10507790104

Committed 22 Aug 2024 12:10PM UTC coverage: 97.925%. First build

Build # 10507790104

Build Type

Pull #146

github

Committed by

web-flow

Commit Message

Merge 92d25e5f3 into cf20e3f45

Pull Request Pull Request #146: Adapt Svärd-Kalisch equations to Serre-Green-Naghdi equations

Run Details

152 of 155 new or added lines in 7 files covered. (98.06%)

1746 of 1783 relevant lines covered (97.92%)

3243389.41 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

99.02

/src/equations/svaerd_kalisch_1d.jl

@doc raw"""
    SvaerdKalischEquations1D(bathymetry_type = bathymetry_variable;
                             gravity_constant, eta0 = 0.0, alpha = 0.0,
                             beta = 0.2308939393939394, gamma = 0.04034343434343434)

Dispersive system by Svärd and Kalisch in one spatial dimension. The equations for variable bathymetry
are given in conservative variables by
```math
\begin{aligned}
  h_t + (hv)_x &= (\hat\alpha(\hat\alpha(h + b)_x)_x)_x,\\
  (hv)_t + (hv^2)_x + gh(h + b)_x &= (\hat\alpha v(\hat\alpha(h + b)_x)_x)_x + (\hat\beta v_x)_{xt} + \frac{1}{2}(\hat\gamma v_x)_{xx} + \frac{1}{2}(\hat\gamma v_{xx})_x,
\end{aligned}
```
where ``\hat\alpha^2 = \alpha\sqrt{gD}D^2``, ``\hat\beta = \beta D^3``, ``\hat\gamma = \gamma\sqrt{gD}D^3``.
The coefficients ``\alpha``, ``\beta`` and ``\gamma`` are provided in dimensionless form and ``D = \eta_0 - b`` is the
still-water depth and `eta0` is the still-water surface (lake-at-rest).
The equations can be rewritten in primitive variables as
```math
\begin{aligned}
  \eta_t + ((\eta + D)v)_x &= (\hat\alpha(\hat\alpha\eta_x)_x)_x,\\
  v_t(\eta + D) - v((\eta + D)v)_x + ((\eta + D)v^2)_x + g(\eta + D)\eta_x &= (\hat\alpha v(\hat\alpha\eta_x)_x)_x - v(\hat\alpha(\hat\alpha\eta_x)_x)_x + (\hat\beta v_x)_{xt} + \frac{1}{2}(\hat\gamma v_x)_{xx} + \frac{1}{2}(\hat\gamma v_{xx})_x.
\end{aligned}
```
The unknown quantities of the Svärd-Kalisch equations are the total water height ``\eta`` 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``.

Currently, the equations only support a general variable bathymetry, see [`bathymetry_variable`](@ref).

`SvärdKalischEquations1D` is an alias for `SvaerdKalischEquations1D`.

The equations by Svärd and Kalisch are presented and analyzed in Svärd and Kalisch (2023).
The semidiscretization implemented here conserves
- the total water (integral of ``h``) as a linear invariant
- the total momentum (integral of ``h v``) as a nonlinear invariant for flat bathymetry
- the total modified energy

for periodic boundary conditions (see Lampert, Ranocha).
Additionally, it is well-balanced for the lake-at-rest stationary solution, see Lampert and Ranocha (2024).

- 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)
- Joshua Lampert, Hendrik Ranocha (2024)
  Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
  [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
"""
struct SvaerdKalischEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <:
       AbstractSvaerdKalischEquations{1, 3}
    bathymetry_type::Bathymetry # type of bathymetry
    gravity::RealT # gravitational constant
    eta0::RealT    # constant still-water surface
    alpha::RealT   # coefficient
    beta::RealT    # coefficient
    gamma::RealT   # coefficient
end

const SvärdKalischEquations1D = SvaerdKalischEquations1D

function SvaerdKalischEquations1D(bathymetry_type = bathymetry_variable;
                                  gravity_constant, eta0 = 0.0, alpha = 0.0,
                                  beta = 0.2308939393939394, gamma = 0.04034343434343434)
    SvaerdKalischEquations1D(bathymetry_type, gravity_constant, eta0, alpha, beta, gamma)
end

varnames(::typeof(prim2prim), ::SvaerdKalischEquations1D) = ("η", "v", "D")
varnames(::typeof(prim2cons), ::SvaerdKalischEquations1D) = ("h", "hv", "b")

"""
    initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, 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::SvaerdKalischEquations1D, mesh)
    g = gravity_constant(equations)
    h0 = 0.8
    A = 0.02
    # omega = 2*pi/(2.02*sqrt(2))
    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - g * k * tanh(k * h0), 1.0) using Roots.jl
    if x < -30.5 * pi / k || x > -8.5 * pi / k
        h = 0.0
    else
        h = A * cos(k * x)
    end
    v = sqrt(g / 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

"""
    initial_condition_manufactured(x, t, equations::SvaerdKalischEquations1D, mesh)

A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref).
"""
function initial_condition_manufactured(x, t,
                                        equations::SvaerdKalischEquations1D,
                                        mesh)
    eta = exp(t) * cospi(2 * (x - 2 * t))
    v = exp(t / 2) * sinpi(2 * (x - t / 2))
    b = -5 - 2 * cospi(2 * x)
    D = equations.eta0 - b
    return SVector(eta, v, D)
end

"""
    source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D, mesh)

A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref).
"""
function source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D)
    g = gravity_constant(equations)
    eta0 = still_water_surface(q, equations)
    alpha = equations.alpha
    beta = equations.beta
    gamma = equations.gamma
    a1 = sinpi(2 * x)
    a2 = cospi(2 * x)
    a3 = sinpi(-t + 2 * x)
    a4 = cospi(-t + 2 * x)
    a5 = sinpi(t - 2 * x)
    a6 = cospi(t - 2 * x)
    a7 = sinpi(-4 * t + 2 * x)
    a8 = exp(t / 2)
    a9 = exp(t)
    a10 = a9 * cospi(-4 * t + 2 * x)
    a11 = eta0 + 2.0 * a2 + 5.0
    a12 = sqrt(g * a11)
    a13 = 0.2 * eta0 + 0.4 * a2 + 1
    a14 = alpha * a12 * a13^2
    a15 = sqrt(a14)
    a16 = -1.0 * pi * a14 * a1 / a11 - 0.8 * pi * alpha * a12 * a13 * a1
    a17 = -20.0 * pi^2 * a15 * a10 - 10.0 * pi * a15 * a16 * a9 * a7 / (a14)
    a18 = -2 * pi * a9 * a7 - 4.0 * pi * a1
    a19 = a10 + 2.0 * a2 + 5.0
    a20 = a18 * a8 * a3 + 2 * pi * a19 * a8 * a4
    a21 = a15 * (40.0 * pi^3 * a15 * a9 * a7 - 40.0 * pi^2 * a15 * a16 * a10 / (a14) -
           20.0 * pi^2 * a15 * a16 * a9 * a1 * a7 / (a14 * a11) -
           16.0 * pi^2 * a15 * a16 * a9 * a1 * a7 / (alpha * a12 * a13^3) -
           10.0 * pi * a15 *
           (-2.0 * pi^2 * a14 * a2 / a11 - 1.6 * pi^2 * alpha * a12 * a13 * a2 +
            3.2 * pi^2 * alpha * a12 * a13 * a1^2 / a11 + 0.56 * pi^2 * alpha * a12 * a1^2) *
           a9 * a7 / (a14) -
           10.0 * pi * a15 * a16^2 * a9 * a7 / (alpha^2 * g * a13^4 * a11))

    dq1 = -5.0 * a21 + a20 + 4 * pi * a9 * a7 + a10 - 5.0 * a15 * a17 * a16 / (a14)

    dq2 = -25.0 * beta * (-2 * pi^2 * a8 * a3 + 4 * pi^3 * a8 * a4) * a13^2 * a11 +
          100.0 * pi * beta * (2 * pi^2 * a8 * a3 + pi * a8 * a4) * a13^2 * a1 +
          40.0 * pi * beta * (2 * pi^2 * a8 * a3 + pi * a8 * a4) * a13 * a11 * a1 -
          2 * pi * g * a19 * a9 * a7 + 100.0 * pi^3 * gamma * a12 * a13^2 * a11 * a8 * a4 -
          300.0 * pi^3 * gamma * a12 * a13^2 * a8 * a1 * a3 -
          80.0 * pi^3 * gamma * a12 * a13 * a11 * a8 * a1 * a3 -
          pi^3 * gamma * a12 *
          (-50.0 * (3.2 * a13 * a2 - 1.28 * a1^2) * a11 * a6 -
           50.0 * (4.0 * a2 / a11 + 0.16 * a1^2 / a13^2) * a13^2 * a11 * a6 -
           200.0 * a13^2 * a11 * a6 - 1200.0 * a13^2 * a1 * a5 - 400.0 * a13^2 * a2 * a6 +
           800.0 * a13^2 * a1^2 * a6 / a11 - 320.0 * a13 * a11 * a1 * a5 +
           960.0 * a13 * a1^2 * a6) * a8 / 2 - 10.0 * pi * a15 * a17 * a8 * a4 -
          2.5 * a21 * a8 * a3 + (5.0 * a21 + 5.0 * a15 * a17 * a16 / (a14)) * a8 * a3 / 2 +
          (a8 * a3 / 2 - pi * a8 * a4) * a19 + a18 * a9 * a3^2 / 2 - (a20) * a8 * a3 / 2 +
          3 * pi * a19 * a9 * a3 * a4 - 2.5 * a15 * a17 * a16 * a8 * a3 / (a14)

    return SVector(dq1, dq2, zero(dq1))
end

function create_cache(mesh, equations::SvaerdKalischEquations1D,
                      solver, initial_condition,
                      ::BoundaryConditionPeriodic,
                      RealT, uEltype)
    D1 = solver.D1
    #  Assume D is independent of time and compute D evaluated at mesh points once.
    D = Array{RealT}(undef, nnodes(mesh))
    x = grid(solver)
    for i in eachnode(solver)
        D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations)
    end
    g = gravity_constant(equations)
    h = ones(RealT, nnodes(mesh))
    hv = zero(h)
    b = zero(h)
    eta_x = zero(h)
    v_x = zero(h)
    alpha_eta_x_x = zero(h)
    y_x = zero(h)
    v_y_x = zero(h)
    yv_x = zero(h)
    vy = zero(h)
    vy_x = zero(h)
    y_v_x = zero(h)
    h_v_x = zero(h)
    hv2_x = zero(h)
    v_xx = zero(h)
    gamma_v_xx_x = zero(h)
    gamma_v_x_xx = zero(h)
    alpha_hat = sqrt.(equations.alpha * sqrt.(g * D) .* D .^ 2)
    beta_hat = equations.beta * D .^ 3
    gamma_hat = equations.gamma * sqrt.(g * D) .* D .^ 3
    tmp2 = zero(h)
    M = mass_matrix(D1)
    M_h = zero(h)
    M_beta = copy(beta_hat)
    scale_by_mass_matrix!(M_beta, D1)
    if D1 isa PeriodicDerivativeOperator ||
       D1 isa UniformPeriodicCoupledOperator
        D1_central = D1
        D1mat = sparse(D1_central)
        minus_MD1betaD1 = D1mat' * Diagonal(M_beta) * D1mat
        system_matrix = let cache = (; M_h, minus_MD1betaD1)
            assemble_system_matrix!(cache, h,
                                    D1, D1mat, equations)
        end
    elseif D1 isa PeriodicUpwindOperators
        D1_central = D1.central
        D1mat_minus = sparse(D1.minus)
        minus_MD1betaD1 = D1mat_minus' * Diagonal(M_beta) * D1mat_minus
        system_matrix = let cache = (; M_h, minus_MD1betaD1)
            assemble_system_matrix!(cache, h,
                                    D1, D1mat_minus, equations)
        end
    else
        @error "unknown type of first-derivative operator: $(typeof(D1))"
    end
    factorization = cholesky(system_matrix)
    cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x,
             alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx,
             gamma_v_xx_x, gamma_v_x_xx,
             alpha_hat, beta_hat, gamma_hat, tmp2, D1_central, M, D1, M_h)
    if D1 isa PeriodicUpwindOperators
        eta_x_upwind = zero(h)
        v_x_upwind = zero(h)
        cache = (; cache..., eta_x_upwind, v_x_upwind)
    end
    return cache
end

function assemble_system_matrix!(cache, h, D1, D1mat,
                                 ::SvaerdKalischEquations1D)
    (; M_h, minus_MD1betaD1) = cache

    @.. M_h = h
    scale_by_mass_matrix!(M_h, 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) + minus_MD1betaD1)
end

# Discretization that conserves
# - the total water (integral of h) as a linear invariant
# - the total momentum (integral of h v) as a nonlinear invariant for flat bathymetry
# - the total modified energy
# for periodic boundary conditions, see
# - Joshua Lampert and Hendrik Ranocha (2024)
#   Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
#   [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
# TODO: Simplify for the case of flat bathymetry and use higher-order operators
function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,
              initial_condition, ::BoundaryConditionPeriodic, source_terms,
              solver, cache)
    (; D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x,
    y_v_x, h_v_x, hv2_x, v_xx, gamma_v_xx_x, gamma_v_x_xx, alpha_hat, gamma_hat,
    tmp1, tmp2, D1_central, M, D1) = cache

    g = gravity_constant(equations)
    eta, v = q.x
    deta, dv, dD = dq.x
    fill!(dD, zero(eltype(dD)))

    @trixi_timeit timer() "deta hyperbolic" begin
        @.. b = equations.eta0 - D
        @.. h = eta - b
        @.. hv = h * v

        if D1 isa PeriodicDerivativeOperator ||
           D1 isa UniformPeriodicCoupledOperator
            mul!(eta_x, D1_central, eta)
            mul!(v_x, D1_central, v)
            @.. tmp1 = alpha_hat * eta_x
            mul!(alpha_eta_x_x, D1_central, tmp1)
            @.. tmp1 = alpha_hat * alpha_eta_x_x
            mul!(y_x, D1_central, tmp1)
            @.. v_y_x = v * y_x
            @.. vy = v * tmp1
            mul!(vy_x, D1_central, vy)
            @.. y_v_x = tmp1 * v_x
            @.. tmp2 = tmp1 - hv
            mul!(deta, D1_central, tmp2)
        elseif D1 isa PeriodicUpwindOperators
            @unpack eta_x_upwind, v_x_upwind = cache
            mul!(eta_x, D1_central, eta)
            mul!(v_x, D1_central, v)
            mul!(eta_x_upwind, D1.plus, eta)
            @.. tmp1 = alpha_hat * eta_x_upwind
            mul!(alpha_eta_x_x, D1_central, tmp1)
            @.. tmp1 = alpha_hat * alpha_eta_x_x
            mul!(y_x, D1.minus, tmp1)
            @.. v_y_x = v * y_x
            @.. vy = v * tmp1
            mul!(vy_x, D1.minus, vy)
            mul!(v_x_upwind, D1.plus, v)
            @.. y_v_x = tmp1 * v_x_upwind
            # deta[:] = D1.minus * tmp1 - D1_central * hv
            mul!(deta, D1.minus, tmp1)
            mul!(deta, D1_central, hv, -1.0, 1.0)
        else
            @error "unknown type of first derivative operator: $(typeof(D1))"
        end
    end

    # split form
    @trixi_timeit timer() "dv hyperbolic" begin
        mul!(h_v_x, D1_central, hv)
        @.. tmp1 = hv * v
        mul!(hv2_x, D1_central, tmp1)
        mul!(v_xx, solver.D2, v)
        @.. tmp1 = gamma_hat * v_xx
        mul!(gamma_v_xx_x, D1_central, tmp1)
        @.. tmp1 = gamma_hat * v_x
        mul!(gamma_v_x_xx, solver.D2, tmp1)
        @.. dv = -(0.5 * (hv2_x + hv * v_x - v * h_v_x) +
                   g * h * eta_x +
                   0.5 * (v_y_x - vy_x - y_v_x) -
                   0.5 * gamma_v_xx_x -
                   0.5 * gamma_v_x_xx)
    end

    # no split form
    #     dv[:] = -(D1_central * (hv .* v) - v .* (D1_central * hv)+
    #               g * h .* eta_x +
    #               vy_x - v_y_c -
    #               0.5 * gamma_v_xx_x -
    #               0.5 * gamma_v_x_xx)

    @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
                                                       solver)
    @trixi_timeit timer() "assemble system matrix" begin
        system_matrix = assemble_system_matrix!(cache, h, D1, D1_central, equations)
    end
    @trixi_timeit timer() "solving elliptic system" begin
        tmp1 .= dv
        solve_system_matrix!(dv, system_matrix, tmp1,
                             equations, D1, cache)
    end

    return nothing
end

@inline function prim2cons(q, equations::SvaerdKalischEquations1D)
    eta, v = q

    b = bathymetry(q, equations)
    h = eta - b
    hv = h * v
    return SVector(h, hv, b)
end

@inline function cons2prim(u, equations::SvaerdKalischEquations1D)
    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, ::SvaerdKalischEquations1D)
    return q[1]
end

@inline function velocity(q, ::SvaerdKalischEquations1D)
    return q[2]
end

@inline function bathymetry(q, equations::SvaerdKalischEquations1D)
    D = q[3]
    return equations.eta0 - D
end

# The modified entropy/energy takes the whole `q` for every point in space
"""
    energy_total_modified(q_global, equations::SvaerdKalischEquations1D, cache)

Return the modified total energy of the primitive variables `q_global` for the
[`SvaerdKalischEquations1D`](@ref). It contains an additional term containing a
derivative compared to the usual [`energy_total`](@ref) modeling
non-hydrostatic contributions. The `energy_total_modified`
is a conserved quantity (for periodic boundary conditions).

It is given by
```math
\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{2} \\hat\\beta 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`.
"""
@inline function energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D,
                                        cache)
    # unpack physical parameters and SBP operator `D1`
    g = gravity_constant(equations)
    (; D1, h, b, v_x, beta_hat) = cache

    # `q_global` is an `ArrayPartition`. It collects the individual arrays for
    # the total water height `eta = h + b` and the velocity `v`.
    eta, v, D = q_global.x
    @.. b = equations.eta0 - D
    @.. h = eta - b

    if D1 isa PeriodicUpwindOperators
        mul!(v_x, D1.minus, v)
    else
        mul!(v_x, D1, v)
    end

    @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 2 * beta_hat * v_x^2
    return e
end

1	@doc raw"""
2	SvaerdKalischEquations1D(bathymetry_type = bathymetry_variable;
3	gravity_constant, eta0 = 0.0, alpha = 0.0,
4	beta = 0.2308939393939394, gamma = 0.04034343434343434)
5
6	Dispersive system by Svärd and Kalisch in one spatial dimension. The equations for variable bathymetry
7	are given in conservative variables by
8	```math
9	\begin{aligned}
10	h_t + (hv)_x &= (\hat\alpha(\hat\alpha(h + b)_x)_x)_x,\\
11	(hv)_t + (hv^2)_x + gh(h + b)_x &= (\hat\alpha v(\hat\alpha(h + b)_x)_x)_x + (\hat\beta v_x)_{xt} + \frac{1}{2}(\hat\gamma v_x)_{xx} + \frac{1}{2}(\hat\gamma v_{xx})_x,
12	\end{aligned}
13	```
14	where ``\hat\alpha^2 = \alpha\sqrt{gD}D^2``, ``\hat\beta = \beta D^3``, ``\hat\gamma = \gamma\sqrt{gD}D^3``.
15	The coefficients ``\alpha``, ``\beta`` and ``\gamma`` are provided in dimensionless form and ``D = \eta_0 - b`` is the
16	still-water depth and `eta0` is the still-water surface (lake-at-rest).
17	The equations can be rewritten in primitive variables as
18	```math
19	\begin{aligned}
20	\eta_t + ((\eta + D)v)_x &= (\hat\alpha(\hat\alpha\eta_x)_x)_x,\\
21	v_t(\eta + D) - v((\eta + D)v)_x + ((\eta + D)v^2)_x + g(\eta + D)\eta_x &= (\hat\alpha v(\hat\alpha\eta_x)_x)_x - v(\hat\alpha(\hat\alpha\eta_x)_x)_x + (\hat\beta v_x)_{xt} + \frac{1}{2}(\hat\gamma v_x)_{xx} + \frac{1}{2}(\hat\gamma v_{xx})_x.
22	\end{aligned}
23	```
24	The unknown quantities of the Svärd-Kalisch equations are the total water height ``\eta`` and the velocity ``v``.
25	The gravitational constant is denoted by `g` and the bottom topography (bathymetry) ``b = \eta_0 - D``.
26	The water height above the bathymetry is therefore given by
27	``h = \eta - \eta_0 + D``.
28
29	Currently, the equations only support a general variable bathymetry, see [`bathymetry_variable`](@ref).
30
31	`SvärdKalischEquations1D` is an alias for `SvaerdKalischEquations1D`.
32
33	The equations by Svärd and Kalisch are presented and analyzed in Svärd and Kalisch (2023).
34	The semidiscretization implemented here conserves
35	- the total water (integral of ``h``) as a linear invariant
36	- the total momentum (integral of ``h v``) as a nonlinear invariant for flat bathymetry
37	- the total modified energy
38
39	for periodic boundary conditions (see Lampert, Ranocha).
40	Additionally, it is well-balanced for the lake-at-rest stationary solution, see Lampert and Ranocha (2024).
41
42	- Magnus Svärd, Henrik Kalisch (2023)
43	A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization
44	[arXiv: 2302.09924](https://arxiv.org/abs/2302.09924)
45	- Joshua Lampert, Hendrik Ranocha (2024)
46	Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
47	[DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
48	"""
49	struct SvaerdKalischEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <:
50	AbstractSvaerdKalischEquations{1, 3}
51	bathymetry_type::Bathymetry # type of bathymetry	40✔
52	gravity::RealT # gravitational constant
53	eta0::RealT # constant still-water surface
54	alpha::RealT # coefficient
55	beta::RealT # coefficient
56	gamma::RealT # coefficient
57	end
58
59	const SvärdKalischEquations1D = SvaerdKalischEquations1D
60
61	function SvaerdKalischEquations1D(bathymetry_type = bathymetry_variable;	96✔
62	gravity_constant, eta0 = 0.0, alpha = 0.0,
63	beta = 0.2308939393939394, gamma = 0.04034343434343434)
64	SvaerdKalischEquations1D(bathymetry_type, gravity_constant, eta0, alpha, beta, gamma)	32✔
65	end
66
67	varnames(::typeof(prim2prim), ::SvaerdKalischEquations1D) = ("η", "v", "D")	1,172✔
68	varnames(::typeof(prim2cons), ::SvaerdKalischEquations1D) = ("h", "hv", "b")	4✔
69
70	"""
71	initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, mesh)
72
73	The initial condition that uses the dispersion relation of the Euler equations
74	to approximate waves generated by a wave maker as it is done by experiments of
75	Dingemans. The topography is a trapezoidal. It is assumed that `equations.eta0 = 0.8`.
76
77	References:
78	- Magnus Svärd, Henrik Kalisch (2023)
79	A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization
80	[arXiv: 2302.09924](https://arxiv.org/abs/2302.09924)
81	- Maarten W. Dingemans (1994)
82	Comparison of computations with Boussinesq-like models and laboratory measurements
83	[link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download)
84	"""
85	function initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, mesh)	41,000✔
86	g = gravity_constant(equations)	41,000✔
87	h0 = 0.8	41,000✔
88	A = 0.02	41,000✔
89	# omega = 2pi/(2.02sqrt(2))
90	k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - g * k * tanh(k * h0), 1.0) using Roots.jl	41,000✔
91	if x < -30.5 * pi / k \|\| x > -8.5 * pi / k	76,640✔
92	h = 0.0	22,680✔
93	else
94	h = A * cos(k * x)	18,320✔
95	end
96	v = sqrt(g / k * tanh(k * h0)) * h / h0	41,000✔
97	if 11.01 <= x && x < 23.04	41,000✔
98	b = 0.6 * (x - 11.01) / (23.04 - 11.01)	2,720✔
99	elseif 23.04 <= x && x < 27.04	38,280✔
100	b = 0.6	860✔
101	elseif 27.04 <= x && x < 33.07	37,420✔
102	b = 0.6 * (33.07 - x) / (33.07 - 27.04)	1,380✔
103	else
104	b = 0.0	36,040✔
105	end
106	eta = h + h0	41,000✔
107	D = equations.eta0 - b	41,000✔
108	return SVector(eta, v, D)	41,000✔
109	end
110
111	"""
112	initial_condition_manufactured(x, t, equations::SvaerdKalischEquations1D, mesh)
113
114	A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref).
115	"""
116	function initial_condition_manufactured(x, t,	17,488✔
117	equations::SvaerdKalischEquations1D,
118	mesh)
119	eta = exp(t) * cospi(2 * (x - 2 * t))	17,488✔
120	v = exp(t / 2) * sinpi(2 * (x - t / 2))	17,488✔
121	b = -5 - 2 * cospi(2 * x)	17,488✔
122	D = equations.eta0 - b	17,488✔
123	return SVector(eta, v, D)	17,488✔
124	end
125
126	"""
127	source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D, mesh)
128
129	A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref).
130	"""
131	function source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D)	93,516,032✔
132	g = gravity_constant(equations)	93,516,032✔
133	eta0 = still_water_surface(q, equations)	93,516,032✔
134	alpha = equations.alpha	93,516,032✔
135	beta = equations.beta	93,516,032✔
136	gamma = equations.gamma	93,516,032✔
137	a1 = sinpi(2 * x)	93,516,032✔
138	a2 = cospi(2 * x)	93,516,032✔
139	a3 = sinpi(-t + 2 * x)	93,516,032✔
140	a4 = cospi(-t + 2 * x)	93,516,032✔
141	a5 = sinpi(t - 2 * x)	93,516,032✔
142	a6 = cospi(t - 2 * x)	93,516,032✔
143	a7 = sinpi(-4 * t + 2 * x)	93,516,032✔
144	a8 = exp(t / 2)	93,516,032✔
145	a9 = exp(t)	93,516,032✔
146	a10 = a9 * cospi(-4 * t + 2 * x)	93,516,032✔
147	a11 = eta0 + 2.0 * a2 + 5.0	93,516,032✔
148	a12 = sqrt(g * a11)	93,516,032✔
149	a13 = 0.2 * eta0 + 0.4 * a2 + 1	93,516,032✔
150	a14 = alpha * a12 * a13^2	93,516,032✔
151	a15 = sqrt(a14)	93,516,032✔
152	a16 = -1.0 * pi * a14 * a1 / a11 - 0.8 * pi * alpha * a12 * a13 * a1	93,516,032✔
153	a17 = -20.0 * pi^2 * a15 * a10 - 10.0 * pi * a15 * a16 * a9 * a7 / (a14)	93,516,032✔
154	a18 = -2 * pi * a9 * a7 - 4.0 * pi * a1	93,516,032✔
155	a19 = a10 + 2.0 * a2 + 5.0	93,516,032✔
156	a20 = a18 * a8 * a3 + 2 * pi * a19 * a8 * a4	93,516,032✔
157	a21 = a15 * (40.0 * pi^3 * a15 * a9 * a7 - 40.0 * pi^2 * a15 * a16 * a10 / (a14) -	93,516,032✔
158	20.0 * pi^2 * a15 * a16 * a9 * a1 * a7 / (a14 * a11) -
159	16.0 * pi^2 * a15 * a16 * a9 * a1 * a7 / (alpha * a12 * a13^3) -
160	10.0 * pi * a15 *
161	(-2.0 * pi^2 * a14 * a2 / a11 - 1.6 * pi^2 * alpha * a12 * a13 * a2 +
162	3.2 * pi^2 * alpha * a12 * a13 * a1^2 / a11 + 0.56 * pi^2 * alpha * a12 * a1^2) *
163	a9 * a7 / (a14) -
164	10.0 * pi * a15 * a16^2 * a9 * a7 / (alpha^2 * g * a13^4 * a11))
165
166	dq1 = -5.0 * a21 + a20 + 4 * pi * a9 * a7 + a10 - 5.0 * a15 * a17 * a16 / (a14)	93,516,032✔
167
168	dq2 = -25.0 * beta * (-2 * pi^2 * a8 * a3 + 4 * pi^3 * a8 * a4) * a13^2 * a11 +	93,516,032✔
169	100.0 * pi * beta * (2 * pi^2 * a8 * a3 + pi * a8 * a4) * a13^2 * a1 +
170	40.0 * pi * beta * (2 * pi^2 * a8 * a3 + pi * a8 * a4) * a13 * a11 * a1 -
171	2 * pi * g * a19 * a9 * a7 + 100.0 * pi^3 * gamma * a12 * a13^2 * a11 * a8 * a4 -
172	300.0 * pi^3 * gamma * a12 * a13^2 * a8 * a1 * a3 -
173	80.0 * pi^3 * gamma * a12 * a13 * a11 * a8 * a1 * a3 -
174	pi^3 * gamma * a12 *
175	(-50.0 * (3.2 * a13 * a2 - 1.28 * a1^2) * a11 * a6 -
176	50.0 * (4.0 * a2 / a11 + 0.16 * a1^2 / a13^2) * a13^2 * a11 * a6 -
177	200.0 * a13^2 * a11 * a6 - 1200.0 * a13^2 * a1 * a5 - 400.0 * a13^2 * a2 * a6 +
178	800.0 * a13^2 * a1^2 * a6 / a11 - 320.0 * a13 * a11 * a1 * a5 +
179	960.0 * a13 * a1^2 * a6) * a8 / 2 - 10.0 * pi * a15 * a17 * a8 * a4 -
180	2.5 * a21 * a8 * a3 + (5.0 * a21 + 5.0 * a15 * a17 * a16 / (a14)) * a8 * a3 / 2 +
181	(a8 * a3 / 2 - pi * a8 * a4) * a19 + a18 * a9 * a3^2 / 2 - (a20) * a8 * a3 / 2 +
182	3 * pi * a19 * a9 * a3 * a4 - 2.5 * a15 * a17 * a16 * a8 * a3 / (a14)
183
184	return SVector(dq1, dq2, zero(dq1))	93,516,032✔
185	end
186
187	function create_cache(mesh, equations::SvaerdKalischEquations1D,	32✔
188	solver, initial_condition,
189	::BoundaryConditionPeriodic,
190	RealT, uEltype)
191	D1 = solver.D1	32✔
192	# Assume D is independent of time and compute D evaluated at mesh points once.
193	D = Array{RealT}(undef, nnodes(mesh))	32✔
194	x = grid(solver)	32✔
195	for i in eachnode(solver)	40✔
196	D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations)	19,168✔
197	end	19,136✔
198	g = gravity_constant(equations)	32✔
199	h = ones(RealT, nnodes(mesh))	32✔
200	hv = zero(h)	32✔
201	b = zero(h)	32✔
202	eta_x = zero(h)	32✔
203	v_x = zero(h)	32✔
204	alpha_eta_x_x = zero(h)	32✔
205	y_x = zero(h)	32✔
206	v_y_x = zero(h)	32✔
207	yv_x = zero(h)	32✔
208	vy = zero(h)	32✔
209	vy_x = zero(h)	32✔
210	y_v_x = zero(h)	32✔
211	h_v_x = zero(h)	32✔
212	hv2_x = zero(h)	32✔
213	v_xx = zero(h)	32✔
214	gamma_v_xx_x = zero(h)	32✔
215	gamma_v_x_xx = zero(h)	32✔
216	alpha_hat = sqrt.(equations.alpha * sqrt.(g * D) .* D .^ 2)	32✔
217	beta_hat = equations.beta * D .^ 3	32✔
218	gamma_hat = equations.gamma * sqrt.(g * D) .* D .^ 3	32✔
219	tmp2 = zero(h)	32✔
220	M = mass_matrix(D1)	32✔
221	M_h = zero(h)	32✔
222	M_beta = copy(beta_hat)	32✔
223	scale_by_mass_matrix!(M_beta, D1)	32✔
224	if D1 isa PeriodicDerivativeOperator \|\|	32✔
225	D1 isa UniformPeriodicCoupledOperator
226	D1_central = D1	28✔
227	D1mat = sparse(D1_central)	28✔
228	minus_MD1betaD1 = D1mat' * Diagonal(M_beta) * D1mat	28✔
229	system_matrix = let cache = (; M_h, minus_MD1betaD1)	28✔
230	assemble_system_matrix!(cache, h,	28✔
231	D1, D1mat, equations)
232	end
233	elseif D1 isa PeriodicUpwindOperators	4✔
234	D1_central = D1.central	4✔
235	D1mat_minus = sparse(D1.minus)	4✔
236	minus_MD1betaD1 = D1mat_minus' * Diagonal(M_beta) * D1mat_minus	4✔
237	system_matrix = let cache = (; M_h, minus_MD1betaD1)	4✔
238	assemble_system_matrix!(cache, h,	4✔
239	D1, D1mat_minus, equations)
240	end
241	else
NEW 242	@error "unknown type of first-derivative operator: $(typeof(D1))"	×
243	end
244	factorization = cholesky(system_matrix)	32✔
245	cache = (; factorization, minus_MD1betaD1, D, h, hv, b, eta_x, v_x,	32✔
246	alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x, y_v_x, h_v_x, hv2_x, v_xx,
247	gamma_v_xx_x, gamma_v_x_xx,
248	alpha_hat, beta_hat, gamma_hat, tmp2, D1_central, M, D1, M_h)
249	if D1 isa PeriodicUpwindOperators	32✔
250	eta_x_upwind = zero(h)	4✔
251	v_x_upwind = zero(h)	4✔
252	cache = (; cache..., eta_x_upwind, v_x_upwind)	4✔
253	end
254	return cache	32✔
255	end
256
257	function assemble_system_matrix!(cache, h, D1, D1mat,	852,166✔
258	::SvaerdKalischEquations1D)
259	(; M_h, minus_MD1betaD1) = cache	852,166✔
260
261	@.. M_h = h	852,166✔
262	scale_by_mass_matrix!(M_h, D1)	852,166✔
263
264	# Floating point errors accumulate a bit and the system matrix
265	# is not necessarily perfectly symmetric but only up to
266	# round-off errors. We wrap it here to avoid issues with the
267	# factorization.
268	return Symmetric(Diagonal(M_h) + minus_MD1betaD1)	852,166✔
269	end
270
271	# Discretization that conserves
272	# - the total water (integral of h) as a linear invariant
273	# - the total momentum (integral of h v) as a nonlinear invariant for flat bathymetry
274	# - the total modified energy
275	# for periodic boundary conditions, see
276	# - Joshua Lampert and Hendrik Ranocha (2024)
277	# Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
278	# [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
279	# TODO: Simplify for the case of flat bathymetry and use higher-order operators
280	function rhs!(dq, q, t, mesh, equations::SvaerdKalischEquations1D,	852,134✔
281	initial_condition, ::BoundaryConditionPeriodic, source_terms,
282	solver, cache)
283	(; D, h, hv, b, eta_x, v_x, alpha_eta_x_x, y_x, v_y_x, yv_x, vy, vy_x,	852,134✔
284	y_v_x, h_v_x, hv2_x, v_xx, gamma_v_xx_x, gamma_v_x_xx, alpha_hat, gamma_hat,
285	tmp1, tmp2, D1_central, M, D1) = cache
286
287	g = gravity_constant(equations)	852,134✔
288	eta, v = q.x	852,134✔
289	deta, dv, dD = dq.x	852,134✔
290	fill!(dD, zero(eltype(dD)))	118,294,528✔
291
292	@trixi_timeit timer() "deta hyperbolic" begin	852,134✔
293	@.. b = equations.eta0 - D	852,134✔
294	@.. h = eta - b	852,134✔
295	@.. hv = h * v	852,134✔
296
297	if D1 isa PeriodicDerivativeOperator \|\|	852,134✔
298	D1 isa UniformPeriodicCoupledOperator
299	mul!(eta_x, D1_central, eta)	851,782✔
300	mul!(v_x, D1_central, v)	851,782✔
301	@.. tmp1 = alpha_hat * eta_x	851,782✔
302	mul!(alpha_eta_x_x, D1_central, tmp1)	851,782✔
303	@.. tmp1 = alpha_hat * alpha_eta_x_x	851,782✔
304	mul!(y_x, D1_central, tmp1)	851,782✔
305	@.. v_y_x = v * y_x	851,782✔
306	@.. vy = v * tmp1	851,782✔
307	mul!(vy_x, D1_central, vy)	851,782✔
308	@.. y_v_x = tmp1 * v_x	851,782✔
309	@.. tmp2 = tmp1 - hv	851,782✔
310	mul!(deta, D1_central, tmp2)	851,782✔
311	elseif D1 isa PeriodicUpwindOperators	352✔
312	@unpack eta_x_upwind, v_x_upwind = cache	352✔
313	mul!(eta_x, D1_central, eta)	352✔
314	mul!(v_x, D1_central, v)	352✔
315	mul!(eta_x_upwind, D1.plus, eta)	352✔
316	@.. tmp1 = alpha_hat * eta_x_upwind	352✔
317	mul!(alpha_eta_x_x, D1_central, tmp1)	352✔
318	@.. tmp1 = alpha_hat * alpha_eta_x_x	352✔
319	mul!(y_x, D1.minus, tmp1)	352✔
320	@.. v_y_x = v * y_x	352✔
321	@.. vy = v * tmp1	352✔
322	mul!(vy_x, D1.minus, vy)	352✔
323	mul!(v_x_upwind, D1.plus, v)	352✔
324	@.. y_v_x = tmp1 * v_x_upwind	352✔
325	# deta[:] = D1.minus * tmp1 - D1_central * hv
326	mul!(deta, D1.minus, tmp1)	352✔
327	mul!(deta, D1_central, hv, -1.0, 1.0)	352✔
328	else
NEW 329	@error "unknown type of first derivative operator: $(typeof(D1))"	×
330	end
331	end
332
333	# split form
334	@trixi_timeit timer() "dv hyperbolic" begin	852,134✔
335	mul!(h_v_x, D1_central, hv)	852,134✔
336	@.. tmp1 = hv * v	852,134✔
337	mul!(hv2_x, D1_central, tmp1)	852,134✔
338	mul!(v_xx, solver.D2, v)	852,134✔
339	@.. tmp1 = gamma_hat * v_xx	852,134✔
340	mul!(gamma_v_xx_x, D1_central, tmp1)	852,134✔
341	@.. tmp1 = gamma_hat * v_x	852,134✔
342	mul!(gamma_v_x_xx, solver.D2, tmp1)	852,134✔
343	@.. dv = -(0.5 * (hv2_x + hv * v_x - v * h_v_x) +	852,134✔
344	g * h * eta_x +
345	0.5 * (v_y_x - vy_x - y_v_x) -
346	0.5 * gamma_v_xx_x -
347	0.5 * gamma_v_x_xx)
348	end
349
350	# no split form
351	# dv[:] = -(D1_central * (hv .* v) - v .* (D1_central * hv)+
352	# g * h .* eta_x +
353	# vy_x - v_y_c -
354	# 0.5 * gamma_v_xx_x -
355	# 0.5 * gamma_v_x_xx)
356
357	@trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,	852,134✔
358	solver)
359	@trixi_timeit timer() "assemble system matrix" begin	852,134✔
360	system_matrix = assemble_system_matrix!(cache, h, D1, D1_central, equations)	852,134✔
361	end
362	@trixi_timeit timer() "solving elliptic system" begin	852,134✔
363	tmp1 .= dv	852,134✔
364	solve_system_matrix!(dv, system_matrix, tmp1,	852,134✔
365	equations, D1, cache)
366	end
367
368	return nothing	852,134✔
369	end
370
371	@inline function prim2cons(q, equations::SvaerdKalischEquations1D)	4✔
372	eta, v = q	4✔
373
374	b = bathymetry(q, equations)	4✔
375	h = eta - b	4✔
376	hv = h * v	4✔
377	return SVector(h, hv, b)	4✔
378	end
379
380	@inline function cons2prim(u, equations::SvaerdKalischEquations1D)	4✔
381	h, hv, b = u	4✔
382
383	eta = h + b	4✔
384	v = hv / h	4✔
385	D = equations.eta0 - b	4✔
386	return SVector(eta, v, D)	4✔
387	end
388
389	@inline function waterheight_total(q, ::SvaerdKalischEquations1D)	330,984✔
390	return q[1]	330,984✔
391	end
392
393	@inline function velocity(q, ::SvaerdKalischEquations1D)	140,600✔
394	return q[2]	140,600✔
395	end
396
397	@inline function bathymetry(q, equations::SvaerdKalischEquations1D)	133,804✔
398	D = q[3]	133,804✔
399	return equations.eta0 - D	133,804✔
400	end
401
402	# The modified entropy/energy takes the whole `q` for every point in space
403	"""
404	energy_total_modified(q_global, equations::SvaerdKalischEquations1D, cache)
405
406	Return the modified total energy of the primitive variables `q_global` for the
407	[`SvaerdKalischEquations1D`](@ref). It contains an additional term containing a
408	derivative compared to the usual [`energy_total`](@ref) modeling
409	non-hydrostatic contributions. The `energy_total_modified`
410	is a conserved quantity (for periodic boundary conditions).
411
412	It is given by
413	```math
414	\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{2} \\hat\\beta v_x^2.
415	```
416
417	`q_global` is a vector of the primitive variables at ALL nodes.
418	`cache` needs to hold the SBP operators used by the `solver`.
419	"""
420	@inline function energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D,	692✔
421	cache)
422	# unpack physical parameters and SBP operator `D1`
423	g = gravity_constant(equations)	692✔
424	(; D1, h, b, v_x, beta_hat) = cache	692✔
425
426	# `q_global` is an `ArrayPartition`. It collects the individual arrays for
427	# the total water height `eta = h + b` and the velocity `v`.
428	eta, v, D = q_global.x	692✔
429	@.. b = equations.eta0 - D	692✔
430	@.. h = eta - b	692✔
431
432	if D1 isa PeriodicUpwindOperators	692✔
433	mul!(v_x, D1.minus, v)	12✔
434	else
435	mul!(v_x, D1, v)	680✔
436	end
437
438	@.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 2 * beta_hat * v_x^2	692✔
439	return e	692✔
440	end

JoshuaLampert / DispersiveShallowWater.jl / 10507790104

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous