From b1cdce79086e2bbecfcb88022dde2f6e4921db74 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 14 May 2024 13:35:04 +0200 Subject: [PATCH 01/19] Add structured mesh support --- .../elixir_euler_free_stream_sc_subcell.jl | 84 +++++ .../elixir_euler_source_terms_sc_subcell.jl | 73 ++++ .../subcell_limiter_idp_correction_2d.jl | 71 ++-- src/solvers/dgsem_structured/dg.jl | 3 + .../dg_2d_subcell_limiters.jl | 111 ++++++ .../dgsem_structured/subcell_limiters_2d.jl | 249 ++++++++++++ .../dgsem_tree/dg_2d_subcell_limiters.jl | 11 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 356 ++++++++++++------ test/test_structured_2d.jl | 53 +++ 9 files changed, 872 insertions(+), 139 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl create mode 100644 examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl create mode 100644 src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl create mode 100644 src/solvers/dgsem_structured/subcell_limiters_2d.jl diff --git a/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl new file mode 100644 index 00000000000..a0ddcb64eb8 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl @@ -0,0 +1,84 @@ + +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]) + +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 + +cells_per_dimension = (16, 16) +mesh = StructuredMesh(cells_per_dimension, mapping, 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()) + +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 diff --git a/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl new file mode 100644 index 00000000000..97c541266fa --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl @@ -0,0 +1,73 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_convergence_test +source_terms = source_terms_convergence_test + +# 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)]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +# Waving flag +f1(s) = SVector(-1.0, s - 1.0) +f2(s) = SVector(1.0, s + 1.0) +f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) +f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) +mapping = Trixi.transfinite_mapping((f1, f2, f3, f4)) + +cells_per_dimension = (16, 16) +mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms) + +############################################################################### +# 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 = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index 6f1723e2a98..9c4a31a515c 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -6,39 +6,60 @@ #! format: noindent function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache) - @unpack inverse_weights = dg.basis - @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes - @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients - @threaded for element in eachelement(dg, cache) # Sign switch as in apply_jacobian! inverse_jacobian = -cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) - # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1} - alpha_flux1 = (1 - alpha1[i, j, element]) * - get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, - element) - alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) * - get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, - j, element) - alpha_flux2 = (1 - alpha2[i, j, element]) * - get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, - element) - alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) * - get_node_vars(antidiffusive_flux2_L, equations, dg, i, - j + 1, element) - - for v in eachvariable(equations) - u[v, i, j, element] += dt * inverse_jacobian * - (inverse_weights[i] * - (alpha_flux1_ip1[v] - alpha_flux1[v]) + - inverse_weights[j] * - (alpha_flux2_jp1[v] - alpha_flux2[v])) - end + perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, cache, + i, j, element) end end return nothing end + +function perform_idp_correction!(u, dt, mesh::StructuredMesh{2}, equations, dg, cache) + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + # Sign switch as in apply_jacobian! + inverse_jacobian = -cache.elements.inverse_jacobian[i, j, element] + + perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, cache, + i, j, element) + end + end + + return nothing +end + +# Function barrier to dispatch outer function by mesh type +@inline function perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, + cache, i, j, element) + (; inverse_weights) = dg.basis + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + (; alpha1, alpha2) = dg.volume_integral.limiter.cache.subcell_limiter_coefficients + + # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1} + alpha_flux1 = (1 - alpha1[i, j, element]) * + get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, element) + alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) * + get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j, + element) + alpha_flux2 = (1 - alpha2[i, j, element]) * + get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, element) + alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) * + get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1, + element) + + for v in eachvariable(equations) + u[v, i, j, element] += dt * inverse_jacobian * + (inverse_weights[i] * + (alpha_flux1_ip1[v] - alpha_flux1[v]) + + inverse_weights[j] * + (alpha_flux2_jp1[v] - alpha_flux2[v])) + end + + return nothing +end end # @muladd diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 00e321fba65..ecb2485c68b 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -85,6 +85,9 @@ include("indicators_1d.jl") include("indicators_2d.jl") include("indicators_3d.jl") +include("subcell_limiters_2d.jl") +include("dg_2d_subcell_limiters.jl") + # Specialized implementations used to improve performance include("dg_2d_compressible_euler.jl") include("dg_3d_compressible_euler.jl") diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl new file mode 100644 index 00000000000..4da402425ea --- /dev/null +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -0,0 +1,111 @@ +# 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 + +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**without non-conservative terms**). +# +# See also `flux_differencing_kernel!`. +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, + mesh::StructuredMesh{2}, + nonconservative_terms::False, equations, + volume_flux, dg::DGSEM, element, cache) + (; contravariant_vectors) = cache.elements + (; weights, derivative_split) = dg.basis + (; flux_temp_threaded) = cache + + flux_temp = flux_temp_threaded[Threads.threadid()] + + # The FV-form fluxes are calculated in a recursive manner, i.e.: + # fhat_(0,1) = w_0 * FVol_0, + # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j, for j=1,...,N-1, + # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i). + + # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated + # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing` + # and saved in in `flux_temp`. + + # Split form volume flux in orientation 1: x direction + flux_temp .= zero(eltype(flux_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + # pull the contravariant vectors in each coordinate direction + Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element) # x direction + + # All diagonal entries of `derivative_split` are zero. Thus, we can skip + # the computation of the diagonal terms. In addition, we use the symmetry + # of the `volume_flux` to save half of the possible two-point flux + # computations. + + # x direction + for ii in (i + 1):nnodes(dg) + u_node_ii = get_node_vars(u, equations, dg, ii, j, element) + # pull the contravariant vectors and compute the average + Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, + element) + Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii) + + # compute the contravariant sharp flux in the direction of the averaged contravariant vector + fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1, + equations, dg, ii, j) + end + end + + # FV-form flux `fhat` in x direction + fhat1_L[:, 1, :] .= zero(eltype(fhat1_L)) + fhat1_L[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_L)) + fhat1_R[:, 1, :] .= zero(eltype(fhat1_R)) + fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) + + for j in eachnode(dg), i in 1:(nnodes(dg) - 1), v in eachvariable(equations) + fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j] + fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j] + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + # pull the contravariant vectors in each coordinate direction + Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element) + + # y direction + for jj in (j + 1):nnodes(dg) + u_node_jj = get_node_vars(u, equations, dg, i, jj, element) + # pull the contravariant vectors and compute the average + Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, + element) + Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj) + # compute the contravariant sharp flux in the direction of the averaged contravariant vector + fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2, + equations, dg, i, jj) + end + end + + # FV-form flux `fhat` in y direction + fhat2_L[:, :, 1] .= zero(eltype(fhat2_L)) + fhat2_L[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_L)) + fhat2_R[:, :, 1] .= zero(eltype(fhat2_R)) + fhat2_R[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_R)) + + for j in 1:(nnodes(dg) - 1), i in eachnode(dg), v in eachvariable(equations) + fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j] + fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1] + end + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem_structured/subcell_limiters_2d.jl b/src/solvers/dgsem_structured/subcell_limiters_2d.jl new file mode 100644 index 00000000000..83a979a16a8 --- /dev/null +++ b/src/solvers/dgsem_structured/subcell_limiters_2d.jl @@ -0,0 +1,249 @@ +# 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::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) + # Get neighboring element ids + left = cache.elements.left_neighbors[1, element] + lower = cache.elements.left_neighbors[2, element] + + if left != 0 + for j in eachnode(dg) + var_left = u[variable, nnodes(dg), j, left] + var_element = u[variable, 1, j, element] + + var_min[1, j, element] = min(var_min[1, j, element], var_left) + var_max[1, j, element] = max(var_max[1, j, element], var_left) + + var_min[nnodes(dg), j, left] = min(var_min[nnodes(dg), j, left], + var_element) + var_max[nnodes(dg), j, left] = max(var_max[nnodes(dg), j, left], + var_element) + end + end + if lower != 0 + for i in eachnode(dg) + var_lower = u[variable, i, nnodes(dg), lower] + var_element = u[variable, i, 1, element] + + var_min[i, 1, element] = min(var_min[i, 1, element], var_lower) + var_max[i, 1, element] = max(var_max[i, 1, element], var_lower) + + var_min[i, nnodes(dg), lower] = min(var_min[i, nnodes(dg), lower], + var_element) + var_max[i, nnodes(dg), lower] = max(var_max[i, nnodes(dg), lower], + var_element) + end + end + end + + # Calc bounds at physical boundaries + if isperiodic(mesh) + return nothing + end + linear_indices = LinearIndices(size(mesh)) + if !isperiodic(mesh, 1) + # - xi direction + for cell_y in axes(mesh, 2) + element = linear_indices[begin, cell_y] + for j in eachnode(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, cache, t, + boundary_conditions[1], Ja1, 1, + mesh, equations, dg, + 1, j, element) + var_outer = u_outer[variable] + + var_min[1, j, element] = min(var_min[1, j, element], var_outer) + var_max[1, j, element] = max(var_max[1, j, element], var_outer) + end + end + # + xi direction + for cell_y in axes(mesh, 2) + element = linear_indices[end, cell_y] + for j in eachnode(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, cache, t, + boundary_conditions[2], Ja1, 2, + mesh, equations, dg, + nnodes(dg), j, element) + var_outer = u_outer[variable] + + var_min[nnodes(dg), j, element] = min(var_min[nnodes(dg), j, element], + var_outer) + var_max[nnodes(dg), j, element] = max(var_max[nnodes(dg), j, element], + var_outer) + end + end + end + if !isperiodic(mesh, 2) + # - eta direction + for cell_x in axes(mesh, 1) + element = linear_indices[cell_x, begin] + for i in eachnode(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, cache, t, + boundary_conditions[3], Ja2, 3, + mesh, equations, dg, + i, 1, element) + var_outer = u_outer[variable] + + var_min[i, 1, element] = min(var_min[i, 1, element], var_outer) + var_max[i, 1, element] = max(var_max[i, 1, element], var_outer) + end + end + # - eta direction + for cell_x in axes(mesh, 1) + element = linear_indices[cell_x, end] + for i in eachnode(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, cache, t, + boundary_conditions[4], Ja2, 4, + mesh, equations, dg, + i, nnodes(dg), element) + var_outer = u_outer[variable] + + var_min[i, nnodes(dg), element] = min(var_min[i, nnodes(dg), element], + var_outer) + var_max[i, nnodes(dg), element] = max(var_max[i, nnodes(dg), element], + var_outer) + end + end + end + + return nothing +end + +function calc_bounds_onesided_interface!(var_minmax, minmax, 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) + # Get neighboring element ids + left = cache.elements.left_neighbors[1, element] + lower = cache.elements.left_neighbors[2, element] + + if left != 0 + for j in eachnode(dg) + var_left = variable(get_node_vars(u, equations, dg, nnodes(dg), j, + left), equations) + var_element = variable(get_node_vars(u, equations, dg, 1, j, element), + equations) + + var_minmax[1, j, element] = minmax(var_minmax[1, j, element], var_left) + var_minmax[nnodes(dg), j, left] = minmax(var_minmax[nnodes(dg), j, + left], var_element) + end + end + if lower != 0 + for i in eachnode(dg) + var_lower = variable(get_node_vars(u, equations, dg, i, nnodes(dg), + lower), equations) + var_element = variable(get_node_vars(u, equations, dg, i, 1, element), + equations) + + var_minmax[i, 1, element] = minmax(var_minmax[i, 1, element], var_lower) + var_minmax[i, nnodes(dg), lower] = minmax(var_minmax[i, nnodes(dg), + lower], + var_element) + end + end + end + + # Calc bounds at physical boundaries + if isperiodic(mesh) + return nothing + end + linear_indices = LinearIndices(size(mesh)) + if !isperiodic(mesh, 1) + # - xi direction + for cell_y in axes(mesh, 2) + element = linear_indices[begin, cell_y] + for j in eachnode(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, cache, t, + boundary_conditions[1], Ja1, 1, + mesh, equations, dg, + 1, j, element) + var_outer = variable(u_outer, equations) + + var_minmax[1, j, element] = minmax(var_minmax[1, j, element], var_outer) + end + end + # + xi direction + for cell_y in axes(mesh, 2) + element = linear_indices[end, cell_y] + for j in eachnode(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, cache, t, + boundary_conditions[2], Ja1, 2, + mesh, equations, dg, + nnodes(dg), j, element) + var_outer = variable(u_outer, equations) + + var_minmax[nnodes(dg), j, element] = minmax(var_minmax[nnodes(dg), j, + element], + var_outer) + end + end + end + if !isperiodic(mesh, 2) + # - eta direction + for cell_x in axes(mesh, 1) + element = linear_indices[cell_x, begin] + for i in eachnode(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, cache, t, + boundary_conditions[3], Ja2, 3, + mesh, equations, dg, + i, 1, element) + var_outer = variable(u_outer, equations) + + var_minmax[i, 1, element] = minmax(var_minmax[i, 1, element], var_outer) + end + end + # + eta direction + for cell_x in axes(mesh, 1) + element = linear_indices[cell_x, end] + for i in eachnode(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, cache, t, + boundary_conditions[4], Ja2, 4, + mesh, equations, dg, + i, nnodes(dg), element) + var_outer = variable(u_outer, equations) + + var_minmax[i, nnodes(dg), element] = minmax(var_minmax[i, nnodes(dg), + element], + var_outer) + end + end + end + + return nothing +end +end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 9af8b65b4cd..ba029667c17 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -5,8 +5,9 @@ @muladd begin #! format: noindent -function create_cache(mesh::TreeMesh{2}, equations, - volume_integral::VolumeIntegralSubcellLimiting, dg::DG, uEltype) +function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + equations, volume_integral::VolumeIntegralSubcellLimiting, + dg::DG, uEltype) cache = create_cache(mesh, equations, VolumeIntegralPureLGLFiniteVolume(volume_integral.volume_flux_fv), dg, uEltype) @@ -56,7 +57,7 @@ function create_cache(mesh::TreeMesh{2}, equations, end function calc_volume_integral!(du, u, - mesh::TreeMesh{2}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache) @@ -70,8 +71,8 @@ function calc_volume_integral!(du, u, end end -@inline function subcell_limiting_kernel!(du, u, - element, mesh::TreeMesh{2}, +@inline function subcell_limiting_kernel!(du, u, element, + mesh::Union{TreeMesh{2}, StructuredMesh{2}}, nonconservative_terms, equations, volume_integral, limiter::SubcellLimiterIDP, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 33ae0599748..7dd0abe3cad 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -275,17 +275,18 @@ end # Local two-sided limiting of conservative variables @inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi) + mesh, _, _, _ = mesh_equations_solver_cache(semi) + for variable in limiter.local_twosided_variables_cons - idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) + idp_local_twosided!(alpha, limiter, u, t, dt, semi, mesh, variable) end return nothing end -@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) +@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, mesh::TreeMesh{2}, + variable) _, _, dg, cache = mesh_equations_solver_cache(semi) - (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes - (; inverse_weights) = dg.basis (; variable_bounds) = limiter.cache.subcell_limiter_coefficients variable_string = string(variable) @@ -296,66 +297,103 @@ end @threaded for element in eachelement(dg, semi.cache) inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) - var = u[variable, i, j, element] - # Real Zalesak type limiter - # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" - # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" - # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is - # for each interface, not each node - - Qp = max(0, (var_max[i, j, element] - var) / dt) - Qm = min(0, (var_min[i, j, element] - var) / dt) - - # Calculate Pp and Pm - # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. - val_flux1_local = inverse_weights[i] * - antidiffusive_flux1_R[variable, i, j, element] - val_flux1_local_ip1 = -inverse_weights[i] * - antidiffusive_flux1_L[variable, i + 1, j, element] - val_flux2_local = inverse_weights[j] * - antidiffusive_flux2_R[variable, i, j, element] - val_flux2_local_jp1 = -inverse_weights[j] * - antidiffusive_flux2_L[variable, i, j + 1, element] - - Pp = max(0, val_flux1_local) + max(0, val_flux1_local_ip1) + - max(0, val_flux2_local) + max(0, val_flux2_local_jp1) - Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + - min(0, val_flux2_local) + min(0, val_flux2_local_jp1) - - Qp = max(0, (var_max[i, j, element] - var) / dt) - Qm = min(0, (var_min[i, j, element] - var) / dt) - - Pp = inverse_jacobian * Pp - Pm = inverse_jacobian * Pm - - # Compute blending coefficient avoiding division by zero - # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) - Qp = abs(Qp) / - (abs(Pp) + eps(typeof(Qp)) * 100 * abs(var_max[i, j, element])) - Qm = abs(Qm) / - (abs(Pm) + eps(typeof(Qm)) * 100 * abs(var_max[i, j, element])) - - # Calculate alpha at nodes - alpha[i, j, element] = max(alpha[i, j, element], 1 - min(1, Qp, Qm)) + idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache, + variable, var_min, var_max, i, j, element) + end + end + + return nothing +end + +@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, + mesh::StructuredMesh{2}, + variable) + _, _, dg, cache = mesh_equations_solver_cache(semi) + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + variable_string = string(variable) + var_min = variable_bounds[Symbol(variable_string, "_min")] + var_max = variable_bounds[Symbol(variable_string, "_max")] + calc_bounds_twosided!(var_min, var_max, variable, u, t, semi) + + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] + idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache, + variable, var_min, var_max, i, j, element) end end return nothing end +# Function barrier to dispatch outer function by mesh type +@inline function idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache, + variable, var_min, var_max, i, j, element) + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + (; inverse_weights) = dg.basis + + var = u[variable, i, j, element] + # Real Zalesak type limiter + # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" + # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" + # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is + # for each interface, not each node + + Qp = max(0, (var_max[i, j, element] - var) / dt) + Qm = min(0, (var_min[i, j, element] - var) / dt) + + # Calculate Pp and Pm + # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. + val_flux1_local = inverse_weights[i] * + antidiffusive_flux1_R[variable, i, j, element] + val_flux1_local_ip1 = -inverse_weights[i] * + antidiffusive_flux1_L[variable, i + 1, j, element] + val_flux2_local = inverse_weights[j] * + antidiffusive_flux2_R[variable, i, j, element] + val_flux2_local_jp1 = -inverse_weights[j] * + antidiffusive_flux2_L[variable, i, j + 1, element] + + Pp = max(0, val_flux1_local) + max(0, val_flux1_local_ip1) + + max(0, val_flux2_local) + max(0, val_flux2_local_jp1) + Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + + min(0, val_flux2_local) + min(0, val_flux2_local_jp1) + + Qp = max(0, (var_max[i, j, element] - var) / dt) + Qm = min(0, (var_min[i, j, element] - var) / dt) + + Pp = inverse_jacobian * Pp + Pm = inverse_jacobian * Pm + + # Compute blending coefficient avoiding division by zero + # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) + Qp = abs(Qp) / + (abs(Pp) + eps(typeof(Qp)) * 100 * abs(var_max[i, j, element])) + Qm = abs(Qm) / + (abs(Pm) + eps(typeof(Qm)) * 100 * abs(var_max[i, j, element])) + + # Calculate alpha at nodes + alpha[i, j, element] = max(alpha[i, j, element], 1 - min(1, Qp, Qm)) + + return nothing +end + ############################################################################## # Local one-sided limiting of nonlinear variables @inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi) + mesh, _, _, _ = mesh_equations_solver_cache(semi) + for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear - idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, min_or_max) + idp_local_onesided!(alpha, limiter, u, t, dt, semi, mesh, + variable, min_or_max) end return nothing end -@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, - min_or_max) +@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, mesh::TreeMesh{2}, + variable, min_or_max) _, equations, dg, cache = mesh_equations_solver_cache(semi) (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_minmax = variable_bounds[Symbol(string(variable), "_", string(min_or_max))] @@ -377,10 +415,36 @@ end return nothing end +@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, + mesh::StructuredMesh{2}, + variable, min_or_max) + _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_minmax = variable_bounds[Symbol(string(variable), "_", string(min_or_max))] + calc_bounds_onesided!(var_minmax, min_or_max, variable, u, t, semi) + + # Perform Newton's bisection method to find new alpha + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] + u_local = get_node_vars(u, equations, dg, i, j, element) + newton_loops_alpha!(alpha, var_minmax[i, j, element], u_local, + i, j, element, variable, min_or_max, + initial_check_local_onesided_newton_idp, + final_check_local_onesided_newton_idp, inverse_jacobian, + dt, equations, dg, cache, limiter) + end + end + + return nothing +end + ############################################################################### # Global positivity limiting @inline function idp_positivity!(alpha, limiter, u, dt, semi) + mesh, _, _, _ = mesh_equations_solver_cache(semi) + # Conservative variables for variable in limiter.positivity_variables_cons @trixi_timeit timer() "conservative variables" idp_positivity_conservative!(alpha, @@ -388,6 +452,7 @@ end u, dt, semi, + mesh, variable) end @@ -397,6 +462,7 @@ end limiter, u, dt, semi, + mesh, variable) end @@ -406,11 +472,9 @@ end ############################################################################### # Global positivity limiting of conservative variables -@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, variable) - mesh, equations, dg, cache = mesh_equations_solver_cache(semi) - (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes - (; inverse_weights) = dg.basis - (; positivity_correction_factor) = limiter +@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, + mesh::TreeMesh{2}, variable) + _, _, dg, cache = mesh_equations_solver_cache(semi) (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_min = variable_bounds[Symbol(string(variable), "_min")] @@ -418,84 +482,158 @@ end @threaded for element in eachelement(dg, semi.cache) inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) - var = u[variable, i, j, element] - if var < 0 - error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") - end + idp_positivity_conservative_inner!(alpha, inverse_jacobian, limiter, u, dt, + dg, cache, variable, var_min, + i, j, element) + end + end - # Compute bound - if limiter.local_twosided && - variable in limiter.local_twosided_variables_cons && - var_min[i, j, element] >= positivity_correction_factor * var - # Local limiting is more restrictive that positivity limiting - # => Skip positivity limiting for this node - continue - end - var_min[i, j, element] = positivity_correction_factor * var - - # Real one-sided Zalesak-type limiter - # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" - # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" - # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is - # for each interface, not each node - Qm = min(0, (var_min[i, j, element] - var) / dt) - - # Calculate Pm - # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. - val_flux1_local = inverse_weights[i] * - antidiffusive_flux1_R[variable, i, j, element] - val_flux1_local_ip1 = -inverse_weights[i] * - antidiffusive_flux1_L[variable, i + 1, j, element] - val_flux2_local = inverse_weights[j] * - antidiffusive_flux2_R[variable, i, j, element] - val_flux2_local_jp1 = -inverse_weights[j] * - antidiffusive_flux2_L[variable, i, j + 1, element] - - Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + - min(0, val_flux2_local) + min(0, val_flux2_local_jp1) - Pm = inverse_jacobian * Pm - - # Compute blending coefficient avoiding division by zero - # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) - Qm = abs(Qm) / (abs(Pm) + eps(typeof(Qm)) * 100) - - # Calculate alpha - alpha[i, j, element] = max(alpha[i, j, element], 1 - Qm) + return nothing +end + +@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, + mesh::Union{StructuredMesh{2}, + P4estMesh{2}}, + variable) + _, _, dg, cache = mesh_equations_solver_cache(semi) + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_min = variable_bounds[Symbol(string(variable), "_min")] + + @threaded for element in eachelement(dg, semi.cache) + for j in eachnode(dg), i in eachnode(dg) + inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] + idp_positivity_conservative_inner!(alpha, inverse_jacobian, limiter, u, dt, + dg, cache, variable, var_min, + i, j, element) end end return nothing end -@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) - _, equations, dg, cache = mesh_equations_solver_cache(semi) +# Function barrier to dispatch outer function by mesh type +@inline function idp_positivity_conservative_inner!(alpha, inverse_jacobian, limiter, u, + dt, dg, cache, variable, var_min, + i, j, element) + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + (; inverse_weights) = dg.basis (; positivity_correction_factor) = limiter + var = u[variable, i, j, element] + if var < 0 + error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") + end + + # Compute bound + if limiter.local_twosided && + variable in limiter.local_twosided_variables_cons && + var_min[i, j, element] >= positivity_correction_factor * var + # Local limiting is more restrictive that positivity limiting + # => Skip positivity limiting for this node + return nothing + end + var_min[i, j, element] = positivity_correction_factor * var + + # Real one-sided Zalesak-type limiter + # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" + # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" + # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is + # for each interface, not each node + Qm = min(0, (var_min[i, j, element] - var) / dt) + + # Calculate Pm + # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. + val_flux1_local = inverse_weights[i] * + antidiffusive_flux1_R[variable, i, j, element] + val_flux1_local_ip1 = -inverse_weights[i] * + antidiffusive_flux1_L[variable, i + 1, j, element] + val_flux2_local = inverse_weights[j] * + antidiffusive_flux2_R[variable, i, j, element] + val_flux2_local_jp1 = -inverse_weights[j] * + antidiffusive_flux2_L[variable, i, j + 1, element] + + Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + + min(0, val_flux2_local) + min(0, val_flux2_local_jp1) + Pm = inverse_jacobian * Pm + + # Compute blending coefficient avoiding division by zero + # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) + Qm = abs(Qm) / (abs(Pm) + eps(typeof(Qm)) * 100) + + # Calculate alpha + alpha[i, j, element] = max(alpha[i, j, element], 1 - Qm) + + return nothing +end + +############################################################################### +# Global positivity limiting of nonlinear variables + +@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, + mesh::TreeMesh{2}, variable) + _, _, dg, cache = mesh_equations_solver_cache(semi) + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_min = variable_bounds[Symbol(string(variable), "_min")] @threaded for element in eachelement(dg, semi.cache) inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) - # Compute bound - u_local = get_node_vars(u, equations, dg, i, j, element) - var = variable(u_local, equations) - if var < 0 - error("Safe low-order method produces negative value for variable $variable. Try a smaller time step.") - end - var_min[i, j, element] = positivity_correction_factor * var + idp_positivity_nonlinear_inner!(alpha, inverse_jacobian, limiter, u, dt, + semi, dg, cache, variable, var_min, + i, j, element) + end + end - # Perform Newton's bisection method to find new alpha - newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, - variable, min, initial_check_nonnegative_newton_idp, - final_check_nonnegative_newton_idp, inverse_jacobian, - dt, equations, dg, cache, limiter) + return nothing +end + +@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, + mesh::StructuredMesh{2}, + variable) + _, _, dg, cache = mesh_equations_solver_cache(semi) + + (; variable_bounds) = limiter.cache.subcell_limiter_coefficients + var_min = variable_bounds[Symbol(string(variable), "_min")] + + @threaded for element in eachelement(dg, semi.cache) + for j in eachnode(dg), i in eachnode(dg) + inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] + idp_positivity_nonlinear_inner!(alpha, inverse_jacobian, limiter, u, dt, + semi, dg, cache, variable, var_min, + i, j, element) end end return nothing end +# Function barrier to dispatch outer function by mesh type +@inline function idp_positivity_nonlinear_inner!(alpha, inverse_jacobian, limiter, u, + dt, semi, dg, cache, variable, var_min, + i, j, element) + _, equations, _, _ = mesh_equations_solver_cache(semi) + + u_local = get_node_vars(u, equations, dg, i, j, element) + var = variable(u_local, equations) + if var < 0 + error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") + end + var_min[i, j, element] = limiter.positivity_correction_factor * var + + # Perform Newton's bisection method to find new alpha + newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, variable, + min, initial_check_nonnegative_newton_idp, + final_check_nonnegative_newton_idp, inverse_jacobian, + dt, equations, dg, cache, limiter) + + return nothing +end + +############################################################################### +# Newton-bisection method + @inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable, min_or_max, initial_check, final_check, inverse_jacobian, dt, equations, dg, cache, diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index f095c97b19e..6875259e788 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -475,6 +475,32 @@ end end end +@trixi_testset "elixir_euler_source_terms_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_sc_subcell.jl"), + l2=[ + 0.00816012943805366, + 0.008658255997419217, + 0.00935190284719103, + 0.027757019482291357, + ], + linf=[ + 0.02722563220635177, + 0.040734034380730755, + 0.03819407626402338, + 0.08080649141597318, + ], + tspan=(0.0, 0.5)) + # 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)) < 10000 + end +end + @trixi_testset "elixir_euler_source_terms_waving_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_waving_flag.jl"), @@ -525,6 +551,33 @@ end end end +# Changing because of disabled bar states +@trixi_testset "elixir_euler_free_stream_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_free_stream_sc_subcell.jl"), + l2=[ + 1.4663777294625118e-15, + 2.320054900530864e-14, + 3.487555722563465e-14, + 2.008802099296406e-14, + ], + linf=[ + 2.3092638912203256e-14, + 2.0623780461193064e-13, + 2.6795232699328153e-13, + 2.362554596402333e-13, + ], + atol=7.0e-13) + # 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)) < 10000 + end +end + @trixi_testset "elixir_euler_free_stream.jl with FluxRotated(flux_lax_friedrichs)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), surface_flux=FluxRotated(flux_lax_friedrichs), From c1af151d6024565828202da0b83725f7decaffa3 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 15 May 2024 11:03:13 +0200 Subject: [PATCH 02/19] Fix non-periodic computation of bounds --- .../dgsem_structured/subcell_limiters_2d.jl | 60 +++++-------------- 1 file changed, 16 insertions(+), 44 deletions(-) diff --git a/src/solvers/dgsem_structured/subcell_limiters_2d.jl b/src/solvers/dgsem_structured/subcell_limiters_2d.jl index 83a979a16a8..99b86a31a67 100644 --- a/src/solvers/dgsem_structured/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_structured/subcell_limiters_2d.jl @@ -57,11 +57,8 @@ 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) - 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, cache, t, - boundary_conditions[1], Ja1, 1, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[1], + cache, t, equations, dg, 1, j, element) var_outer = u_outer[variable] @@ -73,12 +70,8 @@ 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) - 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, cache, t, - boundary_conditions[2], Ja1, 2, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[2], + cache, t, equations, dg, nnodes(dg), j, element) var_outer = u_outer[variable] @@ -94,11 +87,8 @@ 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) - 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, cache, t, - boundary_conditions[3], Ja2, 3, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[3], + cache, t, equations, dg, i, 1, element) var_outer = u_outer[variable] @@ -110,12 +100,8 @@ 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) - 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, cache, t, - boundary_conditions[4], Ja2, 4, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[4], + cache, t, equations, dg, i, nnodes(dg), element) var_outer = u_outer[variable] @@ -178,11 +164,8 @@ 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) - 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, cache, t, - boundary_conditions[1], Ja1, 1, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[1], + cache, t, equations, dg, 1, j, element) var_outer = variable(u_outer, equations) @@ -193,12 +176,8 @@ 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) - 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, cache, t, - boundary_conditions[2], Ja1, 2, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[2], + cache, t, equations, dg, nnodes(dg), j, element) var_outer = variable(u_outer, equations) @@ -213,11 +192,8 @@ 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) - 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, cache, t, - boundary_conditions[3], Ja2, 3, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[3], + cache, t, equations, dg, i, 1, element) var_outer = variable(u_outer, equations) @@ -228,12 +204,8 @@ 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) - 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, cache, t, - boundary_conditions[4], Ja2, 4, - mesh, equations, dg, + u_outer = get_boundary_outer_state(boundary_conditions[4], + cache, t, equations, dg, i, nnodes(dg), element) var_outer = variable(u_outer, equations) From 2ddd78aa0451a00c5929178596ebc711bd19af05 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 15 May 2024 13:23:34 +0200 Subject: [PATCH 03/19] Use local limiting and nonperiodic domain in source terms elixir --- .../elixir_euler_source_terms_sc_subcell.jl | 16 +++++++++++----- test/test_structured_2d.jl | 16 ++++++++-------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl index 97c541266fa..ae413828861 100644 --- a/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl +++ b/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl @@ -10,15 +10,20 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition = initial_condition_convergence_test source_terms = source_terms_convergence_test +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = (x_neg = boundary_condition, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition) + # 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)]) + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) @@ -32,9 +37,10 @@ f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) mapping = Trixi.transfinite_mapping((f1, f2, f3, f4)) cells_per_dimension = (16, 16) -mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = true) +mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions, source_terms = source_terms) ############################################################################### @@ -55,7 +61,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.5) +stepsize_callback = StepsizeCallback(cfl = 0.8) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 6875259e788..ff27f683a40 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -479,16 +479,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_sc_subcell.jl"), l2=[ - 0.00816012943805366, - 0.008658255997419217, - 0.00935190284719103, - 0.027757019482291357, + 7.939813971944535e-5, + 3.6805976512737965e-5, + 7.50142174035037e-5, + 0.0001576144969663713, ], linf=[ - 0.02722563220635177, - 0.040734034380730755, - 0.03819407626402338, - 0.08080649141597318, + 0.0007421863054295486, + 0.00037635704053817776, + 0.0007278708058917616, + 0.001550366737408826, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations From 98c2144ac49bfcbc741c9f782a43c8a79dfb2412 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 15 May 2024 13:24:18 +0200 Subject: [PATCH 04/19] Use local limiting in free stream elixir --- .../elixir_euler_free_stream_sc_subcell.jl | 9 +++++---- test/test_structured_2d.jl | 16 ++++++++-------- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl index a0ddcb64eb8..3f15391049f 100644 --- a/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl +++ b/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl @@ -14,8 +14,9 @@ volume_flux = flux_ranocha polydeg = 3 basis = LobattoLegendreBasis(polydeg) limiter_idp = SubcellLimiterIDP(equations, basis; - positivity_variables_cons = ["rho"], - positivity_variables_nonlinear = [pressure]) + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, + min)]) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, @@ -61,12 +62,12 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval = 10000, +save_solution = SaveSolutionCallback(interval = 100, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.9) +stepsize_callback = StepsizeCallback(cfl = 0.7) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index ff27f683a40..569deffcf43 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -556,16 +556,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_sc_subcell.jl"), l2=[ - 1.4663777294625118e-15, - 2.320054900530864e-14, - 3.487555722563465e-14, - 2.008802099296406e-14, + 2.5709100186298863e-16, + 2.184017068634721e-15, + 2.1781991538183746e-15, + 4.935868571076398e-15, ], linf=[ - 2.3092638912203256e-14, - 2.0623780461193064e-13, - 2.6795232699328153e-13, - 2.362554596402333e-13, + 4.218847493575595e-15, + 2.7602919949742954e-14, + 3.533284775869561e-14, + 5.861977570020827e-14, ], atol=7.0e-13) # Ensure that we do not have excessive memory allocations From ccb833e4e0187abd663ac43776a1bd757dc81589 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 15 May 2024 14:04:49 +0200 Subject: [PATCH 05/19] Remove not needed lines --- src/solvers/dgsem_structured/subcell_limiters_2d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/solvers/dgsem_structured/subcell_limiters_2d.jl b/src/solvers/dgsem_structured/subcell_limiters_2d.jl index 99b86a31a67..5b65a475062 100644 --- a/src/solvers/dgsem_structured/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_structured/subcell_limiters_2d.jl @@ -9,7 +9,6 @@ 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) @@ -120,7 +119,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) # Get neighboring element ids From 003d43c6ce0ebed15ec9a8f6044cd8b7ba9c9e34 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 15 May 2024 14:25:32 +0200 Subject: [PATCH 06/19] Remove P4estMesh --- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 7dd0abe3cad..45791b1a7b0 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -492,8 +492,7 @@ end end @inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, - mesh::Union{StructuredMesh{2}, - P4estMesh{2}}, + mesh::StructuredMesh{2}, variable) _, _, dg, cache = mesh_equations_solver_cache(semi) From e918af0f5ecb82672651e0ce779c7be5bd550d5f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 15 May 2024 15:19:03 +0200 Subject: [PATCH 07/19] Add non-periodic tests with local bounds --- .../elixir_euler_source_terms_sc_subcell.jl | 7 +++- test/test_structured_2d.jl | 33 ++++++++++++++++++- 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl index ae413828861..0efe47ea16c 100644 --- a/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl +++ b/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl @@ -23,7 +23,12 @@ 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 = [], + local_onesided_variables_nonlinear = []) +# Variables for local limiting (`local_twosided_variables_cons` and +# `local_onesided_variables_nonlinear`) are overwritten and used in the tests. + volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 569deffcf43..46c52561740 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -475,7 +475,7 @@ end end end -@trixi_testset "elixir_euler_source_terms_sc_subcell.jl" begin +@trixi_testset "elixir_euler_source_terms_sc_subcell.jl (global bounds)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_sc_subcell.jl"), l2=[ @@ -501,6 +501,37 @@ end end end +@trixi_testset "elixir_euler_source_terms_sc_subcell.jl (local bounds)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_sc_subcell.jl"), + positivity_variables_cons = [], + positivity_variables_nonlinear = [], + local_twosided_variables_cons=["rho"], + local_onesided_variables_nonlinear=[(Trixi.entropy_guermond_etal, min)], + cfl=0.5, + l2=[ + 0.007788373240296921, + 0.006564114774853376, + 0.008411575110395622, + 0.023360113470962406, + ], + linf=[ + 0.033816852778793205, + 0.03938807720446702, + 0.044093017726981376, + 0.08006777547232469, + ], + tspan=(0.0, 0.5)) + # 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)) < 10000 + end +end + @trixi_testset "elixir_euler_source_terms_waving_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_waving_flag.jl"), From 5a7b8ceb28d2a176fc21a2984de07f19b6889f20 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 15 May 2024 15:19:47 +0200 Subject: [PATCH 08/19] fmt --- test/test_structured_2d.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 46c52561740..d6dce99e764 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -504,10 +504,11 @@ end @trixi_testset "elixir_euler_source_terms_sc_subcell.jl (local bounds)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_sc_subcell.jl"), - positivity_variables_cons = [], - positivity_variables_nonlinear = [], + positivity_variables_cons=[], + positivity_variables_nonlinear=[], local_twosided_variables_cons=["rho"], - local_onesided_variables_nonlinear=[(Trixi.entropy_guermond_etal, min)], + local_onesided_variables_nonlinear=[(Trixi.entropy_guermond_etal, + min)], cfl=0.5, l2=[ 0.007788373240296921, From 632a9342ad440098a62d9293b6356c67efff80e5 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 May 2024 09:42:11 +0200 Subject: [PATCH 09/19] Fix test --- test/test_structured_2d.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index d6dce99e764..caec71e954f 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -511,16 +511,16 @@ end min)], cfl=0.5, l2=[ - 0.007788373240296921, - 0.006564114774853376, - 0.008411575110395622, - 0.023360113470962406, + 0.007788374804566746, + 0.006564110929038488, + 0.00841157867996223, + 0.023360113888099734, ], linf=[ - 0.033816852778793205, - 0.03938807720446702, + 0.033816854019882214, + 0.03938807953377843, 0.044093017726981376, - 0.08006777547232469, + 0.08006778793724179, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations From b26a8f24def4cc593befa52108a7f7bfc19709e3 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 May 2024 14:37:58 +0200 Subject: [PATCH 10/19] Use `get_inverse_jacobian` instead of dispatching all routines --- src/solvers/dgsem_structured/dg.jl | 8 + src/solvers/dgsem_tree/dg.jl | 5 + src/solvers/dgsem_tree/subcell_limiters_2d.jl | 361 ++++++------------ 3 files changed, 130 insertions(+), 244 deletions(-) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index ecb2485c68b..5617ae90e3f 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -76,6 +76,14 @@ end end end +@inline function get_inverse_jacobian(inverse_jacobian, + mesh::Union{StructuredMesh, StructuredMeshView, + UnstructuredMesh2D, P4estMesh, + T8codeMesh}, + indices...) + return inverse_jacobian[indices...] +end + include("containers.jl") include("dg_1d.jl") include("dg_2d.jl") diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index ef9a42b4c1a..1c754b172e1 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -42,6 +42,11 @@ function volume_jacobian(element, mesh::TreeMesh, cache) return inv(cache.elements.inverse_jacobian[element])^ndims(mesh) end +@inline function get_inverse_jacobian(inverse_jacobian, mesh::TreeMesh, + indices...) + return inverse_jacobian[last(indices)] +end + # Indicators used for shock-capturing and AMR include("indicators.jl") include("indicators_1d.jl") diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 45791b1a7b0..fe053557503 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -275,18 +275,17 @@ end # Local two-sided limiting of conservative variables @inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi) - mesh, _, _, _ = mesh_equations_solver_cache(semi) - for variable in limiter.local_twosided_variables_cons - idp_local_twosided!(alpha, limiter, u, t, dt, semi, mesh, variable) + idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) end return nothing end -@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, mesh::TreeMesh{2}, - variable) - _, _, dg, cache = mesh_equations_solver_cache(semi) +@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) + mesh, _, dg, cache = mesh_equations_solver_cache(semi) + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + (; inverse_weights) = dg.basis (; variable_bounds) = limiter.cache.subcell_limiter_coefficients variable_string = string(variable) @@ -297,128 +296,69 @@ end @threaded for element in eachelement(dg, semi.cache) inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) - idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache, - variable, var_min, var_max, i, j, element) - end - end - - return nothing -end - -@inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, - mesh::StructuredMesh{2}, - variable) - _, _, dg, cache = mesh_equations_solver_cache(semi) - - (; variable_bounds) = limiter.cache.subcell_limiter_coefficients - variable_string = string(variable) - var_min = variable_bounds[Symbol(variable_string, "_min")] - var_max = variable_bounds[Symbol(variable_string, "_max")] - calc_bounds_twosided!(var_min, var_max, variable, u, t, semi) - - @threaded for element in eachelement(dg, cache) - for j in eachnode(dg), i in eachnode(dg) - inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] - idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache, - variable, var_min, var_max, i, j, element) + inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, element) + var = u[variable, i, j, element] + # Real Zalesak type limiter + # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" + # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" + # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is + # for each interface, not each node + + Qp = max(0, (var_max[i, j, element] - var) / dt) + Qm = min(0, (var_min[i, j, element] - var) / dt) + + # Calculate Pp and Pm + # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. + val_flux1_local = inverse_weights[i] * + antidiffusive_flux1_R[variable, i, j, element] + val_flux1_local_ip1 = -inverse_weights[i] * + antidiffusive_flux1_L[variable, i + 1, j, element] + val_flux2_local = inverse_weights[j] * + antidiffusive_flux2_R[variable, i, j, element] + val_flux2_local_jp1 = -inverse_weights[j] * + antidiffusive_flux2_L[variable, i, j + 1, element] + + Pp = max(0, val_flux1_local) + max(0, val_flux1_local_ip1) + + max(0, val_flux2_local) + max(0, val_flux2_local_jp1) + Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + + min(0, val_flux2_local) + min(0, val_flux2_local_jp1) + + Qp = max(0, (var_max[i, j, element] - var) / dt) + Qm = min(0, (var_min[i, j, element] - var) / dt) + + Pp = inverse_jacobian * Pp + Pm = inverse_jacobian * Pm + + # Compute blending coefficient avoiding division by zero + # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) + Qp = abs(Qp) / + (abs(Pp) + eps(typeof(Qp)) * 100 * abs(var_max[i, j, element])) + Qm = abs(Qm) / + (abs(Pm) + eps(typeof(Qm)) * 100 * abs(var_max[i, j, element])) + + # Calculate alpha at nodes + alpha[i, j, element] = max(alpha[i, j, element], 1 - min(1, Qp, Qm)) end end return nothing end -# Function barrier to dispatch outer function by mesh type -@inline function idp_local_twosided_inner!(alpha, inverse_jacobian, u, dt, dg, cache, - variable, var_min, var_max, i, j, element) - (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes - (; inverse_weights) = dg.basis - - var = u[variable, i, j, element] - # Real Zalesak type limiter - # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" - # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" - # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is - # for each interface, not each node - - Qp = max(0, (var_max[i, j, element] - var) / dt) - Qm = min(0, (var_min[i, j, element] - var) / dt) - - # Calculate Pp and Pm - # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. - val_flux1_local = inverse_weights[i] * - antidiffusive_flux1_R[variable, i, j, element] - val_flux1_local_ip1 = -inverse_weights[i] * - antidiffusive_flux1_L[variable, i + 1, j, element] - val_flux2_local = inverse_weights[j] * - antidiffusive_flux2_R[variable, i, j, element] - val_flux2_local_jp1 = -inverse_weights[j] * - antidiffusive_flux2_L[variable, i, j + 1, element] - - Pp = max(0, val_flux1_local) + max(0, val_flux1_local_ip1) + - max(0, val_flux2_local) + max(0, val_flux2_local_jp1) - Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + - min(0, val_flux2_local) + min(0, val_flux2_local_jp1) - - Qp = max(0, (var_max[i, j, element] - var) / dt) - Qm = min(0, (var_min[i, j, element] - var) / dt) - - Pp = inverse_jacobian * Pp - Pm = inverse_jacobian * Pm - - # Compute blending coefficient avoiding division by zero - # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) - Qp = abs(Qp) / - (abs(Pp) + eps(typeof(Qp)) * 100 * abs(var_max[i, j, element])) - Qm = abs(Qm) / - (abs(Pm) + eps(typeof(Qm)) * 100 * abs(var_max[i, j, element])) - - # Calculate alpha at nodes - alpha[i, j, element] = max(alpha[i, j, element], 1 - min(1, Qp, Qm)) - - return nothing -end - ############################################################################## # Local one-sided limiting of nonlinear variables @inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi) - mesh, _, _, _ = mesh_equations_solver_cache(semi) - for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear - idp_local_onesided!(alpha, limiter, u, t, dt, semi, mesh, - variable, min_or_max) - end - - return nothing -end - -@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, mesh::TreeMesh{2}, - variable, min_or_max) - _, equations, dg, cache = mesh_equations_solver_cache(semi) - (; variable_bounds) = limiter.cache.subcell_limiter_coefficients - var_minmax = variable_bounds[Symbol(string(variable), "_", string(min_or_max))] - calc_bounds_onesided!(var_minmax, min_or_max, variable, u, t, semi) - - # Perform Newton's bisection method to find new alpha - @threaded for element in eachelement(dg, cache) - inverse_jacobian = cache.elements.inverse_jacobian[element] - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - newton_loops_alpha!(alpha, var_minmax[i, j, element], u_local, - i, j, element, variable, min_or_max, - initial_check_local_onesided_newton_idp, - final_check_local_onesided_newton_idp, inverse_jacobian, - dt, equations, dg, cache, limiter) - end + idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, min_or_max) end return nothing end @inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, - mesh::StructuredMesh{2}, variable, min_or_max) - _, equations, dg, cache = mesh_equations_solver_cache(semi) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_minmax = variable_bounds[Symbol(string(variable), "_", string(min_or_max))] calc_bounds_onesided!(var_minmax, min_or_max, variable, u, t, semi) @@ -426,7 +366,8 @@ end # Perform Newton's bisection method to find new alpha @threaded for element in eachelement(dg, cache) for j in eachnode(dg), i in eachnode(dg) - inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] + inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, element) u_local = get_node_vars(u, equations, dg, i, j, element) newton_loops_alpha!(alpha, var_minmax[i, j, element], u_local, i, j, element, variable, min_or_max, @@ -443,8 +384,6 @@ end # Global positivity limiting @inline function idp_positivity!(alpha, limiter, u, dt, semi) - mesh, _, _, _ = mesh_equations_solver_cache(semi) - # Conservative variables for variable in limiter.positivity_variables_cons @trixi_timeit timer() "conservative variables" idp_positivity_conservative!(alpha, @@ -452,7 +391,6 @@ end u, dt, semi, - mesh, variable) end @@ -462,7 +400,6 @@ end limiter, u, dt, semi, - mesh, variable) end @@ -472,164 +409,100 @@ end ############################################################################### # Global positivity limiting of conservative variables -@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, - mesh::TreeMesh{2}, variable) - _, _, dg, cache = mesh_equations_solver_cache(semi) +@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, variable) + mesh, _, dg, cache = mesh_equations_solver_cache(semi) + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes + (; inverse_weights) = dg.basis + (; positivity_correction_factor) = limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_min = variable_bounds[Symbol(string(variable), "_min")] @threaded for element in eachelement(dg, semi.cache) - inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) - idp_positivity_conservative_inner!(alpha, inverse_jacobian, limiter, u, dt, - dg, cache, variable, var_min, - i, j, element) - end - end - - return nothing -end - -@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, - mesh::StructuredMesh{2}, - variable) - _, _, dg, cache = mesh_equations_solver_cache(semi) - - (; variable_bounds) = limiter.cache.subcell_limiter_coefficients - var_min = variable_bounds[Symbol(string(variable), "_min")] + inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, element) + var = u[variable, i, j, element] + if var < 0 + error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") + end - @threaded for element in eachelement(dg, semi.cache) - for j in eachnode(dg), i in eachnode(dg) - inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] - idp_positivity_conservative_inner!(alpha, inverse_jacobian, limiter, u, dt, - dg, cache, variable, var_min, - i, j, element) + # Compute bound + if limiter.local_twosided && + variable in limiter.local_twosided_variables_cons && + var_min[i, j, element] >= positivity_correction_factor * var + # Local limiting is more restrictive that positivity limiting + # => Skip positivity limiting for this node + continue + end + var_min[i, j, element] = positivity_correction_factor * var + + # Real one-sided Zalesak-type limiter + # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" + # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" + # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is + # for each interface, not each node + Qm = min(0, (var_min[i, j, element] - var) / dt) + + # Calculate Pm + # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. + val_flux1_local = inverse_weights[i] * + antidiffusive_flux1_R[variable, i, j, element] + val_flux1_local_ip1 = -inverse_weights[i] * + antidiffusive_flux1_L[variable, i + 1, j, element] + val_flux2_local = inverse_weights[j] * + antidiffusive_flux2_R[variable, i, j, element] + val_flux2_local_jp1 = -inverse_weights[j] * + antidiffusive_flux2_L[variable, i, j + 1, element] + + Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + + min(0, val_flux2_local) + min(0, val_flux2_local_jp1) + Pm = inverse_jacobian * Pm + + # Compute blending coefficient avoiding division by zero + # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) + Qm = abs(Qm) / (abs(Pm) + eps(typeof(Qm)) * 100) + + # Calculate alpha + alpha[i, j, element] = max(alpha[i, j, element], 1 - Qm) end end return nothing end -# Function barrier to dispatch outer function by mesh type -@inline function idp_positivity_conservative_inner!(alpha, inverse_jacobian, limiter, u, - dt, dg, cache, variable, var_min, - i, j, element) - (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes - (; inverse_weights) = dg.basis - (; positivity_correction_factor) = limiter - - var = u[variable, i, j, element] - if var < 0 - error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") - end - - # Compute bound - if limiter.local_twosided && - variable in limiter.local_twosided_variables_cons && - var_min[i, j, element] >= positivity_correction_factor * var - # Local limiting is more restrictive that positivity limiting - # => Skip positivity limiting for this node - return nothing - end - var_min[i, j, element] = positivity_correction_factor * var - - # Real one-sided Zalesak-type limiter - # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" - # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics" - # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is - # for each interface, not each node - Qm = min(0, (var_min[i, j, element] - var) / dt) - - # Calculate Pm - # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. - val_flux1_local = inverse_weights[i] * - antidiffusive_flux1_R[variable, i, j, element] - val_flux1_local_ip1 = -inverse_weights[i] * - antidiffusive_flux1_L[variable, i + 1, j, element] - val_flux2_local = inverse_weights[j] * - antidiffusive_flux2_R[variable, i, j, element] - val_flux2_local_jp1 = -inverse_weights[j] * - antidiffusive_flux2_L[variable, i, j + 1, element] - - Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + - min(0, val_flux2_local) + min(0, val_flux2_local_jp1) - Pm = inverse_jacobian * Pm - - # Compute blending coefficient avoiding division by zero - # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8)) - Qm = abs(Qm) / (abs(Pm) + eps(typeof(Qm)) * 100) - - # Calculate alpha - alpha[i, j, element] = max(alpha[i, j, element], 1 - Qm) - - return nothing -end - ############################################################################### # Global positivity limiting of nonlinear variables -@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, - mesh::TreeMesh{2}, variable) - _, _, dg, cache = mesh_equations_solver_cache(semi) +@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) + mesh, equations, dg, cache = mesh_equations_solver_cache(semi) (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_min = variable_bounds[Symbol(string(variable), "_min")] @threaded for element in eachelement(dg, semi.cache) - inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) - idp_positivity_nonlinear_inner!(alpha, inverse_jacobian, limiter, u, dt, - semi, dg, cache, variable, var_min, - i, j, element) - end - end - - return nothing -end - -@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, - mesh::StructuredMesh{2}, - variable) - _, _, dg, cache = mesh_equations_solver_cache(semi) - - (; variable_bounds) = limiter.cache.subcell_limiter_coefficients - var_min = variable_bounds[Symbol(string(variable), "_min")] + inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, element) + u_local = get_node_vars(u, equations, dg, i, j, element) + var = variable(u_local, equations) + if var < 0 + error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") + end + var_min[i, j, element] = limiter.positivity_correction_factor * var - @threaded for element in eachelement(dg, semi.cache) - for j in eachnode(dg), i in eachnode(dg) - inverse_jacobian = cache.elements.inverse_jacobian[i, j, element] - idp_positivity_nonlinear_inner!(alpha, inverse_jacobian, limiter, u, dt, - semi, dg, cache, variable, var_min, - i, j, element) + # Perform Newton's bisection method to find new alpha + newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, + variable, + min, initial_check_nonnegative_newton_idp, + final_check_nonnegative_newton_idp, inverse_jacobian, + dt, equations, dg, cache, limiter) end end return nothing end -# Function barrier to dispatch outer function by mesh type -@inline function idp_positivity_nonlinear_inner!(alpha, inverse_jacobian, limiter, u, - dt, semi, dg, cache, variable, var_min, - i, j, element) - _, equations, _, _ = mesh_equations_solver_cache(semi) - - u_local = get_node_vars(u, equations, dg, i, j, element) - var = variable(u_local, equations) - if var < 0 - error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") - end - var_min[i, j, element] = limiter.positivity_correction_factor * var - - # Perform Newton's bisection method to find new alpha - newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, variable, - min, initial_check_nonnegative_newton_idp, - final_check_nonnegative_newton_idp, inverse_jacobian, - dt, equations, dg, cache, limiter) - - return nothing -end - ############################################################################### # Newton-bisection method From 8c5047a8fd50dafd10481fe41cc30f771f6833f1 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 May 2024 14:43:20 +0200 Subject: [PATCH 11/19] Simplify `perform_idp_correction!` --- .../subcell_limiter_idp_correction_2d.jl | 76 +++++++------------ 1 file changed, 28 insertions(+), 48 deletions(-) diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index 9c4a31a515c..2b74abdfb40 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -5,61 +5,41 @@ @muladd begin #! format: noindent -function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache) - @threaded for element in eachelement(dg, cache) - # Sign switch as in apply_jacobian! - inverse_jacobian = -cache.elements.inverse_jacobian[element] - - for j in eachnode(dg), i in eachnode(dg) - perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, cache, - i, j, element) - end - end - - return nothing -end +function perform_idp_correction!(u, dt, mesh, equations, dg, cache) + @unpack inverse_weights = dg.basis + @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes + @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients -function perform_idp_correction!(u, dt, mesh::StructuredMesh{2}, equations, dg, cache) @threaded for element in eachelement(dg, cache) for j in eachnode(dg), i in eachnode(dg) # Sign switch as in apply_jacobian! - inverse_jacobian = -cache.elements.inverse_jacobian[i, j, element] - - perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, cache, - i, j, element) + inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian, + mesh, i, j, element) + + # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1} + alpha_flux1 = (1 - alpha1[i, j, element]) * + get_node_vars(antidiffusive_flux1_R, equations, dg, + i, j, element) + alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) * + get_node_vars(antidiffusive_flux1_L, equations, dg, + i + 1, j, element) + alpha_flux2 = (1 - alpha2[i, j, element]) * + get_node_vars(antidiffusive_flux2_R, equations, dg, + i, j, element) + alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) * + get_node_vars(antidiffusive_flux2_L, equations, dg, + i, j + 1, element) + + for v in eachvariable(equations) + u[v, i, j, element] += dt * inverse_jacobian * + (inverse_weights[i] * + (alpha_flux1_ip1[v] - alpha_flux1[v]) + + inverse_weights[j] * + (alpha_flux2_jp1[v] - alpha_flux2[v])) + end end end return nothing end - -# Function barrier to dispatch outer function by mesh type -@inline function perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, - cache, i, j, element) - (; inverse_weights) = dg.basis - (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes - (; alpha1, alpha2) = dg.volume_integral.limiter.cache.subcell_limiter_coefficients - - # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1} - alpha_flux1 = (1 - alpha1[i, j, element]) * - get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, element) - alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) * - get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j, - element) - alpha_flux2 = (1 - alpha2[i, j, element]) * - get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, element) - alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) * - get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1, - element) - - for v in eachvariable(equations) - u[v, i, j, element] += dt * inverse_jacobian * - (inverse_weights[i] * - (alpha_flux1_ip1[v] - alpha_flux1[v]) + - inverse_weights[j] * - (alpha_flux2_jp1[v] - alpha_flux2[v])) - end - - return nothing -end end # @muladd From 80dd9bf82aed6e934c7c4a56b6721729db90c452 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 May 2024 14:47:49 +0200 Subject: [PATCH 12/19] Revert stuff --- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index fe053557503..4095f0853f9 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -294,7 +294,6 @@ end calc_bounds_twosided!(var_min, var_max, variable, u, t, semi) @threaded for element in eachelement(dg, semi.cache) - inverse_jacobian = cache.elements.inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, mesh, i, j, element) @@ -476,6 +475,7 @@ end @inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) + (; positivity_correction_factor) = limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients var_min = variable_bounds[Symbol(string(variable), "_min")] @@ -484,17 +484,18 @@ end for j in eachnode(dg), i in eachnode(dg) inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian, mesh, i, j, element) + + # Compute bound u_local = get_node_vars(u, equations, dg, i, j, element) var = variable(u_local, equations) if var < 0 - error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.") + error("Safe low-order method produces negative value for variable $variable. Try a smaller time step.") end - var_min[i, j, element] = limiter.positivity_correction_factor * var + var_min[i, j, element] = positivity_correction_factor * var # Perform Newton's bisection method to find new alpha newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element, - variable, - min, initial_check_nonnegative_newton_idp, + variable, min, initial_check_nonnegative_newton_idp, final_check_nonnegative_newton_idp, inverse_jacobian, dt, equations, dg, cache, limiter) end From 7102fbdd9e2d88df8c704f56d0239f0fb7f972d6 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 16 May 2024 15:03:38 +0200 Subject: [PATCH 13/19] Remove free stream elixir --- .../elixir_euler_free_stream_sc_subcell.jl | 85 ------------------- test/test_structured_2d.jl | 27 ------ 2 files changed, 112 deletions(-) delete mode 100644 examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl diff --git a/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl deleted file mode 100644 index 3f15391049f..00000000000 --- a/examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl +++ /dev/null @@ -1,85 +0,0 @@ - -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; - local_twosided_variables_cons = ["rho"], - local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, - min)]) - -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 - -cells_per_dimension = (16, 16) -mesh = StructuredMesh(cells_per_dimension, mapping, 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 = 100, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim) - -stepsize_callback = StepsizeCallback(cfl = 0.7) - -callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, - stepsize_callback, - save_solution) - -############################################################################### -# 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 - save_everystep = false, callback = callbacks); -summary_callback() # print the timer summary diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index caec71e954f..3bcf86a00fa 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -583,33 +583,6 @@ end end end -# Changing because of disabled bar states -@trixi_testset "elixir_euler_free_stream_sc_subcell.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_free_stream_sc_subcell.jl"), - l2=[ - 2.5709100186298863e-16, - 2.184017068634721e-15, - 2.1781991538183746e-15, - 4.935868571076398e-15, - ], - linf=[ - 4.218847493575595e-15, - 2.7602919949742954e-14, - 3.533284775869561e-14, - 5.861977570020827e-14, - ], - atol=7.0e-13) - # 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)) < 10000 - end -end - @trixi_testset "elixir_euler_free_stream.jl with FluxRotated(flux_lax_friedrichs)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), surface_flux=FluxRotated(flux_lax_friedrichs), From dc6d5547b0bd6ec4da6c4a063832530123952964 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 22 May 2024 15:20:57 +0200 Subject: [PATCH 14/19] Use sedov blast instead of source term setup; add news --- NEWS.md | 1 + ...ixir_euler_sedov_blast_wave_sc_subcell.jl} | 76 +++++++++++++------ test/test_structured_2d.jl | 50 ++++++------ 3 files changed, 78 insertions(+), 49 deletions(-) rename examples/structured_2d_dgsem/{elixir_euler_source_terms_sc_subcell.jl => elixir_euler_sedov_blast_wave_sc_subcell.jl} (55%) diff --git a/NEWS.md b/NEWS.md index e2902229f71..5de652242f6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,7 @@ for human readability. - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients ([#1812]). - Optional tuple parameter for `GlmSpeedCallback` called `semi_indices` to specify for which semidiscretization of a `SemidiscretizationCoupled` we need to update the GLM speed ([#1835]). - Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh` ([#1792]). +- Add subcell limiting support for `StructuredMesh` ([#1946]). ## Changes when updating to v0.7 from v0.6.x diff --git a/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl similarity index 55% rename from examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl rename to examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 0efe47ea16c..5c11a7d15a7 100644 --- a/examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl +++ b/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -1,14 +1,41 @@ - using OrdinaryDiffEq using Trixi ############################################################################### # semidiscretization of the compressible Euler equations - -equations = CompressibleEulerEquations2D(1.4) - -initial_condition = initial_condition_convergence_test -source_terms = source_terms_convergence_test +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + 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) + # r0 = 0.5 # = more reasonable setup + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) + p0_outer = 1.0e-5 # = true Sedov setup + # p0_outer = 1.0e-3 # = more reasonable 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 boundary_condition = BoundaryConditionDirichlet(initial_condition) boundary_conditions = (x_neg = boundary_condition, @@ -16,42 +43,45 @@ boundary_conditions = (x_neg = boundary_condition, y_neg = boundary_condition, y_pos = boundary_condition) -# Get the DG approximation space 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], - local_twosided_variables_cons = [], - local_onesided_variables_nonlinear = []) -# Variables for local limiting (`local_twosided_variables_cons` and -# `local_onesided_variables_nonlinear`) are overwritten and used in the tests. + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, + min)], + max_iterations_newton = 40, # Default value of 10 iterations is too low to fulfill bounds. + positivity_variables_cons = [], + positivity_variables_nonlinear = []) +# Variables for global limiting (`positivity_variables_cons` and +# `positivity_variables_nonlinear`) are overwritten and used in the tests. volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) -# Waving flag -f1(s) = SVector(-1.0, s - 1.0) -f2(s) = SVector(1.0, s + 1.0) -f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) -f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) -mapping = Trixi.transfinite_mapping((f1, f2, f3, f4)) +# Get the curved quad mesh from a mapping function +# Mapping as described in https://arxiv.org/abs/2012.12040 +function mapping(xi, eta) + y = eta + 0.125 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta)) + + x = xi + 0.125 * (cos(0.5 * pi * xi) * cos(2 * pi * y)) + + return SVector(x, y) +end cells_per_dimension = (16, 16) mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - boundary_conditions = boundary_conditions, - source_terms = source_terms) + boundary_conditions = boundary_conditions) ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 2.0) +tspan = (0.0, 3.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -66,7 +96,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.8) +stepsize_callback = StepsizeCallback(cfl = 0.7) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 3bcf86a00fa..01a1185e70f 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -475,20 +475,20 @@ end end end -@trixi_testset "elixir_euler_source_terms_sc_subcell.jl (global bounds)" begin +@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (local bounds)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_source_terms_sc_subcell.jl"), + "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ - 7.939813971944535e-5, - 3.6805976512737965e-5, - 7.50142174035037e-5, - 0.0001576144969663713, + 0.6337774834710513, + 0.30377119245852724, + 0.3111372568571772, + 1.2976221893997268, ], linf=[ - 0.0007421863054295486, - 0.00037635704053817776, - 0.0007278708058917616, - 0.001550366737408826, + 2.2064877103138207, + 1.541067099687334, + 1.5487587769900337, + 6.271271639873466, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations @@ -501,26 +501,24 @@ end end end -@trixi_testset "elixir_euler_source_terms_sc_subcell.jl (local bounds)" begin +@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (global bounds)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_source_terms_sc_subcell.jl"), - positivity_variables_cons=[], - positivity_variables_nonlinear=[], - local_twosided_variables_cons=["rho"], - local_onesided_variables_nonlinear=[(Trixi.entropy_guermond_etal, - min)], - cfl=0.5, + "elixir_euler_sedov_blast_wave_sc_subcell.jl"), + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure], + local_twosided_variables_cons=[], + local_onesided_variables_nonlinear=[], l2=[ - 0.007788374804566746, - 0.006564110929038488, - 0.00841157867996223, - 0.023360113888099734, + 0.7869912572385168, + 0.39170886758882073, + 0.39613257454431977, + 1.2951760266455101, ], linf=[ - 0.033816854019882214, - 0.03938807953377843, - 0.044093017726981376, - 0.08006778793724179, + 5.156044534854053, + 3.6261667239538986, + 3.1807681416546085, + 6.3028422220287235, ], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations From ee86e27aa9d8c6b635445ebf8da85f082746a7b7 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 22 May 2024 15:43:30 +0200 Subject: [PATCH 15/19] Update dispatching for mesh types --- src/callbacks_stage/subcell_limiter_idp_correction_2d.jl | 2 +- src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl | 8 ++++---- test/test_structured_2d.jl | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index 2b74abdfb40..04d7bf30821 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -5,7 +5,7 @@ @muladd begin #! format: noindent -function perform_idp_correction!(u, dt, mesh, equations, dg, cache) +function perform_idp_correction!(u, dt, mesh::AbstractMesh{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 @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index ba029667c17..ab5182fb389 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::AbstractMesh{2}, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DG, uEltype) cache = create_cache(mesh, equations, @@ -57,7 +57,7 @@ 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, StructuredMesh}, nonconservative_terms, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache) @@ -392,7 +392,7 @@ 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, + u, mesh::AbstractMesh{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 +430,7 @@ 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, + u, mesh::AbstractMesh{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 diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 01a1185e70f..16670db710a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -504,8 +504,8 @@ end @trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (global bounds)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), - positivity_variables_cons = ["rho"], - positivity_variables_nonlinear = [pressure], + positivity_variables_cons=["rho"], + positivity_variables_nonlinear=[pressure], local_twosided_variables_cons=[], local_onesided_variables_nonlinear=[], l2=[ From 7f8f58d27134f209157d691948828abcc1a4516b Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 23 May 2024 15:00:03 +0200 Subject: [PATCH 16/19] Move new tests within test file --- test/test_structured_2d.jl | 112 ++++++++++++++++++------------------- 1 file changed, 56 insertions(+), 56 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 16670db710a..4b3aa5c87e4 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -475,62 +475,6 @@ end end end -@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (local bounds)" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_sedov_blast_wave_sc_subcell.jl"), - l2=[ - 0.6337774834710513, - 0.30377119245852724, - 0.3111372568571772, - 1.2976221893997268, - ], - linf=[ - 2.2064877103138207, - 1.541067099687334, - 1.5487587769900337, - 6.271271639873466, - ], - tspan=(0.0, 0.5)) - # 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)) < 10000 - end -end - -@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (global bounds)" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_sedov_blast_wave_sc_subcell.jl"), - positivity_variables_cons=["rho"], - positivity_variables_nonlinear=[pressure], - local_twosided_variables_cons=[], - local_onesided_variables_nonlinear=[], - l2=[ - 0.7869912572385168, - 0.39170886758882073, - 0.39613257454431977, - 1.2951760266455101, - ], - linf=[ - 5.156044534854053, - 3.6261667239538986, - 3.1807681416546085, - 6.3028422220287235, - ], - tspan=(0.0, 0.5)) - # 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)) < 10000 - end -end - @trixi_testset "elixir_euler_source_terms_waving_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_waving_flag.jl"), @@ -682,6 +626,62 @@ end end end +@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (local bounds)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_sedov_blast_wave_sc_subcell.jl"), + l2=[ + 0.6337774834710513, + 0.30377119245852724, + 0.3111372568571772, + 1.2976221893997268, + ], + linf=[ + 2.2064877103138207, + 1.541067099687334, + 1.5487587769900337, + 6.271271639873466, + ], + tspan=(0.0, 0.5)) + # 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)) < 10000 + end +end + +@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (global bounds)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_sedov_blast_wave_sc_subcell.jl"), + positivity_variables_cons=["rho"], + positivity_variables_nonlinear=[pressure], + local_twosided_variables_cons=[], + local_onesided_variables_nonlinear=[], + l2=[ + 0.7869912572385168, + 0.39170886758882073, + 0.39613257454431977, + 1.2951760266455101, + ], + linf=[ + 5.156044534854053, + 3.6261667239538986, + 3.1807681416546085, + 6.3028422220287235, + ], + tspan=(0.0, 0.5)) + # 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)) < 10000 + end +end + @trixi_testset "elixir_euler_rayleigh_taylor_instability.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_rayleigh_taylor_instability.jl"), From 40f475cf324694dae47122e8dd3d5ec983f6960f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 23 May 2024 15:18:04 +0200 Subject: [PATCH 17/19] Adapt dispatching --- src/callbacks_stage/subcell_bounds_check_2d.jl | 3 ++- src/callbacks_stage/subcell_limiter_idp_correction_2d.jl | 4 +++- src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl | 8 ++++---- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index c81ebc970a0..6f59e44e550 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -5,7 +5,8 @@ @muladd begin #! format: noindent -@inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache, +@inline function check_bounds(u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, + equations, solver, cache, limiter::SubcellLimiterIDP) (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index 04d7bf30821..bf972d8673b 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -5,7 +5,9 @@ @muladd begin #! format: noindent -function perform_idp_correction!(u, dt, mesh::AbstractMesh{2}, equations, dg, cache) +function perform_idp_correction!(u, dt, + mesh::::Union{TreeMesh{2}, StructuredMesh{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 @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index ab5182fb389..2cc5e2cae82 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::AbstractMesh{2}, +function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}}, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DG, uEltype) cache = create_cache(mesh, equations, @@ -57,7 +57,7 @@ function create_cache(mesh::AbstractMesh{2}, end function calc_volume_integral!(du, u, - mesh::Union{TreeMesh, StructuredMesh}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache) @@ -392,7 +392,7 @@ 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::AbstractMesh{2}, + u, mesh::Union{TreeMesh{2}, StructuredMesh{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 +430,7 @@ 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::AbstractMesh{2}, + u, mesh::Union{TreeMesh{2}, StructuredMesh{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 From 12400ee9b07581f2ceccd49bef04c4c1b429c253 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 23 May 2024 15:23:52 +0200 Subject: [PATCH 18/19] Fix typo --- src/callbacks_stage/subcell_limiter_idp_correction_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index bf972d8673b..0eb048ca5b8 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -6,7 +6,7 @@ #! format: noindent function perform_idp_correction!(u, dt, - mesh::::Union{TreeMesh{2}, StructuredMesh{2}}, + mesh::Union{TreeMesh{2}, StructuredMesh{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 From 83a8cef7b9c9d3e477484bd083ad9f586912abeb Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 23 May 2024 15:52:43 +0200 Subject: [PATCH 19/19] Remove not-needed parameters --- src/callbacks_stage/subcell_bounds_check.jl | 12 ++++++------ src/callbacks_stage/subcell_bounds_check_2d.jl | 5 ++--- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl index 3f3e151436f..6268fde6dd7 100644 --- a/src/callbacks_stage/subcell_bounds_check.jl +++ b/src/callbacks_stage/subcell_bounds_check.jl @@ -38,7 +38,7 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage) (; t, iter, alg) = integrator u = wrap_array(u_ode, mesh, equations, solver, cache) - @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache, + @trixi_timeit timer() "check_bounds" check_bounds(u, equations, solver, cache, solver.volume_integral) save_errors = callback.save_errors && (callback.interval > 0) && @@ -48,20 +48,20 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage) (iter + 1) >= integrator.opts.maxiters) # Maximum iterations reached if save_errors @trixi_timeit timer() "save_errors" save_bounds_check_errors(callback.output_directory, - u, t, iter + 1, + t, iter + 1, equations, solver.volume_integral) end end -@inline function check_bounds(u, mesh, equations, solver, cache, +@inline function check_bounds(u, equations, solver, cache, volume_integral::VolumeIntegralSubcellLimiting) - check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter) + check_bounds(u, equations, solver, cache, volume_integral.limiter) end -@inline function save_bounds_check_errors(output_directory, u, t, iter, equations, +@inline function save_bounds_check_errors(output_directory, t, iter, equations, volume_integral::VolumeIntegralSubcellLimiting) - save_bounds_check_errors(output_directory, u, t, iter, equations, + save_bounds_check_errors(output_directory, t, iter, equations, volume_integral.limiter) end diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 6f59e44e550..9c3664ab0c3 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -5,8 +5,7 @@ @muladd begin #! format: noindent -@inline function check_bounds(u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, - equations, solver, cache, +@inline function check_bounds(u, equations, solver, cache, limiter::SubcellLimiterIDP) (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter (; variable_bounds) = limiter.cache.subcell_limiter_coefficients @@ -104,7 +103,7 @@ return nothing end -@inline function save_bounds_check_errors(output_directory, u, time, iter, equations, +@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