7: DG schemes via DGMulti solver
DGMulti is a DG solver that allows meshes with simplex elements. The basic idea and implementation of this solver is explained in section "Meshes". Here, we want to give some examples and a quick overview about the options with DGMulti.
We start with a simple example we already used in the tutorial about flux differencing. There, we implemented a simulation with initial_condition_weak_blast_wave for the 2D compressible Euler equations CompressibleEulerEquations2D and used the DG formulation with flux differencing using volume flux flux_ranocha and surface flux flux_lax_friedrichs.
Here, we want to implement the equivalent example, only now using the DGMulti solver instead of DGSEM.
using Trixi, OrdinaryDiffEq
equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_weak_blast_waveinitial_condition_weak_blast_wave (generic function with 13 methods)To use the Gauss-Lobatto nodes again, we choose approximation_type=SBP() in the solver. Since we want to start with a Cartesian domain discretized with squares, we use the element type Quad().
dg = DGMulti(polydeg = 3,
element_type = Quad(),
approximation_type = SBP(),
surface_flux = flux_lax_friedrichs,
volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
cells_per_dimension = (32, 32)
mesh = DGMultiMesh(dg,
cells_per_dimension, # initial_refinement_level = 5
coordinates_min=(-2.0, -2.0),
coordinates_max=( 2.0, 2.0),
periodicity=true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
boundary_conditions=boundary_condition_periodic)
tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)
alive_callback = AliveCallback(alive_interval=10)
analysis_callback = AnalysisCallback(semi, interval=100, uEltype=real(dg))
callbacks = CallbackSet(analysis_callback, alive_callback);Run the simulation with the same time integration algorithm as before.
sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,
callback=callbacks, save_everystep=false);
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 9.92000000e-07 s
Δt: 0.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: NaN s
PID: Inf s
#DOFs per field: 16384 alloc'd memory: 3103.889 MiB
#elements: 1024
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
Linf error: 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
∑∂S/∂U ⋅ Uₜ : -4.07691781e-17
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 10 │ Δt: 8.7782e-03 │ sim. time: 6.2446e-02 (15.612%) │ run time: 1.4644e-01 s
#timesteps: 20 │ Δt: 1.2426e-02 │ sim. time: 1.7307e-01 (43.267%) │ run time: 2.8826e-01 s
#timesteps: 30 │ Δt: 1.3710e-02 │ sim. time: 3.0622e-01 (76.554%) │ run time: 4.3107e-01 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 37 run time: 5.34031689e-01 s
Δt: 9.57329122e-03 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 4.00000000e-01 (100.000%) time/DOF/rhs!: 8.97765237e-08 s
PID: 9.65199203e-08 s
#DOFs per field: 16384 alloc'd memory: 3108.750 MiB
#elements: 1024
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 6.13970440e-02 4.96516886e-02 4.96432472e-02 2.24583683e-01
Linf error: 2.61815840e-01 2.48009699e-01 2.47516396e-01 9.30972703e-01
∑∂S/∂U ⋅ Uₜ : -2.24491554e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 0.4 Time steps: 37 (accepted), 37 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────using Plots
pd = PlotData2D(sol)
plot(pd)
plot(pd["rho"])
plot!(getmesh(pd))
This simulation is not as fast as the equivalent with TreeMesh since no special optimizations for quads (for instance tensor product structure) have been implemented. Figure 4 in "Efficient implementation of modern entropy stable and kinetic energy preserving discontinuous Galerkin methods for conservation laws" (2021) provides a nice runtime comparison between the different mesh types. On the other hand, the functions are more general and thus we have more option we can choose from.
Simulation with Gauss nodes
For instance, we can change the approximation type of our simulation.
using Trixi, OrdinaryDiffEq
equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_weak_blast_waveinitial_condition_weak_blast_wave (generic function with 13 methods)We now use Gauss nodes instead of Gauss-Lobatto nodes which can be done for the element types Quad() and Hex(). Therefore, we set approximation_type=GaussSBP(). Alternatively, we can use a modal approach using the approximation type Polynomial().
dg = DGMulti(polydeg = 3,
element_type = Quad(),
approximation_type = GaussSBP(),
surface_flux = flux_lax_friedrichs,
volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
cells_per_dimension = (32, 32)
mesh = DGMultiMesh(dg,
cells_per_dimension, # initial_refinement_level = 5
coordinates_min=(-2.0, -2.0),
coordinates_max=( 2.0, 2.0),
periodicity=true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
boundary_conditions=boundary_condition_periodic)
tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)
alive_callback = AliveCallback(alive_interval=10)
analysis_callback = AnalysisCallback(semi, interval=100, uEltype=real(dg))
callbacks = CallbackSet(analysis_callback, alive_callback);
sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,
ode_default_options()..., callback=callbacks);
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 1.65300000e-06 s
Δt: 0.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: NaN s
PID: Inf s
#DOFs per field: 16384 alloc'd memory: 3158.490 MiB
#elements: 1024
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 3.17057848e-16 8.94192483e-18 9.07727588e-18 5.97871667e-16
Linf error: 8.88178420e-16 1.38777878e-16 1.38777878e-16 2.66453526e-15
∑∂S/∂U ⋅ Uₜ : -1.01164379e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 10 │ Δt: 5.0548e-03 │ sim. time: 3.3436e-02 (8.359%) │ run time: 4.7967e-01 s
#timesteps: 20 │ Δt: 8.8313e-03 │ sim. time: 1.0601e-01 (26.503%) │ run time: 9.4567e-01 s
#timesteps: 30 │ Δt: 1.0521e-02 │ sim. time: 2.0438e-01 (51.094%) │ run time: 1.4329e+00 s
#timesteps: 40 │ Δt: 1.1006e-02 │ sim. time: 3.1245e-01 (78.113%) │ run time: 1.9021e+00 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 48 run time: 2.28545686e+00 s
Δt: 1.00192333e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 4.00000000e-01 (100.000%) time/DOF/rhs!: 3.04370725e-07 s
PID: 3.19778189e-07 s
#DOFs per field: 16384 alloc'd memory: 3160.860 MiB
#elements: 1024
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 6.10354479e-02 4.94109888e-02 4.94109888e-02 2.23252171e-01
Linf error: 2.58968601e-01 2.43227118e-01 2.43227118e-01 9.39282114e-01
∑∂S/∂U ⋅ Uₜ : -2.32263587e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 0.4 Time steps: 48 (accepted), 48 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────using Plots
pd = PlotData2D(sol)
plot(pd)
Simulation with triangular elements
Also, we can set another element type. We want to use triangles now.
using Trixi, OrdinaryDiffEq
equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_weak_blast_waveinitial_condition_weak_blast_wave (generic function with 13 methods)Since there is no direct equivalent to Gauss-Lobatto nodes on triangles, the approximation type SBP() now uses Gauss-Lobatto nodes on faces and special nodes in the interior of the triangular. More details can be found in the documentation of StartUpDG.jl.
dg = DGMulti(polydeg = 3,
element_type = Tri(),
approximation_type = SBP(),
surface_flux = flux_lax_friedrichs,
volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
cells_per_dimension = (32, 32)
mesh = DGMultiMesh(dg,
cells_per_dimension, # initial_refinement_level = 5
coordinates_min=(-2.0, -2.0),
coordinates_max=( 2.0, 2.0),
periodicity=true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
boundary_conditions=boundary_condition_periodic)
tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)
alive_callback = AliveCallback(alive_interval=10)
analysis_callback = AnalysisCallback(semi, interval=100, uEltype=real(dg))
callbacks = CallbackSet(analysis_callback, alive_callback);
sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,
ode_default_options()..., callback=callbacks);
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 9.62000000e-07 s
Δt: 0.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: NaN s
PID: Inf s
#DOFs per field: 30720 alloc'd memory: 3198.959 MiB
#elements: 2048
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
Linf error: 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
∑∂S/∂U ⋅ Uₜ : -6.74736856e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 10 │ Δt: 6.6947e-03 │ sim. time: 4.5960e-02 (11.490%) │ run time: 8.4028e-01 s
#timesteps: 20 │ Δt: 1.0524e-02 │ sim. time: 1.3699e-01 (34.249%) │ run time: 1.6419e+00 s
#timesteps: 30 │ Δt: 1.2096e-02 │ sim. time: 2.5197e-01 (62.993%) │ run time: 2.4561e+00 s
#timesteps: 40 │ Δt: 1.2622e-02 │ sim. time: 3.7597e-01 (93.993%) │ run time: 3.2780e+00 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 42 run time: 3.45392812e+00 s
Δt: 1.13592115e-02 └── GC time: 2.48142880e-02 s (0.718%)
sim. time: 4.00000000e-01 (100.000%) time/DOF/rhs!: 2.84264043e-07 s
PID: 2.94181248e-07 s
#DOFs per field: 30720 alloc'd memory: 3179.859 MiB
#elements: 2048
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 6.07149048e-02 4.91826834e-02 4.91814198e-02 2.22089274e-01
Linf error: 2.69551549e-01 2.45792366e-01 2.45935623e-01 9.34191621e-01
∑∂S/∂U ⋅ Uₜ : -1.95711055e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 0.4 Time steps: 42 (accepted), 42 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────using Plots
pd = PlotData2D(sol)
plot(pd)
plot(pd["rho"])
plot!(getmesh(pd))
Triangular meshes on non-Cartesian domains
To use triangular meshes on a non-Cartesian domain, Trixi.jl uses the package StartUpDG.jl. The following example is based on elixir_euler_triangulate_pkg_mesh.jl and uses a pre-defined mesh from StartUpDG.jl.
using Trixi, OrdinaryDiffEqWe want to simulate the smooth initial condition initial_condition_convergence_test with source terms source_terms_convergence_test for the 2D compressible Euler equations.
equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_testsource_terms_convergence_test (generic function with 13 methods)We create the solver DGMulti with triangular elements (Tri()) as before.
dg = DGMulti(polydeg = 3, element_type = Tri(),
approximation_type=Polynomial(),
surface_flux = flux_lax_friedrichs,
volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: ……………………………………………………………… RefElemData{N=3, Polynomial, Tri}. │
│ mortar: …………………………………………………………… nothing │
│ surface integral: ………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: ……………………………………… FluxLaxFriedrichs(max_abs_speed_naive) │
│ volume integral: …………………………………… VolumeIntegralFluxDifferencing │
│ │ volume flux: ………………………………………… flux_ranocha │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘StartUpDG.jl provides for instance a pre-defined Triangulate geometry for a rectangular domain with hole RectangularDomainWithHole. Other pre-defined Triangulate geometries are e.g., SquareDomain, Scramjet, and CircularDomain.
meshIO = StartUpDG.triangulate_domain(StartUpDG.RectangularDomainWithHole());The pre-defined Triangulate geometry in StartUpDG has integer boundary tags. With DGMultiMesh we assign boundary faces based on these integer boundary tags and create a mesh compatible with Trixi.jl.
mesh = DGMultiMesh(dg, meshIO, Dict(:outer_boundary=>1, :inner_boundary=>2))┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DGMultiMesh{2, Trixi.Affine}, │
│ ══════════════════════════════ │
│ number of elements: …………………………… 598 │
│ number of boundaries: ……………………… 2 │
│ │ nfaces on outer_boundary: ……… 52 │
│ │ nfaces on inner_boundary: ……… 4 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = (; :outer_boundary => boundary_condition_convergence_test,
:inner_boundary => boundary_condition_convergence_test)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
source_terms = source_terms,
boundary_conditions = boundary_conditions)
tspan = (0.0, 0.2)
ode = semidiscretize(semi, tspan)
alive_callback = AliveCallback(alive_interval=20)
analysis_callback = AnalysisCallback(semi, interval=200, uEltype=real(dg))
callbacks = CallbackSet(alive_callback, analysis_callback);
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 1.19200000e-06 s
Δt: 1.49245207e-03 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: NaN s
PID: Inf s
#DOFs per field: 5980 alloc'd memory: 3226.828 MiB
#elements: 598
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 3.93036845e-07 3.93036845e-07 3.93036845e-07 1.55064245e-06
Linf error: 4.47141181e-06 4.47141181e-06 4.47141181e-06 1.47933676e-05
∑∂S/∂U ⋅ Uₜ : -1.23444983e-02
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 20 │ Δt: 1.4925e-03 │ sim. time: 2.9849e-02 (14.925%) │ run time: 5.3060e-01 s
#timesteps: 40 │ Δt: 1.4925e-03 │ sim. time: 5.9698e-02 (29.849%) │ run time: 1.0461e+00 s
#timesteps: 60 │ Δt: 1.4925e-03 │ sim. time: 8.9547e-02 (44.774%) │ run time: 1.5613e+00 s
#timesteps: 80 │ Δt: 1.4925e-03 │ sim. time: 1.1940e-01 (59.698%) │ run time: 2.0860e+00 s
#timesteps: 100 │ Δt: 1.4925e-03 │ sim. time: 1.4925e-01 (74.623%) │ run time: 2.6126e+00 s
#timesteps: 120 │ Δt: 1.4925e-03 │ sim. time: 1.7909e-01 (89.547%) │ run time: 3.1286e+00 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 0.2 Time steps: 135 (accepted), 135 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 135 run time: 3.51601787e+00 s
Δt: 1.14223386e-05 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 2.00000000e-01 (100.000%) time/DOF/rhs!: 8.61506334e-07 s
PID: 8.67846343e-07 s
#DOFs per field: 5980 alloc'd memory: 3227.301 MiB
#elements: 598
Variable: rho rho_v1 rho_v2 rho_e
L2 error: 4.43758204e-06 3.71131543e-06 4.42671638e-06 9.85291883e-06
Linf error: 6.10809367e-05 4.53698746e-05 6.62331341e-05 1.25625347e-04
∑∂S/∂U ⋅ Uₜ : -1.03046341e-02
────────────────────────────────────────────────────────────────────────────────────────────────────using Plots
pd = PlotData2D(sol)
plot(pd["rho"])
plot!(getmesh(pd))
For more information, please have a look in the StartUpDG.jl documentation.
Package versions
These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["Trixi", "StartUpDG", "OrdinaryDiffEq", "Plots"],
mode=PKGMODE_MANIFEST)Julia Version 1.9.4
Commit 8e5136fa297 (2023-11-14 08:46 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 7763 64-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-14.0.6 (ORCJIT, znver3)
Threads: 1 on 4 virtual cores
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/Trixi.jl/Trixi.jl/docs/Manifest.toml`
⌅ [1dea7af3] OrdinaryDiffEq v6.66.0
[91a5bcdd] Plots v1.40.4
[472ebc20] StartUpDG v0.17.7
[a7f1ee26] Trixi v0.7.8 `~/work/Trixi.jl/Trixi.jl`
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`This page was generated using Literate.jl.