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 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:
\[\begin{aligned}
+ \eta_t + ((\eta + D)v)_x - \frac{1}{6}(D^2\eta_{xt})_x &= 0,\\
+ v_t + g\eta_x + \left(\frac{1}{2}v^2\right)_x - \frac{1}{6}(D^2v_t)_{xx} &= 0.
+\end{aligned}\]
Here, $\eta = h + b$ describes the total water height, $h$ the water height above the bottom topography (bathymetry), $b = -D$ the bathymetry and $v$ the velocity in horizontal direction. The gravitational acceleration is denoted as $g$. A sketch of the water height and the bathymetry can be found below.
In order to conduct a numerical simulation with DispersiveShallowWater.jl, we perform the following steps.
First, we load the necessary libraries:
using DispersiveShallowWater, OrdinaryDiffEq
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 travelling wave that moves towards a beach, which is modeled by a linearly increasing bathymetry.
equations = BBMBBMVariableEquations1D(gravity_constant = 9.81)
+
+function initial_condition_shoaling(x, t, equations::BBMBBMVariableEquations1D, mesh)
+ A = 0.07 # amplitude of wave
+ x0 = -30 # initial center
+ eta = A * exp(-0.1*(x - x0)^2)
+ v = 0
+ D = x <= 0.0 ? 0.7 : 0.7 - 1/50 * x
+ return SVector(eta, v, D)
+end
+
+initial_condition = initial_condition_shoaling
+boundary_conditions = boundary_condition_periodic
+
+coordinates_min = -130.0
+coordinates_max = 20.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N + 1)
Mesh1D{Float64}
+ xmin: -130.0
+ xmax: 20.0
+ N: 513
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 since, up to now, DispersiveShallowWater.jl only provides support for this type of boundary conditions. Lastly, we define the physical domain as the interval from -130 to 20 and we choose 512 intermediate nodes. The mesh is homogeneous, i.e. the distance between each two nodes is constant. We choose the left boundary very far to the left in order to avoid an interaction of the left- and right-travelling waves.
In the next step, we build a Semidiscretization that bundles all ingredients for the spatial discretization of the model. Especially, we need to define a Solver. 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.
solver = Solver(mesh, 4)
+
+semi = Semidiscretization(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions)
Semidiscretization
+ #spatial dimensions: 1
+ mesh: Mesh1D{Float64} with length 513
+ equations: BBMBBMVariableEquations1D
+ initial condition: initial_condition_shoaling
+ boundary condition: boundary_condition_periodic
Finally, we put the mesh, the equations, the initial_condition, the solver and the boundary_conditions together in a semidiscretization semi.
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 by calling semidiscretize on the semidiscretization and the time span. Additionally, we can analyze the numerical solution using an AnalysisCallback. 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 waterheight $\eta$ waterheight_total, the velocity and the entropy are computed and saved for each time step. The total waterheight and the total velocity are linear invariants of the BBM-BBM equations, i.e. they do not change over time. The total entropy
\[\mathcal E(t; \eta, v) = \frac{1}{2}\int_\Omega g\eta^2 + (\eta + D)v^2\textrm{d}x\]
is a nonlinear invariant and should be constant over time as well. During the simulation, the AnalysisCallback will print the results to the terminal.
Finally, the ode can be solved using the interface from 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.
tspan = (0.0, 25.0)
+ode = semidiscretize(semi, tspan)
+analysis_callback = AnalysisCallback(semi; interval = 10,
+ extra_analysis_errors = (:conservation_error,),
+ extra_analysis_integrals = (waterheight_total,
+ velocity, entropy),
+ io = devnull)
+callbacks = CallbackSet(analysis_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+ save_everystep = false, callback = callbacks, saveat = saveat)
retcode: Success
+Interpolation: 1st order linear
+t: 100-element Vector{Float64}:
+ 0.0
+ 0.25252525252525254
+ 0.5050505050505051
+ 0.7575757575757576
+ 1.0101010101010102
+ 1.2626262626262625
+ 1.5151515151515151
+ 1.7676767676767677
+ 2.0202020202020203
+ 2.272727272727273
+ ⋮
+ 22.97979797979798
+ 23.232323232323232
+ 23.484848484848484
+ 23.737373737373737
+ 23.98989898989899
+ 24.242424242424242
+ 24.494949494949495
+ 24.747474747474747
+ 25.0
+u: 100-element Vector{Vector{Float64}}:
+ [0.0, 0.0, 0.7, 0.0, 0.0, 0.7, 0.0, 0.0, 0.7, 0.0 … 0.32339181286549706, 1.1159911034969113e-106, 0.0, 0.31754385964912274, 6.2568988186290324e-108, 0.0, 0.31169590643274847, 3.4485093610035348e-109, 0.0, 0.30584795321637426]
+ [-3.3321047588718723e-56, 8.251130546822754e-57, 0.7, 3.4075338894055293e-56, 3.9326309362743164e-59, 0.7, -1.283173375043596e-56, 4.5802968391757424e-57, 0.7, 9.844834652703703e-57 … 0.32339181286549706, 9.30667650388066e-55, -7.283148374436256e-55, 0.31754385964912274, -3.350796824759337e-55, 2.838498358666995e-55, 0.31169590643274847, 9.268569900755344e-56, -9.95556903369101e-56, 0.30584795321637426]
+ [-4.8547983995290075e-53, 1.8098135295902446e-53, 0.7, 4.842606456334392e-53, -1.4594298990367477e-55, 0.7, -1.8544463592440024e-53, 9.787993172206424e-54, 0.7, 1.4193303481545844e-53 … 0.32339181286549706, 1.2619638760700319e-51, -1.4765571825182888e-51, 0.31754385964912274, -4.641540581198118e-52, 5.871522647401302e-52, 0.31169590643274847, 1.320914893327927e-52, -2.0982370814323856e-52, 0.30584795321637426]
+ [-1.161242193741655e-50, 5.352623097482111e-51, 0.7, 1.1320524430148154e-50, -1.0908865403610494e-52, 0.7, -4.4040845541169946e-51, 2.825787874193396e-51, 0.7, 3.3635723672500614e-51 … 0.32339181286549706, 2.818933981521643e-49, -4.048343871083631e-49, 0.31754385964912274, -1.058197448349471e-49, 1.6416275820272595e-49, 0.31169590643274847, 3.094178030291858e-50, -5.973671280788524e-50, 0.30584795321637426]
+ [-1.1119854014450027e-48, 5.886890378685159e-49, 0.7, 1.0619537996508051e-48, -1.8735941788709885e-50, 0.7, -4.1915194291215604e-49, 3.0416513299400544e-49, 0.7, 3.1953508449033937e-49 … 0.32339181286549706, 2.533649020285857e-47, -4.1512247516289456e-47, 0.31754385964912274, -9.69339932925456e-48, 1.7143529346264834e-47, 0.31169590643274847, 2.906332549421794e-48, -6.344239666763591e-48, 0.30584795321637426]
+ [-6.211498239967898e-47, 3.6607311826978715e-47, 0.7, 5.815475691322198e-47, -1.5891572584843106e-48, 0.7, -2.3283193351518106e-47, 1.8525174302835778e-47, 0.7, 1.7719355506070908e-47 … 0.32339181286549706, 1.3288099374184154e-45, -2.4065744350964748e-45, 0.31754385964912274, -5.181054795720947e-46, 1.0122888587233171e-45, 0.31169590643274847, 1.5926499782855233e-46, -3.8101344190304886e-46, 0.30584795321637426]
+ [-2.3254636832401748e-45, 1.492195217448078e-45, 0.7, 2.137213732722148e-45, -8.159201413073166e-47, 0.7, -8.674202199248216e-46, 7.40703410470262e-46, 0.7, 6.591215847050772e-46 … 0.32339181286549706, 4.681077910562855e-44, -9.168674170174745e-44, 0.31754385964912274, -1.8589374303101976e-44, 3.925994514103777e-44, 0.31169590643274847, 5.853556853116871e-45, -1.5021291360961555e-44, 0.30584795321637426]
+ [-6.392817767444815e-44, 4.4001360552137954e-44, 0.7, 5.774217233847734e-44, -2.8902522521195595e-45, 0.7, -2.3744584776238255e-44, 2.1453566106837627e-44, 0.7, 1.8017470562625006e-44 … 0.32339181286549706, 1.2132383779023626e-42, -2.5326738369849148e-42, 0.31754385964912274, -4.904477852717804e-43, 1.1034157858015082e-42, 0.31169590643274847, 1.58074539515666e-43, -4.289517530474068e-43, 0.30584795321637426]
+ [-1.3787503725949053e-42, 1.0073886682961005e-42, 0.7, 1.225119365614405e-42, -7.708120841618133e-44, 0.7, -5.102248914644414e-43, 4.830053503339048e-43, 0.7, 3.866658258223505e-43 … 0.32339181286549706, 2.470475164410606e-41, -5.441183704745909e-41, 0.31754385964912274, -1.0162156416809096e-41, 2.411037908305353e-41, 0.31169590643274847, 3.3505723972756644e-42, -9.519811676617257e-42, 0.30584795321637426]
+ [-2.4329942874765706e-41, 1.872657301401079e-41, 0.7, 2.1286192694904517e-41, -1.6338489485120413e-42, 0.7, -8.975456699039445e-42, 8.838532478005529e-42, 0.7, 6.79390435768371e-42 … 0.32339181286549706, 4.120411812198934e-40, -9.504003188795472e-40, 0.31754385964912274, -1.7241821494301705e-40, 4.281941434433367e-40, 0.31169590643274847, 5.812886428139709e-41, -1.7167524470673666e-40, 0.30584795321637426]
+ ⋮
+ [-0.002176856805330325, -0.005471508691304251, 0.7, -0.0010991531445136095, -0.0029871877626790164, 0.7, -0.0005923337764093438, -0.0016674303283173346, 0.7, 0.0007288095726285177 … 0.32339181286549706, -0.002446501805835471, -0.01637377384830738, 0.31754385964912274, -0.0016159836137552256, -0.013551798527814577, 0.31169590643274847, -0.002406693562270947, -0.015385993847199556, 0.30584795321637426]
+ [-0.002505875471498055, -0.008500215882787399, 0.7, -0.0018252401083375593, -0.008052197913994339, 0.7, -0.0012170710172801332, -0.0067102232275051904, 0.7, -0.0004346836902188495 … 0.32339181286549706, -0.002185674622026492, -0.015967209051955054, 0.31754385964912274, -0.0019264030106736325, -0.013402096254005758, 0.31169590643274847, -0.002238819035873452, -0.014632665884650377, 0.30584795321637426]
+ [-0.0018622554232660814, -0.008813488862025683, 0.7, -0.0021123681485556265, -0.009894126107164319, 0.7, -0.0020128406188384905, -0.0090495784363708, 0.7, -0.0018802163331252746 … 0.32339181286549706, -0.0014215010258198128, -0.009963901856291011, 0.31754385964912274, -0.0019272324089571029, -0.012910684832531511, 0.31169590643274847, -0.001510589219542568, -0.013900945522401677, 0.30584795321637426]
+ [-0.001133675860455097, -0.0065355391224715086, 0.7, -0.0019792195079004647, -0.0076241143777409925, 0.7, -0.0023947609189684006, -0.008186003096997648, 0.7, -0.002719468672374811 … 0.32339181286549706, -0.00034224853703392327, -0.004573764010269953, 0.31754385964912274, -0.0011555137032242004, -0.011452364518071355, 0.31169590643274847, -0.000895344278176558, -0.012890538549893145, 0.30584795321637426]
+ [-0.000851135833074273, -0.0035764780554477377, 0.7, -0.0015164577577977954, -0.0038311634821191435, 0.7, -0.0019345932262210196, -0.005644616064316189, 0.7, -0.0024032659014425048 … 0.32339181286549706, 0.00048676953199567566, -0.003240993175422954, 0.31754385964912274, -0.0001314980022209947, -0.007826563056479848, 0.31169590643274847, -0.0006596159967522925, -0.009519268041658644, 0.30584795321637426]
+ [-0.0007926983849507133, -0.0017254112742891728, 0.7, -0.0009183634647563134, -0.0015781479578823652, 0.7, -0.0009349159821018389, -0.0034054992530152357, 0.7, -0.0013306484115036508 … 0.32339181286549706, 0.0007085992217668003, -0.003227914524735484, 0.31754385964912274, 0.000367272135332628, -0.0027193442234356067, 0.31169590643274847, -0.0005220550795835639, -0.004148042975382936, 0.30584795321637426]
+ [-0.0005045949122518838, -0.001062607734299834, 0.7, -0.00039778847670194907, -0.0013768855203502293, 0.7, -0.00012800109379717544, -0.0022001075693594907, 0.7, -0.0004075931121329941 … 0.32339181286549706, 0.0005311117049730989, -0.0007705176971119265, 0.31754385964912274, 0.0002934328821294129, 0.001461899074030676, 0.31169590643274847, -0.00014917439661158696, 0.0005201307641518416, 0.30584795321637426]
+ [7.803777012574535e-5, -0.00042936266828309447, 0.7, -4.958365526502639e-5, -0.0013994492190678646, 0.7, 8.636575719676265e-5, -0.0013855062117349463, 0.7, -0.0001392132859487371 … 0.32339181286549706, 0.00044778090104512366, 0.0034101223194062317, 0.31754385964912274, 0.00020695567860438865, 0.0031178119822827246, 0.31169590643274847, 0.0003893407334412704, 0.0027561528810216467, 0.30584795321637426]
+ [0.0006073616394551456, 0.0008907054710588314, 0.7, 0.000180731216408728, -8.677905121445732e-5, 0.7, -2.611552969243042e-5, -0.00015434992446160915, 0.7, -0.0002524639408819109 … 0.32339181286549706, 0.0006577431971269266, 0.005711584342451313, 0.31754385964912274, 0.00043018536651269755, 0.0030669676280136887, 0.31169590643274847, 0.0007358879405425351, 0.0033152700567688075, 0.30584795321637426]
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 NamedTuples by errors(analysis_callback) and integrals(analysis_callback).
After running the simulation, the results can be visualized using 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.
using Plots
+
+plot(semi => sol)
+savefig("shoaling_solution.png")
"/home/runner/work/DispersiveShallowWater.jl/DispersiveShallowWater.jl/docs/build/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 or waterheight_total if one only wants to plot the total waterheight. The resulting plot will have one subplot for each of the returned variables of the conversion variable. By default, the conversion function is just prim2prim, 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.
anim = @animate for step in 1:length(sol.u)
+ plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), ylims = [(-0.8, 0.1)])
+end
+gif(anim, "shoaling_solution.gif", fps = 25)
It is also possible to plot the solution variables at a fixed spatial point over time by calling plot(semi => sol, x) for some x-value, see plot_examples.jl for some examples.
Often, it is interesting to have a look at how the quantities that are recorded by the AnalysisCallback evolve in time. To this end, you can plot the AnalysisCallback by
plot(analysis_callback)
+savefig("analysis_callback.png")
"/home/runner/work/DispersiveShallowWater.jl/DispersiveShallowWater.jl/docs/build/analysis_callback.png"
This creates the following figure:
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.
To obtain entropy-conserving time-stepping schemes DispersiveShallowWater.jl uses the relaxation method introduced in and further developed in . The relaxation method is implemented as a RelaxationCallback, 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:
analysis_callback = AnalysisCallback(semi; interval = 10,
+ extra_analysis_errors = (:conservation_error,),
+ extra_analysis_integrals = (waterheight_total,
+ velocity, entropy),
+ io = devnull)
+relaxation_callback = RelaxationCallback(invariant = entropy)
+callbacks = CallbackSet(relaxation_callback, analysis_callback)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+ save_everystep = false, callback = callbacks, saveat = saveat)
retcode: Success
+Interpolation: 1st order linear
+t: 100-element Vector{Float64}:
+ 0.0
+ 0.25252525252525254
+ 0.5050505050505051
+ 0.7575757575757576
+ 1.0101010101010102
+ 1.2626262626262625
+ 1.5151515151515151
+ 1.7676767676767677
+ 2.0202020202020203
+ 2.272727272727273
+ ⋮
+ 22.97979797979798
+ 23.232323232323232
+ 23.484848484848484
+ 23.737373737373737
+ 23.98989898989899
+ 24.242424242424242
+ 24.494949494949495
+ 24.747474747474747
+ 25.0
+u: 100-element Vector{Vector{Float64}}:
+ [0.0, 0.0, 0.7, 0.0, 0.0, 0.7, 0.0, 0.0, 0.7, 0.0 … 0.32339181286549706, 1.1159911034969113e-106, 0.0, 0.31754385964912274, 6.2568988186290324e-108, 0.0, 0.31169590643274847, 3.4485093610035348e-109, 0.0, 0.30584795321637426]
+ [-3.3297849150159086e-56, 8.241895836085858e-57, 0.7, 3.4052185300811195e-56, 3.8881748563163253e-59, 0.7, -1.282263578596884e-56, 4.575170044505072e-57, 0.7, 9.837918819156205e-57 … 0.32339181286549706, 9.30063321540606e-55, -7.276268643806448e-55, 0.31754385964912274, -3.348542619905238e-55, 2.8354727931795833e-55, 0.31169590643274847, 9.26217699429647e-56, -9.944702313348365e-56, 0.30584795321637426]
+ [-4.854358932980195e-53, 1.8096131337839353e-53, 0.7, 4.842178696468332e-53, -1.4590119227273747e-55, 0.7, -1.8542798548245878e-53, 9.786942140740117e-54, 0.7, 1.419203146786536e-53 … 0.32339181286549706, 1.2618571098416748e-51, -1.476408228991929e-51, 0.31754385964912274, -4.641140459485744e-52, 5.870913979781286e-52, 0.31169590643274847, 1.3207978093786314e-52, -2.098014486656577e-52, 0.30584795321637426]
+ [-1.1611677105185517e-50, 5.352253124846145e-51, 0.7, 1.131980612153987e-50, -1.0907865764699485e-52, 0.7, -4.4038029803197864e-51, 2.8255949510378287e-51, 0.7, 3.363357527736327e-51 … 0.32339181286549706, 2.8187589821588535e-49, -4.048074914318266e-49, 0.31754385964912274, -1.058131124345205e-49, 1.6415173993859553e-49, 0.31169590643274847, 3.0939815656145766e-50, -5.973266555947453e-50, 0.30584795321637426]
+ [-1.1119326231997787e-48, 5.886595840555343e-49, 0.7, 1.0619038621697752e-48, -1.8734843160733145e-50, 0.7, -4.19132101207301e-49, 3.0415006396698293e-49, 0.7, 3.1951997058944475e-49 … 0.32339181286549706, 2.5335322329688205e-47, -4.1510237734450563e-47, 0.31754385964912274, -9.692948588306525e-48, 1.7142692229097626e-47, 0.31169590643274847, 2.906195831886557e-48, -6.34392741624436e-48, 0.30584795321637426]
+ [-6.211340532217373e-47, 3.660634557587833e-47, 0.7, 5.815329227299939e-47, -1.589110730509207e-48, 0.7, -2.3282603489608144e-47, 1.852468935079328e-47, 0.7, 1.7718906905846917e-47 … 0.32339181286549706, 1.3287770760107856e-45, -2.4065127076113577e-45, 0.31754385964912274, -5.180925644645881e-46, 1.0122626968239563e-45, 0.31169590643274847, 1.5926098618954286e-46, -3.8100352622923995e-46, 0.30584795321637426]
+ [-2.3254161163064466e-45, 1.4921634451642905e-45, 0.7, 2.137170436417654e-45, -8.159008862774204e-47, 0.7, -8.674025199829319e-46, 7.406877933416167e-46, 0.7, 6.591081457003216e-46 … 0.32339181286549706, 4.680985227093627e-44, -9.16848575783732e-44, 0.31754385964912274, -1.858900258859365e-44, 3.9259130714242533e-44, 0.31169590643274847, 5.853438288559902e-45, -1.5020977044370005e-44, 0.30584795321637426]
+ [-6.392704749915101e-44, 4.400055673573975e-44, 0.7, 5.774116043216231e-44, -2.890195018350533e-45, 0.7, -2.3744165857848013e-44, 2.1453177616145644e-44, 0.7, 1.801715290597514e-44 … 0.32339181286549706, 1.2132175751434895e-42, -2.5326290720182764e-42, 0.31754385964912274, -4.904392969886902e-43, 1.1033961094737897e-42, 0.31169590643274847, 1.580717704226583e-43, -4.28944041824579e-43, 0.30584795321637426]
+ [-1.37872581855927e-42, 1.0073701184950383e-42, 0.7, 1.2250977603393081e-42, -7.707966800413413e-44, 0.7, -5.102158236014128e-43, 4.829965430378285e-43, 0.7, 3.8665895903999155e-43 … 0.32339181286549706, 2.470432699751572e-41, -5.441087292069507e-41, 0.31754385964912274, -1.0161979839696196e-41, 2.4109947428435356e-41, 0.31169590643274847, 3.3505133570864116e-42, -9.519639622274105e-42, 0.30584795321637426]
+ [-2.4329533363631523e-41, 1.8726248499178522e-41, 0.7, 2.128583769263703e-41, -1.633818571864118e-42, 0.7, -8.975305888394391e-42, 8.838380696071164e-42, 0.7, 6.79379028093039e-42 … 0.32339181286549706, 4.120344807669663e-40, -9.503844505374361e-40, 0.31754385964912274, -1.7241538134168129e-40, 4.281869218700696e-40, 0.31169590643274847, 5.812789588337876e-41, -1.7167232265535456e-40, 0.30584795321637426]
+ ⋮
+ [-0.002165983808856068, -0.005462316236984185, 0.7, -0.0010949955191127826, -0.0029946780624175756, 0.7, -0.0005920726738235036, -0.001663002733240682, 0.7, 0.0007264876839319572 … 0.32339181286549706, -0.002445583550959634, -0.016317015455508743, 0.31754385964912274, -0.0016212354480480459, -0.0135289108928958, 0.31169590643274847, -0.002396885241065637, -0.015332363345196748, 0.30584795321637426]
+ [-0.002494650957585022, -0.008474305024845712, 0.7, -0.0018211108261937223, -0.008025725779320049, 0.7, -0.001219316764107308, -0.006688548166333462, 0.7, -0.00043839507595232945 … 0.32339181286549706, -0.0021804302393484075, -0.015901863182964018, 0.31754385964912274, -0.001926215965289191, -0.013401193628170173, 0.31169590643274847, -0.0022324484169123932, -0.014629354103317616, 0.30584795321637426]
+ [-0.0018597237172829416, -0.00878590816730313, 0.7, -0.002108029863372187, -0.009851920630150503, 0.7, -0.0020098381159790635, -0.009018710861238517, 0.7, -0.0018764289682748508 … 0.32339181286549706, -0.0014174867196700058, -0.009955159836300121, 0.31754385964912274, -0.0019200772426444747, -0.012904053282721307, 0.31169590643274847, -0.0015113103521350384, -0.013907495809161632, 0.30584795321637426]
+ [-0.00113828874159783, -0.0065262409142084485, 0.7, -0.001975563757344359, -0.0076033730942781244, 0.7, -0.0023866475979837714, -0.0081656639530465, 0.7, -0.0027089237075481136 … 0.32339181286549706, -0.00034409204632864453, -0.00461356538213619, 0.31754385964912274, -0.001150846246748618, -0.011431378783412934, 0.31169590643274847, -0.0008995601062787645, -0.012874967436299414, 0.30584795321637426]
+ [-0.0008547943895035903, -0.0035873298420343543, 0.7, -0.0015148921303085032, -0.0038438470097323266, 0.7, -0.0019284466945888319, -0.005644129387757497, 0.7, -0.002394651334213295 … 0.32339181286549706, 0.00048076743678606326, -0.0032655471937228053, 0.31754385964912274, -0.00013500790692118908, -0.007808778243089342, 0.31169590643274847, -0.0006613057774831235, -0.009496236940962872, 0.30584795321637426]
+ [-0.0007911392786587517, -0.00173935695552966, 0.7, -0.0009192247258541208, -0.0016012733935531832, 0.7, -0.0009360248152808115, -0.003415976453312026, 0.7, -0.0013305081833950056 … 0.32339181286549706, 0.0007042413951191663, -0.0032091589879661687, 0.31754385964912274, 0.0003609870405987438, -0.002724084533864073, 0.31169590643274847, -0.0005197291955476894, -0.004146557210079369, 0.30584795321637426]
+ [-0.0005009119017044968, -0.0010643539943987947, 0.7, -0.00039971046012839053, -0.0013829707632917734, 0.7, -0.00013392608631433202, -0.002206735492919327, 0.7, -0.0004136793726674155 … 0.32339181286549706, 0.0005317240886545098, -0.0007458554006872201, 0.31754385964912274, 0.0002918660500146642, 0.001438538950550272, 0.31169590643274847, -0.0001466728723783037, 0.0004996358322929233, 0.30584795321637426]
+ [7.877033748192486e-5, -0.00042112363962542665, 0.7, -5.0759247674105346e-5, -0.0013861889826008388, 0.7, 8.248268091035848e-5, -0.0013829603946877696, 0.7, -0.00014394433252075652 … 0.32339181286549706, 0.0004509502003486682, 0.003401019776759001, 0.31754385964912274, 0.00020988107240593255, 0.0030986111294804813, 0.31169590643274847, 0.00038828103079828485, 0.0027368433078830117, 0.30584795321637426]
+ [0.0006041217258655958, 0.0008950650498073837, 0.7, 0.00018078926007732303, -7.445962324102933e-5, 0.7, -2.451652846529534e-5, -0.0001493914009530353, 0.7, -0.00025139259773707763 … 0.32339181286549706, 0.0006585897415158146, 0.00567992343085633, 0.31754385964912274, 0.0004319778826393798, 0.003066796203847099, 0.31169590643274847, 0.0007320680612231083, 0.003311140575127919, 0.30584795321637426]
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.
plot(analysis_callback)
+savefig("analysis_callback_relaxation.png")
"/home/runner/work/DispersiveShallowWater.jl/DispersiveShallowWater.jl/docs/build/analysis_callback_relaxation.png"
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, which needs to be imported first:
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.
mesh = Mesh1D(coordinates_min, coordinates_max, N)
+accuracy_order = 4
+D_legendre = legendre_derivative_operator(-1.0, 1.0, accuracy_order)
+uniform_mesh = UniformPeriodicMesh1D(mesh.xmin, mesh.xmax, div(mesh.N, accuracy_order))
SummationByPartsOperators.UniformPeriodicMesh1D{Float64} with 128 cells in (-130.0, 20.0)
Upwind DG operators in negative, central and positive operators can be obtained by couple_discontinuously
central = couple_discontinuously(D_legendre, uniform_mesh)
+minus = couple_discontinuously(D_legendre, uniform_mesh, Val(:minus))
+plus = couple_discontinuously(D_legendre, uniform_mesh, Val(:plus))
+D1 = PeriodicUpwindOperators(minus, central, plus)
Upwind SBP first-derivative operators of order 3 on a grid in [-129.99999999999997, 19.999999999999996] using 512 nodes
+and coefficients of Module
In order to still have an entropy-conserving semidiscretization the second-derivative SBP operator needs to be
using SparseArrays: sparse
+D2 = sparse(plus) * sparse(minus)
512×512 SparseArrays.SparseMatrixCSC{Float64, Int64} with 3072 stored entries:
+⎡⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⎤
+⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
+⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
+⎣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦
The Solver object can now be created by passing the two SBP operators to the constructor, which, in turn, can be used to construct a Semidiscretization:
solver = Solver(D1, D2)
+semi = Semidiscretization(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions)
Semidiscretization
+ #spatial dimensions: 1
+ mesh: Mesh1D{Float64} with length 512
+ equations: BBMBBMVariableEquations1D
+ initial condition: initial_condition_shoaling
+ boundary condition: boundary_condition_periodic
As before, we can run the simulation by
analysis_callback = AnalysisCallback(semi; interval = 10,
+ extra_analysis_errors = (:conservation_error,),
+ extra_analysis_integrals = (waterheight_total,
+ velocity, entropy),
+ io = devnull)
+relaxation_callback = RelaxationCallback(invariant = entropy)
+callbacks = CallbackSet(relaxation_callback, analysis_callback)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+ save_everystep = false, callback = callbacks, saveat = saveat)
+anim = @animate for step in 1:length(sol.u)
+ plot(semi => sol, plot_initial = true, conversion = waterheight_total, step = step, xlim = (-50, 20), ylims = [(-0.8, 0.1)])
+end
+gif(anim, "shoaling_solution_dg.gif", fps = 25)
For more details see also the documentation of SummationByPartsOperators.jl
Some more examples sorted by the simulated equations can be found in the examples/ subdirectory. Especially, in examples/svaerd_kalisch_1d/ you can find Julia scripts that solve the SvaerdKalischEquations1D 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.
More examples, especially focussing on plotting, can be found in the scripts create_figures.jl and plot_examples.jl.
