From 3205c8fed7acdb190f76d333fc07723a040cbecc Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 11 Apr 2024 12:02:04 +0200 Subject: [PATCH 01/17] add multilayer SWE in 2D --- ...xir_shallowwater_multilayer_convergence.jl | 61 ++ ...lixir_shallowwater_multilayer_dam_break.jl | 92 ++ ...r_shallowwater_multilayer_well_balanced.jl | 79 ++ ...xir_shallowwater_multilayer_convergence.jl | 62 ++ ...lixir_shallowwater_multilayer_dam_break.jl | 97 ++ ...r_shallowwater_multilayer_well_balanced.jl | 77 ++ src/TrixiShallowWater.jl | 2 +- src/equations/equations.jl | 1 + src/equations/shallow_water_multilayer_2d.jl | 829 ++++++++++++++++++ test/test_tree_2d.jl | 248 ++++++ test/test_unit.jl | 34 +- test/test_unstructured_2d.jl | 160 ++++ 12 files changed, 1739 insertions(+), 3 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl create mode 100644 src/equations/shallow_water_multilayer_2d.jl diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl new file mode 100644 index 0000000..391eef6 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -0,0 +1,61 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with three layers + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 10.0, + rhos = (0.9, 1.0, 1.1)) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.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, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +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, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8, + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl new file mode 100644 index 0000000..2dbb108 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -0,0 +1,92 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations for a dam break test with a +# discontinuous bottom topography function to test entropy conservation + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0, + rhos = (0.9, 0.95, 1.0)) + +# This academic testcase sets up a discontinuous bottom topography +# function and initial condition to test entropy conservation. + +function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 0.3 * exp(-0.5 * ((x[1])^2 + (x[2])^2)) + + if x[1] < 0.0 + H = SVector(1.0, 0.8, 0.6) + else + H = SVector(0.9, 0.7, 0.5) + b += 0.1 + end + + v1 = zero(H) + v2 = zero(H) + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_dam_break + +boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break) + +############################################################################### +# 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 = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.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, 2.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_total, + energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + 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 diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl new file mode 100644 index 0000000..1a944e1 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -0,0 +1,79 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a bottom topography function +# to test well-balancedness + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, + rhos = (0.7, 0.8, 0.9, 1.0)) + +# An initial condition with constant total water height, zero velocities and a bottom topography to +# test well-balancedness +function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations2D) + H = SVector(0.6, 0.55, 0.5, 0.45) + v1 = zero(H) + v2 = zero(H) + b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + + -0.5)^2)) + 1) : 0.0) + + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_well_balanced + +############################################################################### +# 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 = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +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, + 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 diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl new file mode 100644 index 0000000..708e7a2 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -0,0 +1,62 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a periodic +# bottom topography function (set in the initial conditions) + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 10.0, + rhos = (0.9, 1.0, 1.1)) + +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 = 8, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form convergence test on a periodic domain + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semidiscretization 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 = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +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, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8, + save_everystep = false, callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl new file mode 100644 index 0000000..2545288 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -0,0 +1,97 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations for a dam break test with a +# discontinuous bottom topography function to test entropy conservation + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0, + rhos = (0.9, 0.95, 1.0)) + +# This academic testcase sets up a discontinuous bottom topography +# function and initial condition to test entropy conservation. + +function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 0.3 * exp(-0.5 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2)) + + if x[1] < sqrt(2) / 2 + H = SVector(1.0, 0.8, 0.6) + else + H = SVector(0.9, 0.7, 0.5) + b += 0.1 + end + + v1 = zero(H) + v2 = zero(H) + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_dam_break + +boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break) + +############################################################################### +# 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 = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = false) + +# Boundary conditions +boundary_condition = Dict(:Top => boundary_condition_slip_wall, + :Left => boundary_condition_slip_wall, + :Right => boundary_condition_slip_wall, + :Bottom => boundary_condition_slip_wall) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, + solver, boundary_conditions = boundary_condition) + +############################################################################### +# ODE solver + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_total, + energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + 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 diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl new file mode 100644 index 0000000..4c4763b --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -0,0 +1,77 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a bottom topography function +# to test well-balancedness + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, + rhos = (0.7, 0.8, 0.9, 1.0)) + +# An initial condition with constant total water height, zero velocities and a bottom topography to +# test well-balancedness +function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations2D) + H = SVector(0.6, 0.55, 0.5, 0.45) + v1 = zero(H) + v2 = zero(H) + b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + + -0.5)^2)) + 1) : 0.0) + + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_well_balanced + +############################################################################### +# 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 = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form convergence test on a periodic domain + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +############################################################################### +# ODE solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +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, + 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 diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index ac74b97..c3b0147 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -20,7 +20,7 @@ include("solvers/indicators.jl") # Export types/functions that define the public API of TrixiShallowWater.jl export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, - ShallowWaterMultiLayerEquations1D + ShallowWaterMultiLayerEquations1D, ShallowWaterMultiLayerEquations2D export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle, flux_hll_chen_noelle, diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a36f993..89713b9 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -17,6 +17,7 @@ include("shallow_water_two_layer_2d.jl") abstract type AbstractShallowWaterMultiLayerEquations{NDIMS, NVARS, NLAYERS} <: Trixi.AbstractEquations{NDIMS, NVARS} end include("shallow_water_multilayer_1d.jl") +include("shallow_water_multilayer_2d.jl") """ eachlayer(equations::AbstractShallowWaterMultiLayerEquations) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl new file mode 100644 index 0000000..73dea6d --- /dev/null +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -0,0 +1,829 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + ShallowWaterMultiLayerEquations2D(gravity, H0, rhos) + +Multi-Layer Shallow Water equations (MLSWE) in two space dimension. The equations are given by +```math +\left\{ + \begin{aligned} + &\partial_t h_m + \partial_x h_mv_m = 0,\\ + &\partial h_mv1_m + \partial_x h_mv1_m^2 + \partial_y h_mv1_mv2_m = -gh_m\partial_x \bigg(b + \sum\limits_{k\geq j}h_k + \sum\limits_{k\ + ISBN: 978-0-12-088759-0 +""" + +struct ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT <: Real} <: + AbstractShallowWaterMultiLayerEquations{2, NVARS, NLAYERS} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + rhos::SVector{NLAYERS, RealT} # Vector of layer densities + + function ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}(gravity::RealT, + H0::RealT, + rhos::SVector{NLAYERS, + RealT}) where { + NVARS, + NLAYERS, + RealT <: + Real + } + # Ensure that layer densities are all positive and in increasing order + issorted(rhos) || + throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])")) + min(rhos...) > 0 || throw(ArgumentError("densities must be positive")) + + new(gravity, H0, rhos) + end +end + +# Allow for flexibility to set the gravitational constant within an elixir depending on the +# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. +# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" +# well-balancedness test cases. +function ShallowWaterMultiLayerEquations2D(; gravity_constant, + H0 = zero(gravity_constant), rhos) + + # Promote all variables to a common type + _rhos = promote(rhos...) + RealT = promote_type(eltype(_rhos), eltype(gravity_constant), eltype(H0)) + __rhos = SVector(map(RealT, _rhos)) + + # Extract number of layers and variables + NLAYERS = length(rhos) + NVARS = 3 * NLAYERS + 1 + + return ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}(gravity_constant, + H0, __rhos) +end + +@inline function Base.real(::ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}) where { + NVARS, + NLAYERS, + RealT <: + Real + } + RealT +end + +Trixi.have_nonconservative_terms(::ShallowWaterMultiLayerEquations2D) = True() +function Trixi.varnames(::typeof(cons2cons), + equations::ShallowWaterMultiLayerEquations2D) + waterheight = ntuple(n -> "h" * string(n), Val(nlayers(equations))) + momentum_1 = ntuple(n -> "h" * string(n) * "_v1", Val(nlayers(equations))) + momentum_2 = ntuple(n -> "h" * string(n) * "_v2", Val(nlayers(equations))) + return (waterheight..., momentum_1..., momentum_2..., "b") +end + +# We use the total layer heights, H = ∑h + b as primitive variables for easier visualization and setting initial +# conditions +function Trixi.varnames(::typeof(cons2prim), + equations::ShallowWaterMultiLayerEquations2D) + total_layer_height = ntuple(n -> "H" * string(n), Val(nlayers(equations))) + velocity_1 = ntuple(n -> "v" * string(n) * "_1", Val(nlayers(equations))) + velocity_2 = ntuple(n -> "v" * string(n) * "_2", Val(nlayers(equations))) + return (total_layer_height..., velocity_1..., velocity_2..., "b") +end + +# TODO: This needs to be redone +# # Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations2D) + +A smooth initial condition for a three-layer configuration used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) (and +[`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function Trixi.initial_condition_convergence_test(x, t, + equations::ShallowWaterMultiLayerEquations2D) + # Check if the equation has been set up with three layers + if nlayers(equations) != 3 + throw(ArgumentError("This initial condition is only valid for three layers")) + end + + # Some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] + ω = 2.0 * pi * sqrt(2.0) + + H3 = 1.5 + 0.1 * cos(ω * x[1] + t) + 0.1 * cos(ω * x[2] + t) + H2 = 2.0 + 0.1 * sin(ω * x[1] + t) + 0.1 * sin(ω * x[2] + t) + H1 = 4.0 + 0.1 * cos(ω * x[1] + t) + 0.1 * cos(ω * x[2] + t) + H = (H1, H2, H3) + v1 = (0.8, 0.8, 0.8) + v2 = (1.0, 1.0, 1.0) + + b = 1.0 + 0.1 * cos(ω * x[1]) + 0.1 * cos(ω * x[2]) + + return prim2cons(SVector(H..., v1..., v2..., b), equations) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterMultiLayerEquations2D) + +Source terms used for convergence tests with a three-layer configuration in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) +in non-periodic domains). +""" +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterMultiLayerEquations2D) + # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because + # this manufactured solution velocity is taken to be constant + ω = 2 * pi * sqrt(2.0) + + du1 = -0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - + 0.1cos(t + x[1] * ω) - 0.1cos(t + x[2] * ω) * ω + + 0.8(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + du2 = 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω) + 0.1cos(t + x[2] * ω) * ω + + 0.8(0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 0.1sin(t + x[2] * ω) * ω + + du3 = -0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + 0.1sin(x[2] * ω) * ω + + 0.8(0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + du4 = 0.8(-0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - + 0.1cos(t + x[1] * ω)) + + 0.8(-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + 0.8^2 * (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + + 10.0(2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - + 0.1sin(t + x[2] * ω) + 0.1cos(t + x[1] * ω)) * + (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + + (2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω + + du5 = 0.8(0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω)) + + 0.8(0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (-0.1sin(t + x[1] * ω) * ω + + 0.9(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω)) + + du6 = 0.8(-0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω)) + + 0.8(0.1sin(x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * + (0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * (-0.1sin(x[1] * ω) * ω + + 0.9 / 1.1 * (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + + 1.0 / 1.1 * (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω)) + + du7 = -0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - + 0.1cos(t + x[1] * ω) - 0.1cos(t + x[2] * ω) * ω + + 0.8(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + 10.0(2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - + 0.1sin(t + x[2] * ω) + 0.1cos(t + x[1] * ω)) * + (-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + (2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω + + du8 = 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω) + 0.1cos(t + x[2] * ω) * ω + + 0.8(0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 0.1sin(t + x[2] * ω) * ω + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (0.9(-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω) + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) + + du9 = -0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + 0.1sin(x[2] * ω) * ω + + 0.8(0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * + (0.9 / 1.1 * (-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + 1.0 / 1.1 * (0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) - + 0.1sin(x[2] * ω) * ω) + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * (0.1sin(x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, zero(eltype(u))) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::ShallowWaterMultiLayerEquations2D) +Create a boundary state by reflecting the normal velocity component and keep +the tangential velocity component unchanged. The boundary water height is taken from +the internal value. +For details see Section 9.2.5 of the book: +- Eleuterio F. Toro (2001) + Shock-Capturing Methods for Free-Surface Shallow Flows + 1st edition + ISBN 0471987662 +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::ShallowWaterMultiLayerEquations2D) + # normalize the outward pointing direction + normal = normal_direction / norm(normal_direction) + + # Extract internal values + h = waterheight(u_inner, equations) + h_v1, h_v2 = momentum(u_inner, equations) + b = u_inner[end] + + # u_boundary = SVector(h..., -hv..., b) + + # compute the normal velocity + u_normal = normal[1] * h_v1 + normal[2] * h_v2 + + # create the "external" boundary solution state + u_boundary = SVector(h..., + (h_v1 - 2.0 * u_normal * normal[1])..., + (h_v2 - 2.0 * u_normal * normal[2])..., + b) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) + + return flux +end + +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::ShallowWaterMultiLayerEquations2D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::ShallowWaterMultiLayerEquations2D) + # Extract internal values + h = waterheight(u_inner, equations) + h_v1, h_v2 = momentum(u_inner, equations) + b = u_inner[end] + + ## get the appropriate normal vector from the orientation + if orientation == 1 + u_boundary = SVector(h..., -h_v1..., h_v2..., b) + else # orientation == 2 + u_boundary = SVector(h..., h_v1..., -h_v2..., b) + end + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + + return flux +end + +# Calculate 1D advective portion of the flux for a single point +# Note, the bottom topography has no flux +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterMultiLayerEquations2D) + + # Extract waterheight and momentum and compute velocities + h_v1, h_v2 = momentum(u, equations) + v1, v2 = velocity(u, equations) + + # Initialize flux vector + f = zero(MVector{3 * nlayers(equations) + 1, real(equations)}) + + # Calculate fluxes in each layer + # The momentum flux simplifies as the pressure is included in the nonconservative term + if orientation == 1 + for i in eachlayer(equations) + f_h = h_v1[i] + f_hv1 = f_h * v1[i] + f_hv2 = f_h * v2[i] + setlayer!(f, f_h, f_hv1, f_hv2, i, equations) + end + else # orientation == 2 + for i in eachlayer(equations) + f_h = h_v2[i] + f_hv1 = f_h * v1[i] + f_hv2 = f_h * v2[i] + setlayer!(f, f_h, f_hv1, f_hv2, i, equations) + end + end + + return SVector(f) +end + +# Calculate 1D flux for a single point in the normal direction +# Note, this directional vector is not normalized and the bottom topography has no flux +@inline function Trixi.flux(u, normal_direction::AbstractVector, + equations::ShallowWaterMultiLayerEquations2D) + # Extract waterheight and momentum and compute velocities + h = waterheight(u, equations) + v1, v2 = velocity(u, equations) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + h_v_normal = h .* v_normal + + # Initialize flux vector + f = zero(MVector{3 * nlayers(equations) + 1, real(equations)}) + + # Calculate fluxes in each layer + # The momentum flux simplifies as the pressure is included in the nonconservative term + for i in eachlayer(equations) + f_h = h_v_normal[i] + f_hv1 = f_h * v1[i] + f_hv2 = f_h * v2[i] + setlayer!(f, f_h, f_hv1, f_hv2, i, equations) + end + + return SVector(f) +end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterMultiLayerEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterMultiLayerEquations2D) + +Non-symmetric path-conservative two-point flux discretizing the nonconservative (source) term +that contains the gradients of the bottom topography and waterheights from the coupling between layers +and the nonconservative pressure formulation [`ShallowWaterMultiLayerEquations2D`](@ref). + +When the bottom topography is nonzero this scheme will be well-balanced when used with [`flux_ersing_etal`](@ref). + +In the two-layer setting this combination is equivalent to the fluxes in: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) +""" +@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterMultiLayerEquations2D) + # Pull the necessary left and right state information + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) + b_rr = u_rr[end] + b_ll = u_ll[end] + + # Compute the jumps + h_jump = h_rr - h_ll + b_jump = b_rr - b_ll + g = equations.gravity + + # Initialize flux vector + f = zero(MVector{3 * nlayers(equations) + 1, real(equations)}) + + # Compute the nonconservative flux in each layer + # where f_hv[i] = g * h[i] * (b + ∑h[k] + ∑σ[k] * h[k])_x and σ[k] = ρ[k] / ρ[i] denotes the + # density ratio of different layers + for i in eachlayer(equations) + f_hv = g * h_ll[i] * b_jump + for j in eachlayer(equations) + if j < i + f_hv += g * h_ll[i] * + (equations.rhos[j] / equations.rhos[i] * h_jump[j]) + else # (i Date: Thu, 11 Apr 2024 13:14:08 +0200 Subject: [PATCH 02/17] fix doc references --- src/equations/shallow_water_multilayer_2d.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index 73dea6d..366f2d9 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -51,7 +51,6 @@ A good introduction for the MLSWE is available in Chapter 12 of the book: \ ISBN: 978-0-12-088759-0 """ - struct ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT <: Real} <: AbstractShallowWaterMultiLayerEquations{2, NVARS, NLAYERS} gravity::RealT # gravitational constant @@ -131,7 +130,7 @@ end A smooth initial condition for a three-layer configuration used for convergence tests in combination with [`source_terms_convergence_test`](@ref) (and -[`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +[`BoundaryConditionDirichlet`](@extref) in non-periodic domains). """ function Trixi.initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations2D) @@ -160,7 +159,7 @@ end Source terms used for convergence tests with a three-layer configuration in combination with [`initial_condition_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) +(and [`BoundaryConditionDirichlet`](@extref) in non-periodic domains). """ @inline function Trixi.source_terms_convergence_test(u, x, t, @@ -494,7 +493,7 @@ When the bottom topography is nonzero this scheme will be well-balanced when use nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain an entropy stable formulation the `surface_flux` can be set as -[`FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()), flux_nonconservative_ersing_etal`](@ref). +`FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()), flux_nonconservative_ersing_etal`. In the two-layer setting this combination is equivalent to the fluxes in: - Patrick Ersing, Andrew R. Winters (2023) From 929881f057303d9e6a09333f90698c29fa73b880 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 11 Apr 2024 13:15:41 +0200 Subject: [PATCH 03/17] fix typo and formatting --- src/equations/shallow_water_multilayer_2d.jl | 2 +- test/test_unstructured_2d.jl | 106 +++++++++++++++++-- 2 files changed, 98 insertions(+), 10 deletions(-) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index 366f2d9..fcf5364 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -20,7 +20,7 @@ Multi-Layer Shallow Water equations (MLSWE) in two space dimension. The equation ``` where ``m = 1, 2, ..., M`` is the layer index and the unknown variables are the water height ``h`` and -the velocities ``v1, v2`` in both spatial dimenions . Furthermore, ``g`` denotes the gravitational +the velocities ``v1, v2`` in both spatial dimensions . Furthermore, ``g`` denotes the gravitational constant, ``b(x)`` the bottom topography and ``\rho_m`` the m-th layer density, that must be chosen such that ``\rho_1 < \rho_2 < ... < \rho_M``, to ensure that different layers are ordered from top to bottom, with increasing density. diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index df79deb..a5ae122 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -431,8 +431,30 @@ end # 2LSWE @trixi_testset "elixir_shallowwater_multilayer_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), - l2 = [1.8656084144845292e-5, 1.917796508073001e-5, 1.7454512418793462e-5, 0.00014048665094953486, 3.1328989570708975e-5, 3.512474430566952e-5, 0.00014804378967592363, 3.2299027874400884e-5, 3.881063343808959e-5, 2.393199280374949e-7], - linf = [0.00013984354089258133, 0.00011014730412189921, 0.00010077985342477058, 0.0011054081397450233, 0.00022137639447733504, 0.0002451045776211691, 0.0009651462415440903, 0.00026431961702122475, 0.0002728818589384785, 1.03643740911874e-6], + l2=[ + 1.8656084144845292e-5, + 1.917796508073001e-5, + 1.7454512418793462e-5, + 0.00014048665094953486, + 3.1328989570708975e-5, + 3.512474430566952e-5, + 0.00014804378967592363, + 3.2299027874400884e-5, + 3.881063343808959e-5, + 2.393199280374949e-7, + ], + linf=[ + 0.00013984354089258133, + 0.00011014730412189921, + 0.00010077985342477058, + 0.0011054081397450233, + 0.00022137639447733504, + 0.0002451045776211691, + 0.0009651462415440903, + 0.00026431961702122475, + 0.0002728818589384785, + 1.03643740911874e-6, + ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -447,8 +469,30 @@ end # 2LSWE @trixi_testset "elixir_shallowwater_multilayer_convergence.jl with flux_lax_friedrichs" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), - l2 = [1.0098144596783888e-5, 8.667598940295975e-6, 7.5615805463112364e-6, 1.8150527998748446e-5, 6.935960840543087e-6, 7.335149149911856e-6, 2.403682140927019e-5, 9.684693895783513e-6, 8.27535107484968e-6, 2.393199280374949e-7], - linf = [6.62853303028399e-5, 5.144519565108974e-5, 3.6540170641030656e-5, 0.00013228821066646468, 4.67241008220709e-5, 4.623162831440819e-5, 0.00015925845769571012, 7.204370895297352e-5, 4.566973688335807e-5, 1.03643740911874e-6], + l2=[ + 1.0098144596783888e-5, + 8.667598940295975e-6, + 7.5615805463112364e-6, + 1.8150527998748446e-5, + 6.935960840543087e-6, + 7.335149149911856e-6, + 2.403682140927019e-5, + 9.684693895783513e-6, + 8.27535107484968e-6, + 2.393199280374949e-7, + ], + linf=[ + 6.62853303028399e-5, + 5.144519565108974e-5, + 3.6540170641030656e-5, + 0.00013228821066646468, + 4.67241008220709e-5, + 4.623162831440819e-5, + 0.00015925845769571012, + 7.204370895297352e-5, + 4.566973688335807e-5, + 1.03643740911874e-6, + ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), tspan=(0.0, 0.25)) @@ -555,8 +599,30 @@ end # 2LSWE @trixi_testset "elixir_shallowwater_multilayer_dam_break.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), - l2 = [0.008264086240458098, 0.008904024684630463, 0.015040884852261934, 0.0068380853981447965, 0.006803190764735809, 0.007681935438533932, 0.00014126104190197007, 0.00013800738645477302, 0.00015518676788954237, 0.0062500008375575705], - linf = [0.0271281524526652, 0.038446397522241216, 0.10189873392468171, 0.02279784901518067, 0.02344053175169211, 0.0384823318069141, 0.0015613916350439184, 0.0015293762865072726, 0.0020984165518842324, 0.05808926980412543], + l2=[ + 0.008264086240458098, + 0.008904024684630463, + 0.015040884852261934, + 0.0068380853981447965, + 0.006803190764735809, + 0.007681935438533932, + 0.00014126104190197007, + 0.00013800738645477302, + 0.00015518676788954237, + 0.0062500008375575705, + ], + linf=[ + 0.0271281524526652, + 0.038446397522241216, + 0.10189873392468171, + 0.02279784901518067, + 0.02344053175169211, + 0.0384823318069141, + 0.0015613916350439184, + 0.0015293762865072726, + 0.0020984165518842324, + 0.05808926980412543, + ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -571,9 +637,31 @@ end # 2LSWE @trixi_testset "elixir_shallowwater_multilayer_dam_break.jl with flux_lax_friedrichs" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), - l2 = [0.008289485483638842, 0.008744140233824497, 0.014162544835509572, 0.0067108677626928565, 0.006699962946263383, 0.0075441143283640315, 0.00012772255623776208, 0.0001243555119128215, 0.00013480820924093037, 0.0062500008375575705], - linf = [0.024574998802901732, 0.024386474982493633, 0.09088941363545519, 0.01828534534858913, 0.017626892238781267, 0.02546243728251761, 0.0007231987972067275, 0.0007131387402676157, 0.0007959906366350036, 0.05808926980412543], - surface_flux=(flux_lax_friedrichs, + l2=[ + 0.008289485483638842, + 0.008744140233824497, + 0.014162544835509572, + 0.0067108677626928565, + 0.006699962946263383, + 0.0075441143283640315, + 0.00012772255623776208, + 0.0001243555119128215, + 0.00013480820924093037, + 0.0062500008375575705, + ], + linf=[ + 0.024574998802901732, + 0.024386474982493633, + 0.09088941363545519, + 0.01828534534858913, + 0.017626892238781267, + 0.02546243728251761, + 0.0007231987972067275, + 0.0007131387402676157, + 0.0007959906366350036, + 0.05808926980412543, + ], + surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations From 33263c9ba0de4a9c0ff597cef5447738418fd454 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 11 Apr 2024 13:19:58 +0200 Subject: [PATCH 04/17] soften tolerance to fix test --- test/test_tree_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index 95f1d87..349d068 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -511,7 +511,8 @@ end # 2LSWE ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), - tspan=(0.0, 0.25)) + tspan=(0.0, 0.25), + atol=1e-11) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 9270816fc503ac9270d29c2952630383ebaad892 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 11 Apr 2024 15:31:22 +0200 Subject: [PATCH 05/17] add wall_bc to dam_break_test to increase coverage --- ...lixir_shallowwater_multilayer_dam_break.jl | 7 +- src/equations/shallow_water_multilayer_2d.jl | 5 +- test/test_tree_2d.jl | 72 +++++++++---------- 3 files changed, 42 insertions(+), 42 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl index 2dbb108..883f404 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -32,7 +32,7 @@ end initial_condition = initial_condition_dam_break -boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break) +boundary_conditions = boundary_condition_slip_wall ############################################################################### # Get the DG approximation space @@ -50,10 +50,11 @@ coordinates_max = (1.0, 1.0) mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 4, n_cells_max = 10_000, - periodicity = true) + periodicity = false) # Create the semi discretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) ############################################################################### # ODE solver diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index fcf5364..e57f7ee 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -123,14 +123,13 @@ function Trixi.varnames(::typeof(cons2prim), return (total_layer_height..., velocity_1..., velocity_2..., "b") end -# TODO: This needs to be redone # # Set initial conditions at physical location `x` for time `t` """ initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations2D) A smooth initial condition for a three-layer configuration used for convergence tests in combination with [`source_terms_convergence_test`](@ref) (and -[`BoundaryConditionDirichlet`](@extref) in non-periodic domains). +[`Trixi.BoundaryConditionDirichlet`](@extref) in non-periodic domains). """ function Trixi.initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations2D) @@ -159,7 +158,7 @@ end Source terms used for convergence tests with a three-layer configuration in combination with [`initial_condition_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet`](@extref) +(and [`Trixi.BoundaryConditionDirichlet`](@extref) in non-periodic domains). """ @inline function Trixi.source_terms_convergence_test(u, x, t, diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index 349d068..666b1eb 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -617,27 +617,27 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), l2=[ - 0.009265145927444818, - 0.009576578670098401, - 0.018038886862235676, - 0.008005440938803041, - 0.007959348712974561, - 0.011311110764427515, - 7.142111268477517e-5, - 7.111110836039663e-5, - 9.79599676952869e-5, + 0.006943714920464853, + 0.007232726823994391, + 0.011985810753673342, + 0.005778397278079132, + 0.005736810674427555, + 0.006946616973320985, + 6.143372810603547e-5, + 6.116048140831271e-5, + 7.85179348120683e-5, 0.005455457843589704, ], linf=[ - 0.026783140584084097, - 0.026443746256489098, - 0.09689330703943652, - 0.02878552200736042, - 0.027788281211946812, - 0.03634279241261087, - 0.0003437550406588945, - 0.0003388927276266049, - 0.0006183912668599117, + 0.0264647222938314, + 0.02644739838686727, + 0.09689908398594305, + 0.021624857501702666, + 0.02083541825119268, + 0.029889274751729353, + 0.0003438446930977006, + 0.00033898478545986335, + 0.0006204362806211161, 0.05596747200337038, ], tspan=(0.0, 0.25)) @@ -655,27 +655,27 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), l2=[ - 0.009274909889999189, - 0.009671864215718903, - 0.018865905911242695, - 0.0078068599882339705, - 0.007804963380095172, - 0.011176901444985472, - 6.699629339166091e-5, - 6.713946069649228e-5, - 9.306402661769737e-5, + 0.006846433642482368, + 0.007157521853802785, + 0.012242149312850049, + 0.00562742302611423, + 0.005619979904755309, + 0.006853160451226712, + 5.8363153191008005e-5, + 5.848194044876578e-5, + 7.584883097507706e-5, 0.005455457843589704, ], linf=[ - 0.028399154794939735, - 0.03232097498575831, - 0.1400666269752741, - 0.023201818785649527, - 0.022579330938837905, - 0.03808019720850861, - 0.00032068413360869893, - 0.0003181092631009146, - 0.00042369446429336467, + 0.026631073359625668, + 0.025615528815443434, + 0.0817112552789852, + 0.020564308208119463, + 0.020040028098637117, + 0.024332298183399693, + 0.0003206729530360014, + 0.00031809818924231675, + 0.00042369787238329603, 0.05596747200337038, ], surface_flux=(flux_lax_friedrichs, From 0e2f1011b556da1c22b8b6a1b07b10957e06020d Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 11 Apr 2024 16:30:03 +0200 Subject: [PATCH 06/17] fix comments, update convergence test values --- ...lixir_shallowwater_multilayer_dam_break.jl | 2 +- ...r_shallowwater_multilayer_well_balanced.jl | 2 +- ...xir_shallowwater_multilayer_convergence.jl | 4 +- ...r_shallowwater_multilayer_well_balanced.jl | 2 +- test/test_unstructured_2d.jl | 80 +++++++++---------- 5 files changed, 45 insertions(+), 45 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl index 883f404..af064ef 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -43,7 +43,7 @@ solver = DGSEM(polydeg = 3, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### -# Get the TreeMesh and setup a periodic mesh +# Get the TreeMesh and setup a non-periodic mesh with wall boundary conditions coordinates_min = (-1.0, -1.0) coordinates_max = (1.0, 1.0) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl index 1a944e1..2600aed 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -4,7 +4,7 @@ using Trixi using TrixiShallowWater ############################################################################### -# Semidiscretization of the two-layer shallow water equations with a bottom topography function +# Semidiscretization of the multilayer shallow water equations with a bottom topography function # to test well-balancedness equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl index 708e7a2..a9662ee 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -4,7 +4,7 @@ using Trixi using TrixiShallowWater ############################################################################### -# Semidiscretization of the two-layer shallow water equations with a periodic +# Semidiscretization of the multilayer shallow water equations with a periodic # bottom topography function (set in the initial conditions) equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 10.0, @@ -17,7 +17,7 @@ initial_condition = initial_condition_convergence_test volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg = 8, surface_flux = surface_flux, +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl index 4c4763b..99161e9 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -4,7 +4,7 @@ using Trixi using TrixiShallowWater ############################################################################### -# Semidiscretization of the two-layer shallow water equations with a bottom topography function +# Semidiscretization of the multilayer shallow water equations with a bottom topography function # to test well-balancedness equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index a5ae122..1976aab 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -432,28 +432,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 1.8656084144845292e-5, - 1.917796508073001e-5, - 1.7454512418793462e-5, - 0.00014048665094953486, - 3.1328989570708975e-5, - 3.512474430566952e-5, - 0.00014804378967592363, - 3.2299027874400884e-5, - 3.881063343808959e-5, - 2.393199280374949e-7, + 0.0004897160984647814, + 0.0005096737100476838, + 0.000386126716413923, + 0.0037553984372061602, + 0.0008470734515158057, + 0.000914391804695984, + 0.0039026859922448687, + 0.0008458945196610379, + 0.0009797369075336638, + 1.4701175107332196e-5, ], linf=[ - 0.00013984354089258133, - 0.00011014730412189921, - 0.00010077985342477058, - 0.0011054081397450233, - 0.00022137639447733504, - 0.0002451045776211691, - 0.0009651462415440903, - 0.00026431961702122475, - 0.0002728818589384785, - 1.03643740911874e-6, + 0.003156020123921799, + 0.0023782353180252236, + 0.002165877463488175, + 0.027549055479562767, + 0.007426673513840409, + 0.0057007815151305374, + 0.03394812114803547, + 0.00528271152947668, + 0.008033142366791368, + 5.280913235339302e-5, ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations @@ -470,28 +470,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 1.0098144596783888e-5, - 8.667598940295975e-6, - 7.5615805463112364e-6, - 1.8150527998748446e-5, - 6.935960840543087e-6, - 7.335149149911856e-6, - 2.403682140927019e-5, - 9.684693895783513e-6, - 8.27535107484968e-6, - 2.393199280374949e-7, + 0.0002738452813944285, + 0.00024728982452184266, + 0.00018186742817211048, + 0.000515064164935216, + 0.00020716566255560428, + 0.0001967259461755342, + 0.0006765281760969105, + 0.00026967321034573937, + 0.0002097997387723772, + 1.4701175107334422e-5, ], linf=[ - 6.62853303028399e-5, - 5.144519565108974e-5, - 3.6540170641030656e-5, - 0.00013228821066646468, - 4.67241008220709e-5, - 4.623162831440819e-5, - 0.00015925845769571012, - 7.204370895297352e-5, - 4.566973688335807e-5, - 1.03643740911874e-6, + 0.0016143270108415209, + 0.0013659806938262076, + 0.0008667976189581372, + 0.0024878154391543283, + 0.0007750025239923741, + 0.0012310851890787178, + 0.0031150595992053276, + 0.0015746347140027095, + 0.0011906648185069368, + 5.280913235350404e-5, ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), From 2eec19df7e184318796af5c0c14f2be8673918ab Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 18 Apr 2024 09:42:27 +0200 Subject: [PATCH 07/17] address comment changes from code review --- src/equations/shallow_water_multilayer_1d.jl | 1 + src/equations/shallow_water_multilayer_2d.jl | 7 ++++--- src/equations/shallow_water_two_layer_1d.jl | 1 + src/equations/shallow_water_two_layer_2d.jl | 1 + src/equations/shallow_water_wet_dry_1d.jl | 1 + src/equations/shallow_water_wet_dry_2d.jl | 1 + 6 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 328ea85..0b391b3 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -104,6 +104,7 @@ end end Trixi.have_nonconservative_terms(::ShallowWaterMultiLayerEquations1D) = True() + function Trixi.varnames(::typeof(cons2cons), equations::ShallowWaterMultiLayerEquations1D) waterheight = ntuple(n -> "h" * string(n), Val(nlayers(equations))) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index e57f7ee..16dd026 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -105,6 +105,7 @@ end end Trixi.have_nonconservative_terms(::ShallowWaterMultiLayerEquations2D) = True() + function Trixi.varnames(::typeof(cons2cons), equations::ShallowWaterMultiLayerEquations2D) waterheight = ntuple(n -> "h" * string(n), Val(nlayers(equations))) @@ -123,7 +124,7 @@ function Trixi.varnames(::typeof(cons2prim), return (total_layer_height..., velocity_1..., velocity_2..., "b") end -# # Set initial conditions at physical location `x` for time `t` +# Set initial conditions at physical location `x` for time `t` """ initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations2D) @@ -324,7 +325,7 @@ Should be used together with [`TreeMesh`](@ref). return flux end -# Calculate 1D advective portion of the flux for a single point +# Calculate 2D advective portion of the flux for a single point # Note, the bottom topography has no flux @inline function Trixi.flux(u, orientation::Integer, equations::ShallowWaterMultiLayerEquations2D) @@ -357,7 +358,7 @@ end return SVector(f) end -# Calculate 1D flux for a single point in the normal direction +# Calculate 2D flux for a single point in the normal direction # Note, this directional vector is not normalized and the bottom topography has no flux @inline function Trixi.flux(u, normal_direction::AbstractVector, equations::ShallowWaterMultiLayerEquations2D) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index b3c2b0e..7e6276c 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -84,6 +84,7 @@ function ShallowWaterTwoLayerEquations1D(; gravity_constant, end Trixi.have_nonconservative_terms(::ShallowWaterTwoLayerEquations1D) = True() + function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations1D) ("h_upper", "h_v_upper", "h_lower", "h_v_lower", "b") diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index feb7952..7b8fd6f 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -98,6 +98,7 @@ function ShallowWaterTwoLayerEquations2D(; gravity_constant, end Trixi.have_nonconservative_terms(::ShallowWaterTwoLayerEquations2D) = True() + function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations2D) ("h_upper", "h_v1_upper", "h_v2_upper", "h_lower", "h_v1_lower", "h_v2_lower", "b") end diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index c395786..bb5e726 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -91,6 +91,7 @@ function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_con end Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry1D) = True() + function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquationsWetDry1D) ("h", "h_v", "b") end diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index 32c7e81..2d59294 100644 --- a/src/equations/shallow_water_wet_dry_2d.jl +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -94,6 +94,7 @@ function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_con end Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry2D) = True() + Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquationsWetDry2D) = ("h", "h_v1", "h_v2", "b") # Note, we use the total water height, H = h + b, as the first primitive variable for easier From db5f66057376d0da6aed06cea5cf94b80760c8d0 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 22 Apr 2024 13:06:55 +0200 Subject: [PATCH 08/17] set true discontinuities in 2D --- ...lixir_shallowwater_multilayer_dam_break.jl | 41 ++++++++++ ...lixir_shallowwater_multilayer_dam_break.jl | 49 +++++++++++- test/test_tree_2d.jl | 80 +++++++++---------- test/test_unstructured_2d.jl | 80 +++++++++---------- 4 files changed, 168 insertions(+), 82 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl index af064ef..ad61a69 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -60,6 +60,47 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) +############################################################################### +#= +Workaround for TreeMesh2D to set true discontinuities for debugging and testing. +Essentially, this is a slight augmentation of the `compute_coefficients` where the `x` node values +passed here are slightly perturbed in order to set a true discontinuity that avoids the doubled +value of the LGL nodes at a particular element interface. +=# + +# Point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# Reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for i in eachnode(semi.solver), j in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, + semi.solver, i, j, element) + # Changing the node positions passed to the initial condition by the minimum + # amount possible with the current type of floating point numbers allows setting + # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump` + # works if the jump location `x_jump` is at the position of an interface. + if i == 1 && j == 1 # bottom left corner + x_node = SVector(nextfloat(x_node[1]), nextfloat(x_node[2])) + elseif i == 1 && j == nnodes(semi.solver) # top left corner + x_node = SVector(nextfloat(x_node[1]), prevfloat(x_node[2])) + elseif i == nnodes(semi.solver) && j == 1 # bottom right corner + x_node = SVector(prevfloat(x_node[1]), nextfloat(x_node[2])) + elseif i == nnodes(semi.solver) && j == nnodes(semi.solver) # top right corner + x_node = SVector(prevfloat(x_node[1]), prevfloat(x_node[2])) + elseif i == 1 # left boundary + x_node = SVector(nextfloat(x_node[1]), x_node[2]) + elseif j == 1 # bottom boundary + x_node = SVector(x_node[1], nextfloat(x_node[2])) + elseif i == nnodes(semi.solver) # right boundary + x_node = SVector(prevfloat(x_node[1]), x_node[2]) + elseif j == nnodes(semi.solver) # top boundary + x_node = SVector(x_node[1], prevfloat(x_node[2])) + end + + u_node = initial_condition_dam_break(x_node, first(tspan), equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end ############################################################################### # Callbacks diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl index 2545288..4f1f1bb 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -10,8 +10,10 @@ using TrixiShallowWater equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0, rhos = (0.9, 0.95, 1.0)) -# This academic testcase sets up a discontinuous bottom topography -# function and initial condition to test entropy conservation. +# This test case uses a special work around to setup a truly discontinuous bottom topography +# function and initial condition for this academic testcase of entropy conservation. First, a +# dummy initial_condition_dam_break is introduced to create the semidiscretization. Then the initial +# condition is reset with the true discontinuous values from initial_condition_discontinuous_dam_break. function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D) # Bottom topography @@ -65,6 +67,49 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, tspan = (0.0, 0.5) ode = semidiscretize(semi, tspan) +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# test case and initial condition. +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! + +function initial_condition_discontinuous_dam_break(x, t, element_id, + equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 0.3 * exp(-0.5 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2)) + + # Left side of discontinuity + IDs = [1, 2, 5, 6, 9, 10, 13, 14] + if element_id in IDs + H = SVector(1.0, 0.8, 0.6) + # Right side of discontinuity + else + H = SVector(0.9, 0.7, 0.5) + b += 0.1 + end + + v1 = zero(H) + v2 = zero(H) + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for j in eachnode(semi.solver), i in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, + semi.solver, i, j, element) + u_node = initial_condition_discontinuous_dam_break(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + ############################################################################### # Callbacks diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index 666b1eb..ff0a752 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -617,28 +617,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), l2=[ - 0.006943714920464853, - 0.007232726823994391, - 0.011985810753673342, - 0.005778397278079132, - 0.005736810674427555, - 0.006946616973320985, - 6.143372810603547e-5, - 6.116048140831271e-5, - 7.85179348120683e-5, - 0.005455457843589704, + 0.006867196239306527, + 0.007212095601190101, + 0.014309127811972761, + 0.005801387118671292, + 0.0057555558162007536, + 0.006951246408270715, + 6.084817535430214e-5, + 6.068733435363927e-5, + 7.807791370377477e-5, + 0.003857583749542185, ], linf=[ - 0.0264647222938314, - 0.02644739838686727, - 0.09689908398594305, - 0.021624857501702666, - 0.02083541825119268, - 0.029889274751729353, - 0.0003438446930977006, - 0.00033898478545986335, - 0.0006204362806211161, - 0.05596747200337038, + 0.028792046975677804, + 0.036503344705121676, + 0.19267105917517297, + 0.028685238833612656, + 0.02761561086221269, + 0.029892140357308587, + 0.0004158435510054702, + 0.0004119057098728111, + 0.0007709896609887244, + 0.10000011323773067, ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations @@ -655,28 +655,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), l2=[ - 0.006846433642482368, - 0.007157521853802785, - 0.012242149312850049, - 0.00562742302611423, - 0.005619979904755309, - 0.006853160451226712, - 5.8363153191008005e-5, - 5.848194044876578e-5, - 7.584883097507706e-5, - 0.005455457843589704, + 0.007162130568658371, + 0.007507487998680027, + 0.01500636533827742, + 0.005593279150854283, + 0.005603720137202845, + 0.0068519457928639385, + 5.748086634734491e-5, + 5.7787474070077486e-5, + 7.492163249697633e-5, + 0.003857583749542185, ], linf=[ - 0.026631073359625668, - 0.025615528815443434, - 0.0817112552789852, - 0.020564308208119463, - 0.020040028098637117, - 0.024332298183399693, - 0.0003206729530360014, - 0.00031809818924231675, - 0.00042369787238329603, - 0.05596747200337038, + 0.037711432402825124, + 0.0427565001934887, + 0.17275386567546575, + 0.025348821041845597, + 0.02434895938442919, + 0.03396954814003092, + 0.0003026537474088657, + 0.00029947595638190504, + 0.00037146307629417057, + 0.10000011323773067, ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 1976aab..758b0ff 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -600,28 +600,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), l2=[ - 0.008264086240458098, - 0.008904024684630463, - 0.015040884852261934, - 0.0068380853981447965, - 0.006803190764735809, - 0.007681935438533932, - 0.00014126104190197007, - 0.00013800738645477302, - 0.00015518676788954237, - 0.0062500008375575705, + 0.008133096076522545, + 0.009059379390936938, + 0.0156472431836693, + 0.006787446774526842, + 0.0067495775839911685, + 0.007737394569507473, + 0.0001835435993951076, + 0.00018347939156198144, + 0.00023500246849234892, + 0.004003203849568449, ], linf=[ - 0.0271281524526652, - 0.038446397522241216, - 0.10189873392468171, - 0.02279784901518067, - 0.02344053175169211, - 0.0384823318069141, - 0.0015613916350439184, - 0.0015293762865072726, - 0.0020984165518842324, - 0.05808926980412543, + 0.0268742067416933, + 0.05243441379266703, + 0.1590063807340898, + 0.028581575095277572, + 0.024486565205936787, + 0.03953111184090555, + 0.002625120329243771, + 0.0026387938388929785, + 0.0038007822188701103, + 0.10000000026183736, ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations @@ -638,28 +638,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), l2=[ - 0.008289485483638842, - 0.008744140233824497, - 0.014162544835509572, - 0.0067108677626928565, - 0.006699962946263383, - 0.0075441143283640315, - 0.00012772255623776208, - 0.0001243555119128215, - 0.00013480820924093037, - 0.0062500008375575705, + 0.0089213866827554, + 0.009502261925352482, + 0.017690193299797062, + 0.006845800468099525, + 0.006849466485572439, + 0.00761550757687342, + 0.0004491107702429648, + 0.00043570105077582146, + 0.0005832036370132536, + 0.004003203849568449, ], linf=[ - 0.024574998802901732, - 0.024386474982493633, - 0.09088941363545519, - 0.01828534534858913, - 0.017626892238781267, - 0.02546243728251761, - 0.0007231987972067275, - 0.0007131387402676157, - 0.0007959906366350036, - 0.05808926980412543, + 0.04276512135712615, + 0.054061877851283885, + 0.18896205511083025, + 0.018902398329971263, + 0.017901194463979406, + 0.025135202096835247, + 0.004065137704454743, + 0.0038290081477974037, + 0.005930216968236107, + 0.10000000026183736, ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), From 6fc27e2dab0532d44c931e57f2ba2fa6fb56f36b Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Thu, 25 Apr 2024 10:45:02 +0200 Subject: [PATCH 09/17] Apply suggestions from code review remove implicit multiplication; fix comment Co-authored-by: Andrew Winters --- ...lixir_shallowwater_multilayer_dam_break.jl | 3 +- src/equations/shallow_water_multilayer_2d.jl | 158 +++++++++--------- 2 files changed, 80 insertions(+), 81 deletions(-) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl index 4f1f1bb..5fb3f91 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -85,8 +85,7 @@ function initial_condition_discontinuous_dam_break(x, t, element_id, IDs = [1, 2, 5, 6, 9, 10, 13, 14] if element_id in IDs H = SVector(1.0, 0.8, 0.6) - # Right side of discontinuity - else + else # Right side of discontinuity H = SVector(0.9, 0.7, 0.5) b += 0.1 end diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index 16dd026..ab3009c 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -168,85 +168,85 @@ in non-periodic domains). # this manufactured solution velocity is taken to be constant ω = 2 * pi * sqrt(2.0) - du1 = -0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - - 0.1cos(t + x[1] * ω) - 0.1cos(t + x[2] * ω) * ω + - 0.8(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) - - 0.1sin(t + x[2] * ω) * ω - - du2 = 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + - 0.1cos(t + x[1] * ω) + 0.1cos(t + x[2] * ω) * ω + - 0.8(0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + - 0.1sin(t + x[2] * ω) * ω - - du3 = -0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + 0.1sin(x[2] * ω) * ω + - 0.8(0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) - - 0.1sin(t + x[2] * ω) * ω - - du4 = 0.8(-0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - - 0.1cos(t + x[1] * ω)) + - 0.8(-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + - 0.8^2 * (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + - 10.0(2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - - 0.1sin(t + x[2] * ω) + 0.1cos(t + x[1] * ω)) * - (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + - (2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + - 0.1cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω - - du5 = 0.8(0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + - 0.1cos(t + x[1] * ω)) + - 0.8(0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) + - 0.8^2 * (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + - 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + - 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * - (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + - 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + - 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * - (-0.1sin(t + x[1] * ω) * ω + - 0.9(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω)) - - du6 = 0.8(-0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω)) + - 0.8(0.1sin(x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + - 0.8^2 * (0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) + - 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + - 0.1cos(t + x[1] * ω)) * - (0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) + - 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + - 0.1cos(t + x[1] * ω)) * (-0.1sin(x[1] * ω) * ω + - 0.9 / 1.1 * (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + - 1.0 / 1.1 * (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω)) - - du7 = -0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - - 0.1cos(t + x[1] * ω) - 0.1cos(t + x[2] * ω) * ω + - 0.8(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) - - 0.1sin(t + x[2] * ω) * ω + - 10.0(2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - - 0.1sin(t + x[2] * ω) + 0.1cos(t + x[1] * ω)) * - (-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + - (2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + - 0.1cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω - - du8 = 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + - 0.1cos(t + x[1] * ω) + 0.1cos(t + x[2] * ω) * ω + - 0.8(0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + - 0.1sin(t + x[2] * ω) * ω + - 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + - 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * - (0.9(-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) - - 0.1sin(t + x[2] * ω) * ω) + - 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + - 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * - (0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) - - du9 = -0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + 0.1sin(x[2] * ω) * ω + - 0.8(0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) - - 0.1sin(t + x[2] * ω) * ω + - 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + - 0.1cos(t + x[1] * ω)) * - (0.9 / 1.1 * (-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + - 1.0 / 1.1 * (0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) - - 0.1sin(x[2] * ω) * ω) + - 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + - 0.1cos(t + x[1] * ω)) * (0.1sin(x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + du1 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω) + + du2 = (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 0.1 * sin(t + x[2] * ω) * ω) + + du3 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + 0.1 * sin(x[2] * ω) * ω + + 0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω) + + du4 = (0.8 * (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω)) + + 0.8 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + 0.8^2 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + 10.0 * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω) + + du5 = (0.8 * (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) + + 0.8 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(t + x[1] * ω) * ω + + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω))) + + du6 = (0.8 * (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω)) + + 0.8 * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(x[1] * ω) * ω + + 0.9 / 1.1 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + 1.0 / 1.1 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω))) + + du7 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω + + 10.0 * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω) + + du8 = (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 0.1 * sin(t + x[2] * ω) * ω + + 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.9 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω)) + + du9 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + 0.1 * sin(x[2] * ω) * ω + + 0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω + + 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * + (0.9 / 1.1 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + 1.0 / 1.1 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) - + 0.1 * sin(x[2] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω)) return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, zero(eltype(u))) end From 592a91a7e9d6cf56a98d16ee7b52a397fed004b6 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 25 Apr 2024 13:26:10 +0200 Subject: [PATCH 10/17] soften tolerances for macOS --- test/test_structured_2d.jl | 3 ++- test/test_tree_1d.jl | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index fac1fe6..abb0d92 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -83,7 +83,8 @@ isdir(outdir) && rm(outdir, recursive = true) 4.4849667259339357e-14, 1.9999999999999993, ], - tspan=(0.0, 0.25)) + tspan=(0.0, 0.25), + atol = 1e-12) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 2d19361..45838a6 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -145,7 +145,7 @@ isdir(outdir) && rm(outdir, recursive = true) ], tspan=(0.0, 0.25), # Soften the tolerance as test results vary between different CPUs - atol=1000 * eps()) + atol=2000 * eps()) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From eb55c008eaff63b0d49026e1157cac842c426a21 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 25 Apr 2024 13:27:36 +0200 Subject: [PATCH 11/17] apply formatter --- src/equations/shallow_water_multilayer_2d.jl | 166 +++++++++++-------- 1 file changed, 93 insertions(+), 73 deletions(-) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index ab3009c..61dfaad 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -168,85 +168,105 @@ in non-periodic domains). # this manufactured solution velocity is taken to be constant ω = 2 * pi * sqrt(2.0) - du1 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) - - 0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω + - 0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) - - 0.1 * sin(t + x[2] * ω) * ω) + du1 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω) du2 = (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) + - 0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω + - 0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + - 0.1 * sin(t + x[2] * ω) * ω) - - du3 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + 0.1 * sin(x[2] * ω) * ω + - 0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) - - 0.1 * sin(t + x[2] * ω) * ω) - - du4 = (0.8 * (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) - - 0.1 * cos(t + x[1] * ω)) + - 0.8 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + - 0.8^2 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + - 10.0 * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - - 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * - (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + - (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω) - - du5 = (0.8 * (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) + - 0.8 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) + - 0.8^2 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * - (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * - (-0.1 * sin(t + x[1] * ω) * ω + - 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω))) + 0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 0.1 * sin(t + x[2] * ω) * ω) + + du3 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + + 0.1 * sin(x[2] * ω) * ω + + 0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω) + + du4 = (0.8 * (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω)) + + 0.8 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + 0.8^2 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + 10.0 * + (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω) + + du5 = (0.8 * (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) + + 0.8 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 10.0 * + (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 10.0 * + (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(t + x[1] * ω) * ω + + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω))) du6 = (0.8 * (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω)) + - 0.8 * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + - 0.8^2 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) * - (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(x[1] * ω) * ω + - 0.9 / 1.1 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + - 1.0 / 1.1 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω))) - - du7 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) - - 0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω + - 0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) - - 0.1 * sin(t + x[2] * ω) * ω + - 10.0 * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - - 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * - (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + - (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω) + 0.8 * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + + 10.0 * + (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + + 10.0 * + (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(x[1] * ω) * ω + + 0.9 / 1.1 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + 1.0 / 1.1 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω))) + + du7 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω + + 10.0 * + (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω) du8 = (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) + - 0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω + - 0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + - 0.1 * sin(t + x[2] * ω) * ω + - 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * - (0.9 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) - - 0.1 * sin(t + x[2] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * - (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω)) - - du9 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + 0.1 * sin(x[2] * ω) * ω + - 0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) - - 0.1 * sin(t + x[2] * ω) * ω + - 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) * - (0.9 / 1.1 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + - 1.0 / 1.1 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) - - 0.1 * sin(x[2] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω)) + 0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω + + 0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 0.1 * sin(t + x[2] * ω) * ω + + 10.0 * + (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.9 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω) + + 10.0 * + (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω)) + + du9 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + + 0.1 * sin(x[2] * ω) * ω + + 0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω + + 10.0 * + (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * + (0.9 / 1.1 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + 1.0 / 1.1 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) - + 0.1 * sin(x[2] * ω) * ω) + + 10.0 * + (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω)) return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, zero(eltype(u))) end From 5145e7e627c719d5dc95ce66cf2466fab857d8a6 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 25 Apr 2024 13:30:44 +0200 Subject: [PATCH 12/17] format again --- test/test_structured_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index abb0d92..eff0178 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -84,7 +84,7 @@ isdir(outdir) && rm(outdir, recursive = true) 1.9999999999999993, ], tspan=(0.0, 0.25), - atol = 1e-12) + atol=1e-12) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 7f00ff711295156813d77f36c85f8a43e0817643 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 26 Apr 2024 11:31:33 +0200 Subject: [PATCH 13/17] update 2D convergence test --- ...xir_shallowwater_multilayer_convergence.jl | 14 ++-- ...xir_shallowwater_multilayer_convergence.jl | 9 ++- src/equations/shallow_water_multilayer_2d.jl | 61 ++++++++------ test/test_tree_2d.jl | 80 +++++++++---------- test/test_unstructured_2d.jl | 80 +++++++++---------- 5 files changed, 130 insertions(+), 114 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl index 391eef6..23252fa 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -6,7 +6,7 @@ using TrixiShallowWater ############################################################################### # Semidiscretization of the multilayer shallow water equations with three layers -equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 10.0, +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.1, rhos = (0.9, 1.0, 1.1)) initial_condition = initial_condition_convergence_test @@ -25,8 +25,8 @@ solver = DGSEM(polydeg = 3, coordinates_min = (0.0, 0.0) coordinates_max = (sqrt(2.0), sqrt(2.0)) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 4, - n_cells_max = 10_000, + initial_refinement_level = 2, + n_cells_max = 100_000, periodicity = true) # create the semi discretization object @@ -50,12 +50,16 @@ save_solution = SaveSolutionCallback(interval = 500, save_initial_solution = true, save_final_solution = true) -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) ############################################################################### # run the simulation # use a Runge-Kutta method with automatic (error based) time step size control -sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8, +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 diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl index a9662ee..08a9283 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -7,7 +7,7 @@ using TrixiShallowWater # Semidiscretization of the multilayer shallow water equations with a periodic # bottom topography function (set in the initial conditions) -equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 10.0, +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.1, rhos = (0.9, 1.0, 1.1)) initial_condition = initial_condition_convergence_test @@ -50,13 +50,16 @@ save_solution = SaveSolutionCallback(interval = 500, save_initial_solution = true, save_final_solution = true) -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) ############################################################################### # run the simulation # use a Runge-Kutta method with automatic (error based) time step size control -sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8, +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 diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index 61dfaad..d9798a1 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -140,7 +140,7 @@ function Trixi.initial_condition_convergence_test(x, t, end # Some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] - ω = 2.0 * pi * sqrt(2.0) + ω = pi * sqrt(2.0) H3 = 1.5 + 0.1 * cos(ω * x[1] + t) + 0.1 * cos(ω * x[2] + t) H2 = 2.0 + 0.1 * sin(ω * x[1] + t) + 0.1 * sin(ω * x[2] + t) @@ -166,7 +166,8 @@ in non-periodic domains). equations::ShallowWaterMultiLayerEquations2D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocity is taken to be constant - ω = 2 * pi * sqrt(2.0) + ω = pi * sqrt(2.0) + g = equations.gravity du1 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) - @@ -189,10 +190,12 @@ in non-periodic domains). 0.1 * cos(t + x[1] * ω)) + 0.8 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + 0.8^2 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + - 10.0 * + g * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - - 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + 0.1 * g * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω) @@ -202,24 +205,26 @@ in non-periodic domains). 0.1 * cos(t + x[1] * ω)) + 0.8 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) + 0.8^2 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + - 10.0 * + g * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * - (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + - 10.0 * - (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * + 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(t + x[1] * ω) * ω + - 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω))) + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) + + g * + (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + + 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω)) du6 = (0.8 * (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω)) + 0.8 * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + 0.8^2 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + - 10.0 * + g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + - 10.0 * + g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(x[1] * ω) * ω + @@ -231,39 +236,43 @@ in non-periodic domains). 0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω + 0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) - 0.1 * sin(t + x[2] * ω) * ω + - 10.0 * + 0.1 * g * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - - 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * - (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + + 0.1 * sin(t + x[2] * ω) + + 0.1 * cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω + + g * (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + - 0.1 * cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω)) du8 = (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) + 0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω + 0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + 0.1 * sin(t + x[2] * ω) * ω + - 10.0 * + g * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * - (0.9 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) - - 0.1 * sin(t + x[2] * ω) * ω) + - 10.0 * + 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω)) * + (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) + + g * (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + - 0.1 * sin(t + x[2] * ω) - 0.1 * cos(t + x[1] * ω)) * - (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω)) + 0.1 * sin(t + x[2] * ω) - + 0.1 * cos(t + x[1] * ω)) * + (0.9 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) - + 0.1 * sin(t + x[2] * ω) * ω)) du9 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) + 0.1 * sin(x[2] * ω) * ω + 0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) - 0.1 * sin(t + x[2] * ω) * ω + - 10.0 * + g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * (0.9 / 1.1 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) + 1.0 / 1.1 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) - 0.1 * sin(x[2] * ω) * ω) + - 10.0 * + g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) + 0.1 * cos(t + x[1] * ω)) * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω)) diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index ff0a752..18a780a 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -448,28 +448,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.00013147471504344405, - 0.00015514765626121514, - 0.000189385179009786, - 0.00011608534902216235, - 0.00015078355278842716, - 0.00014740625262810186, - 0.00011588655436490155, - 0.00018765407379989162, - 0.00016853529382291021, - 1.2548293535443184e-5, + 0.001578480912165503, + 0.0010725791806397446, + 0.00036117546296274316, + 0.0010107766709183063, + 0.0009612673392275463, + 0.00026390576962071974, + 0.001140711809558893, + 0.001197469124790338, + 0.00031384421384170294, + 0.00019675440964325044, ], linf=[ - 0.0007914697417876759, - 0.0006348222705068185, - 0.001002926474891086, - 0.0006466327311540621, - 0.0009859041580843608, - 0.0008083921716079412, - 0.0006016100163761529, - 0.0010877289257850142, - 0.0010435037400048364, - 3.639351911033373e-5, + 0.006611273262461914, + 0.0047987596637860674, + 0.001476163166610922, + 0.004378481797734368, + 0.004059620398096986, + 0.00107253722812789, + 0.004650692998510841, + 0.005014864957408438, + 0.0013175223929244861, + 0.0004374891172380657, ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations @@ -486,28 +486,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.00015249970137499414, - 9.061891979391998e-5, - 9.396591776903845e-5, - 0.00011273701474177127, - 6.279405808220591e-5, - 5.3495033579633576e-5, - 0.00010797255635235481, - 7.167327550380933e-5, - 7.40293792497647e-5, - 1.2548293535443157e-5, + 0.0009657284251663235, + 0.0007574316192008629, + 0.00014884847789419973, + 0.0008365221427534766, + 0.0005829573732731771, + 0.00013759996082561183, + 0.0009795167874289885, + 0.0007208939038902295, + 0.00015863343571957187, + 0.00019675440964325044, ], linf=[ - 0.0007475490650712402, - 0.00040410931043857734, - 0.0004618601405063094, - 0.0005819120302106295, - 0.00033956825032638305, - 0.00016617565116200383, - 0.0005672335445119359, - 0.0003107057954721548, - 0.00038873185787224873, - 3.639351911033373e-5, + 0.005343578840408814, + 0.005546044302341735, + 0.0005110705217063471, + 0.004897047461771331, + 0.004047517471484768, + 0.0005991894625433924, + 0.006191481703377466, + 0.005110806058533313, + 0.000775668423696807, + 0.0004374891172380657, ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 758b0ff..8bdd0bd 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -432,28 +432,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0004897160984647814, - 0.0005096737100476838, - 0.000386126716413923, - 0.0037553984372061602, - 0.0008470734515158057, - 0.000914391804695984, - 0.0039026859922448687, - 0.0008458945196610379, - 0.0009797369075336638, - 1.4701175107332196e-5, + 2.845107067734722e-5, + 2.541332828758267e-5, + 1.3776844294886904e-5, + 3.21641609844227e-5, + 1.9645259635774905e-5, + 1.2494077534998188e-5, + 3.842168759551814e-5, + 2.3940963660476904e-5, + 1.5539403052020932e-5, + 7.244891628173627e-7, ], linf=[ - 0.003156020123921799, - 0.0023782353180252236, - 0.002165877463488175, - 0.027549055479562767, - 0.007426673513840409, - 0.0057007815151305374, - 0.03394812114803547, - 0.00528271152947668, - 0.008033142366791368, - 5.280913235339302e-5, + 0.0001686597682639679, + 0.00016943475945224717, + 0.00010161333860780886, + 0.00018234203945666216, + 0.0001440335843042595, + 7.821014430642315e-5, + 0.0002159224957090089, + 0.00016273539529820802, + 9.798918563769243e-5, + 4.040896422807805e-6, ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations @@ -470,28 +470,28 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0002738452813944285, - 0.00024728982452184266, - 0.00018186742817211048, - 0.000515064164935216, - 0.00020716566255560428, - 0.0001967259461755342, - 0.0006765281760969105, - 0.00026967321034573937, - 0.0002097997387723772, - 1.4701175107334422e-5, + 1.918097575328871e-5, + 1.502019981850744e-5, + 5.042301180852648e-6, + 1.9602948809781437e-5, + 1.1406394406547229e-5, + 5.066224008439183e-6, + 2.3354649343891286e-5, + 1.433469494607502e-5, + 6.046742009426357e-6, + 7.244891628173627e-7, ], linf=[ - 0.0016143270108415209, - 0.0013659806938262076, - 0.0008667976189581372, - 0.0024878154391543283, - 0.0007750025239923741, - 0.0012310851890787178, - 0.0031150595992053276, - 0.0015746347140027095, - 0.0011906648185069368, - 5.280913235350404e-5, + 0.0001267907408812885, + 0.00013376989668978378, + 3.3611013632417475e-5, + 0.00013760680376173617, + 0.00011659018473553218, + 3.15196814599239e-5, + 0.00017067740318088553, + 0.00012830136581865048, + 3.6074428482746335e-5, + 4.040896422807805e-6, ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), From 737ebcc72b9b58b7b1c8c0a22275e75c4d6d4a95 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 26 Apr 2024 11:44:40 +0200 Subject: [PATCH 14/17] update 1D convergence test --- ...xir_shallowwater_multilayer_convergence.jl | 14 +++-- ...xir_shallowwater_multilayer_convergence.jl | 1 - ...xir_shallowwater_multilayer_convergence.jl | 4 +- src/equations/shallow_water_multilayer_1d.jl | 44 ++++++++------- test/test_tree_1d.jl | 56 +++++++++---------- 5 files changed, 63 insertions(+), 56 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl index 8a9384f..05af3f3 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -6,7 +6,7 @@ using TrixiShallowWater ############################################################################### # Semidiscretization of the multilayer shallow water equations with three layers -equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 10.0, +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.1, rhos = (0.9, 1.0, 1.1)) initial_condition = initial_condition_convergence_test @@ -25,7 +25,7 @@ solver = DGSEM(polydeg = 3, coordinates_min = 0.0 coordinates_max = sqrt(2.0) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 4, + initial_refinement_level = 2, n_cells_max = 10_000, periodicity = true) @@ -50,12 +50,16 @@ save_solution = SaveSolutionCallback(interval = 500, save_initial_solution = true, save_final_solution = true) -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) ############################################################################### # run the simulation -# use a Runge-Kutta method with automatic (error based) time step size control -sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8, +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 diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl index 23252fa..7b070fb 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -58,7 +58,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav ############################################################################### # run the simulation -# use a Runge-Kutta method with automatic (error based) time step size control 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); diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl index 08a9283..1a9fec4 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -52,12 +52,12 @@ save_solution = SaveSolutionCallback(interval = 500, stepsize_callback = StepsizeCallback(cfl = 0.9) -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) ############################################################################### # run the simulation -# use a Runge-Kutta method with automatic (error based) time step size control 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); diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 0b391b3..ba6e1ef 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -138,7 +138,7 @@ function Trixi.initial_condition_convergence_test(x, t, end # Some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] - ω = 2.0 * pi * sqrt(2.0) + ω = pi * sqrt(2.0) H1 = 4.0 + 0.1 * cos(ω * x[1] + t) H2 = 2.0 + 0.1 * sin(ω * x[1] + t) @@ -146,7 +146,7 @@ function Trixi.initial_condition_convergence_test(x, t, v1 = 0.9 v2 = 1.0 v3 = 1.1 - b = 1.0 + 0.1 * cos(2.0 * ω * x[1]) + b = 1.0 + 0.1 * cos(ω * x[1]) return prim2cons(SVector(H1, H2, H3, v1, v2, v3, b), equations) end @@ -163,39 +163,43 @@ in non-periodic domains). equations::ShallowWaterMultiLayerEquations1D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocity is taken to be constant - ω = 2 * pi * sqrt(2.0) + ω = pi * sqrt(2.0) + g = equations.gravity du1 = -0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω) + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + du2 = 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω) + 0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω + du3 = -0.1 * sin(t + x[1] * ω) + - 1.1 * (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0 * x[1] * ω) * ω) + 1.1 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + du4 = 0.9 * (-0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) + 0.81 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + - 10.0 * (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + g * (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + - (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + 0.1 * g * (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω + du5 = 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω) + 0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω + - 10.0 * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) * + g * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + g * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(t + x[1] * ω) * ω + - 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) + - 10.0 * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) * - (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) - du6 = -0.11 * sin(t + x[1] * ω) + - 1.21 * (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0x[1] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(2.0 * x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * - (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0 * x[1] * ω) * ω) + - 10.0 * (0.5 - 0.1 * cos(2.0 * x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * - (0.8181818181818181 * - (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + - 0.9090909090909091 * - (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) - - 0.2 * sin(2.0 * x[1] * ω) * ω) + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) + + du6 = -0.11000000000000001 * sin(t + x[1] * ω) + + 1.2100000000000002 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + + g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(x[1] * ω) * ω + + 1.0 / 1.1 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + + 0.9 / 1.1 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) + + g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) return SVector(du1, du2, du3, du4, du5, du6, zero(eltype(u))) end diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 45838a6..34a2585 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -480,22 +480,22 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0024190350732389444, - 0.0011553525849735248, - 0.001427790649735925, - 0.0006301962973580928, - 0.00047816381506058254, - 0.0005518957180338281, - 0.00013912637728693124, + 0.0012190787287581911, + 0.0007900811065315603, + 0.00024265888138436312, + 0.00044695521604641624, + 0.0009421189161161149, + 0.0001681269882698457, + 0.000139126377287091, ], linf=[ - 0.011265033642474886, - 0.005221927831678352, - 0.007355554947688914, - 0.0021140505291696865, - 0.001562655691777548, - 0.002409077213895827, - 0.0002187445586192549, + 0.0035625489595796367, + 0.0025099269565310167, + 0.000677952137713711, + 0.0019817824651813254, + 0.002749046992753801, + 0.0006393958093637853, + 0.00021874455861881081, ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations @@ -512,22 +512,22 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0019005751369889307, - 0.0006416934635148739, - 0.001323748067661235, - 0.0002514848406032196, - 0.00019551914601560465, - 0.00029396199135642574, - 0.00013912637728693124, + 0.0007413378395106811, + 0.0005405760373610748, + 8.772497321219398e-5, + 0.0006969915241438401, + 0.0004924779641660143, + 0.00011816455915431224, + 0.000139126377287091, ], linf=[ - 0.005947927726240643, - 0.0028032845734302647, - 0.003944642659854947, - 0.0009866355929912807, - 0.0009301080749861135, - 0.0012826136652659414, - 0.0002187445586192549, + 0.002657085141947846, + 0.0028102251075707296, + 0.00022855706065982861, + 0.003310637567035979, + 0.002490976099376596, + 0.0005211086991627756, + 0.00021874455861881081, ], surface_flux=(flux_lax_friedrichs, flux_nonconservative_ersing_etal), From ade337759d875a3b3b37ee0b1bc8fc8bb643c096 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 26 Apr 2024 16:05:21 +0200 Subject: [PATCH 15/17] add changes from code review --- src/equations/shallow_water_multilayer_1d.jl | 4 ++-- src/equations/shallow_water_multilayer_2d.jl | 2 -- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index ba6e1ef..138768f 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -192,8 +192,8 @@ in non-periodic domains). (-0.1 * sin(t + x[1] * ω) * ω + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) - du6 = -0.11000000000000001 * sin(t + x[1] * ω) + - 1.2100000000000002 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + + du6 = -0.11 * sin(t + x[1] * ω) + + 1.21 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) + g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * (-0.1 * sin(x[1] * ω) * ω + 1.0 / 1.1 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index d9798a1..da26309 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -305,8 +305,6 @@ For details see Section 9.2.5 of the book: h_v1, h_v2 = momentum(u_inner, equations) b = u_inner[end] - # u_boundary = SVector(h..., -hv..., b) - # compute the normal velocity u_normal = normal[1] * h_v1 + normal[2] * h_v2 From 1b897af4519889eafc17622de1d49678a8080583 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 26 Apr 2024 16:14:24 +0200 Subject: [PATCH 16/17] remove empty line --- .../elixir_shallowwater_multilayer_dam_break.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl index 5fb3f91..5e41ef1 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -75,7 +75,6 @@ ode = semidiscretize(semi, tspan) # In contrast to the usual signature of initial conditions, this one get passed the # `element_id` explicitly. In particular, this initial conditions works as intended # only for the specific mesh loaded above! - function initial_condition_discontinuous_dam_break(x, t, element_id, equations::ShallowWaterMultiLayerEquations2D) # Bottom topography From ccda9560a7c49fb730cb677f5c4c658bed1870c1 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 26 Apr 2024 16:32:11 +0200 Subject: [PATCH 17/17] update formatting of bottom topography function --- .../elixir_shallowwater_multilayer_well_balanced.jl | 4 ++-- .../elixir_shallowwater_multilayer_well_balanced.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl index 2600aed..3069112 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -17,8 +17,8 @@ function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayer v1 = zero(H) v2 = zero(H) b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? - 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + - -0.5)^2)) + 1) : 0.0) + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] - + 0.5)^2)) + 1) : 0.0) return prim2cons(SVector(H..., v1..., v2..., b), equations) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl index 99161e9..986f236 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -17,8 +17,8 @@ function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayer v1 = zero(H) v2 = zero(H) b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? - 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + - -0.5)^2)) + 1) : 0.0) + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] - + 0.5)^2)) + 1) : 0.0) return prim2cons(SVector(H..., v1..., v2..., b), equations)