6: Subcell limiting with the IDP Limiter

In the previous tutorial, the element-wise limiting with IndicatorHennemannGassner and VolumeIntegralShockCapturingHG was explained. This tutorial contains a short introduction to the idea and implementation of subcell shock capturing approaches in Trixi.jl, which is also based on the DGSEM scheme in flux differencing formulation. Trixi.jl contains the a-posteriori invariant domain-preserving (IDP) limiter which was introduced by Pazner (2020) and Rueda-Ramírez, Pazner, Gassner (2022). It is a flux-corrected transport-type (FCT) limiter and is implemented using SubcellLimiterIDP and VolumeIntegralSubcellLimiting. Since it is an a-posteriori limiter you have to apply a correction stage after each Runge-Kutta stage. This is done by passing the stage callback SubcellLimiterIDPCorrection to the time integration method.

Time integration method

As mentioned before, the IDP limiting is an a-posteriori limiter. Its limiting process guarantees the target bounds for an explicit (forward) Euler time step. To still achieve a high-order approximation, the implementation uses strong-stability preserving (SSP) Runge-Kutta methods, which can be written as convex combinations of forward Euler steps. As such, they preserve the convexity of convex functions and functionals, such as the TVD semi-norm and the maximum principle in 1D, for instance.

Since IDP/FCT limiting procedure operates on independent forward Euler steps, its a-posteriori correction stage is implemented as a stage callback that is triggered after each forward Euler step in an SSP Runge-Kutta method. Unfortunately, the solve(...) routines in OrdinaryDiffEq.jl, typically employed for time integration in Trixi.jl, do not support this type of stage callback.

Therefore, subcell limiting with the IDP limiter requires the use of a Trixi-intern time integration SSPRK method called with

Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...)

Right now, only the canonical three-stage, third-order SSPRK method (Shu-Osher) Trixi.SimpleSSPRK33 is implemented.

IDP Limiting

The implementation of the invariant domain preserving (IDP) limiting approach (SubcellLimiterIDP) is based on Pazner (2020) and Rueda-Ramírez, Pazner, Gassner (2022). It supports several types of limiting which are enabled by passing parameters individually.

Global bounds

If enabled, the global bounds enforce physical admissibility conditions, such as non-negativity of variables. This can be done for conservative variables, where the limiter is of a one-sided Zalesak-type (Zalesak, 1979), and general non-linear variables, where a Newton-bisection algorithm is used to enforce the bounds.

The Newton-bisection algorithm is an iterative method and requires some parameters. It uses a fixed maximum number of iteration steps (max_iterations_newton = 10) and relative/absolute tolerances (newton_tolerances = (1.0e-12, 1.0e-14)). The given values are sufficient in most cases and therefore used as default. If the implemented bounds checking functionality indicates problems with the limiting (see below) the Newton method with the chosen parameters might not manage to converge. If so, adapting the mentioned parameters helps fix that. Additionally, there is the parameter gamma_constant_newton, which can be used to scale the antidiffusive flux for the computation of the blending coefficients of nonlinear variables. The default value is 2 * ndims(equations), as it was shown by Pazner (2020) [Section 4.2.2.] that this value guarantees the fulfillment of bounds for a forward-Euler increment.

Very small non-negative values can be an issue as well. That's why we use an additional correction factor in the calculation of the global bounds,

\[u^{new} \geq \beta * u^{FV}.\]

By default, $\beta$ (named positivity_correction_factor) is set to 0.1 which works properly in most of the tested setups.

Conservative variables

The procedure to enforce global bounds for a conservative variables is as follows: If you want to guarantee non-negativity for the density of the compressible Euler equations, you pass the specific quantity name of the conservative variable.

using Trixi
equations = CompressibleEulerEquations2D(1.4)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ CompressibleEulerEquations2D                                                                     │
│ ════════════════════════════                                                                     │
│ #variables: ………………………………………………… 4                                                                │
│ │ variable 1: …………………………………………… rho                                                              │
│ │ variable 2: …………………………………………… rho_v1                                                           │
│ │ variable 3: …………………………………………… rho_v2                                                           │
│ │ variable 4: …………………………………………… rho_e                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

The quantity name of the density is rho which is how we enable its limiting.

positivity_variables_cons = ["rho"]
1-element Vector{String}:
 "rho"

The quantity names are passed as a vector to allow several quantities. This is used, for instance, if you want to limit the density of two different components using the multicomponent compressible Euler equations.

equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
                                                       gas_constants = (0.287, 1.578))
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ CompressibleEulerMulticomponentEquations2D                                                       │
│ ══════════════════════════════════════════                                                       │
│ #variables: ………………………………………………… 5                                                                │
│ │ variable 1: …………………………………………… rho_v1                                                           │
│ │ variable 2: …………………………………………… rho_v2                                                           │
│ │ variable 3: …………………………………………… rho_e                                                            │
│ │ variable 4: …………………………………………… rho1                                                             │
│ │ variable 5: …………………………………………… rho2                                                             │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

Then, we just pass both quantity names.

positivity_variables_cons = ["rho1", "rho2"]
2-element Vector{String}:
 "rho1"
 "rho2"

Alternatively, it is possible to all limit all density variables with a general command using

positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations)]
2-element Vector{String}:
 "rho1"
 "rho2"

Non-linear variables

To allow limitation for all possible non-linear variables, including variables defined on-the-fly, you can directly pass the function that computes the quantity for which you want to enforce positivity. For instance, if you want to enforce non-negativity for the pressure, do as follows.

positivity_variables_nonlinear = [pressure]
1-element Vector{typeof(pressure)}:
 pressure (generic function with 15 methods)

Local bounds

Second, Trixi.jl supports the limiting with local bounds for conservative variables using a two-sided Zalesak-type limiter (Zalesak, 1979) and for general non-linear variables using a one-sided Newton-bisection algorithm. They allow to avoid spurious oscillations within the global bounds and to improve the shock-capturing capabilities of the method. The corresponding numerical admissibility conditions are frequently formulated as local maximum or minimum principles. The local bounds are computed using the maximum and minimum values of all local neighboring nodes. Within this calculation we use the low-order FV solution values for each node.

As for the limiting with global bounds you are passing the quantity names of the conservative variables you want to limit. So, to limit the density with lower and upper local bounds pass the following.

local_twosided_variables_cons = ["rho"]
1-element Vector{String}:
 "rho"

To limit non-linear variables locally, pass the variable function combined with the requested bound (min or max) as a tuple. For instance, to impose a lower local bound on the modified specific entropy Trixi.entropy_guermond_etal, use

local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, min)]
1-element Vector{Tuple{typeof(Trixi.entropy_guermond_etal), typeof(min)}}:
 (Trixi.entropy_guermond_etal, min)

Exemplary simulation

How to set up a simulation using the IDP limiting becomes clearer when looking at an exemplary setup. This will be a simplified version of tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl. Since the setup is mostly very similar to a pure DGSEM setup as in tree_2d_dgsem/elixir_euler_blast_wave.jl, the equivalent parts are used without any explanation here.

using OrdinaryDiffEq
using Trixi

equations = CompressibleEulerEquations2D(1.4)

function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
    # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
    # Set up polar coordinates
    inicenter = SVector(0.0, 0.0)
    x_norm = x[1] - inicenter[1]
    y_norm = x[2] - inicenter[2]
    r = sqrt(x_norm^2 + y_norm^2)
    phi = atan(y_norm, x_norm)
    sin_phi, cos_phi = sincos(phi)

    # Calculate primitive variables
    rho = r > 0.5 ? 1.0 : 1.1691
    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
    p = r > 0.5 ? 1.0E-3 : 1.245

    return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_blast_wave;

Since the surface integral is equal for both the DG and the subcell FV method, the limiting is applied only in the volume integral.

Note, that the DG method is based on the flux differencing formulation. Hence, you have to use a two-point flux, such as flux_ranocha, flux_shima_etal, flux_chandrashekar or flux_kennedy_gruber, for the DG volume flux.

surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
flux_ranocha (generic function with 9 methods)

The limiter is implemented within SubcellLimiterIDP. It always requires the parameters equations and basis. With additional parameters (described above or listed in the docstring) you can specify and enable additional limiting options. Here, the simulation should contain local limiting for the density using lower and upper bounds.

basis = LobattoLegendreBasis(3)
limiter_idp = SubcellLimiterIDP(equations, basis;
                                local_twosided_variables_cons = ["rho"])
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SubcellLimiterIDP                                                                                │
│ ═════════════════                                                                                │
│ Limiter: …………………………………………………………                                                                  │
│ : …………………………………………………………………………… Local two-sided limiting for conservative variables [1]          │
│ Local bounds: …………………………………………… FV solution                                                      │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

The initialized limiter is passed to VolumeIntegralSubcellLimiting in addition to the volume fluxes of the low-order and high-order scheme.

volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                volume_flux_dg = volume_flux,
                                                volume_flux_fv = surface_flux)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ VolumeIntegralSubcellLimiting                                                                    │
│ ═════════════════════════════                                                                    │
│ volume flux DG: ……………………………………… flux_ranocha                                                     │
│ volume flux FV: ……………………………………… FluxLaxFriedrichs(max_abs_speed_naive)                           │
│ limiter: ………………………………………………………… SubcellLimiterIDP                                                │
│ │ Limiter: ……………………………………………………                                                                  │
│ │ : ……………………………………………………………………… Local two-sided limiting for conservative variables [1]          │
│ │ Local bounds: ……………………………………… FV solution                                                      │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

Then, the volume integral is passed to solver as it is done for the standard flux-differencing DG scheme or the element-wise limiting.

solver = DGSEM(basis, surface_flux, volume_integral)
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64}                                                                                      │
│ ═══════════                                                                                      │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3)                         │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3)                      │
│ surface integral: ………………………………… SurfaceIntegralWeakForm                                          │
│ │ surface flux: ……………………………………… FluxLaxFriedrichs(max_abs_speed_naive)                           │
│ volume integral: …………………………………… VolumeIntegralSubcellLimiting                                    │
│ │ volume flux DG: ………………………………… flux_ranocha                                                     │
│ │ volume flux FV: ………………………………… FluxLaxFriedrichs(max_abs_speed_naive)                           │
│ │ limiter: …………………………………………………… SubcellLimiterIDP                                                │
│ │ │ Limiter: ………………………………………………                                                                  │
│ │ │ : ………………………………………………………………… Local two-sided limiting for conservative variables [1]          │
│ │ │ Local bounds: ………………………………… FV solution                                                      │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
coordinates_min = (-2.0, -2.0)
coordinates_max = (2.0, 2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level = 5,
                n_cells_max = 10_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

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

summary_callback = SummaryCallback()

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

alive_callback = AliveCallback(analysis_interval = analysis_interval)

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

stepsize_callback = StepsizeCallback(cfl = 0.3)

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

As explained above, the IDP limiter works a-posteriori and requires the additional use of a correction stage implemented with the stage callback SubcellLimiterIDPCorrection. This callback is passed within a tuple to the time integration method.

stage_callbacks = (SubcellLimiterIDPCorrection(),)
(SubcellLimiterIDPCorrection(),)

Moreover, as mentioned before as well, simulations with subcell limiting require a Trixi-intern SSPRK time integration methods with passed stage callbacks and a Trixi-intern Trixi.solve(...) routine.

sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                  callback = callbacks);
summary_callback() # print the timer summary

████████╗██████╗ ██╗██╗  ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
   ██║   ██████╔╝██║ ╚███╔╝ ██║
   ██║   ██╔══██╗██║ ██╔██╗ ██║
   ██║   ██║  ██║██║██╔╝ ██╗██║
   ╚═╝   ╚═╝  ╚═╝╚═╝╚═╝  ╚═╝╚═╝

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic                                                                     │
│ ════════════════════════════                                                                     │
│ #spatial dimensions: ………………………… 2                                                                │
│ mesh: ………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2}} with length 1365                │
│ equations: …………………………………………………… CompressibleEulerEquations2D                                     │
│ initial condition: ……………………………… initial_condition_blast_wave                                     │
│ boundary conditions: ………………………… Trixi.BoundaryConditionPeriodic                                  │
│ source terms: …………………………………………… nothing                                                          │
│ solver: …………………………………………………………… DG                                                               │
│ total #DOFs per field: …………………… 16384                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{2, Trixi.SerialTree{2}}                                                                 │
│ ════════════════════════════════                                                                 │
│ center: …………………………………………………………… [0.0, 0.0]                                                       │
│ length: …………………………………………………………… 4.0                                                              │
│ periodicity: ……………………………………………… (true, true)                                                     │
│ current #cells: ……………………………………… 1365                                                             │
│ #leaf-cells: ……………………………………………… 1024                                                             │
│ maximum #cells: ……………………………………… 10000                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ CompressibleEulerEquations2D                                                                     │
│ ════════════════════════════                                                                     │
│ #variables: ………………………………………………… 4                                                                │
│ │ variable 1: …………………………………………… rho                                                              │
│ │ variable 2: …………………………………………… rho_v1                                                           │
│ │ variable 3: …………………………………………… rho_v2                                                           │
│ │ variable 4: …………………………………………… rho_e                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64}                                                                                      │
│ ═══════════                                                                                      │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3)                         │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3)                      │
│ surface integral: ………………………………… SurfaceIntegralWeakForm                                          │
│ │ surface flux: ……………………………………… FluxLaxFriedrichs(max_abs_speed_naive)                           │
│ volume integral: …………………………………… VolumeIntegralSubcellLimiting                                    │
│ │ volume flux DG: ………………………………… flux_ranocha                                                     │
│ │ volume flux FV: ………………………………… FluxLaxFriedrichs(max_abs_speed_naive)                           │
│ │ limiter: …………………………………………………… SubcellLimiterIDP                                                │
│ │ │ Limiter: ………………………………………………                                                                  │
│ │ │ : ………………………………………………………………… Local two-sided limiting for conservative variables [1]          │
│ │ │ Local bounds: ………………………………… FV solution                                                      │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AnalysisCallback                                                                                 │
│ ════════════════                                                                                 │
│ interval: ……………………………………………………… 1000                                                             │
│ analyzer: ……………………………………………………… LobattoLegendreAnalyzer{Float64}(polydeg=6)                      │
│ │ error 1: …………………………………………………… l2_error                                                         │
│ │ error 2: …………………………………………………… linf_error                                                       │
│ │ integral 1: …………………………………………… entropy_timederivative                                           │
│ save analysis to file: …………………… no                                                               │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AliveCallback                                                                                    │
│ ═════════════                                                                                    │
│ interval: ……………………………………………………… 100                                                              │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SaveSolutionCallback                                                                             │
│ ════════════════════                                                                             │
│ interval: ……………………………………………………… 1000                                                             │
│ solution variables: …………………………… cons2prim                                                        │
│ save initial solution: …………………… yes                                                              │
│ save final solution: ………………………… yes                                                              │
│ output directory: ………………………………… /home/runner/work/Trixi.jl/Trixi.jl/docs/build/tutorials/out     │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback                                                                                 │
│ ════════════════                                                                                 │
│ CFL number: ………………………………………………… 0.3                                                              │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration                                                                                 │
│ ════════════════                                                                                 │
│ Start time: ………………………………………………… 0.0                                                              │
│ Final time: ………………………………………………… 2.0                                                              │
│ time integrator: …………………………………… SimpleSSPRK33                                                    │
│ adaptive: ……………………………………………………… false                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information                                                                          │
│ ═══════════════════════                                                                          │
│ #threads: ……………………………………………………… 1                                                                │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'CompressibleEulerEquations2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                  0                run time:       7.52000000e-07 s
 Δt:             1.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:       1202.326 MiB
 #elements:                1024

 Variable:       rho              rho_v1           rho_v2           rho_e
 L2 error:       6.25621384e-03   5.88786362e-03   5.81457821e-03   1.15827059e-01
 Linf error:     1.06470791e-01   2.46283676e-01   1.37585923e-01   1.97119199e+00
 ∑∂S/∂U ⋅ Uₜ :  -6.22526192e+02
────────────────────────────────────────────────────────────────────────────────────────────────────

#timesteps:    100 │ Δt: 2.9061e-03 │ sim. time: 2.9394e-01 (14.697%)  │ run time: 1.6542e+00 s
#timesteps:    200 │ Δt: 3.2316e-03 │ sim. time: 6.0126e-01 (30.063%)  │ run time: 2.7737e+00 s
#timesteps:    300 │ Δt: 3.8243e-03 │ sim. time: 9.5591e-01 (47.795%)  │ run time: 3.9175e+00 s
#timesteps:    400 │ Δt: 4.3979e-03 │ sim. time: 1.3674e+00 (68.371%)  │ run time: 5.0887e+00 s
#timesteps:    500 │ Δt: 5.0302e-03 │ sim. time: 1.8395e+00 (91.977%)  │ run time: 6.2695e+00 s

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'CompressibleEulerEquations2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                532                run time:       6.65722482e+00 s
 Δt:             2.58201930e-03                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      2.00000000e+00 (100.000%)     time/DOF/rhs!:  1.92363713e-07 s
                                               PID:            2.54345909e-07 s
 #DOFs per field:         16384                alloc'd memory:       1215.900 MiB
 #elements:                1024

 Variable:       rho              rho_v1           rho_v2           rho_e
 L2 error:       1.01030662e+00   3.05177209e-01   3.05188326e-01   6.85712113e-01
 Linf error:     3.29436672e+00   1.45717455e+00   1.45659687e+00   2.92881060e+00
 ∑∂S/∂U ⋅ Uₜ :  -2.91165098e+01
────────────────────────────────────────────────────────────────────────────────────────────────────

────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished.  Final time: 2.0  Time steps: 532 (accepted), 532 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────

────────────────────────────────────────────────────────────────────────────────
           Trixi.jl                    Time                    Allocations
                              ───────────────────────   ────────────────────────
      Tot / % measured:            6.67s / 100.0%           21.1MiB /  99.5%

Section               ncalls     time    %tot     avg     alloc    %tot      avg
────────────────────────────────────────────────────────────────────────────────
main loop                  1    6.15s   92.2%   6.15s   8.68MiB   41.2%  8.68MiB
  rhs!                 1.60k    5.03s   75.5%  3.15ms   9.33KiB    0.0%    5.98B
    volume integral    1.60k    4.48s   67.2%  2.81ms     0.00B    0.0%    0.00B
    interface flux     1.60k    290ms    4.3%   182μs     0.00B    0.0%    0.00B
    prolong2interf...  1.60k    122ms    1.8%  76.2μs     0.00B    0.0%    0.00B
    surface integral   1.60k    102ms    1.5%  63.8μs     0.00B    0.0%    0.00B
    reset ∂u/∂t        1.60k   18.4ms    0.3%  11.5μs     0.00B    0.0%    0.00B
    Jacobian           1.60k   17.6ms    0.3%  11.0μs     0.00B    0.0%    0.00B
    ~rhs!~             1.60k   2.88ms    0.0%  1.81μs   9.33KiB    0.0%    5.98B
    prolong2mortars    1.60k    131μs    0.0%  82.1ns     0.00B    0.0%    0.00B
    prolong2bounda...  1.60k    111μs    0.0%  69.7ns     0.00B    0.0%    0.00B
    mortar flux        1.60k   78.4μs    0.0%  49.1ns     0.00B    0.0%    0.00B
    boundary flux      1.60k   54.0μs    0.0%  33.8ns     0.00B    0.0%    0.00B
    source terms       1.60k   48.2μs    0.0%  30.2ns     0.00B    0.0%    0.00B
  a posteriori lim...  1.60k    965ms   14.5%   604μs    214KiB    1.0%     137B
    blending factors   1.60k    656ms    9.8%   411μs    213KiB    1.0%     137B
      local twosided   1.60k    587ms    8.8%   368μs    212KiB    1.0%     136B
      ~blending fa...  1.60k   63.9ms    1.0%  40.0μs   1.47KiB    0.0%    0.94B
      reset alpha      1.60k   5.01ms    0.1%  3.14μs     0.00B    0.0%    0.00B
    ~a posteriori ...  1.60k    309ms    4.6%   193μs      752B    0.0%    0.47B
  ~main loop~              1   91.2ms    1.4%  91.2ms   2.31MiB   11.0%  2.31MiB
  calculate dt           532   51.0ms    0.8%  96.0μs     0.00B    0.0%    0.00B
  analyze solution         1   6.66ms    0.1%  6.66ms   4.64MiB   22.0%  4.64MiB
  I/O                      1   2.02ms    0.0%  2.02ms   1.52MiB    7.2%  1.52MiB
    save solution          1   2.01ms    0.0%  2.01ms   1.51MiB    7.2%  1.51MiB
    ~I/O~                  1   5.22μs    0.0%  5.22μs   3.94KiB    0.0%  3.94KiB
    get node varia...      1    862ns    0.0%   862ns     0.00B    0.0%    0.00B
    save mesh              1    221ns    0.0%   221ns     0.00B    0.0%    0.00B
    get element va...      1    211ns    0.0%   211ns     0.00B    0.0%    0.00B
I/O                        2    514ms    7.7%   257ms   7.60MiB   36.1%  3.80MiB
  ~I/O~                    2    512ms    7.7%   256ms   6.09MiB   29.0%  3.05MiB
  save solution            1   2.43ms    0.0%  2.43ms   1.51MiB    7.2%  1.51MiB
  get element vari...      1   1.02μs    0.0%  1.02μs     0.00B    0.0%    0.00B
  get node variables       1    942ns    0.0%   942ns     0.00B    0.0%    0.00B
  save mesh                1   60.0ns    0.0%  60.0ns     0.00B    0.0%    0.00B
analyze solution           1   6.36ms    0.1%  6.36ms   4.76MiB   22.6%  4.76MiB
calculate dt               1   96.6μs    0.0%  96.6μs     0.00B    0.0%    0.00B
────────────────────────────────────────────────────────────────────────────────

Visualization

As for a standard simulation in Trixi.jl, it is possible to visualize the solution using the plot routine from Plots.jl.

using Plots
plot(sol)
Example block output

To get an additional look at the amount of limiting that is used, you can use the visualization approach using the SaveSolutionCallback, Trixi2Vtk and ParaView. More details about this procedure can be found in the visualization documentation.

With that implementation and the standard procedure used for Trixi2Vtk you get the following dropdown menu in ParaView.

ParaView_Dropdownmenu

The resulting visualization of the density and the limiting parameter then looks like this. blast_wave_paraview

You can see that the limiting coefficient does not lie in the interval [0,1] because Trixi2Vtk interpolates all quantities to regular nodes by default. You can disable this functionality with reinterpolate=false within the call of trixi2vtk(...) and get the following visualization. blast_wave_paraview_reinterpolate=false

Bounds checking

Subcell limiting is based on the fulfillment of target bounds - either global or local. Although the implementation works and has been thoroughly tested, there are some cases where these bounds are not met. For instance, the deviations could be in machine precision, which is not problematic. Larger deviations can be cause by too large time-step sizes (which can be easily fixed by reducing the CFL number), specific boundary conditions or source terms. Insufficient parameters for the Newton-bisection algorithm can also be a reason when limiting non-linear variables. There are described above.

In many cases, it is reasonable to monitor the bounds deviations. Because of that, Trixi.jl supports a bounds checking routine implemented using the stage callback BoundsCheckCallback. It checks all target bounds for fulfillment in every RK stage. If added to the tuple of stage callbacks like

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback())

and passed to the time integration method, a summary is added to the final console output. For the given example, this summary shows that all bounds are met at all times.

────────────────────────────────────────────────────────────────────────────────────────────────────
Maximum deviation from bounds:
────────────────────────────────────────────────────────────────────────────────────────────────────
rho:
- lower bound: 0.0
- upper bound: 0.0
────────────────────────────────────────────────────────────────────────────────────────────────────

Moreover, it is also possible to monitor the bounds deviations incurred during the simulations. To do that use the parameter save_errors = true, such that the instant deviations are written to deviations.txt in output_directory every interval time steps.

BoundsCheckCallback(save_errors = true, output_directory = "out", interval = 100)

Then, for the given example the deviations file contains all daviations for the current timestep and simulation time.

iter, simu_time, rho_min, rho_max
100, 0.29103427131404924, 0.0, 0.0
200, 0.5980281923063808, 0.0, 0.0
300, 0.9520853560765293, 0.0, 0.0
400, 1.3630295622683186, 0.0, 0.0
500, 1.8344999624013498, 0.0, 0.0
532, 1.9974179806990118, 0.0, 0.0

This page was generated using Literate.jl.