Suggestions on what to implement next

hello everyone!

I have recently (almost a year ago!) started working in the field of physical oceanography - by training I am a computational chemist with a lot of experience in high performance computing. I started this job because in the end "it's all code and maths, right?"

I am quite stubborn so I don't like writing code if I don't understand the underlying physics so I've gone deep into oceanography and the physics/math behind it (oh lord). Then I believe that I don't understand things until I've coded them up. In order to not go crazy I decided to start with a simple 2D shallow water solver over a structured Cartesian grid - simple! The other requirement that I placed myself is that everything will be GPU native.

I've written an HLL Riemann solver with minmod reconstruction. The timestepping is done via strong stablity preserving (SSP) RK2. I also have fully balanced wetting/drying.

I added 3D physics via a hydrostatic and also a non-hydrostatic approach (the non-hydrostatic was a fucking nightmare). I use sigma layers for the hydrostatic with baroclinic forcing.

The non-hydrostatic is done via an awful conjugate gradient Poisson solver with barotropic/baroclinic mode splitting. I spent an entire week getting this working and it was god awful. It is also slow.

All the code is parallelized on GPUs with MPI for multi-node, multi-GPU scaling. So far I've gone up to (an ideal dam break simulation) on up to 32 GPUs and it scales pretty nicely. For real case scenarios it is a bit trickier.

For 3D physics I've implemented salinity transport, a linear equation of state (yay simple), a baroclinic pressure gradient and a pressure Jacobian because I was getting some mass instabilities. I can also add Coriolis, per layer bottom friction, and horizontal diffusion for salinity. The goal of all of these was simulating river plumes.

My only real world test case so far was simulating Sydney Harbour and bits of the Paramatta river (this was still all 2d, so no fresh water intake). I validated the 2D physics by reproducing the tide measurements at the Fort Denison measuring gauge. I was able to simulate an entire month around March from last year and the measurements look good. I can account for forcing from wind, and things like pressure and rain but this run was just tides obtained from the EOT20 model. I've got now the TPXO ones which should make the sim more accurate. This Sydney Harbour simulation is a domain 41km by 25km resolved at a 25x25m grid (around 2.3 million little cells woo). A 50x50m resolution failed to capture some physics so I had to go down. It takes around 24 minutes to simulate 3 days of Sydney Harbour with 2D physics. If I enable 3D layering the compute time multiples a bit, I haven't optimized this code because I wanted salinity to add freshwater intake from the Paramatta river.

So this said and done, I'd like to ask for input from the community on "what is an interesting thing to look at?", for now this solver is constrained for coastal, estuarine simulations because, we just because the shallow water equations will die at non shallow water!

The code is still not public, but will be once it is "stable". It is all written in Fortran.

I am thinking on implementing unstructured grids and do some weird nesting with structured in the simple parts and unstructured in the inside, coastal domains where accuracy is needed. Then deal with the communication at the boundary. But that is for later.