PlotData2D non-periodic TreeMesh

[maxbertrand1996]

When working with the two dimensional shallow water equations and visualizing the bottom topography b with PlotData2D, the resulting plot can look not as expected.

Taking for example this code

using Plots
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the shallow water equations with a discontinuous
# bottom topography function

equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=3.25)

# An initial condition with constant total water height and zero velocities to test well-balancedness.
# Note, this routine is used to compute errors in the analysis callback but the initialization is
# overwritten by `initial_condition_discontinuous_well_balancedness` below.
function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D)
  # Set the background values
  H = equations.H0
  v1 = 0.0
  v2 = 0.0
  # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
  x1, x2 = x
  b = (  1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) )
       + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) )
  return prim2cons(SVector(H, v1, v2, b), equations)
end

initial_condition = initial_condition_well_balancedness

###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
solver = DGSEM(polydeg=4, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# Get the TreeMesh and setup a periodic mesh

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=10_000)

# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solver

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# Callbacks

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(lake_at_rest_error,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true)

stepsize_callback = StepsizeCallback(cfl=3.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

pd = PlotData2D(sol)

pyplot()
surface(pd["b"], title="Bottom topography")

where at the end the bottom topography function

\begin{aligned}
b(x,y) = (1.5 / \exp( 0.5 \cdot ((x - 1.0)^2 + (y - 1.0)^2) ) + 0.75 / \exp( 0.5 \cdot ((x + 1.0)^2 + (y + 1.0)^2) ) )
\end{aligned}

of the solution is displayed via the following surface plot.

You would not expect these 'walls' at the borders of the plot, but rather the follwoing:

It seems that PlotData2D always assumes that the data is periodic.

[ranocha]

From @sloede on Slack:

for TreeMesh and PlotData2D, the visualization routines [..] just assume a periodic mesh in all directions

From the docstring for cell2node:

Convert cell-centered values to node-centered values by averaging over all four neighbors and making use of the periodicity of the solution

(highlight is mine) https://github.com/trixi-framework/Trixi.jl/blob/b3e4f1542e33921af390641d0090e61c773cb10c/src/visualization/utilities.jl#L445-L496 cell2node is called from get_data_2d, which is called from PlotData2DCartesian, which is ultimately called by PlotData2D when called with a TreeMesh mesh: https://github.com/trixi-framework/Trixi.jl/blob/b3e4f1542e33921af390641d0090e61c773cb10c/src/visualization/types.jl#L239 Thus, I think someone would need to polish this up or at the very least add a comment somewhere more prominent (e.g., in the visualization docs and the docstring for PlotData2D)

[jlchan]

@maxbertrand1996 out of curiosity, does using pd = PlotData2DTriangulated(sol) produce a different plot?

[maxbertrand1996]

It says that PlotData2DTriangulated is not defined

ERROR: LoadError: UndefVarError: PlotData2DTriangulated not defined
Stacktrace:
 [1] top-level scope
   @ C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\examples\tree_2d_dgsem\PlotData2D_error.jl:82
 [2] include(fname::String)
   @ Base.MainInclude .\client.jl:451
 [3] top-level scope
   @ REPL[1]:1

Note: I updated my local repo just before this run to the latest version of Trixi.jl

[jlchan]

Oops, you may need Trixi.PlotData2DTriangulated because it's not exported directly.

[sloede]

It is always helpful to show the command + error. Have you tried Trixi.PlotData2DTriangulated? EDIT: Sorry, I guess @jlchan was just a little too quick for me 😬

[maxbertrand1996]

using Plots
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the shallow water equations with a discontinuous
# bottom topography function

equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=3.25)

# An initial condition with constant total water height and zero velocities to test well-balancedness.
# Note, this routine is used to compute errors in the analysis callback but the initialization is
# overwritten by `initial_condition_discontinuous_well_balancedness` below.
function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D)
  # Set the background values
  H = equations.H0
  v1 = 0.0
  v2 = 0.0
  # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
  x1, x2 = x
  b = (  1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) )
       + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) )
  return prim2cons(SVector(H, v1, v2, b), equations)
end

initial_condition = initial_condition_well_balancedness

###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
solver = DGSEM(polydeg=4, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# Get the TreeMesh and setup a periodic mesh

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=10_000)

# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solver

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# Callbacks

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(lake_at_rest_error,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true)

stepsize_callback = StepsizeCallback(cfl=3.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

pd = Trixi.PlotData2DTriangulated(sol)

pyplot()
surface(pd["b"], title="Bottom topography")

This results in a heatmap and not a surface plot.

[maxbertrand1996]

But at least the boundaries seem more reasonable

[jlchan]

Ah right. We support surface plots in Makie.jl, but unfortunately when using Plots.jl we only have access to heatmaps due to the use of TriplotRecipes.jl.

[sloede]

@maxbertrand1996 have you tried using Makie, e.g., via the GLMakie package?

[maxbertrand1996]

Isn't this only supprorted for UnstructuredMesh2D? See documentation

I am using TreeMesh

[maxbertrand1996]

I still tried it using this code

using GLMakie
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the shallow water equations with a discontinuous
# bottom topography function

equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=3.25)

# An initial condition with constant total water height and zero velocities to test well-balancedness.
# Note, this routine is used to compute errors in the analysis callback but the initialization is
# overwritten by `initial_condition_discontinuous_well_balancedness` below.
function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D)
  # Set the background values
  H = equations.H0
  v1 = 0.0
  v2 = 0.0
  # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
  x1, x2 = x
  b = (  1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) )
       + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) )
  return prim2cons(SVector(H, v1, v2, b), equations)
end

initial_condition = initial_condition_well_balancedness

###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
solver = DGSEM(polydeg=4, surface_flux=surface_flux,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# Get the TreeMesh and setup a periodic mesh

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=10_000)

# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solver

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# Callbacks

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(lake_at_rest_error,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true)

stepsize_callback = StepsizeCallback(cfl=3.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

pd = Trixi.PlotData2DTriangulated(sol)

iplot(sol)

And got this error message

ERROR: LoadError: ArgumentError: range step cannot be zero
Stacktrace:
  [1] (::Colon)(start::Float32, step::Float32, stop::Float32)
    @ Base .\twiceprecision.jl:401
  [2] get_minor_tickvalues(i::IntervalsBetween, scale::Function, tickvalues::Vector{Float32}, vmin::Float32, vmax::Float32)
    @ Makie.MakieLayout C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\lineaxis.jl:606
  [3] (::Makie.MakieLayout.var"#182#213"{Observable{Vector{Float32}}, Observable{Any}, Attributes})(tickvalues::Vector{Float32}, minorticks::IntervalsBetween)
    @ Makie.MakieLayout C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\lineaxis.jl:239
  [4] (::Observables.OnUpdate{Makie.MakieLayout.var"#182#213"{Observable{Vector{Float32}}, Observable{Any}, Attributes}, Tuple{Observable{Vector{Float32}}, Observable{Any}}})(#unused#::Vector{Float32})
    @ Observables C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:334
  [5] #invokelatest#2
    @ .\essentials.jl:716 [inlined]
  [6] invokelatest
    @ .\essentials.jl:714 [inlined]
  [7] notify
    @ C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:88 [inlined]
  [8] setindex!(observable::Observable{Vector{Float32}}, val::Vector{Float64})
    @ Observables C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:248
  [9] (::Makie.MakieLayout.var"#178#209"{Observable{Vector{AbstractString}}, Observable{Vector{Point{2, Float32}}}, Observable{Vector{Float32}}, Observable{Tuple{Float32, Tuple{Float32, Float32}, Bool}}, Observable{Any}, Attributes})(tickvalues_labels_unfiltered::Tuple{Vector{Float64}, Vector{String}}, reversed::Bool)
    @ Makie.MakieLayout C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\lineaxis.jl:214
 [10] (::Observables.OnUpdate{Makie.MakieLayout.var"#178#209"{Observable{Vector{AbstractString}}, Observable{Vector{Point{2, Float32}}}, Observable{Vector{Float32}}, Observable{Tuple{Float32, Tuple{Float32, Float32}, Bool}}, Observable{Any}, Attributes}, Tuple{Observable{Tuple{Vector{Float64}, Vector{String}}}, Observable{Any}}})(#unused#::Tuple{Vector{Float64}, Vector{String}})
    @ Observables C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:334
 [11] #invokelatest#2
    @ .\essentials.jl:716 [inlined]
 [12] invokelatest
    @ .\essentials.jl:714 [inlined]
 [13] notify
    @ C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:88 [inlined]
 [14] setindex!(observable::Observable{Tuple{Vector{Float64}, Vector{String}}}, val::Tuple{Vector{Float64}, Vector{String}})
    @ Observables C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:248
 [15] (::Observables.MapUpdater{Makie.MakieLayout.var"#177#208", Tuple{Vector{Float64}, Vector{String}}})(::Tuple{Float32, Tuple{Float32, Float32}, Bool}, ::Vararg{Any})
    @ Observables C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:372
 [16] (::Observables.OnUpdate{Observables.MapUpdater{Makie.MakieLayout.var"#177#208", Tuple{Vector{Float64}, Vector{String}}}, Tuple{Observable{Tuple{Float32, Tuple{Float32, Float32}, Bool}}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}}})(#unused#::Tuple{Float32, Float32})
    @ Observables C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:334
 [17] #invokelatest#2
    @ .\essentials.jl:716 [inlined]
 [18] invokelatest
    @ .\essentials.jl:714 [inlined]
 [19] notify
    @ C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:88 [inlined]
 [20] setindex!(observable::Observable{Any}, val::Tuple{Float32, Float32})
    @ Observables C:\Users\maxbe\.julia\packages\Observables\OFj0u\src\Observables.jl:248
 [21] Makie.MakieLayout.LineAxis(parent::Scene; kwargs::Base.Pairs{Symbol, Any, NTuple{35, Symbol}, NamedTuple{(:endpoints, :flipped, :limits, :ticklabelalign, :label, :labelpadding, :labelvisible, :labelsize, :labelcolor, :labelfont, :ticklabelfont, :ticks, :tickformat, :ticklabelsize, :ticklabelsvisible, :ticksize, :ticksvisible, :ticklabelpad, :tickalign, :ticklabelrotation, :tickwidth, :tickcolor, :spinewidth, :ticklabelspace, :ticklabelcolor, :spinecolor, :spinevisible, :flip_vertical_label, :minorticksvisible, :minortickalign, :minorticksize, :minortickwidth, :minortickcolor, :minorticks, :scale), Tuple{Observable{Tuple{Point{2, Float32}, Point{2, Float32}}}, Observable{Any}, Observable{Tuple{Float32, Float32}}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Symbol, Bool, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}, Observable{Any}}}})
    @ Makie.MakieLayout C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\lineaxis.jl:372
 [22] layoutable(::Type{Colorbar}, fig_or_scene::Figure; bbox::Nothing, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:limits, :colormap), Tuple{Observable{Tuple{Float32, Float32}}, Symbol}}})
    @ Makie.MakieLayout C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\layoutables\colorbar.jl:330
 [23] #_layoutable#11
    @ C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\layoutables.jl:69 [inlined]
 [24] _layoutable(::Type{Colorbar}, ::GridPosition; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:limits, :colormap), Tuple{Observable{Tuple{Float32, Float32}}, Symbol}}})
    @ Makie.MakieLayout C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\layoutables.jl:64
 [25] #_#9
    @ C:\Users\maxbe\.julia\packages\Makie\xbI6d\src\makielayout\layoutables.jl:49 [inlined]
 [26] iplot(pd::Trixi.PlotData2DTriangulated{StructArray{SVector{4, Float64}, 2, NTuple{4, Matrix{Float64}}, Int64}, Matrix{Float64}, Matrix{Float64}, StructArray{SVector{4, Float64}, 2, NTuple{4, Matrix{Float64}}, Int64}, SVector{4, String}, Matrix{Int32}}; plot_mesh::Bool, show_axis::Bool, colormap::Symbol, variable_to_plot_in::Int64)
    @ Trixi C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\src\visualization\recipes_makie.jl:106
in expression starting at C:\Users\maxbe\Repos\maxbertrand1996\Trixi.jl\examples\tree_2d_dgsem\PlotData2D_error.jl:84

[sloede]

Hm, dunno why this is happening. If I do

using Trixi
trixi_include(default_example())
using GLMakie
iplot(sol)

then I get this:

[maxbertrand1996]

Does it also work, if you execute my code?

[sloede]

Does it also work, if you execute my code?

Hm, no. But there it already doesn't work for plot(sol), I don't even have to go to iplot(sol) :-/

@jlchan I can reproduce the issue locally. It seems to be caused by faulty input to https://github.com/trixi-framework/Trixi.jl/blob/b3e4f1542e33921af390641d0090e61c773cb10c/src/visualization/recipes_makie.jl#L251 although I cannot figure out what the issue is (it is a "normal" TreeMesh2D without AMR etc. as far as I can tell). Do you have any ideas? The only thing I can think of is that there are abnormal triangles or something similar that cannot be handled properly (e.g., three points almost on a single line). But here I am at the end of my capacity...

[sloede]

OK, getting closer. If I comment the line I just referenced https://github.com/trixi-framework/Trixi.jl/blob/b3e4f1542e33921af390641d0090e61c773cb10c/src/visualization/recipes_makie.jl#L251 then plot(sol) works. If I also comment https://github.com/trixi-framework/Trixi.jl/blob/b3e4f1542e33921af390641d0090e61c773cb10c/src/visualization/recipes_makie.jl#L106 iplot(sol) works. However, I believe here is the problem: Because H does not vary, i.e.,

julia> pd = PlotData2D(sol)
PlotData2DCartesian{Vector{Float64},Vector{Matrix{Float64}},SVector{4, String},Vector{Float64}}(<x>, <y>, <data>, <variable_names>, <mesh_vertices_x>, <mesh_vertices_y>)

julia> mini, maxi = extrema(pd.data[1])
(3.2499999999999942, 3.2500000000000053)

julia> mini ≈ maxi
true

the Colorbar function fails, which seems to be a known issue for Makie: https://github.com/MakieOrg/Makie.jl./issues/1400

[jlchan]

Thanks for looking into that. That Makie colorbar issue is rather annoying.

What do you think a good fix is? I am thinking we could turn the colorbar on/off via keyword option...

[sloede]

I thought something like this could be a possible fix, i.e., checking for approximate equality of the extrema and then just adding a offset in this case:

  mini, maxi = extrema(solution_z)
  if isapprox(mini, maxi)
    _one = one(eltype(solution_z))
    Makie.Colorbar(fig[1, 3], limits = Makie.@lift(extrema($solution_z) .+ (-_one, _one)), colormap=colormap)
  else
    Makie.Colorbar(fig[1, 3], limits = Makie.@lift(extrema($solution_z)), colormap=colormap)
  end

However, this does not work, since in the function where this is executed, solution_z is not actually available - the @lift macro just creates closures apparently, thus I do not know how to implement such logic based on actual solution values.

I also tried to see if we can get by using just a "small" offset, e.g.,

Makie.Colorbar(fig[1, 3], limits = Makie.@lift(extrema($solution_z) .+ (-1.0e-8, 1.0e-8)), colormap=colormap)

but as you can see, small is not really small (1e-8 might not be negligible in certain situations).

@jlchan would you know how to implement this using Makie's recipe style?

[jlchan]

I wonder if we can include the extreme call inside the @lift closure? We could scale by 1e-8 * the magnitude of the solution.

I can take a look at this later today and see if this works.

[maxbertrand1996]

Sorry I was not contributing to the discussion the past days. But is there anything I can do for now? I do not have any experience with the Makie code.

[jlchan]

Apologies - it's been busy, and I just got a chance to look again at this today. I'm adding a PR which should fix this issue.

[jlchan]

https://github.com/trixi-framework/Trixi.jl/pull/1259 should fix this with a scaling.

[maxbertrand1996]

Can I close this issue now?

[ranocha]

I guess only iplot is fixed, not your original issue? So we should keep it open

[maxbertrand1996]

Yes, you are right.

[jlchan]

I posted a workaround in https://github.com/trixi-framework/Trixi.jl/issues/1265. However, we should leave this issue open since the original problem of Plot2DCartesian is still not resolved.