Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Implement positivty limiter with numbers for conservative variables #113

Merged
merged 3 commits into from
Jun 19, 2023
Merged

Implement positivty limiter with numbers for conservative variables #113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
374fccf
Implement positivty limiter with numbers of cons vars
bennibolm Jun 14, 2023
5f21d5d
Merge branch 'IDPlimiting-positivity-cons' into IDPlimiting-positivit…
bennibolm Jun 16, 2023
43de3b7
Merge branch 'IDPlimiting-positivity-cons' into IDPlimiting-positivit…
bennibolm Jun 16, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
5 changes: 1 addition & 4 deletions examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
21 changes: 11 additions & 10 deletions src/solvers/dgsem_tree/indicators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,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)
Expand All @@ -240,25 +240,26 @@ 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
Comment on lines +243 to +245
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure if a tuple or an array is more useful. The length can vary depending on the example, which is why I now use an array. Is there still the possibility to use a tuple with undefined length?
The non-linear variables (in latter PRs) are still passed as functions in a tuple. But maybe the same type (tuple/array) should be used here (?).

cache::Cache
positivity_correction_factor::RealT
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)
Expand All @@ -270,7 +271,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, ")")
Expand All @@ -288,11 +289,11 @@ 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 *
"" => " "^11 *
"and positivity correction factor $(indicator.positivity_correction_factor)",
]
end
Expand Down
16 changes: 5 additions & 11 deletions src/solvers/dgsem_tree/indicators_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,7 +252,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
Expand All @@ -270,19 +270,13 @@ 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)
antidiffusive_flux1[variable, i, j, element]
val_flux1_local_ip1 = -inverse_weights[i] *
variable(get_node_vars(antidiffusive_flux1, equations,
dg, i + 1, j, element),
equations)
antidiffusive_flux1[variable, i + 1, j, element]
val_flux2_local = inverse_weights[j] *
variable(get_node_vars(antidiffusive_flux2, equations, dg,
i, j, element), equations)
antidiffusive_flux2[variable, i, j, element]
val_flux2_local_jp1 = -inverse_weights[j] *
variable(get_node_vars(antidiffusive_flux2, equations,
dg, i, j + 1, element),
equations)
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)
Expand Down