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

[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 https://github.com/CliMA/Oceananigans.jl/discussions/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 https://github.com/CliMA/Oceananigans.jl/discussions/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 :joy:

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.