Shallow water equations with DifferentialEquations.jl



I am really interested in using the DifferentialEquations.jl framework to solve the 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? :slight_smile:


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 :slight_smile: )


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.


What’s the full name of this? (I should probably become more familiar with fluid dynamics PDEs…)


Here: . It’s a nice book, but I couldn’t say I’ve read it all.

It goes with the Clawpack software:


Hi @ChrisRackauckas,

I appreciate if you can skim through this quick overview on FVM for the equations:

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 :slight_smile:


Thanks @mauro3 for the reference, super valuable :slight_smile:


@juliohm: Have a look at . 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.


See also ANUGA which is a FVM to model Tsunami inundations (and more). It’s a Python program that resorts to C to do the heavy work.


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)).


@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: