Add ImmersedPressureSolver with FFTBasedPoissonSolver as preconditioner for complex domains

[xkykai]

Using the FFTBasedPoissonSolver in complex domains on the ImmersedBoundaryGrid leads to resulting solution that either

1. Does not conserve tracers in the flow; or
2. Leads to divergent flow along the immersed boundary.

On the other hand, using the PreconditionedConjugateGradientSolver on complex domains leads to exact solution at a much slower speed on the GPU. This PR implements an ImmersedPressureSolver that uses the FFTBasedPoissonSolver as a preconditioner, then iterates PreconditionedConjugateGradientSolver to arrive at the correct solution. This approach is much faster than using a vanilla PreconditionedConjugateGradientSolver as the number of iterations scale much more slowly with grid size on the GPU.

This PR will include the ImmersedPressureSolver, a few validation examples for benchmarking purposes and a number of example use cases where this solver would be useful.

@glwagner @simone-silvestri

[navidcy]

can we name this something bit more descriptive and less referential to MITgcm?

[xkykai]

Something like SymmetricLocalPreconditioner? I'm not entirely sure, could use some ideas too!

[navidcy]

Well I that sounds much better because at least it tries to describe something and not just refer to some other model.

I'm not sure what this preconditioner is doing so that's why I'm not in position to propose something. I don't even know what local refers to. But definitely a SomethingThatDescribesMePreconditioner is better than SimilarToSomeOtherModelOrSomebodyElsesPreconditioner. That's what I wanted to convey with my comment.

[glwagner]

I believe this preconditioner is some kind of perturbation about a diagonal preconditioner, is that right? DiagonallyDominantPreconditioner or something. In fact shoudln't we also delete the wrong version that's already present @simone-silvestri ?

[simone-silvestri]

Exactly, it should be the same one but corrected to be 3D and with the correct denominator (we should correct also the 2D version of it)

[glwagner]

Ok so maybe DiagonallyDominantThreeDimensionalPreconditioner and also TwoDimensional.

[tomchor]

I'm not sure what the state of this PR is, but last I checked it was working, but was about 5 times slower than the FFT-based solver, no?

Coming from the discussion in #3529, I'd argue that it'd be useful to include this solver in the main code, even if it's a bit slow for now. Then it'd be easy to check the influence of near-immersed-boundary divergence when investigating results such as #3529.

Furthermore, I can see situations where I'd be willing to spend 5 times longer running simulations in order to be more sure of the results. Oceananigans can be so fast on GPUs that a 5x slower code might still be as good as a more traditional CPU code running the same problem 😂

What do you guys think?

[glwagner]

I'm not sure what the state of this PR is, but last I checked it was working, but was about 5 times slower than the FFT-based solver, no?

Coming from the discussion in #3529, I'd argue that it'd be useful to include this solver in the main code, even if it's a bit slow for now. Then it'd be easy to check the influence of near-immersed-boundary divergence when investigating results such as #3529.

Furthermore, I can see situations where I'd be willing to spend 5 times longer running simulations in order to be more sure of the results. Oceananigans is can be so fast on GPUs that a 5x slower code might still be as good as a more traditional CPU code running the same problem 😂

What do you guys think?

There's no reason to wait on merging this as far as I'm concerned. @xkykai and @simone-silvestri can perhaps comment on their strategy / plan.