Add support for P4estMesh to subcell-limiting (#121)

* Add support for P4estMesh to `subcell-limiting` * Add elixir with supersonic flow * Extract inner code of `characteristic_boundary_value_function` * Add support for `BoundaryConditionDirichlet` for P4estMesh * Fix calculation of `vn` * Remove default routine for `get_boundary_outer_state` * Revise calculation of bar states * Fix allocations with `foreach(enumerate(...))` * Fix normalization in `get_boundary_outer_state` * Add TODO note for future * Remove unnecessary code; format * Adapt test errors after normalizing normal vectors * Add dispatch for equations * Fix normalization; Is done inside boundary_value_function * Fix bug in calculation of bar state at interfaces for P4estMesh * Add mesh as parameter of `get_boundary_outer_state` - Fixes bug that P4estMesh calls the right `characteristic_boundary_state` routine * Add p4est double mach elixir * Fix error * Remove comment * Add test for supersonic flow * fmt * Add comments * Update src/equations/compressible_euler_2d.jl Co-authored-by: Andrés Rueda-Ramírez <[email protected]> * fmt * Fix error --------- Co-authored-by: Andrés Rueda-Ramírez <[email protected]>
trixi-framework · Dec 5, 2023 · 1a09eaa · 1a09eaa
1 parent 647ead5
commit 1a09eaa
Show file tree

Hide file tree

Showing 18 changed files with 1,066 additions and 30 deletions.
diff --git a/examples/p4est_2d_dgsem/elixir_euler_double_mach.jl b/examples/p4est_2d_dgsem/elixir_euler_double_mach.jl
@@ -0,0 +1,167 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_double_mach_reflection(x, t, equations::CompressibleEulerEquations2D)
+
+Compressible Euler setup for a double Mach reflection problem.
+Involves strong shock interactions as well as steady / unsteady flow structures.
+Also exercises special boundary conditions along the bottom of the domain that is a mixture of
+Dirichlet and slip wall.
+See Section IV c on the paper below for details.
+
+- Paul Woodward and Phillip Colella (1984)
+  The Numerical Simulation of Two-Dimensional Fluid Flows with Strong Shocks.
+  [DOI: 10.1016/0021-9991(84)90142-6](https://doi.org/10.1016/0021-9991(84)90142-6)
+"""
+@inline function initial_condition_double_mach_reflection(x, t,
+                                                          equations::CompressibleEulerEquations2D)
+    if x[1] < 1 / 6 + (x[2] + 20 * t) / sqrt(3)
+        phi = pi / 6
+        sin_phi, cos_phi = sincos(phi)
+
+        rho = 8.0
+        v1 = 8.25 * cos_phi
+        v2 = -8.25 * sin_phi
+        p = 116.5
+    else
+        rho = 1.4
+        v1 = 0.0
+        v2 = 0.0
+        p = 1.0
+    end
+
+    prim = SVector(rho, v1, v2, p)
+    return prim2cons(prim, equations)
+end
+
+initial_condition = initial_condition_double_mach_reflection
+
+boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition)
+
+# Special mixed boundary condition type for the :y_neg of the domain.
+# It is charachteristic-based when x < 1/6 and a slip wall when x >= 1/6
+# Note: Only for P4estMesh
+@inline function boundary_condition_mixed_characteristic_wall(u_inner,
+                                                              normal_direction::AbstractVector,
+                                                              x, t, surface_flux_function,
+                                                              equations::CompressibleEulerEquations2D)
+    if x[1] < 1 / 6
+        # From the BoundaryConditionCharacteristic
+        # get the external state of the solution
+        u_boundary = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
+                                                                  u_inner,
+                                                                  normal_direction,
+                                                                  x, t,
+                                                                  equations)
+        # Calculate boundary flux
+        flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+    else # x[1] >= 1 / 6
+        # Use the free slip wall BC otherwise
+        flux = boundary_condition_slip_wall(u_inner, normal_direction, x, t,
+                                            surface_flux_function, equations)
+    end
+
+    return flux
+end
+
+# Note: Only for P4estMesh
+@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
+                                                boundary_condition::typeof(boundary_condition_mixed_characteristic_wall),
+                                                normal_direction::AbstractVector, direction,
+                                                mesh::P4estMesh{2},
+                                                equations::CompressibleEulerEquations2D,
+                                                dg, indices...)
+    x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)
+    if x[1] < 1 / 6 # BoundaryConditionCharacteristic
+        u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
+                                                               u_inner,
+                                                               normal_direction,
+                                                               x, t, equations)
+
+    else # if x[1] >= 1 / 6 # boundary_condition_slip_wall
+        factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3])
+        u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction
+
+        u_outer = SVector(u_inner[1],
+                          u_inner[2] - 2.0 * u_normal[1],
+                          u_inner[3] - 2.0 * u_normal[2],
+                          u_inner[4])
+    end
+    return u_outer
+end
+
+boundary_conditions = Dict(:y_neg => boundary_condition_mixed_characteristic_wall,
+                           :y_pos => boundary_condition_inflow_outflow,
+                           :x_pos => boundary_condition_inflow_outflow,
+                           :x_neg => boundary_condition_inflow_outflow)
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                local_minmax_variables_cons = ["rho"],
+                                spec_entropy = true,
+                                positivity_correction_factor = 0.1,
+                                max_iterations_newton = 100,
+                                bar_states = true)
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+initial_refinement_level = 4
+trees_per_dimension = (4 * 2^initial_refinement_level, 2^initial_refinement_level)
+coordinates_min = (0.0, 0.0)
+coordinates_max = (4.0, 1.0)
+mesh = P4estMesh(trees_per_dimension, polydeg = polydeg,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 initial_refinement_level = 0,
+                 periodicity = false)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     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,
+                        stepsize_callback,
+                        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
diff --git a/examples/p4est_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl
@@ -0,0 +1,91 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_constant
+
+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_correction_factor = 0.1,
+                                spec_entropy = false,
+                                bar_states = true)
+
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D.
+# This particular mesh is unstructured in the yz-plane, but extruded in x-direction.
+# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded
+# in x-direction to ensure free stream preservation on a non-conforming mesh.
+# See https://doi.org/10.1007/s10915-018-00897-9, Section 6.
+
+# Mapping as described in https://arxiv.org/abs/2012.12040, but reduced to 2D
+function mapping(xi_, eta_)
+    # Transform input variables between -1 and 1 onto [0,3]
+    xi = 1.5 * xi_ + 1.5
+    eta = 1.5 * eta_ + 1.5
+
+    y = eta + 3 / 8 * (cos(1.5 * pi * (2 * xi - 3) / 3) *
+                       cos(0.5 * pi * (2 * eta - 3) / 3))
+
+    x = xi + 3 / 8 * (cos(0.5 * pi * (2 * xi - 3) / 3) *
+                      cos(2 * pi * (2 * y - 3) / 3))
+
+    return SVector(x, y)
+end
+
+trees_per_dimension = (16, 16)
+
+# Create P4estMesh with 16 x 16 trees and 16 x 16 elements
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 mapping = mapping,
+                 initial_refinement_level = 0, periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.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 = 10000,
+                                     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,
+                        stepsize_callback,
+                        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
+                  save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary