Subcell positivity IDP limiting for conservative variables

examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl

@@ -0,0 +1,91 @@
+
+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)
+indicator_sc = IndicatorIDP(equations, basis;
+                            positivity=true, variables_cons=(Trixi.density,), positivity_correction_factor=0.5)
+volume_integral = VolumeIntegralShockCapturingSubcell(indicator_sc;
+                                                      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 = (AntidiffusiveStage(),)
+
+sol = Trixi.solve(ode; alg=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_sc_subcell.jl

@@ -0,0 +1,143 @@
+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)
+
+density1(u, equations::CompressibleEulerMulticomponentEquations2D) = u[1+3]
+density2(u, equations::CompressibleEulerMulticomponentEquations2D) = u[2+3]
+
+indicator_sc = IndicatorIDP(equations, basis;
+                            positivity=true,
+                            variables_cons=(density1, density2))
+
+volume_integral=VolumeIntegralShockCapturingSubcell(indicator_sc; 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 = (AntidiffusiveStage(),)
+
+sol = Trixi.solve(ode; alg=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

@@ -201,6 +201,7 @@ export DG,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
+       VolumeIntegralShockCapturingSubcell, IndicatorIDP,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
@@ -231,7 +232,7 @@ export ControllerThreeLevel, ControllerThreeLevelCombined,
        IndicatorLöhner, IndicatorLoehner, IndicatorMax,
        IndicatorNeuralNetwork, NeuralNetworkPerssonPeraire, NeuralNetworkRayHesthaven, NeuralNetworkCNN
 
-export PositivityPreservingLimiterZhangShu
+export PositivityPreservingLimiterZhangShu, AntidiffusiveStage
 
 export trixi_include, examples_dir, get_examples, default_example,
        default_example_unstructured, ode_default_options

src/callbacks_stage/antidiffusive_stage.jl

@@ -0,0 +1,73 @@
+# 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
+
+
+"""
+    AntidiffusiveStage()
+
+Perform antidiffusive stage for IDP limiting.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+"""
+struct AntidiffusiveStage end
+
+function (antidiffusive_stage!::AntidiffusiveStage)(u_ode, integrator, stage)
+
+  antidiffusive_stage!(u_ode, integrator.p, integrator.t, integrator.dt, integrator.p.solver.volume_integral)
+end
+
+(::AntidiffusiveStage)(u_ode, semi, t, dt, volume_integral::AbstractVolumeIntegral) = nothing
+
+function (antidiffusive_stage!::AntidiffusiveStage)(u_ode, semi, t, dt, volume_integral::VolumeIntegralShockCapturingSubcell)
+
+  @trixi_timeit timer() "antidiffusive_stage!" antidiffusive_stage!(u_ode, semi, t, dt, volume_integral.indicator)
+end
+
+function (antidiffusive_stage!::AntidiffusiveStage)(u_ode, semi, t, dt, indicator::IndicatorIDP)
+  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+
+  u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+  @trixi_timeit timer() "alpha calculation" semi.solver.volume_integral.indicator(u, semi, solver, t, dt)
+
+  perform_idp_correction(u, dt, mesh, equations, solver, cache)
+
+  return nothing
+end
+
+@inline function perform_idp_correction(u, dt, mesh::TreeMesh2D, equations, dg, cache)
+  @unpack inverse_weights = dg.basis
+  @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.ContainerAntidiffusiveFlux2D
+  @unpack alpha1, alpha2 = dg.volume_integral.indicator.cache.ContainerShockCapturingIndicator
+
+  # Loop over blended DG-FV elements
+  @threaded for element in eachelement(dg, cache)
+    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.0 - alpha1[i,   j, element]) * get_node_vars(antidiffusive_flux1, equations, dg, i,   j, element)
+      alpha_flux1_ip1 = (1.0 - alpha1[i+1, j, element]) * get_node_vars(antidiffusive_flux1, equations, dg, i+1, j, element)
+      alpha_flux2     = (1.0 - alpha2[i,   j, element]) * get_node_vars(antidiffusive_flux2, equations, dg, i,   j, element)
+      alpha_flux2_jp1 = (1.0 - 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
+
+init_callback(callback::AntidiffusiveStage, semi) = nothing
+
+finalize_callback(antidiffusive_stage!::AntidiffusiveStage, semi) = nothing
+
+
+end # @muladd

src/callbacks_stage/callbacks_stage.jl

@@ -6,6 +6,7 @@
 
 
 include("positivity_zhang_shu.jl")
+include("antidiffusive_stage.jl")
 
 
 end # @muladd

src/solvers/dg.jl

@@ -172,6 +172,47 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralPureLGLFi
 end
 
 
+"""
+    VolumeIntegralShockCapturingSubcell(indicator;
+                                        volume_flux_dg, volume_flux_fv)
+
+A shock-capturing volume integral type for DG methods based on a subcell blending approach
+with a low-order FV method from the preprint paper
+- Rueda-Ramírez, Pazner, Gassner (2022)
+  "Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods"
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+
+See also: [`VolumeIntegralShockCapturingHG`](@ref)
+"""
+struct VolumeIntegralShockCapturingSubcell{VolumeFluxDG, VolumeFluxFV, Indicator} <: AbstractVolumeIntegral
+  volume_flux_dg::VolumeFluxDG
+  volume_flux_fv::VolumeFluxFV
+  indicator::Indicator
+end
+
+function VolumeIntegralShockCapturingSubcell(indicator; volume_flux_dg,
+                                                        volume_flux_fv)
+  VolumeIntegralShockCapturingSubcell{typeof(volume_flux_dg), typeof(volume_flux_fv), typeof(indicator)}(
+    volume_flux_dg, volume_flux_fv, indicator)
+end
+
+function Base.show(io::IO, mime::MIME"text/plain", integral::VolumeIntegralShockCapturingSubcell)
+  @nospecialize integral # reduce precompilation time
+
+  if get(io, :compact, false)
+    show(io, integral)
+  else
+    summary_header(io, "VolumeIntegralShockCapturingSubcell")
+    summary_line(io, "volume flux DG", integral.volume_flux_dg)
+    summary_line(io, "volume flux FV", integral.volume_flux_fv)
+    summary_line(io, "indicator", integral.indicator |> typeof |> nameof)
+    show(increment_indent(io), mime, integral.indicator)
+    summary_footer(io)
+  end
+end
+
 # TODO: FD. Should this definition live in a different file because it is
 # not strictly a DG method?
 """