Implement IDP mortars with limiting

examples/tree_2d_dgsem/elixir_euler_density_wave_nonconforming_idp_mortars.jl

@@ -0,0 +1,76 @@
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_density_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 5
+basis = LobattoLegendreBasis(polydeg)
+
+# No limiting is needed in the volume integral
+# It is still required to use `VolumeIntegralSubcellLimiting` in order to use `MortarIDP`.
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = [],
+                                positivity_variables_nonlinear = [])
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+mortar = MortarIDP(equations, basis;
+                   positivity_variables_cons = ["rho"])
+solver = DGSEM(basis, surface_flux, volume_integral, mortar)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+refinement_patches = ((type = "box", coordinates_min = (-0.5, -0.5),
+                       coordinates_max = (0.0, 1.0)),
+                      (type = "box", coordinates_min = (0.5, -1.0),
+                       coordinates_max = (1.0, 0.5)))
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                refinement_patches = refinement_patches,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+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 = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = (:limiting_coefficient,))
+
+stepsize_callback = StepsizeCallback(cfl = 0.4)
+
+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
+                  callback = callbacks);

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl → examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_sc_subcell.jl

@@ -30,25 +30,20 @@ function initial_condition_kelvin_helmholtz_instability(x, t,
 end
 initial_condition = initial_condition_kelvin_helmholtz_instability
 
-# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
-# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
-# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
-# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
-# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
-# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
-# `StepsizeCallback` (CFL-Condition) and less diffusion.
-surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
+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])
+                                positivity_variables_nonlinear = [pressure],
+                                local_twosided_variables_cons = []) # required for testing
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)
-solver = DGSEM(basis, surface_flux, volume_integral)
+mortar = MortarIDP(equations, basis; positivity_variables_cons = ["rho"])
+solver = DGSEM(basis, surface_flux, volume_integral, mortar)
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)
@@ -79,12 +74,27 @@ save_solution = SaveSolutionCallback(interval = 100,
 save_restart = SaveRestartCallback(interval = 1000,
                                    save_final_restart = true)
 
+amr_indicator = IndicatorHennemannGassner(semi,
+                                          alpha_max = 1.0,
+                                          alpha_min = 0.0001,
+                                          alpha_smooth = false,
+                                          variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 4,
+                                      med_level = 0, med_threshold = 0.0003,
+                                      max_level = 6, max_threshold = 0.003)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 1,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
 stepsize_callback = StepsizeCallback(cfl = 0.7)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        stepsize_callback,
-                        save_restart, save_solution)
+                        amr_callback,
+                        save_restart, save_solution,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
@@ -95,5 +105,4 @@ 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
-                  ode_default_options()...,
                   callback = callbacks);

src/Trixi.jl

@@ -265,7 +265,7 @@ export DG,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
-       MortarL2
+       MortarL2, MortarIDP
 
 export reconstruction_O2_inner, reconstruction_O2_full,
        reconstruction_constant,

src/callbacks_stage/subcell_limiter_idp_correction.jl

@@ -56,6 +56,19 @@ function (limiter!::SubcellLimiterIDPCorrection)(u_ode, semi, t, dt,
 
     perform_idp_correction!(u, dt, mesh, equations, solver, cache)
 
+    if solver.mortar isa Trixi.LobattoLegendreMortarIDP
+        @trixi_timeit timer() "mortar blending factors" calc_mortar_limiting_factor!(u,
+                                                                                     semi,
+                                                                                     t,
+                                                                                     dt)
+
+        @trixi_timeit timer() "mortar correction" perform_idp_mortar_correction(u, dt,
+                                                                                mesh,
+                                                                                equations,
+                                                                                solver,
+                                                                                cache)
+    end
+
     return nothing
 end
 

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -67,4 +67,112 @@ function perform_idp_correction!(u, dt,
 
     return nothing
 end
+
+function perform_idp_mortar_correction(u, dt, mesh::TreeMesh{2}, equations, dg, cache)
+    (; orientations, limiting_factor) = cache.mortars
+
+    (; surface_flux_values) = cache.elements
+    (; surface_flux_values_high_order) = cache.antidiffusive_fluxes
+    (; boundary_interpolation) = dg.basis
+
+    for mortar in eachmortar(dg, cache)
+        if isapprox(limiting_factor[mortar], one(eltype(limiting_factor)))
+            continue
+        end
+
+        if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+            if orientations[mortar] == 1
+                direction_small = 1
+                direction_large = 2
+            else
+                direction_small = 3
+                direction_large = 4
+            end
+            # In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
+            # This sign switch is directly applied to the boundary interpolation factors here.
+            factor_small = boundary_interpolation[1, 1]
+            factor_large = -boundary_interpolation[nnodes(dg), 2]
+        else # large_sides[mortar] == 2 -> small elements on left side
+            if orientations[mortar] == 1
+                direction_small = 2
+                direction_large = 1
+            else
+                direction_small = 4
+                direction_large = 3
+            end
+            # In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
+            # This sign switch is directly applied to the boundary interpolation factors here.
+            factor_large = boundary_interpolation[1, 1]
+            factor_small = -boundary_interpolation[nnodes(dg), 2]
+        end
+
+        for i in eachnode(dg)
+            if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (1, i)
+                    indices_large = (nnodes(dg), i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, 1)
+                    indices_large = (i, nnodes(dg))
+                end
+            else # large_sides[mortar] == 2 -> small elements on left side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (nnodes(dg), i)
+                    indices_large = (1, i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, nnodes(dg))
+                    indices_large = (i, 1)
+                end
+            end
+
+            # small elements
+            for small_element_index in 1:2
+                small_element = cache.mortars.neighbor_ids[small_element_index, mortar]
+                inverse_jacobian_small = get_inverse_jacobian(cache.elements.inverse_jacobian,
+                                                              mesh, indices_small...,
+                                                              small_element)
+
+                flux_small_high_order = get_node_vars(surface_flux_values_high_order,
+                                                      equations, dg,
+                                                      i, direction_small, small_element)
+                flux_small_low_order = get_node_vars(surface_flux_values, equations, dg,
+                                                     i, direction_small, small_element)
+                flux_difference_small = factor_small *
+                                        (flux_small_high_order .- flux_small_low_order)
+
+                multiply_add_to_node_vars!(u,
+                                           dt * inverse_jacobian_small *
+                                           (1 - limiting_factor[mortar]),
+                                           flux_difference_small, equations, dg,
+                                           indices_small..., small_element)
+            end
+
+            # large element
+            large_element = cache.mortars.neighbor_ids[3, mortar]
+            inverse_jacobian_large = get_inverse_jacobian(cache.elements.inverse_jacobian,
+                                                          mesh, indices_large...,
+                                                          large_element)
+
+            flux_large_high_order = get_node_vars(surface_flux_values_high_order,
+                                                  equations, dg, i, direction_large,
+                                                  large_element)
+            flux_large_low_order = get_node_vars(surface_flux_values, equations, dg, i,
+                                                 direction_large, large_element)
+            flux_difference_large = factor_large *
+                                    (flux_large_high_order .- flux_large_low_order)
+
+            multiply_add_to_node_vars!(u,
+                                       dt * inverse_jacobian_large *
+                                       (1 - limiting_factor[mortar]),
+                                       flux_difference_large, equations, dg,
+                                       indices_large..., large_element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd