# Shallow water equations with DifferentialEquations.jl
**URL:** https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691
**Category:** Numerics
**Tags:** question, diffeq
**Created:** [March 15, 2017, 11:01pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691 "2017-03-15T23:01:46Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [March 15, 2017, 11:01pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/1 "2017-03-15T23:01:46Z")
I am really interested in using the `DifferentialEquations.jl` framework to solve the [shallow water equations](https://en.wikipedia.org/wiki/Shallow_water_equations) over a 2D domain efficiently instead of coding my own solution.
How would you formulate the problem in `DifferentialEquations.jl` ignoring all frictional, viscous and Coriolis terms, and only considering gravity? The corresponding system of equations for this set of assumptions in conservative form is in the Wikipedia link above.
In terms of the implementation: given a 2D map (or Julia array) of the landscape `H` and a map indicating the depth of water `h_t` on top of the landscape at time `t`, I would like to advance in time and obtain `h_{t+1}` iteratively until reach a time span [0,T].
I have a finite volume solution that is working for now, but I am not happy reinventing the wheel. It would be great to compare it with the solution obtained with `DifferentialEquations.jl` and use that instead.
@ChrisRackauckas can you give a hand? 🙂
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 16, 2017, 12:31am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/2 "2017-03-16T00:31:35Z")
I am not familiar with this PDE, so if you could fill me in on algorithms usually used to solve these?
Normally, they correspond directly to some spatial + timestepper method. The way this is normally done, what you do is the spatial discretization part (discretize every spatial derivative, and the variables) and turn it into an ODE, and then send that ODE to some DifferentialEquations.jl solver. For example: “Crank-Nicholson” is just doing a spatial discretization and then using the `Trapezoid()` method on the resulting ODEs.
However, you need to make sure that the spatial discretization of the PDE is sensible, which is why you should check the standard PDE methods first.
(Also, there are some people interested in PDE discretizations for GSoC. So share some references on this PDE, and I’ll add it to the list of possible PDEs to tackle 🙂 )
---
### Author:  [@mauro3](https://discourse.julialang.org/u/mauro3)
#### Post date: [March 16, 2017, 8:15am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/3 "2017-03-16T08:15:18Z")
It being a hyperbolic problem the FVM should be well suited. Once you got the spatial discretisation, then using an explicit ODE time stepper should work best. The Leveque book suggests: at least second order and not a multistep method.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 16, 2017, 8:22am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/4 "2017-03-16T08:22:47Z")
> [@mauro3](#):
>
> The Leveque book suggests
What’s the full name of this? (I should probably become more familiar with fluid dynamics PDEs…)
---
### Author:  [@mauro3](https://discourse.julialang.org/u/mauro3)
#### Post date: [March 16, 2017, 8:26am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/5 "2017-03-16T08:26:03Z")
Here: [Finite Volume Methods for Hyperbolic Problems - Randall J. LeVeque - Google Books](https://books.google.ch/books?id=O_ZjpMSZiwoC&dq=Finite+volume+methods+leveque&lr=&source=gbs_navlinks_s) . It’s a nice book, but I couldn’t say I’ve read it all.
It goes with the Clawpack software: [http://www.clawpack.org/](http://www.clawpack.org/)
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [March 16, 2017, 4:24pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/6 "2017-03-16T16:24:08Z")
Hi @ChrisRackauckas,
I appreciate if you can skim through this quick overview on FVM for the equations: [http://caminos.udc.es/info/asignaturas/201/Finite%20Volumes%20IMWE.pdf](http://caminos.udc.es/info/asignaturas/201/Finite%20Volumes%20IMWE.pdf)
What I have currently for space discretizattion is a backward difference in space that propagates the fluxes properly in the direction of the velocity field (it is called the upwind scheme) given that the time step is small enough (satisfies the CLF condition described in the document).
I think my implementation is working fine, but I would like to be able to switch from one PDE to another quickly in the future and play with different solvers, in other words, enjoy the goodies of `DifferentialEquations.jl` 🙂
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [March 16, 2017, 4:24pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/7 "2017-03-16T16:24:59Z")
Thanks @mauro3 for the reference, super valuable 🙂
---
### Author:  [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)
#### Post date: [March 23, 2017, 6:33pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/8 "2017-03-23T18:33:06Z")
@juliohm: Have a look at [http://dedalus-project.org/](http://dedalus-project.org/) . Dedalus is a Python-based system that takes a PDE in symbolic form, discretizes it into a time-stepping system, and then does the time-stepping with standard C libraries (FFTW, LAPACK, MPI). It’s not Julia, but it might do what you want with minimal effort on your part.
BTW, I am also quite interested in seeing how Julia works for PDE simulation, in my case Navier-Stokes in simple geometries (pipes, channels, cubes), with distributed-memory parallelization. I hope to implement a Fourier-Chebyshev channel-flow code or isotropic turbulence in a triply periodic cube. But I believe some fundamental work needs to happen to enable this, namely getting distributed-memory FFTW working in Julia and coordinated with DistributedArrays.
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [March 23, 2017, 6:41pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/9 "2017-03-23T18:41:54Z")
See also [ANUGA](https://anuga.anu.edu.au/) which is a FVM to model Tsunami inundations (and more). It’s a Python program that resorts to C to do the heavy work.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 23, 2017, 6:53pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/10 "2017-03-23T18:53:35Z")
> [@John\_Gibson](#):
>
> @juliohm: Have a look at [http://dedalus-project.org/](http://dedalus-project.org/) . Dedalus is a Python-based system that takes a PDE in symbolic form, discretizes it into a time-stepping system, and then does the time-stepping with standard C libraries (FFTW, LAPACK, MPI). It’s not Julia, but it might do what you want with minimal effort on your part.
That looks similar to some of what the ApproxFun+DifferentialEquations setup has going on (well, is working towards). I think the right way to their kind of functionality is to make a distributed `Fun` type which uses what you’re talking about, and then pairing that with the timestepping routines on `Fun` types. The timestepping part is coming along just fine. But maybe @dlfivefifty has something to say about possible parallelism of `Fun`s?
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [March 23, 2017, 7:53pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/11 "2017-03-23T19:53:22Z")
> [@ChrisRackauckas](#):
>
> That looks similar to some of what the ApproxFun+DifferentialEquations setup has going on (well, is working towards). I think the right way to their kind of functionality is to make a distributed Fun type which uses what you’re talking about, and then pairing that with the timestepping routines on Fun types. The timestepping part is coming along just fine. But maybe @dlfivefifty has something to say about possible parallelism of Funs?
Dedalus is a great project (a colleague of mine Geoff Vasil is one of the people behind it) that discretises in a very similar way to ApproxFun: using Chebyshev polynomials and their derivative. Hopefully in a few years ApproxFun + DifferentialEquations will be able to compete.
As for parallelism: the primary parallelism in Dedalus is to decouple the Fourier modes (assuming you are using the Fourier basis in one of the directions) to handle the linear, implicit part of the equation. This is possible since Fourier gives you diagonal differential operators. Then parallelised FFT and transposes allow you to handle the nonlinear, explicit part by transforming back to value space.
This is potentially doable in ApproxFun. The main part would be to recognise that a partial differential operator is diagonal in one dimension and then auto-parrallelise the solve. This would all happen in the initial factorisation of the linear part (so perhaps a flag `qrfact(A;parallelize=true)`).
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [March 23, 2017, 9:59pm UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/12 "2017-03-23T21:59:21Z")
@John_Gibson thanks for sharing, when it comes to Python, there are plenty of packages. First time I heard of Dedalus, bookmarking it here. I think the most famous in Python is the Fenics project:
[https://fenicsproject.org/](https://fenicsproject.org/)
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [May 10, 2019, 4:31am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/13 "2019-05-10T04:31:50Z")
Hey @juliohm. In the end did you manage to use DiffEq.jl?
I am also looking at dedalus or fenics to work with the Shallow Water system.
Did DiffEq + ApproxFun evolve nicely?
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [May 10, 2019, 5:19am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/14 "2019-05-10T05:19:22Z")
There hasn’t been much work on DiffEq + ApproxFun, though @ChrisRackauckas did demonstrate some examples at JuliaCon 2018, his talk is up on YouTube.
It’s a project just waiting for the right PHD student, really.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 10, 2019, 6:51am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/15 "2019-05-10T06:51:21Z")
Well the issue is still in broadcast, but now DiffEq is only using broadcast internally (yay [http://juliadiffeq.org/2019/05/09/GPU.html](http://juliadiffeq.org/2019/05/09/GPU.html)) so the MWE is now pretty simple:
```julia
using ApproxFun
x = Fun(identity,-1..1)
f = exp(x)
g = f/sqrt(2-x^2)
h = @. f + 0.1g
length(f.coefficients) #18
length(g.coefficients) # 51
length(h.coefficients) # 39
ll = deepcopy(f.coefficients)
resize!(ll,51)
_h = Fun(f.space,ll .+ 0.1.*g.coefficients)
__h = f + 0.1g
using Plots
plot(_h)
plot!(__h)
length(__h.coefficients) #51
```
While it was helpful in this case, it’s not in all cases:
```julia
f = exp(x)
g = f/sqrt(2-x^2)
h = @. f + 0.1g # boom
```
This behavior of refitting in each broadcast can really cause the coefficient sizes to grow, when in reality it doesn’t need to be more than the maximum of the coefficients of any other Fun. This was the issue pre-v1.0, and now it seems to be more contained. Now I think we just need a way to say “prune coefficients at some constant maximum number of coefficients, or to a tolerance of epsilon (equal to the ODE solver tolerance)”.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 10, 2019, 7:19am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/16 "2019-05-10T07:19:09Z")
A better display is then this integration example:
```julia
using ApproxFun, OrdinaryDiffEq
x = Fun(identity,-1..1)
f = exp(x)
g = f/sqrt(2-x^2)
h = @. f + 0.1g
function ff(u,p,t)
2u - u^2
end
prob = ODEProblem(ff,h,(0.0,1.0))
OrdinaryDiffEq.recursive_bottom_eltype(a::ApproxFun.Fun) = eltype(a.coefficients)
OrdinaryDiffEq.recursive_unitless_bottom_eltype(a::ApproxFun.Fun) = eltype(a.coefficients)
OrdinaryDiffEq.recursive_unitless_eltype(a::ApproxFun.Fun) = eltype(a.coefficients)
isnan(h) # :(
Base.isnan(a::ApproxFun.Fun) = any(isnan,a.coefficients)
sol = DiffEqBase.__solve(prob,Tsit5(),adaptive=false,dt=0.1)
@show sol[end].coefficients
sol = DiffEqBase.__solve(prob,Tsit5())
@show sol[end].coefficients
```
So just a few reasonable overloads (could be added via Requires) and that’s probably the simplest adaptive time adaptive space semilinear PDE solver script you’ve ever seen 🙂. The coefficients lists look like:
```julia
x = Fun(Chebyshev(-1..1),[1.37369, 1.23314, 0.312644, 0.057217, 0.0102633, 0.00210412, 0.000680313, 0.000220778, 9.27862e-5, 3.23932e-5, 1.40827e-5, 4.98013e-6, 2.19425e-6, 7.81442e-7, 3.47373e-7, 1.24314e-7, 5.56183e-8, 1.99758e-8, 8.98101e-9, 3.23453e-9, 1.45984e-9, 5.26919e-10, 2.38555e-10, 8.62588e-11, 3.91528e-11, 1.41782e-11, 6.44926e-12, 2.33834e-12, 1.06549e-12, 3.86616e-13, 1.76285e-13,
6.3722e-14, 2.87327e-14, 1.01941e-14, 4.5861e-15, 1.61385e-15, 6.88532e-16, 2.68315e-16, 1.32806e-16])
(sol[end]).coefficients = [1.74261, 0.397775, -0.0454502, -0.00270729, 0.00176207, -0.000219944, 6.33866e-5, -1.62488e-5, 1.19699e-5, -3.06074e-6, 1.83144e-6, -4.60346e-7, 2.83213e-7, -7.24161e-8, 4.46337e-8, -1.15028e-8, 7.09702e-9, -1.8429e-9, 1.13747e-9, -2.97332e-10, 1.83457e-10,
-4.82406e-11, 2.97425e-11, -7.863e-12, 4.84288e-12, -1.28656e-12, 7.91322e-13, -2.11464e-13, 1.29487e-13, -3.47097e-14, 2.12311e-14, -5.75179e-15, 3.54605e-15, -1.21433e-15, 9.9087e-16, -6.27891e-16, 2.29117e-16, -2.48411e-16, 1.45315e-16, 4.15729e-16, -7.54313e-18, -1.87555e-18, -5.55151e-19, -1.89897e-19, -7.16383e-20, -2.71441e-20, -1.05819e-20, -4.04815e-21, -1.58455e-21, -6.07215e-22, -2.3798e-22, -9.11932e-23, -3.5739e-23, -1.36772e-23, -5.35425e-24, -2.04399e-24, -7.98466e-25, -3.03684e-25, -1.18248e-25, -4.47412e-26, -1.73416e-26, -6.51539e-27, -2.5094e-27, -9.33774e-28, -3.56495e-28, -1.309e-28, -4.93575e-29, -1.77803e-29, -6.58306e-30, -2.30334e-30, -8.28578e-31, -2.77257e-31, -9.58546e-32, -2.97471e-32, -9.46569e-33, -2.61974e-33, -6.36441e-34, 1.57045e-35, 5.02227e-36, 1.70695e-36, 6.18296e-37, 2.30489e-37, 8.70625e-38, 3.28727e-38, 1.24436e-38, 4.69274e-39, 1.77251e-39, 6.6648e-40, 2.50898e-40, 9.39579e-41, 3.52145e-41, 1.31187e-41, 4.88924e-42, 1.80954e-42, 6.69697e-43, 2.45846e-43, 9.01968e-44, 3.27759e-44, 1.18944e-44, 4.26708e-45, 1.52726e-45, 5.3895e-46, 1.89486e-46, 6.54383e-47, 2.24694e-47, 7.53617e-48, 2.50516e-48, 8.0617e-49, 2.55813e-49, 7.74247e-50, 2.28975e-50, 6.23784e-51, 1.61996e-51, 3.54841e-52, 5.93231e-53]
(sol[end]).coefficients = [1.74262, 0.397778, -0.0454475, -0.00270537, 0.00176325, -0.000219333, 6.35796e-5, -1.62111e-5, 1.19832e-5, -3.05395e-6, 1.83367e-6, -4.59648e-7, 2.83495e-7, -7.23072e-8, 4.46756e-8, -1.14851e-8, 7.10432e-9, -1.84003e-9, 1.13863e-9, -2.96853e-10, 1.83652e-10, -4.81613e-11, 2.97752e-11, -7.85006e-12, 4.84817e-12, -1.28459e-12, 7.92376e-13, -2.11071e-13, 1.29953e-13, -3.44992e-14, 2.15547e-14, -5.61981e-15, 3.54937e-15, -1.37966e-15, 2.62586e-16, 3.8377e-17, 1.03977e-16, -3.55172e-17, -4.26953e-16, -7.84177e-17, 1.85717e-16, -1.83744e-16, -5.06391e-19, -1.732e-19, -6.53278e-20, -2.4752e-20, -9.64887e-21, -3.69133e-21, -1.44485e-21, -5.53705e-22, -2.17004e-22, -8.31592e-23, -3.25898e-23, -1.24727e-23, -4.8826e-24, -1.86403e-24, -7.28156e-25, -2.76957e-25, -1.07839e-25, -4.08055e-26, -1.58159e-26, -5.94256e-27, -2.28874e-27, -8.51733e-28, -3.2517e-28, -1.19409e-28, -4.50246e-29, -1.62213e-29, -6.00592e-30, -2.10174e-30, -7.56086e-31, -2.53061e-31, -8.74967e-32, -2.71649e-32, -8.64636e-33,
-2.39494e-33, -5.82687e-34, 1.3202e-35, 4.2246e-36, 1.43632e-36, 5.20266e-37, 1.93934e-37, 7.32492e-38, 2.76571e-38, 1.04692e-38, 3.94828e-39, 1.49132e-39, 5.60774e-40, 2.11105e-40, 7.90595e-41, 2.96308e-41, 1.1039e-41, 4.11417e-42, 1.52275e-42, 5.63557e-43, 2.06892e-43, 7.59046e-44,
2.75837e-44, 1.00101e-44, 3.59124e-45, 1.28535e-45, 4.53604e-46, 1.59476e-46, 5.5077e-47, 1.89109e-47, 6.34292e-48, 2.10837e-48, 6.78502e-49, 2.15277e-49, 6.51563e-50, 1.9265e-50, 5.24791e-51, 1.36221e-51, 2.98239e-52, 4.97929e-53]
```
and that’s what I mean by the cutoff issue, since the ODE solver tolerance is just a relative tolerance of 1e-3 and abstol of 1e-6, so it should really look like:
```julia
(sol[end]).coefficients = [1.74262, 0.397778, -0.0454475, -0.00270537, 0.00176325, -0.000219333, 6.35796e-5, -1.62111e-5, 1.19832e-5, -3.05395e-6, 1.83367e-6]
```
so if there was a constructor that could autochop:
```julia
function ff(u,p,t)
chop!(u,1e-6)
2u - u^2
end
prob = ODEProblem(ff,h,(0.0,1.0))
sol = DiffEqBase.__solve(prob,Tsit5())
@show sol[end].coefficients
```
this would have an efficient solution
```julia
(sol[end]).coefficients = [1.74262, 0.397778, -0.0454475, -0.00270537, 0.00176325, -0.000219333, 6.35794e-5, -1.62109e-5, 1.19852e-5, -3.06974e-6, 1.88027e-6, -4.7376e-7,
3.15883e-7]
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 10, 2019, 7:27am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/17 "2019-05-10T07:27:19Z")
Oh and BTW, implicit methods via iteration work:
```julia
sol = DiffEqBase.__solve(prob,KenCarp4(nlsolve=NLFunctional()))
@show sol[end].coefficients
sol = DiffEqBase.__solve(prob,KenCarp4(nlsolve=NLAnderson()))
sol = DiffEqBase.__solve(prob,KenCarp4(nlsolve=NLNewton()))
```
but using quasi-Newton will fail because ApproxFun doesn’t expose an array interface on the coefficients:
```julia
h[2] # h.coefficients[2]?
length(h) # length(h.coefficients)?
```
if that was exposed, then the nonlinear solver would automatically work like a pseudospectral method by solving the ODE update equation in coefficient space, assuming the Jacobian calculation plays nicely. So I just need to coordinate with @dlfivefifty to get these things acting enough like arrays for them to work, or enough like scalars for them to work.
---
### Author:  [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)
#### Post date: [May 10, 2019, 7:36am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/18 "2019-05-10T07:36:34Z")
The long term goal is to make ApproxFun a thin wrapper around a continuum array: [https://github.com/JuliaApproximation/ContinuumArrays.jl](https://github.com/JuliaApproximation/ContinuumArrays.jl)
Then all the Operator terminology becomes (possibly infinite via InfiniteArrays.jl) Array terminology, making it a lot easier to describe basic operations like restrict to a subspace, mass matrices, etc. In particular, it gives a nice language for generating structured and parallelised arrays.
We have the proof of concept implemented for linear spline FEM and (infinite) ultraspherical methods, see the tests. Once this is done syncing with DiffEq should be trivial.
However, I’m busy working on papers at the moment so am actively trying to avoid working on ContinuumArrays.jl. But perhaps I’ll have a software relapse and show that it can be linked in with DiffEq.jl.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 10, 2019, 7:36am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/19 "2019-05-10T07:36:44Z")
L-stable integrators will work with the right derivative definition:
```julia
OrdinaryDiffEq.jacobian(ff,x::ApproxFun.Fun,diff_type) = 0.0 # should be dff/du
OrdinaryDiffEq.derivative(ff,x::Number,integrator) = 0.0 # should be dff/dt
OrdinaryDiffEq._vec(x::ApproxFun.Fun) = x
OrdinaryDiffEq._reshape(x::ApproxFun.Fun,siz) = x
sol = DiffEqBase.__solve(prob,Rosenbrock23())
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 10, 2019, 7:37am UTC](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691/20 "2019-05-10T07:37:37Z")
> [@dlfivefifty](#):
>
> The long term goal is to make ApproxFun a thin wrapper around a continuum array: [GitHub - JuliaApproximation/ContinuumArrays.jl: A package for representing quasi arrays with continuous indices](https://github.com/JuliaApproximation/ContinuumArrays.jl)
>
> Then all the Operator terminology becomes (possibly infinite via InfiniteArrays.jl) Array terminology, making it a lot easier to describe basic operations like restrict to a subspace, mass matrices, etc. In particular, it gives a nice language for generating structured and parallelised arrays.
>
> We have the proof of concept implemented for linear spline FEM and (infinite) ultraspherical methods, see the tests. Once this is done syncing with DiffEq should be trivial.
Yeah that sounds good, so I can wait 🙂
[Next page](https://discourse.julialang.org/t/shallow-water-equations-with-differentialequations-jl/2691.md?page=2)