From 374fccfc7347a7c0ee17692d88d1f3073d5dccaa Mon Sep 17 00:00:00 2001 From: Benjamin Bolm Date: Wed, 14 Jun 2023 11:07:46 +0200 Subject: [PATCH] Implement positivty limiter with numbers of cons vars --- .github/workflows/SpellCheck.yml | 2 +- Project.toml | 4 ++-- .../elixir_euler_shockcapturing_subcell.jl | 2 +- ...lixir_eulermulti_shock_bubble_sc_subcell.jl | 5 +---- src/auxiliary/precompile.jl | 2 +- src/callbacks_step/analysis_dgmulti.jl | 8 ++++---- src/solvers/dgsem_tree/indicators.jl | 18 +++++++++--------- src/solvers/dgsem_tree/indicators_1d.jl | 2 +- src/solvers/dgsem_tree/indicators_2d.jl | 10 +++++----- src/time_integration/methods_SSP.jl | 1 - 10 files changed, 25 insertions(+), 29 deletions(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 9e27ce20d0..c4ab3a9855 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v3 - name: Check spelling - uses: crate-ci/typos@v1.14.11 + uses: crate-ci/typos@v1.14.12 diff --git a/Project.toml b/Project.toml index e8922d2639..9d51e4dcff 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.28-pre" +version = "0.5.29-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" @@ -73,7 +73,7 @@ Requires = "1.1" SciMLBase = "1.90" Setfield = "0.8, 1" SimpleUnPack = "1.1" -StartUpDG = "0.16" +StartUpDG = "0.17" Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8" StaticArrayInterface = "1.4" StaticArrays = "1" diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 8061a4006f..528b6486ea 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -39,7 +39,7 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) indicator_sc = IndicatorIDP(equations, basis; - positivity=true, variables_cons=(Trixi.density,), positivity_correction_factor=0.5) + positivity=true, variables_cons=[1], positivity_correction_factor=0.5) volume_integral = VolumeIntegralSubcellLimiting(indicator_sc; volume_flux_dg=volume_flux, volume_flux_fv=surface_flux) diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_sc_subcell.jl index 25b4b79241..02e4d8c15a 100644 --- a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_sc_subcell.jl @@ -83,12 +83,9 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) -density1(u, equations::CompressibleEulerMulticomponentEquations2D) = u[1+3] -density2(u, equations::CompressibleEulerMulticomponentEquations2D) = u[2+3] - indicator_sc = IndicatorIDP(equations, basis; positivity=true, - variables_cons=(density1, density2)) + variables_cons=[(i+3 for i in eachcomponent(equations))...]) volume_integral=VolumeIntegralSubcellLimiting(indicator_sc; volume_flux_dg=volume_flux, volume_flux_fv=surface_flux) diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 0695e72efa..3b1cb58e14 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -146,7 +146,7 @@ function _precompile_manual_() function equations_types_1d(RealT) ( LinearScalarAdvectionEquation1D{RealT}, - HyperbolicDiffusionEquation1D{RealT}, + HyperbolicDiffusionEquations1D{RealT}, CompressibleEulerEquations1D{RealT}, IdealGlmMhdEquations1D{RealT}, ) diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index 2fbd8eda87..18640c9379 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/src/callbacks_step/analysis_dgmulti.jl @@ -155,10 +155,10 @@ function integrate(func::typeof(enstrophy), u, gradient_z_quadrature_values = local_gradient_quadrature_values[3][Threads.threadid()] # interpolate to quadrature on each element - apply_to_each_field(mul_by(dg.basis.Vq), u_quadrature_values, view(u, :, e)) - apply_to_each_field(mul_by(dg.basis.Vq), gradient_x_quadrature_values, view(gradients_x, :, e)) - apply_to_each_field(mul_by(dg.basis.Vq), gradient_y_quadrature_values, view(gradients_y, :, e)) - apply_to_each_field(mul_by(dg.basis.Vq), gradient_z_quadrature_values, view(gradients_z, :, e)) + apply_to_each_field(mul_by!(dg.basis.Vq), u_quadrature_values, view(u, :, e)) + apply_to_each_field(mul_by!(dg.basis.Vq), gradient_x_quadrature_values, view(gradients_x, :, e)) + apply_to_each_field(mul_by!(dg.basis.Vq), gradient_y_quadrature_values, view(gradients_y, :, e)) + apply_to_each_field(mul_by!(dg.basis.Vq), gradient_z_quadrature_values, view(gradients_z, :, e)) # integrate over the element for i in eachindex(u_quadrature_values) diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index e0c4706f7b..3424658222 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -215,7 +215,7 @@ end """ IndicatorIDP(equations::AbstractEquations, basis; positivity=false, - variables_cons=(), + variables_cons=[], positivity_correction_factor=0.1) Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref) @@ -237,9 +237,9 @@ The bounds are calculated using the low-order FV solution. The positivity limite !!! warning "Experimental implementation" This is an experimental feature and may change in future releases. """ -struct IndicatorIDP{RealT<:Real, LimitingVariablesCons, Cache} <: AbstractIndicator +struct IndicatorIDP{RealT<:Real, Cache} <: AbstractIndicator positivity::Bool - variables_cons::LimitingVariablesCons # Positivity of conservative variables + variables_cons::Vector{Int} # Impose positivity for conservative variables cache::Cache positivity_correction_factor::RealT end @@ -247,14 +247,14 @@ end # this method is used when the indicator is constructed as for shock-capturing volume integrals function IndicatorIDP(equations::AbstractEquations, basis; positivity=false, - variables_cons=(), + variables_cons=[], positivity_correction_factor=0.1) number_bounds = positivity * length(variables_cons) cache = create_cache(IndicatorIDP, equations, basis, number_bounds) - IndicatorIDP{typeof(positivity_correction_factor), typeof(variables_cons), typeof(cache)}(positivity, - variables_cons, cache, positivity_correction_factor) + IndicatorIDP{typeof(positivity_correction_factor), typeof(cache)}(positivity, variables_cons, + cache, positivity_correction_factor) end function Base.show(io::IO, indicator::IndicatorIDP) @@ -266,7 +266,7 @@ function Base.show(io::IO, indicator::IndicatorIDP) print(io, "No limiter selected => pure DG method") else print(io, "limiter=(") - positivity && print(io, "Positivity with variables $(indicator.variables_cons)") + positivity && print(io, "positivity") print(io, "), ") end print(io, ")") @@ -284,9 +284,9 @@ function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorIDP) else setup = ["limiter" => ""] if positivity - string = "Positivity with variables $(indicator.variables_cons))" + string = "positivity with conservative variables $(indicator.variables_cons)" setup = [setup..., "" => string] - setup = [setup..., "" => " "^14 * "and positivity correction factor $(indicator.positivity_correction_factor)"] + setup = [setup..., "" => " "^11 * "and positivity correction factor $(indicator.positivity_correction_factor)"] end end summary_box(io, "IndicatorIDP", setup) diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index c1a8816124..7086d77a1a 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -332,7 +332,7 @@ function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkRayHesthaven})( for direction in eachdirection(mesh.tree) if !has_any_neighbor(mesh.tree, cell_id, direction) - neighbor_ids[direction] = element_id + neighbor_ids[direction] = element continue end if has_neighbor(mesh.tree, cell_id, direction) diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 33fed020e6..9fd8c9fa2a 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -242,7 +242,7 @@ 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 = variable(get_node_vars(u, equations, dg, i, j, element), equations) + var = u[variable, i, j, element] if var < 0 error("Safe $variable is not safe. element=$element, node: $i $j, value=$var") end @@ -259,10 +259,10 @@ end # Calculate Pm # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. - val_flux1_local = inverse_weights[i] * variable(get_node_vars(antidiffusive_flux1, equations, dg, i, j, element), equations) - val_flux1_local_ip1 = -inverse_weights[i] * variable(get_node_vars(antidiffusive_flux1, equations, dg, i+1, j, element), equations) - val_flux2_local = inverse_weights[j] * variable(get_node_vars(antidiffusive_flux2, equations, dg, i, j, element), equations) - val_flux2_local_jp1 = -inverse_weights[j] * variable(get_node_vars(antidiffusive_flux2, equations, dg, i, j+1, element), equations) + val_flux1_local = inverse_weights[i] * antidiffusive_flux1[variable, i, j, element] + val_flux1_local_ip1 = -inverse_weights[i] * antidiffusive_flux1[variable, i+1, j, element] + val_flux2_local = inverse_weights[j] * antidiffusive_flux2[variable, i, j, element] + val_flux2_local_jp1 = -inverse_weights[j] * antidiffusive_flux2[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) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index f0786634d5..82f2b251b3 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -132,7 +132,6 @@ function solve(ode::ODEProblem, alg=SimpleSSPRK33()::SimpleAlgorithmSSP; end function solve!(integrator::SimpleIntegratorSSP) - @unpack indicator = integrator.p.solver.volume_integral @unpack prob = integrator.sol @unpack alg = integrator t_end = last(prob.tspan)