Skip to content

Commit 5fe8999

Browse files
patrickersingandrewwinters5000amrueda
authored
Subcell-limiting support for SWE (#111)
* add nonconservative sc-fluxes for multilayer-swe * separate velocity desingularization callback, specialize subcell limiter * add specialized two-sided limiter that uses the total water height * add elixirs for testing * fix documentation link * increase allowed allocations for tests with custom time integrators * increase atol for two-sided limiter * apply formatter * remove custom time integrator workaround for the positivity limiter * fix issues with p4est and subcell limiting * add boundary condition to blast wave elixir to increase coverage * use the same name for the subcell fluxes implementation * create news entry * Revert "use the same name for the subcell fluxes implementation" This reverts commit c64084e. * Apply suggestions from code review Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se> * rename elixirs * fix reference for `dg_2d_subcell_limiters.jl` * add specialized two-sided limiter for p4est mesh * update test values after bug fix * add p4est blast wave test and fix bug in p4est implementation * add p4est test with wall BC * set periodicity=true in p4est test * Update examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_sc_subcell.jl Co-authored-by: Andrés Rueda-Ramírez <aruedara@uni-koeln.de> * Update src/equations/shallow_water_multilayer_2d.jl Co-authored-by: Andrés Rueda-Ramírez <aruedara@uni-koeln.de> * refactor computation of limiting coefficients * update argument types for calc_bounds_twosided! * add subcell-limiting upstream test * fix elixir path for upstream test * Typo fix in wet dry elixir --------- Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se> Co-authored-by: Andrés Rueda-Ramírez <aruedara@uni-koeln.de>
1 parent f95d39b commit 5fe8999

17 files changed

+1674
-104
lines changed

‎NEWS.md‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,19 @@ TrixiShallowWater.jl follows the interpretation of
55
used in the Julia ecosystem. Notable changes will be documented in this file
66
for human readability.
77

8+
## Changes in the v0.2 lifecycle
9+
10+
#### Added
11+
- Introduce node-wise limiting functionality for the `ShallowWaterMultiLayerEquations2D`. ([#111])
12+
13+
#### Changed
14+
- Velocity desingularization procedure has been moved into a distinct `VelocityDesingularization`
15+
stage callback. ([#111])
16+
17+
#### Deprecated
18+
19+
#### Removed
20+
821
## Changes when updating to v0.2 from v0.1.x
922

1023
#### Added

‎examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_blast_wet_dry_sc_subcell.jl‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
2+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
3+
using Trixi
4+
using TrixiShallowWater
5+
6+
###############################################################################
7+
# Semidiscretization of the multilayer shallow water equations with a single layer and a bottom
8+
# topography function for a blast wave test with a wet/dry front with discontinuous initial conditions.
9+
# This example combines the subcell limiter with a velocity desingularization callback to ensure
10+
# robustness at wet/dry fronts.
11+
12+
equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.45,
13+
rhos = (1.0),
14+
threshold_desingularization = 1e-6)
15+
16+
function initial_condition_blast_wave(x, t, equations::ShallowWaterMultiLayerEquations2D)
17+
# Set up polar coordinates
18+
RealT = eltype(x)
19+
inicenter = SVector(convert(RealT, 2.0), convert(RealT, 2.0))
20+
x_norm = x[1] - inicenter[1]
21+
y_norm = x[2] - inicenter[2]
22+
r = sqrt(x_norm^2 + y_norm^2)
23+
phi = atan(y_norm, x_norm)
24+
sin_phi, cos_phi = sincos(phi)
25+
26+
# Calculate primitive variables
27+
H = r > 0.5f0 ? 2.0f0 : 4.0f0
28+
v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
29+
v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi
30+
31+
b = min(sqrt((2 - x[1])^2 + (2 - x[2])^2)^2, 4.0)
32+
33+
H = max(H, b + equations.threshold_limiter)
34+
return prim2cons(SVector(H, v1, v2, b), equations)
35+
end
36+
37+
initial_condition = initial_condition_blast_wave
38+
39+
boundary_condition = Dict(:x_neg => boundary_condition_slip_wall,
40+
:y_neg => boundary_condition_slip_wall,
41+
:x_pos => boundary_condition_slip_wall,
42+
:y_pos => boundary_condition_slip_wall)
43+
44+
###############################################################################
45+
# Get the DG approximation space
46+
47+
polydeg = 4
48+
volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
49+
surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
50+
DissipationLocalLaxFriedrichs()),
51+
hydrostatic_reconstruction_ersing_etal),
52+
FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal_local_jump,
53+
hydrostatic_reconstruction_ersing_etal))
54+
55+
basis = LobattoLegendreBasis(polydeg)
56+
limiter_idp = SubcellLimiterIDP(equations, basis;
57+
positivity_variables_cons = ["h1"],
58+
local_twosided_variables_cons = ["h1"],
59+
positivity_correction_factor = 0.5,)
60+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
61+
volume_flux_dg = volume_flux,
62+
volume_flux_fv = surface_flux)
63+
solver = DGSEM(basis, surface_flux, volume_integral)
64+
65+
###############################################################################
66+
# Get the P4estMesh
67+
68+
coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
69+
coordinates_max = (4.0, 4.0) # maximum coordinates (max(x), max(y))
70+
71+
trees_per_dimension = (1, 1)
72+
73+
mesh = P4estMesh(trees_per_dimension, polydeg = 3,
74+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
75+
initial_refinement_level = 5,
76+
periodicity = false)
77+
78+
# Create the semi discretization object
79+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
80+
boundary_conditions = boundary_condition)
81+
82+
###############################################################################
83+
# ODE solver
84+
85+
tspan = (0.0, 1.0)
86+
ode = semidiscretize(semi, tspan)
87+
88+
summary_callback = SummaryCallback()
89+
90+
analysis_interval = 100
91+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
92+
analysis_polydeg = polydeg,
93+
extra_analysis_errors = (:conservation_error,))
94+
95+
stepsize_callback = StepsizeCallback(cfl = 0.25)
96+
97+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
98+
99+
save_solution = SaveSolutionCallback(dt = 0.1,
100+
save_initial_solution = true,
101+
save_final_solution = true,
102+
extra_node_variables = (:limiting_coefficient,))
103+
104+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
105+
stepsize_callback)
106+
107+
###############################################################################
108+
109+
# Setup the stage callbacks
110+
stage_limiter! = VelocityDesingularization()
111+
stage_callbacks = (SubcellLimiterIDPCorrection(), stage_limiter!)
112+
113+
# run the simulation
114+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
115+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
116+
ode_default_options()...,
117+
callback = callbacks, adaptive = false);

‎examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_convergence_sc_subcell.jl‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
using TrixiShallowWater
4+
5+
###############################################################################
6+
# Semidiscretization of the multilayer shallow water equations with three layers
7+
8+
equations = ShallowWaterMultiLayerEquations2D(gravity = 1.1,
9+
rhos = (0.9, 1.0, 1.1))
10+
11+
initial_condition = initial_condition_convergence_test
12+
13+
###############################################################################
14+
# Get the DG approximation space
15+
16+
polydeg = 3
17+
volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
18+
surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
19+
20+
basis = LobattoLegendreBasis(polydeg)
21+
limiter_idp = SubcellLimiterIDP(equations, basis;)
22+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
23+
volume_flux_dg = volume_flux,
24+
volume_flux_fv = surface_flux)
25+
solver = DGSEM(basis, surface_flux, volume_integral)
26+
27+
###############################################################################
28+
# Get the P4estMesh and setup a periodic mesh
29+
30+
coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
31+
coordinates_max = (sqrt(2.0), sqrt(2.0)) # maximum coordinates (max(x), max(y))
32+
33+
# Create P4estMesh with 4 x 4 elements
34+
trees_per_dimension = (2, 2)
35+
mesh = P4estMesh(trees_per_dimension, polydeg = 3,
36+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
37+
initial_refinement_level = 1,
38+
periodicity = true)
39+
40+
# create the semi discretization object
41+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
42+
source_terms = source_terms_convergence_test)
43+
44+
###############################################################################
45+
# ODE solvers, callbacks etc.
46+
47+
tspan = (0.0, 0.5)
48+
ode = semidiscretize(semi, tspan)
49+
50+
summary_callback = SummaryCallback()
51+
52+
analysis_interval = 500
53+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
54+
55+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
56+
57+
save_solution = SaveSolutionCallback(interval = 500,
58+
save_initial_solution = true,
59+
save_final_solution = true)
60+
61+
stepsize_callback = StepsizeCallback(cfl = 0.25)
62+
63+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
64+
stepsize_callback)
65+
66+
###############################################################################
67+
# run the simulation
68+
69+
stage_callbacks = (SubcellLimiterIDPCorrection(),)
70+
71+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
72+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
73+
ode_default_options()...,
74+
callback = callbacks, adaptive = false);

‎examples/tree_2d_dgsem/elixir_shallowwater_multilayer_blast_wet_dry_sc_subcell.jl‎

Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
2+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
3+
using Trixi
4+
using TrixiShallowWater
5+
6+
###############################################################################
7+
# Semidiscretization of the multilayer shallow water equations with a single layer and a bottom
8+
# topography function for a blast wave test with a wet/dry front with discontinuous initial conditions.
9+
# This example combines the subcell limiter with a velocity desingularization callback to ensure
10+
# robustness at wet/dry fronts.
11+
12+
equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.45,
13+
rhos = (1.0),
14+
threshold_desingularization = 1e-6)
15+
16+
function initial_condition_blast_wave(x, t, equations::ShallowWaterMultiLayerEquations2D)
17+
# Set up polar coordinates
18+
RealT = eltype(x)
19+
inicenter = SVector(convert(RealT, 2.0), convert(RealT, 2.0))
20+
x_norm = x[1] - inicenter[1]
21+
y_norm = x[2] - inicenter[2]
22+
r = sqrt(x_norm^2 + y_norm^2)
23+
phi = atan(y_norm, x_norm)
24+
sin_phi, cos_phi = sincos(phi)
25+
26+
# Calculate primitive variables
27+
H = r > 0.5f0 ? 2.0f0 : 4.0f0
28+
v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
29+
v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi
30+
31+
b = min(sqrt((2 - x[1])^2 + (2 - x[2])^2)^2, 4.0)
32+
33+
H = max(H, b + equations.threshold_limiter)
34+
return prim2cons(SVector(H, v1, v2, b), equations)
35+
end
36+
37+
initial_condition = initial_condition_blast_wave
38+
39+
###############################################################################
40+
# Get the DG approximation space
41+
42+
polydeg = 4
43+
volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
44+
surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
45+
DissipationLocalLaxFriedrichs()),
46+
hydrostatic_reconstruction_ersing_etal),
47+
FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal_local_jump,
48+
hydrostatic_reconstruction_ersing_etal))
49+
50+
basis = LobattoLegendreBasis(polydeg)
51+
limiter_idp = SubcellLimiterIDP(equations, basis;
52+
positivity_variables_cons = ["h1"],
53+
local_twosided_variables_cons = ["h1"],
54+
positivity_correction_factor = 0.5,)
55+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
56+
volume_flux_dg = volume_flux,
57+
volume_flux_fv = surface_flux)
58+
solver = DGSEM(basis, surface_flux, volume_integral)
59+
60+
###############################################################################
61+
# Get the TreeMesh
62+
63+
coordinates_min = (0.0, 0.0)
64+
coordinates_max = (4.0, 4.0)
65+
mesh = TreeMesh(coordinates_min, coordinates_max,
66+
initial_refinement_level = 5,
67+
n_cells_max = 10_000,
68+
periodicity = false)
69+
70+
# Create the semi discretization object
71+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
72+
boundary_conditions = BoundaryConditionDirichlet(initial_condition))
73+
74+
###############################################################################
75+
# ODE solver
76+
77+
tspan = (0.0, 1.0)
78+
ode = semidiscretize(semi, tspan)
79+
80+
summary_callback = SummaryCallback()
81+
82+
analysis_interval = 100
83+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
84+
analysis_polydeg = polydeg,
85+
extra_analysis_errors = (:conservation_error,))
86+
87+
stepsize_callback = StepsizeCallback(cfl = 0.25)
88+
89+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
90+
91+
save_solution = SaveSolutionCallback(dt = 0.1,
92+
save_initial_solution = true,
93+
save_final_solution = true,
94+
extra_node_variables = (:limiting_coefficient,))
95+
96+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
97+
stepsize_callback)
98+
99+
###############################################################################
100+
101+
# Setup the stage callbacks
102+
stage_limiter! = VelocityDesingularization()
103+
stage_callbacks = (SubcellLimiterIDPCorrection(), stage_limiter!)
104+
105+
# run the simulation
106+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
107+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
108+
ode_default_options()...,
109+
callback = callbacks, adaptive = false);

‎examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence_sc_subcell.jl‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
2+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
3+
using Trixi
4+
using TrixiShallowWater
5+
6+
###############################################################################
7+
# Semidiscretization of the multilayer shallow water equations with three layers
8+
9+
equations = ShallowWaterMultiLayerEquations2D(gravity = 1.1,
10+
rhos = (0.9, 1.0, 1.1))
11+
12+
initial_condition = initial_condition_convergence_test
13+
14+
###############################################################################
15+
# Get the DG approximation space
16+
17+
polydeg = 4
18+
volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
19+
surface_flux = (FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()),
20+
flux_nonconservative_ersing_etal_local_jump)
21+
22+
basis = LobattoLegendreBasis(polydeg)
23+
limiter_idp = SubcellLimiterIDP(equations, basis;)
24+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
25+
volume_flux_dg = volume_flux,
26+
volume_flux_fv = surface_flux)
27+
solver = DGSEM(basis, surface_flux, volume_integral)
28+
29+
###############################################################################
30+
# Get the TreeMesh and setup a periodic mesh
31+
32+
coordinates_min = (0.0, 0.0)
33+
coordinates_max = (sqrt(2.0), sqrt(2.0))
34+
mesh = TreeMesh(coordinates_min, coordinates_max,
35+
initial_refinement_level = 2,
36+
n_cells_max = 100_000,
37+
periodicity = true)
38+
39+
# create the semi discretization object
40+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
41+
source_terms = source_terms_convergence_test)
42+
43+
###############################################################################
44+
# ODE solvers, callbacks etc.
45+
46+
tspan = (0.0, 1.0)
47+
ode = semidiscretize(semi, tspan)
48+
49+
summary_callback = SummaryCallback()
50+
51+
analysis_interval = 500
52+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
53+
54+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
55+
56+
save_solution = SaveSolutionCallback(interval = 500,
57+
save_initial_solution = true,
58+
save_final_solution = true)
59+
60+
stepsize_callback = StepsizeCallback(cfl = 0.3)
61+
62+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
63+
stepsize_callback)
64+
65+
###############################################################################
66+
# run the simulation
67+
68+
stage_callbacks = (SubcellLimiterIDPCorrection(),)
69+
70+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
71+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
72+
ode_default_options()...,
73+
callback = callbacks, adaptive = false);

0 commit comments

Comments
 (0)