diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index c10f4bc8..526960d8 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -47,7 +47,7 @@ jobs: - '1.9' os: - ubuntu-latest - - macos-13 + - macos-latest - windows-latest arch: - x64 diff --git a/.github/workflows/DocPreviewCleanup.yml b/.github/workflows/DocPreviewCleanup.yml new file mode 100644 index 00000000..5be23b96 --- /dev/null +++ b/.github/workflows/DocPreviewCleanup.yml @@ -0,0 +1,33 @@ +name: Doc Preview Cleanup + +on: + pull_request: + types: [closed] + +# Ensure that only one "Doc Preview Cleanup" workflow is force pushing at a time +concurrency: + group: doc-preview-cleanup + cancel-in-progress: false + +jobs: + doc-preview-cleanup: + runs-on: ubuntu-latest + permissions: + contents: write + steps: + - name: Checkout gh-pages branch + uses: actions/checkout@v4 + with: + ref: gh-pages + - name: Delete preview and history + push changes + run: | + if [ -d "${preview_dir}" ]; then + git config user.name "Documenter.jl" + git config user.email "documenter@juliadocs.github.io" + git rm -rf "${preview_dir}" + git commit -m "delete preview" + git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree}) + git push --force origin gh-pages-new:gh-pages + fi + env: + preview_dir: previews/PR${{ github.event.number }} diff --git a/README.md b/README.md index 0cbbeee4..67e1e569 100644 --- a/README.md +++ b/README.md @@ -16,8 +16,8 @@ To date, it provides provably conservative, entropy-conserving and well-balanced * the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924), * the [Serre-Green-Naghdi equations](https://arxiv.org/abs/2408.02665). -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. +The semidiscretizations are based on summation-by-parts (SBP) operators, which are implemented in [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/). +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. A more detailed documentation can be found [online](https://JoshuaLampert.github.io./DispersiveShallowWater.jl/stable/). diff --git a/docs/make.jl b/docs/make.jl index 33b6133a..8503ee40 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,6 +2,31 @@ using Documenter using DispersiveShallowWater using TrixiBase +# Dynamically replace all files in subdirectories of the source directory to include all files in these subdirectories +# This way they don't need to be listed explicitly +EQUATIONS_FILES_TO_BE_INSERTED = joinpath.(Ref("equations"), + readdir(joinpath(dirname(@__DIR__), "src", + "equations"))) +CALLBACKS_STEP_FILES_TO_BE_INSERTED = joinpath.(Ref("callbacks_step"), + readdir(joinpath(dirname(@__DIR__), "src", + "callbacks_step"))) + +ref_path = joinpath(@__DIR__, "src", "ref.md") +lines = readlines(ref_path) +open(ref_path, "w") do io + for line in lines + if contains(line, "EQUATIONS_FILES_TO_BE_INSERTED") + line = replace(line, + "EQUATIONS_FILES_TO_BE_INSERTED" => EQUATIONS_FILES_TO_BE_INSERTED) + end + if contains(line, "CALLBACKS_STEP_FILES_TO_BE_INSERTED") + line = replace(line, + "CALLBACKS_STEP_FILES_TO_BE_INSERTED" => CALLBACKS_STEP_FILES_TO_BE_INSERTED) + end + println(io, line) + end +end + # Define module-wide setups such that the respective modules are available in doctests DocMeta.setdocmeta!(DispersiveShallowWater, :DocTestSetup, :(using DispersiveShallowWater); recursive = true) @@ -32,4 +57,5 @@ makedocs(; deploydocs(; repo = "github.com/JoshuaLampert/DispersiveShallowWater.jl", - devbranch = "main") + devbranch = "main", + push_preview = true) diff --git a/docs/src/development.md b/docs/src/development.md index 5e984c20..9299bf9d 100644 --- a/docs/src/development.md +++ b/docs/src/development.md @@ -48,3 +48,6 @@ julia --project=docs --color=yes docs/make.jl ``` The resulting `.html` files can then be found in `docs/build/` and you can look at them by opening them in a browser. +For pull requests from the main repository (i.e. not from a fork), the documentation is automatically built and can +be previewed under `https://joshualampert.github.io/DispersiveShallowWater.jl/previews/PRXXX/` where `XXX` is the number +of the pull request. diff --git a/docs/src/index.md b/docs/src/index.md index a66df7af..a695ec55 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -10,13 +10,14 @@ [![DOI](https://zenodo.org/badge/635090135.svg)](https://zenodo.org/doi/10.5281/zenodo.10034636) [**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 for two dispersive shallow water models: +To date, it provides provably conservative, entropy-conserving and well-balanced numerical schemes for some dispersive shallow water models: * the [BBM-BBM equations with varying bottom topography](https://iopscience.iop.org/article/10.1088/1361-6544/ac3c29), -* the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924). +* the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924), +* the [Serre-Green-Naghdi equations](https://arxiv.org/abs/2408.02665). -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. +The semidiscretizations are based on summation-by-parts (SBP) operators, which are implemented in [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/). +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 diff --git a/docs/src/overview.md b/docs/src/overview.md index fc14950d..30cb9eb2 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -2,7 +2,11 @@ ## 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 ``q = (\eta, v)`` they can be written as: +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 ``q = (\eta, v)`` they can be written as: ```math \begin{aligned} @@ -11,7 +15,10 @@ In this tutorial we describe how to numerically solve the BBM-BBM (Benjamin-Bona \end{aligned} ``` -Here, ``\eta = h + b`` describes the total water height, ``h`` the water height above the bottom topography (bathymetry), ``b = \eta_0 - D`` the bathymetry and ``v`` the velocity in horizontal direction. Here, ``\eta_0`` is a reference water height also called still water height. In the case of the BBM-BBM equations, ``\eta_0`` is usually taken to be 0. The gravitational acceleration is denoted as ``g``. A sketch of the water height and the bathymetry can be found below. +Here, ``\eta = h + b`` describes the total water height, ``h`` the water height above the bottom topography (bathymetry), ``b = \eta_0 - D`` +the bathymetry and ``v`` the velocity in horizontal direction. Here, ``\eta_0`` is a reference water height also called still water height. +In the case of the BBM-BBM equations, ``\eta_0`` is usually taken to be 0. The gravitational acceleration is denoted as ``g``. A sketch of +the water height and the bathymetry can be found below. ![water height and bathymetry](bathymetry.png) @@ -25,7 +32,9 @@ 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 the following example, the initial condition describes a traveling wave that moves towards a beach, which is modeled by a linearly increasing bathymetry. +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 the following example, the initial condition +describes a traveling wave that moves towards a beach, which is modeled by a linearly increasing bathymetry. ```@example overview equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) @@ -48,13 +57,26 @@ 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 gravitational 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 solution to be compared with the numerical solution. Otherwise, you can just keep the time variable unused. 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`. +The first line specifies that we want to solve the BBM-BBM equations with variable bathymetry using a gravitational 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 solution to be compared with the numerical solution. Otherwise, you can just keep the time variable unused. 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`. -Next, we choose periodic boundary conditions. DispersiveShallowWater.jl also supports reflecting boundary conditions for the [`BBMBBMEquations1D`](@ref) and [`BBMBBMVariableEquations1D`](@ref), see [`boundary_condition_reflecting`](@ref). 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-traveling waves. +Next, we choose periodic boundary conditions. DispersiveShallowWater.jl also supports reflecting boundary conditions for the [`BBMBBMEquations1D`](@ref) +and [`BBMBBMVariableEquations1D`](@ref), see [`boundary_condition_reflecting`](@ref). 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-traveling 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). Note that for non-periodic boundary conditions, the solver also needs to be created with non-periodic operators, see, e.g. [examples/bbm\_bbm\_1d/bbm\_bbm\_1d\_basic\_reflecting.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/examples/bbm_bbm_1d/bbm_bbm_1d_basic_reflecting.jl). +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 (SBP) 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). Note that for non-periodic boundary conditions, the solver also needs to be created with non-periodic +operators, see, e.g. [examples/bbm\_bbm\_1d/bbm\_bbm\_1d\_basic\_reflecting.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/examples/bbm_bbm_1d/bbm_bbm_1d_basic_reflecting.jl). ```@example overview solver = Solver(mesh, 4) @@ -66,7 +88,15 @@ Finally, we put the `mesh`, the `equations`, the `initial_condition`, the `solve ## 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`` error and ``L^\infty`` error of the different solution's variables compared to the initial condition (or, if available, at the same time analytical solution). Additional errors can be passed by the keyword argument `extra_analysis_errors`. 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 water height ``\eta`` [`waterheight_total`](@ref), the [`velocity`](@ref) and the [`entropy`](@ref) are computed and saved for each time step. The total water height and the total velocity are linear invariants of the BBM-BBM equations, i.e. they do not change over time. The total entropy +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`` error and ``L^\infty`` error of the different solution's variables compared to the initial condition (or, if available, +at the same time analytical solution). Additional errors can be passed by the keyword argument `extra_analysis_errors`. 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 water height ``\eta`` [`waterheight_total`](@ref), +the [`velocity`](@ref) and the [`entropy`](@ref) are computed and saved for each time step. The total water height 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 @@ -74,7 +104,9 @@ Once we have obtained a semidiscretization, we can solve the resulting system of 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). This means, we can specify a time-stepping scheme we want to use the tolerances for the adaptive time-stepping and the time values, where the solution values should be saved. In this case, we use the adaptive explicit Runge-Kutta method `Tsit5` by Tsitouras of order 5(4). Here, we save the solution at 100 equidistant points. +Finally, the `ode` can be `solve`d using the interface from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl). This means, we can specify a time-stepping +scheme we want to use the tolerances for the adaptive time-stepping and the time values, where the solution values should be saved. In this case, we use the adaptive +explicit Runge-Kutta method `Tsit5` by Tsitouras of order 5(4). Here, we save the solution at 100 equidistant points. ```@example overview tspan = (0.0, 25.0) @@ -91,11 +123,13 @@ 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). +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. +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 @@ -107,7 +141,12 @@ nothing # hide ![shoaling solution](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 values of the primitive variables `q` at one node, and the `equations` as input and should return an `SVector` of any length as output. For a user defined conversion function, 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 water height. The resulting plot will have one subplot for each of the returned variables of the conversion variable. By default, the conversion function is just [`prim2prim`](@ref), i.e. the identity. +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 values +of the primitive variables `q` at one node, and the `equations` as input and should return an `SVector` of any length as output. For a user defined conversion function, +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 water height. The resulting plot will have one subplot for each of the returned +variables of the conversion variable. By default, the conversion function is just [`prim2prim`](@ref), i.e. the identity. 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. @@ -121,8 +160,9 @@ nothing # hide ![shoaling solution](shoaling_solution.gif) -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 [plot_examples.jl](https://github.com/JoshuaLampert/2023-master-thesis/blob/main/code/plot_examples.jl) from -the reproducibility repository of the master thesis of Joshua Lampert for some examples. +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 +[plot_examples.jl](https://github.com/JoshuaLampert/2023-master-thesis/blob/main/code/plot_examples.jl) from the reproducibility repository of the master +thesis of Joshua Lampert 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 this end, you can `plot` the `AnalysisCallback` by @@ -136,11 +176,17 @@ This creates the following figure: ![analysis callback](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`` and rises 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. +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`` and rises 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 -To 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`: +To 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, @@ -154,7 +200,8 @@ 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. +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. @@ -169,13 +216,16 @@ nothing # hide ## [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: +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. +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) @@ -232,9 +282,13 @@ For more details see also the [documentation of SummationByPartsOperators.jl](ht ## 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 to these equations. Attention must be paid 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). +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 to these equations. +Attention must be paid 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 scripts [create_figures.jl](https://github.com/JoshuaLampert/2023-master-thesis/blob/main/code/create_figures.jl) and [plot_examples.jl](https://github.com/JoshuaLampert/2023-master-thesis/blob/main/code/plot_examples.jl) from the reproducibility repository of the master thesis of Joshua Lampert. +More examples, especially focussing on plotting, can be found in the scripts [create_figures.jl](https://github.com/JoshuaLampert/2023-master-thesis/blob/main/code/create_figures.jl) +and [plot_examples.jl](https://github.com/JoshuaLampert/2023-master-thesis/blob/main/code/plot_examples.jl) from the reproducibility repository of the master thesis of Joshua Lampert. ## References diff --git a/docs/src/ref.md b/docs/src/ref.md index ab8144e4..713fc9a4 100644 --- a/docs/src/ref.md +++ b/docs/src/ref.md @@ -6,4 +6,54 @@ CurrentModule = DispersiveShallowWater ```@autodocs Modules = [DispersiveShallowWater] +Pages = ["DispersiveShallowWater.jl"] +``` + +## Equations + +```@autodocs +Modules = [DispersiveShallowWater] +Pages = EQUATIONS_FILES_TO_BE_INSERTED +``` + +## Mesh + +```@autodocs +Modules = [DispersiveShallowWater] +Pages = ["mesh.jl"] +``` + +## Boundary conditions + +```@autodocs +Modules = [DispersiveShallowWater] +Pages = ["boundary_conditions.jl"] +``` + +## Solver + +```@autodocs +Modules = [DispersiveShallowWater] +Pages = ["solver.jl"] +``` + +## Semidiscretization + +```@autodocs +Modules = [DispersiveShallowWater] +Pages = ["semidiscretization.jl"] +``` + +## Callbacks + +```@autodocs +Modules = [DispersiveShallowWater] +Pages = CALLBACKS_STEP_FILES_TO_BE_INSERTED +``` + +## Utilities + +```@autodocs +Modules = [DispersiveShallowWater] +Pages = ["util.jl"] ``` diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 761ca562..294bb2e5 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -1,3 +1,18 @@ +""" + DispersiveShallowWater + +**DispersiveShallowWater.jl** is a Julia package that implements structure-preserving numerical methods for dispersive shallow water models. +It provides provably conservative, entropy-conserving, and well-balanced numerical schemes for some dispersive shallow water models. + +The semidiscretizations are based on summation-by-parts (SBP) operators, which are implemented in +[SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl). 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 provided by DispersiveShallowWater.jl. + +See also: [DispersiveShallowWater.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl) +""" module DispersiveShallowWater using BandedMatrices: BandedMatrix diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index 846a2ad5..e90ee896 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -161,7 +161,7 @@ Dingemans. The topography is a trapezoidal. !!! warning "Translation of water height" The initial condition for the water height is translated to be around 0, which is needed for the simulation because the `BBMBBMVariableEquations1D` are only implemented - for ``\eta_0 = 0``. + for ``\\eta_0 = 0``. References: - Magnus Svärd, Henrik Kalisch (2023) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ab4edced..3ceae561 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -42,7 +42,7 @@ Common choices of the `conversion_function` are [`prim2prim`](@ref) and """ function varnames end -""" +@doc raw""" AbstractShallowWaterEquations{NDIMS, NVARS} An abstract supertype of all equation system that contain the classical diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index f07dc6fd..0a5949ea 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -15,7 +15,7 @@ where ``\hat\alpha^2 = \alpha\sqrt{gD}D^2``, ``\hat\beta = \beta D^3``, ``\hat\g 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,\\ + \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} ``` @@ -337,9 +337,12 @@ end energy_total_modified(q_global, equations::SvaerdKalischEquations1D, cache) Return the modified total energy of the primitive variables `q_global` 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. +[`SvaerdKalischEquations1D`](@ref). It contains an additional term containing a +derivative compared to the usual [`energy_total`](@ref). The `energy_total_modified` +is a conserved quantity of the Svärd-Kalisch equations given by +```math +\\frac{1}{2} g h^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 first-derivative SBP operator `D1`. diff --git a/src/solver.jl b/src/solver.jl index a510b884..98f2b14d 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -8,7 +8,7 @@ abstract type AbstractSolver end """ Solver -A `struct` that holds the summation by parts (SBP) operators that are used for the spatial discretization. +A `struct` that holds the summation-by-parts (SBP) operators that are used for the spatial discretization. """ struct Solver{RealT <: Real, FirstDerivative <: AbstractDerivativeOperator{RealT}, @@ -36,7 +36,7 @@ end """ Solver(mesh, accuracy_order) -Create a solver, where the summation by parts (SBP) operators are of order `accuracy_order` and +Create a solver, where the summation-by-parts (SBP) operators are of order `accuracy_order` and associated to the `mesh`. """ function Solver(mesh, accuracy_order) diff --git a/test/test_bbm_bbm_1d.jl b/test/test_bbm_bbm_1d.jl index 236acafe..6fad0a9e 100644 --- a/test/test_bbm_bbm_1d.jl +++ b/test/test_bbm_bbm_1d.jl @@ -102,7 +102,9 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d") cons_error=[4.272316442782926e-11 0.5469460931577768], change_waterheight=4.272316442782926e-11, change_velocity=0.5469460931577768, - change_entropy=130.69415963528348) + change_entropy=130.69415963528348, + atol=1e-11, + atol_ints=1e-11) # in order to make CI pass @test_allocations(semi, sol)