Add subcell positivity limiting of non-linear variables

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl

@@ -0,0 +1,91 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
+
+A version of the classical Kelvin-Helmholtz instability based on
+- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
+  A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
+  of the Euler Equations
+  [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
+"""
+function initial_condition_kelvin_helmholtz_instability(x, t,
+                                                        equations::CompressibleEulerEquations2D)
+    # change discontinuity to tanh
+    # typical resolution 128^2, 256^2
+    # domain size is [-1,+1]^2
+    slope = 15
+    amplitude = 0.02
+    B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
+    rho = 0.5 + 0.75 * B
+    v1 = 0.5 * (B - 1)
+    v2 = 0.1 * sin(2 * pi * x[1])
+    p = 1.0
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_kelvin_helmholtz_instability
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = ["rho"],
+                                positivity_variables_nonlinear = [pressure])
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.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, 3.7)
+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 = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+save_restart = SaveRestartCallback(interval = 1000,
+                                   save_final_restart = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback,
+                        save_restart, save_solution)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
+
+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

examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl

@@ -22,7 +22,7 @@ function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
     r = sqrt(x[1]^2 + x[2]^2)
 
     pmax = 10.0
-    pmin = 1.0
+    pmin = 0.01
     rhomax = 1.0
     rhomin = 0.01
     if r <= 0.09
@@ -52,7 +52,8 @@ basis = LobattoLegendreBasis(3)
 
 limiter_idp = SubcellLimiterIDP(equations, basis;
                                 positivity_variables_cons = ["rho"],
-                                positivity_correction_factor = 0.5)
+                                positivity_variables_nonlinear = [pressure],
+                                positivity_correction_factor = 0.1)
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)
@@ -84,7 +85,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-cfl = 0.5
+cfl = 0.4
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
 glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)

src/callbacks_stage/subcell_bounds_check.jl

@@ -97,6 +97,9 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
                 end
                 print(f, ", " * string(variables[v]) * "_min")
             end
+            for variable in limiter.positivity_variables_nonlinear
+                print(f, ", " * string(variable) * "_min")
+            end
         end
         println(f)
     end
@@ -140,6 +143,11 @@ end
             println(string(variables[v]) * ":\n- positivity: ",
                     idp_bounds_delta[Symbol(string(v), "_min")][2])
         end
+        for variable in limiter.positivity_variables_nonlinear
+            variable_string = string(variable)
+            println(variable_string * ":\n- positivity: ",
+                    idp_bounds_delta[Symbol(variable_string, "_min")][2])
+        end
     end
     println("─"^100 * "\n")
 

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -48,6 +48,19 @@
             end
             deviation[2] = max(deviation[2], deviation[1])
         end
+        for variable in limiter.positivity_variables_nonlinear
+            key = Symbol(string(variable), "_min")
+            deviation = idp_bounds_delta[key]
+            for element in eachelement(solver, cache), j in eachnode(solver),
+                i in eachnode(solver)
+
+                var = variable(get_node_vars(u, equations, solver, i, j, element),
+                               equations)
+                deviation[1] = max(deviation[1],
+                                   variable_bounds[key][i, j, element] - var)
+            end
+            deviation[2] = max(deviation[2], deviation[1])
+        end
     end
     if save_errors
         # Print to output file
@@ -67,6 +80,10 @@
                     end
                     print(f, ", ", idp_bounds_delta[Symbol(string(v), "_min")][1])
                 end
+                for variable in limiter.positivity_variables_nonlinear
+                    print(f, ", ",
+                          idp_bounds_delta[Symbol(string(variable), "_min")][1])
+                end
             end
             println(f)
         end

src/equations/compressible_euler_2d.jl

@@ -1404,6 +1404,17 @@ end
     return SVector(w1, w2, w3, w4)
 end
 
+# Transformation from conservative variables u to d(p)/d(u)
+@inline function pressure(u, equations::CompressibleEulerEquations2D, derivative::True)
+    rho, rho_v1, rho_v2, rho_e = u
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v_square = v1^2 + v2^2
+
+    return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0)
+end
+
 @inline function entropy2cons(w, equations::CompressibleEulerEquations2D)
     # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD
     # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1)
@@ -1428,6 +1439,14 @@ end
     return SVector(rho, rho_v1, rho_v2, rho_e)
 end
 
+@inline function is_valid_state(cons, equations::CompressibleEulerEquations2D)
+    p = pressure(cons, equations)
+    if cons[1] <= 0.0 || p <= 0.0
+        return false
+    end
+    return true
+end
+
 # Convert primitive to conservative variables
 @inline function prim2cons(prim, equations::CompressibleEulerEquations2D)
     rho, v1, v2, p = prim

src/equations/ideal_glm_mhd_2d.jl

@@ -1062,6 +1062,19 @@ end
     return SVector(w1, w2, w3, w4, w5, w6, w7, w8, w9)
 end
 
+# Transformation from conservative variables u to d(p)/d(u)
+@inline function pressure(u, equations::IdealGlmMhdEquations2D, derivative::True)
+    rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v3 = rho_v3 / rho
+    v_square = v1^2 + v2^2 + v3^2
+
+    return (equations.gamma - 1.0) *
+           SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi)
+end
+
 # Convert entropy variables to conservative variables
 @inline function entropy2cons(w, equations::IdealGlmMhdEquations2D)
     w1, w2, w3, w4, w5, w6, w7, w8, w9 = w
@@ -1089,6 +1102,14 @@ end
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
 end
 
+@inline function is_valid_state(cons, equations::IdealGlmMhdEquations2D)
+    p = pressure(cons, equations)
+    if cons[1] <= 0.0 || p <= 0.0
+        return false
+    end
+    return true
+end
+
 # Convert primitive to conservative variables
 @inline function prim2cons(prim, equations::IdealGlmMhdEquations2D)
     rho, v1, v2, v3, p, B1, B2, B3, psi = prim

src/solvers/dgsem_tree/subcell_limiters.jl

@@ -16,18 +16,28 @@ end
     SubcellLimiterIDP(equations::AbstractEquations, basis;
                       local_minmax_variables_cons = String[],
                       positivity_variables_cons = String[],
-                      positivity_correction_factor = 0.1)
+                      positivity_variables_nonlinear = [],
+                      positivity_correction_factor = 0.1,
+                      max_iterations_newton = 10,
+                      newton_tolerances = (1.0e-12, 1.0e-14),
+                      gamma_constant_newton = 2 * ndims(equations))
 
 Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref)
 including:
 - Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`)
-- Positivity limiting for conservative variables (`positivity_variables_cons`)
+- Positivity limiting for conservative variables (`positivity_variables_cons`) and non-linear variables
+(`positivity_variables_nonlinear`)
 
 Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]`
-and `positivity_variables_cons = ["rho"]`.
+and `positivity_variables_cons = ["rho"]`. For non-linear variables the specific functions are
+passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`.
 
 The bounds are calculated using the low-order FV solution. The positivity limiter uses
 `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`.
+The limiting of non-linear variables uses a Newton-bisection method with a maximum of
+`max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances`
+and a gamma constant of `gamma_constant_newton` (`gamma_constant_newton>=2*d`,
+where `d = #dimensions`).
 
 !!! note
     This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together.
@@ -45,22 +55,32 @@ The bounds are calculated using the low-order FV solution. The positivity limite
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.
 """
-struct SubcellLimiterIDP{RealT <: Real, Cache} <: AbstractSubcellLimiter
+struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <:
+       AbstractSubcellLimiter
     local_minmax::Bool
     local_minmax_variables_cons::Vector{Int}                   # Local mininum/maximum principles for conservative variables
     positivity::Bool
     positivity_variables_cons::Vector{Int}                     # Positivity for conservative variables
+    positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables
     positivity_correction_factor::RealT
     cache::Cache
+    max_iterations_newton::Int
+    newton_tolerances::Tuple{RealT, RealT}  # Relative and absolute tolerances for Newton's method
+    gamma_constant_newton::RealT            # Constant for the subcell limiting of convex (nonlinear) constraints
 end
 
 # this method is used when the limiter is constructed as for shock-capturing volume integrals
 function SubcellLimiterIDP(equations::AbstractEquations, basis;
                            local_minmax_variables_cons = String[],
                            positivity_variables_cons = String[],
-                           positivity_correction_factor = 0.1)
+                           positivity_variables_nonlinear = [],
+                           positivity_correction_factor = 0.1,
+                           max_iterations_newton = 10,
+                           newton_tolerances = (1.0e-12, 1.0e-14),
+                           gamma_constant_newton = 2 * ndims(equations))
     local_minmax = (length(local_minmax_variables_cons) > 0)
-    positivity = (length(positivity_variables_cons) > 0)
+    positivity = (length(positivity_variables_cons) +
+                  length(positivity_variables_nonlinear) > 0)
 
     local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons,
                                                        equations)
@@ -80,13 +100,20 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis;
             bound_keys = (bound_keys..., Symbol(string(v), "_min"))
         end
     end
+    for variable in positivity_variables_nonlinear
+        bound_keys = (bound_keys..., Symbol(string(variable), "_min"))
+    end
 
     cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys)
 
     SubcellLimiterIDP{typeof(positivity_correction_factor),
+                      typeof(positivity_variables_nonlinear),
                       typeof(cache)}(local_minmax, local_minmax_variables_cons_,
                                      positivity, positivity_variables_cons_,
-                                     positivity_correction_factor, cache)
+                                     positivity_variables_nonlinear,
+                                     positivity_correction_factor, cache,
+                                     max_iterations_newton, newton_tolerances,
+                                     gamma_constant_newton)
 end
 
 function Base.show(io::IO, limiter::SubcellLimiterIDP)
@@ -124,7 +151,7 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP)
                 ]
             end
             if positivity
-                string = "positivity for conservative variables $(limiter.positivity_variables_cons)"
+                string = "positivity for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)"
                 setup = [setup..., "" => string]
                 setup = [
                     setup...,