From eb9606f2292ecbb4240fd0af902479d8e287f0a9 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 15 Aug 2024 16:59:44 +0200 Subject: [PATCH 01/34] variable bathymetry Serre-Green-Naghdi equations --- src/DispersiveShallowWater.jl | 2 +- src/equations/equations.jl | 27 +- src/equations/serre_green_naghdi_1d.jl | 869 ++++++++++++++++++++ src/equations/serre_green_naghdi_flat_1d.jl | 456 ---------- test/test_serre_green_naghdi_1d.jl | 120 +++ 5 files changed, 1016 insertions(+), 458 deletions(-) create mode 100644 src/equations/serre_green_naghdi_1d.jl delete mode 100644 src/equations/serre_green_naghdi_flat_1d.jl diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 761ca562..671acaaa 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -65,7 +65,7 @@ export Semidiscretization, semidiscretize, grid export boundary_condition_periodic, boundary_condition_reflecting -export bathymetry_flat +export bathymetry_flat, bathymetry_mild_slope, bathymetry_variable export initial_condition_convergence_test, initial_condition_manufactured, source_terms_manufactured, diff --git a/src/equations/equations.jl b/src/equations/equations.jl index eba21306..465e8f3a 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -325,14 +325,39 @@ Default analysis integrals used by the [`AnalysisCallback`](@ref). default_analysis_integrals(::AbstractEquations) = Symbol[] abstract type AbstractBathymetry end + struct BathymetryFlat <: AbstractBathymetry end """ bathymetry_flat = DispersiveShallowWater.BathymetryFlat() A singleton struct indicating a flat bathymetry. + +See also [`bathymetry_mild_slope`](@ref) and [`bathymetry_variable`](@ref). """ const bathymetry_flat = BathymetryFlat() +struct BathymetryMildSlope <: AbstractBathymetry end +""" + bathymetry_mild_slope = DispersiveShallowWater.BathymetryMildSlope() + +A singleton struct indicating a variable bathymetry with mild-slope +approximation. + +See also [`bathymetry_flat`](@ref) and [`bathymetry_variable`](@ref). +""" +const bathymetry_mild_slope = BathymetryMildSlope() + +struct BathymetryVariable <: AbstractBathymetry end +""" + bathymetry_variable = DispersiveShallowWater.BathymetryVariable() + +A singleton struct indicating a variable bathymetry (without +mild-slope approximation). + +See also [`bathymetry_flat`](@ref) and [`bathymetry_mild_slope`](@ref). +""" +const bathymetry_variable = BathymetryVariable() + # BBM-BBM equations abstract type AbstractBBMBBMEquations{NDIMS, NVARS} <: AbstractShallowWaterEquations{NDIMS, NVARS} end @@ -347,4 +372,4 @@ include("svaerd_kalisch_1d.jl") # Serre-Green-Naghdi equations abstract type AbstractSerreGreenNaghdiEquations{NDIMS, NVARS} <: AbstractShallowWaterEquations{NDIMS, NVARS} end -include("serre_green_naghdi_flat_1d.jl") +include("serre_green_naghdi_1d.jl") diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl new file mode 100644 index 00000000..4deee588 --- /dev/null +++ b/src/equations/serre_green_naghdi_1d.jl @@ -0,0 +1,869 @@ +@doc raw""" + SerreGreenNaghdiEquations1D(bathymetry = bathymetry_flat; + gravity_constant, eta0 = 0.0) + +Serre-Green-Naghdi system in one spatial dimension. +The equations for flat bathymetry are given by +```math +\begin{aligned} + h_t + (h v)_x &= 0,\\ + h v_t - \frac{1}{3} (h^3 v_{tx})_x + \frac{1}{2} g (h^2)_x + \frac{1}{2} h (v^2)_x + p_x &= 0,\\ + p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx}. +\end{aligned} +``` +The unknown quantities of the Serre-Green-Naghdi equations are the +total water height ``\eta = h + b`` and the velocity ``v``. +The gravitational constant is denoted by `g` and the bottom topography +(bathymetry) ``b = \eta_0 - D``. The water height above the bathymetry +is therefore given by ``h = \eta - \eta_0 + D``. The Serre-Green-Naghdi +equations are only implemented for ``\eta_0 = 0``. +The total water height is therefore given by ``\eta = h + b``. + +Three types of variable `bathymetry` are supported: +- [`bathymetry_flat`](@ref): flat bathymetry (typically ``b = 0`` everywhere) +- [`bathymetry_mild_slope`](@ref): variable bathymetry with mild-slope approximation +- [`bathymetry_variable`](@ref): general variable bathymetry + +For the mild-slope approximation, the Serre-Green-Naghdi equations are +```math +\begin{aligned} + h_t + (h v)_x &= 0,\\ + h v_t - \frac{1}{3} (h^3 v_{tx})_x + \frac{1}{2} (h^2 b_x u_t)_x - \frac{1}{2} h^2 b_x u_{tx} + \frac{3}{4} h b_x^2 u_t + + \frac{1}{2} g (h^2)_x + g h b_x + \frac{1}{2} h (v^2)_x + + p_x + \frac{3}{2} \frac{p}{h} b_x &= 0,\\ + p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx} + + \frac{1}{2} h^2 v (b_x v)_x. +\end{aligned} +``` +For the general case of variable vathymetry without mild-slope +approximation, the Serre-Green-Naghdi equations are +```math +\begin{aligned} + h_t + (h v)_x &= 0,\\ + h v_t - \frac{1}{3} (h^3 v_{tx})_x + \frac{1}{2} (h^2 b_x u_t)_x - \frac{1}{2} h^2 b_x u_{tx} + h b_x^2 u_t + + \frac{1}{2} g (h^2)_x + g h b_x + \frac{1}{2} h (v^2)_x + + p_x + \frac{3}{2} \frac{p}{h} b_x + \psi b_x &= 0,\\ + p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx} + + \frac{1}{2} h^2 v (b_x v)_x, + \psi &= \frac{1}{4} h v (b_x v)_x. +\end{aligned} +``` + +References for the Serre-Green-Naghdi system can be found in +- Serre (1953) + Contribution â l'étude des écoulements permanents et variables dans les canaux + [DOI: 10.1051/lhb/1953034](https://doi.org/10.1051/lhb/1953034) +- Green and Naghdi (1976) + A derivation of equations for wave propagation in water of variable depth + [DOI: 10.1017/S0022112076002425](https://doi.org/10.1017/S0022112076002425) + +The semidiscretization implemented here conserves +- the total water mass (integral of h) as a linear invariant +- the total momentum (integral of h v) as a nonlinear invariant +- the total energy + +for periodic boundary conditions, see +- Hendrik Ranocha and Mario Ricchiuto (2024) + Structure-preserving approximations of the Serre-Green-Naghdi + equations in standard and hyperbolic form + [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) +""" +struct SerreGreenNaghdiEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <: + AbstractSerreGreenNaghdiEquations{1, 3} + bathymetry::Bathymetry # type of bathymetry + gravity::RealT # gravitational constant + eta0::RealT # constant still-water surface +end + +function SerreGreenNaghdiEquations1D(bathymetry = bathymetry_flat; + gravity_constant, eta0 = 0.0) + eta0 == 0.0 || + @warn "The still-water surface needs to be 0 for the Serre-Green-Naghdi equations" + SerreGreenNaghdiEquations1D(bathymetry, gravity_constant, eta0) +end + +varnames(::typeof(prim2prim), ::SerreGreenNaghdiEquations1D) = ("η", "v", "D") +varnames(::typeof(prim2cons), ::SerreGreenNaghdiEquations1D) = ("h", "hv", "b") + +""" + initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEquations1D, mesh) + +A soliton solution used for convergence tests in a periodic domain. +""" +function initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEquations1D, + mesh) + g = equations.gravity + + # setup parameters data + h1 = 1.0 + h2 = 1.2 + c = sqrt(g * h2) + + x_t = mod(x - c * t - xmin(mesh), xmax(mesh) - xmin(mesh)) + xmin(mesh) + + h = h1 + (h2 - h1) * sech(x_t / 2 * sqrt(3 * (h2 - h1) / (h1^2 * h2)))^2 + v = c * (1 - h1 / h) + + return SVector(h, v, 0) +end + +# flat bathymetry +function create_cache(mesh, + equations::SerreGreenNaghdiEquations1D{BathymetryFlat}, + solver, + initial_condition, + ::BoundaryConditionPeriodic, + RealT, uEltype) + D = solver.D1 + + # create temporary storage + h = ones(RealT, nnodes(mesh)) + h_x = zero(h) + v_x = zero(h) + h2_x = zero(h) + hv_x = zero(h) + v2_x = zero(h) + h2_v_vx_x = zero(h) + h_vx_x = zero(h) + p_x = zero(h) + tmp = zero(h) + M_h = zero(h) + M_h3_3 = zero(h) + + if D isa PeriodicUpwindOperators + v_x_upwind = zero(h) + + Dmat_minus = sparse(D.minus) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + @. M_h = h + scale_by_mass_matrix!(M_h, D) + @. M_h3_3 = (1 / 3) * h^3 + scale_by_mass_matrix!(M_h3_3, D) + system_matrix = Symmetric(Diagonal(M_h) + + + Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) + factorization = cholesky(system_matrix) + + cache = (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_x, tmp, + M_h, M_h3_3, + D, Dmat_minus, factorization) + else + if D isa FourierDerivativeOperator + Dmat = Matrix(D) + + cache = (; h_x, v_x, h2_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_x, tmp, + M_h, M_h3_3, + D, Dmat) + else + Dmat = sparse(D) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + @. M_h = h + scale_by_mass_matrix!(M_h, D) + @. M_h3_3 = (1 / 3) * h^3 + scale_by_mass_matrix!(M_h3_3, D) + system_matrix = Symmetric(Diagonal(M_h) + + + Dmat' * Diagonal(M_h3_3) * Dmat) + + factorization = cholesky(system_matrix) + + cache = (; h_x, v_x, h2_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_x, tmp, + M_h, M_h3_3, + D, Dmat, factorization) + end + end + + return cache +end + +# variable bathymetry +function create_cache(mesh, + equations::SerreGreenNaghdiEquations1D, + solver, + initial_condition, + ::BoundaryConditionPeriodic, + RealT, uEltype) + D = solver.D1 + + # create temporary storage + h = ones(RealT, nnodes(mesh)) + h_x = zero(h) + b = zero(h) + b_x = zero(h) + v_x = zero(h) + h_hpb_x = zero(h) + hv_x = zero(h) + v2_x = zero(h) + h2_v_vx_x = zero(h) + h_vx_x = zero(h) + p_h = zero(h) + p_x = zero(h) + tmp = zero(h) + M_h_p_h_bx2 = zero(h) + M_h3_3 = zero(h) + M_h2_bx = zero(h) + + if D isa PeriodicUpwindOperators + v_x_upwind = zero(h) + p_0 = zero(h) + + Dmat_minus = sparse(D.minus) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + if equations.bathymetry isa BathymetryMildSlope + factor = 0.75 + elseif equations.bathymetry isa BathymetryVariable + factor = 1.0 + else + throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry)")) + end + @. M_h_p_h_bx2 = h + factor * h * b_x^2 + scale_by_mass_matrix!(M_h_p_h_bx2, D) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D) + @. M_h2_bx = 0.5 * h^2 * b_x + scale_by_mass_matrix!(M_h2_bx, D) + system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus + - Diagonal(M_h2_bx)) + - Diagonal(M_h2_bx) * Dmat_minus) + + factorization = cholesky(system_matrix) + + cache = (; h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D, Dmat_minus, factorization) + else + if D isa FourierDerivativeOperator + Dmat = Matrix(D) + + cache = (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D, Dmat) + else + Dmat = sparse(D) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + if equations.bathymetry isa BathymetryMildSlope + factor = 0.75 + elseif equations.bathymetry isa BathymetryVariable + factor = 1.0 + else + throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry)")) + end + @. M_h_p_h_bx2 = h + factor * h * b_x^2 + scale_by_mass_matrix!(M_h_p_h_bx2, D) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D) + @. M_h2_bx = 0.5 * h^2 * b_x + scale_by_mass_matrix!(M_h2_bx, D) + system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + + Dmat' * (Diagonal(M_h3_3) * Dmat + - Diagonal(M_h2_bx)) + - Diagonal(M_h2_bx) * Dmat) + + factorization = cholesky(system_matrix) + + cache = (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D, Dmat, factorization) + end + end + + if equations.bathymetry isa BathymetryVariable + ψ = zero(h) + cache = (; cache..., ψ) + end + + return cache +end + +# Discretization that conserves +# - the total water mass (integral of h) as a linear invariant +# - the total momentum (integral of h v) as a nonlinear invariant +# - the total energy +# for periodic boundary conditions, see +# - Hendrik Ranocha and Mario Ricchiuto (2024) +# Structure-preserving approximations of the Serre-Green-Naghdi +# equations in standard and hyperbolic form +# [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) +# TODO: Implement source terms +function rhs!(dq, q, t, mesh, + equations::SerreGreenNaghdiEquations1D, + initial_condition, + ::BoundaryConditionPeriodic, + source_terms::Nothing, + solver, cache) + if cache.D isa PeriodicUpwindOperators + rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry) + else + rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry) + end + + return nothing +end + +function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat) + # Unpack physical parameters and SBP operator `D` as well as the + # SBP operator in sparse matrix form `Dmat` + g = gravity_constant(equations) + (; D, Dmat) = cache + + # `q` and `dq` are `ArrayPartition`s. They collect the individual + # arrays for the water height `h` and the velocity `v`. + h, v = q.x + dh, dv = dq.x + + @trixi_timeit timer() "hyperbolic terms" begin + # Compute all derivatives required below + (; h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, + h_vx_x, p_x, tmp, M_h, M_h3_3) = cache + + mul!(h_x, D, h) + mul!(v_x, D, v) + @. tmp = h^2 + mul!(h2_x, D, tmp) + @. tmp = h * v + mul!(hv_x, D, tmp) + @. tmp = v^2 + mul!(v2_x, D, tmp) + + @. tmp = h^2 * v * v_x + mul!(h2_v_vx_x, D, tmp) + @. tmp = h * v_x + mul!(h_vx_x, D, tmp) + inv6 = 1 / 6 + @. tmp = (0.5 * h^2 * (h * v_x + h_x * v) * v_x + - + inv6 * h * h2_v_vx_x + - + inv6 * h^2 * v * h_vx_x) + mul!(p_x, D, tmp) + + # Plain: h_t + (h v)_x = 0 + # + # Split form for energy conservation: + # h_t + h_x v + h v_x = 0 + @. dh = -(h_x * v + h * v_x) + + # Plain: h v_t + ... = 0 + # + # Split form for energy conservation: + @. tmp = -(g * h2_x - g * h * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + + p_x) + end + + @trixi_timeit timer() "assembling elliptic operator" begin + # The code below is equivalent to + # dv .= (Diagonal(h) - Dmat * Diagonal(1/3 .* h.^3) * Dmat) \ tmp + # but faster since the symbolic factorization is reused. + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off + # errors. We wrap it here to avoid issues with the factorization. + @. M_h = h + scale_by_mass_matrix!(M_h, D) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D) + system_matrix = Symmetric(Diagonal(M_h) + + + Dmat' * Diagonal(M_h3_3) * Dmat) + end + + @trixi_timeit timer() "solving elliptic system" begin + if issparse(system_matrix) + (; factorization) = cache + cholesky!(factorization, system_matrix; check = false) + if issuccess(factorization) + scale_by_mass_matrix!(tmp, D) + dv .= factorization \ tmp + else + # The factorization may fail if the time step is too large + # and h becomes negative. + fill!(dv, NaN) + end + else + factorization = cholesky!(system_matrix) + scale_by_mass_matrix!(tmp, D) + ldiv!(dv, factorization, tmp) + end + end + + return nothing +end + +function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat) + # Unpack physical parameters and SBP operator `D` as well as the + # SBP upwind operator in sparse matrix form `Dmat_minus` + g = gravity_constant(equations) + (; Dmat_minus) = cache + D_upwind = cache.D + D = D_upwind.central + + # `q` and `dq` are `ArrayPartition`s. They collect the individual + # arrays for the water height `h` and the velocity `v`. + h, v = q.x + dh, dv, dD = dq.x + fill!(dD, zero(eltype(dD))) + + @trixi_timeit timer() "hyperbolic terms" begin + # Compute all derivatives required below + (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_x, tmp, + M_h, M_h3_3) = cache + + mul!(h_x, D, h) + mul!(v_x, D, v) + mul!(v_x_upwind, D_upwind.minus, v) + @. tmp = h^2 + mul!(h2_x, D, tmp) + @. tmp = h * v + mul!(hv_x, D, tmp) + @. tmp = v^2 + mul!(v2_x, D, tmp) + + @. tmp = h^2 * v * v_x + mul!(h2_v_vx_x, D, tmp) + @. tmp = h * v_x + mul!(h_vx_x, D, tmp) + # p_+ + @. tmp = 0.5 * h^2 * (h * v_x + h_x * v) * v_x_upwind + mul!(p_x, D_upwind.plus, tmp) + # p_0 + minv6 = -1 / 6 + @. tmp = minv6 * (h * h2_v_vx_x + + + h^2 * v * h_vx_x) + mul!(p_x, D, tmp, 1.0, 1.0) + + # Plain: h_t + (h v)_x = 0 + # + # Split form for energy conservation: + # h_t + h_x v + h v_x = 0 + @. dh = -(h_x * v + h * v_x) + + # Plain: h v_t + ... = 0 + # + # Split form for energy conservation: + @. tmp = -(g * h2_x - g * h * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + + p_x) + end + + @trixi_timeit timer() "assembling elliptic operator" begin + # The code below is equivalent to + # dv .= (Diagonal(h) - Dmat_plus * Diagonal(1/3 .* h.^3) * Dmat_minus) \ tmp + # but faster since the symbolic factorization is reused. + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off errors. + # We wrap it here to avoid issues with the factorization. + @. M_h = h + scale_by_mass_matrix!(M_h, D) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D) + system_matrix = Symmetric(Diagonal(M_h) + + + Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) + end + + @trixi_timeit timer() "solving elliptic system" begin + (; factorization) = cache + cholesky!(factorization, system_matrix; check = false) + if issuccess(factorization) + scale_by_mass_matrix!(tmp, D) + dv .= factorization \ tmp + else + # The factorization may fail if the time step is too large + # and h becomes negative. + fill!(dv, NaN) + end + end + + return nothing +end + +function rhs_sgn_central!(dq, q, equations, source_terms, cache, + ::Union{BathymetryMildSlope, BathymetryVariable}) + # Unpack physical parameters and SBP operator `D` as well as the + # SBP operator in sparse matrix form `Dmat` + g = gravity_constant(equations) + (; D, Dmat) = cache + + # `q` and `dq` are `ArrayPartition`s. They collect the individual + # arrays for the water height `h` and the velocity `v`. + h, v = q.x + dh, dv, dD = dq.x + fill!(dD, zero(eltype(dD))) + + @trixi_timeit timer() "hyperbolic terms" begin + # Compute all derivatives required below + (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + if equations.bathymetry isa BathymetryVariable + (; ψ) = cache + end + + @. b = equations.eta0 - q.x[3] + mul!(b_x, D, b) + + mul!(h_x, D, h) + mul!(v_x, D, v) + @. tmp = h * (h + b) + mul!(h_hpb_x, D, tmp) + @. tmp = h * v + mul!(hv_x, D, tmp) + @. tmp = v^2 + mul!(v2_x, D, tmp) + + @. tmp = h^2 * v * v_x + mul!(h2_v_vx_x, D, tmp) + @. tmp = h * v_x + mul!(h_vx_x, D, tmp) + inv6 = 1 / 6 + @. p_h = ( 0.5 * h * (h * v_x + h_x * v) * v_x + - inv6 * h2_v_vx_x + - inv6 * h * v * h_vx_x) + @. tmp = h * b_x * v^2 + mul!(p_x, D, tmp) + @. p_h += 0.25 * p_x + if equations.bathymetry isa BathymetryVariable + @. ψ = 0.125 * p_x + end + @. tmp = b_x * v + mul!(p_x, D, tmp) + @. p_h += 0.25 * h * v * p_x + if equations.bathymetry isa BathymetryVariable + @. ψ += 0.125 * h * v * p_x + end + @. p_h = p_h - 0.25 * (h_x * v + h * v_x) * b_x * v + if equations.bathymetry isa BathymetryVariable + @. ψ = ψ - 0.125 * (h_x * v + h * v_x) * b_x * v + end + @. tmp = p_h * h + mul!(p_x, D, tmp) + + # Plain: h_t + (h v)_x = 0 + # + # Split form for energy conservation: + # h_t + h_x v + h v_x = 0 + @. dh = -(h_x * v + h * v_x) + + # Plain: h v_t + ... = 0 + # + # Split form for energy conservation: + @. tmp = -( g * h_hpb_x - g * (h + b) * h_x + + 0.5 * h * v2_x + - 0.5 * v^2 * h_x + + 0.5 * hv_x * v + - 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) + if equations.bathymetry isa BathymetryVariable + @. tmp = tmp - ψ * b_x + end + end + + @trixi_timeit timer() "assembling elliptic operator" begin + if equations.bathymetry isa BathymetryMildSlope + factor = 0.75 + else # equations.bathymetry isa BathymetryVariable + factor = 1.0 + end + # The code below is equivalent to + # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp + # but faster since the symbolic factorization is reused. + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off + # errors. We wrap it here to avoid issues with the factorization. + @. M_h_p_h_bx2 = h + factor * h * b_x^2 + scale_by_mass_matrix!(M_h_p_h_bx2, D) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D) + @. M_h2_bx = 0.5 * h^2 * b_x + scale_by_mass_matrix!(M_h2_bx, D) + system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + + Dmat' * (Diagonal(M_h3_3) * Dmat + - Diagonal(M_h2_bx)) + - Diagonal(M_h2_bx) * Dmat) + end + + @trixi_timeit timer() "solving elliptic system" begin + if issparse(system_matrix) + (; factorization) = cache + cholesky!(factorization, system_matrix; check = false) + if issuccess(factorization) + scale_by_mass_matrix!(tmp, D) + dv .= factorization \ tmp + else + # The factorization may fail if the time step is too large + # and h becomes negative. + fill!(dv, NaN) + end + else + factorization = cholesky!(system_matrix) + scale_by_mass_matrix!(tmp, D) + ldiv!(dv, factorization, tmp) + end + end + + return nothing +end + +function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, + ::Union{BathymetryMildSlope, BathymetryVariable}) + # Unpack physical parameters and SBP operator `D` as well as the + # SBP operator in sparse matrix form `Dmat` + g = gravity_constant(equations) + (; Dmat_minus) = cache + D_upwind = cache.D + D = D_upwind.central + + # `q` and `dq` are `ArrayPartition`s. They collect the individual + # arrays for the water height `h` and the velocity `v`. + h, v, b = q.x + dh, dv, db = dq.x + fill!(db, zero(eltype(db))) + + @trixi_timeit timer() "hyperbolic terms" begin + # Compute all derivatives required below + (; h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + if equations.bathymetry isa BathymetryVariable + (; ψ) = cache + end + + @. b = equations.eta0 - q.x[3] + mul!(b_x, D, b) + + mul!(h_x, D, h) + mul!(v_x, D, v) + mul!(v_x_upwind, D_upwind.minus, v) + @. tmp = h * (h + b) + mul!(h_hpb_x, D, tmp) + @. tmp = h * v + mul!(hv_x, D, tmp) + @. tmp = v^2 + mul!(v2_x, D, tmp) + + @. tmp = h^2 * v * v_x + mul!(h2_v_vx_x, D, tmp) + @. tmp = h * v_x + mul!(h_vx_x, D, tmp) + # p_0 + minv6 = -1 / 6 + @. p_h = minv6 * ( h2_v_vx_x + + h * v * h_vx_x) + @. tmp = h * b_x * v^2 + mul!(p_x, D, tmp) + @. p_h += 0.25 * p_x + if equations.bathymetry isa BathymetryVariable + @. ψ = 0.125 * p_x + end + @. tmp = b_x * v + mul!(p_x, D, tmp) + @. p_h += 0.25 * h * v * p_x + if equations.bathymetry isa BathymetryVariable + @. ψ += 0.125 * h * v * p_x + end + @. p_0 = p_h * h + mul!(p_x, D, p_0) + # p_+ + @. tmp = ( 0.5 * h * (h * v_x + h_x * v) * v_x_upwind + - 0.25 * (h_x * v + h * v_x) * b_x * v) + if equations.bathymetry isa BathymetryVariable + @. ψ = 0.125 * (h_x * v + h * v_x) * b_x * v + end + @. p_h = p_h + tmp + @. tmp = tmp * h + mul!(p_x, D_upwind.plus, tmp, 1.0, 1.0) + + # Plain: h_t + (h v)_x = 0 + # + # Split form for energy conservation: + # h_t + h_x v + h v_x = 0 + @. dh = -(h_x * v + h * v_x) + + # Plain: h v_t + ... = 0 + # + # Split form for energy conservation: + @. tmp = -( g * h_hpb_x - g * (h + b) * h_x + + 0.5 * h * v2_x + - 0.5 * v^2 * h_x + + 0.5 * hv_x * v + - 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) + if equations.bathymetry isa BathymetryVariable + @. tmp = tmp - ψ * b_x + end + end + + @trixi_timeit timer() "assembling elliptic operator" begin + if equations.bathymetry isa BathymetryMildSlope + factor = 0.75 + else # equations.bathymetry isa BathymetryVariable + factor = 1.0 + end + # The code below is equivalent to + # dv .= (Diagonal(h .+ 0.75 .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp + # but faster since the symbolic factorization is reused. + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off + # errors. We wrap it here to avoid issues with the factorization. + @. M_h_p_h_bx2 = h + factor * h * b_x^2 + scale_by_mass_matrix!(M_h_p_h_bx2, D) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D) + @. M_h2_bx = 0.5 * h^2 * b_x + scale_by_mass_matrix!(M_h2_bx, D) + system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus + - Diagonal(M_h2_bx)) + - Diagonal(M_h2_bx) * Dmat_minus) + end + + @trixi_timeit timer() "solving elliptic system" begin + if issparse(system_matrix) + (; factorization) = cache + cholesky!(factorization, system_matrix; check = false) + if issuccess(factorization) + scale_by_mass_matrix!(tmp, D) + dv .= factorization \ tmp + else + # The factorization may fail if the time step is too large + # and h becomes negative. + fill!(dv, NaN) + end + else + factorization = cholesky!(system_matrix) + scale_by_mass_matrix!(tmp, D) + ldiv!(dv, factorization, tmp) + end + end + + return nothing +end + +@inline function prim2cons(q, equations::SerreGreenNaghdiEquations1D) + h = waterheight(q, equations) + v = velocity(q, equations) + b = bathymetry(q, equations) + + hv = h * v + return SVector(h, hv, b) +end + +@inline function cons2prim(u, equations::SerreGreenNaghdiEquations1D) + h, hv, b = u + + eta = h + b + v = hv / h + D = equations.eta0 - b + return SVector(eta, v, D) +end + +@inline function waterheight_total(q, equations::SerreGreenNaghdiEquations1D) + return q[1] +end + +@inline function velocity(q, equations::SerreGreenNaghdiEquations1D) + return q[2] +end + +@inline function bathymetry(q, equations::SerreGreenNaghdiEquations1D) + D = q[3] + return equations.eta0 - D +end + +@inline function waterheight(q, equations::SerreGreenNaghdiEquations1D) + return waterheight_total(q, equations) - bathymetry(q, equations) +end + +# The entropy/energy takes the whole `q` for every point in space +""" + energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache) + +Return the modified total energy of the primitive variables `q_global` for the +[`SerreGreenNaghdiEquations1D`](@ref). +It contains an additional term containing a +derivative compared to the usual [`energy_total`](@ref) modeling +non-hydrostatic contributions. The [`energy_total_modified`](@ref) +is a conserved quantity (for periodic boundary conditions). + +For a [`bathymetry_flat`](@ref) the total energy is given by +```math +\\frac{1}{2} g h^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2. +``` + +`q_global` is a vector of the primitive variables at ALL nodes. +`cache` needs to hold the SBP operators used by the `solver`. +""" +function energy_total_modified(q_global, + equations::SerreGreenNaghdiEquations1D, + cache) + # unpack physical parameters and SBP operator `D` + g = equations.gravity + (; D, v_x) = cache + + # `q_global` is an `ArrayPartition`. It collects the individual arrays for + # the water height `h` and the velocity `v`. + h, v = q_global.x + + N = length(v) + e = zeros(eltype(q_global), N) + + # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 v_x^2 + if D isa PeriodicUpwindOperators + mul!(v_x, D.minus, v) + else + mul!(v_x, D, v) + end + + @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h^3 * v_x^2 + + return e +end diff --git a/src/equations/serre_green_naghdi_flat_1d.jl b/src/equations/serre_green_naghdi_flat_1d.jl deleted file mode 100644 index 0858679f..00000000 --- a/src/equations/serre_green_naghdi_flat_1d.jl +++ /dev/null @@ -1,456 +0,0 @@ -@doc raw""" - SerreGreenNaghdiEquations1D(; gravity_constant, eta0 = 0.0) - -Serre-Green-Naghdi system in one spatial dimension. -The equations for flat bathymetry are given by -```math -\begin{aligned} - h_t + (h v)_x &= 0,\\ - h v_t - \frac{1}{3} (h^3 v_{tx})_x + \frac{1}{2} g (h^2)_x + \frac{1}{2} h (v^2)_x + p_x &= 0,\\ - p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx}. -\end{aligned} -``` -The unknown quantities of the Serre-Green-Naghdi equations are the -total water height ``\eta = h + b`` and the velocity ``v``. -The gravitational constant is denoted by `g` and the bottom topography -(bathymetry) ``b = \eta_0 - D``. The water height above the bathymetry -is therefore given by ``h = \eta - \eta_0 + D``. The Serre-Green-Naghdi -equations are only implemented for ``\eta_0 = 0``. -The total water height is therefore given by ``\eta = h + b``. - -References for the Serre-Green-Naghdi system can be found in -- Serre (1953) - Contribution â l'étude des écoulements permanents et variables dans les canaux - [DOI: 10.1051/lhb/1953034](https://doi.org/10.1051/lhb/1953034) -- Green and Naghdi (1976) - A derivation of equations for wave propagation in water of variable depth - [DOI: 10.1017/S0022112076002425](https://doi.org/10.1017/S0022112076002425) - -The semidiscretization implemented here conserves -- the total water mass (integral of h) as a linear invariant -- the total momentum (integral of h v) as a nonlinear invariant -- the total energy - -for periodic boundary conditions, see -- Hendrik Ranocha and Mario Ricchiuto (2024) - Structure-preserving approximations of the Serre-Green-Naghdi - equations in standard and hyperbolic form - [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) -""" -struct SerreGreenNaghdiEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <: - AbstractSerreGreenNaghdiEquations{1, 3} - bathymetry::Bathymetry # type of bathymetry - gravity::RealT # gravitational constant - eta0::RealT # constant still-water surface -end - -function SerreGreenNaghdiEquations1D(; gravity_constant, eta0 = 0.0) - eta0 == 0.0 || - @warn "The still-water surface needs to be 0 for the Serre-Green-Naghdi equations" - SerreGreenNaghdiEquations1D(bathymetry_flat, gravity_constant, eta0) -end - -varnames(::typeof(prim2prim), ::SerreGreenNaghdiEquations1D) = ("η", "v", "D") -varnames(::typeof(prim2cons), ::SerreGreenNaghdiEquations1D) = ("h", "hv", "b") - -""" - initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEquations1D, mesh) - -A soliton solution used for convergence tests in a periodic domain. -""" -function initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEquations1D, - mesh) - g = equations.gravity - - # setup parameters data - h1 = 1.0 - h2 = 1.2 - c = sqrt(g * h2) - - x_t = mod(x - c * t - xmin(mesh), xmax(mesh) - xmin(mesh)) + xmin(mesh) - - h = h1 + (h2 - h1) * sech(x_t / 2 * sqrt(3 * (h2 - h1) / (h1^2 * h2)))^2 - v = c * (1 - h1 / h) - - return SVector(h, v, 0) -end - -function create_cache(mesh, - equations::SerreGreenNaghdiEquations1D{BathymetryFlat}, - solver, - initial_condition, - ::BoundaryConditionPeriodic, - RealT, uEltype) - D = solver.D1 - - # create temporary storage - h = ones(RealT, nnodes(mesh)) - h_x = zero(h) - v_x = zero(h) - h2_x = zero(h) - hv_x = zero(h) - v2_x = zero(h) - h2_v_vx_x = zero(h) - h_vx_x = zero(h) - p_x = zero(h) - tmp = zero(h) - M_h = zero(h) - M_h3_3 = zero(h) - - if D isa PeriodicUpwindOperators - v_x_upwind = zero(h) - - Dmat_minus = sparse(D.minus) - - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) - factorization = cholesky(system_matrix) - - cache = (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_x, tmp, - M_h, M_h3_3, - D, Dmat_minus, factorization) - else - if D isa FourierDerivativeOperator - Dmat = Matrix(D) - - cache = (; h_x, v_x, h2_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_x, tmp, - M_h, M_h3_3, - D, Dmat) - else - Dmat = sparse(D) - - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat' * Diagonal(M_h3_3) * Dmat) - factorization = cholesky(system_matrix) - - cache = (; h_x, v_x, h2_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_x, tmp, - M_h, M_h3_3, - D, Dmat, factorization) - end - end - - return cache -end - -# Discretization that conserves -# - the total water mass (integral of h) as a linear invariant -# - the total momentum (integral of h v) as a nonlinear invariant -# - the total energy -# for periodic boundary conditions, see -# - Hendrik Ranocha and Mario Ricchiuto (2024) -# Structure-preserving approximations of the Serre-Green-Naghdi -# equations in standard and hyperbolic form -# [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) -# TODO: Implement source terms -# TODO: Implement variable bathymetry -function rhs!(dq, q, t, mesh, - equations::SerreGreenNaghdiEquations1D, - initial_condition, - ::BoundaryConditionPeriodic, - source_terms::Nothing, - solver, cache) - if cache.D isa PeriodicUpwindOperators - rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry) - else - rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry) - end - - return nothing -end - -function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP operator in sparse matrix form `Dmat` - g = gravity_constant(equations) - (; D, Dmat) = cache - - # `q` and `dq` are `ArrayPartition`s. They collect the individual - # arrays for the water height `h` and the velocity `v`. - h, v = q.x - dh, dv = dq.x - - @trixi_timeit timer() "hyperbolic terms" begin - # Compute all derivatives required below - (; h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, - h_vx_x, p_x, tmp, M_h, M_h3_3) = cache - - mul!(h_x, D, h) - mul!(v_x, D, v) - @. tmp = h^2 - mul!(h2_x, D, tmp) - @. tmp = h * v - mul!(hv_x, D, tmp) - @. tmp = v^2 - mul!(v2_x, D, tmp) - - @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) - @. tmp = h * v_x - mul!(h_vx_x, D, tmp) - inv6 = 1 / 6 - @. tmp = (0.5 * h^2 * (h * v_x + h_x * v) * v_x - - - inv6 * h * h2_v_vx_x - - - inv6 * h^2 * v * h_vx_x) - mul!(p_x, D, tmp) - - # Plain: h_t + (h v)_x = 0 - # - # Split form for energy conservation: - # h_t + h_x v + h v_x = 0 - @. dh = -(h_x * v + h * v_x) - - # Plain: h v_t + ... = 0 - # - # Split form for energy conservation: - @. tmp = -(g * h2_x - g * h * h_x - + - 0.5 * h * v2_x - - - 0.5 * v^2 * h_x - + - 0.5 * hv_x * v - - - 0.5 * h * v * v_x - + - p_x) - end - - @trixi_timeit timer() "assembling elliptic operator" begin - # The code below is equivalent to - # dv .= (Diagonal(h) - Dmat * Diagonal(1/3 .* h.^3) * Dmat) \ tmp - # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat' * Diagonal(M_h3_3) * Dmat) - end - - @trixi_timeit timer() "solving elliptic system" begin - if issparse(system_matrix) - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) - end - else - factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D) - ldiv!(dv, factorization, tmp) - end - end - - return nothing -end - -function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP upwind operator in sparse matrix form `Dmat_minus` - g = gravity_constant(equations) - (; Dmat_minus) = cache - D_upwind = cache.D - D = D_upwind.central - - # `q` and `dq` are `ArrayPartition`s. They collect the individual - # arrays for the water height `h` and the velocity `v`. - h, v = q.x - dh, dv = dq.x - - @trixi_timeit timer() "hyperbolic terms" begin - # Compute all derivatives required below - (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_x, tmp, - M_h, M_h3_3) = cache - - mul!(h_x, D, h) - mul!(v_x, D, v) - mul!(v_x_upwind, D_upwind.minus, v) - @. tmp = h^2 - mul!(h2_x, D, tmp) - @. tmp = h * v - mul!(hv_x, D, tmp) - @. tmp = v^2 - mul!(v2_x, D, tmp) - - @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) - @. tmp = h * v_x - mul!(h_vx_x, D, tmp) - # p_+ - @. tmp = 0.5 * h^2 * (h * v_x + h_x * v) * v_x_upwind - mul!(p_x, D_upwind.plus, tmp) - # p_0 - minv6 = -1 / 6 - @. tmp = minv6 * (h * h2_v_vx_x - + - h^2 * v * h_vx_x) - mul!(p_x, D, tmp, 1.0, 1.0) - - # Plain: h_t + (h v)_x = 0 - # - # Split form for energy conservation: - # h_t + h_x v + h v_x = 0 - @. dh = -(h_x * v + h * v_x) - - # Plain: h v_t + ... = 0 - # - # Split form for energy conservation: - @. tmp = -(g * h2_x - g * h * h_x - + - 0.5 * h * v2_x - - - 0.5 * v^2 * h_x - + - 0.5 * hv_x * v - - - 0.5 * h * v * v_x - + - p_x) - end - - @trixi_timeit timer() "assembling elliptic operator" begin - # The code below is equivalent to - # dv .= (Diagonal(h) - Dmat_plus * Diagonal(1/3 .* h.^3) * Dmat_minus) \ tmp - # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off errors. - # We wrap it here to avoid issues with the factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) - end - - @trixi_timeit timer() "solving elliptic system" begin - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) - end - end - - return nothing -end - -@inline function prim2cons(q, equations::SerreGreenNaghdiEquations1D) - h = waterheight(q, equations) - v = velocity(q, equations) - b = bathymetry(q, equations) - - hv = h * v - return SVector(h, hv, b) -end - -@inline function cons2prim(u, equations::SerreGreenNaghdiEquations1D) - h, hv, b = u - - eta = h + b - v = hv / h - D = equations.eta0 - b - return SVector(eta, v, D) -end - -@inline function waterheight_total(q, equations::SerreGreenNaghdiEquations1D) - return q[1] -end - -@inline function velocity(q, equations::SerreGreenNaghdiEquations1D) - return q[2] -end - -@inline function bathymetry(q, equations::SerreGreenNaghdiEquations1D) - D = q[3] - return equations.eta0 - D -end - -@inline function waterheight(q, equations::SerreGreenNaghdiEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - -# The entropy/energy takes the whole `q` for every point in space -""" - energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache) - -Return the modified total energy of the primitive variables `q_global` for the -[`SerreGreenNaghdiEquations1D`](@ref). -It contains an additional term containing a -derivative compared to the usual [`energy_total`](@ref) modeling -non-hydrostatic contributions. The [`energy_total_modified`](@ref) -is a conserved quantity (for periodic boundary conditions). - -For a [`bathymetry_flat`](@ref) the total energy is given by -```math -\\frac{1}{2} g h^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2. -``` - -`q_global` is a vector of the primitive variables at ALL nodes. -`cache` needs to hold the SBP operators used by the `solver`. -""" -function energy_total_modified(q_global, - equations::SerreGreenNaghdiEquations1D, - cache) - # unpack physical parameters and SBP operator `D` - g = equations.gravity - (; D, v_x) = cache - - # `q_global` is an `ArrayPartition`. It collects the individual arrays for - # the water height `h` and the velocity `v`. - h, v = q_global.x - - N = length(v) - e = zeros(eltype(q_global), N) - - # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 v_x^2 - if D isa PeriodicUpwindOperators - mul!(v_x, D.minus, v) - else - mul!(v_x, D, v) - end - - @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h^3 * v_x^2 - - return e -end diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index c5bb1e55..9de8d5f5 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -21,6 +21,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, 550_000) end + @trixi_testset "serre_green_naghdi_soliton.jl with bathymetry_mild_slope" begin + # same values as serre_green_naghdi_soliton.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton.jl"), + equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + tspan=(0.0, 0.1), + l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], + linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], + cons_error=[0.0, 8.174581012099225e-10, 0.0], + change_waterheight=0.0, + change_entropy_modified=-3.1093350116861984e-11) + + @test_allocations(semi, sol, 800_000) + end + + @trixi_testset "serre_green_naghdi_soliton.jl with bathymetry_variable" begin + # same values as serre_green_naghdi_soliton.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton.jl"), + equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + tspan=(0.0, 0.1), + l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], + linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], + cons_error=[0.0, 8.174581012099225e-10, 0.0], + change_waterheight=0.0, + change_entropy_modified=-3.1093350116861984e-11) + + @test_allocations(semi, sol, 800_000) + end + @trixi_testset "serre_green_naghdi_soliton_fourier.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), @@ -34,6 +64,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=450_000) end + @trixi_testset "serre_green_naghdi_soliton_fourier.jl with bathymetry_mild_slope" begin + # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_fourier.jl"), + equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + tspan=(0.0, 0.1), + l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], + linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.627409566637652e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=-3.097966327914037e-11) + + @test_allocations(semi, sol, allocs=850_000) + end + + @trixi_testset "serre_green_naghdi_soliton_fourier.jl with bathymetry_variable" begin + # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_fourier.jl"), + equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + tspan=(0.0, 0.1), + l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], + linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.627409566637652e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=-3.097966327914037e-11) + + @test_allocations(semi, sol, allocs=850_000) + end + @trixi_testset "serre_green_naghdi_soliton_upwind.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), @@ -47,6 +107,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=500_000) end + @trixi_testset "serre_green_naghdi_soliton_upwind.jl with bathymetry_mild_slope" begin + # same as serre_green_naghdi_soliton_upwind.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_upwind.jl"), + tspan=(0.0, 0.1), + equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], + linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], + cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], + change_waterheight=4.263256414560601e-14, + change_entropy_modified=-3.1036506698001176e-11) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_soliton_upwind.jl with bathymetry_variable" begin + # same as serre_green_naghdi_soliton_upwind.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_upwind.jl"), + tspan=(0.0, 0.1), + equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], + linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], + cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], + change_waterheight=4.263256414560601e-14, + change_entropy_modified=-3.1036506698001176e-11) + + @test_allocations(semi, sol, allocs=750_000) + end + @trixi_testset "serre_green_naghdi_soliton_relaxation.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), @@ -59,6 +149,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=450_000) end + + @trixi_testset "serre_green_naghdi_soliton_relaxation.jl with bathymetry_mild_slope" begin + # same as serre_green_naghdi_soliton_relaxation.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_relaxation.jl"), + tspan=(0.0, 0.1), + equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], + linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=0.0) + + @test_allocations(semi, sol, allocs=850_000) + end + + @trixi_testset "serre_green_naghdi_soliton_relaxation.jl with bathymetry_variable" begin + # same as serre_green_naghdi_soliton_relaxation.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_relaxation.jl"), + tspan=(0.0, 0.1), + equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], + linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=0.0) + + @test_allocations(semi, sol, allocs=850_000) + end end end # module From 6e1b1e6f7c38952f812002c6a13895559f718d66 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 15 Aug 2024 17:40:30 +0200 Subject: [PATCH 02/34] Dingemans --- .gitignore | 1 + .../serre_green_naghdi_dingemans.jl | 49 ++++++++++ src/equations/serre_green_naghdi_1d.jl | 89 ++++++++++++++++--- test/test_serre_green_naghdi_1d.jl | 29 ++++++ 4 files changed, 157 insertions(+), 11 deletions(-) create mode 100644 examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl diff --git a/.gitignore b/.gitignore index c1dd5bb5..ffeb3567 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ out*/ /docs/build/ run run/* +**/*.code-workspace diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl new file mode 100644 index 00000000..854a8879 --- /dev/null +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl @@ -0,0 +1,49 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: upwind_operators, periodic_derivative_operator + +############################################################################### +# Semidiscretization of the Svärd-Kalisch equations + +bathymetry = bathymetry_variable # or bathymetry_mild_slope +equations = SerreGreenNaghdiEquations1D(bathymetry; + gravity_constant = 9.81) + +initial_condition = initial_condition_dingemans +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -140.0 +coordinates_max = 100.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic upwind SBP operators of accuracy order 4 +# (note that only central-type with even order are used) +D1 = upwind_operators(periodic_derivative_operator; + derivative_order = 1, accuracy_order = 3, + xmin = xmin(mesh), xmax = xmax(mesh), + N = nnodes(mesh)) +solver = Solver(D1, nothing) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 70.0) +ode = semidiscretize(semi, tspan) +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + momentum, + entropy, + entropy_modified)) + +callbacks = CallbackSet(analysis_callback, summary_callback) + +saveat = range(tspan..., length = 500) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 4deee588..d10836db 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -77,11 +77,19 @@ end function SerreGreenNaghdiEquations1D(bathymetry = bathymetry_flat; gravity_constant, eta0 = 0.0) - eta0 == 0.0 || - @warn "The still-water surface needs to be 0 for the Serre-Green-Naghdi equations" SerreGreenNaghdiEquations1D(bathymetry, gravity_constant, eta0) end +function get_name(equations::SerreGreenNaghdiEquations1D) + name = equations |> typeof |> nameof |> string + bathymetry = equations.bathymetry + if bathymetry isa BathymetryFlat + return name + else # variable bathymetry + return name * "-" * string(nameof(typeof(bathymetry))) + end +end + varnames(::typeof(prim2prim), ::SerreGreenNaghdiEquations1D) = ("η", "v", "D") varnames(::typeof(prim2cons), ::SerreGreenNaghdiEquations1D) = ("h", "hv", "b") @@ -107,6 +115,46 @@ function initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEqu return SVector(h, v, 0) end +""" + initial_condition_dingemans(x, t, equations::SerreGreenNaghdiEquations1D, mesh) + +The initial condition that uses the dispersion relation of the Euler equations +to approximate waves generated by a wave maker as it is done by experiments of +Dingemans. The topography is a trapezoidal. It is assumed that `equations.eta0 = 0.8`. + +References: +- Magnus Svärd, Henrik Kalisch (2023) + A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization + [arXiv: 2302.09924](https://arxiv.org/abs/2302.09924) +- Maarten W. Dingemans (1994) + Comparison of computations with Boussinesq-like models and laboratory measurements + [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download) +""" +function initial_condition_dingemans(x, t, equations::SerreGreenNaghdiEquations1D, mesh) + h0 = 0.8 + A = 0.02 + # omega = 2*pi/(2.02*sqrt(2)) + k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * h0), 1.0) using Roots.jl + if x < -34.5 * pi / k || x > -8.5 * pi / k + h = 0.0 + else + h = A * cos(k * x) + end + v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0 + if 11.01 <= x && x < 23.04 + b = 0.6 * (x - 11.01) / (23.04 - 11.01) + elseif 23.04 <= x && x < 27.04 + b = 0.6 + elseif 27.04 <= x && x < 33.07 + b = 0.6 * (33.07 - x) / (33.07 - 27.04) + else + b = 0.0 + end + eta = h + h0 + D = equations.eta0 - b + return SVector(eta, v, D) +end + # flat bathymetry function create_cache(mesh, equations::SerreGreenNaghdiEquations1D{BathymetryFlat}, @@ -245,7 +293,7 @@ function create_cache(mesh, factorization = cholesky(system_matrix) - cache = (; h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, + cache = (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx, D, Dmat_minus, factorization) @@ -253,7 +301,7 @@ function create_cache(mesh, if D isa FourierDerivativeOperator Dmat = Matrix(D) - cache = (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx, D, Dmat) @@ -285,7 +333,7 @@ function create_cache(mesh, factorization = cholesky(system_matrix) - cache = (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx, D, Dmat, factorization) @@ -531,13 +579,13 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - h, v = q.x + eta, v = q.x dh, dv, dD = dq.x fill!(dD, zero(eltype(dD))) @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below - (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache if equations.bathymetry isa BathymetryVariable @@ -545,6 +593,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end @. b = equations.eta0 - q.x[3] + @. h = eta - b mul!(b_x, D, b) mul!(h_x, D, h) @@ -662,13 +711,13 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - h, v, b = q.x + eta, v, b = q.x dh, dv, db = dq.x fill!(db, zero(eltype(db))) @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below - (; h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, + (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache if equations.bathymetry isa BathymetryVariable @@ -676,6 +725,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end @. b = equations.eta0 - q.x[3] + @. h = eta - b mul!(b_x, D, b) mul!(h_x, D, h) @@ -856,14 +906,31 @@ function energy_total_modified(q_global, N = length(v) e = zeros(eltype(q_global), N) - # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 v_x^2 + # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 w^2 + # and + 1/8 h (v b_x)^2 for full bathymetry without mild-slope approximation if D isa PeriodicUpwindOperators mul!(v_x, D.minus, v) else mul!(v_x, D, v) end - @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h^3 * v_x^2 + if equations.bathymetry isa BathymetryFlat + b_x = cache.tmp + fill!(b_x, zero(eltype(b_x))) + else + (; b, b_x) = cache + @. b = equations.eta0 - q_global.x[3] + if D isa PeriodicUpwindOperators + mul!(b_x, D.central, b) + else + mul!(b_x, D, b) + end + end + + @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2 + if equations.bathymetry isa BathymetryVariable + @. e += 1 / 8 * h * (v * b_x)^2 + end return e end diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index 9de8d5f5..367fc0bf 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -179,6 +179,35 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=850_000) end + + @trixi_testset "serre_green_naghdi_dingemans.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_dingemans.jl"), + tspan=(0.0, 1.0), + l2=[0.2463295492331802, 0.805357513755897, 0.0], + linf=[0.036265772921736605, 0.11763218152350845, 0.0], + cons_error=[5.684341886080802e-14, 3.509473432346808e-5, 0.0], + change_waterheight=5.684341886080802e-14, + change_entropy=1.9489684518703143e-5, + change_entropy_modified=-1.2177679309388623e-8) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_dingemans.jl with bathymetry_mild_slope" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_dingemans.jl"), + tspan=(0.0, 1.0), + equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + l2=[0.2463295492331802, 0.805357513755897, 0.0], + linf=[0.036265772921736605, 0.11763218152350845, 0.0], + cons_error=[5.684341886080802e-14, 3.509473432346624e-5, 0.0], + change_waterheight=5.684341886080802e-14, + change_entropy=1.9489684518703143e-5, + change_entropy_modified=-1.2177679309388623e-8) + + @test_allocations(semi, sol, allocs=750_000) + end end end # module From 37d42016c3f706785a5caa871ebc9d39bf202991 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 15 Aug 2024 17:51:37 +0200 Subject: [PATCH 03/34] format --- src/equations/serre_green_naghdi_1d.jl | 116 +++++++++++++++---------- test/test_serre_green_naghdi_1d.jl | 27 ++++-- 2 files changed, 88 insertions(+), 55 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index d10836db..7ecd3d1d 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -287,24 +287,27 @@ function create_cache(mesh, @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus - - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat_minus) + + + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus + - + Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx) * Dmat_minus) factorization = cholesky(system_matrix) cache = (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, - M_h_p_h_bx2, M_h3_3, M_h2_bx, - D, Dmat_minus, factorization) + h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D, Dmat_minus, factorization) else if D isa FourierDerivativeOperator Dmat = Matrix(D) cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_h, p_x, tmp, - M_h_p_h_bx2, M_h3_3, M_h2_bx, - D, Dmat) + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D, Dmat) else Dmat = sparse(D) @@ -327,16 +330,19 @@ function create_cache(mesh, @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + Dmat' * (Diagonal(M_h3_3) * Dmat - - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat) + + + Dmat' * (Diagonal(M_h3_3) * Dmat + - + Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx) * Dmat) factorization = cholesky(system_matrix) cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_h, p_x, tmp, - M_h_p_h_bx2, M_h3_3, M_h2_bx, - D, Dmat, factorization) + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D, Dmat, factorization) end end @@ -586,8 +592,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_h, p_x, tmp, - M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache if equations.bathymetry isa BathymetryVariable (; ψ) = cache end @@ -610,9 +616,11 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, @. tmp = h * v_x mul!(h_vx_x, D, tmp) inv6 = 1 / 6 - @. p_h = ( 0.5 * h * (h * v_x + h_x * v) * v_x - - inv6 * h2_v_vx_x - - inv6 * h * v * h_vx_x) + @. p_h = (0.5 * h * (h * v_x + h_x * v) * v_x + - + inv6 * h2_v_vx_x + - + inv6 * h * v * h_vx_x) @. tmp = h * b_x * v^2 mul!(p_x, D, tmp) @. p_h += 0.25 * p_x @@ -641,13 +649,17 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -( g * h_hpb_x - g * (h + b) * h_x - + 0.5 * h * v2_x - - 0.5 * v^2 * h_x - + 0.5 * hv_x * v - - 0.5 * h * v * v_x - + p_x - + 1.5 * p_h * b_x) + @. tmp = -(g * h_hpb_x - g * (h + b) * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) if equations.bathymetry isa BathymetryVariable @. tmp = tmp - ψ * b_x end @@ -673,9 +685,12 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + Dmat' * (Diagonal(M_h3_3) * Dmat - - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat) + + + Dmat' * (Diagonal(M_h3_3) * Dmat + - + Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx) * Dmat) end @trixi_timeit timer() "solving elliptic system" begin @@ -718,8 +733,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, - h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, - M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache if equations.bathymetry isa BathymetryVariable (; ψ) = cache end @@ -744,8 +759,9 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, mul!(h_vx_x, D, tmp) # p_0 minv6 = -1 / 6 - @. p_h = minv6 * ( h2_v_vx_x - + h * v * h_vx_x) + @. p_h = minv6 * (h2_v_vx_x + + + h * v * h_vx_x) @. tmp = h * b_x * v^2 mul!(p_x, D, tmp) @. p_h += 0.25 * p_x @@ -761,8 +777,9 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, @. p_0 = p_h * h mul!(p_x, D, p_0) # p_+ - @. tmp = ( 0.5 * h * (h * v_x + h_x * v) * v_x_upwind - - 0.25 * (h_x * v + h * v_x) * b_x * v) + @. tmp = (0.5 * h * (h * v_x + h_x * v) * v_x_upwind + - + 0.25 * (h_x * v + h * v_x) * b_x * v) if equations.bathymetry isa BathymetryVariable @. ψ = 0.125 * (h_x * v + h * v_x) * b_x * v end @@ -779,13 +796,17 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -( g * h_hpb_x - g * (h + b) * h_x - + 0.5 * h * v2_x - - 0.5 * v^2 * h_x - + 0.5 * hv_x * v - - 0.5 * h * v * v_x - + p_x - + 1.5 * p_h * b_x) + @. tmp = -(g * h_hpb_x - g * (h + b) * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) if equations.bathymetry isa BathymetryVariable @. tmp = tmp - ψ * b_x end @@ -811,9 +832,12 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus - - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat_minus) + + + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus + - + Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx) * Dmat_minus) end @trixi_timeit timer() "solving elliptic system" begin diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index 367fc0bf..3113d0d2 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -25,7 +25,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton.jl"), - equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; + gravity_constant = 9.81), tspan=(0.0, 0.1), l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], @@ -40,7 +41,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton.jl"), - equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_variable; + gravity_constant = 9.81), tspan=(0.0, 0.1), l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], @@ -68,7 +70,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), - equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; + gravity_constant = 9.81), tspan=(0.0, 0.1), l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], @@ -83,7 +86,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), - equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_variable; + gravity_constant = 9.81), tspan=(0.0, 0.1), l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], @@ -112,7 +116,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), tspan=(0.0, 0.1), - equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; + gravity_constant = 9.81), l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], @@ -127,7 +132,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), tspan=(0.0, 0.1), - equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_variable; + gravity_constant = 9.81), l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], @@ -155,7 +161,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), tspan=(0.0, 0.1), - equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; + gravity_constant = 9.81), l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], @@ -170,7 +177,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), tspan=(0.0, 0.1), - equations = SerreGreenNaghdiEquations1D(bathymetry_variable; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_variable; + gravity_constant = 9.81), l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], @@ -198,7 +206,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_dingemans.jl"), tspan=(0.0, 1.0), - equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81), + equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; + gravity_constant = 9.81), l2=[0.2463295492331802, 0.805357513755897, 0.0], linf=[0.036265772921736605, 0.11763218152350845, 0.0], cons_error=[5.684341886080802e-14, 3.509473432346624e-5, 0.0], From 2ae501231a2af4c1988aed59f049902608701b5b Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 17:06:21 +0200 Subject: [PATCH 04/34] various fixes --- src/equations/serre_green_naghdi_1d.jl | 80 ++++++++++++++++---------- 1 file changed, 50 insertions(+), 30 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 7ecd3d1d..e86ed55f 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -166,6 +166,7 @@ function create_cache(mesh, # create temporary storage h = ones(RealT, nnodes(mesh)) + b = zero(h) h_x = zero(h) v_x = zero(h) h2_x = zero(h) @@ -196,7 +197,7 @@ function create_cache(mesh, Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) factorization = cholesky(system_matrix) - cache = (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, + cache = (; h, b, h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, D, Dmat_minus, factorization) @@ -204,7 +205,7 @@ function create_cache(mesh, if D isa FourierDerivativeOperator Dmat = Matrix(D) - cache = (; h_x, v_x, h2_x, hv_x, v2_x, + cache = (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, D, Dmat) @@ -225,7 +226,7 @@ function create_cache(mesh, factorization = cholesky(system_matrix) - cache = (; h_x, v_x, h2_x, hv_x, v2_x, + cache = (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, D, Dmat, factorization) @@ -347,8 +348,8 @@ function create_cache(mesh, end if equations.bathymetry isa BathymetryVariable - ψ = zero(h) - cache = (; cache..., ψ) + psi = zero(h) + cache = (; cache..., psi) end return cache @@ -387,14 +388,18 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - h, v = q.x - dh, dv = dq.x + eta, v = q.x + dh, dv, dD = dq.x # dh = deta since b is constant in time + fill!(dD, zero(eltype(dD))) @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below - (; h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, + (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache + @. b = equations.eta0 - q.x[3] + @. h = eta - b + mul!(h_x, D, h) mul!(v_x, D, v) @. tmp = h^2 @@ -487,16 +492,19 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - h, v = q.x - dh, dv, dD = dq.x + eta, v = q.x + dh, dv, dD = dq.x # dh = deta since b is constant in time fill!(dD, zero(eltype(dD))) @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below - (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, + (; h, b, h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache + @. b = equations.eta0 - q.x[3] + @. h = eta - b + mul!(h_x, D, h) mul!(v_x, D, v) mul!(v_x_upwind, D_upwind.minus, v) @@ -586,7 +594,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. eta, v = q.x - dh, dv, dD = dq.x + dh, dv, dD = dq.x # dh = deta since b is constant in time fill!(dD, zero(eltype(dD))) @trixi_timeit timer() "hyperbolic terms" begin @@ -595,7 +603,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, h2_v_vx_x, h_vx_x, p_h, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache if equations.bathymetry isa BathymetryVariable - (; ψ) = cache + (; psi) = cache end @. b = equations.eta0 - q.x[3] @@ -625,17 +633,17 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, mul!(p_x, D, tmp) @. p_h += 0.25 * p_x if equations.bathymetry isa BathymetryVariable - @. ψ = 0.125 * p_x + @. psi = 0.125 * p_x end @. tmp = b_x * v mul!(p_x, D, tmp) @. p_h += 0.25 * h * v * p_x if equations.bathymetry isa BathymetryVariable - @. ψ += 0.125 * h * v * p_x + @. psi += 0.125 * h * v * p_x end @. p_h = p_h - 0.25 * (h_x * v + h * v_x) * b_x * v if equations.bathymetry isa BathymetryVariable - @. ψ = ψ - 0.125 * (h_x * v + h * v_x) * b_x * v + @. psi = psi - 0.125 * (h_x * v + h * v_x) * b_x * v end @. tmp = p_h * h mul!(p_x, D, tmp) @@ -661,7 +669,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, + p_x + 1.5 * p_h * b_x) if equations.bathymetry isa BathymetryVariable - @. tmp = tmp - ψ * b_x + @. tmp = tmp - psi * b_x end end @@ -726,9 +734,9 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - eta, v, b = q.x - dh, dv, db = dq.x - fill!(db, zero(eltype(db))) + eta, v = q.x + dh, dv, dD = dq.x # dh = deta since b is constant in time + fill!(dD, zero(eltype(dD))) @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below @@ -736,7 +744,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache if equations.bathymetry isa BathymetryVariable - (; ψ) = cache + (; psi) = cache end @. b = equations.eta0 - q.x[3] @@ -766,13 +774,13 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, mul!(p_x, D, tmp) @. p_h += 0.25 * p_x if equations.bathymetry isa BathymetryVariable - @. ψ = 0.125 * p_x + @. psi = 0.125 * p_x end @. tmp = b_x * v mul!(p_x, D, tmp) @. p_h += 0.25 * h * v * p_x if equations.bathymetry isa BathymetryVariable - @. ψ += 0.125 * h * v * p_x + @. psi += 0.125 * h * v * p_x end @. p_0 = p_h * h mul!(p_x, D, p_0) @@ -781,7 +789,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, - 0.25 * (h_x * v + h * v_x) * b_x * v) if equations.bathymetry isa BathymetryVariable - @. ψ = 0.125 * (h_x * v + h * v_x) * b_x * v + @. psi = psi - 0.125 * (h_x * v + h * v_x) * b_x * v end @. p_h = p_h + tmp @. tmp = tmp * h @@ -808,7 +816,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, + p_x + 1.5 * p_h * b_x) if equations.bathymetry isa BathymetryVariable - @. tmp = tmp - ψ * b_x + @. tmp = tmp - psi * b_x end end @@ -910,7 +918,15 @@ is a conserved quantity (for periodic boundary conditions). For a [`bathymetry_flat`](@ref) the total energy is given by ```math -\\frac{1}{2} g h^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2. +\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2. +``` +For a [`bathymetry_mild_slope`](@ref) the total energy is given by +```math +\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h (-h v_x + 1.5 v b_x)^2. +``` +For a [`bathymetry_variable`](@ref) the total energy has the additional term +```math ++ \\frac{1}{8} h (v b_x)^2. ``` `q_global` is a vector of the primitive variables at ALL nodes. @@ -921,11 +937,15 @@ function energy_total_modified(q_global, cache) # unpack physical parameters and SBP operator `D` g = equations.gravity - (; D, v_x) = cache + (; D, h, b, v_x) = cache # `q_global` is an `ArrayPartition`. It collects the individual arrays for - # the water height `h` and the velocity `v`. - h, v = q_global.x + # the total water height `eta = h + b` and the velocity `v`. + eta, v = q_global.x + let D = q_global.x[3] + @. b = equations.eta0 - D + @. h = eta - b + end N = length(v) e = zeros(eltype(q_global), N) @@ -951,7 +971,7 @@ function energy_total_modified(q_global, end end - @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2 + @. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2 if equations.bathymetry isa BathymetryVariable @. e += 1 / 8 * h * (v * b_x)^2 end From f32caeed7b1e4428f6e5931154c3728acc29168c Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 18:39:10 +0200 Subject: [PATCH 05/34] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index e86ed55f..36abba78 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -59,10 +59,11 @@ References for the Serre-Green-Naghdi system can be found in The semidiscretization implemented here conserves - the total water mass (integral of h) as a linear invariant -- the total momentum (integral of h v) as a nonlinear invariant -- the total energy +- the total momentum (integral of h v) as a nonlinear invariant if the bathymetry is constant +- the total modified energy -for periodic boundary conditions, see +for periodic boundary conditions (see Ranocha and Ricchiuto (2024)). +Additionally, it is well-balanced for the lake-at-rest stationary solution, see - Hendrik Ranocha and Mario Ricchiuto (2024) Structure-preserving approximations of the Serre-Green-Naghdi equations in standard and hyperbolic form From 6a4e907f27b560cac0868ab64c441b7f85af0fd3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 18:36:19 +0200 Subject: [PATCH 06/34] move function waterheight --- src/equations/bbm_bbm_1d.jl | 4 ---- src/equations/bbm_bbm_variable_bathymetry_1d.jl | 4 ---- src/equations/equations.jl | 14 +++++++++----- src/equations/serre_green_naghdi_1d.jl | 4 ---- src/equations/svaerd_kalisch_1d.jl | 4 ---- 5 files changed, 9 insertions(+), 21 deletions(-) diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 594607db..18e1ffc2 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -288,8 +288,4 @@ end return equations.eta0 - equations.D end -@inline function waterheight(q, equations::BBMBBMEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - @inline entropy(q, equations::BBMBBMEquations1D) = energy_total(q, equations) diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index 846a2ad5..8c20f39c 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -374,10 +374,6 @@ end return equations.eta0 - D end -@inline function waterheight(q, equations::BBMBBMVariableEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - @inline entropy(q, equations::BBMBBMVariableEquations1D) = energy_total(q, equations) # Calculate the error for the "lake-at-rest" test case where eta should diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 465e8f3a..949e3df3 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -96,8 +96,8 @@ function cons2prim end waterheight_total(q, equations) Return the total waterheight of the primitive variables `q` for a given set of -`equations`, i.e., the [`waterheight`](@ref) plus the -[`bathymetry`](@ref). +`equations`, i.e., the [`waterheight`](@ref) ``h`` plus the +[`bathymetry`](@ref) ``b``. `q` is a vector of the primitive variables at a single node, i.e., a vector of the correct length `nvariables(equations)`. @@ -107,15 +107,19 @@ function waterheight_total end varnames(::typeof(waterheight_total), equations) = ("η",) """ - waterheight(q, equations) + waterheight(q, equations::AbstractShallowWaterEquations) Return the waterheight of the primitive variables `q` for a given set of -`equations`, i.e., the waterheight above the bathymetry. +`equations`, i.e., the waterheight ``h`` above the bathymetry ``b``. `q` is a vector of the primitive variables at a single node, i.e., a vector of the correct length `nvariables(equations)`. + +See also [`waterheight_total`](@ref), [`bathymetry`](@ref). """ -function waterheight end +@inline function waterheight(q, equations::AbstractShallowWaterEquations) + return waterheight_total(q, equations) - bathymetry(q, equations) +end varnames(::typeof(waterheight), equations) = ("h",) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 36abba78..f7c39ee3 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -902,10 +902,6 @@ end return equations.eta0 - D end -@inline function waterheight(q, equations::SerreGreenNaghdiEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - # The entropy/energy takes the whole `q` for every point in space """ energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index f07dc6fd..34637c4e 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -326,10 +326,6 @@ end return equations.eta0 - D end -@inline function waterheight(q, equations::SvaerdKalischEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - @inline entropy(u, equations::SvaerdKalischEquations1D) = energy_total(u, equations) # The modified entropy/total energy takes the whole `q` for every point in space From 22a354a3851332f6a74f3e62477de3b20075ba42 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 18:41:18 +0200 Subject: [PATCH 07/34] bathymetry -> bathymetry_type --- src/equations/serre_green_naghdi_1d.jl | 60 +++++++++++++------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index f7c39ee3..5036d837 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -71,7 +71,7 @@ Additionally, it is well-balanced for the lake-at-rest stationary solution, see """ struct SerreGreenNaghdiEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <: AbstractSerreGreenNaghdiEquations{1, 3} - bathymetry::Bathymetry # type of bathymetry + bathymetry_type::Bathymetry # type of bathymetry gravity::RealT # gravitational constant eta0::RealT # constant still-water surface end @@ -83,11 +83,11 @@ end function get_name(equations::SerreGreenNaghdiEquations1D) name = equations |> typeof |> nameof |> string - bathymetry = equations.bathymetry - if bathymetry isa BathymetryFlat + bathymetry_type = equations.bathymetry_type + if bathymetry_type isa BathymetryFlat return name - else # variable bathymetry - return name * "-" * string(nameof(typeof(bathymetry))) + else # variable bathymetry_type + return name * "-" * string(nameof(typeof(bathymetry_type))) end end @@ -274,12 +274,12 @@ function create_cache(mesh, # is not necessarily perfectly symmetric but only up to # round-off errors. We wrap it here to avoid issues with the # factorization. - if equations.bathymetry isa BathymetryMildSlope + if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 - elseif equations.bathymetry isa BathymetryVariable + elseif equations.bathymetry_type isa BathymetryVariable factor = 1.0 else - throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry)")) + throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry_type)")) end @. M_h_p_h_bx2 = h + factor * h * b_x^2 scale_by_mass_matrix!(M_h_p_h_bx2, D) @@ -317,12 +317,12 @@ function create_cache(mesh, # is not necessarily perfectly symmetric but only up to # round-off errors. We wrap it here to avoid issues with the # factorization. - if equations.bathymetry isa BathymetryMildSlope + if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 - elseif equations.bathymetry isa BathymetryVariable + elseif equations.bathymetry_type isa BathymetryVariable factor = 1.0 else - throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry)")) + throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry_type)")) end @. M_h_p_h_bx2 = h + factor * h * b_x^2 scale_by_mass_matrix!(M_h_p_h_bx2, D) @@ -348,7 +348,7 @@ function create_cache(mesh, end end - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable psi = zero(h) cache = (; cache..., psi) end @@ -373,9 +373,9 @@ function rhs!(dq, q, t, mesh, source_terms::Nothing, solver, cache) if cache.D isa PeriodicUpwindOperators - rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry) + rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type) else - rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry) + rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type) end return nothing @@ -603,7 +603,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable (; psi) = cache end @@ -633,17 +633,17 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, @. tmp = h * b_x * v^2 mul!(p_x, D, tmp) @. p_h += 0.25 * p_x - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. psi = 0.125 * p_x end @. tmp = b_x * v mul!(p_x, D, tmp) @. p_h += 0.25 * h * v * p_x - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. psi += 0.125 * h * v * p_x end @. p_h = p_h - 0.25 * (h_x * v + h * v_x) * b_x * v - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. psi = psi - 0.125 * (h_x * v + h * v_x) * b_x * v end @. tmp = p_h * h @@ -669,15 +669,15 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, 0.5 * h * v * v_x + p_x + 1.5 * p_h * b_x) - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. tmp = tmp - psi * b_x end end @trixi_timeit timer() "assembling elliptic operator" begin - if equations.bathymetry isa BathymetryMildSlope + if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 - else # equations.bathymetry isa BathymetryVariable + else # equations.bathymetry_type isa BathymetryVariable factor = 1.0 end # The code below is equivalent to @@ -744,7 +744,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable (; psi) = cache end @@ -774,13 +774,13 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, @. tmp = h * b_x * v^2 mul!(p_x, D, tmp) @. p_h += 0.25 * p_x - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. psi = 0.125 * p_x end @. tmp = b_x * v mul!(p_x, D, tmp) @. p_h += 0.25 * h * v * p_x - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. psi += 0.125 * h * v * p_x end @. p_0 = p_h * h @@ -789,7 +789,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, @. tmp = (0.5 * h * (h * v_x + h_x * v) * v_x_upwind - 0.25 * (h_x * v + h * v_x) * b_x * v) - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. psi = psi - 0.125 * (h_x * v + h * v_x) * b_x * v end @. p_h = p_h + tmp @@ -816,15 +816,15 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, 0.5 * h * v * v_x + p_x + 1.5 * p_h * b_x) - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. tmp = tmp - psi * b_x end end @trixi_timeit timer() "assembling elliptic operator" begin - if equations.bathymetry isa BathymetryMildSlope + if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 - else # equations.bathymetry isa BathymetryVariable + else # equations.bathymetry_type isa BathymetryVariable factor = 1.0 end # The code below is equivalent to @@ -955,7 +955,7 @@ function energy_total_modified(q_global, mul!(v_x, D, v) end - if equations.bathymetry isa BathymetryFlat + if equations.bathymetry_type isa BathymetryFlat b_x = cache.tmp fill!(b_x, zero(eltype(b_x))) else @@ -969,7 +969,7 @@ function energy_total_modified(q_global, end @. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2 - if equations.bathymetry isa BathymetryVariable + if equations.bathymetry_type isa BathymetryVariable @. e += 1 / 8 * h * (v * b_x)^2 end From 59b2a9d140118d8aea9b5817916690cb6542b059 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 19:01:14 +0200 Subject: [PATCH 08/34] add bathymetry_type to elixirs and well-balanced test --- ...bm_variable_bathymetry_1d_well_balanced.jl | 1 - .../serre_green_naghdi_dingemans.jl | 4 +- .../serre_green_naghdi_soliton.jl | 4 +- .../serre_green_naghdi_soliton_fourier.jl | 4 +- .../serre_green_naghdi_soliton_relaxation.jl | 4 +- .../serre_green_naghdi_soliton_upwind.jl | 4 +- .../serre_green_naghdi_well_balanced.jl | 66 +++++++++++++++++++ .../bbm_bbm_variable_bathymetry_1d.jl | 7 -- src/equations/equations.jl | 13 ++++ src/equations/svaerd_kalisch_1d.jl | 7 -- test/test_serre_green_naghdi_1d.jl | 58 +++++++++++----- 11 files changed, 133 insertions(+), 39 deletions(-) create mode 100644 examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl index ab640e15..08b8d26c 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl @@ -57,7 +57,6 @@ analysis_callback = AnalysisCallback(semi; interval = 10, extra_analysis_integrals = (waterheight_total, velocity, entropy, lake_at_rest_error)) -# Always put relaxation_callback before analysis_callback to guarantee conservation of the invariant callbacks = CallbackSet(analysis_callback, summary_callback) dt = 0.5 diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl index 854a8879..f424fa4c 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl @@ -5,8 +5,8 @@ using SummationByPartsOperators: upwind_operators, periodic_derivative_operator ############################################################################### # Semidiscretization of the Svärd-Kalisch equations -bathymetry = bathymetry_variable # or bathymetry_mild_slope -equations = SerreGreenNaghdiEquations1D(bathymetry; +bathymetry_type = bathymetry_variable # or bathymetry_mild_slope +equations = SerreGreenNaghdiEquations1D(bathymetry_type; gravity_constant = 9.81) initial_condition = initial_condition_dingemans diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl index 033197d7..bfcaaddf 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl @@ -4,7 +4,9 @@ using DispersiveShallowWater ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl index d22196f9..b8f4bf7f 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl @@ -5,7 +5,9 @@ using SummationByPartsOperators: fourier_derivative_operator ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl index e095c9ff..d484a733 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl @@ -5,7 +5,9 @@ using SummationByPartsOperators: fourier_derivative_operator ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl index 8bb14375..cd1c3eb7 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl @@ -6,7 +6,9 @@ using SparseArrays: sparse ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl new file mode 100644 index 00000000..38e21c67 --- /dev/null +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl @@ -0,0 +1,66 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: upwind_operators, periodic_derivative_operator + +############################################################################### +# Semidiscretization of the BBM-BBM equations + +bathymetry_type = bathymetry_variable # or bathymetry_mild_slope +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 1.0, eta0 = 2.0) + +# Setup a truly discontinuous bottom topography function for this academic +# testcase of well-balancedness. The errors from the analysis callback are +# not important but the error for this lake-at-rest test case +# `∫|η-η₀|` should be around machine roundoff. +function initial_condition_discontinuous_well_balancedness(x, t, + equations::SerreGreenNaghdiEquations1D, + mesh) + # Set the background values + eta = equations.eta0 + v = 0.0 + D = equations.eta0 - 1.0 + + # Setup a discontinuous bottom topography + if x >= 0.5 && x <= 0.75 + D = equations.eta0 - 1.5 - 0.5 * sinpi(2.0 * x) + end + + return SVector(eta, v, D) +end + +initial_condition = initial_condition_discontinuous_well_balancedness +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -1.0 +coordinates_max = 1.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic upwind SBP operators +D1 = upwind_operators(periodic_derivative_operator, + derivative_order = 1, + accuracy_order = 3, # the resulting central operators have order 4 + xmin = xmin(mesh), xmax = xmax(mesh), + N = nnodes(mesh)) +solver = Solver(D1, nothing) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 30.0) +ode = semidiscretize(semi, tspan) +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + momentum, + entropy_modified, + lake_at_rest_error)) +callbacks = CallbackSet(analysis_callback, summary_callback) + +sol = solve(ode, Tsit5(), dt = 0.5, adaptive = false, callback = callbacks) diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index 8c20f39c..ffc1e1fb 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -375,10 +375,3 @@ end end @inline entropy(q, equations::BBMBBMVariableEquations1D) = energy_total(q, equations) - -# Calculate the error for the "lake-at-rest" test case where eta should -# be a constant value over time -@inline function lake_at_rest_error(q, equations::BBMBBMVariableEquations1D) - eta, _, _ = q - return abs(equations.eta0 - eta) -end diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 949e3df3..5b902a75 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -188,6 +188,19 @@ of the correct length `nvariables(equations)`. """ function bathymetry end +""" + lake_at_rest_error(q, equations::AbstractShallowWaterEquations) + +Calculate the error for the "lake-at-rest" test case where the +[`waterheight_total`](@ref) ``\\eta = h + b`` should +be a constant value over time (given by the value ``\\eta_0`` passed to the +`equations` when constructing them). +""" +@inline function lake_at_rest_error(q, equations::AbstractShallowWaterEquations) + eta = waterheight_total(q, equations) + return abs(equations.eta0 - eta) +end + """ entropy(q, equations) diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 34637c4e..6c904fcc 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -361,10 +361,3 @@ is a conserved quantity of the Svärd-Kalisch equations. end return e_modified end - -# Calculate the error for the "lake-at-rest" test case where eta should -# be a constant value over time -@inline function lake_at_rest_error(u, equations::SvaerdKalischEquations1D) - eta, _, _ = u - return abs(equations.eta0 - eta) -end diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index 3113d0d2..b4c7fb06 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -25,8 +25,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton.jl"), - equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; - gravity_constant = 9.81), + bathymetry_type = bathymetry_mild_slope, tspan=(0.0, 0.1), l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], @@ -41,8 +40,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton.jl"), - equations=SerreGreenNaghdiEquations1D(bathymetry_variable; - gravity_constant = 9.81), + bathymetry_type = bathymetry_variable, tspan=(0.0, 0.1), l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], @@ -70,8 +68,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), - equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; - gravity_constant = 9.81), + bathymetry_type = bathymetry_mild_slope, tspan=(0.0, 0.1), l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], @@ -86,8 +83,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), - equations=SerreGreenNaghdiEquations1D(bathymetry_variable; - gravity_constant = 9.81), + bathymetry_type = bathymetry_variable, tspan=(0.0, 0.1), l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], @@ -116,8 +112,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), tspan=(0.0, 0.1), - equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; - gravity_constant = 9.81), + bathymetry_type = bathymetry_mild_slope, l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], @@ -132,8 +127,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), tspan=(0.0, 0.1), - equations=SerreGreenNaghdiEquations1D(bathymetry_variable; - gravity_constant = 9.81), + bathymetry_type = bathymetry_variable, l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], @@ -161,8 +155,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), tspan=(0.0, 0.1), - equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; - gravity_constant = 9.81), + bathymetry_type = bathymetry_mild_slope, l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], @@ -177,8 +170,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), tspan=(0.0, 0.1), - equations=SerreGreenNaghdiEquations1D(bathymetry_variable; - gravity_constant = 9.81), + bathymetry_type = bathymetry_variable, l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], @@ -188,6 +180,37 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=850_000) end + @trixi_testset "serre_green_naghdi_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_well_balanced.jl"), + tspan=(0.0, 2.0), + l2=[0, 0, 0], + linf=[0, 0, 0], + cons_error=[0, 0, 0], + change_waterheight=0.0, + change_momentum=0.0, + change_entropy_modified=0.0, + lake_at_rest=0.0) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_well_balanced.jl with bathymetry_mild_slope" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_well_balanced.jl"), + bathymetry_type = bathymetry_mild_slope, + tspan=(0.0, 2.0), + l2=[0, 0, 0], + linf=[0, 0, 0], + cons_error=[0, 0, 0], + change_waterheight=0.0, + change_momentum=0.0, + change_entropy_modified=0.0, + lake_at_rest=0.0) + + @test_allocations(semi, sol, allocs=750_000) + end + @trixi_testset "serre_green_naghdi_dingemans.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_dingemans.jl"), @@ -206,8 +229,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_dingemans.jl"), tspan=(0.0, 1.0), - equations=SerreGreenNaghdiEquations1D(bathymetry_mild_slope; - gravity_constant = 9.81), + bathymetry_type = bathymetry_mild_slope, l2=[0.2463295492331802, 0.805357513755897, 0.0], linf=[0.036265772921736605, 0.11763218152350845, 0.0], cons_error=[5.684341886080802e-14, 3.509473432346624e-5, 0.0], From 05c25cd1655774901af0fea8ac267aac2cb6d134 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 19:05:08 +0200 Subject: [PATCH 09/34] move comments --- src/equations/serre_green_naghdi_1d.jl | 56 +++++++++++++------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 5036d837..ab6f78a1 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -185,14 +185,14 @@ function create_cache(mesh, Dmat_minus = sparse(D.minus) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. @. M_h = h scale_by_mass_matrix!(M_h, D) @. M_h3_3 = (1 / 3) * h^3 scale_by_mass_matrix!(M_h3_3, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. system_matrix = Symmetric(Diagonal(M_h) + Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) @@ -213,14 +213,14 @@ function create_cache(mesh, else Dmat = sparse(D) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. @. M_h = h scale_by_mass_matrix!(M_h, D) @. M_h3_3 = (1 / 3) * h^3 scale_by_mass_matrix!(M_h3_3, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. system_matrix = Symmetric(Diagonal(M_h) + Dmat' * Diagonal(M_h3_3) * Dmat) @@ -270,10 +270,6 @@ function create_cache(mesh, Dmat_minus = sparse(D.minus) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 elseif equations.bathymetry_type isa BathymetryVariable @@ -288,6 +284,10 @@ function create_cache(mesh, scale_by_mass_matrix!(M_h3_3, D) @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus @@ -313,10 +313,6 @@ function create_cache(mesh, else Dmat = sparse(D) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 elseif equations.bathymetry_type isa BathymetryVariable @@ -331,6 +327,10 @@ function create_cache(mesh, scale_by_mass_matrix!(M_h3_3, D) @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + Dmat' * (Diagonal(M_h3_3) * Dmat @@ -448,14 +448,14 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla # The code below is equivalent to # dv .= (Diagonal(h) - Dmat * Diagonal(1/3 .* h.^3) * Dmat) \ tmp # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. @. M_h = h scale_by_mass_matrix!(M_h, D) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 scale_by_mass_matrix!(M_h3_3, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off + # errors. We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h) + Dmat' * Diagonal(M_h3_3) * Dmat) @@ -556,14 +556,14 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat # The code below is equivalent to # dv .= (Diagonal(h) - Dmat_plus * Diagonal(1/3 .* h.^3) * Dmat_minus) \ tmp # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off errors. - # We wrap it here to avoid issues with the factorization. @. M_h = h scale_by_mass_matrix!(M_h, D) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 scale_by_mass_matrix!(M_h3_3, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off errors. + # We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h) + Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) @@ -683,9 +683,6 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, # The code below is equivalent to # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. @. M_h_p_h_bx2 = h + factor * h * b_x^2 scale_by_mass_matrix!(M_h_p_h_bx2, D) inv3 = 1 / 3 @@ -693,6 +690,9 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, scale_by_mass_matrix!(M_h3_3, D) @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off + # errors. We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + Dmat' * (Diagonal(M_h3_3) * Dmat @@ -830,9 +830,6 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, # The code below is equivalent to # dv .= (Diagonal(h .+ 0.75 .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. @. M_h_p_h_bx2 = h + factor * h * b_x^2 scale_by_mass_matrix!(M_h_p_h_bx2, D) inv3 = 1 / 3 @@ -840,6 +837,9 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, scale_by_mass_matrix!(M_h3_3, D) @. M_h2_bx = 0.5 * h^2 * b_x scale_by_mass_matrix!(M_h2_bx, D) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to round-off + # errors. We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus From 853662ebcf7b0f276ef94b6592f3951aaaee91c9 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 19:05:32 +0200 Subject: [PATCH 10/34] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index ab6f78a1..9c97b583 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -681,7 +681,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, factor = 1.0 end # The code below is equivalent to - # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp + # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp # but faster since the symbolic factorization is reused. @. M_h_p_h_bx2 = h + factor * h * b_x^2 scale_by_mass_matrix!(M_h_p_h_bx2, D) From dc5ed5ad753250f8cbc6703c43e01ba595eb0af2 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 16 Aug 2024 19:06:41 +0200 Subject: [PATCH 11/34] format --- test/test_serre_green_naghdi_1d.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index b4c7fb06..2b2099c6 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -25,7 +25,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton.jl"), - bathymetry_type = bathymetry_mild_slope, + bathymetry_type=bathymetry_mild_slope, tspan=(0.0, 0.1), l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], @@ -40,7 +40,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton.jl"), - bathymetry_type = bathymetry_variable, + bathymetry_type=bathymetry_variable, tspan=(0.0, 0.1), l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], @@ -68,7 +68,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), - bathymetry_type = bathymetry_mild_slope, + bathymetry_type=bathymetry_mild_slope, tspan=(0.0, 0.1), l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], @@ -83,7 +83,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), - bathymetry_type = bathymetry_variable, + bathymetry_type=bathymetry_variable, tspan=(0.0, 0.1), l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], @@ -112,7 +112,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), tspan=(0.0, 0.1), - bathymetry_type = bathymetry_mild_slope, + bathymetry_type=bathymetry_mild_slope, l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], @@ -127,7 +127,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), tspan=(0.0, 0.1), - bathymetry_type = bathymetry_variable, + bathymetry_type=bathymetry_variable, l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], @@ -155,7 +155,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), tspan=(0.0, 0.1), - bathymetry_type = bathymetry_mild_slope, + bathymetry_type=bathymetry_mild_slope, l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], @@ -170,7 +170,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), tspan=(0.0, 0.1), - bathymetry_type = bathymetry_variable, + bathymetry_type=bathymetry_variable, l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], @@ -198,7 +198,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @trixi_testset "serre_green_naghdi_well_balanced.jl with bathymetry_mild_slope" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_well_balanced.jl"), - bathymetry_type = bathymetry_mild_slope, + bathymetry_type=bathymetry_mild_slope, tspan=(0.0, 2.0), l2=[0, 0, 0], linf=[0, 0, 0], @@ -229,7 +229,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_dingemans.jl"), tspan=(0.0, 1.0), - bathymetry_type = bathymetry_mild_slope, + bathymetry_type=bathymetry_mild_slope, l2=[0.2463295492331802, 0.805357513755897, 0.0], linf=[0.036265772921736605, 0.11763218152350845, 0.0], cons_error=[5.684341886080802e-14, 3.509473432346624e-5, 0.0], From 1740e5aea64c76eb0a667592e50193000e9af6b6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:14:56 +0200 Subject: [PATCH 12/34] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 9c97b583..3a540e26 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -1,5 +1,5 @@ @doc raw""" - SerreGreenNaghdiEquations1D(bathymetry = bathymetry_flat; + SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_flat; gravity_constant, eta0 = 0.0) Serre-Green-Naghdi system in one spatial dimension. @@ -19,7 +19,7 @@ is therefore given by ``h = \eta - \eta_0 + D``. The Serre-Green-Naghdi equations are only implemented for ``\eta_0 = 0``. The total water height is therefore given by ``\eta = h + b``. -Three types of variable `bathymetry` are supported: +Three types of variable `bathymetry_type` are supported: - [`bathymetry_flat`](@ref): flat bathymetry (typically ``b = 0`` everywhere) - [`bathymetry_mild_slope`](@ref): variable bathymetry with mild-slope approximation - [`bathymetry_variable`](@ref): general variable bathymetry @@ -76,9 +76,9 @@ struct SerreGreenNaghdiEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Re eta0::RealT # constant still-water surface end -function SerreGreenNaghdiEquations1D(bathymetry = bathymetry_flat; +function SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_flat; gravity_constant, eta0 = 0.0) - SerreGreenNaghdiEquations1D(bathymetry, gravity_constant, eta0) + SerreGreenNaghdiEquations1D(bathymetry_type, gravity_constant, eta0) end function get_name(equations::SerreGreenNaghdiEquations1D) From 711864f8014a10909c9ad7ab48d26f8e919116c0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:20:12 +0200 Subject: [PATCH 13/34] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 3a540e26..1c1c5c98 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -613,7 +613,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, mul!(h_x, D, h) mul!(v_x, D, v) - @. tmp = h * (h + b) + @. tmp = h * eta mul!(h_hpb_x, D, tmp) @. tmp = h * v mul!(hv_x, D, tmp) From b1ba5e67c950a5dd6131ba433fe7114742d17181 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:20:54 +0200 Subject: [PATCH 14/34] Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 1c1c5c98..1b7b1be6 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -947,7 +947,7 @@ function energy_total_modified(q_global, N = length(v) e = zeros(eltype(q_global), N) - # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 w^2 + # 1/2 g eta^2 + 1/2 h v^2 + 1/6 h^3 w^2 # and + 1/8 h (v b_x)^2 for full bathymetry without mild-slope approximation if D isa PeriodicUpwindOperators mul!(v_x, D.minus, v) From 5578e539e745bde2eb0272c89423f14cd74753fd Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:30:32 +0200 Subject: [PATCH 15/34] use D1 for SBP op. --- src/equations/serre_green_naghdi_1d.jl | 275 +++++++++++++------------ 1 file changed, 138 insertions(+), 137 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 1b7b1be6..d5725f1c 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -121,7 +121,8 @@ end The initial condition that uses the dispersion relation of the Euler equations to approximate waves generated by a wave maker as it is done by experiments of -Dingemans. The topography is a trapezoidal. It is assumed that `equations.eta0 = 0.8`. +Dingemans for a flow over a trapezoidal bathymetry. +It is assumed that `equations.eta0 = 0.8`. References: - Magnus Svärd, Henrik Kalisch (2023) @@ -163,7 +164,7 @@ function create_cache(mesh, initial_condition, ::BoundaryConditionPeriodic, RealT, uEltype) - D = solver.D1 + D1 = solver.D1 # create temporary storage h = ones(RealT, nnodes(mesh)) @@ -180,57 +181,57 @@ function create_cache(mesh, M_h = zero(h) M_h3_3 = zero(h) - if D isa PeriodicUpwindOperators + if D1 isa PeriodicUpwindOperators v_x_upwind = zero(h) - Dmat_minus = sparse(D.minus) + D1mat_minus = sparse(D1.minus) @. M_h = h - scale_by_mass_matrix!(M_h, D) + scale_by_mass_matrix!(M_h, D1) @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to # round-off errors. We wrap it here to avoid issues with the # factorization. system_matrix = Symmetric(Diagonal(M_h) + - Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) + D1mat_minus' * Diagonal(M_h3_3) * D1mat_minus) factorization = cholesky(system_matrix) cache = (; h, b, h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, - D, Dmat_minus, factorization) + D1, D1mat_minus, factorization) else - if D isa FourierDerivativeOperator - Dmat = Matrix(D) + if D1 isa FourierDerivativeOperator + D1mat = Matrix(D1) cache = (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, - D, Dmat) + D1, D1mat) else - Dmat = sparse(D) + D1mat = sparse(D1) @. M_h = h - scale_by_mass_matrix!(M_h, D) + scale_by_mass_matrix!(M_h, D1) @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to # round-off errors. We wrap it here to avoid issues with the # factorization. system_matrix = Symmetric(Diagonal(M_h) + - Dmat' * Diagonal(M_h3_3) * Dmat) + D1mat' * Diagonal(M_h3_3) * D1mat) factorization = cholesky(system_matrix) cache = (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, - D, Dmat, factorization) + D1, D1mat, factorization) end end @@ -244,7 +245,7 @@ function create_cache(mesh, initial_condition, ::BoundaryConditionPeriodic, RealT, uEltype) - D = solver.D1 + D1 = solver.D1 # create temporary storage h = ones(RealT, nnodes(mesh)) @@ -264,11 +265,11 @@ function create_cache(mesh, M_h3_3 = zero(h) M_h2_bx = zero(h) - if D isa PeriodicUpwindOperators + if D1 isa PeriodicUpwindOperators v_x_upwind = zero(h) p_0 = zero(h) - Dmat_minus = sparse(D.minus) + D1mat_minus = sparse(D1.minus) if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 @@ -278,40 +279,40 @@ function create_cache(mesh, throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry_type)")) end @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D) + scale_by_mass_matrix!(M_h_p_h_bx2, D1) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D) + scale_by_mass_matrix!(M_h2_bx, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to # round-off errors. We wrap it here to avoid issues with the # factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + - Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus + D1mat_minus' * (Diagonal(M_h3_3) * D1mat_minus - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat_minus) + Diagonal(M_h2_bx) * D1mat_minus) factorization = cholesky(system_matrix) cache = (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx, - D, Dmat_minus, factorization) + D1, D1mat_minus, factorization) else - if D isa FourierDerivativeOperator - Dmat = Matrix(D) + if D1 isa FourierDerivativeOperator + D1mat = Matrix(D1) cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx, - D, Dmat) + D1, D1mat) else - Dmat = sparse(D) + D1mat = sparse(D1) if equations.bathymetry_type isa BathymetryMildSlope factor = 0.75 @@ -321,30 +322,30 @@ function create_cache(mesh, throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry_type)")) end @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D) + scale_by_mass_matrix!(M_h_p_h_bx2, D1) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D) + scale_by_mass_matrix!(M_h2_bx, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to # round-off errors. We wrap it here to avoid issues with the # factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + - Dmat' * (Diagonal(M_h3_3) * Dmat + D1mat' * (Diagonal(M_h3_3) * D1mat - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat) + Diagonal(M_h2_bx) * D1mat) factorization = cholesky(system_matrix) cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_h, p_x, tmp, M_h_p_h_bx2, M_h3_3, M_h2_bx, - D, Dmat, factorization) + D1, D1mat, factorization) end end @@ -372,7 +373,7 @@ function rhs!(dq, q, t, mesh, ::BoundaryConditionPeriodic, source_terms::Nothing, solver, cache) - if cache.D isa PeriodicUpwindOperators + if cache.D1 isa PeriodicUpwindOperators rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type) else rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type) @@ -382,14 +383,14 @@ function rhs!(dq, q, t, mesh, end function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP operator in sparse matrix form `Dmat` + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP operator in sparse matrix form `D1mat` g = gravity_constant(equations) - (; D, Dmat) = cache + (; D1, D1mat) = cache # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - eta, v = q.x + eta, v, D = q.x dh, dv, dD = dq.x # dh = deta since b is constant in time fill!(dD, zero(eltype(dD))) @@ -398,29 +399,29 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache - @. b = equations.eta0 - q.x[3] + @. b = equations.eta0 - D @. h = eta - b - mul!(h_x, D, h) - mul!(v_x, D, v) + mul!(h_x, D1, h) + mul!(v_x, D1, v) @. tmp = h^2 - mul!(h2_x, D, tmp) + mul!(h2_x, D1, tmp) @. tmp = h * v - mul!(hv_x, D, tmp) + mul!(hv_x, D1, tmp) @. tmp = v^2 - mul!(v2_x, D, tmp) + mul!(v2_x, D1, tmp) @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) + mul!(h2_v_vx_x, D1, tmp) @. tmp = h * v_x - mul!(h_vx_x, D, tmp) + mul!(h_vx_x, D1, tmp) inv6 = 1 / 6 @. tmp = (0.5 * h^2 * (h * v_x + h_x * v) * v_x - inv6 * h * h2_v_vx_x - inv6 * h^2 * v * h_vx_x) - mul!(p_x, D, tmp) + mul!(p_x, D1, tmp) # Plain: h_t + (h v)_x = 0 # @@ -446,19 +447,19 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla @trixi_timeit timer() "assembling elliptic operator" begin # The code below is equivalent to - # dv .= (Diagonal(h) - Dmat * Diagonal(1/3 .* h.^3) * Dmat) \ tmp + # dv .= (Diagonal(h) - D1mat * Diagonal(1/3 .* h.^3) * D1mat) \ tmp # but faster since the symbolic factorization is reused. @. M_h = h - scale_by_mass_matrix!(M_h, D) + scale_by_mass_matrix!(M_h, D1) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to round-off # errors. We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h) + - Dmat' * Diagonal(M_h3_3) * Dmat) + D1mat' * Diagonal(M_h3_3) * D1mat) end @trixi_timeit timer() "solving elliptic system" begin @@ -466,7 +467,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla (; factorization) = cache cholesky!(factorization, system_matrix; check = false) if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) + scale_by_mass_matrix!(tmp, D1) dv .= factorization \ tmp else # The factorization may fail if the time step is too large @@ -475,7 +476,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla end else factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D) + scale_by_mass_matrix!(tmp, D1) ldiv!(dv, factorization, tmp) end end @@ -484,16 +485,16 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla end function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP upwind operator in sparse matrix form `Dmat_minus` + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP upwind operator in sparse matrix form `D1mat_minus` g = gravity_constant(equations) - (; Dmat_minus) = cache - D_upwind = cache.D - D = D_upwind.central + (; D1mat_minus) = cache + D1_upwind = cache.D1 + D1 = D1_upwind.central # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - eta, v = q.x + eta, v, D = q.x dh, dv, dD = dq.x # dh = deta since b is constant in time fill!(dD, zero(eltype(dD))) @@ -503,32 +504,32 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache - @. b = equations.eta0 - q.x[3] + @. b = equations.eta0 - D @. h = eta - b - mul!(h_x, D, h) - mul!(v_x, D, v) - mul!(v_x_upwind, D_upwind.minus, v) + mul!(h_x, D1, h) + mul!(v_x, D1, v) + mul!(v_x_upwind, D1_upwind.minus, v) @. tmp = h^2 - mul!(h2_x, D, tmp) + mul!(h2_x, D1, tmp) @. tmp = h * v - mul!(hv_x, D, tmp) + mul!(hv_x, D1, tmp) @. tmp = v^2 - mul!(v2_x, D, tmp) + mul!(v2_x, D1, tmp) @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) + mul!(h2_v_vx_x, D1, tmp) @. tmp = h * v_x - mul!(h_vx_x, D, tmp) + mul!(h_vx_x, D1, tmp) # p_+ @. tmp = 0.5 * h^2 * (h * v_x + h_x * v) * v_x_upwind - mul!(p_x, D_upwind.plus, tmp) + mul!(p_x, D1_upwind.plus, tmp) # p_0 minv6 = -1 / 6 @. tmp = minv6 * (h * h2_v_vx_x + h^2 * v * h_vx_x) - mul!(p_x, D, tmp, 1.0, 1.0) + mul!(p_x, D1, tmp, 1.0, 1.0) # Plain: h_t + (h v)_x = 0 # @@ -554,26 +555,26 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat @trixi_timeit timer() "assembling elliptic operator" begin # The code below is equivalent to - # dv .= (Diagonal(h) - Dmat_plus * Diagonal(1/3 .* h.^3) * Dmat_minus) \ tmp + # dv .= (Diagonal(h) - D1mat_plus * Diagonal(1/3 .* h.^3) * D1mat_minus) \ tmp # but faster since the symbolic factorization is reused. @. M_h = h - scale_by_mass_matrix!(M_h, D) + scale_by_mass_matrix!(M_h, D1) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to round-off errors. # We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h) + - Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) + D1mat_minus' * Diagonal(M_h3_3) * D1mat_minus) end @trixi_timeit timer() "solving elliptic system" begin (; factorization) = cache cholesky!(factorization, system_matrix; check = false) if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) + scale_by_mass_matrix!(tmp, D1) dv .= factorization \ tmp else # The factorization may fail if the time step is too large @@ -587,14 +588,14 @@ end function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::Union{BathymetryMildSlope, BathymetryVariable}) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP operator in sparse matrix form `Dmat` + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP operator in sparse matrix form `D1mat` g = gravity_constant(equations) - (; D, Dmat) = cache + (; D1, D1mat) = cache # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - eta, v = q.x + eta, v, D = q.x dh, dv, dD = dq.x # dh = deta since b is constant in time fill!(dD, zero(eltype(dD))) @@ -607,23 +608,23 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, (; psi) = cache end - @. b = equations.eta0 - q.x[3] + @. b = equations.eta0 - D @. h = eta - b - mul!(b_x, D, b) + mul!(b_x, D1, b) - mul!(h_x, D, h) - mul!(v_x, D, v) + mul!(h_x, D1, h) + mul!(v_x, D1, v) @. tmp = h * eta - mul!(h_hpb_x, D, tmp) + mul!(h_hpb_x, D1, tmp) @. tmp = h * v - mul!(hv_x, D, tmp) + mul!(hv_x, D1, tmp) @. tmp = v^2 - mul!(v2_x, D, tmp) + mul!(v2_x, D1, tmp) @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) + mul!(h2_v_vx_x, D1, tmp) @. tmp = h * v_x - mul!(h_vx_x, D, tmp) + mul!(h_vx_x, D1, tmp) inv6 = 1 / 6 @. p_h = (0.5 * h * (h * v_x + h_x * v) * v_x - @@ -631,13 +632,13 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, - inv6 * h * v * h_vx_x) @. tmp = h * b_x * v^2 - mul!(p_x, D, tmp) + mul!(p_x, D1, tmp) @. p_h += 0.25 * p_x if equations.bathymetry_type isa BathymetryVariable @. psi = 0.125 * p_x end @. tmp = b_x * v - mul!(p_x, D, tmp) + mul!(p_x, D1, tmp) @. p_h += 0.25 * h * v * p_x if equations.bathymetry_type isa BathymetryVariable @. psi += 0.125 * h * v * p_x @@ -647,7 +648,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, @. psi = psi - 0.125 * (h_x * v + h * v_x) * b_x * v end @. tmp = p_h * h - mul!(p_x, D, tmp) + mul!(p_x, D1, tmp) # Plain: h_t + (h v)_x = 0 # @@ -681,25 +682,25 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, factor = 1.0 end # The code below is equivalent to - # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp + # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp # but faster since the symbolic factorization is reused. @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D) + scale_by_mass_matrix!(M_h_p_h_bx2, D1) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D) + scale_by_mass_matrix!(M_h2_bx, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to round-off # errors. We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + - Dmat' * (Diagonal(M_h3_3) * Dmat + D1mat' * (Diagonal(M_h3_3) * D1mat - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat) + Diagonal(M_h2_bx) * D1mat) end @trixi_timeit timer() "solving elliptic system" begin @@ -707,7 +708,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, (; factorization) = cache cholesky!(factorization, system_matrix; check = false) if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) + scale_by_mass_matrix!(tmp, D1) dv .= factorization \ tmp else # The factorization may fail if the time step is too large @@ -716,7 +717,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end else factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D) + scale_by_mass_matrix!(tmp, D1) ldiv!(dv, factorization, tmp) end end @@ -726,16 +727,16 @@ end function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::Union{BathymetryMildSlope, BathymetryVariable}) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP operator in sparse matrix form `Dmat` + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP operator in sparse matrix form `D1mat` g = gravity_constant(equations) - (; Dmat_minus) = cache - D_upwind = cache.D - D = D_upwind.central + (; D1mat_minus) = cache + D1_upwind = cache.D1 + D1 = D1_upwind.central # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - eta, v = q.x + eta, v, D = q.x dh, dv, dD = dq.x # dh = deta since b is constant in time fill!(dD, zero(eltype(dD))) @@ -748,43 +749,43 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, (; psi) = cache end - @. b = equations.eta0 - q.x[3] + @. b = equations.eta0 - D @. h = eta - b - mul!(b_x, D, b) + mul!(b_x, D1, b) - mul!(h_x, D, h) - mul!(v_x, D, v) - mul!(v_x_upwind, D_upwind.minus, v) + mul!(h_x, D1, h) + mul!(v_x, D1, v) + mul!(v_x_upwind, D1_upwind.minus, v) @. tmp = h * (h + b) - mul!(h_hpb_x, D, tmp) + mul!(h_hpb_x, D1, tmp) @. tmp = h * v - mul!(hv_x, D, tmp) + mul!(hv_x, D1, tmp) @. tmp = v^2 - mul!(v2_x, D, tmp) + mul!(v2_x, D1, tmp) @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) + mul!(h2_v_vx_x, D1, tmp) @. tmp = h * v_x - mul!(h_vx_x, D, tmp) + mul!(h_vx_x, D1, tmp) # p_0 minv6 = -1 / 6 @. p_h = minv6 * (h2_v_vx_x + h * v * h_vx_x) @. tmp = h * b_x * v^2 - mul!(p_x, D, tmp) + mul!(p_x, D1, tmp) @. p_h += 0.25 * p_x if equations.bathymetry_type isa BathymetryVariable @. psi = 0.125 * p_x end @. tmp = b_x * v - mul!(p_x, D, tmp) + mul!(p_x, D1, tmp) @. p_h += 0.25 * h * v * p_x if equations.bathymetry_type isa BathymetryVariable @. psi += 0.125 * h * v * p_x end @. p_0 = p_h * h - mul!(p_x, D, p_0) + mul!(p_x, D1, p_0) # p_+ @. tmp = (0.5 * h * (h * v_x + h_x * v) * v_x_upwind - @@ -794,7 +795,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end @. p_h = p_h + tmp @. tmp = tmp * h - mul!(p_x, D_upwind.plus, tmp, 1.0, 1.0) + mul!(p_x, D1_upwind.plus, tmp, 1.0, 1.0) # Plain: h_t + (h v)_x = 0 # @@ -828,25 +829,25 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, factor = 1.0 end # The code below is equivalent to - # dv .= (Diagonal(h .+ 0.75 .* h .* b_x.^2) - Dmat * (Diagonal(1/3 .* h.^3) * Dmat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * Dmat) \ tmp + # dv .= (Diagonal(h .+ 0.75 .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp # but faster since the symbolic factorization is reused. @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D) + scale_by_mass_matrix!(M_h_p_h_bx2, D1) inv3 = 1 / 3 @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) + scale_by_mass_matrix!(M_h3_3, D1) @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D) + scale_by_mass_matrix!(M_h2_bx, D1) # Floating point errors accumulate a bit and the system matrix # is not necessarily perfectly symmetric but only up to round-off # errors. We wrap it here to avoid issues with the factorization. system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + - Dmat_minus' * (Diagonal(M_h3_3) * Dmat_minus + D1mat_minus' * (Diagonal(M_h3_3) * D1mat_minus - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * Dmat_minus) + Diagonal(M_h2_bx) * D1mat_minus) end @trixi_timeit timer() "solving elliptic system" begin @@ -854,7 +855,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, (; factorization) = cache cholesky!(factorization, system_matrix; check = false) if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) + scale_by_mass_matrix!(tmp, D1) dv .= factorization \ tmp else # The factorization may fail if the time step is too large @@ -863,7 +864,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end else factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D) + scale_by_mass_matrix!(tmp, D1) ldiv!(dv, factorization, tmp) end end @@ -932,9 +933,9 @@ For a [`bathymetry_variable`](@ref) the total energy has the additional term function energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache) - # unpack physical parameters and SBP operator `D` + # unpack physical parameters and SBP operator `D1` g = equations.gravity - (; D, h, b, v_x) = cache + (; D1, h, b, v_x) = cache # `q_global` is an `ArrayPartition`. It collects the individual arrays for # the total water height `eta = h + b` and the velocity `v`. @@ -949,10 +950,10 @@ function energy_total_modified(q_global, # 1/2 g eta^2 + 1/2 h v^2 + 1/6 h^3 w^2 # and + 1/8 h (v b_x)^2 for full bathymetry without mild-slope approximation - if D isa PeriodicUpwindOperators - mul!(v_x, D.minus, v) + if D1 isa PeriodicUpwindOperators + mul!(v_x, D1.minus, v) else - mul!(v_x, D, v) + mul!(v_x, D1, v) end if equations.bathymetry_type isa BathymetryFlat @@ -961,10 +962,10 @@ function energy_total_modified(q_global, else (; b, b_x) = cache @. b = equations.eta0 - q_global.x[3] - if D isa PeriodicUpwindOperators - mul!(b_x, D.central, b) + if D1 isa PeriodicUpwindOperators + mul!(b_x, D1.central, b) else - mul!(b_x, D, b) + mul!(b_x, D1, b) end end From c1cbda49bba8992b6df2a751cae857752a680d5c Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:32:27 +0200 Subject: [PATCH 16/34] (h + b) -> eta --- src/equations/serre_green_naghdi_1d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index d5725f1c..8ce9839c 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -659,7 +659,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -(g * h_hpb_x - g * (h + b) * h_x + @. tmp = -(g * h_hpb_x - g * eta * h_x + 0.5 * h * v2_x - @@ -756,7 +756,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, mul!(h_x, D1, h) mul!(v_x, D1, v) mul!(v_x_upwind, D1_upwind.minus, v) - @. tmp = h * (h + b) + @. tmp = h * eta mul!(h_hpb_x, D1, tmp) @. tmp = h * v mul!(hv_x, D1, tmp) @@ -806,7 +806,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, # Plain: h v_t + ... = 0 # # Split form for energy conservation: - @. tmp = -(g * h_hpb_x - g * (h + b) * h_x + @. tmp = -(g * h_hpb_x - g * eta * h_x + 0.5 * h * v2_x - From 33115e34646ba68dc8e6a4537337f07ee8f2c437 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:34:03 +0200 Subject: [PATCH 17/34] explain mild slope --- src/equations/equations.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 5b902a75..19b6e3aa 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -358,7 +358,8 @@ struct BathymetryMildSlope <: AbstractBathymetry end bathymetry_mild_slope = DispersiveShallowWater.BathymetryMildSlope() A singleton struct indicating a variable bathymetry with mild-slope -approximation. +approximation. Typically, this means that some terms like ``b_x^2`` +are neglected. See also [`bathymetry_flat`](@ref) and [`bathymetry_variable`](@ref). """ From e321a8df549122a1685f3ddb0c16567d02dd8bc9 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:37:57 +0200 Subject: [PATCH 18/34] -= --- src/equations/serre_green_naghdi_1d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 8ce9839c..96d62e61 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -645,7 +645,7 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end @. p_h = p_h - 0.25 * (h_x * v + h * v_x) * b_x * v if equations.bathymetry_type isa BathymetryVariable - @. psi = psi - 0.125 * (h_x * v + h * v_x) * b_x * v + @. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v end @. tmp = p_h * h mul!(p_x, D1, tmp) @@ -791,7 +791,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, - 0.25 * (h_x * v + h * v_x) * b_x * v) if equations.bathymetry_type isa BathymetryVariable - @. psi = psi - 0.125 * (h_x * v + h * v_x) * b_x * v + @. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v end @. p_h = p_h + tmp @. tmp = tmp * h From 3e633a5f412b151ec231e41e233ce9ad2b29d3f3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:38:31 +0200 Subject: [PATCH 19/34] format --- src/equations/serre_green_naghdi_1d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 96d62e61..acb0f22e 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -335,8 +335,8 @@ function create_cache(mesh, system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + D1mat' * (Diagonal(M_h3_3) * D1mat - - - Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx)) - Diagonal(M_h2_bx) * D1mat) @@ -697,8 +697,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) + D1mat' * (Diagonal(M_h3_3) * D1mat - - - Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx)) - Diagonal(M_h2_bx) * D1mat) end From 6008569a87ae2eee0ced568d9f4eb4b859d3be4e Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 06:48:36 +0200 Subject: [PATCH 20/34] update Dingemans test values --- test/test_serre_green_naghdi_1d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index 2b2099c6..a508c5dd 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -217,7 +217,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") tspan=(0.0, 1.0), l2=[0.2463295492331802, 0.805357513755897, 0.0], linf=[0.036265772921736605, 0.11763218152350845, 0.0], - cons_error=[5.684341886080802e-14, 3.509473432346808e-5, 0.0], + cons_error=[5.684341886080802e-14, 3.509473432334318e-5, 0.0], change_waterheight=5.684341886080802e-14, change_entropy=1.9489684518703143e-5, change_entropy_modified=-1.2177679309388623e-8) @@ -232,7 +232,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") bathymetry_type=bathymetry_mild_slope, l2=[0.2463295492331802, 0.805357513755897, 0.0], linf=[0.036265772921736605, 0.11763218152350845, 0.0], - cons_error=[5.684341886080802e-14, 3.509473432346624e-5, 0.0], + cons_error=[5.684341886080802e-14, 3.5094734323342747e-5, 0.0], change_waterheight=5.684341886080802e-14, change_entropy=1.9489684518703143e-5, change_entropy_modified=-1.2177679309388623e-8) From 5a4dbb2265b568b459e21df11cd7128e61045391 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 10:18:27 +0200 Subject: [PATCH 21/34] fix bug --- src/equations/serre_green_naghdi_1d.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index acb0f22e..cb40dadf 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -265,6 +265,18 @@ function create_cache(mesh, M_h3_3 = zero(h) M_h2_bx = zero(h) + # b_x appears in the system matrix, so it must not be zero to get the + # correct pattern of the factorization. Just setting it to unity does + # also not work since some terms cancel. Thus, b_x needs to be + # "variable enough" to get the correct pattern of non-zero entries in + # the system matrix and the factorization. + # TODO: This is a hack and should be improved. It would be ideal if we + # had access to the initial condition and the exact value of the + # bathymetry here. + let x = grid(D1) + @. b_x = x^3 + end + if D1 isa PeriodicUpwindOperators v_x_upwind = zero(h) p_0 = zero(h) From d5bb6e690ce7226775ca331075a52a4e3d6d809b Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 10:35:12 +0200 Subject: [PATCH 22/34] add new test used for benchmarks --- .../serre_green_naghdi_conservation.jl | 69 +++++++++++++++++++ test/test_serre_green_naghdi_1d.jl | 25 +++++++ 2 files changed, 94 insertions(+) create mode 100644 examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl new file mode 100644 index 00000000..2cd858e8 --- /dev/null +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl @@ -0,0 +1,69 @@ +# This elixir contains an aritficial setup that can be used to check the +# conservation properties of the equations and numerical methods as well as +# a possible directional bias (if the velocity is set to zero). See +# - Hendrik Ranocha and Mario Ricchiuto (2024) +# Structure-preserving approximations of the Serre-Green-Naghdi +# equations in standard and hyperbolic form +# [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) + +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: upwind_operators, periodic_derivative_operator + +############################################################################### +# Semidiscretization of the Svärd-Kalisch equations + +bathymetry_type = bathymetry_variable # or bathymetry_mild_slope +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81, + eta0 = 1.0) + +function initial_condition_conservation_test(x, t, + equations::SerreGreenNaghdiEquations1D, + mesh) + eta = 1 + exp(-x^2) + v = 1.0e-2 # set this to zero to test a directional bias + b = 0.25 * cospi(x / 75) + + D = equations.eta0 - b + return SVector(eta, v, D) +end + +# create homogeneous mesh +coordinates_min = -150.0 +coordinates_max = +150.0 +N = 1_000 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 2 +accuracy_order = 2 +solver = Solver(mesh, accuracy_order) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, + initial_condition_conservation_test, solver; + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 35.0) +ode = semidiscretize(semi, tspan) + +# The callbacks support an additional `io` argument to write output to a file +# or any other IO stream. The default is stdout. We use this here to enable +# setting it to `devnull` to benchmark the full simulation including the time +# to compute the errors etc. but without the time to write the output to the +# terminal. +io = stdout +summary_callback = SummaryCallback(io) +analysis_callback = AnalysisCallback(semi; interval = 10, io, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + momentum, + entropy, + entropy_modified)) + +callbacks = CallbackSet(analysis_callback, summary_callback) + +sol = solve(ode, Tsit5(); + save_everystep = false, callback = callbacks) diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index a508c5dd..1ca2b549 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -239,6 +239,31 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=750_000) end + + @trixi_testset "serre_green_naghdi_conservation.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_conservation.jl"), + l2=[1.3655498085989206, 2.3967486930606716, 0.0], + linf=[1.001076318001934, 0.8052527556023067, 0.0], + cons_error=[0.0, 0.0002674927404067162, 0.0], + change_entropy=-0.0584189861183404, + change_entropy_modified=0.059273537492344985) + + @test_allocations(semi, sol, allocs=900_000) + end + + @trixi_testset "serre_green_naghdi_conservation.jl with bathymetry_mild_slope" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_conservation.jl"), + bathymetry_type=bathymetry_mild_slope, + l2=[1.3655493671985637, 2.3967828251339003, 0.0], + linf=[1.001075913983051, 0.8052680970114169, 0.0], + cons_error=[1.1368683772161603e-13, 0.00026407261543415217, 0.0], + change_entropy=-0.058352284294869605, + change_entropy_modified=0.05927339747017868) + + @test_allocations(semi, sol, allocs=900_000) + end end end # module From d27cfd7287ea6e3a3cadee87d68f871c890b00f4 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 10:36:10 +0200 Subject: [PATCH 23/34] fix typo --- .../serre_green_naghdi_1d/serre_green_naghdi_conservation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl index 2cd858e8..51dd0897 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl @@ -1,4 +1,4 @@ -# This elixir contains an aritficial setup that can be used to check the +# This elixir contains an artificial setup that can be used to check the # conservation properties of the equations and numerical methods as well as # a possible directional bias (if the velocity is set to zero). See # - Hendrik Ranocha and Mario Ricchiuto (2024) From 3bfc324b7b1e6625890974712b4e4a81eb3dfcb4 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 10:52:35 +0200 Subject: [PATCH 24/34] Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- .../serre_green_naghdi_1d/serre_green_naghdi_conservation.jl | 2 +- examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl index 51dd0897..709f3958 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl @@ -11,7 +11,7 @@ using DispersiveShallowWater using SummationByPartsOperators: upwind_operators, periodic_derivative_operator ############################################################################### -# Semidiscretization of the Svärd-Kalisch equations +# Semidiscretization of the Serre-Green-Naghdi equations bathymetry_type = bathymetry_variable # or bathymetry_mild_slope equations = SerreGreenNaghdiEquations1D(bathymetry_type; diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl index f424fa4c..0596f421 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl @@ -3,7 +3,7 @@ using DispersiveShallowWater using SummationByPartsOperators: upwind_operators, periodic_derivative_operator ############################################################################### -# Semidiscretization of the Svärd-Kalisch equations +# Semidiscretization of the Serre-Green-Naghdi equations bathymetry_type = bathymetry_variable # or bathymetry_mild_slope equations = SerreGreenNaghdiEquations1D(bathymetry_type; From e6734145a790251b4cab102b48108292d72ca15f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 10:53:53 +0200 Subject: [PATCH 25/34] Update examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- .../serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl index 38e21c67..aac19ffb 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl @@ -3,7 +3,7 @@ using DispersiveShallowWater using SummationByPartsOperators: upwind_operators, periodic_derivative_operator ############################################################################### -# Semidiscretization of the BBM-BBM equations +# Semidiscretization of the Serre-Green-Naghdi equations bathymetry_type = bathymetry_variable # or bathymetry_mild_slope equations = SerreGreenNaghdiEquations1D(bathymetry_type; From c7ff29a129e56ec43da45bb35944eb98caa9473c Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 11:22:56 +0200 Subject: [PATCH 26/34] assemble_system_matrix! --- src/equations/serre_green_naghdi_1d.jl | 233 +++++++++---------------- 1 file changed, 82 insertions(+), 151 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index cb40dadf..c56ea60d 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -186,17 +186,10 @@ function create_cache(mesh, D1mat_minus = sparse(D1.minus) - @. M_h = h - scale_by_mass_matrix!(M_h, D1) - @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - system_matrix = Symmetric(Diagonal(M_h) - + - D1mat_minus' * Diagonal(M_h3_3) * D1mat_minus) + system_matrix = let cache = (; M_h, M_h3_3) + assemble_system_matrix!(cache, h, + D1, D1mat_minus, equations) + end factorization = cholesky(system_matrix) cache = (; h, b, h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, @@ -214,18 +207,10 @@ function create_cache(mesh, else D1mat = sparse(D1) - @. M_h = h - scale_by_mass_matrix!(M_h, D1) - @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - system_matrix = Symmetric(Diagonal(M_h) - + - D1mat' * Diagonal(M_h3_3) * D1mat) - + system_matrix = let cache = (; M_h, M_h3_3) + assemble_system_matrix!(cache, h, + D1, D1mat, equations) + end factorization = cholesky(system_matrix) cache = (; h, b, h_x, v_x, h2_x, hv_x, v2_x, @@ -273,6 +258,7 @@ function create_cache(mesh, # TODO: This is a hack and should be improved. It would be ideal if we # had access to the initial condition and the exact value of the # bathymetry here. + @show initial_condition let x = grid(D1) @. b_x = x^3 end @@ -283,32 +269,10 @@ function create_cache(mesh, D1mat_minus = sparse(D1.minus) - if equations.bathymetry_type isa BathymetryMildSlope - factor = 0.75 - elseif equations.bathymetry_type isa BathymetryVariable - factor = 1.0 - else - throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry_type)")) + system_matrix = let cache = (; M_h_p_h_bx2, M_h3_3, M_h2_bx) + assemble_system_matrix!(cache, h, b_x, + D1, D1mat_minus, equations) end - @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D1) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + - D1mat_minus' * (Diagonal(M_h3_3) * D1mat_minus - - - Diagonal(M_h2_bx)) - - - Diagonal(M_h2_bx) * D1mat_minus) - factorization = cholesky(system_matrix) cache = (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, @@ -326,32 +290,10 @@ function create_cache(mesh, else D1mat = sparse(D1) - if equations.bathymetry_type isa BathymetryMildSlope - factor = 0.75 - elseif equations.bathymetry_type isa BathymetryVariable - factor = 1.0 - else - throw(ArgumentError("Unsupported bathymetry type $(equations.bathymetry_type)")) + system_matrix = let cache = (; M_h_p_h_bx2, M_h3_3, M_h2_bx) + assemble_system_matrix!(cache, h, b_x, + D1, D1mat, equations) end - @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D1) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + - D1mat' * (Diagonal(M_h3_3) * D1mat - - - Diagonal(M_h2_bx)) - - - Diagonal(M_h2_bx) * D1mat) - factorization = cholesky(system_matrix) cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, @@ -369,6 +311,49 @@ function create_cache(mesh, return cache end +function assemble_system_matrix!(cache, h, D1, D1mat, + equations::SerreGreenNaghdiEquations1D{BathymetryFlat}) + (; M_h, M_h3_3) = cache + + @. M_h = h + scale_by_mass_matrix!(M_h, D1) + @. M_h3_3 = (1 / 3) * h^3 + scale_by_mass_matrix!(M_h3_3, D1) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + return Symmetric(Diagonal(M_h) + D1mat' * Diagonal(M_h3_3) * D1mat) +end + +# variable bathymetry +function assemble_system_matrix!(cache, h, b_x, D1, D1mat, + equations::SerreGreenNaghdiEquations1D) + (; M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + + if equations.bathymetry_type isa BathymetryMildSlope + factor = 0.75 + elseif equations.bathymetry_type isa BathymetryVariable + factor = 1.0 + end + @. M_h_p_h_bx2 = h + factor * h * b_x^2 + scale_by_mass_matrix!(M_h_p_h_bx2, D1) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D1) + @. M_h2_bx = 0.5 * h^2 * b_x + scale_by_mass_matrix!(M_h2_bx, D1) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + return Symmetric(Diagonal(M_h_p_h_bx2) + + D1mat' * (Diagonal(M_h3_3) * D1mat + - Diagonal(M_h2_bx)) + - Diagonal(M_h2_bx) * D1mat) +end + # Discretization that conserves # - the total water mass (integral of h) as a linear invariant # - the total momentum (integral of h v) as a nonlinear invariant @@ -457,21 +442,13 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla p_x) end + # The code below is equivalent to + # dv .= (Diagonal(h) - D1mat * Diagonal(1/3 .* h.^3) * D1mat) \ tmp + # but faster since the symbolic factorization is reused. @trixi_timeit timer() "assembling elliptic operator" begin - # The code below is equivalent to - # dv .= (Diagonal(h) - D1mat * Diagonal(1/3 .* h.^3) * D1mat) \ tmp - # but faster since the symbolic factorization is reused. - @. M_h = h - scale_by_mass_matrix!(M_h, D1) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. - system_matrix = Symmetric(Diagonal(M_h) - + - D1mat' * Diagonal(M_h3_3) * D1mat) + system_matrix = assemble_system_matrix!(cache, h, + D1, D1mat, + equations) end @trixi_timeit timer() "solving elliptic system" begin @@ -565,21 +542,13 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat p_x) end + # The code below is equivalent to + # dv .= (Diagonal(h) - D1mat_plus * Diagonal(1/3 .* h.^3) * D1mat_minus) \ tmp + # but faster since the symbolic factorization is reused. @trixi_timeit timer() "assembling elliptic operator" begin - # The code below is equivalent to - # dv .= (Diagonal(h) - D1mat_plus * Diagonal(1/3 .* h.^3) * D1mat_minus) \ tmp - # but faster since the symbolic factorization is reused. - @. M_h = h - scale_by_mass_matrix!(M_h, D1) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off errors. - # We wrap it here to avoid issues with the factorization. - system_matrix = Symmetric(Diagonal(M_h) - + - D1mat_minus' * Diagonal(M_h3_3) * D1mat_minus) + system_matrix = assemble_system_matrix!(cache, h, + D1, D1mat_minus, + equations) end @trixi_timeit timer() "solving elliptic system" begin @@ -687,32 +656,13 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end end + # The code below is equivalent to + # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp + # but faster since the symbolic factorization is reused. @trixi_timeit timer() "assembling elliptic operator" begin - if equations.bathymetry_type isa BathymetryMildSlope - factor = 0.75 - else # equations.bathymetry_type isa BathymetryVariable - factor = 1.0 - end - # The code below is equivalent to - # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp - # but faster since the symbolic factorization is reused. - @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D1) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. - system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + - D1mat' * (Diagonal(M_h3_3) * D1mat - - - Diagonal(M_h2_bx)) - - - Diagonal(M_h2_bx) * D1mat) + system_matrix = assemble_system_matrix!(cache, h, b_x, + D1, D1mat, + equations) end @trixi_timeit timer() "solving elliptic system" begin @@ -834,32 +784,13 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end end + # The code below is equivalent to + # dv .= (Diagonal(h .+ 0.75 .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp + # but faster since the symbolic factorization is reused. @trixi_timeit timer() "assembling elliptic operator" begin - if equations.bathymetry_type isa BathymetryMildSlope - factor = 0.75 - else # equations.bathymetry_type isa BathymetryVariable - factor = 1.0 - end - # The code below is equivalent to - # dv .= (Diagonal(h .+ 0.75 .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp - # but faster since the symbolic factorization is reused. - @. M_h_p_h_bx2 = h + factor * h * b_x^2 - scale_by_mass_matrix!(M_h_p_h_bx2, D1) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D1) - @. M_h2_bx = 0.5 * h^2 * b_x - scale_by_mass_matrix!(M_h2_bx, D1) - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. - system_matrix = Symmetric(Diagonal(M_h_p_h_bx2) - + - D1mat_minus' * (Diagonal(M_h3_3) * D1mat_minus - - - Diagonal(M_h2_bx)) - - - Diagonal(M_h2_bx) * D1mat_minus) + system_matrix = assemble_system_matrix!(cache, h, b_x, + D1, D1mat_minus, + equations) end @trixi_timeit timer() "solving elliptic system" begin From 0062720c18c8e60a7e1970f333fc5ec52f373497 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 11:27:35 +0200 Subject: [PATCH 27/34] solve_system_matrix! --- src/equations/serre_green_naghdi_1d.jl | 88 +++++++++----------------- 1 file changed, 30 insertions(+), 58 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index c56ea60d..931c719d 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -354,6 +354,28 @@ function assemble_system_matrix!(cache, h, b_x, D1, D1mat, - Diagonal(M_h2_bx) * D1mat) end +function solve_system_matrix!(dv, system_matrix, rhs, + ::SerreGreenNaghdiEquations1D, + D1, cache) + if issparse(system_matrix) + (; factorization) = cache + cholesky!(factorization, system_matrix; check = false) + if issuccess(factorization) + scale_by_mass_matrix!(rhs, D1) + dv .= factorization \ rhs + else + # The factorization may fail if the time step is too large + # and h becomes negative. + fill!(dv, NaN) + end + else + factorization = cholesky!(system_matrix) + scale_by_mass_matrix!(rhs, D1) + ldiv!(dv, factorization, rhs) + end + return nothing +end + # Discretization that conserves # - the total water mass (integral of h) as a linear invariant # - the total momentum (integral of h v) as a nonlinear invariant @@ -452,22 +474,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla end @trixi_timeit timer() "solving elliptic system" begin - if issparse(system_matrix) - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D1) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) - end - else - factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D1) - ldiv!(dv, factorization, tmp) - end + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) end return nothing @@ -552,16 +560,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat end @trixi_timeit timer() "solving elliptic system" begin - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D1) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) - end + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) end return nothing @@ -666,22 +666,8 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - if issparse(system_matrix) - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D1) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) - end - else - factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D1) - ldiv!(dv, factorization, tmp) - end + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) end return nothing @@ -794,22 +780,8 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end @trixi_timeit timer() "solving elliptic system" begin - if issparse(system_matrix) - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D1) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) - end - else - factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D1) - ldiv!(dv, factorization, tmp) - end + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) end return nothing From 5b2bc46a81e3355daad44c38bf2566ec3e772962 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 11:30:33 +0200 Subject: [PATCH 28/34] Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 931c719d..fcd13807 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -349,9 +349,12 @@ function assemble_system_matrix!(cache, h, b_x, D1, D1mat, # round-off errors. We wrap it here to avoid issues with the # factorization. return Symmetric(Diagonal(M_h_p_h_bx2) - + D1mat' * (Diagonal(M_h3_3) * D1mat - - Diagonal(M_h2_bx)) - - Diagonal(M_h2_bx) * D1mat) + + + D1mat' * (Diagonal(M_h3_3) * D1mat + - + Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx) * D1mat) end function solve_system_matrix!(dv, system_matrix, rhs, From a0cdae432ab5ffaadea2ca854983bc76e4e0f102 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 11:37:04 +0200 Subject: [PATCH 29/34] update Dingemans --- .../serre_green_naghdi_dingemans.jl | 4 ++-- test/test_serre_green_naghdi_1d.jl | 20 +++++++++---------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl index 0596f421..ab4e0692 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl @@ -13,8 +13,8 @@ initial_condition = initial_condition_dingemans boundary_conditions = boundary_condition_periodic # create homogeneous mesh -coordinates_min = -140.0 -coordinates_max = 100.0 +coordinates_min = -138.0 +coordinates_max = 46.0 N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index 1ca2b549..f231ad7f 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -215,12 +215,12 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_dingemans.jl"), tspan=(0.0, 1.0), - l2=[0.2463295492331802, 0.805357513755897, 0.0], - linf=[0.036265772921736605, 0.11763218152350845, 0.0], - cons_error=[5.684341886080802e-14, 3.509473432334318e-5, 0.0], + l2=[0.24638385233253557, 0.8055038575889898, 0.0], + linf=[0.03627087716483823, 0.11899056073196576, 0.0], + cons_error=[5.684341886080802e-14, 3.57950443142954e-5, 0.0], change_waterheight=5.684341886080802e-14, - change_entropy=1.9489684518703143e-5, - change_entropy_modified=-1.2177679309388623e-8) + change_entropy=2.4420619411102962e-5, + change_entropy_modified=-1.0404164640931413e-8) @test_allocations(semi, sol, allocs=750_000) end @@ -230,12 +230,12 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") "serre_green_naghdi_dingemans.jl"), tspan=(0.0, 1.0), bathymetry_type=bathymetry_mild_slope, - l2=[0.2463295492331802, 0.805357513755897, 0.0], - linf=[0.036265772921736605, 0.11763218152350845, 0.0], - cons_error=[5.684341886080802e-14, 3.5094734323342747e-5, 0.0], + l2=[0.24638385233253557, 0.8055038575889898, 0.0], + linf=[0.03627087716483823, 0.11899056073196576, 0.0], + cons_error=[5.684341886080802e-14, 3.579504431429478e-5, 0.0], change_waterheight=5.684341886080802e-14, - change_entropy=1.9489684518703143e-5, - change_entropy_modified=-1.2177679309388623e-8) + change_entropy=2.4420619411102962e-5, + change_entropy_modified=-1.0404164640931413e-8) @test_allocations(semi, sol, allocs=750_000) end From a898e7ed1695434ca4c2e48b69da78ae2a1a6106 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 13:05:59 +0200 Subject: [PATCH 30/34] Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index fcd13807..578a7ddc 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -258,7 +258,6 @@ function create_cache(mesh, # TODO: This is a hack and should be improved. It would be ideal if we # had access to the initial condition and the exact value of the # bathymetry here. - @show initial_condition let x = grid(D1) @. b_x = x^3 end From 7e28939d2f1b43d23d99b18ed22925879d942ba0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 13:07:12 +0200 Subject: [PATCH 31/34] Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/equations/serre_green_naghdi_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 578a7ddc..a5d34a81 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -773,7 +773,7 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, end # The code below is equivalent to - # dv .= (Diagonal(h .+ 0.75 .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp + # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp # but faster since the symbolic factorization is reused. @trixi_timeit timer() "assembling elliptic operator" begin system_matrix = assemble_system_matrix!(cache, h, b_x, From 59f77296741f4e0d103aecc9f9eb40e0761ea5f4 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 14:55:50 +0200 Subject: [PATCH 32/34] fix docstring of energy_total_modified --- src/equations/serre_green_naghdi_1d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index a5d34a81..a3387ab8 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -831,15 +831,15 @@ derivative compared to the usual [`energy_total`](@ref) modeling non-hydrostatic contributions. The [`energy_total_modified`](@ref) is a conserved quantity (for periodic boundary conditions). -For a [`bathymetry_flat`](@ref) the total energy is given by +For a [`bathymetry_flat`](@ref) the total modified energy is given by ```math \\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2. ``` -For a [`bathymetry_mild_slope`](@ref) the total energy is given by +For a [`bathymetry_mild_slope`](@ref) the total modified energy is given by ```math \\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h (-h v_x + 1.5 v b_x)^2. ``` -For a [`bathymetry_variable`](@ref) the total energy has the additional term +For a [`bathymetry_variable`](@ref) the total modified energy has the additional term ```math + \\frac{1}{8} h (v b_x)^2. ``` From 0dfab87b95e88f5904f11b8fe41890d78effa3b8 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 20:06:42 +0200 Subject: [PATCH 33/34] adapt tolerances for macOS CI --- test/test_serre_green_naghdi_1d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index f231ad7f..71efeab7 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -247,7 +247,8 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") linf=[1.001076318001934, 0.8052527556023067, 0.0], cons_error=[0.0, 0.0002674927404067162, 0.0], change_entropy=-0.0584189861183404, - change_entropy_modified=0.059273537492344985) + change_entropy_modified=0.059273537492344985, + atol_ints=4e-9) # to make CI pass @test_allocations(semi, sol, allocs=900_000) end From d4f9dd6d9c891859bf7ae937d14b141ddd8eee7e Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 20:09:23 +0200 Subject: [PATCH 34/34] not only integrals --- test/test_serre_green_naghdi_1d.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index 71efeab7..d0679264 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -248,6 +248,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") cons_error=[0.0, 0.0002674927404067162, 0.0], change_entropy=-0.0584189861183404, change_entropy_modified=0.059273537492344985, + atol=1e-11, # to make CI pass atol_ints=4e-9) # to make CI pass @test_allocations(semi, sol, allocs=900_000)