Shallow water visualization

[jlchan]

A few users have requested code to plot water height on top of bathymetry. There is a known bug for the default plotting behavior using TreeMesh (see https://github.com/trixi-framework/Trixi.jl/issues/1256), and iplot with Makie doesn't currently support plotting multiple fields on the same plot.

Until these issues are addressed, here is a workaround for creating such a plot:

using GLMakie # use CairoMakie for non-interactive plots (i.e., saving pictures)

# create triangulated plot data for Makie
pd = Trixi.PlotData2DTriangulated(sol)

# plot bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["b"])
fig, ax, plt = Makie.mesh(plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
                          colormap=Trixi.default_Makie_colormap())

# plot water as blue on top of the bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["H"])
Makie.mesh!(ax, plotting_mesh; color=getindex.(plotting_mesh.position, 3), colormap=:blues)

fig # displays the plot

The complete example is below:

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

using GLMakie # use CairoMakie for non-interactive plots (i.e., saving pictures)

# create triangulated plot data for Makie
pd = Trixi.PlotData2DTriangulated(sol)

# plot bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["b"])
fig, ax, plt = Makie.mesh(plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
                          colormap=Trixi.default_Makie_colormap())

# plot water as blue on top of the bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["H"])
Makie.mesh!(ax, plotting_mesh; color=getindex.(plotting_mesh.position, 3), colormap=:blues)

fig # displays the plot

@maxbertrand1996 @svengoldberg this should help in your thesis preparations.

[maxbertrand1996]

If you want to control the opacity of the water surface, the package Colors has to be added to the code posted by @jlchan .

using GLMakie # use CairoMakie for non-interactive plots (i.e., saving pictures)

# create triangulated plot data for Makie
pd = Trixi.PlotData2DTriangulated(sol)

# plot bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["b"])
fig, ax, plt = Makie.mesh(plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
                          colormap=Trixi.default_Makie_colormap())

# plot water as blue on top of the bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["H"])
Makie.mesh!(ax, plotting_mesh; color=getindex.(plotting_mesh.position, 3), colormap=:blues)

fig # displays the plot

Afterwards the colormap attribute has du be extended by RGBAf0.(Colors.color.(to_colormap(:blues)), 0.6). This leads to the following extension to the previous code

using GLMakie # use CairoMakie for non-interactive plots (i.e., saving pictures)
using Colors

# create triangulated plot data for Makie
pd = Trixi.PlotData2DTriangulated(sol)

# plot bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["b"])
fig, ax, plt = Makie.mesh(plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
                          colormap=Trixi.default_Makie_colormap())

# plot water as blue on top of the bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["H"])
Makie.mesh!(ax, plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
            colormap = RGBAf0.(Colors.color.(to_colormap(:blues)), 0.6))

fig # displays the plot

The number in (Colors.color.(to_colormap(:blues)), 0.6) controls the opacity and can be set to a value between 1.0 and 0.0, where 1.0 is full opacity and 0.0 is no opacity.

[sloede]

@maxbertrand1996 Interesting! Can you maybe post an example with and without opacity?

[maxbertrand1996]

Sure

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

using GLMakie # use CairoMakie for non-interactive plots (i.e., saving pictures)
using Colors

# create triangulated plot data for Makie
pd = Trixi.PlotData2DTriangulated(sol)

# plot bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["b"])
fig, ax, plt = Makie.mesh(plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
                          colormap=Trixi.default_Makie_colormap())

# plot water as blue on top of the bathymetry
plotting_mesh = Trixi.global_plotting_triangulation_makie(pd["H"])
Makie.mesh!(ax, plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
            colormap = RGBAf0.(Colors.color.(to_colormap(:blues)), 1.0))

fig # displays the plot

This is the resulting plot

Setting the opacity value to 0.1 gives the following plot

If you set the opacity value to 0.0, you would not see a water surface at all.

[maxbertrand1996]

A little addendum. Makie does no longer support RGBAf0. Now RGBAf has to be used.

This means that the last section of the code

Makie.mesh!(ax, plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
           colormap = RGBAf0.(Colors.color.(to_colormap(:blues)), 1.0))

fig # displays the plot

now has to be

Makie.mesh!(ax, plotting_mesh; color=getindex.(plotting_mesh.position, 3), 
            colormap = RGBAf.(Colors.color.(to_colormap(:blues)), 1.0))

fig # displays the plot

[jlchan]

Thanks for catching that. Does this break the existing iplot functionality in Trixi?

[maxbertrand1996]

The default example still works with iplot