Unify implementation VolumeIntegralShockCapturingHG and VolumeIntegralShockCapturingRRG

NEWS.md

@@ -14,6 +14,12 @@ This is useful for (locally) diffusion-dominated problems.
 This enables in particular adaptive mesh refinement for that solver-mesh combination ([#2712]).
 - Added functionality to `ScalarPlotData2D` allowing visualization a field provided by a user-defined scalar function ([#2796]).
 - Added `NonIdealCompressibleEuler2D` ([#2768]).
+- Generalization of `VolumeIntegralShockCapturingHG` and `VolumeIntegralShockCapturingRRG` to support different volume integrals on the 
+  non-stabilized and stabilized elements/cells.
+  The generalized volume integral is called `VolumeIntegralShockCapturingHGType` and takes the three keyword arguments `volume_integral_default`,
+  `volume_integral_blend_high_order`, and `volume_integral_blend_low_order` besides the usual `indicator` argument.
+  In particular, `volume_integral_default` may be e.g.  `VolumeIntegralWeakForm` or `VolumeIntegralAdaptive`, i.e.,
+  the non-stabilized elements/cells are no longer restricted to `VolumeIntegralFluxDifferencing` only ([#2802]).
 
 ## Changes when updating to v0.15 from v0.14.x
 

examples/p4est_2d_dgsem/elixir_navierstokes_RAE2822airfoil_separation.jl

@@ -0,0 +1,138 @@
+using OrdinaryDiffEqSSPRK
+using OrdinaryDiffEqCore: PIDController
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+gamma() = 1.4
+equations = CompressibleEulerEquations2D(gamma())
+
+Re() = 6.5 * 10^6
+airfoil_chord_length() = 1.0
+
+# See https://www1.grc.nasa.gov/wp-content/uploads/case_c2.1.pdf or
+# https://cfd.ku.edu/hiocfd/case_c2.2.html
+U_inf() = 0.734 # Mach_inf = 1.0
+rho_inf() = gamma() # => p_inf = 1.0
+
+mu() = rho_inf() * U_inf() * airfoil_chord_length() / Re()
+prandtl_number() = 0.71
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+                                                          Prandtl = prandtl_number())
+
+aoa() = deg2rad(2.79) # 2.79 Degree angle of attack
+
+@inline function initial_condition_mach085_flow(x, t, equations)
+    v1 = 0.73312995164809
+    v2 = 0.03572777625978245
+
+    prim = SVector(1.4, v1, v2, 1.0)
+    return prim2cons(prim, equations)
+end
+initial_condition = initial_condition_mach085_flow
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+shock_indicator = IndicatorHennemannGassner(equations, basis,
+                                            alpha_max = 1.0,
+                                            alpha_min = 0.001,
+                                            alpha_smooth = true,
+                                            variable = density_pressure)
+
+# In non-blended/limited regions, we use the cheaper weak form volume integral
+volume_integral_default = VolumeIntegralWeakForm()
+
+surface_flux = flux_hllc
+volume_flux = flux_ranocha
+# For the blended/limited regions, we need to supply high-order and low-order volume integrals.
+volume_integral_blend_high_order = VolumeIntegralFluxDifferencing(volume_flux)
+volume_integral_blend_low_order = VolumeIntegralPureLGLFiniteVolumeO2(basis;
+                                                                      volume_flux_fv = surface_flux,
+                                                                      reconstruction_mode = reconstruction_O2_inner,
+                                                                      slope_limiter = minmod)
+
+volume_integral = VolumeIntegralShockCapturingHGType(shock_indicator;
+                                                     volume_integral_default = volume_integral_default,
+                                                     volume_integral_blend_high_order = volume_integral_blend_high_order,
+                                                     volume_integral_blend_low_order = volume_integral_blend_low_order)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+###############################################################################
+
+# mesh downloaded from https://cfd.ku.edu/hiocfd/rae2822/
+mesh_file = Trixi.download("https://gist.githubusercontent.com/DanielDoehring/373727b8bc43e4aaeb63a6fcea77f098/raw/99cfd7c6b35df1a28d11db71be4b7702522cc84f/rae2822_level3.inp",
+                           joinpath(@__DIR__, "rae2822_level3.inp"))
+
+boundary_symbols = [:WallBoundary, :FarfieldBoundary]
+mesh = P4estMesh{2}(mesh_file, boundary_symbols = boundary_symbols)
+
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+
+velocity_bc_airfoil = NoSlip((x, t, equations) -> SVector(0.0, 0.0))
+heat_bc = Adiabatic((x, t, equations) -> 0.0)
+boundary_condition_airfoil = BoundaryConditionNavierStokesWall(velocity_bc_airfoil, heat_bc)
+
+boundary_conditions_hyp = (; FarfieldBoundary = boundary_condition_free_stream,
+                           WallBoundary = boundary_condition_slip_wall)
+
+boundary_conditions_para = (; FarfieldBoundary = boundary_condition_free_stream,
+                            WallBoundary = boundary_condition_airfoil)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions_hyp,
+                                                                    boundary_conditions_para))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+t_c = airfoil_chord_length() / U_inf() # convective time
+tspan = (0.0, 25 * t_c)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+save_sol_interval = 50_000
+save_solution = SaveSolutionCallback(interval = save_sol_interval,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+force_boundary_names = (:WallBoundary,)
+drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
+                                           DragCoefficientPressure2D(aoa(), rho_inf(),
+                                                                     U_inf(),
+                                                                     airfoil_chord_length()))
+lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
+                                           LiftCoefficientPressure2D(aoa(), rho_inf(),
+                                                                     U_inf(),
+                                                                     airfoil_chord_length()))
+
+analysis_callback = AnalysisCallback(semi, interval = save_sol_interval,
+                                     output_directory = "out",
+                                     analysis_errors = Symbol[],
+                                     analysis_integrals = (drag_coefficient,
+                                                           lift_coefficient))
+
+alive_callback = AliveCallback(alive_interval = 500)
+
+save_restart = SaveRestartCallback(interval = save_sol_interval,
+                                   save_final_restart = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution, save_restart)
+
+###############################################################################
+# run the simulation
+
+ode_algorithm = SSPRK43(thread = Trixi.True())
+
+time_int_tol = 1e-4
+sol = solve(ode, ode_algorithm;
+            abstol = time_int_tol, reltol = time_int_tol, dt = 1e-6,
+            maxiters = Inf, # long simulation
+            controller = PIDController(0.55, -0.27, 0.05), # optimized for SSPRK43
+            ode_default_options()..., callback = callbacks)

examples/tree_3d_dgsem/elixir_euler_sedov_blast_weak_form_sc.jl

@@ -0,0 +1,122 @@
+using OrdinaryDiffEqSSPRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+"""
+    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+
+The Sedov blast wave setup based on example 35.1.4 from Flash
+- https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
+"""
+function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D)
+    # Set up polar coordinates
+    RealT = eltype(x)
+    inicenter = SVector(0, 0, 0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
+
+    # Setup based on example 35.1.4 in https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
+    r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    E = 1
+    nu = 3 # dims
+    p0_inner = 3 * (equations.gamma - 1) * E / ((nu + 1) * convert(RealT, pi) * r0^nu)
+    p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup
+
+    # Calculate primitive variables
+    rho = 1
+    v1 = 0
+    v2 = 0
+    v3 = 0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+initial_condition = initial_condition_sedov_blast_wave
+
+basis = LobattoLegendreBasis(3)
+
+# Use standard Hennemann-Gassner a-priori shock & blending indicator
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+
+# In non-blended/limited regions, we use the cheaper weak form volume integral
+volume_integral_default = VolumeIntegralWeakForm()
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_chandrashekar
+# For the blended/limited regions, we need to supply high-order and low-order volume integrals.
+volume_integral_blend_high_order = VolumeIntegralFluxDifferencing(volume_flux)
+volume_integral_blend_low_order = VolumeIntegralPureLGLFiniteVolume(surface_flux)
+
+volume_integral = VolumeIntegralShockCapturingHGType(indicator_sc;
+                                                     volume_integral_default = volume_integral_default,
+                                                     volume_integral_blend_high_order = volume_integral_blend_high_order,
+                                                     volume_integral_blend_low_order = volume_integral_blend_low_order)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-2.0, -2.0, -2.0)
+coordinates_max = (2.0, 2.0, 2.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 5,
+                n_cells_max = 1_000_000,
+                periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(alive_interval = 20)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_indicator = IndicatorHennemannGassner(semi,
+                                          alpha_max = 1.0,
+                                          alpha_min = 0.0,
+                                          alpha_smooth = false,
+                                          variable = density_pressure)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 2,
+                                      max_level = 6, max_threshold = 0.0003)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 2,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = false)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        amr_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK54(thread = Trixi.True());
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -269,6 +269,7 @@ export DG,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, VolumeIntegralShockCapturingRRG,
+       VolumeIntegralShockCapturingHGType,
        VolumeIntegralAdaptive, IndicatorEntropyChange,
        IndicatorHennemannGassner,
        VolumeIntegralUpwind,