diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 2819d0b1..d55df443 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -5,17 +5,25 @@ on: - main tags: ['*'] paths-ignore: + - 'CITATION.bib' - 'LICENSE.md' - 'README.md' + - '.zenodo.json' - '.github/workflows/CompatHelper.yml' + - '.github/workflows/Documenter.yml' + - '.github/workflows/Format-check.yml' - '.github/workflows/TagBot.yml' - '.github/workflows/SpellCheck.yml' - 'docs/**' pull_request: paths-ignore: + - 'CITATION.bib' - 'LICENSE.md' - 'README.md' + - '.zenodo.json' - '.github/workflows/CompatHelper.yml' + - '.github/workflows/Documenter.yml' + - '.github/workflows/Format-check.yml' - '.github/workflows/TagBot.yml' - '.github/workflows/SpellCheck.yml' - 'docs/**' diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index de582b4f..afb0d62f 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -6,6 +6,7 @@ on: - 'main' tags: '*' paths-ignore: + - '.zenodo.json' - '.github/workflows/CI.yml' - '.github/workflows/CompatHelper.yml' - '.github/workflows/TagBot.yml' diff --git a/.github/workflows/Format-check.yml b/.github/workflows/Format-check.yml index aa6e808e..240bdf53 100644 --- a/.github/workflows/Format-check.yml +++ b/.github/workflows/Format-check.yml @@ -25,7 +25,7 @@ jobs: - name: Install JuliaFormatter and format run: | julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))' - julia -e 'using JuliaFormatter; format(["src", "test", "examples", "visualization"], verbose = true)' + julia -e 'using JuliaFormatter; format(["src", "test", "examples"], verbose = true)' - name: Format check run: | julia -e ' diff --git a/CITATION.bib b/CITATION.bib new file mode 100644 index 00000000..e7f70104 --- /dev/null +++ b/CITATION.bib @@ -0,0 +1,9 @@ +@misc{lampert2023dispersive, + title={{D}ispersive{S}hallow{W}ater.jl: {S}tructure-preserving numerical + methods for dispersive shallow water models}, + author={Lampert, Joshua}, + year={2023}, + month={10}, + howpublished={\url{https://github.com/JoshuaLampert/DispersiveShallowWater.jl}}, + doi={10.5281/zenodo.10034636} +} diff --git a/Project.toml b/Project.toml index 4e18eea3..060118ea 100644 --- a/Project.toml +++ b/Project.toml @@ -21,12 +21,15 @@ SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" [compat] DiffEqBase = "6.121" Interpolations = "0.14" +LinearAlgebra = "1" PolynomialBases = "0.4.15" +Printf = "1" RecipesBase = "1.1" Reexport = "1.0" Roots = "2.0.17" SciMLBase = "1.90, 2" SimpleUnPack = "1.1" +SparseArrays = "1" StaticArrays = "1" SummationByPartsOperators = "0.5.41" julia = "1.8" diff --git a/README.md b/README.md index 3470042c..16e9e07d 100644 --- a/README.md +++ b/README.md @@ -56,9 +56,24 @@ for more details. Other examples can be found in the subdirectory [examples/](ht 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"))`. +## Referencing + +You can directly refer to DispersiveShallowWater.jl as +```bibtex +@misc{lampert2023dispersive, + title={{D}ispersive{S}hallow{W}ater.jl: {S}tructure-preserving numerical + methods for dispersive shallow water models}, + author={Lampert, Joshua}, + year={2023}, + month={10}, + howpublished={\url{https://github.com/JoshuaLampert/DispersiveShallowWater.jl}}, + doi={10.5281/zenodo.10034636} +} +``` + ## 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). +The package is developed and maintained by Joshua Lampert (University of Hamburg). 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/). diff --git a/docs/make.jl b/docs/make.jl index cd24ef9f..fe83b895 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -18,7 +18,6 @@ makedocs(; pages = [ "Home" => "index.md", "Overview" => "overview.md", - "Reproduce figures" => "reproduce.md", "Reference" => "ref.md", "License" => "license.md", ]) diff --git a/docs/src/index.md b/docs/src/index.md index 73e7dbf5..966d634a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,7 +18,7 @@ The semidiscretizations are based on summation by parts (SBP) operators, which a 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 +## Installation If you have not yet installed Julia, then 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. DispersiveShallowWater.jl is a registered Julia package. Therefore, you can install it by executing the following commands from the Julia REPL @@ -34,7 +34,7 @@ which can be installed running julia> Pkg.add("SummationByPartsOperators") ``` -# Usage +## Usage In the Julia REPL, first load the package DispersiveShallowWater.jl ```julia @@ -55,12 +55,27 @@ for more details. Other examples can be found in the subdirectory [examples/](ht 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 +## Referencing -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). +You can directly refer to DispersiveShallowWater.jl as +```bibtex +@misc{lampert2023dispersive, + title={{D}ispersive{S}hallow{W}ater.jl: {S}tructure-preserving numerical + methods for dispersive shallow water models}, + author={Lampert, Joshua}, + year={2023}, + month={10}, + howpublished={\url{https://github.com/JoshuaLampert/DispersiveShallowWater.jl}}, + doi={10.5281/zenodo.10034636} +} +``` + +## Authors + +The package is developed and maintained by Joshua Lampert (University of Hamburg). 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 +## 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/overview.md b/docs/src/overview.md index 00eac065..ef1d1a5b 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -117,7 +117,8 @@ 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 [plot_examples.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/visualization/plot_examples.jl) 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 @@ -224,7 +225,7 @@ For more details see also the [documentation of SummationByPartsOperators.jl](ht 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/DispersiveShallowWater.jl/blob/main/visualization/create_figures.jl) and [plot_examples.jl](https://github.com/JoshuaLampert/DispersiveShallowWater.jl/blob/main/visualization/plot_examples.jl). +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/reproduce.md b/docs/src/reproduce.md deleted file mode 100644 index 235be748..00000000 --- a/docs/src/reproduce.md +++ /dev/null @@ -1,19 +0,0 @@ -# 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()) -``` -Note that for one figure [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) is needed, so download Trixi.jl first: - -```julia -julia> using Pkg -julia> Pkg.add("Trixi") -``` -Executing the 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/visualization/create_figures.jl), modify it accordingly and run it by: - -```julia -julia> include("create_figures.jl") -``` diff --git a/visualization/bbm_bbm_1d_widestencil.jl b/examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl similarity index 73% rename from visualization/bbm_bbm_1d_widestencil.jl rename to examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl index 8dd057b0..c051cb3e 100644 --- a/visualization/bbm_bbm_1d_widestencil.jl +++ b/examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl @@ -1,36 +1,33 @@ using OrdinaryDiffEq using DispersiveShallowWater -using SummationByPartsOperators: periodic_derivative_operator -using SparseArrays: sparse ############################################################################### # Semidiscretization of the BBM-BBM equations -equations = BBMBBMEquations1D(gravity_constant = 9.81, D = 2.0) +equations = BBMBBMEquations1D(gravity_constant = 9.81, D = 4.0) -# initial_condition_convergence_test needs periodic boundary conditions -initial_condition = initial_condition_convergence_test +initial_condition = initial_condition_manufactured +source_terms = source_terms_manufactured boundary_conditions = boundary_condition_periodic # create homogeneous mesh -coordinates_min = -35.0 -coordinates_max = 35.0 +coordinates_min = 0.0 +coordinates_max = 1.0 N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 accuracy_order = 4 -D1 = periodic_derivative_operator(1, accuracy_order, mesh.xmin, mesh.xmax, mesh.N) -D2 = sparse(D1)^2 -solver = Solver(D1, D2) +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, - boundary_conditions = boundary_conditions) + boundary_conditions = boundary_conditions, + source_terms = source_terms) ############################################################################### # Create `ODEProblem` and run the simulation -tspan = (0.0, 30.0) +tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) analysis_callback = AnalysisCallback(semi; interval = 10, extra_analysis_errors = (:conservation_error,), diff --git a/visualization/bbm_bbm_variable_bathymetry_1d_widestencil.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl similarity index 76% rename from visualization/bbm_bbm_variable_bathymetry_1d_widestencil.jl rename to examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl index 8f62c419..12204888 100644 --- a/visualization/bbm_bbm_variable_bathymetry_1d_widestencil.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl @@ -1,36 +1,33 @@ using OrdinaryDiffEq using DispersiveShallowWater -using SummationByPartsOperators: periodic_derivative_operator -using SparseArrays: sparse ############################################################################### # Semidiscretization of the BBM-BBM equations equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) -# initial_condition_convergence_test needs periodic boundary conditions -initial_condition = initial_condition_convergence_test +initial_condition = initial_condition_manufactured +source_terms = source_terms_manufactured boundary_conditions = boundary_condition_periodic # create homogeneous mesh -coordinates_min = -35.0 -coordinates_max = 35.0 +coordinates_min = 0.0 +coordinates_max = 1.0 N = 512 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 accuracy_order = 4 -D1 = periodic_derivative_operator(1, accuracy_order, mesh.xmin, mesh.xmax, mesh.N) -D2 = sparse(D1)^2 -solver = Solver(D1, D2) +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, - boundary_conditions = boundary_conditions) + boundary_conditions = boundary_conditions, + source_terms = source_terms) ############################################################################### # Create `ODEProblem` and run the simulation -tspan = (0.0, 30.0) +tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) analysis_callback = AnalysisCallback(semi; interval = 10, extra_analysis_errors = (:conservation_error,), diff --git a/visualization/bbm_bbm_variable_bathymetry_1d_upwind.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl similarity index 64% rename from visualization/bbm_bbm_variable_bathymetry_1d_upwind.jl rename to examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl index 3ab09684..781d4bf3 100644 --- a/visualization/bbm_bbm_variable_bathymetry_1d_upwind.jl +++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_manufactured.jl @@ -1,38 +1,36 @@ using OrdinaryDiffEq using DispersiveShallowWater -using SummationByPartsOperators: upwind_operators, periodic_derivative_operator -using SparseArrays: sparse ############################################################################### -# Semidiscretization of the BBM-BBM equations +# Semidiscretization of the Svärd-Kalisch equations -equations = BBMBBMVariableEquations1D(gravity_constant = 9.81) +equations = SvaerdKalischEquations1D(gravity_constant = 1.0, eta0 = 0.0, + alpha = 0.0004040404040404049, + beta = 0.49292929292929294, + gamma = 0.15707070707070708) -# initial_condition_convergence_test needs periodic boundary conditions -initial_condition = initial_condition_convergence_test +initial_condition = initial_condition_manufactured +source_terms = source_terms_manufactured boundary_conditions = boundary_condition_periodic # create homogeneous mesh -coordinates_min = -35.0 -coordinates_max = 35.0 -N = 512 +coordinates_min = 0.0 +coordinates_max = 1.0 +N = 128 mesh = Mesh1D(coordinates_min, coordinates_max, N) # create solver with periodic SBP operators of accuracy order 4 accuracy_order = 4 -D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, - accuracy_order = accuracy_order, xmin = mesh.xmin, xmax = mesh.xmax, - N = mesh.N) -D2 = sparse(D1.plus) * sparse(D1.minus) -solver = Solver(D1, D2) +solver = Solver(mesh, accuracy_order) # semidiscretization holds all the necessary data structures for the spatial discretization semi = Semidiscretization(mesh, equations, initial_condition, solver, - boundary_conditions = boundary_conditions) + boundary_conditions = boundary_conditions, + source_terms = source_terms) ############################################################################### # Create `ODEProblem` and run the simulation -tspan = (0.0, 30.0) +tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) analysis_callback = AnalysisCallback(semi; interval = 10, extra_analysis_errors = (:conservation_error,), diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 4d74c399..77b98726 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -24,8 +24,8 @@ import SummationByPartsOperators: grid, xmin, xmax include("boundary_conditions.jl") include("mesh.jl") -include("solver.jl") include("equations/equations.jl") +include("solver.jl") include("semidiscretization.jl") include("callbacks_step/callbacks_step.jl") include("visualization.jl") @@ -48,7 +48,9 @@ export Semidiscretization, semidiscretize, grid export boundary_condition_periodic -export initial_condition_convergence_test, initial_condition_dingemans +export initial_condition_convergence_test, + initial_condition_manufactured, source_terms_manufactured, + initial_condition_dingemans export AnalysisCallback, RelaxationCallback export tstops, errors, integrals diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 8e8046b5..fa5844aa 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -53,6 +53,41 @@ function initial_condition_convergence_test(x, t, equations::BBMBBMEquations1D, return SVector(eta, v) end +""" + initial_condition_manufactured(x, t, equations::BBMBBMEquations1D, mesh) + +A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref). +""" +function initial_condition_manufactured(x, t, + equations::BBMBBMEquations1D, + mesh) + eta = exp(t) * cospi(2 * (x - 2 * t)) + v = exp(t / 2) * sinpi(2 * (x - t / 2)) + return SVector(eta, v) +end + +""" + source_terms_manufactured(q, x, t, equations::BBMBBMEquations1D, mesh) + +A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref). +""" +function source_terms_manufactured(q, x, t, equations::BBMBBMEquations1D) + g = equations.gravity + D = equations.D + a3 = cospi(t - 2 * x) + a4 = sinpi(t - 2 * x) + a5 = sinpi(2 * t - 4 * x) + a6 = sinpi(4 * t - 2 * x) + a7 = cospi(4 * t - 2 * x) + dq1 = -2 * pi^2 * D^2 * (4 * pi * a6 - a7) * exp(t) / 3 + 2 * pi * D * exp(t / 2) * a3 - + 2 * pi * exp(3 * t / 2) * a4 * a6 + 2 * pi * exp(3 * t / 2) * a3 * a7 - + 4 * pi * exp(t) * a6 + exp(t) * a7 + dq2 = -pi^2 * D^2 * (a4 + 2 * pi * a3) * exp(t / 2) / 3 + 2 * pi * g * exp(t) * a6 - + exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5 + + return SVector(dq1, dq2) +end +# function create_cache(mesh, equations::BBMBBMEquations1D, solver, @@ -61,15 +96,14 @@ function create_cache(mesh, uEltype) if solver.D1 isa PeriodicDerivativeOperator || solver.D1 isa UniformPeriodicCoupledOperator - invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ Matrix(solver.D1) + invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2)) elseif solver.D1 isa PeriodicUpwindOperators - invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ - Matrix(solver.D1.central) + invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2)) else @error "unknown type of first-derivative operator" end - tmp1 = Array{RealT}(undef, nnodes(mesh)) - return (invImD2_D = invImD2_D, tmp1 = tmp1) + tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback` + return (invImD2 = invImD2, tmp1 = tmp1) end # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see @@ -77,8 +111,8 @@ end # 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!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_condition, - ::BoundaryConditionPeriodic, solver, cache) - @unpack invImD2_D, tmp1 = cache + ::BoundaryConditionPeriodic, source_terms, solver, cache) + @unpack invImD2 = cache q = wrap_array(u_ode, mesh, equations, solver) dq = wrap_array(du_ode, mesh, equations, solver) @@ -89,11 +123,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_cond dv = view(dq, 2, :) # energy and mass conservative semidiscretization - @. tmp1 = -(equations.D * v + eta * v) - mul!(deta, invImD2_D, tmp1) + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator + deta[:] = -solver.D1 * (equations.D * v + eta .* v) + dv[:] = -solver.D1 * (equations.gravity * eta + 0.5 * v .^ 2) + elseif solver.D1 isa PeriodicUpwindOperators + deta[:] = -solver.D1.central * (equations.D * v + eta .* v) + dv[:] = -solver.D1.central * (equations.gravity * eta + 0.5 * v .^ 2) + else + @error "unknown type of first-derivative operator" + end + + calc_sources!(dq, q, t, source_terms, equations, solver) - @. tmp1 = -(equations.gravity * eta + 0.5 * v^2) - mul!(dv, invImD2_D, tmp1) + deta[:] = invImD2 * deta + dv[:] = invImD2 * dv return nothing end diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index d69c10bf..ad63e6a0 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -41,8 +41,7 @@ For details see Example 5 in Section 3 from (here adapted for dimensional equati Exact Traveling-Wave Solutions to Bidirectional Wave Equations [DOI: 10.1023/A:1026667903256](https://doi.org/10.1023/A:1026667903256) """ -function initial_condition_convergence_test(x, - t, +function initial_condition_convergence_test(x, t, equations::BBMBBMVariableEquations1D, mesh) g = equations.gravity @@ -58,12 +57,55 @@ function initial_condition_convergence_test(x, return SVector(eta, v, D) end +""" + initial_condition_manufactured(x, t, equations::BBMBBMVariableEquations1D, mesh) + +A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref). +""" +function initial_condition_manufactured(x, t, + equations::BBMBBMVariableEquations1D, + mesh) + eta = exp(t) * cospi(2 * (x - 2 * t)) + v = exp(t / 2) * sinpi(2 * (x - t / 2)) + D = 5.0 + 2.0 * cospi(2 * x) + return SVector(eta, v, D) +end + +""" + source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D, mesh) + +A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref). +""" +function source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D) + g = equations.gravity + a1 = cospi(2 * x) + a2 = sinpi(2 * x) + a3 = cospi(t - 2 * x) + a4 = sinpi(t - 2 * x) + a5 = sinpi(2 * t - 4 * x) + a6 = sinpi(4 * t - 2 * x) + a7 = cospi(4 * t - 2 * x) + dq1 = -2 * pi^2 * (4 * pi * a6 - a7) * (2 * a1 + 5)^2 * exp(t) / 3 + + 8 * pi^2 * (a6 + 4 * pi * a7) * (2 * a1 + 5) * exp(t) * a2 / 3 + + 2 * pi * (2 * a1 + 5) * exp(t / 2) * a3 - 2 * pi * exp(3 * t / 2) * a4 * a6 + + 2 * pi * exp(3 * t / 2) * a3 * a7 + 4 * pi * exp(t / 2) * a2 * a4 - + 4 * pi * exp(t) * a6 + exp(t) * a7 + dq2 = 2 * pi * g * exp(t) * a6 - + pi^2 * + (8 * (2 * pi * a4 - a3) * (2 * a1 + 5) * a2 + + (a4 + 2 * pi * a3) * (2 * a1 + 5)^2 + + 4 * (a4 + 2 * pi * a3) * (16 * sinpi(x)^4 - 26 * sinpi(x)^2 + 7)) * exp(t / 2) / + 3 - exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5 + + return SVector(dq1, dq2, 0.0) +end + """ initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh) The initial condition that uses the dispersion relation of the Euler equations to approximate waves generated by a wave maker as it is done by experiments of -Dingemans. The topography is a trapesoidal. +Dingemans. The topography is a trapezoidal. References: - Magnus Svärd, Henrik Kalisch (2023) @@ -84,16 +126,14 @@ function initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, h = A * cos(k * x) end v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0 - if x < 11.01 || x >= 33.07 - b = 0.0 - elseif 11.01 <= x && x < 23.04 + if 11.01 <= x && x < 23.04 b = 0.6 * (x - 11.01) / (23.04 - 11.01) elseif 23.04 <= x && x < 27.04 b = 0.6 elseif 27.04 <= x && x < 33.07 b = 0.6 * (33.07 - x) / (33.07 - 27.04) else - error("should not happen") + b = 0.0 end eta = h + eta0 D = -b @@ -102,7 +142,7 @@ end function create_cache(mesh, equations::BBMBBMVariableEquations1D, - solver::Solver, + solver, initial_condition, RealT, uEltype) @@ -115,18 +155,16 @@ function create_cache(mesh, K = Diagonal(D .^ 2) if solver.D1 isa PeriodicDerivativeOperator || solver.D1 isa UniformPeriodicCoupledOperator - invImDKD_D = (I - 1 / 6 * sparse(solver.D1) * K * sparse(solver.D1)) \ - Matrix(solver.D1) - invImD2K_D = (I - 1 / 6 * sparse(solver.D2) * K) \ Matrix(solver.D1) + invImDKD = inv(I - 1 / 6 * Matrix(solver.D1) * K * Matrix(solver.D1)) + invImD2K = inv(I - 1 / 6 * Matrix(solver.D2) * K) elseif solver.D1 isa PeriodicUpwindOperators - invImDKD_D = (I - 1 / 6 * sparse(solver.D1.minus) * K * sparse(solver.D1.plus)) \ - Matrix(solver.D1.minus) - invImD2K_D = (I - 1 / 6 * sparse(solver.D2) * K) \ Matrix(solver.D1.plus) + invImDKD = inv(I - 1 / 6 * Matrix(solver.D1.minus) * K * Matrix(solver.D1.plus)) + invImD2K = inv(I - 1 / 6 * Matrix(solver.D2) * K) else @error "unknown type of first-derivative operator" end - tmp1 = Array{RealT}(undef, nnodes(mesh)) - return (invImDKD_D = invImDKD_D, invImD2K_D = invImD2K_D, tmp1 = tmp1) + tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback` + return (invImDKD = invImDKD, invImD2K = invImD2K, tmp1 = tmp1) end # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see @@ -135,9 +173,9 @@ end # [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119) # Here, adapted for spatially varying bathymetry. function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D, - initial_condition, - ::BoundaryConditionPeriodic, solver, cache) - @unpack invImDKD_D, invImD2K_D, tmp1 = cache + initial_condition, ::BoundaryConditionPeriodic, source_terms, + solver, cache) + @unpack invImDKD, invImD2K = cache q = wrap_array(u_ode, mesh, equations, solver) dq = wrap_array(du_ode, mesh, equations, solver) @@ -150,11 +188,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D, dD = view(dq, 3, :) fill!(dD, zero(eltype(dD))) - @. tmp1 = -(D * v + eta * v) - mul!(deta, invImDKD_D, tmp1) + if solver.D1 isa PeriodicDerivativeOperator || + solver.D1 isa UniformPeriodicCoupledOperator + deta[:] = -solver.D1 * (D .* v + eta .* v) + dv[:] = -solver.D1 * (equations.gravity * eta + 0.5 * v .^ 2) + elseif solver.D1 isa PeriodicUpwindOperators + deta[:] = -solver.D1.minus * (D .* v + eta .* v) + dv[:] = -solver.D1.plus * (equations.gravity * eta + 0.5 * v .^ 2) + else + @error "unknown type of first-derivative operator" + end + + calc_sources!(dq, q, t, source_terms, equations, solver) - @. tmp1 = -(equations.gravity * eta + 0.5 * v^2) - mul!(dv, invImD2K_D, tmp1) + deta[:] = invImDKD * deta + dv[:] = invImD2K * dv return nothing end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 3a53f376..cf15d74e 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -49,7 +49,7 @@ varnames(::typeof(prim2cons), ::SvaerdKalischEquations1D) = ("h", "hv", "b") The initial condition that uses the dispersion relation of the Euler equations to approximate waves generated by a wave maker as it is done by experiments of -Dingemans. The topography is a trapesoidal. It is assumed that `equations.eta0 = 0.8`. +Dingemans. The topography is a trapezoidal. It is assumed that `equations.eta0 = 0.8`. References: - Magnus Svärd, Henrik Kalisch (2023) @@ -70,25 +70,69 @@ function initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, h = A * cos(k * x) end v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0 - if x < 11.01 || x >= 33.07 - b = 0.0 - elseif 11.01 <= x && x < 23.04 + if 11.01 <= x && x < 23.04 b = 0.6 * (x - 11.01) / (23.04 - 11.01) elseif 23.04 <= x && x < 27.04 b = 0.6 elseif 27.04 <= x && x < 33.07 b = 0.6 * (33.07 - x) / (33.07 - 27.04) else - error("should not happen") + b = 0.0 end eta = h + eta0 D = -b return SVector(eta, v, D) end +""" + initial_condition_manufactured(x, t, equations::SvaerdKalischEquations1D, mesh) + +A smooth manufactured solution in combination with [`source_terms_manufactured`](@ref). +""" +function initial_condition_manufactured(x, t, + equations::SvaerdKalischEquations1D, + mesh) + eta = exp(t) * cospi(2 * (x - 2 * t)) + v = exp(t / 2) * sinpi(2 * (x - t / 2)) + D = 3.0 + return SVector(eta, v, D) +end + +""" + source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D, mesh) + +A smooth manufactured solution in combination with [`initial_condition_manufactured`](@ref). +""" +function source_terms_manufactured(q, x, t, equations::SvaerdKalischEquations1D) + g = equations.gravity + D = q[3, 1] # D is constant, thus simply take the first entry + eta0 = equations.eta0 + alpha = equations.alpha + beta = equations.beta + gamma = equations.gamma + a1 = cospi(t - 2 * x) + a2 = sinpi(t - 2 * x) + a3 = cospi(4 * t - 2 * x) + a4 = sinpi(4 * t - 2 * x) + + dq1 = 8 * pi^3 * alpha * sqrt(g * (D + eta0)) * (D + eta0)^2 * exp(t) * a4 + + 2 * pi * (D + exp(t) * a3) * exp(t / 2) * a1 - 2 * pi * exp(3 * t / 2) * a2 * a4 - + 4 * pi * exp(t) * a4 + exp(t) * a3 + dq2 = 2 * pi * D * g * exp(t) * a4 - D * exp(t / 2) * a2 / 2 - + pi * D * exp(t / 2) * a1 - 2 * pi * D * exp(t) * a2 * a1 + + 8 * pi^3 * alpha * (D + eta0)^2 * sqrt(D * g + eta0 * g) * exp(3 * t / 2) * a1 * + a3 - 2 * pi^2 * beta * (D + eta0)^3 * exp(t / 2) * a2 - + 4 * pi^3 * beta * (D + eta0)^3 * exp(t / 2) * a1 + + 2 * pi * g * exp(2 * t) * a4 * a3 + + 8.0 * pi^3 * gamma * (D + eta0)^3 * sqrt(D * g + eta0 * g) * exp(t / 2) * a1 - + exp(3 * t / 2) * a2 * a3 / 2 - pi * exp(3 * t / 2) * a1 * a3 - + 2 * pi * exp(2 * t) * a2 * a1 * a3 + return SVector(dq1, dq2, 0.0) +end +# function create_cache(mesh, equations::SvaerdKalischEquations1D, - solver::Solver, + solver, initial_condition, RealT, uEltype) @@ -125,8 +169,8 @@ end # Discretization that conserves the mass (for eta and v) and is energy-bounded for periodic boundary conditions function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, - initial_condition, - ::BoundaryConditionPeriodic, solver, cache) + initial_condition, ::BoundaryConditionPeriodic, source_terms, + solver, cache) @unpack hmD1betaD1, D1betaD1, d, h, hv, alpha_hat, beta_hat, gamma_hat, tmp1, tmp2, D1_central = cache q = wrap_array(u_ode, mesh, equations, solver) dq = wrap_array(du_ode, mesh, equations, solver) @@ -139,7 +183,7 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, dD = view(dq, 3, :) fill!(dD, zero(eltype(dD))) - @. h = eta + D + h = eta .+ D hv = h .* v if solver.D1 isa PeriodicDerivativeOperator || @@ -166,18 +210,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D, hmD1betaD1 = Diagonal(h) - D1betaD1 # split form - tmp2 = -(0.5 * (D1_central * (hv .* v) + hv .* D1v - v .* (D1_central * hv)) + - equations.gravity * h .* D1eta + - 0.5 * (vD1y - D1vy - yD1v) - - 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - - 0.5 * solver.D2 * (gamma_hat .* D1v)) - # not split form - # tmp2 = -(D1_central * (hv .* v) - v .* (D1_central * hv)+ - # equations.gravity * h .* D1eta + - # vD1y - D1vy - - # 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - - # 0.5 * solver.D2 * (gamma_hat .* D1v)) - dv[:] = hmD1betaD1 \ tmp2 + dv[:] = -(0.5 * (D1_central * (hv .* v) + hv .* D1v - v .* (D1_central * hv)) + + equations.gravity * h .* D1eta + + 0.5 * (vD1y - D1vy - yD1v) - + 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - + 0.5 * solver.D2 * (gamma_hat .* D1v)) + + # no split form + # dv[:] = -(D1_central * (hv .* v) - v .* (D1_central * hv)+ + # equations.gravity * h .* D1eta + + # vD1y - D1vy - + # 0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) - + # 0.5 * solver.D2 * (gamma_hat .* D1v)) + + calc_sources!(dq, q, t, source_terms, equations, solver) + dv[:] = hmD1betaD1 \ dv return nothing end diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl index 05aadf7f..eda6bdb5 100644 --- a/src/semidiscretization.jl +++ b/src/semidiscretization.jl @@ -6,7 +6,7 @@ A `struct` containing everything needed to describe a spatial semidiscretization of an equation. """ struct Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions, - Solver, Cache} + SourceTerms, Solver, Cache} mesh::Mesh equations::Equations @@ -15,30 +15,34 @@ struct Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions, initial_condition::InitialCondition boundary_conditions::BoundaryConditions + source_terms::SourceTerms solver::Solver cache::Cache function Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions, - Solver, + SourceTerms, Solver, Cache}(mesh::Mesh, equations::Equations, initial_condition::InitialCondition, boundary_conditions::BoundaryConditions, + source_terms::SourceTerms, solver::Solver, cache::Cache) where {Mesh, Equations, InitialCondition, - BoundaryConditions, Solver, - Cache} + BoundaryConditions, SourceTerms, + Solver, Cache} @assert ndims(mesh) == ndims(equations) @assert xmin(mesh) == xmin(solver.D1) @assert xmax(mesh) == xmax(solver.D1) @assert nnodes(mesh) == length(grid(solver)) - new(mesh, equations, initial_condition, boundary_conditions, solver, cache) + new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver, + cache) end end """ Semidiscretization(mesh, equations, initial_condition, solver; + source_terms=nothing, boundary_conditions=boundary_condition_periodic, RealT=real(solver), uEltype=RealT, @@ -47,6 +51,7 @@ end Construct a semidiscretization of a PDE. """ function Semidiscretization(mesh, equations, initial_condition, solver; + source_terms = nothing, boundary_conditions = boundary_condition_periodic, # `RealT` is used as real type for node locations etc. # while `uEltype` is used as element type of solutions etc. @@ -57,12 +62,10 @@ function Semidiscretization(mesh, equations, initial_condition, solver; initial_cache...) Semidiscretization{typeof(mesh), typeof(equations), typeof(initial_condition), - typeof(boundary_conditions), typeof(solver), typeof(cache)}(mesh, - equations, - initial_condition, - boundary_conditions, - solver, - cache) + typeof(boundary_conditions), typeof(source_terms), + typeof(solver), typeof(cache)}(mesh, equations, initial_condition, + boundary_conditions, source_terms, + solver, cache) end function Base.show(io::IO, semi::Semidiscretization) @@ -73,6 +76,7 @@ function Base.show(io::IO, semi::Semidiscretization) print(io, ", ", semi.equations) print(io, ", ", semi.initial_condition) print(io, ", ", semi.boundary_conditions) + print(io, ", ", semi.source_terms) print(io, ", ", semi.solver) print(io, ", cache(") for (idx, key) in enumerate(keys(semi.cache)) @@ -93,7 +97,8 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization) println(io, " mesh: ", semi.mesh) println(io, " equations: ", get_name(semi.equations)) println(io, " initial condition: ", semi.initial_condition) - print(io, " boundary condition: ", semi.boundary_conditions) + println(io, " boundary condition: ", semi.boundary_conditions) + print(io, " source terms: ", semi.source_terms) end end @@ -196,10 +201,11 @@ function wrap_array(u_ode, semi::Semidiscretization) end function rhs!(du_ode, u_ode, semi::Semidiscretization, t) - @unpack mesh, equations, initial_condition, boundary_conditions, solver, cache = semi + @unpack mesh, equations, initial_condition, boundary_conditions, solver, source_terms, cache = semi - rhs!(du_ode, u_ode, t, mesh, equations, initial_condition, boundary_conditions, solver, - cache) + rhs!(du_ode, u_ode, t, mesh, equations, initial_condition, boundary_conditions, + source_terms, + solver, cache) return nothing end diff --git a/src/solver.jl b/src/solver.jl index 06e8a6d8..1bc2d645 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -114,3 +114,20 @@ function calc_error_norms(u_ode, t, initial_condition, mesh::Mesh1D, equations, end return l2_error, linf_error end + +function calc_sources!(dq, q, t, source_terms::Nothing, + equations::AbstractEquations{1}, solver::Solver) + return nothing +end + +function calc_sources!(dq, q, t, source_terms, + equations::AbstractEquations{1}, solver::Solver) + x = grid(solver) + for i in eachnode(solver) + local_source = source_terms(view(q, :, i), x[i], t, equations) + for v in eachvariable(equations) + dq[v, i] += local_source[v] + end + end + return nothing +end diff --git a/src/util.jl b/src/util.jl index df7d7378..c15072af 100644 --- a/src/util.jl +++ b/src/util.jl @@ -48,22 +48,6 @@ function default_example() "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, "visualization", "create_figures.jl") -end - # Note: We can't call the method below `DispersiveShallowWater.include` since that is created automatically # inside `module DispersiveShallowWater` to `include` source files and evaluate them within the global scope # of `DispersiveShallowWater`. However, users will want to evaluate in the global scope of `Main` or something @@ -96,6 +80,16 @@ julia> redirect_stdout(devnull) do ``` """ function trixi_include(mod::Module, example::AbstractString; kwargs...) + # Check that all kwargs exist as assignments + code = read(example, String) + expr = Meta.parse("begin \n$code \nend") + expr = insert_maxiters(expr) + + for (key, val) in kwargs + # This will throw an error when `key` is not found + find_assignment(expr, key) + end + Base.include(ex -> replace_assignments(insert_maxiters(ex); kwargs...), mod, example) end @@ -282,6 +276,7 @@ end function find_assignment(expr, destination) # declare result to be able to assign to it in the closure local result + found = false # find explicit and keyword assignments walkexpr(expr) do x @@ -289,12 +284,17 @@ function find_assignment(expr, destination) if (x.head === Symbol("=") || x.head === :kw) && x.args[1] === Symbol(destination) result = x.args[2] + found = true # dump(x) end end return x end + if !found + throw(ArgumentError("assignment `$destination` not found in expression")) + end + result end diff --git a/test/Project.toml b/test/Project.toml index d699460a..cdb84b7a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -7,7 +7,9 @@ SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -Aqua = "0.7" +Aqua = "0.7, 0.8" OrdinaryDiffEq = "6.49.1" Plots = "1.16" +SparseArrays = "1" SummationByPartsOperators = "0.5.41" +Test = "1" diff --git a/test/test_bbm_bbm_1d.jl b/test/test_bbm_bbm_1d.jl index 22236e2c..f7a72458 100644 --- a/test/test_bbm_bbm_1d.jl +++ b/test/test_bbm_bbm_1d.jl @@ -41,6 +41,19 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d") change_velocity=-4.547473508864641e-13, change_entropy=0.0) end + + @trixi_testset "bbm_bbm_1d_manufactured" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_1d_manufactured.jl"), + tspan=(0.0, 1.0), + l2=[4.365176233405813e-9 6.7151849982388e-10], + linf=[6.226559268185383e-9 9.698699621196738e-10], + cons_error=[3.873483998828586e-12 2.2986355942039745e-11], + change_waterheight=3.873483998828586e-12, + change_velocity=2.2986355942039745e-11, + change_entropy=17.387441847193436, + atol=1e-10, + atol_ints=1e-10) # in order to make CI pass + end end end # module diff --git a/test/test_bbm_bbm_variable_bathymetry_1d.jl b/test/test_bbm_bbm_variable_bathymetry_1d.jl index 8bcea168..7910ab09 100644 --- a/test/test_bbm_bbm_variable_bathymetry_1d.jl +++ b/test/test_bbm_bbm_variable_bathymetry_1d.jl @@ -54,7 +54,22 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d") cons_error=[7.194245199571014e-14 1.2388735870505485e-12 0.0], change_waterheight=-1.5765166949677223e-13, change_velocity=1.8791999206735355e-13, - change_entropy=2.1316282072803006e-14) + change_entropy=2.1316282072803006e-14, + atol=1e-11) # in order to make CI pass) + end + + @trixi_testset "bbm_bbm_variable_bathymetry_1d_manufactured" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "bbm_bbm_variable_bathymetry_1d_manufactured.jl"), + tspan=(0.0, 1.0), + l2=[6.867287425380051e-9 3.446178245128195e-9 0.0], + linf=[1.1667709465257303e-8 5.917459633408839e-9 0.0], + cons_error=[1.359075785939412e-11 3.8711139735371144e-13 0.0], + change_waterheight=-1.359075785939412e-11, + change_velocity=-3.8711139735371144e-13, + change_entropy=17.81701226932122, + atol=1e-10, + atol_ints=1e-11) # in order to make CI pass end @trixi_testset "bbm_bbm_variable_bathymetry_1d_well_balanced" begin diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl index 935e22b4..96bc6969 100644 --- a/test/test_svaerd_kalisch_1d.jl +++ b/test/test_svaerd_kalisch_1d.jl @@ -8,6 +8,18 @@ include("test_util.jl") EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d") @testset "SvaerdKalisch1D" begin + @trixi_testset "svaerd_kalisch_1d_manufactured" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "svaerd_kalisch_1d_manufactured.jl"), + tspan=(0.0, 0.1), + l2=[7.467887060263923e-5 2.7796353838948894e-8 0.0], + linf=[0.0001613144395267163 4.344495230235168e-8 0.0], + cons_error=[2.3635607360183997e-16 8.084235123776567e-10 0.0], + change_waterheight=-2.3635607360183997e-16, + change_entropy=0.1342289500320556, + atol=1e-4) # in order to make CI pass + end + @trixi_testset "svaerd_kalisch_1d_dingemans" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "svaerd_kalisch_1d_dingemans.jl"), diff --git a/test/test_unit.jl b/test/test_unit.jl index ebe64afc..90fa26a7 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -216,15 +216,13 @@ using SparseArrays: sparse, SparseMatrixCSC @testset "util" begin @test_nowarn get_examples() - @test_nowarn DispersiveShallowWater.path_create_figures() @test_nowarn trixi_include(default_example(), tspan = (0.0, 0.1)) accuracy_orders = [2, 4, 6] for accuracy_order in accuracy_orders - eoc_mean_values, _ = convergence_test(default_example(), 2, initial_N = 128, + eoc_mean_values, _ = convergence_test(default_example(), 2, N = 512, tspan = (0.0, 1.0), accuracy_order = accuracy_order) - show(eoc_mean_values) @test isapprox(eoc_mean_values[:l2][1], accuracy_order, atol = 0.5) @test isapprox(eoc_mean_values[:linf][2], accuracy_order, atol = 0.5) @test isapprox(eoc_mean_values[:l2][1], accuracy_order, atol = 0.5) diff --git a/test/test_util.jl b/test/test_util.jl index 8a356847..81027d90 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -41,8 +41,8 @@ macro test_trixi_include(example, args...) if (arg.head == :(=) && !(arg.args[1] in (:l2, :linf, :cons_error, :change_waterheight, :change_velocity, :change_momentum, :change_entropy, - :change_entropy_modified, :atol, :rtol, :atol_ints, - :rtol_ints))) + :change_entropy_modified, :lake_at_rest, + :atol, :rtol, :atol_ints, :rtol_ints))) push!(kwargs, Pair(arg.args...)) end end diff --git a/visualization/create_figures.jl b/visualization/create_figures.jl deleted file mode 100644 index 873cc933..00000000 --- a/visualization/create_figures.jl +++ /dev/null @@ -1,735 +0,0 @@ -using DispersiveShallowWater -using SummationByPartsOperators: SummationByPartsOperators, - periodic_derivative_operator, - upwind_operators, - legendre_derivative_operator, - legendre_second_derivative_operator, - UniformPeriodicMesh1D, - couple_discontinuously, - couple_continuously -using Trixi: PlotData1D -using SparseArrays: sparse -using Plots -using LaTeXStrings - -const OUT = "out/" -ispath(OUT) || mkpath(OUT) -const VISUALIZATION_DIR = pkgdir(DispersiveShallowWater, "visualization") -const EXAMPLES_DIR_BBMBBM = joinpath(examples_dir(), "bbm_bbm_1d") -const EXAMPLES_DIR_BBMBBM_VARIABLE = joinpath(examples_dir(), - "bbm_bbm_variable_bathymetry_1d") -const EXAMPLES_DIR_SVAERD_KALISCH = joinpath(examples_dir(), "svaerd_kalisch_1d") - -# Chapter 2 -# Plot of bathymetry and waterheight -function fig_1() - L = 1.0 - n = 100 - x = LinRange(0.0, L, n) - fontsize = 20 - - # just pick some function for b and eta that look nice - H = 1.012 - - b(x) = x * cos.(3 * pi * x) + H - plot(x, b, color = :gray, fill = (0, 0.8, :gray), fillstyle = :/, linewidth = 3, - legend = nothing, ticks = nothing, border = :none) - - eta(x) = x / (x^2 + 1) * sin(2 * pi * x) + H + 1.5 - plot!(x, eta, color = :blue, fill = (b.(x), 0.4, :blue), linewidth = 3) - - x1 = 0.2 - plot!([x1, x1], [b(x1), eta(x1)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), - annotation = (x1 - 0.08, (eta(x1) + b(x1)) / 2, text(L"h(t, x)", fontsize)), - linewidth = 2) - x2 = 0.4 - plot!([x2, x2], [0.0, b(x2)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), - annotation = (x2 + 0.06, b(x2) / 2, text(L"b(x)", fontsize)), linewidth = 2) - x3 = 0.8 - plot!([x3, x3], [0.0, eta(x3)], line = (Plots.Arrow(:open, :both, 2.5, 2.0), :black), - annotation = (x3 - 0.08, eta(x3) / 2, text(L"\eta(t, x)", fontsize)), - linewidth = 2) - - savefig(joinpath(OUT, "bathymetry.pdf")) -end - -# Plot of dispersion relations -function fig_2() - linewidth = 2 - markersize = 5 - - h0 = 1.0 - g = 1.0 - c0 = sqrt(g * h0) - - k = 0.01:0.5:(8 * pi) - k_zoom = 0.01:0.3:pi - ylim = (0.0, 1.1) - - omega_euler(k) = sqrt(g * k) * sqrt(tanh(h0 * k)) - c_euler(k) = omega_euler(k) / k - plot(k, c_euler.(k) ./ c0, label = "Euler", ylim = ylim, xguide = L"k", - yguide = L"c/c_0", linewidth = linewidth, markershape = :circle, - markersize = markersize) - plot!(k_zoom, c_euler.(k_zoom) ./ c0, ylim = (0.54, 1.0), - inset = bbox(0.35, 0.1, 0.35, 0.3), subplot = 2, legend = nothing, - linewidth = linewidth, markershape = :circle, markersize = markersize, - framestyle = :box) - - function plot_dispersion_relation(omega, label, markershape) - c(k) = omega(k) / k - plot!(k, c.(k) ./ c0, label = label, linewidth = linewidth, - markershape = markershape, markersize = markersize) - plot!(k_zoom, c.(k_zoom) ./ c0, subplot = 2, linewidth = linewidth, - markershape = markershape, markersize = markersize) - end - - omega_bbmbbm_(k, d0) = sqrt(g * h0) * k / (1 + 1 / 6 * (d0 * k)^2) - omega_bbmbbm(k) = omega_bbmbbm_(k, h0) - plot_dispersion_relation(omega_bbmbbm, "BBM-BBM", :cross) - - alpha_set1 = -1 / 3 * c0 * h0^2 - beta_set1 = 0.0 * h0^3 - gamma_set1 = 0.0 * c0 * h0^3 - - alpha_set2 = 0.0004040404040404049 * c0 * h0^2 - beta_set2 = 0.49292929292929294 * h0^3 - gamma_set2 = 0.15707070707070708 * c0 * h0^3 - - alpha_set3 = 0.0 * c0 * h0^2 - beta_set3 = 0.27946992481203003 * h0^3 - gamma_set3 = 0.0521077694235589 * c0 * h0^3 - - alpha_set4 = 0.0 * c0 * h0^2 - beta_set4 = 0.2308939393939394 * h0^3 - gamma_set4 = 0.04034343434343434 * c0 * h0^3 - - function char_equation(alpha, beta, gamma, k) - a = (1 + beta / h0 * k^2) - b = (-alpha - beta * alpha / h0 * k^2 - gamma / h0) * k^3 - c = -g * h0 * k^2 + gamma * alpha / h0 * k^6 - omega1 = (-b + sqrt(b^2 - 4 * a * c)) / (2 * a) - # omega2 = (-b - sqrt(b^2 - 4*a*c))/(2*a) - return omega1 - end - - omega_set1(k) = char_equation(alpha_set1, beta_set1, gamma_set1, k) - plot_dispersion_relation(omega_set1, "S.-K. set 1", :rtriangle) - - omega_set2(k) = char_equation(alpha_set2, beta_set2, gamma_set2, k) - plot_dispersion_relation(omega_set2, "S.-K. set 2", :star5) - - omega_set3(k) = char_equation(alpha_set3, beta_set3, gamma_set3, k) - plot_dispersion_relation(omega_set3, "S.-K. set 3", :star8) - - omega_set4(k) = char_equation(alpha_set4, beta_set4, gamma_set4, k) - plot_dispersion_relation(omega_set4, "S.-K. set 4", :diamond) - - # Plot box - plot!([0.0, pi], [0.54, 0.54], color = :black, label = :none) - plot!([0.0, pi], [1.0, 1.0], color = :black, label = :none) - plot!([0.0, 0.0], [0.54, 1.0], color = :black, label = :none) - plot!([pi, pi], [0.54, 1.0], color = :black, label = :none) - - # Plot connecting lines - plot!([pi, 6.8], [0.54, 0.629], color = :black, label = :none) - plot!([pi, 6.8], [1, 1.01], color = :black, label = :none) - - savefig(joinpath(OUT, "dispersion_relations.pdf")) -end - -# Chapter 4 -# Chapter 4.1 Soliton - -const OUT_SOLITON = joinpath(OUT, "soliton") -ispath(OUT_SOLITON) || mkpath(OUT_SOLITON) - -# Plot convergence orders for baseline and relaxation -function fig_3() - tspan = (0.0, 10.0) - accuracy_orders = [2, 4, 6, 8] - iters = [4, 4, 4, 3] - initial_Ns = [128, 128, 128, 128] - - all_Ns = minimum(initial_Ns) * 2 .^ (0:(maximum(iters) - 1)) - - linewidth = 2 - markersize = 5 - markershapes = [:circle, :star5, :star8, :rtriangle] - plot(label = :none, xscale = :log2, yscale = :log10, xlabel = "N", ylim = (1e-5, 1e2), - ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2", - legend = :bottomleft, layout = 2) - - # left subplot: baseline - for i in 1:length(accuracy_orders) - Ns = initial_Ns[i] * 2 .^ (0:(iters[i] - 1)) - _, errormatrix = convergence_test(joinpath(EXAMPLES_DIR_BBMBBM, - "bbm_bbm_1d_basic.jl"), - iters[i]; N = initial_Ns[i], tspan = tspan, - accuracy_order = accuracy_orders[i]) - # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 - l2_err = sum(errormatrix[:l2], dims = 2) - eocs = log.(l2_err[2:end] ./ l2_err[1:(end - 1)]) ./ log(0.5) - eoc_mean = round(sum(eocs) / length(eocs), digits = 2) - plot!(Ns, l2_err, label = "p = $(accuracy_orders[i]), EOC: $eoc_mean", - markershape = markershapes[i], linewidth = linewidth, markersize = markersize, - subplot = 1) - end - - # right subplot: relaxation - for i in 1:length(accuracy_orders) - Ns = initial_Ns[i] * 2 .^ (0:(iters[i] - 1)) - _, errormatrix = convergence_test(joinpath(EXAMPLES_DIR_BBMBBM, - "bbm_bbm_1d_relaxation.jl"), - iters[i]; N = initial_Ns[i], tspan = tspan, - accuracy_order = accuracy_orders[i]) - # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 - l2_err = sum(errormatrix[:l2], dims = 2) - eocs = log.(l2_err[2:end] ./ l2_err[1:(end - 1)]) ./ log(0.5) - eoc_mean = round(sum(eocs) / length(eocs), digits = 2) - plot!(Ns, l2_err, label = "p = $(accuracy_orders[i]), EOC: $eoc_mean", - markershape = markershapes[i], linewidth = linewidth, markersize = markersize, - subplot = 2) - end - xticks!(all_Ns, string.(all_Ns)) - savefig(joinpath(OUT_SOLITON, "orders.pdf")) -end - -# Plot errors, change of invariants, and solution at final time for baseline and relaxation -function fig_4_5_6() - linewidth = 2 - linestyles = [:dash, :dot] - - g = 9.81 - D = 2.0 - c = 5 / 2 * sqrt(g * D) - xmin = -35.0 - xmax = 35.0 - tspan = (0.0, 50 * (xmax - xmin) / c) - N = 512 - accuracy_order = 8 - - # baseline - trixi_include(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_basic.jl"), - gravity_constant = g, D = D, coordinates_min = xmin, - coordinates_max = xmax, tspan = tspan, N = N, - accuracy_order = accuracy_order) - p1 = plot(analysis_callback, title = "", label_extension = "baseline", - linestyles = [:solid :dash :dot], - linewidth = linewidth, layout = 2, subplot = 1) - p2 = plot(analysis_callback, title = "", what = (:errors,), - label_extension = "baseline", linestyle = linestyles[1], - linewidth = linewidth, - ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2", - exclude = [:conservation_error]) - p3 = plot(semi => sol, label = "baseline", plot_initial = true, plot_bathymetry = false, - linestyle = linestyles[1], linewidth = linewidth, plot_title = "", title = "", - ylims = [(-8, 3) (-1, 40)]) - x = DispersiveShallowWater.grid(semi) - q = DispersiveShallowWater.wrap_array(sol.u[end], semi) - plot!(p3, x, view(q, 1, :), inset = (1, bbox(0.11, 0.6, 0.35, 0.32)), subplot = 3, - xlim = (-20, -10), - ylim = (-0.05, 0.05), legend = nothing, linewidth = linewidth, - linestyle = linestyles[1], - color = 2, - tickfontsize = 5, yticks = [0.04, 0.0, -0.04], xticks = [-20, -15, -10], - plot_initial = true, plot_bathymetry = false, framestyle = :box) - q_exact = DispersiveShallowWater.wrap_array(DispersiveShallowWater.compute_coefficients(initial_condition, - tspan[2], - semi), - semi) - plot!(p3, x, view(q_exact, 1, :), subplot = 3, legend = nothing, linewidth = linewidth, - linestyle = :solid, color = 1) - - # relaxation - trixi_include(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"), - gravity_constant = g, D = D, coordinates_min = xmin, - coordinates_max = xmax, tspan = tspan, N = N, - accuracy_order = accuracy_order) - plot!(p1, analysis_callback, title = "", label_extension = "relaxation", - linestyles = [:solid :dash :dot], - linewidth = linewidth, subplot = 2) - plot!(p2, analysis_callback, title = "", what = (:errors,), - label_extension = "relaxation", linestyle = linestyles[2], linewidth = linewidth, - ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2", - exclude = [:conservation_error]) - plot!(p3, semi => sol, plot_bathymetry = false, label = "relaxation", - linestyle = linestyles[2], - linewidth = linewidth, plot_title = "", title = "", color = 3) - x = DispersiveShallowWater.grid(semi) - q = DispersiveShallowWater.wrap_array(sol.u[end], semi) - plot!(p3, x, view(q, 1, :), subplot = 3, legend = nothing, linewidth = linewidth, - linestyle = linestyles[2], color = 3) - - # Plot box - plot!(p3, [-20, -10], [-0.1, -0.1], color = :black, label = :none) - plot!(p3, [-20, -10], [0.1, 0.1], color = :black, label = :none) - plot!(p3, [-20, -20], [-0.1, 0.1], color = :black, label = :none) - plot!(p3, [-10, -10], [-0.1, 0.1], color = :black, label = :none) - - # Plot connecting lines - plot!(p3, [-20, -29], [-0.1, -3.6], color = :black, label = :none) - plot!(p3, [-10, -3.15], [-0.1, -3.6], color = :black, label = :none) - - savefig(p1, joinpath(OUT_SOLITON, "invariants.pdf")) - savefig(p2, joinpath(OUT_SOLITON, "errors.pdf")) - savefig(p3, joinpath(OUT_SOLITON, "solution.pdf")) -end - -# Plot errors for narrow-stencil, wide-stencil and upwind operators (all using relaxation) -function fig_7() - linewidth = 2 - linestyles = [:solid, :dash, :dot, :dashdot] - - g = 9.81 - D = 2.0 - c = 5 / 2 * sqrt(g * D) - xmin = -35.0 - xmax = 35.0 - tspan = (0.0, 15 * (xmax - xmin) / c) - N = 512 - accuracy_order = 8 - - plot() - - D1 = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N) - D2 = sparse(D1)^2 - solver_widestencil = Solver(D1, D2) - - D1 = periodic_derivative_operator(1, accuracy_order, xmin, xmax, N) - D2 = periodic_derivative_operator(2, accuracy_order, xmin, xmax, N) - solver_narrowstencil = Solver(D1, D2) - - D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, - accuracy_order = accuracy_order, xmin = xmin, xmax = xmax, - N = N) - D2 = sparse(D1.plus) * sparse(D1.minus) - solver_upwind = Solver(D1, D2) - solvers = [ - solver_narrowstencil, - solver_narrowstencil, - solver_widestencil, - solver_upwind, - ] - labels = [ - "narrow-stencil", - "narrow-stencil in velocity equation", - "wide-stencil", - "upwind", - ] - examples = [joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"), - joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_relaxation.jl"), - joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_relaxation.jl"), - joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_relaxation.jl")] - - for (i, solver) in enumerate(solvers) - trixi_include(examples[i], - gravity_constant = g, D = D, coordinates_min = xmin, - coordinates_max = xmax, tspan = tspan, N = N, - accuracy_order = accuracy_order, solver = solver) - plot!(analysis_callback, title = "", what = (:errors,), - label_extension = labels[i], linestyle = linestyles[i], - linewidth = linewidth, - ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2", - exclude = [:conservation_error, :linf_error]) - end - - savefig(joinpath(OUT_SOLITON, "errors_stencils.pdf")) -end - -# Compare the orders of narrow-stencil, wide-stencil and upwind SBP operators -# Not used in the thesis, but nonetheless interesting -function fig_orders_different_stencils() - tspan = (0.0, 10.0) - xmin = -35.0 - xmax = 35.0 - accuracy_orders = [2, 4, 6, 8] - iters = [4, 4, 4, 3] - initial_Ns = [128, 128, 128, 128] - - all_Ns = minimum(initial_Ns) * 2 .^ (0:(maximum(iters) - 1)) - - linewidth = 2 - markersize = 5 - markershapes = [:circle, :star5, :star8, :rtriangle] - plot(label = :none, xscale = :log2, yscale = :log10, xlabel = "N", ylim = (1e-4, 1e3), - ylabel = L"\Vert\eta - \eta_{ana}\Vert_2 + \Vert v - v_{ana}\Vert_2", - legend = :topright, layout = (1, 3)) - - # put examples in separate files since the different solvers cannot be set with the convergence_test - # because they depend on N - examples = [ - joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, "bbm_bbm_variable_bathymetry_1d_basic.jl"), - joinpath(VISUALIZATION_DIR, "bbm_bbm_variable_bathymetry_1d_widestencil.jl"), - joinpath(VISUALIZATION_DIR, "bbm_bbm_variable_bathymetry_1d_upwind.jl")] - for (j, example) in enumerate(examples) - for i in 1:length(accuracy_orders) - Ns = initial_Ns[i] * 2 .^ (0:(iters[i] - 1)) - _, errormatrix = convergence_test(example, - iters[i]; N = initial_Ns[i], tspan = tspan, - accuracy_order = accuracy_orders[i], - coordinates_min = xmin, - coordinates_max = xmax) - # Use sum over all L^2-errors for all variables, i.e. ||η - η_ana||_2 + ||v - v_ana||_2 - l2_err = sum(errormatrix[:l2], dims = 2) - eocs = log.(l2_err[2:end] ./ l2_err[1:(end - 1)]) ./ log(0.5) - eoc_mean = round(sum(eocs) / length(eocs), digits = 2) - plot!(Ns, l2_err, label = "p = $(accuracy_orders[i]), EOC: $eoc_mean", - markershape = markershapes[i], linewidth = linewidth, - markersize = markersize, - subplot = j) - end - end - - xticks!(all_Ns, string.(all_Ns)) - plot!(top_margin = 3 * Plots.mm, subplot = 1) - savefig(joinpath(OUT_SOLITON, "orders_stencils.pdf")) -end - -# Chapter 4.2 Lake-at-rest -const OUT_LAKEATREST = joinpath(OUT, "lake_at_rest") -ispath(OUT_LAKEATREST) || mkpath(OUT_LAKEATREST) - -# Lake-at-rest error for long-time simulation with discontinuous bottom -function fig_8() - linewidth = 2 - N = 100 - accuracy_order = 4 - xmin = -1.0 - xmax = 1.0 - tspan = (0.0, 100.0) - D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, - accuracy_order = accuracy_order, xmin = xmin, xmax = xmax, - N = N) - D2 = sparse(D1.plus) * sparse(D1.minus) - solver = Solver(D1, D2) - trixi_include(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_well_balanced.jl"); - N = N, tspan = tspan, solver = solver, dt = 0.5) - plot(analysis_callback, exclude = [:waterheight_total, :velocity, :entropy], - label_extension = "BBM-BBM", plot_title = "", title = "", - ylabel = "lake-at-rest error", linewidth = linewidth) - - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_well_balanced.jl"); - N = N, tspan = tspan, solver = solver, dt = 0.003) - plot!(analysis_callback, exclude = [:waterheight_total, :momentum, :entropy], - label_extension = "Svärd-Kalisch", plot_title = "", title = "", - ylabel = "lake-at-rest error", linestyle = :dash, linewidth = linewidth) - savefig(joinpath(OUT_LAKEATREST, "lake_at_rest_error_discontinuous.pdf")) -end - -# Plot of condition number of matrix that needs to be inverted for the Svärd-Kalisch equations for different order of accuracy -# Not used in the thesis, but nonetheless interesting -function fig_condition_number() - xmin = -1.0 - xmax = 1.0 - accuracy_orders = [2, 4, 6, 8] - eta0 = 2.0 - beta = 0.49292929292929294 - Ns = 10:10:500 - conds = Array{Float64}(undef, length(Ns)) - plot(xlabel = "N", ylabel = L"cond_2") - for accuracy_order in accuracy_orders - for (i, N) in enumerate(Ns) - D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1, - accuracy_order = accuracy_order, xmin = xmin, xmax = xmax, - N = N) - eta = fill(eta0, N) - D = fill(-1.0, N) - for (i, x) in enumerate(SummationByPartsOperators.grid(D1)) - if x >= 0.5 && x <= 0.75 - D[i] = -1.5 - 0.5 * sinpi(2.0 * x) - end - end - d = eta0 .+ D - beta_hat = beta * d .^ 3 - - D1betaD1 = sparse(D1.plus) * Diagonal(beta_hat) * sparse(D1.minus) - hmD1betaD1 = Diagonal(eta .+ D) - D1betaD1 - conds[i] = cond(Matrix(hmD1betaD1)) - end - plot!(Ns, conds, label = "p = $accuracy_order", linewidth = 2, linestyle = :auto) - end - savefig(joinpath(OUT_LAKEATREST, "condition_numbers.pdf")) -end - -# Chapter 4.3 Dingemans -const OUT_DINGEMANS = joinpath(OUT, "dingemans") -ispath(OUT_DINGEMANS) || mkpath(OUT_DINGEMANS) - -# Plot of total waterheight for different models at different points in time -function fig_9() - linewidth = 3 - fontsize = 16 - linestyles = [:solid, :dash, :dot] - - N = 512 - steps = [100, 200, 300, 500] - xlims_zoom = [(-25, 0), (5, 30), (20, 45), (-100, -75)] - ylim_zoom = (0.75, 0.85) - - trixi_include(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_dingemans.jl"); - N = N) - plot(layout = (2, 2), ylim = (-0.05, 0.86), size = (1200, 800), - titlefontsize = fontsize) - for (i, step) in enumerate(steps) - plot!(semi => sol, step = step, conversion = waterheight_total, label = "BBM-BBM", - subplot = i, plot_title = "", linewidth = linewidth, legend = :none, - guidefontsize = fontsize, tickfontsize = fontsize, linestyle = linestyles[1]) - plot!(semi => sol, step = step, inset = (i, bbox(0.1, 0.2, 0.6, 0.5)), - conversion = waterheight_total, linewidth = linewidth, legend = :none, - framestyle = :box, xlim = xlims_zoom[i], ylim = ylim_zoom, - subplot = length(steps) + i, plot_title = "", title = "", xguide = "", - yguide = "", linestyle = linestyles[1]) - end - - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, "svaerd_kalisch_1d_dingemans.jl"); - N = N) - for (i, step) in enumerate(steps) - plot!(semi => sol, step = step, conversion = waterheight_total, - label = "Svärd-Kalisch (set 3)", subplot = i, plot_bathymetry = false, - plot_title = "", linewidth = linewidth, legend = :none, - guidefontsize = fontsize, tickfontsize = fontsize, color = 2, - linestyle = linestyles[2]) - plot!(semi => sol, step = step, conversion = waterheight_total, - linewidth = linewidth, legend = :none, framestyle = :box, - xlim = xlims_zoom[i], ylim = ylim_zoom, subplot = length(steps) + i, - plot_title = "", title = "", xguide = "", yguide = "", color = 2, - linestyle = linestyles[2]) - end - - trixi_include("elixir_shallowwater_1d_dingemans.jl") - for (i, step) in enumerate(steps) - pd = PlotData1D(sol.u[step], semi) - plot!(pd["H"], label = "Shallow water", subplot = i, - title = "t = $(round(sol.t[step], digits = 2))", plot_title = "", - linewidth = linewidth, legend = :none, guidefontsize = fontsize, - tickfontsize = fontsize, color = 3, linestyle = linestyles[3]) - plot!(pd["H"], linewidth = linewidth, legend = :none, - framestyle = :box, xlim = xlims_zoom[i], ylim = ylim_zoom, - subplot = length(steps) + i, plot_title = "", title = "", xguide = "", - yguide = "", color = 3, linestyle = linestyles[3]) - end - - # dirty hack to have one legend for all subplots - plot!(subplot = 3, legend_column = 2, bottom_margin = 22 * Plots.mm, - legend = (0.7, -0.34), legendfontsize = 12) - plot!(left_margin = 5 * Plots.mm) - - # plot boxes - for i in 1:length(steps) - plot!([xlims_zoom[i][1], xlims_zoom[i][2]], [ylim_zoom[1], ylim_zoom[1]], - color = :black, label = :none, subplot = i, linewidth = 2) - plot!([xlims_zoom[i][1], xlims_zoom[i][2]], [ylim_zoom[2], ylim_zoom[2]], - color = :black, label = :none, subplot = i, linewidth = 2) - plot!([xlims_zoom[i][1], xlims_zoom[i][1]], [ylim_zoom[1], ylim_zoom[2]], - color = :black, label = :none, subplot = i, linewidth = 2) - plot!([xlims_zoom[i][2], xlims_zoom[i][2]], [ylim_zoom[1], ylim_zoom[2]], - color = :black, label = :none, subplot = i, linewidth = 2) - end - # plot connecting lines - upper_corners = [[-119.5, 0.68], [-9.5, 0.68]] - for i in 1:length(steps) - plot!([xlims_zoom[i][1], upper_corners[1][1]], [ylim_zoom[1], upper_corners[1][2]], - color = :black, label = :none, subplot = i, linewidth = 2) - plot!([xlims_zoom[i][2], upper_corners[2][1]], [ylim_zoom[1], upper_corners[2][2]], - color = :black, label = :none, subplot = i, linewidth = 2) - end - savefig(joinpath(OUT_DINGEMANS, "waterheight_over_time.pdf")) -end - -# Plot of total waterheight for Svärd-Kalisch equations at different points in space and different orders of accuracy -function fig_10() - ylim = (0.75, 0.85) - yticks = [0.76, 0.78, 0.8, 0.82, 0.84] - x_values = [3.04, 9.44, 20.04, 26.04, 30.44, 37.04] - tlims = [ - (15.0, 45.0), - (19.0, 48.0), - (25.0, 52.0), - (30.0, 60.0), - (33.0, 61.0), - (35.0, 65.0), - ] - plot(layout = (3, 2)) - - N = 512 - tspan = (0.0, 70.0) - saveat = range(tspan..., length = 1000) - accuracy_orders = [2, 4, 6] - linestyles = [:solid, :dash, :dot] - - for (i, accuracy_order) in enumerate(accuracy_orders) - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans.jl"); - N = N, tspan = tspan, accuracy_order = accuracy_order, - saveat = saveat) - for (j, x) in enumerate(x_values) - index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x)) - title = "x = $(round(DispersiveShallowWater.grid(semi)[index], digits = 4))" - plot!(semi => sol, x, conversion = waterheight_total, subplot = j, - xlim = tlims[j], ylim = ylim, plot_title = "", title = title, - legend = nothing, yticks = yticks, linewidth = 2, titlefontsize = 10, - label = "p = $accuracy_order ", linestyle = linestyles[i]) - end - end - - plot!(subplot = 5, legend = (0.82, -1.0), legend_column = 3, legendfontsize = 8, - bottom_margin = 10 * Plots.mm) - savefig(joinpath(OUT_DINGEMANS, "waterheight_at_x_accuracy_order.pdf")) -end - -# Plots of total waterheight for Svärd-Kalisch equations at different points in space and different types of solvers -function fig_11() - ylim = (0.75, 0.85) - yticks = [0.76, 0.78, 0.8, 0.82, 0.84] - x_values = [3.04, 9.44, 20.04, 26.04, 30.44, 37.04] - tlims = [ - (15.0, 45.0), - (19.0, 48.0), - (25.0, 52.0), - (30.0, 60.0), - (33.0, 61.0), - (35.0, 65.0), - ] - plot(layout = (3, 2)) - - N = 512 - tspan = (0.0, 70.0) - saveat = range(tspan..., length = 1000) - accuracy_order = 4 - linestyles = [:solid, :dash, :dot] - - coordinates_min = -138.0 - coordinates_max = 46.0 - p = 3 # N needs to be divisible by p + 1 - D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) - uniform_mesh = UniformPeriodicMesh1D(coordinates_min, coordinates_max, div(N, p + 1)) - D1 = couple_discontinuously(D_legendre, uniform_mesh) - D_pl = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus)) - D_min = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus)) - D2 = sparse(D_pl) * sparse(D_min) - solver_DG = Solver(D1, D2) - - p = 4 # N needs to be divisible by p - D_legendre = legendre_derivative_operator(-1.0, 1.0, p + 1) - uniform_mesh = UniformPeriodicMesh1D(coordinates_min, coordinates_max, div(N, p)) - D1 = couple_continuously(D_legendre, uniform_mesh) - D2_legendre = legendre_second_derivative_operator(-1.0, 1.0, p + 1) - D2 = couple_continuously(D2_legendre, uniform_mesh) - solver_CG = Solver(D1, D2) - - solvers = [solver_DG, :none, solver_CG] - labels = ["DG ", "FD ", "CG "] - - for (i, solver) in enumerate(solvers) - if solver == :none - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans.jl"); - N = N, tspan = tspan, accuracy_order = accuracy_order, - saveat = saveat) - else - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans.jl"); - N = N, tspan = tspan, accuracy_order = accuracy_order, - saveat = saveat, solver = solvers[i]) - end - for (j, x) in enumerate(x_values) - index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x)) - title = "x = $(round(DispersiveShallowWater.grid(semi)[index], digits = 4))" - plot!(semi => sol, x, conversion = waterheight_total, subplot = j, - xlim = tlims[j], ylim = ylim, plot_title = "", title = title, - legend = nothing, yticks = yticks, linewidth = 2, titlefontsize = 10, - label = labels[i], linestyle = linestyles[i]) - end - end - - plot!(subplot = 5, legend = (0.86, -1.0), legend_column = 3, legendfontsize = 8, - bottom_margin = 10 * Plots.mm) - savefig(joinpath(OUT_DINGEMANS, "waterheight_at_x_solver_types.pdf")) -end - -# Plot solution at different points in space and invariants for entropy conservative and dissipative schemes -function fig_12_13() - ylim = (0.75, 0.85) - yticks = [0.76, 0.78, 0.8, 0.82, 0.84] - x_values = [3.04, 9.44, 20.04, 26.04, 30.44, 37.04] - tlims = [ - (15.0, 45.0), - (19.0, 48.0), - (25.0, 52.0), - (30.0, 60.0), - (33.0, 61.0), - (35.0, 65.0), - ] - p1 = plot(layout = (3, 2)) - - N = 512 - tspan = (0.0, 70.0) - saveat = range(tspan..., length = 1000) - accuracy_order = 4 - linestyles = [:solid, :dash, :dot] - linewidth = 2 - titlefontsize = 10 - - labels = ["EC baseline", "EC relaxation ", "ED upwind"] - - function plot_at_x(semi, sol, i) - for (j, x) in enumerate(x_values) - index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x)) - title = "x = $(round(DispersiveShallowWater.grid(semi)[index], digits = 4))" - plot!(p1, semi => sol, x, conversion = waterheight_total, subplot = j, - xlim = tlims[j], ylim = ylim, plot_title = "", title = title, - legend = nothing, yticks = yticks, linewidth = linewidth, - titlefontsize = titlefontsize, - label = labels[i], linestyle = linestyles[i]) - end - end - - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, "svaerd_kalisch_1d_dingemans.jl"); - N = N, tspan = tspan, accuracy_order = accuracy_order, saveat = saveat) - plot_at_x(semi, sol, 1) - p2 = plot(analysis_callback, title = labels[1], legend = :none, - linestyles = [:solid :dash :dot], - linewidth = linewidth, layout = 3, subplot = 1, titlefontsize = titlefontsize) - - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans_relaxation.jl"); - N = N, tspan = tspan, accuracy_order = accuracy_order, saveat = saveat) - plot_at_x(semi, sol, 2) - plot!(p2, analysis_callback, title = labels[2], legend = :none, - linestyles = [:solid :dash :dot], - linewidth = linewidth, subplot = 2, titlefontsize = titlefontsize) - - trixi_include(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans_upwind.jl"); - N = N, tspan = tspan, accuracy_order = accuracy_order, saveat = saveat) - plot_at_x(semi, sol, 3) - plot!(p2, analysis_callback, title = labels[3], legend = :none, - linestyles = [:solid :dash :dot], - linewidth = linewidth, subplot = 3, titlefontsize = titlefontsize) - - plot!(p1, subplot = 5, legend = (0.55, -1.1), legend_column = 3, legendfontsize = 8, - bottom_margin = 10 * Plots.mm) - plot!(p2, subplot = 3, legend = (1.3, 0.6), legendfontsize = 8) - savefig(p1, joinpath(OUT_DINGEMANS, "waterheight_at_x_ec.pdf")) - savefig(p2, joinpath(OUT_DINGEMANS, "invariants_ec.pdf")) -end - -fig_1() -fig_2() -fig_3() -fig_4_5_6() -fig_7() -# fig_orders_different_stencils() -fig_8() -# fig_condition_number() -fig_9() -fig_10() -fig_11() -fig_12_13() diff --git a/visualization/elixir_shallowwater_1d_dingemans.jl b/visualization/elixir_shallowwater_1d_dingemans.jl deleted file mode 100644 index 87e83522..00000000 --- a/visualization/elixir_shallowwater_1d_dingemans.jl +++ /dev/null @@ -1,59 +0,0 @@ -using Trixi -using OrdinaryDiffEq - -equations = ShallowWaterEquations1D(gravity_constant = 9.81, H0 = 0.8) - -function initial_condition_dingemans_trixi(x, t, equations::ShallowWaterEquations1D) - eta0 = 0.8 - A = 0.02 - # omega = 2*pi/(2.02*sqrt(2)) - k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * eta0), 1.0) using Roots.jl - if x[1] < -30.5 * pi / k || x[1] > -8.5 * pi / k - h = 0.0 - else - h = A * cos(k * x[1]) - end - v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0 - if x[1] < 11.01 || x[1] >= 33.07 - b = 0.0 - elseif 11.01 <= x[1] && x[1] < 23.04 - b = 0.6 * (x[1] - 11.01) / (23.04 - 11.01) - elseif 23.04 <= x[1] && x[1] < 27.04 - b = 0.6 - elseif 27.04 <= x[1] && x[1] < 33.07 - b = 0.6 * (33.07 - x[1]) / (33.07 - 27.04) - else - error("should not happen") - end - eta = h + eta0 - D = -b - return Trixi.prim2cons(SVector(eta, v, b), equations) -end - -initial_condition = initial_condition_dingemans_trixi - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) -solver = DGSEM(polydeg = 3, surface_flux = surface_flux, - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) - -coordinates_min = -138.0 -coordinates_max = 46.0 -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 7, - n_cells_max = 10_000) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) -tspan = (0.0, 70.0) -ode = Trixi.semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() - -analysis_callback = Trixi.AnalysisCallback(semi, interval = 100) - -callbacks = CallbackSet(summary_callback, analysis_callback) - -saveat = range(tspan..., length = 500) -sol = solve(ode, Tsit5(), reltol = 1e-7, abstol = 1e-7, - save_everystep = false, callback = callbacks, saveat = saveat); -summary_callback() # print the timer summary diff --git a/visualization/plot_examples.jl b/visualization/plot_examples.jl deleted file mode 100644 index 58e2adb8..00000000 --- a/visualization/plot_examples.jl +++ /dev/null @@ -1,211 +0,0 @@ -# This script solves all the examples in `examples_dir()` and creates a .gif of the solution -# as well as some additional plots. All plots are saved in a directory `out/` and -# subdirectories therein. - -using DispersiveShallowWater -using Plots - -# Use a macro to avoid world age issues when defining new initial conditions etc. -# inside an example. -macro plot_example(filename, args...) - local ylims_gif = get_kwarg(args, :ylims_gif, nothing) - local ylims_x = get_kwarg(args, :ylims_x, nothing) - local x_values = get_kwarg(args, :x_values, []) - local tlims = get_kwarg(args, :tlims, []) - - local kwargs = Pair{Symbol, Any}[] - for arg in args - if (arg.head == :(=) && !(arg.args[1] in (:ylims_gif, :ylims_x, :x_values, :tlims))) - push!(kwargs, Pair(arg.args...)) - end - end - - quote - trixi_include(joinpath(examples_dir(), $filename); $kwargs...) - elixirname = splitext(basename($filename))[1] - outdir = joinpath("out", dirname($filename), elixirname) - ispath(outdir) || mkpath(outdir) - - # Plot solution - anim = @animate for step in 1:length(sol.u) - plot(semi => sol, plot_initial = true, step = step, ylims = $ylims_gif) - end - gif(anim, joinpath(outdir, "solution.gif"), fps = 25) - - # Plot error in invariants - plot(analysis_callback) - savefig(joinpath(outdir, "invariants.pdf")) - - # Plot at different x values over time - @assert size($x_values) == size($tlims) - for (i, x) in enumerate($x_values) - plot(semi => sol, x, xlim = $tlims[i], ylims = $ylims_x) - savefig(joinpath(outdir, "solution_at_x_" * string(x) * ".pdf")) - end - end -end - -# Get the first value assigned to `keyword` in `args` and return `default_value` -# if there are no assignments to `keyword` in `args`. -function get_kwarg(args, keyword, default_value) - val = default_value - for arg in args - if arg.head == :(=) && arg.args[1] == keyword - val = arg.args[2] - break - end - end - return val -end - -const EXAMPLES_DIR_BBMBBM = "bbm_bbm_1d" -const EXAMPLES_DIR_BBMBBM_VARIABLE = "bbm_bbm_variable_bathymetry_1d" -const EXAMPLES_DIR_SVAERD_KALISCH = "svaerd_kalisch_1d" - -############################################################################### -# Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions -# using periodic SBP operators -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_basic.jl"), - ylims_gif=[(-8, 4) :auto], tspan=(0.0, 50.0)) - -############################################################################### -# Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions -# using discontinuously coupled Legendre SBP operators -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_dg.jl"), - ylims_gif=[(-4, 2) :auto]) - -############################################################################### -# Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions -# using periodic SBP operators and relaxation, is energy-conservative -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM, "bbm_bbm_1d_relaxation.jl"), - ylims_gif=[(-8, 4) (-10, 30)], - tspan=(0.0, 30.0)) - -############################################################################### -# Travelling wave solution for one-dimensional BBM-BBM equations with periodic boundary conditions -# using periodic SBP operators. Uses the BBM-BBM equations with variable bathymetry, but sets the bathymetry -# as a constant. Should give the same result as "bbm_bbm_1d_basic.jl" -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_basic.jl"), - ylims_gif=[(-8, 4) :auto], - tspan=(0.0, 50.0)) - -############################################################################### -# One-dimensional BBM-BBM equations with a Gaussian bump as initial condition for the water height -# and initially still water. The bathymetry is a sine function. Relaxation is used, so the solution -# is energy-conservative. Uses periodic finite difference SBP operators. -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_relaxation.jl"), - ylims_gif=[(-1.5, 6.0) (-10.0, 10.0)], - tspan=(0.0, 10.0)) - -############################################################################### -# One-dimensional BBM-BBM equations with a Gaussian bump as initial condition for the water height -# and initially still water. The bathymetry is a sine function. Relaxation is used, so the solution -# is energy-conservative. Uses upwind discontinuously coupled SBP operators. -@plot_elixir(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_dg_upwind_relaxation.jl"), - ylims_gif=[(-1.5, 6.0) (-10.0, 10.0)], - tspan=(0.0, 10.0)) - -############################################################################### -# One-dimensional BBM-BBM equations with a Gaussian bump as initial condition for the water height -# and initially still water. The bathymetry is a sine function. Relaxation is used, so the solution -# is energy-conservative. Uses periodic finite difference discontinuously coupled SBP operators. -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_upwind_relaxation.jl"), - ylims_gif=[(-1.5, 6.0) (-10.0, 10.0)], - tspan=(0.0, 10.0)) - -############################################################################### -# One-dimensional BBM-BBM equations with a constant water height -# and initially still water. The bathymetry is discontinuous. Relaxation is used, so the solution -# is energy-conservative. Uses periodic finite difference SBP operators. The solution should be -# (exactly) constant in time. -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_well_balanced.jl"), - ylims_gif=[(2.0 - 1e-3, 2.0 + 1e-3) (-1e-3, 1e-3)], - tspan=(0.0, 10.0)) - -############################################################################### -# One-dimensional BBM-BBM equations with initial condition that models -# a wave make. This setup comes from experiments by W. M. Dingemans. -@plot_example(joinpath(EXAMPLES_DIR_BBMBBM_VARIABLE, - "bbm_bbm_variable_bathymetry_1d_dingemans.jl"), - ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)], - ylims_x=[:auto :auto], - x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04], - tlims=[ - (15.0, 45.0), - (19.0, 48.0), - (25.0, 52.0), - (30.0, 60.0), - (33.0, 61.0), - (35.0, 65.0), - ], - tspan=(0.0, 70.0)) - -############################################################################### -# One-dimensional equations from Svärd and Kalisch with initial condition that models -# a wave make. This setup comes from experiments by W. M. Dingemans. -@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans.jl"), - ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)], - ylims_x=[:auto :auto], - x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04], - tlims=[ - (15.0, 45.0), - (19.0, 48.0), - (25.0, 52.0), - (30.0, 60.0), - (33.0, 61.0), - (35.0, 65.0), - ], - tspan=(0.0, 70.0)) - -############################################################################### -# One-dimensional equations from Svärd and Kalisch with initial condition that models -# a wave make. This setup comes from experiments by W. M. Dingemans. -@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans_upwind.jl"), - ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)], - ylims_x=[:auto :auto], - x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04], - tlims=[ - (15.0, 45.0), - (19.0, 48.0), - (25.0, 52.0), - (30.0, 60.0), - (33.0, 61.0), - (35.0, 65.0), - ], - tspan=(0.0, 70.0)) - -############################################################################### -# One-dimensional equations from Svärd and Kalisch with initial condition that models -# a wave make. This setup comes from experiments by W. M. Dingemans. Relaxation is used -# to preserve the modified entropy. -@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_dingemans_relaxation.jl"), - ylims_gif=[(-0.1, 0.9) (-0.3, 0.3)], - ylims_x=[:auto :auto], - x_values=[3.04, 9.44, 20.04, 26.04, 30.44, 37.04], - tlims=[ - (15.0, 45.0), - (19.0, 48.0), - (25.0, 52.0), - (30.0, 60.0), - (33.0, 61.0), - (35.0, 65.0), - ], - tspan=(0.0, 70.0)) - -############################################################################### -# One-dimensional Svärd-Kalisch equations with a constant water height -# and initially still water. The bathymetry is discontinuous. Relaxation is used, so the solution -# is energy-conservative. Uses periodic finite difference SBP operators. The solution should be -# (exactly) constant in time. -@plot_example(joinpath(EXAMPLES_DIR_SVAERD_KALISCH, - "svaerd_kalisch_1d_well_balanced.jl"), - ylims=[(2.0 - 1e-3, 2.0 + 1e-3) (-1e-3, 1e-3)], - tspan=(0.0, 10.0))