From 2da08630c97d3808fac9b52661fac32421098f26 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Fri, 24 May 2024 05:45:31 +0200 Subject: [PATCH] Add subcell limiting support for StructuredMesh (#1946) * Add structured mesh support * Fix non-periodic computation of bounds * Use local limiting and nonperiodic domain in source terms elixir * Use local limiting in free stream elixir * Remove not needed lines * Remove P4estMesh * Add non-periodic tests with local bounds * fmt * Fix test * Use `get_inverse_jacobian` instead of dispatching all routines * Simplify `perform_idp_correction!` * Revert stuff * Remove free stream elixir * Use sedov blast instead of source term setup; add news * Update dispatching for mesh types * Move new tests within test file * Adapt dispatching * Fix typo * Remove not-needed parameters --- NEWS.md | 1 + ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 114 +++++++++ src/callbacks_stage/subcell_bounds_check.jl | 12 +- .../subcell_bounds_check_2d.jl | 4 +- .../subcell_limiter_idp_correction_2d.jl | 27 ++- src/solvers/dgsem_structured/dg.jl | 3 + .../dg_2d_subcell_limiters.jl | 111 +++++++++ .../dgsem_structured/subcell_limiters_2d.jl | 220 ++++++++++++++++++ .../dgsem_tree/dg_2d_subcell_limiters.jl | 15 +- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 31 ++- test/test_structured_2d.jl | 56 +++++ 11 files changed, 557 insertions(+), 37 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_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/NEWS.md b/NEWS.md index ebc8d9cda39..ecbd70ce472 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,7 @@ for human readability. - Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh` ([#1792]). - New time integrator `PairedExplicitRK2`, implementing the second-order paired explicit Runge-Kutta method with [Convex.jl](https://github.com/jump-dev/Convex.jl) and [ECOS.jl](https://github.com/jump-dev/ECOS.jl) ([#1908]) +- 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_sedov_blast_wave_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl new file mode 100644 index 00000000000..5c11a7d15a7 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -0,0 +1,114 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +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, + x_pos = boundary_condition, + y_neg = boundary_condition, + y_pos = boundary_condition) + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + local_twosided_variables_cons = ["rho"], + local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, + min)], + max_iterations_newton = 40, # Default 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) + +# 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) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.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, + 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_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 c81ebc970a0..9c3664ab0c3 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -5,7 +5,7 @@ @muladd begin #! format: noindent -@inline function check_bounds(u, mesh::AbstractMesh{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 @@ -103,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 diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index 6f1723e2a98..0eb048ca5b8 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -5,29 +5,32 @@ @muladd begin #! format: noindent -function perform_idp_correction!(u, dt, mesh::TreeMesh2D, 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 @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) + # Sign switch as in apply_jacobian! + 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) + 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) + 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) + 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) + 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 * diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 424ed5e1e7e..5617ae90e3f 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -93,6 +93,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..5b65a475062 --- /dev/null +++ b/src/solvers/dgsem_structured/subcell_limiters_2d.jl @@ -0,0 +1,220 @@ +# 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 + + # 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) + u_outer = get_boundary_outer_state(boundary_conditions[1], + cache, t, 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) + u_outer = get_boundary_outer_state(boundary_conditions[2], + cache, t, 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) + u_outer = get_boundary_outer_state(boundary_conditions[3], + cache, t, 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) + u_outer = get_boundary_outer_state(boundary_conditions[4], + cache, t, 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 + + # 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) + u_outer = get_boundary_outer_state(boundary_conditions[1], + cache, t, 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) + u_outer = get_boundary_outer_state(boundary_conditions[2], + cache, t, 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) + u_outer = get_boundary_outer_state(boundary_conditions[3], + cache, t, 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) + u_outer = get_boundary_outer_state(boundary_conditions[4], + cache, t, 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..2cc5e2cae82 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) @@ -391,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::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 @@ -429,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::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 diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 33ae0599748..4095f0853f9 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -283,7 +283,7 @@ end end @inline function idp_local_twosided!(alpha, limiter, u, t, dt, semi, variable) - _, _, dg, cache = mesh_equations_solver_cache(semi) + 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 @@ -294,8 +294,9 @@ 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) var = u[variable, i, j, element] # Real Zalesak type limiter # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids" @@ -354,17 +355,18 @@ end return nothing end -@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, variable, - min_or_max) - _, equations, dg, cache = mesh_equations_solver_cache(semi) +@inline function idp_local_onesided!(alpha, limiter, u, t, dt, semi, + variable, min_or_max) + 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) # 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) + 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, @@ -407,7 +409,7 @@ 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) + 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 @@ -416,8 +418,9 @@ end 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) + 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.") @@ -467,16 +470,21 @@ end return nothing end +############################################################################### +# Global positivity limiting of nonlinear variables + @inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable) - _, equations, dg, cache = mesh_equations_solver_cache(semi) + 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")] @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) + # Compute bound u_local = get_node_vars(u, equations, dg, i, j, element) var = variable(u_local, equations) @@ -496,6 +504,9 @@ end 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..4b3aa5c87e4 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -626,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"),