diff --git a/README.md b/README.md index a9246ff8..d48980fb 100644 --- a/README.md +++ b/README.md @@ -39,8 +39,12 @@ The result can be visualized by using the package Plots.jl julia> using Plots julia> plot(semi => sol) ``` -The command `plot` expects a `Pair` consisting of a `Semidiscretization` and an `ODESolution`. The visualization can also be customized, see the [documentation](https://JoshuaLampert.github.io./DispersiveShallowWater.jl/stable/) for more details. Other examples can be found in the subdirectory [examples/](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/tree/main/examples). A list of all examples is returned by running `get_examples()`. You can pass the filename of one of the examples or your own simulation file to `include` in order to run it, e.g., `include(joinpath(examples_dir(), "svaerd_kalisch_1d", "svaerd_kalisch_1d_dingemans_relaxation.jl"))`. +The command `plot` expects a `Pair` consisting of a `Semidiscretization` and an `ODESolution`. The visualization can also be customized, see the [documentation](https://JoshuaLampert.github.io./DispersiveShallowWater.jl/stable/overview.html#Visualize-results) for more details. Other examples can be found in the subdirectory [examples/](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/tree/main/examples). A list of all examples is returned by running `get_examples()`. You can pass the filename of one of the examples or your own simulation file to `include` in order to run it, e.g., `include(joinpath(examples_dir(), "svaerd_kalisch_1d", "svaerd_kalisch_1d_dingemans_relaxation.jl"))`. # Authors The package is developed and maintained by Joshua Lampert and was initiated as part of the master thesis "Structure-preserving Numerical Methods for Dispersive Shallow Water Models" (2023). Some parts of this repository are based on parts of [Dispersive-wave-schemes-notebooks. A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations](https://github.com/ranocha/Dispersive-wave-schemes-notebooks) by Hendrik Ranocha, Dimitrios Mitsotakis and David Ketcheson. The code structure is inspired by [Trixi.jl](https://github.com/trixi-framework/Trixi.jl/). + +# License and contributing + +DispersiveShallowWater.jl is published under the MIT license (see [License](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/LICENSE)). We are pleased to accept contributions from everyone, preferably in the form of a PR. diff --git a/docs/make.jl b/docs/make.jl index 98bcf4a4..899655f4 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,26 +1,28 @@ using DispersiveShallowWater using Documenter -DocMeta.setdocmeta!(DispersiveShallowWater, :DocTestSetup, :(using DispersiveShallowWater); recursive=true) +# Define module-wide setups such that the respective modules are available in doctests +DocMeta.setdocmeta!(DispersiveShallowWater, :DocTestSetup, :(using DispersiveShallowWater); + recursive = true) makedocs(; - modules=[DispersiveShallowWater], - authors="Joshua Lampert ", - repo="https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/{commit}{path}#{line}", - sitename="DispersiveShallowWater.jl", - format=Documenter.HTML(; - prettyurls=get(ENV, "CI", "false") == "true", - canonical="https://JoshuaLampert/DispersiveShallowWater.jl", - edit_link="main", - assets=String[], - ), - pages=[ - "Home" => "index.md", - "Reference" => "ref.md", - ], -) + modules = [DispersiveShallowWater], + authors = "Joshua Lampert ", + repo = "https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/{commit}{path}#{line}", + sitename = "DispersiveShallowWater.jl", + format = Documenter.HTML(; + prettyurls = get(ENV, "CI", "false") == "true", + canonical = "https://JoshuaLampert.github.io/DispersiveShallowWater.jl/stable", + edit_link = "main", + assets = String[]), + pages = [ + "Home" => "index.md", + "Overview" => "overview.md", + "Reproduce figures" => "reproduce.md", + "Reference" => "ref.md", + "License" => "license.md", + ]) -#deploydocs(; -# repo="github.com/JoshuaLampert/DispersiveShallowWater.jl", -# devbranch="main", -#) +deploydocs(; + repo = "github.com/JoshuaLampert/DispersiveShallowWater.jl", + devbranch = "main") diff --git a/docs/src/index.md b/docs/src/index.md index bca66269..bfdd601e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,9 +1,50 @@ # DispersiveShallowWater.jl -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://trixi-framework.github.io/DispersiveShallowWater.jl/stable/) -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://trixi-framework.github.io/DispersiveShallowWater.jl/dev/) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JoshuaLampert.github.io/DispersiveShallowWater.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JoshuaLampert.github.io/DispersiveShallowWater.jl/dev/) [![Build Status](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) +[![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT) -[**DispersiveShallowWater.jl**](https://github.com/JoshuaLampert/DispersiveShallowWater.jl) is a numerical simulation framework for dispersive shallow water models. The main objective lies in structure-preserving numerical methods. -This is the reproducibality repository for the master thesis "Structure-preserving Numerical Methods for Dispersive Shallow Water Models" by Joshua Lampert (2023). +[**DispersiveShallowWater.jl**](https://github.com/JoshuaLampert/DispersiveShallowWater.jl) is a [Julia](https://julialang.org/) package that implements structure-preserving numerical methods for dispersive shallow water models. To date, it provides provably conservative, entropy-conserving and well-balanced numerical schemes of the [BBM-BBM equations with varying bottom topography](https://iopscience.iop.org/article/10.1088/1361-6544/ac3c29), and the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924). The semidiscretizations are based on summation by parts (SBP) operators, which are implemented in [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/). In order to obtain fully discrete schemes, the time integration methods from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) are used to solve the resulting ordinary differential equations. Fully discrete entropy-conservative methods can be obtained by using the [relaxation method](https://epubs.siam.org/doi/10.1137/19M1263662) provided by DispersiveShallowWater.jl. +# Installation + +If you have not yet installed Julia, you first need to [download Julia](https://julialang.org/downloads/). Please [follow the instructions for your operating system](https://julialang.org/downloads/platform/). DispersiveShallowWater.jl works with Julia v1.8 and newer. You can install DispersiveShallowWater.jl by executing the following commands from the Julia REPL +```julia +julia> using Pkg + +julia> Pkg.add(url="https://github.com/JoshuaLampert/DispersiveShallowWater.jl") + +julia> Pkg.add(["OrdinaryDiffEq", "Plots"]) +``` +The last command installs also the package "OrdinaryDiffEq.jl" used for time-integration and "Plots.jl" to visualize the results. If you want to use other SBP operators than the default operators that DispersiveShallowWater.jl uses, you also need SummationByPartsOperators.jl, which can be installed running +```julia +julia> Pkg.add("SummationByPartsOperators") +``` + +# Usage + +In the Julia REPL, first load the package DispersiveShallowWater.jl +```julia +julia> using DispersiveShallowWater +``` + +You can run a basic simulation that solves the BBM-BBM equations by executing +```julia +julia> include(default_example()); +``` +The result can be visualized by using the package Plots.jl +```julia +julia> using Plots +julia> plot(semi => sol) +``` +The command `plot` expects a `Pair` consisting of a [`Semidiscretization`](@ref) and an `ODESolution`. The visualization can also be customized, see the [documentation](@ref visualize_results) for more details. Other examples can be found in the subdirectory [examples/](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/tree/main/examples). A list of all examples is returned by running [`get_examples()`](@ref). You can pass the filename of one of the examples or your own simulation file to `include` in order to run it, e.g., `include(joinpath(examples_dir(), "svaerd_kalisch_1d", "svaerd_kalisch_1d_dingemans_relaxation.jl"))`. + +# Authors + +The package is developed and maintained by Joshua Lampert and was initiated as part of the master thesis "Structure-preserving Numerical Methods for Dispersive Shallow Water Models" (2023). Some parts of this repository are based on parts of [Dispersive-wave-schemes-notebooks. A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations](https://github.com/ranocha/Dispersive-wave-schemes-notebooks) by Hendrik Ranocha, Dimitrios Mitsotakis and David Ketcheson. The code structure is inspired by [Trixi.jl](https://github.com/trixi-framework/Trixi.jl/). + +# License and contributing + +DispersiveShallowWater.jl is published under the MIT license (see [License](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/LICENSE)). We are pleased to accept contributions from everyone, preferably in the form of a PR. diff --git a/docs/src/license.md b/docs/src/license.md new file mode 100644 index 00000000..922bfc18 --- /dev/null +++ b/docs/src/license.md @@ -0,0 +1,23 @@ +# License + +MIT License + +Copyright (c) 2023-present Joshua Lampert + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/docs/src/overview.md b/docs/src/overview.md new file mode 100644 index 00000000..7bb0fb77 --- /dev/null +++ b/docs/src/overview.md @@ -0,0 +1,240 @@ +# Running a simulation + +## Introduction + +In this tutorial we describe how to numerically solve the BBM-BBM (Benjamin-Bona-Mahony) equations with variable bottom topography in one dimension, which has been proposed in [^IsrawiKalischKatsaounisMitsotakis2021] for two spatial dimensions. The equations describe a dispersive shallow water model, i.e. they extend the well-known shallow water equations in the sense that dispersion is modeled. The shallow water equations are a system of first order hyperbolic partial differential equations that can be written in the form of a balance law. In contrast, the BBM-BBM equations additionally include third-order mixed derivatives. In primitive variables ``(\eta, v)`` they can be written as: + +```math +\begin{aligned} + \eta_t + ((\eta + D)v)_x - \frac{1}{6}(D^2\eta_{xt})_x &= 0,\\ + v_t + g\eta_x + \left(\frac{1}{2}v^2\right)_x - \frac{1}{6}(D^2v_t)_{xx} &= 0. +\end{aligned} +``` + +Here, ``\eta = h + b`` describes the total water height, ``h`` the water height above the bottom topography (bathymetry), ``b = -D`` the bathymetry and ``v`` the velocity in horizontal direction. The gravitational acceleration is denoted as ``g``. A sketch of the water height and the bathymetry can be found below. + +![](bathymetry.png) + +In order to conduct a numerical simulation with DispersiveShallowWater.jl, we perform the following steps. + +First, we load the necessary libraries: + +```@example overview +using DispersiveShallowWater, OrdinaryDiffEq +``` + +## Define physical setup + +As a first step of a numerical simulation, we define the physical setup we want to solve. This includes the set of equations, potentially including physical parameters, initial and boundary conditions as well as the domain. In DispersiveShallowWater.jl this can be done as follows: + +```@example overview +equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) + +function initial_condition_shoaling(x, t, equations::BBMBBMVariableEquations1D, mesh) + A = 0.07 # amplitude of wave + x0 = -30 # initial center + eta = A * exp(-0.1*(x - x0)^2) + v = 0 + D = x <= 0.0 ? 0.7 : 0.7 - 1/50 * x + return SVector(eta, v, D) +end + +initial_condition = initial_condition_shoaling +boundary_conditions = boundary_condition_periodic + +coordinates_min = -130.0 +coordinates_max = 20.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N + 1) +``` + +The first line specifies that we want to solve the BBM-BBM equations with variable bathymetry using a gravitaional acceleration ``g = 9.81``. Afterwards, we define the initial condition, which is described as a function with the spatial variable `x`, the time `t`, the `equations` and a `mesh` as parameters. If an analytical solution is available, the time variable `t` can be used, and the initial condition can serve as an analytical to compare the numerical solution with. An initial condition in DispersiveShallowWater.jl is supposed to return an `SVector` holding the values for each of the unknown variables. Since the bathymetry is treated as a variable (with time derivative 0) for convenience, we need to provide the value for the primitive variables `eta` and `v` as well as for `D`. Otherwise, you can just keep the time variable unused. In this case, the initial condition describes a travelling wave that moves towards a beach, which is modeled by a linearly increasing bathymetry. + +Next, we choose periodic boundary conditions since, up to now, DispersiveShallowWater.jl only provides support for this type of boundary conditions. Lastly, we define the physical domain as the interval from -130 to 20 and we choose 512 intermediate nodes. The mesh is homogeneous, i.e. the distance between each two nodes is constant. We choose the left boundary very far to the left in order to avoid an interaction of the left- and right-travelling waves. + +## Define numerical solver + +In the next step, we build a [`Semidiscretization`](@ref) that bundles all ingredients for the spatial discretization of the model. Especially, we need to define a [`Solver`](@ref). The simplest way to define a solver is to call the constructor by providing the mesh and a desired order of accuracy. In the following example, we use an accuracy order of 4. The default constructor simply creates periodic first- and second-derivative central finite difference summation by parts operators of the provided order of accuracy. How to use other summation by parts operators, is described in the section on [how to customize the solver](@ref customize_solver). + +```@example overview +solver = Solver(mesh, 4) + +semi = Semidiscretization(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) +``` + +Finally, we put the `mesh`, the `equations`, the `initial_condition`, the `solver` and the `boundary_conditions` together in a semidiscretization `semi`. + +## Solve system of ordinary differential equations + +Once we have obtained a semidiscretization, we can solve the resulting system of ordinary differential equations. To do so, we specify the time interval that we want to simulate and obtain an `ODEProblem` from the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/) by calling [`semidiscretize`](@ref) on the semidiscretization and the time span. Additionally, we can analyze the numerical solution using an [`AnalysisCallback`](@ref). The analysis includes computing the ``L_2`` and ``L_\infty`` errors of the different variables of the solution compared to the initial condition (or, if available, at the same time analytical solution) and potentially additional errors and quantities. Additional errors can be passed by the keyword argument `extra_analysis_errors` and additional integral quantities that should be analyzed can be passed by keyword argument `extra_analysis_integrals`. In this example we pass the `conservation_error`, which computes the temporal change of the total amount (i.e. integral) of the different variables over time. In addition, the integrals of the total waterheight ``\eta`` [`waterheight_total`](@ref), the [`velocity`](@ref) and the [`entropy`](@ref) are computed and saved for each time step. The total waterheight and the total velocity are linear invariants of the BBM-BBM equations, i.e. they do not change over time. The total entropy + +```math +\mathcal E(t; \eta, v) = \frac{1}{2}\int_\Omega g\eta^2 + (\eta + D)v^2\textrm{d}x +``` + +is a nonlinear invariant and should be constant over time as well. During the simulation, the `AnalysisCallback` will print the results to the terminal. + +Finally, the `ode` can be `solve`d using the interface from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), i.e. we can specify a time-stepping scheme we want to use (in this case we use the adaptive explicit Runge-Kutta method `Tsit5` by Tsitouras of order 5(4)), the tolerances for the adaptive time-stepping and the time values, where the solution values should be saved. Here, we choose 100 equidistant points. + +```@example overview +tspan = (0.0, 25.0) +ode = semidiscretize(semi, tspan) +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +callbacks = CallbackSet(analysis_callback) + +saveat = range(tspan..., length = 100) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) +``` + +After solving the equations, `sol` contains the solution for each of the three variables at every spatial point for each of the 100 points in time. The errors and integrals recorded by the `AnalysisCallback` can be obtained as `NamedTuple`s by [`errors(analysis_callback)`](@ref) and [`integrals(analysis_callback)`](@ref). + +## [Visualize results](@id visualize_results) + +After running the simulation, the results can be visualized using [Plots.jl](https://github.com/JuliaPlots/Plots.jl), which needs to be imported first. Then, we can plot the solution at the final time by calling `plot` on a `Pair` of the `Semidiscretization` and the corresponding `ODESolution` `sol`. The result is depicted in the following picture. + +```@example overview +using Plots + +plot(semi => sol) +savefig("shoaling_solution.png") +``` + +![](shoaling_solution.png) + +By default, this will plot the bathymetry, but not the initial (analytical) solution. You can adjust this by passing the boolean values `plot_bathymetry` (if `true` always plot to first subplot) and `plot_initial`. You can also provide a `conversion` function that converts the solution. A conversion function should take the primitive variables `q` at one node and the `equations` as input and should return an `SVector` of any length. There should also exist a function `varnames(conversion, equations)` that returns a `Tuple` of the variable names used for labelling. The conversion function can, e.g., be [`prim2cons`](@ref) or [`waterheight_total`](@ref) if one only wants to plot the total waterheight. The resulting plot will have one subplot for each of the returned variables of the conversion variable. By default, this is just [`prim2prim`](@ref), i.e. the identity. The limits of the y-axis can be set by the keyword argument `yli` (**note**: `ylim` and similar are already occupied by Plots.jl and do not work with different y limits for multiple subplots), which should be given as a vector of `Tuple`s with the same length as there are subplots. Plotting an animation over time can, e.g., be done by the following command, which uses `step` to plot the solution at a specific time step. + +```@example overview +anim = @animate for step in 1:length(sol.u) + plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), yli = [(-0.8, 0.1)]) +end +gif(anim, "shoaling_solution.gif", fps = 25) +``` + +It is also possible to plot the solution variables at a fixed spatial point over time by calling `plot(semi => sol, x)` for some `x`-value, see [create_figures.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/create_figures.jl) for some examples. + +Often, it is interesting to have a look at how the quantities that are recorded by the `AnalysisCallback` evolve in time. To do this end, you can `plot` the `AnalysisCallback` by + +```@example overview +plot(analysis_callback) +savefig("analysis_callback.png") +``` + +This creates the following figure. + +![](analysis_callback.png) + +You can see that the linear invariants ``\int_\Omega\eta\textrm{d}x`` and ``\int_\Omega v\textrm{d}x`` are indeed conserved exactly. The entropy, however, starts growing at around ``t = 17`` up to approximately 5e-5. This is because of the fact that, during the time integration, a nonlinear invariant is not necessarily conserved, even if the semidiscretization conserves the quantity exactly. How to obtain a fully-discrete structure-preserving numerical scheme is explained in the following section. + +## Use entropy-conserving time integration + +Two obtain entropy-conserving time-stepping schemes DispersiveShallowWater.jl uses the relaxation method introduced in [^Ketcheson2019] and further developed in [^RanochaSayyariDalcinParsaniKetcheson2020]. The relaxation method is implemented as a [`RelaxationCallback`](@ref), which takes a function representing the conserved quantity as the keyword argument `invariant`. Therefore, we can run the same example as above, but using relaxation on the entropy by simply adding another callback to the `CallbackSet`: + +```@example overview +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +relaxation_callback = RelaxationCallback(invariant = entropy) +callbacks = CallbackSet(relaxation_callback, analysis_callback) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) +``` + +When you use both, an `AnalysisCallback` and a `RelaxationCallback`, note that the `relaxation_callback` needs to come first inside the `CallbackSet` as it needs to be invoked prior to the `analysis_callback`, such that the `analysis_callback` analyzes the solution with the already updated values. + +Plotting the `analysis_callback` again, we can see that now also the `entropy` is conserved up to machine precision. + +```@example overview +plot(analysis_callback) +savefig("analysis_callback_relaxation.png") +``` + +![](analysis_callback_relaxation.png) + + +## [Customize solver](@id customize_solver) + +In the semidiscretization created above, we used the default SBP operators, which are periodic finite difference operators. Using different SBP operators for the semidiscretization can be done leveraging [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/), which needs to be imported first: + +```@example overview +using SummationByPartsOperators: legendre_derivative_operator, UniformPeriodicMesh1D, couple_discontinuously, PeriodicUpwindOperators +``` + +As an example, let us create a semidiscretization based on discontinuous Galerkin (DG) upwind operators. A semidiscretization implemented in DispersiveShallowWater.jl needs one first-derivative and one second-derivative SBP operator. To build the first-derivative operator, we first create a `LegendreDerivativeOperator` with polynomial degree 3 on a reference element `[-1.0, 1.0]` and a `UniformPeriodicMesh1D` for the coupling. + +```@example overview +mesh = Mesh1D(coordinates_min, coordinates_max, N) +accuracy_order = 4 +Dop = legendre_derivative_operator(-1.0, 1.0, accuracy_order) +sbp_mesh = UniformPeriodicMesh1D(mesh.xmin, mesh.xmax, div(mesh.N, accuracy_order)) +``` + +Upwind DG operators in negative, central and positive operators can be obtained by `couple_discontinuously` + +```@example overview +central = couple_discontinuously(Dop, sbp_mesh) +minus = couple_discontinuously(Dop, sbp_mesh, Val(:minus)) +plus = couple_discontinuously(Dop, sbp_mesh, Val(:plus)) +D1 = PeriodicUpwindOperators(minus, central, plus) +``` + +In order to still have an entropy-conserving semidiscretization the second-derivative SBP operator needs to be + +```@example overview +using SparseArrays: sparse +D2 = sparse(plus) * sparse(minus) +``` + +The [`Solver`](@ref) object can now be created by passing the two SBP operators to the constructor, which, in turn, can be used to construct a `Semidiscretization`: + +```@example overview +solver = Solver(D1, D2) +semi = Semidiscretization(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) +``` + +As before, we can run the simulation by + +```@example overview +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + velocity, entropy)) +relaxation_callback = RelaxationCallback(invariant = entropy) +callbacks = CallbackSet(relaxation_callback, analysis_callback) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) +anim = @animate for step in 1:length(sol.u) + plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), yli = [(-0.8, 0.1)]) +end +gif(anim, "shoaling_solution_dg.gif", fps = 25) +``` + +For more details see also the [documentation of SummationByPartsOperators.jl](https://ranocha.de/SummationByPartsOperators.jl/stable/) + +## Additional resources + +Some more examples sorted by the simulated equations can be found in the [examples/](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/tree/main/examples) subdirectory. Especially, in [examples/svaerd\_kalisch\_1d/](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/tree/main/examples/svaerd_kalisch_1d) you can find Julia scripts that solve the [`SvaerdKalischEquations1D`](@ref) that were not covered in this tutorial. The same steps as described above, however, apply in the same way also to these equations. Attention must be taken for these equations because they do not conserve the classical total entropy ``\mathcal E``, but a modified entropy ``\hat{\mathcal E}``, available as [`entropy_modified`](@ref). + +More examples, especially focussing on plotting can be found in the script [create_figures.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/create_figures.jl). + +## References + +[^IsrawiKalischKatsaounisMitsotakis2021]: + Israwi, Kalisch, Katsaounis, Mitsotakis (2021). + A regularized shallow-water waves system with slip-wall boundary conditions in a basin: theory and numerical analysis. + [DOI: 10.1088/1361-6544/ac3c29](https://doi.org/10.1088/1361-6544/ac3c29) + +[^Ketcheson2019]: + Ketcheson (2019): + Relaxation Runge-Kutta Methods: Conservation and stability for Inner-Product Norms. + [DOI: 10.1137/19M1263662](https://doi.org/10.1137/19M1263662) + +[^RanochaSayyariDalcinParsaniKetcheson2020]: + Ranocha, Sayyari, Dalcin, Parsani, Ketcheson (2020): + Relaxation Runge–Kutta Methods: Fully-Discrete Explicit Entropy-Stable Schemes for the Compressible Euler and Navier–Stokes Equations + [DOI: 10.1137/19M1263480](https://doi.org/10.1137/19M1263480) + diff --git a/docs/src/reproduce.md b/docs/src/reproduce.md new file mode 100644 index 00000000..54ac9e17 --- /dev/null +++ b/docs/src/reproduce.md @@ -0,0 +1,12 @@ +# How to reproduce the figures + +In order to reproduce all figures used in the master thesis "Structure-preserving Numerical Methods for Dispersive Shallow Water Models" (2023) by Joshua Lampert execute the file located at the path [`DispersiveShallowWater.path_create_figures()`](@ref). From the Julia REPL, this can be done by: + +```julia +julia> using DispersiveShallowWater +julia> include(DispersiveShallowWater.path_create_figures()) +``` +Executing this script may take a while. It will generate a folder `out/` with certain subfolders containing the figures. If you want to modify the plots or only produce a subset of plots, you can download the script [`create_figures.jl`](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/create_figures.jl), modify accordingly and run it by +```julia +julia> include("create_figures.jl") +``` diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl index a73292ed..5409edbf 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl @@ -21,8 +21,8 @@ mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver accuracy_order = 4 -Dop = legendre_derivative_operator(mesh.xmin, mesh.xmax, accuracy_order) -sbp_mesh = UniformPeriodicMesh1D(-1.0, 1.0, div(mesh.N, accuracy_order)) +Dop = legendre_derivative_operator(-1.0, 1.0, accuracy_order) +sbp_mesh = UniformPeriodicMesh1D(mesh.xmin, mesh.xmax, div(mesh.N, accuracy_order)) central = couple_discontinuously(Dop, sbp_mesh) minus = couple_discontinuously(Dop, sbp_mesh, Val(:minus)) plus = couple_discontinuously(Dop, sbp_mesh, Val(:plus)) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index be3b1320..5c709ec8 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -216,7 +216,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, integrator, semi) println() println("─"^100) - println("Simulation running '", get_name(equations), "' with '", semi.initial_condition, "'") + println("Simulation running '", get_name(equations), "' with '", semi.initial_condition, + "'") println("─"^100) println(" #timesteps: " * @sprintf("% 14d", iter) * " " * diff --git a/src/callbacks_step/relaxation.jl b/src/callbacks_step/relaxation.jl index 489fa988..686025a6 100644 --- a/src/callbacks_step/relaxation.jl +++ b/src/callbacks_step/relaxation.jl @@ -6,7 +6,7 @@ Use a relaxation method in time in order to exactly preserve the (nonlinear) `invariant` is `invariant = entropy`. Reference -- Hendrik Ranocha, Mohammed Sayyari, Lisandro Dalcin, Matteo Parsani, David I. Ketcheson (2019) +- Hendrik Ranocha, Mohammed Sayyari, Lisandro Dalcin, Matteo Parsani, David I. Ketcheson (2020) Relaxation Runge–Kutta Methods: Fully-Discrete Explicit Entropy-Stable Schemes for the Compressible Euler and Navier–Stokes Equations [DOI: 10.1137/19M1263480](https://doi.org/10.1137/19M1263480) diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 0646c44f..c0411f5a 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -4,8 +4,8 @@ BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension. The equations are given by ```math \begin{aligned} - \frac{\partial\eta}{\partial t} + \frac{\partial}{\partial x}((\eta + D)v) - \frac{1}{6}D^2\frac{\partial}{\partial t}\frac{\partial^2}{\partial x^2}\eta &= 0,\\ - \frac{\partial v}{\partial t} + g\frac{\partial\eta}{\partial x} + \frac{\partial}{\partial x}\left(\frac{1}{2}v^2\right) - \frac{1}{6}D^2\frac{\partial}{\partial t}\frac{\partial^2}{\partial x^2}v &= 0. + \eta_t + ((\eta + D)v)_x - \frac{1}{6}D^2\eta_{xxt} &= 0,\\ + v_t + g\eta_x + \left(\frac{1}{2}v^2\right)_x - \frac{1}{6}D^2v_{xxt} &= 0. \end{aligned} ``` The unknown quantities of the BBM-BBM equations are the total water height ``\eta`` and the velocity ``v``. diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index 22958171..44e0587b 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -1,11 +1,11 @@ @doc raw""" - BBMBBMVariableEquations1D(gravity, eta0 = 0.0) + BBMBBMVariableEquations1D(gravity, eta0 = 1.0) BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension with spatially varying bathymetry. The equations are given by ```math \begin{aligned} - \frac{\partial\eta}{\partial t} + \frac{\partial}{\partial x}((\eta + D)v) - \frac{1}{6}\frac{\partial}{\partial x}(D^2\frac{\partial}{\partial t}\frac{\partial}{\partial x}\eta) &= 0,\\ - \frac{\partial v}{\partial t} + g\frac{\partial\eta}{\partial x} + \frac{\partial}{\partial x}\left(\frac{1}{2}v^2\right) - \frac{1}{6}\frac{\partial^2}{\partial x^2}(D^2\frac{\partial}{\partial t}v) &= 0. + \eta_t + ((\eta + D)v)_x - \frac{1}{6}(D^2\eta_{xt})_x &= 0,\\ + v_t + g\eta_x + \left(\frac{1}{2}v^2\right)_x - \frac{1}{6}(D^2v_t)_{xx} &= 0. \end{aligned} ``` The unknown quantities of the BBM-BBM equations are the total water height ``\eta`` and the velocity ``v``. diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 3fde07fd..22437d9a 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -1,20 +1,20 @@ @doc raw""" - SvärdKalischEquations1D(gravity, eta0 = 0.0, alpha = 0.0, beta = 0.2308939393939394, gamma = 0.04034343434343434) + SvärdKalischEquations1D(gravity, eta0 = 1.0, alpha = 0.0, beta = 0.2308939393939394, gamma = 0.04034343434343434) Dispersive system by Svärd and Kalisch in one spatial dimension with spatially varying bathymetry. The equations 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). + (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}. + 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 = -D``. The water height above the bathymetry is therefore given by @@ -193,50 +193,50 @@ end return SVector(eta, v, D) end -@inline function waterheight_total(u, equations::SvaerdKalischEquations1D) - return u[1] +@inline function waterheight_total(q, equations::SvaerdKalischEquations1D) + return q[1] end -@inline function velocity(u, equations::SvaerdKalischEquations1D) - return u[2] +@inline function velocity(q, equations::SvaerdKalischEquations1D) + return q[2] end -@inline function bathymetry(u, equations::SvaerdKalischEquations1D) - return -u[3] +@inline function bathymetry(q, equations::SvaerdKalischEquations1D) + return -q[3] end -@inline function waterheight(u, equations::SvaerdKalischEquations1D) - return waterheight_total(u, equations) - bathymetry(u, equations) +@inline function waterheight(q, equations::SvaerdKalischEquations1D) + return waterheight_total(q, equations) - bathymetry(q, equations) end -@inline function energy_total(u, equations::SvaerdKalischEquations1D) - eta, v, D = u +@inline function energy_total(q, equations::SvaerdKalischEquations1D) + eta, v, D = q e = 0.5 * (equations.gravity * eta^2 + (D + eta) * v^2) return e end @inline entropy(u, equations::SvaerdKalischEquations1D) = energy_total(u, equations) -# The modified entropy/total energy takes the whole `u` for every point in space +# The modified entropy/total energy takes the whole `q` for every point in space """ - energy_total_modified(u, equations::SvaerdKalischEquations1D, cache) + energy_total_modified(q, equations::SvaerdKalischEquations1D, cache) -Return the modified total energy of the conserved variables `u` for the +Return the modified total energy of the primitive variables `q` for the `SvaerdKalischEquations1D`. It contains an additional term containing a derivative compared to the usual `energy_total`. The `energy_total_modified` is a conserved quantity of the Svärd-Kalisch equations. -`u` is a vector of the conserved variables at ALL nodes, i.e., a matrix +`q` is a vector of the primitive variables at ALL nodes, i.e., a matrix of the correct length `nvariables(equations)` as first dimension and the number of nodes as length of the second dimension. `cache` needs to hold the first-derivative SBP operator `D1`. """ -@inline function energy_total_modified(u, equations::SvaerdKalischEquations1D, cache) - e_modified = zeros(eltype(u), size(u, 2)) +@inline function energy_total_modified(q, equations::SvaerdKalischEquations1D, cache) + e_modified = zeros(eltype(q), size(q, 2)) # Need to compute new beta_hat, do not use the old one from the `cache` - eta = view(u, 1, :) - v = view(u, 2, :) - D = view(u, 3, :) + eta = view(q, 1, :) + v = view(q, 2, :) + D = view(q, 3, :) beta_hat = equations.beta * (eta .+ D) .^ 3 if cache.D1 isa PeriodicDerivativeOperator tmp = 0.5 * beta_hat .* ((cache.D1 * v) .^ 2) @@ -245,21 +245,25 @@ number of nodes as length of the second dimension. else @error "unknown type of first-derivative operator" end - for i in 1:size(u, 2) - e_modified[i] = energy_total(view(u, :, i), equations) + tmp[i] + for i in 1:size(q, 2) + e_modified[i] = energy_total(view(q, :, i), equations) + tmp[i] end return e_modified end +varnames(::typeof(energy_total_modified), equations) = ("e_modified",) + """ - entropy_modified(u, equations::SvaerdKalischEquations1D, cache) + entropy_modified(q, equations::SvaerdKalischEquations1D, cache) Alias for [`energy_total_modified`](@ref). """ -@inline function entropy_modified(u, equations::SvaerdKalischEquations1D, cache) - energy_total_modified(u, equations, cache) +@inline function entropy_modified(q, equations::SvaerdKalischEquations1D, cache) + energy_total_modified(q, equations, cache) end +varnames(::typeof(entropy_modified), equations) = ("U_modified",) + # 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) diff --git a/src/util.jl b/src/util.jl index 4273a626..08e6314e 100644 --- a/src/util.jl +++ b/src/util.jl @@ -44,7 +44,24 @@ See also [`examples_dir`](@ref) and [`get_examples`](@ref). Copied from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl). """ function default_example() - joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d", "bbm_bbm_variable_bathymetry_1d_basic.jl") + joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d", + "bbm_bbm_variable_bathymetry_1d_basic.jl") +end + +""" + path_create_figures() + +Return the path to the file that creates all figures used in the master thesis "Structure-preserving +Numerical Methods for Dispersive Shallow Water Model" (2023). Executing this julia script may take a +while. + +# Examples +```@example +include(DispersiveShallowWater.path_create_figures()) +``` +""" +function path_create_figures() + pkgdir(DispersiveShallowWater, "create_figures.jl") end # Note: We can't call the method below `DispersiveShallowWater.include` since that is created automatically @@ -62,7 +79,8 @@ providing examples with sensible default values for users. Before replacing assignments in `example`, the keyword argument `maxiters` is inserted into calls to `solve` and `DispersiveShallowWater.solve` with it's default value used in the SciML ecosystem -for ODEs, see https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Miscellaneous. +for ODEs, see the "Miscellaneous" section of the +[documentation](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). Copied from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).