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

Shallow water exner equations 1D #51

Merged
merged 13 commits into from
Aug 15, 2024
6 changes: 4 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,18 @@ version = "0.1.0-pre"
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
andrewwinters5000 marked this conversation as resolved.
Show resolved Hide resolved
Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[compat]
LinearAlgebra = "1"
MuladdMacro = "0.2.2"
Printf = "1"
Roots = "2.1.6"
Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8"
StaticArrays = "1"
Trixi = "0.7, 0.8"
julia = "1.8"
Printf = "1"
116 changes: 116 additions & 0 deletions examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@

using OrdinaryDiffEq
using Trixi
using TrixiShallowWater
using Roots

###############################################################################
# Semidiscretization of the shallow water exner equations for a channel flow problem
# with sediment transport

equations = ShallowWaterExnerEquations1D(gravity_constant = 9.81, rho_f = 1.0, rho_s = 1.0,
porosity = 0.4,
sediment_model = GrassModel(A_g = 0.01))

# Initial condition for for a channel flow problem over a sediment hump.
# An asymptotic solution based on the method of characteristics was derived under a rigid-lid
# approximation in chapter 3.5.1 of the thesis:
# - Justin Hudson (2001)
# "Numerical Techniques for Morphodynamic Modelling"
# [PhD Thesis, University of Reading](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=f78dbae9cfbb12ae823975c6ce9d2585b40417ba)
function initial_condition_channel(x, t,

Check warning on line 21 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L21

Added line #L21 was not covered by tests
equations::ShallowWaterExnerEquations1D)
(; sediment_model, porosity_inv) = equations

Check warning on line 23 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L23

Added line #L23 was not covered by tests

H_ref = 10.0 # Reference water level
hv = 10.0

Check warning on line 26 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L25-L26

Added lines #L25 - L26 were not covered by tests

# Use the method of characteristics to find the asymptotic solution for the bed height, see
# Eq. 3.16 in the reference.
# First use an iterative method to predict x0
f(x0) = x[1] - x0 -

Check warning on line 31 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L31

Added line #L31 was not covered by tests
sediment_model.A_g * porosity_inv * 3 * hv^3 * t *
(H_ref - sinpi((x0 - 300) / 200)^2)^(-4)
fx = Roots.ZeroProblem(f, 400.0)
x0 = Roots.solve(fx, atol = 0.0, rtol = 0.0)

Check warning on line 35 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L34-L35

Added lines #L34 - L35 were not covered by tests

# If the result is outside 300 < x < 500 the result is invalid and instead compute x0 from
if x0 > 500 || x0 < 300
x0 = x[1] -

Check warning on line 39 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L38-L39

Added lines #L38 - L39 were not covered by tests
sediment_model.A_g * porosity_inv * 3 * hv^3 * t * H_ref^(-4)
end

# Compute the sediment and water height
300 < x0 < 500 ? h_b = 0.1 + sinpi((x0 - 300) / 200)^2 : h_b = 0.1
h = H_ref - h_b

Check warning on line 45 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L44-L45

Added lines #L44 - L45 were not covered by tests

return SVector(h, hv, h_b)

Check warning on line 47 in examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl#L47

Added line #L47 was not covered by tests
end

initial_condition = initial_condition_channel

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

volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
surface_flux = (FluxPlusDissipation(flux_ersing_etal, dissipation_roe),
flux_nonconservative_ersing_etal)

basis = LobattoLegendreBasis(6)

indicator_sc = IndicatorHennemannGassner(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = water_sediment_height)
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 periodic mesh

coordinates_min = 0.0
coordinates_max = 1000.0
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 4,
n_cells_max = 10_000,
periodicity = true)

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

###############################################################################
# ODE solver

tspan = (0.0, 30_000.0)
ode = semidiscretize(semi, tspan)

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

summary_callback = SummaryCallback()

analysis_interval = 10000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

stepsize_callback = StepsizeCallback(cfl = 0.8)

save_solution = SaveSolutionCallback(dt = 10_000.0,
save_initial_solution = true,
save_final_solution = true)

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

###############################################################################
# run the simulation
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
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
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
using OrdinaryDiffEq
using Trixi
using TrixiShallowWater

###############################################################################
# Semidiscretization of the SWE-Exner equations with source terms for convergence testing

# Equations with Grass model
equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5,
rho_s = 1.0, porosity = 0.5,
friction = ManningFriction(n = 0.0),
sediment_model = GrassModel(A_g = 0.01))

initial_condition = initial_condition_convergence_test

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

volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)

solver = DGSEM(polydeg = 4,
surface_flux = surface_flux,
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

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

coordinates_min = 0.0
coordinates_max = sqrt(2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 2,
n_cells_max = 10_000,
periodicity = true)

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

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

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

summary_callback = SummaryCallback()

analysis_interval = 200
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

stepsize_callback = StepsizeCallback(cfl = 1.0)

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

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

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

sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
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
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
using OrdinaryDiffEq
using Trixi
using TrixiShallowWater

###############################################################################
# Semidiscretization of the SWE-Exner equations with source terms for convergence testing

# Equations with Meyer-Peter-Mueller model
equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5,
rho_s = 1.0, porosity = 0.5,
friction = ManningFriction(n = 0.01),
sediment_model = MeyerPeterMueller(theta_c = 0.0,
d_s = 1e-3))

initial_condition = initial_condition_convergence_test

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

volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)

solver = DGSEM(polydeg = 4,
surface_flux = surface_flux,
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

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

coordinates_min = 0.0
coordinates_max = sqrt(2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 2,
n_cells_max = 10_000,
periodicity = true)

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

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

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

summary_callback = SummaryCallback()

analysis_interval = 200
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

stepsize_callback = StepsizeCallback(cfl = 1.0)

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

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

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

sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
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
104 changes: 104 additions & 0 deletions examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@

using OrdinaryDiffEq
using Trixi
using TrixiShallowWater

###############################################################################
# Semidiscretization of the SWE-Exner equations with a discontinuous sediment bed
# to test well-balancedness

# Equations with Meyer-Peter-Mueller sedimentation model
equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, H0 = 1.0,
rho_f = 0.5, rho_s = 1.0, porosity = 0.5,
friction = ManningFriction(n = 0.01),
sediment_model = MeyerPeterMueller(theta_c = 0.047,
d_s = 1e-3))
patrickersing marked this conversation as resolved.
Show resolved Hide resolved

function initial_condition_steady_state(x, t,

Check warning on line 17 in examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl#L17

Added line #L17 was not covered by tests
equations::ShallowWaterExnerEquations1D)
hv = 0.0

Check warning on line 19 in examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl#L19

Added line #L19 was not covered by tests

# discontinuous sediment bed
if -2.0 < x[1] < 0.0
h_s = 0.3

Check warning on line 23 in examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl#L22-L23

Added lines #L22 - L23 were not covered by tests
else
h_s = 0.1

Check warning on line 25 in examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl#L25

Added line #L25 was not covered by tests
end

h = 1.0 - h_s

Check warning on line 28 in examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl#L28

Added line #L28 was not covered by tests

return SVector(h, hv, h_s)

Check warning on line 30 in examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl

View check run for this annotation

Codecov / codecov/patch

examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl#L30

Added line #L30 was not covered by tests
end

initial_condition = initial_condition_steady_state

boundary_condition = boundary_condition_slip_wall

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

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

basis = LobattoLegendreBasis(3)

indicator_sc = IndicatorHennemannGassner(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = water_sediment_height)
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

coordinates_min = -2.0
coordinates_max = 2.0
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 4,
n_cells_max = 10_000,
periodicity = false)

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

###############################################################################
# ODE solver

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

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

summary_callback = SummaryCallback()

analysis_interval = 10000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
extra_analysis_integrals = (lake_at_rest_error,))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

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

stepsize_callback = StepsizeCallback(cfl = 1.0)

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

###############################################################################
# run the simulation
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
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
Loading
Loading