10813110923

Committed 11 Sep 2024 01:53PM UTC coverage: 97.939% (+0.2%) from 97.78%

Build # 10813110923

Build Type

Pull #150

github

Committed by

web-flow

Commit Message

Merge 1e0179980 into 11355c045

Pull Request Pull Request #150: Add `BBMEquation1D`

Run Details

91 of 92 new or added lines in 4 files covered. (98.91%)

2 existing lines in 1 file now uncovered.

1711 of 1747 relevant lines covered (97.94%)

3372397.2 hits per line

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

98.51

/src/equations/bbm_1d.jl

# TODO: Use this nondimensionalized version or use dimensional version including gravity and depth?
@doc raw"""
    BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0)

BBM (Benjamin–Bona–Mahony) equation in one spatial dimension. The equations are given by
```math
\begin{aligned}
  \eta_t + \eta\eta_x + \eta_x - \eta_{xxt} &= 0.
\end{aligned}
```

The unknown quantity of the BBM equation is the total water height ``\eta``. The water height
is given by ``h = \eta - \eta_0``, where ``\eta_0`` is the constant still-water surface.

Currently, the equations only support a flat bathymetry ``b = 0``, see [`bathymetry_flat`](@ref).

The BBM equation is first described in Benjamin, Bona, and Mahony (1972).
The semidiscretization implemented here is developed in Ranocha, Mitsotakis, and Ketcheson (2020).
It conserves
- the total water mass (integral of ``h``) as a linear invariant
- a quadratic invariant (integral of ``\eta^2 + \eta_x^2``), which is called here [`energy_total_modified`](@ref)
  (and [`entropy_modified`](@ref)) because it contains derivatives of the solution

for periodic boundary conditions.

- Thomas B. Benjamin, Jerry L. Bona and John J. Mahony (1972)
  Model equations for long waves in nonlinear dispersive systems
  [DOI: 10.1098/rsta.1972.0032](https://doi.org/10.1098/rsta.1972.0032)
- Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
  A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
  [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
"""
struct BBMEquation1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <:
       AbstractBBMEquation{1, 1}
    bathymetry_type::Bathymetry # type of bathymetry
    eta0::RealT # constant still-water surface
end

function BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0)
    BBMEquation1D(bathymetry_type, eta0)
end

"""
    initial_condition_convergence_test(x, t, equations::BBMEquation1D, mesh)

A travelling-wave solution used for convergence tests in a periodic domain.

See section 4.1.3 in (there is an error in paper: it should be sech^2 instead of cosh):
- Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
  A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
  [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
"""
function initial_condition_convergence_test(x, t, equations::BBMEquation1D, mesh)
    c = 1.2
    A = 3 * (c - 1)
    K = 0.5 * sqrt(1 - 1 / c)
    x_t = mod(x - c * t - xmin(mesh), xmax(mesh) - xmin(mesh)) + xmin(mesh)
    eta = A * sech(K * x_t)^2
    return SVector(eta)
end

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

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

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

A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref).
"""
function source_terms_manufactured(q, x, t, equations::BBMEquation1D)
    a3 = cospi(t - 2 * x)
    a4 = sinpi(t - 2 * x)
    a5 = sinpi(2 * t - 4 * x)
    a14 = exp(t / 2)
    dq1 = -2 * pi^2 * (a4 + 2 * pi * a3) * a14 - a14 * a4 / 2 + pi * a14 * a3 -
          pi * exp(t) * a5

    return SVector(dq1)
end

function create_cache(mesh, equations::BBMEquation1D,
                      solver, initial_condition,
                      ::BoundaryConditionPeriodic,
                      RealT, uEltype)
    invImD2 = lu(I - sparse(solver.D2))
    eta2 = zeros(RealT, nnodes(mesh))
    eta2_x = zero(eta2)
    eta_x = zero(eta2)
    etaeta_x = zero(eta2)
    cache = (; invImD2, eta2, eta2_x, eta_x, etaeta_x, solver.D1)
    if solver.D1 isa PeriodicUpwindOperators
        eta_x_upwind = zero(eta2)
        cache = (; cache..., eta_x_upwind)
    end
    return cache
end

# Discretization that conserves the mass for eta and the modified energy for periodic boundary conditions, see
# - Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
#   A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
#   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
function rhs!(dq, q, t, mesh, equations::BBMEquation1D, initial_condition,
              ::BoundaryConditionPeriodic, source_terms, solver, cache)
    (; invImD2, eta2, eta2_x, eta_x, etaeta_x) = cache
    if solver.D1 isa PeriodicUpwindOperators
        (; eta_x_upwind) = cache
    end

    eta, = q.x
    deta, = dq.x
    # energy and mass conservative semidiscretization
    @trixi_timeit timer() "hyperbolic" begin
        @.. eta2 = eta^2
        if solver.D1 isa PeriodicDerivativeOperator ||
           solver.D1 isa UniformPeriodicCoupledOperator ||
           solver.D1 isa FourierDerivativeOperator
            mul!(eta2_x, solver.D1, eta2)
            mul!(eta_x, solver.D1, eta)
            @.. etaeta_x = eta * eta_x
            # split form
            @.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x)
        elseif solver.D1 isa PeriodicUpwindOperators
            mul!(eta2_x, solver.D1.central, eta2)
            mul!(eta_x_upwind, solver.D1.minus, eta)
            mul!(eta_x, solver.D1.central, eta)
            @.. etaeta_x = eta * eta_x
            # split form
            @.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x_upwind)
        else
            @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
        end
    end

    @trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,
                                                       solver)

    @trixi_timeit timer() "elliptic" begin
        solve_system_matrix!(deta, invImD2, equations)
    end
    return nothing
end

"""
    energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)

Return the modified total energy `e` of the primitive variables `q_global` for the
[`BBMEquation1D`](@ref). The `energy_total_modified`
is a conserved quantity (for periodic boundary conditions).

It is given by
```math
\\frac{1}{2} \\eta^2 + \\frac{1}{2} \\eta_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`.

See also [`energy_total_modified`](@ref).
"""
function energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)
    eta, = q_global.x

    (; D1, eta_x) = cache
    if D1 isa PeriodicUpwindOperators
        mul!(eta_x, D1.central, eta)
    else
        mul!(eta_x, D1, eta)
    end

    @.. e = 0.5 * (eta^2 + eta_x^2)
    return nothing
end

"""
    invariant_cubic(q, equations::BBMEquation1D)

Return the cubic invariant of the primitive variables `q` for the
[`BBMEquation1D`](@ref). The cubic invariant is given by
```math
(\\eta + 1)^3.
```
"""
function invariant_cubic(q, equations::BBMEquation1D)
    eta = waterheight_total(q, equations)
    return (eta + 1)^3
end

varnames(::typeof(invariant_cubic), equations) = ("(η + 1)³",)

1	# TODO: Use this nondimensionalized version or use dimensional version including gravity and depth?
2	@doc raw"""
3	BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0)
4
5	BBM (Benjamin–Bona–Mahony) equation in one spatial dimension. The equations are given by
6	```math
7	\begin{aligned}
8	\eta_t + \eta\eta_x + \eta_x - \eta_{xxt} &= 0.
9	\end{aligned}
10	```
11
12	The unknown quantity of the BBM equation is the total water height ``\eta``. The water height
13	is given by ``h = \eta - \eta_0``, where ``\eta_0`` is the constant still-water surface.
14
15	Currently, the equations only support a flat bathymetry ``b = 0``, see [`bathymetry_flat`](@ref).
16
17	The BBM equation is first described in Benjamin, Bona, and Mahony (1972).
18	The semidiscretization implemented here is developed in Ranocha, Mitsotakis, and Ketcheson (2020).
19	It conserves
20	- the total water mass (integral of ``h``) as a linear invariant
21	- a quadratic invariant (integral of ``\eta^2 + \eta_x^2``), which is called here [`energy_total_modified`](@ref)
22	(and [`entropy_modified`](@ref)) because it contains derivatives of the solution
23
24	for periodic boundary conditions.
25
26	- Thomas B. Benjamin, Jerry L. Bona and John J. Mahony (1972)
27	Model equations for long waves in nonlinear dispersive systems
28	[DOI: 10.1098/rsta.1972.0032](https://doi.org/10.1098/rsta.1972.0032)
29	- Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
30	A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
31	[DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
32	"""
33	struct BBMEquation1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <:
34	AbstractBBMEquation{1, 1}
35	bathymetry_type::Bathymetry # type of bathymetry	25✔
36	eta0::RealT # constant still-water surface
37	end
38
39	function BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0)	40✔
40	BBMEquation1D(bathymetry_type, eta0)	20✔
41	end
42
43	"""
44	initial_condition_convergence_test(x, t, equations::BBMEquation1D, mesh)
45
46	A travelling-wave solution used for convergence tests in a periodic domain.
47
48	See section 4.1.3 in (there is an error in paper: it should be sech^2 instead of cosh):
49	- Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
50	A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
51	[DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
52	"""
53	function initial_condition_convergence_test(x, t, equations::BBMEquation1D, mesh)	196,608✔
54	c = 1.2	196,608✔
55	A = 3 * (c - 1)	196,608✔
56	K = 0.5 * sqrt(1 - 1 / c)	196,608✔
57	x_t = mod(x - c * t - xmin(mesh), xmax(mesh) - xmin(mesh)) + xmin(mesh)	393,184✔
58	eta = A * sech(K * x_t)^2	196,608✔
59	return SVector(eta)	196,608✔
60	end
61
62	"""
63	initial_condition_manufactured(x, t, equations::BBMEquation1D, mesh)
64
65	A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref).
66	"""
67	function initial_condition_manufactured(x, t,	16,384✔
68	equations::BBMEquation1D,
69	mesh)
70	eta = exp(t / 2) * sinpi(2 * (x - t / 2))	16,384✔
71	return SVector(eta)	16,384✔
72	end
73
74	"""
75	source_terms_manufactured(q, x, t, equations::BBMEquation1D, mesh)
76
77	A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref).
78	"""
79	function source_terms_manufactured(q, x, t, equations::BBMEquation1D)	434,176✔
80	a3 = cospi(t - 2 * x)	434,176✔
81	a4 = sinpi(t - 2 * x)	434,176✔
82	a5 = sinpi(2 * t - 4 * x)	434,176✔
83	a14 = exp(t / 2)	434,176✔
84	dq1 = -2 * pi^2 * (a4 + 2 * pi * a3) * a14 - a14 * a4 / 2 + pi * a14 * a3 -	434,176✔
85	pi * exp(t) * a5
86
87	return SVector(dq1)	434,176✔
88	end
89
90	function create_cache(mesh, equations::BBMEquation1D,	24✔
91	solver, initial_condition,
92	::BoundaryConditionPeriodic,
93	RealT, uEltype)
94	invImD2 = lu(I - sparse(solver.D2))	24✔
95	eta2 = zeros(RealT, nnodes(mesh))	24✔
96	eta2_x = zero(eta2)	24✔
97	eta_x = zero(eta2)	24✔
98	etaeta_x = zero(eta2)	24✔
99	cache = (; invImD2, eta2, eta2_x, eta_x, etaeta_x, solver.D1)	24✔
100	if solver.D1 isa PeriodicUpwindOperators	24✔
101	eta_x_upwind = zero(eta2)	8✔
102	cache = (; cache..., eta_x_upwind)	8✔
103	end
104	return cache	24✔
105	end
106
107	# Discretization that conserves the mass for eta and the modified energy for periodic boundary conditions, see
108	# - Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
109	# A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
110	# [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
111	function rhs!(dq, q, t, mesh, equations::BBMEquation1D, initial_condition,	22,108✔
112	::BoundaryConditionPeriodic, source_terms, solver, cache)
113	(; invImD2, eta2, eta2_x, eta_x, etaeta_x) = cache	22,108✔
114	if solver.D1 isa PeriodicUpwindOperators	22,108✔
115	(; eta_x_upwind) = cache	5,520✔
116	end
117
118	eta, = q.x	22,108✔
119	deta, = dq.x	22,108✔
120	# energy and mass conservative semidiscretization
121	@trixi_timeit timer() "hyperbolic" begin	22,108✔
122	@.. eta2 = eta^2	22,108✔
123	if solver.D1 isa PeriodicDerivativeOperator \|\|	22,108✔
124	solver.D1 isa UniformPeriodicCoupledOperator \|\|
125	solver.D1 isa FourierDerivativeOperator
126	mul!(eta2_x, solver.D1, eta2)	16,588✔
127	mul!(eta_x, solver.D1, eta)	16,588✔
128	@.. etaeta_x = eta * eta_x	16,588✔
129	# split form
130	@.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x)	16,588✔
131	elseif solver.D1 isa PeriodicUpwindOperators	5,520✔
132	mul!(eta2_x, solver.D1.central, eta2)	5,520✔
133	mul!(eta_x_upwind, solver.D1.minus, eta)	5,520✔
134	mul!(eta_x, solver.D1.central, eta)	5,520✔
135	@.. etaeta_x = eta * eta_x	5,520✔
136	# split form
137	@.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x_upwind)	5,520✔
138	else
NEW 139	@error "unknown type of first-derivative operator: $(typeof(solver.D1))"	×
140	end
141	end
142
143	@trixi_timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations,	22,108✔
144	solver)
145
146	@trixi_timeit timer() "elliptic" begin	22,108✔
147	solve_system_matrix!(deta, invImD2, equations)	22,108✔
148	end
149	return nothing	22,108✔
150	end
151
152	"""
153	energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)
154
155	Return the modified total energy `e` of the primitive variables `q_global` for the
156	[`BBMEquation1D`](@ref). The `energy_total_modified`
157	is a conserved quantity (for periodic boundary conditions).
158
159	It is given by
160	```math
161	\\frac{1}{2} \\eta^2 + \\frac{1}{2} \\eta_x^2.
162	```
163
164	`q_global` is a vector of the primitive variables at ALL nodes.
165	`cache` needs to hold the SBP operators used by the `solver`.
166
167	See also [`energy_total_modified`](@ref).
168	"""
169	function energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)	9,522✔
170	eta, = q_global.x	9,522✔
171
172	(; D1, eta_x) = cache	9,522✔
173	if D1 isa PeriodicUpwindOperators	9,522✔
174	mul!(eta_x, D1.central, eta)	92✔
175	else
176	mul!(eta_x, D1, eta)	9,430✔
177	end
178
179	@.. e = 0.5 * (eta^2 + eta_x^2)	9,522✔
180	return nothing	9,522✔
181	end
182
183	"""
184	invariant_cubic(q, equations::BBMEquation1D)
185
186	Return the cubic invariant of the primitive variables `q` for the
187	[`BBMEquation1D`](@ref). The cubic invariant is given by
188	```math
189	(\\eta + 1)^3.
190	```
191	"""
192	function invariant_cubic(q, equations::BBMEquation1D)	188,420✔
193	eta = waterheight_total(q, equations)	188,420✔
194	return (eta + 1)^3	188,420✔
195	end
196
197	varnames(::typeof(invariant_cubic), equations) = ("(η + 1)³",)	4✔

JoshuaLampert / DispersiveShallowWater.jl / 10813110923

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