From 2ac203e7c52182e8cbd7841d58ecb1e285c96322 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Fri, 30 Aug 2024 12:19:16 +0200 Subject: [PATCH] Add subcell limiting support for P4estMesh (#1954) * First part of P4estMesh support * Adapt `get_boundary_outer_state` to allow supersoniv cylinder eilixir * Add test * Rename routine * Rename other routine * Adapt elixir and complete testing * Adapt tests * Adapt cfl in elixir * Fix test * Add sedov elixir to test periodic boundaries * fmt * Activate test for coverage run * Remove extra lines of code * Implement suggestions * Adapt parameters to reduce bounds checking errors * Implement first suggestions * Remove `mesh` from `get_boundary_outer_state` * Use `foreach_enumerate` to remove allocations * Move `get_boundary_outer_state` to elixir --------- Co-authored-by: Michael Schlottke-Lakemper --- ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 102 +++++++ ...ir_euler_supersonic_cylinder_sc_subcell.jl | 162 +++++++++++ src/callbacks_stage/subcell_bounds_check.jl | 42 +++ .../subcell_bounds_check_2d.jl | 47 +--- .../subcell_limiter_idp_correction_2d.jl | 3 +- src/solvers/dgsem_p4est/dg.jl | 2 + .../dgsem_p4est/subcell_limiters_2d.jl | 257 ++++++++++++++++++ .../dg_2d_subcell_limiters.jl | 2 +- .../dgsem_structured/subcell_limiters_2d.jl | 62 +++-- .../dgsem_tree/dg_2d_subcell_limiters.jl | 37 ++- src/solvers/dgsem_tree/subcell_limiters.jl | 29 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 14 +- test/test_p4est_2d.jl | 52 ++++ 13 files changed, 717 insertions(+), 94 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl create mode 100644 examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl create mode 100644 src/solvers/dgsem_p4est/subcell_limiters_2d.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl new file mode 100644 index 00000000000..6a22aa262ef --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -0,0 +1,102 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +The Sedov blast wave setup based on Flash +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 +""" +function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # 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) + + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 + r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) + p0_outer = 1.0e-5 # = true Sedov setup + + # Calculate primitive variables + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_sedov_blast_wave + +# Get the DG approximation space +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, + min)], + max_iterations_newton = 40, # Default parameters are not sufficient to fulfill bounds properly. + newton_tolerances = (1.0e-14, 1.0e-15)) + +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) + +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg = polydeg, initial_refinement_level = 2, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 300 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 300, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) + +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_supersonic_cylinder_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl new file mode 100644 index 00000000000..78df12d60dd --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl @@ -0,0 +1,162 @@ +# Channel flow around a cylinder at Mach 3 +# +# Boundary conditions are supersonic Mach 3 inflow at the left portion of the domain +# and supersonic outflow at the right portion of the domain. The top and bottom of the +# channel as well as the cylinder are treated as Euler slip wall boundaries. +# This flow results in strong shock reflections / interactions as well as Kelvin-Helmholtz +# instabilities at later times as two Mach stems form above and below the cylinder. +# +# For complete details on the problem setup see Section 5.7 of the paper: +# - Jean-Luc Guermond, Murtazo Nazarov, Bojan Popov, and Ignacio Tomas (2018) +# Second-Order Invariant Domain Preserving Approximation of the Euler Equations using Convex Limiting. +# [DOI: 10.1137/17M1149961](https://doi.org/10.1137/17M1149961) +# +# Keywords: supersonic flow, shock capturing, AMR, unstructured curved mesh, positivity preservation, compressible Euler, 2D + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +@inline function initial_condition_mach3_flow(x, t, equations::CompressibleEulerEquations2D) + # set the freestream flow parameters + rho_freestream = 1.4 + v1 = 3.0 + v2 = 0.0 + p_freestream = 1.0 + + prim = SVector(rho_freestream, v1, v2, p_freestream) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach3_flow + +# Supersonic inflow boundary condition. +# Calculate the boundary flux entirely from the external solution state, i.e., set +# external solution state values for everything entering the domain. +@inline function boundary_condition_supersonic_inflow(u_inner, + normal_direction::AbstractVector, + x, t, surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach3_flow(x, t, equations) + flux = Trixi.flux(u_boundary, normal_direction, equations) + + return flux +end + +# For subcell limiting, the calculation of local bounds for non-periodic domains requires the +# boundary outer state. Those functions return the boundary value for a specific boundary condition +# at time `t`, for the node with spatial indices `indices` and the given `normal_direction`. +# only for P4estMesh{2} +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_supersonic_inflow), + normal_direction::AbstractVector, + equations, dg, cache, + indices...) + x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) + + return initial_condition_mach3_flow(x, t, equations) +end + +# Supersonic outflow boundary condition. +# Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow +# except all the solution state values are set from the internal solution as everything leaves the domain +@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + flux = Trixi.flux(u_inner, normal_direction, equations) + + return flux +end + +# only for P4estMesh{2} +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_outflow), + normal_direction::AbstractVector, + equations, dg, cache, + indices...) + return u_inner +end + +# only for P4estMesh{2} +@inline function Trixi.get_boundary_outer_state(u_inner, t, + boundary_condition::typeof(boundary_condition_slip_wall), + normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D, + dg, cache, indices...) + factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3]) + u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction + + return SVector(u_inner[1], + u_inner[2] - 2 * u_normal[1], + u_inner[3] - 2 * u_normal[2], + u_inner[4]) +end + +boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall, + :Circle => boundary_condition_slip_wall, + :Top => boundary_condition_slip_wall, + :Right => boundary_condition_outflow, + :Left => boundary_condition_supersonic_inflow) + +volume_flux = flux_ranocha_turbo +surface_flux = flux_lax_friedrichs +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + local_twosided_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure], + local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, + min)], + max_iterations_newton = 50) # Default value of 10 iterations is too low to fulfill bounds. + +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp", + joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp")) + +mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 0) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +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.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) + +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/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 6268fde6dd7..15231b25922 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -171,5 +171,47 @@ end return nothing end +@inline function save_bounds_check_errors(output_directory, time, iter, equations, + limiter::SubcellLimiterIDP) + (; local_twosided, positivity, local_onesided) = limiter + (; idp_bounds_delta_local) = limiter.cache + + # Print to output file + open(joinpath(output_directory, "deviations.txt"), "a") do f + print(f, iter, ", ", time) + if local_twosided + for v in limiter.local_twosided_variables_cons + v_string = string(v) + print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")], + ", ", idp_bounds_delta_local[Symbol(v_string, "_max")]) + end + end + if local_onesided + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear + key = Symbol(string(variable), "_", string(min_or_max)) + print(f, ", ", idp_bounds_delta_local[key]) + end + end + if positivity + for v in limiter.positivity_variables_cons + if v in limiter.local_twosided_variables_cons + continue + end + print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")]) + end + for variable in limiter.positivity_variables_nonlinear + print(f, ", ", idp_bounds_delta_local[Symbol(string(variable), "_min")]) + end + end + println(f) + end + # Reset local maximum deviations + for (key, _) in idp_bounds_delta_local + idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key])) + end + + return nothing +end + include("subcell_bounds_check_2d.jl") end # @muladd diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 9c3664ab0c3..8ac4c74b518 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -5,8 +5,8 @@ @muladd begin #! format: noindent -@inline function check_bounds(u, equations, solver, cache, - limiter::SubcellLimiterIDP) +@inline function check_bounds(u, equations::AbstractEquations{2}, # only works for 2D + solver, cache, limiter::SubcellLimiterIDP) (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache @@ -102,47 +102,4 @@ return nothing end - -@inline function save_bounds_check_errors(output_directory, time, iter, equations, - limiter::SubcellLimiterIDP) - (; local_twosided, positivity, local_onesided) = limiter - (; idp_bounds_delta_local) = limiter.cache - - # Print to output file - open(joinpath(output_directory, "deviations.txt"), "a") do f - print(f, iter, ", ", time) - if local_twosided - for v in limiter.local_twosided_variables_cons - v_string = string(v) - print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")], - ", ", idp_bounds_delta_local[Symbol(v_string, "_max")]) - end - end - if local_onesided - for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear - key = Symbol(string(variable), "_", string(min_or_max)) - print(f, ", ", idp_bounds_delta_local[key]) - end - end - if positivity - for v in limiter.positivity_variables_cons - if v in limiter.local_twosided_variables_cons - continue - end - print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")]) - end - for variable in limiter.positivity_variables_nonlinear - print(f, ", ", idp_bounds_delta_local[Symbol(string(variable), "_min")]) - end - end - println(f) - end - - # Reset local maximum deviations - for (key, _) in idp_bounds_delta_local - idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key])) - end - - return nothing -end end # @muladd diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index 0eb048ca5b8..e3b5616d5e5 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -6,7 +6,8 @@ #! format: noindent function perform_idp_correction!(u, dt, - mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + P4estMesh{2}}, equations, dg, cache) @unpack inverse_weights = dg.basis @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 10cc075089c..8197ad4a2d0 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -53,4 +53,6 @@ include("dg_2d_parabolic.jl") include("dg_3d.jl") include("dg_3d_parabolic.jl") include("dg_parallel.jl") + +include("subcell_limiters_2d.jl") end # @muladd diff --git a/src/solvers/dgsem_p4est/subcell_limiters_2d.jl b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl new file mode 100644 index 00000000000..95267f387f9 --- /dev/null +++ b/src/solvers/dgsem_p4est/subcell_limiters_2d.jl @@ -0,0 +1,257 @@ +# 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 +#! format: noindent + +function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, + mesh::P4estMesh{2}) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + + (; neighbor_ids, node_indices) = cache.interfaces + index_range = eachnode(dg) + + # Calc bounds at interfaces and periodic boundaries + for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + + # Create the local i,j indexing + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + i_secondary = i_secondary_start + j_secondary = j_secondary_start + + for node in eachnode(dg) + var_primary = u[variable, i_primary, j_primary, primary_element] + var_secondary = u[variable, i_secondary, j_secondary, secondary_element] + + var_min[i_primary, j_primary, primary_element] = min(var_min[i_primary, + j_primary, + primary_element], + var_secondary) + var_max[i_primary, j_primary, primary_element] = max(var_max[i_primary, + j_primary, + primary_element], + var_secondary) + + var_min[i_secondary, j_secondary, secondary_element] = min(var_min[i_secondary, + j_secondary, + secondary_element], + var_primary) + var_max[i_secondary, j_secondary, secondary_element] = max(var_max[i_secondary, + j_secondary, + secondary_element], + var_primary) + + # Increment primary element indices + i_primary += i_primary_step + j_primary += j_primary_step + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + # Calc bounds at physical boundaries + calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t, + boundary_conditions, + mesh, equations, dg, cache) + + return nothing +end + +@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t, + boundary_conditions::BoundaryConditionPeriodic, + mesh::P4estMesh{2}, + equations, dg, cache) + return nothing +end + +@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t, + boundary_conditions, + mesh::P4estMesh{2}, + equations, dg, cache) + (; boundary_condition_types, boundary_indices) = boundary_conditions + (; contravariant_vectors) = cache.elements + + (; boundaries) = cache + index_range = eachnode(dg) + + foreach_enumerate(boundary_condition_types) do (i, boundary_condition) + for boundary in boundary_indices[i] + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], + index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], + index_range) + + i_node = i_node_start + j_node = j_node_start + for i in eachnode(dg) + normal_direction = get_normal_direction(direction, + contravariant_vectors, + i_node, j_node, element) + + u_inner = get_node_vars(u, equations, dg, i_node, j_node, element) + + u_outer = get_boundary_outer_state(u_inner, t, boundary_condition, + normal_direction, + equations, dg, cache, + i_node, j_node, element) + var_outer = u_outer[variable] + + var_min[i_node, j_node, element] = min(var_min[i_node, j_node, element], + var_outer) + var_max[i_node, j_node, element] = max(var_max[i_node, j_node, element], + var_outer) + + i_node += i_node_step + j_node += j_node_step + end + end + end + + return nothing +end + +function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, semi, + mesh::P4estMesh{2}) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi + + (; neighbor_ids, node_indices) = cache.interfaces + index_range = eachnode(dg) + index_end = last(index_range) + + # Calc bounds at interfaces and periodic boundaries + for interface in eachinterface(dg, cache) + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + + # Create the local i,j indexing + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + i_secondary = i_secondary_start + j_secondary = j_secondary_start + + for node in eachnode(dg) + var_primary = variable(get_node_vars(u, equations, dg, i_primary, j_primary, + primary_element), equations) + var_secondary = variable(get_node_vars(u, equations, dg, i_secondary, + j_secondary, secondary_element), + equations) + + var_minmax[i_primary, j_primary, primary_element] = minmax(var_minmax[i_primary, + j_primary, + primary_element], + var_secondary) + var_minmax[i_secondary, j_secondary, secondary_element] = minmax(var_minmax[i_secondary, + j_secondary, + secondary_element], + var_primary) + + # Increment primary element indices + i_primary += i_primary_step + j_primary += j_primary_step + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + # Calc bounds at physical boundaries + calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t, + boundary_conditions, + mesh, equations, dg, cache) + + return nothing +end + +@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t, + boundary_conditions::BoundaryConditionPeriodic, + mesh::P4estMesh{2}, + equations, dg, cache) + return nothing +end + +@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t, + boundary_conditions, + mesh::P4estMesh{2}, + equations, dg, cache) + (; boundary_condition_types, boundary_indices) = boundary_conditions + (; contravariant_vectors) = cache.elements + + (; boundaries) = cache + index_range = eachnode(dg) + + foreach_enumerate(boundary_condition_types) do (i, boundary_condition) + for boundary in boundary_indices[i] + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], + index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], + index_range) + + i_node = i_node_start + j_node = j_node_start + for node in eachnode(dg) + normal_direction = get_normal_direction(direction, + contravariant_vectors, + i_node, j_node, element) + + u_inner = get_node_vars(u, equations, dg, i_node, j_node, element) + + u_outer = get_boundary_outer_state(u_inner, t, boundary_condition, + normal_direction, + equations, dg, cache, + i_node, j_node, element) + var_outer = variable(u_outer, equations) + + var_minmax[i_node, j_node, element] = minmax(var_minmax[i_node, j_node, + element], + var_outer) + + i_node += i_node_step + j_node += j_node_step + end + end + end + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index e10a5295459..b8ea14bf762 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -10,7 +10,7 @@ # # See also `flux_differencing_kernel!`. @inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, - mesh::StructuredMesh{2}, + mesh::Union{StructuredMesh{2}, P4estMesh{2}}, nonconservative_terms::False, equations, volume_flux, dg::DGSEM, element, cache) (; contravariant_vectors) = cache.elements diff --git a/src/solvers/dgsem_structured/subcell_limiters_2d.jl b/src/solvers/dgsem_structured/subcell_limiters_2d.jl index 5b65a475062..823c4817229 100644 --- a/src/solvers/dgsem_structured/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_structured/subcell_limiters_2d.jl @@ -9,6 +9,7 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, mesh::StructuredMesh{2}) _, equations, dg, cache = mesh_equations_solver_cache(semi) (; boundary_conditions) = semi + (; contravariant_vectors) = cache.elements # Calc bounds at interfaces and periodic boundaries for element in eachelement(dg, cache) @@ -56,8 +57,11 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, for cell_y in axes(mesh, 2) element = linear_indices[begin, cell_y] for j in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[1], - cache, t, equations, dg, + Ja1 = get_contravariant_vector(1, contravariant_vectors, 1, j, element) + u_inner = get_node_vars(u, equations, dg, 1, j, element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[1], Ja1, 1, + equations, dg, cache, 1, j, element) var_outer = u_outer[variable] @@ -69,8 +73,12 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, for cell_y in axes(mesh, 2) element = linear_indices[end, cell_y] for j in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[2], - cache, t, equations, dg, + Ja1 = get_contravariant_vector(1, contravariant_vectors, nnodes(dg), j, + element) + u_inner = get_node_vars(u, equations, dg, nnodes(dg), j, element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[2], Ja1, 2, + equations, dg, cache, nnodes(dg), j, element) var_outer = u_outer[variable] @@ -86,8 +94,11 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, for cell_x in axes(mesh, 1) element = linear_indices[cell_x, begin] for i in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[3], - cache, t, equations, dg, + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, 1, element) + u_inner = get_node_vars(u, equations, dg, i, 1, element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[3], Ja2, 3, + equations, dg, cache, i, 1, element) var_outer = u_outer[variable] @@ -99,8 +110,12 @@ function calc_bounds_twosided_interface!(var_min, var_max, variable, u, t, semi, for cell_x in axes(mesh, 1) element = linear_indices[cell_x, end] for i in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[4], - cache, t, equations, dg, + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, nnodes(dg), + element) + u_inner = get_node_vars(u, equations, dg, i, nnodes(dg), element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[4], Ja2, 4, + equations, dg, cache, i, nnodes(dg), element) var_outer = u_outer[variable] @@ -119,6 +134,7 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem mesh::StructuredMesh{2}) _, equations, dg, cache = mesh_equations_solver_cache(semi) (; boundary_conditions) = semi + (; contravariant_vectors) = cache.elements # Calc bounds at interfaces and periodic boundaries for element in eachelement(dg, cache) @@ -163,8 +179,11 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem for cell_y in axes(mesh, 2) element = linear_indices[begin, cell_y] for j in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[1], - cache, t, equations, dg, + Ja1 = get_contravariant_vector(1, contravariant_vectors, 1, j, element) + u_inner = get_node_vars(u, equations, dg, 1, j, element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[1], Ja1, 1, + equations, dg, cache, 1, j, element) var_outer = variable(u_outer, equations) @@ -175,8 +194,12 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem for cell_y in axes(mesh, 2) element = linear_indices[end, cell_y] for j in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[2], - cache, t, equations, dg, + Ja1 = get_contravariant_vector(1, contravariant_vectors, nnodes(dg), j, + element) + u_inner = get_node_vars(u, equations, dg, nnodes(dg), j, element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[2], Ja1, 2, + equations, dg, cache, nnodes(dg), j, element) var_outer = variable(u_outer, equations) @@ -191,8 +214,11 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem for cell_x in axes(mesh, 1) element = linear_indices[cell_x, begin] for i in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[3], - cache, t, equations, dg, + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, 1, element) + u_inner = get_node_vars(u, equations, dg, i, 1, element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[3], Ja2, 3, + equations, dg, cache, i, 1, element) var_outer = variable(u_outer, equations) @@ -203,8 +229,12 @@ function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u, t, sem for cell_x in axes(mesh, 1) element = linear_indices[cell_x, end] for i in eachnode(dg) - u_outer = get_boundary_outer_state(boundary_conditions[4], - cache, t, equations, dg, + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, nnodes(dg), + element) + u_inner = get_node_vars(u, equations, dg, i, nnodes(dg), element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[4], Ja2, 4, + equations, dg, cache, i, nnodes(dg), element) var_outer = variable(u_outer, equations) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 998790ef8e1..c27327be8ad 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -5,7 +5,7 @@ @muladd begin #! format: noindent -function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}}, +function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}}, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DG, uEltype) cache = create_cache(mesh, equations, @@ -57,7 +57,8 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}}, end function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + P4estMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache) @@ -72,7 +73,8 @@ function calc_volume_integral!(du, u, end @inline function subcell_limiting_kernel!(du, u, element, - mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + P4estMesh{2}}, nonconservative_terms, equations, volume_integral, limiter::SubcellLimiterIDP, dg::DGSEM, cache) @@ -392,7 +394,9 @@ end # Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems. @inline function calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fstar1_L, fstar1_R, fstar2_L, fstar2_R, - u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + P4estMesh{2}}, nonconservative_terms::False, equations, limiter::SubcellLimiterIDP, dg, element, cache) @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes @@ -430,7 +434,9 @@ end # Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems. @inline function calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fstar1_L, fstar1_R, fstar2_L, fstar2_R, - u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + P4estMesh{2}}, nonconservative_terms::True, equations, limiter::SubcellLimiterIDP, dg, element, cache) @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes @@ -466,17 +472,24 @@ end end """ - get_boundary_outer_state(boundary_condition::BoundaryConditionDirichlet, - cache, t, equations, dg, indices...) -For subcell limiting, the calculation of local bounds for non-periodic domains require the boundary -outer state. This function returns the boundary value at time `t` and for node with spatial -indices `indices`. + get_boundary_outer_state(u_inner, t, + boundary_condition::BoundaryConditionDirichlet, + orientation_or_normal, direction, + equations, dg, cache, indices...) +For subcell limiting, the calculation of local bounds for non-periodic domains requires the boundary +outer state. This function returns the boundary value for [`BoundaryConditionDirichlet`](@ref) at +time `t` and for node with spatial indices `indices` at the boundary with `orientation_or_normal` +and `direction`. + +Should be used together with [`TreeMesh`](@ref) or [`StructuredMesh`](@ref). !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. """ -@inline function get_boundary_outer_state(boundary_condition::BoundaryConditionDirichlet, - cache, t, equations, dg, indices...) +@inline function get_boundary_outer_state(u_inner, t, + boundary_condition::BoundaryConditionDirichlet, + orientation_or_normal, direction, + equations, dg, cache, indices...) (; node_coordinates) = cache.elements x = get_node_coords(node_coordinates, equations, dg, indices...) diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl index 17c9488316d..f560018f461 100644 --- a/src/solvers/dgsem_tree/subcell_limiters.jl +++ b/src/solvers/dgsem_tree/subcell_limiters.jl @@ -193,28 +193,27 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP) else setup = ["Limiter" => ""] if local_twosided - setup = [ - setup..., - "" => "Local two-sided limiting for conservative variables $(limiter.local_twosided_variables_cons)", - ] + push!(setup, + "" => "Local two-sided limiting for conservative variables $(limiter.local_twosided_variables_cons)") end if positivity - string = "Positivity limiting for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)" - setup = [setup..., "" => string] - setup = [ - setup..., - "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)", - ] + if !isempty(limiter.positivity_variables_cons) + string = "conservative variables $(limiter.positivity_variables_cons)" + push!(setup, "" => "Positivity limiting for " * string) + end + if !isempty(limiter.positivity_variables_nonlinear) + string = "$(limiter.positivity_variables_nonlinear)" + push!(setup, "" => "Positivity limiting for " * string) + end + push!(setup, + "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)") end if local_onesided for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear - setup = [setup..., "" => "Local $min_or_max limiting for $variable"] + push!(setup, "" => "Local $min_or_max limiting for $variable") end end - setup = [ - setup..., - "Local bounds" => "FV solution", - ] + push!(setup, "Local bounds" => "FV solution") end summary_box(io, "SubcellLimiterIDP", setup) end diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index b185cd4c540..5dae5798a83 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -155,8 +155,11 @@ end index = reverse(index) boundary_index += 2 end - u_outer = get_boundary_outer_state(boundary_conditions[boundary_index], - cache, t, equations, dg, + u_inner = get_node_vars(u, equations, dg, index..., element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[boundary_index], + orientation, boundary_index, + equations, dg, cache, index..., element) var_outer = u_outer[variable] @@ -258,8 +261,11 @@ end index = reverse(index) boundary_index += 2 end - u_outer = get_boundary_outer_state(boundary_conditions[boundary_index], - cache, t, equations, dg, + u_inner = get_node_vars(u, equations, dg, index..., element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[boundary_index], + orientation, boundary_index, + equations, dg, cache, index..., element) var_outer = variable(u_outer, equations) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 95c69b86954..aad6337ffcd 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -274,6 +274,32 @@ end end end +@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_sedov_blast_wave_sc_subcell.jl"), + l2=[ + 0.4573787784168518, + 0.28520972760728397, + 0.28527281808006966, + 1.2881460122982442, + ], + linf=[ + 1.644411040701827, + 1.6743368119653912, + 1.6760847977977988, + 6.268843623142863, + ], + tspan=(0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_euler_sedov.jl with HLLC Flux" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[ @@ -462,6 +488,32 @@ end end end +@trixi_testset "elixir_euler_supersonic_cylinder_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_supersonic_cylinder_sc_subcell.jl"), + l2=[ + 0.11085870166618325, + 0.23309905989870722, + 0.13505351590735631, + 0.7932047512585592, + ], + linf=[ + 2.9808773737943564, + 4.209364526217892, + 6.265341002817672, + 24.077904874883338, + ], + tspan=(0.0, 0.02)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end +end + @trixi_testset "elixir_euler_NACA6412airfoil_mach2.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA6412airfoil_mach2.jl"), l2=[