Direct 2D solves with FFT-based solver for rigid lid and implicit free surface

[glwagner]

HydrostaticFreeSurfaceModel has two options for free_surface: ExplicitFreeSurface and ImplicitFreeSurface.

Unfortunately, ImplicitFreeSurface using a PreconditionedConjugateGradientSolver incurs a significant cost (even relative to the cost of a three-dimensional direct non hydrostatic pressure solve!), while ExplicitFreeSurface severely limits time-step with realistic gravitational_acceleration and velocity scales.

For grids that are regular and rectilinear in the horizontal, there is a cure: we can implement a direct solve using the FFT based pressure solver either of the 2D barotropic poisson pressure equation, or for the 2D implicit free surface step.

It will be interesting to compare the total cost of the barotropic integral + 2D solve, which requires integrals to extra the barotropic transport, with the 3D solve used by IncompressibleModel. The vertical FFT is asymptotically more costly than the couple vertical integrals required for a barotropic pressure solve; but on the other hand the FFTs are highly optimized.

[navidcy]

Interesting. How other models usually do this?

[glwagner]

MITgcm I believe uses the same two-dimensional preconditioned conjugate gradient solver for the rigid lid case as for the implicit free surface case.

Many ocean models often have a split explicit method for stepping forward the free surface so there is no elliptic solve. (Killworth et al. 1991 discusses the disadvantages of either implicit free surface or rigid lid for ocean modeling with realistically complex geometries / coastlines.)

This optimization really applies just to hydrostatic models on regular grids (no horizontal stretching). Much of the time I think it would be preferable to use the nonhydrostatic model for this case, since our FFT solver is so fast that the price paid is utterly minor. Yet with an immersed boundary (and perhaps only with a non-grid-fitted immersed boundary), there are some lingering issues that we haven't resolved about whether the FFT solver can be used as is while maintaining mass conservation. The hydrostatic solver uses a vertical integral of the continuity equation and is thus far more straightforward to maintain incompressibility with non-grid-fitted boundaries. So there is a little corner case in which you might want this feature. It's also obviously useful for testing the hydrostatic model.

[jmbeckers]

If the FFT-based solver (on a regular grid) is so fast, maybe there is a way to use it as an (optional) preconditionner when the grid is (weakly) stretched?

[christophernhill]

@jmbeckers yes that could be possible, but it needs some code.

@navidcy @glwagner and @jmbeckers two other related points.

1. I tried something in a gist here ( https://gist.github.com/christophernhill/d1f07cec81c19df9903b18461ba38250 ) that seems to suggest that if we accept we are only really interested in finding a potential whose "difference" (weighted in some way) makes our source terms non-divergent then it may be possible to use FFT approach with a stretched grid? In other words we aren't so interested in the pressure, more the gradient of the pressure.
   I am not sure I believe the result, since it seems like someone would have done it 40+ years ago if this was the case. On the other hand the gist does seem to work, and I can't see what I am missing! I would be interested to chat about whether there is something I am overlooking.

2. the iterative solver should be useful for many problems. If it isn't then it is probably something that needs fixing in the code - more than the approach itself. For a rigid lid, however, it may be slow to converge. Is there a reason to use a rigid lid for hydrostatic problems?

[jmbeckers]

@jmbeckers yes that could be possible, but it needs some code.

Yes I'm fully aware of that. My point is that if you are ready to make some coding to include the FFT-based solver into the shallow water model, maybe the alternative is to include it rather as a preconditionner into the iterative solver (if that is not much more work than the direct implementation). Should the grid be regular, invoking the iterative solver with FFT preconditionner should then converge in one iteration (i.e. equivalent to direkt solver) and hopefully (?) with weakly stretched grids convergence could be fast?
Just speculating ...

[glwagner]

@jmbeckers I think that's a nice idea. If we have an implementation of a direct 2D solve for the implicit free surface step I don't think it would be hard to also use it as a preconditioner in the conjugate gradient algorithm.

@christophernhill still thinking about the stretched FFT solver...

[jmbeckers]

1. I tried something in a gist here ( https://gist.github.com/christophernhill/d1f07cec81c19df9903b18461ba38250 ) that seems to suggest  that if we accept we are only really interested in finding a potential whose "difference" (weighted in some way) makes our source terms non-divergent then it may be possible to use FFT approach with a stretched grid?  In other words we aren't so interested in the pressure, more the gradient of the pressure.
   I am not sure I believe the result, since it seems like someone would have done it 40+ years ago if this was the case. On the other hand the gist does seem to work, and I can't see what I am missing! I would be interested to chat about whether there is something I am overlooking.

That is intriguing. If I understand correctly, the test code indeed provides a divergence-free field by adding a finite-difference to the velocity. I wonder if on a stretched grid the approach makes sure the added force is actually a pressure gradient (i.e. derives from a potential). In other words, maybe on a stretched grid the code tested might lead to energy inconsistencies ?

[christophernhill]

1. I tried something in a gist here ( https://gist.github.com/christophernhill/d1f07cec81c19df9903b18461ba38250 ) that seems to suggest  that if we accept we are only really interested in finding a potential whose "difference" (weighted in some way) makes our source terms non-divergent then it may be possible to use FFT approach with a stretched grid?  In other words we aren't so interested in the pressure, more the gradient of the pressure.
   I am not sure I believe the result, since it seems like someone would have done it 40+ years ago if this was the case. On the other hand the gist does seem to work, and I can't see what I am missing! I would be interested to chat about whether there is something I am overlooking.

That is intriguing. If I understand correctly, the test code indeed provides a divergence-free field by adding a finite-difference to the velocity. I wonder if on a stretched grid the approach makes sure the added force is actually a pressure gradient (i.e. derives from a potential). In other words, maybe on a stretched grid the code tested might lead to energy inconsistencies ?

@jmbeckers I am going to try and experiment a bit more next week - was tinkering with other stuff this last week. I did also remember that I think you can introduce any non-divergent stream function in to the Poisson solution in multiple dimensions and it would still be a solution, so that may also cause some trouble in practice and may be another manifestation of possible energetic inconsistencies!

[glwagner]

We have a direct implicit free surface solver now.