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.

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

Add wet/dry functionality for the Multilayer SWE in 1D #38

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
126 changes: 126 additions & 0 deletions examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,126 @@

using OrdinaryDiffEq
using Trixi
using TrixiShallowWater

###############################################################################
# Semidiscretization of the shallow water equations

# By passing only a single value for rhos, the system recovers the standard shallow water equations
equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.812, rhos = 1.0)
patrickersing marked this conversation as resolved.
Show resolved Hide resolved

"""
initial_condition_beach(x, t, equations:: ShallowWaterMultiLayerEquationsWetDry1D)

Initial condition to simulate a wave running towards a beach and crashing. Difficult test
including both wetting and drying in the domain using slip wall boundary conditions.
The bottom topography is altered to be differentiable on the domain [0,8] and
differs from the reference below.

The water height and speed functions used here, are adapted from the initial condition
found in section 5.2 of the paper:
- Andreas Bollermann, Sebastian Noelle, Maria Lukáčová-Medvidová (2011)
Finite volume evolution Galerkin methods for the shallow water equations with dry beds\n
[DOI: 10.4208/cicp.220210.020710a](https://dx.doi.org/10.4208/cicp.220210.020710a)
"""
function initial_condition_beach(x, t, equations::ShallowWaterMultiLayerEquations1D)
D = 1
delta = 0.02
gamma = sqrt((3 * delta) / (4 * D))
x_a = sqrt((4 * D) / (3 * delta)) * acosh(sqrt(20))

f = D + 40 * delta * sech(gamma * (8 * x[1] - x_a))^2

# steep curved beach
b = 0.01 + 99 / 409600 * 4^x[1]

if x[1] >= 6
H = b
v = 0.0
else
H = f
v = sqrt(equations.gravity / D) * H
end

# It is mandatory to shift the water level at dry areas to make sure the water height h
# stays positive. The system would not be stable for h set to a hard 0 due to division by h in
# the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
# with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above
# for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0.
# This default value can be changed within the constructor call depending on the simulation setup.
H = max(H, b + equations.threshold_limiter)
return prim2cons(SVector(H, v, b), equations)
end

initial_condition = initial_condition_beach
boundary_condition = boundary_condition_slip_wall

###############################################################################
# Get the DG approximation space

volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
DissipationLocalLaxFriedrichs()),
hydrostatic_reconstruction_ersing_etal),
FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal,
hydrostatic_reconstruction_ersing_etal))

basis = LobattoLegendreBasis(3)

indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = waterheight_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)

solver = DGSEM(basis, surface_flux, volume_integral)

###############################################################################
# Create the TreeMesh for the domain [0, 8]

coordinates_min = 0.0
coordinates_max = 8.0

mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 7,
n_cells_max = 10_000,
periodicity = false)

# create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_condition)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
save_analysis = false,
extra_analysis_integrals = (energy_kinetic,
energy_internal))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(dt = 0.5,
save_initial_solution = true,
save_final_solution = true)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)

stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))

###############################################################################
# run the simulation

sol = solve(ode, SSPRK43(stage_limiter!);
ode_default_options()..., callback = callbacks);

summary_callback() # print the timer summary
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@

using OrdinaryDiffEq
using Trixi
using TrixiShallowWater

###############################################################################
# Semidiscretization of the multilayer shallow water equations for a dam break
# test over a dry domain with a discontinuous bottom topography function

equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
rhos = (0.8, 0.85, 0.9, 0.95, 1.0))

# Initial condition of a dam break with a discontinuous water heights and bottom topography.
# To test the wet/dry functionality this test case considers a dam break over a dry domain with
# multiple layers.
# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh
# refinement level is changed the initial condition below may need changed as well to
# ensure that the discontinuities lie on an element interface.
function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations1D)
# Set the discontinuity
if x[1] <= 10.0
H = [3.5, 3.0, 2.5, 2.0, 1.5]
b = 0.0
else
# Right side of the domain is dry
b = 0.5
H = [b, b, b, b, b]
end

# It is mandatory to shift the water level at dry areas to make sure the water height h
# stays positive. The system would not be stable for h set to a hard 0 due to division by h in
# the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
# with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above
# for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0.
# This default value can be changed within the constructor call depending on the simulation setup.
for i in reverse(eachlayer(equations))
if i == nlayers(equations)
H[i] = max(H[i], b + equations.threshold_limiter)
else
H[i] = max(H[i], H[i + 1] + equations.threshold_limiter)
end
end

# Initialize zero velocity
v = zero(H)

return prim2cons(SVector(H..., v..., b), equations)
end

initial_condition = initial_condition_dam_break

###############################################################################
# Get the DG approximation space

volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
DissipationLocalLaxFriedrichs()),
hydrostatic_reconstruction_ersing_etal),
FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal,
hydrostatic_reconstruction_ersing_etal))

basis = LobattoLegendreBasis(3)

indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = waterheight_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)

solver = DGSEM(basis, surface_flux, volume_integral)

###############################################################################
# Get the TreeMesh and setup a non-periodic mesh

coordinates_min = 0.0
coordinates_max = 20.0
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 5,
n_cells_max = 10000,
periodicity = false)

boundary_condition = boundary_condition_slip_wall

# create the semidiscretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_condition)

###############################################################################
# ODE solvers

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# Callbacks

summary_callback = SummaryCallback()

analysis_interval = 500
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
save_analysis = false)

stepsize_callback = StepsizeCallback(cfl = 0.2)
patrickersing marked this conversation as resolved.
Show resolved Hide resolved

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 500,
save_initial_solution = true,
save_final_solution = true)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
stepsize_callback, save_solution)

stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))

###############################################################################
# run the simulation

sol = solve(ode, SSPRK43(stage_limiter!), dt = 1.0,
save_everystep = false, callback = callbacks, adaptive = false);
summary_callback() # print the timer summary
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@

using OrdinaryDiffEq
using Trixi
using TrixiShallowWater

###############################################################################
# Semidiscretization of the shallow water equations

# By passing only a single value for rhos, the system recovers the standard shallow water equations
equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, rhos = 1.0)
patrickersing marked this conversation as resolved.
Show resolved Hide resolved

"""
initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterMultiLayerEquations1D)

Well-known initial condition to test the [`hydrostatic_reconstruction_ersing_etal`](@ref) and its
wet-dry mechanics. This test has analytical solutions. The initial condition is defined by the
analytical solution at time t=0. The bottom topography defines a bowl and the water level is given
by an oscillating lake.

The original test and its analytical solution in two dimensions were first presented in
- William C. Thacker (1981)
Some exact solutions to the nonlinear shallow-water wave equations
[DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882).

The particular setup below is taken from Section 6.2 of
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018)
An entropy stable discontinuous Galerkin method for the shallow water equations on
curvilinear meshes with wet/dry fronts accelerated by GPUs
[DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038).
"""
function initial_condition_parabolic_bowl(x, t,
equations::ShallowWaterMultiLayerEquations1D)
a = 1
h_0 = 0.1
sigma = 0.5
ω = sqrt(2 * equations.gravity * h_0) / a

v = -sigma * ω * sin(ω * t)

b = h_0 * x[1]^2 / a^2

H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) - sigma) + h_0

#=
It is mandatory to shift the water level at dry areas to make sure the water height h
stays positive. The system would not be stable for h set to a hard 0 due to division by h in
the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above
for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0.
This default value can be changed within the constructor call depending on the simulation setup.
=#
H = max(H, b + equations.threshold_limiter)
return prim2cons(SVector(H, v, b), equations)
end

initial_condition = initial_condition_parabolic_bowl

###############################################################################
# Get the DG approximation space

volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
DissipationLocalLaxFriedrichs()),
hydrostatic_reconstruction_ersing_etal),
FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal,
hydrostatic_reconstruction_ersing_etal))

basis = LobattoLegendreBasis(5)

indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = waterheight_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)

solver = DGSEM(basis, surface_flux, volume_integral)

###############################################################################
# Create the TreeMesh for the domain [-2, 2]

coordinates_min = -2.0
coordinates_max = 2.0

mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 6,
n_cells_max = 10_000)

# create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
save_analysis = false,
extra_analysis_integrals = (energy_kinetic,
energy_internal))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 1000,
save_initial_solution = true,
save_final_solution = true)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)

stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))

###############################################################################
# run the simulation

sol = solve(ode, SSPRK43(stage_limiter!);
patrickersing marked this conversation as resolved.
Show resolved Hide resolved
ode_default_options()..., callback = callbacks);

summary_callback() # print the timer summary
Loading
Loading