Subcell positivity IDP limiting for conservative variables

examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl

@@ -0,0 +1,92 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+"""
+    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+
+A medium blast wave (modified to lower density and higher pressure) taken from
+- Sebastian Hennemann, Gregor J. Gassner (2020)
+  A provably entropy stable subcell shock capturing approach for high order split form DG
+  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
+"""
+function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+  # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> modified to lower density, higher pressure
+  # 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         "normal" medium blast wave
+  rho = r > 0.5 ? 0.1 : 0.2691            # rho = r > 0.5 ? 1 : 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-1 : 1.245          # p   = r > 0.5 ? 1.0E-3 : 1.245
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux  = flux_ranocha
+basis = LobattoLegendreBasis(3)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons=[1],
+                                positivity_correction_factor=0.5)
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg=volume_flux,
+                                                volume_flux_fv=surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+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=100_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=0.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(),)
+
+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
+                  save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl

@@ -0,0 +1,140 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler multicomponent equations
+
+# 1) Dry Air  2) Helium + 28% Air
+equations = CompressibleEulerMulticomponentEquations2D(gammas        = (1.4, 1.648),
+                                                       gas_constants = (0.287, 1.578))
+
+"""
+    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})
+
+A shock-bubble testcase for multicomponent Euler equations
+- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
+  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
+  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
+"""
+function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})
+  # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
+  # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
+  # typical domain is rectangular, we change it to a square, as Trixi can only do squares
+  @unpack gas_constants = equations
+
+  # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
+  delta   = 0.03
+
+  # Region I
+  rho1_1  = delta
+  rho2_1  = 1.225 * gas_constants[1]/gas_constants[2] - delta
+  v1_1    = zero(delta)
+  v2_1    = zero(delta)
+  p_1     = 101325
+
+  # Region II
+  rho1_2  = 1.225-delta
+  rho2_2  = delta
+  v1_2    = zero(delta)
+  v2_2    = zero(delta)
+  p_2     = 101325
+
+  # Region III
+  rho1_3  = 1.6861 - delta
+  rho2_3  = delta
+  v1_3    = -113.5243
+  v2_3    = zero(delta)
+  p_3     = 159060
+
+  # Set up Region I & II:
+  inicenter = SVector(zero(delta), zero(delta))
+  x_norm = x[1] - inicenter[1]
+  y_norm = x[2] - inicenter[2]
+  r = sqrt(x_norm^2 + y_norm^2)
+
+  if (x[1] > 0.50)
+    # Set up Region III
+    rho1    = rho1_3
+    rho2    = rho2_3
+    v1      = v1_3
+    v2      = v2_3
+    p       = p_3
+  elseif (r < 0.25)
+    # Set up Region I
+    rho1    = rho1_1
+    rho2    = rho2_1
+    v1      = v1_1
+    v2      = v2_1
+    p       = p_1
+  else
+    # Set up Region II
+    rho1    = rho1_2
+    rho2    = rho2_2
+    v1      = v1_2
+    v2      = v2_2
+    p       = p_2
+  end
+
+  return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
+end
+initial_condition = initial_condition_shock_bubble
+
+surface_flux        = flux_lax_friedrichs
+volume_flux         = flux_ranocha
+basis               = LobattoLegendreBasis(3)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons=[(i+3 for i in eachcomponent(equations))...])
+
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg=volume_flux,
+                                                volume_flux_fv=surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min     = (-2.25, -2.225)
+coordinates_max     = ( 2.20,  2.225)
+mesh                = TreeMesh(coordinates_min, coordinates_max,
+                               initial_refinement_level=3,
+                               n_cells_max=1_000_000)
+
+semi                = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan               = (0.0, 0.01)
+ode                 = semidiscretize(semi, tspan)
+
+summary_callback    = SummaryCallback()
+
+analysis_interval   = 300
+analysis_callback   = AnalysisCallback(semi, interval=analysis_interval,
+                                       extra_analysis_integrals=(Trixi.density,))
+
+alive_callback      = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution       = SaveSolutionCallback(interval=300,
+                                           save_initial_solution=true,
+                                           save_final_solution=true,
+                                           solution_variables=cons2prim)
+
+stepsize_callback   = StepsizeCallback(cfl=0.9)
+
+callbacks           = CallbackSet(summary_callback,
+                                  analysis_callback,
+                                  alive_callback,
+                                  save_solution,
+                                  stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(),)
+
+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
+                  save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -119,10 +119,10 @@ include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
 include("semidiscretization/semidiscretization_coupled.jl")
+include("time_integration/time_integration.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
-include("time_integration/time_integration.jl")
 
 # `trixi_include` and special elixirs such as `convergence_test`
 include("auxiliary/special_elixirs.jl")
@@ -215,6 +215,7 @@ export DG,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
+       VolumeIntegralSubcellLimiting, SubcellLimiterIDP,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
@@ -248,7 +249,7 @@ export ControllerThreeLevel, ControllerThreeLevelCombined,
        IndicatorNeuralNetwork, NeuralNetworkPerssonPeraire, NeuralNetworkRayHesthaven,
        NeuralNetworkCNN
 
-export PositivityPreservingLimiterZhangShu
+export PositivityPreservingLimiterZhangShu, SubcellLimiterIDPCorrection
 
 export trixi_include, examples_dir, get_examples, default_example,
        default_example_unstructured, ode_default_options

src/callbacks_stage/callbacks_stage.jl

@@ -6,4 +6,5 @@
 #! format: noindent
 
 include("positivity_zhang_shu.jl")
+include("subcell_limiter_idp_correction.jl")
 end # @muladd

src/callbacks_stage/subcell_limiter_idp_correction.jl

@@ -0,0 +1,68 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+"""
+    SubcellLimiterIDPCorrection()
+
+Perform antidiffusive correction stage for the a posteriori IDP limiter [`SubcellLimiterIDP`](@ref)
+called with [`VolumeIntegralSubcellLimiting`](@ref).
+
+!!! note
+    This callback and the actual limiter [`SubcellLimiterIDP`](@ref) only work together.
+    This is not a replacement but a necessary addition.
+
+## References
+
+- Rueda-Ramírez, Pazner, Gassner (2022)
+  Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods
+  [DOI: 10.1016/j.compfluid.2022.105627](https://doi.org/10.1016/j.compfluid.2022.105627)
+- Pazner (2020)
+  Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting
+  [DOI: 10.1016/j.cma.2021.113876](https://doi.org/10.1016/j.cma.2021.113876)
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+"""
+struct SubcellLimiterIDPCorrection end
+
+function (limiter!::SubcellLimiterIDPCorrection)(u_ode,
+                                                 integrator::Trixi.SimpleIntegratorSSP,
+                                                 stage)
+    limiter!(u_ode, integrator.p, integrator.t, integrator.dt,
+             integrator.p.solver.volume_integral)
+end
+
+function (limiter!::SubcellLimiterIDPCorrection)(u_ode, semi, t, dt,
+                                                 volume_integral::VolumeIntegralSubcellLimiting)
+    @trixi_timeit timer() "a posteriori limiter" limiter!(u_ode, semi, t, dt,
+                                                          volume_integral.limiter)
+end
+
+function (limiter!::SubcellLimiterIDPCorrection)(u_ode, semi, t, dt,
+                                                 limiter::SubcellLimiterIDP)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+    # Calculate blending factor alpha in [0,1]
+    # f_ij = alpha_ij * f^(FV)_ij + (1 - alpha_ij) * f^(DG)_ij
+    #      = f^(FV)_ij + (1 - alpha_ij) * f^(antidiffusive)_ij
+    @trixi_timeit timer() "blending factors" solver.volume_integral.limiter(u, semi,
+                                                                            solver, t,
+                                                                            dt)
+
+    perform_idp_correction!(u, dt, mesh, equations, solver, cache)
+
+    return nothing
+end
+
+init_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing
+
+finalize_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing
+
+include("subcell_limiter_idp_correction_2d.jl")
+end # @muladd

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -0,0 +1,44 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@inline function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache)
+    @unpack inverse_weights = dg.basis
+    @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.container_antidiffusive_flux
+    @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.container_subcell_limiter
+
+    @threaded for element in eachelement(dg, cache)
+        # Sign switch as in apply_jacobian!
+        inverse_jacobian = -cache.elements.inverse_jacobian[element]
+
+        for j in eachnode(dg), i in eachnode(dg)
+            # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1}
+            alpha_flux1 = (1 - alpha1[i, j, element]) *
+                          get_node_vars(antidiffusive_flux1, equations, dg, i, j,
+                                        element)
+            alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) *
+                              get_node_vars(antidiffusive_flux1, equations, dg, i + 1,
+                                            j, element)
+            alpha_flux2 = (1 - alpha2[i, j, element]) *
+                          get_node_vars(antidiffusive_flux2, equations, dg, i, j,
+                                        element)
+            alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) *
+                              get_node_vars(antidiffusive_flux2, equations, dg, i,
+                                            j + 1, element)
+
+            for v in eachvariable(equations)
+                u[v, i, j, element] += dt * inverse_jacobian *
+                                       (inverse_weights[i] *
+                                        (alpha_flux1_ip1[v] - alpha_flux1[v]) +
+                                        inverse_weights[j] *
+                                        (alpha_flux2_jp1[v] - alpha_flux2[v]))
+            end
+        end
+    end
+
+    return nothing
+end
+end # @muladd