8: 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_wave
initial_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:       7.82000000e-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:       1140.831 MiB
 #elements:               16384

 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ₜ :  -8.96621432e-17
────────────────────────────────────────────────────────────────────────────────────────────────────

#timesteps:     10 │ Δt: 8.7782e-03 │ sim. time: 6.2446e-02 (15.612%)  │ run time: 1.4625e-01 s
#timesteps:     20 │ Δt: 1.2426e-02 │ sim. time: 1.7307e-01 (43.267%)  │ run time: 2.8822e-01 s
#timesteps:     30 │ Δt: 1.3710e-02 │ sim. time: 3.0622e-01 (76.554%)  │ run time: 4.4698e-01 s

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 37                run time:       5.50386058e-01 s
 Δt:             9.57329122e-03                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      4.00000000e-01 (100.000%)     time/DOF/rhs!:  9.24932783e-08 s
                                               PID:            9.94337898e-08 s
 #DOFs per field:         16384                alloc'd memory:       1150.638 MiB
 #elements:               16384

 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)
Example block output
plot(pd["rho"])
plot!(getmesh(pd))
Example block output

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_wave
initial_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.04200000e-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:        813.325 MiB
 #elements:               16384

 Variable:       rho              rho_v1           rho_v2           rho_e
 L2 error:       3.31396240e-15   1.31101638e-16   1.30776731e-16   8.29084420e-15
 Linf error:     7.32747196e-15   3.13037310e-15   3.20427569e-15   1.99840144e-14
 ∑∂S/∂U ⋅ Uₜ :  -1.01164379e-01
────────────────────────────────────────────────────────────────────────────────────────────────────

#timesteps:     10 │ Δt: 5.0548e-03 │ sim. time: 3.3436e-02 (8.359%)   │ run time: 4.7266e-01 s
#timesteps:     20 │ Δt: 8.8313e-03 │ sim. time: 1.0601e-01 (26.503%)  │ run time: 9.4942e-01 s
#timesteps:     30 │ Δt: 1.0521e-02 │ sim. time: 2.0438e-01 (51.094%)  │ run time: 1.4109e+00 s
#timesteps:     40 │ Δt: 1.1006e-02 │ sim. time: 3.1245e-01 (78.113%)  │ run time: 1.8761e+00 s

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 48                run time:       2.25618988e+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.01897874e-07 s
                                               PID:            3.15617011e-07 s
 #DOFs per field:         16384                alloc'd memory:        820.837 MiB
 #elements:               16384

 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)
Example block output

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_wave
initial_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:       7.81000000e-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:        857.188 MiB
 #elements:               30720

 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ₜ :  -7.93719284e-03
────────────────────────────────────────────────────────────────────────────────────────────────────

#timesteps:     10 │ Δt: 6.6587e-03 │ sim. time: 4.5551e-02 (11.388%)  │ run time: 7.7240e-01 s
#timesteps:     20 │ Δt: 1.0494e-02 │ sim. time: 1.3622e-01 (34.054%)  │ run time: 1.5280e+00 s
#timesteps:     30 │ Δt: 1.2082e-02 │ sim. time: 2.5106e-01 (62.764%)  │ run time: 2.3098e+00 s
#timesteps:     40 │ Δt: 1.2625e-02 │ sim. time: 3.7504e-01 (93.761%)  │ run time: 3.1020e+00 s

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'CompressibleEulerEquations2D' with DGMulti(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                 42                run time:       3.26829446e+00 s
 Δt:             1.22969880e-02                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      4.00000000e-01 (100.000%)     time/DOF/rhs!:  2.71147187e-07 s
                                               PID:            2.78351510e-07 s
 #DOFs per field:         30720                alloc'd memory:        871.241 MiB
 #elements:               30720

 Variable:       rho              rho_v1           rho_v2           rho_e
 L2 error:       6.07273896e-02   4.91888840e-02   4.91932276e-02   2.22135410e-01
 Linf error:     2.69688605e-01   2.45714816e-01   2.45652668e-01   9.33596710e-01
 ∑∂S/∂U ⋅ Uₜ :  -1.95902350e-03
────────────────────────────────────────────────────────────────────────────────────────────────────

────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished.  Final time: 0.4  Time steps: 42 (accepted), 42 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
using Plots
pd = PlotData2D(sol)
plot(pd)
Example block output
plot(pd["rho"])
plot!(getmesh(pd))
Example block output

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, OrdinaryDiffEq

We 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_test
source_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:       8.91000000e-07 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:        860.435 MiB
 #elements:                5980

 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.4721e-01 s
#timesteps:     40 │ Δt: 1.4925e-03 │ sim. time: 5.9698e-02 (29.849%)  │ run time: 1.0818e+00 s
#timesteps:     60 │ Δt: 1.4925e-03 │ sim. time: 8.9547e-02 (44.774%)  │ run time: 1.6161e+00 s
#timesteps:     80 │ Δt: 1.4925e-03 │ sim. time: 1.1940e-01 (59.698%)  │ run time: 2.1502e+00 s
#timesteps:    100 │ Δt: 1.4925e-03 │ sim. time: 1.4925e-01 (74.623%)  │ run time: 2.6851e+00 s
#timesteps:    120 │ Δt: 1.4925e-03 │ sim. time: 1.7909e-01 (89.547%)  │ run time: 3.2197e+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.62148036e+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.85740292e-07 s
                                               PID:            8.94290938e-07 s
 #DOFs per field:          5980                alloc'd memory:        867.745 MiB
 #elements:                5980

 Variable:       rho              rho_v1           rho_v2           rho_e
 L2 error:       4.43758204e-06   3.71131542e-06   4.42671638e-06   9.85291883e-06
 Linf error:     6.10809368e-05   4.53698746e-05   6.62331342e-05   1.25625347e-04
 ∑∂S/∂U ⋅ Uₜ :  -1.03046341e-02
────────────────────────────────────────────────────────────────────────────────────────────────────
using Plots
pd = PlotData2D(sol)
plot(pd["rho"])
plot!(getmesh(pd))
Example block output

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.10.6
Commit 67dffc4a8ae (2024-10-28 12:23 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-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (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.9
  [472ebc20] StartUpDG v1.1.5
  [a7f1ee26] Trixi v0.9.7 `~/work/Trixi.jl/Trixi.jl`
Info Packages marked with  have new versions available and may be upgradable.

This page was generated using Literate.jl.