A significant computational bottleneck for some problems with NonhydrostaticModel (mainly complex domains and multi-GPU situations) is the cost to solve the pressure Poisson equation.

We've long been aware of Le and Moin 1991, which claims to allow us to reduce the number of Poisson solves to 1 per RK3 time-step. I can't really remember why we haven't pursued this --- perhaps some skepticism that it would work? cc @simone-silvestri @xkykai @tomchor

More recently, an algorithm deemed (FastRK3) has been defined by Aithal and Ferrante (2020), which is similar (but perhaps more clearly explained than Le and Moin 1991).

The gist is that FastRK3 modifies the momentum tendencies to use a pressure "estimate" rather than the exact pressure for the first two substeps. It requires storing an additional variable (the pressure from the previous timestep), and also starting time-stepping with a normal RK3 step (ie the method is not self-starting). But all of that is not too hard...

I think we should consider implementing FastRK3, especially for multi-GPU.