Add tests for `ConjugateGradientPoissonSolver` with `PerturbationAdvectionOpenBoundaryConditions` in `ImmersedBoundaryGrid`s

[tomchor]

Originally this PR did what is described below. Then https://github.com/CliMA/Oceananigans.jl/pull/4582 and https://github.com/CliMA/Oceananigans.jl/pull/4554 were merged, which did part of what was originally proposed here. Now this PR will only add tests to ensure that things keep working with the ConjugateGradientPoissonSolver when using PerturbationAdvectionOpenBoundaryConditions in ImmersedBoundaryGrids. At least in simple cases.

---

The goal of this PR is to fix the initialization of PerturbationAdvectionOpenBoundaryCondition. I was playing around with it and did notice that the indices at i=1 weren't receiving the external boundary condition (the one we get through getbc()) which led to, so this PR tests a solution for that. For now this fix is just an ifelse statement, but maybe we can find a better one.

I hadn't opened a PR earlier because I was trying to figure out why the implementation of the open boundary condition was different on for left boundary than for the right boundaries (right boundaries seem to work fine). It was unclear to be where the asymmetry was coming from. Then @jagoosw's comment comment earlier mentioned that the left boundary in stepped but apparently the right boundary is not.

As a benchmark, this MWE is currently failing on main, but is running on this branch:

using Oceananigans
using Oceananigans.BoundaryConditions: PerturbationAdvectionOpenBoundaryCondition
using Oceananigans.Solvers: ConjugateGradientPoissonSolver

Lx, Lz = 10, 3
Nx = Nz = 8

grid_base = RectilinearGrid(topology = (Bounded, Flat, Bounded), size = (Nx, Nz), x = (0, Lx), z = (0, Lz))
flat_bottom(x) = 1
grid = ImmersedBoundaryGrid(grid_base, PartialCellBottom(flat_bottom))
U = 1
inflow_timescale = 0.0
outflow_timescale = Inf
u_boundaries = FieldBoundaryConditions(
    west   = PerturbationAdvectionOpenBoundaryCondition(U; inflow_timescale, outflow_timescale),
    east   = PerturbationAdvectionOpenBoundaryCondition(U; inflow_timescale, outflow_timescale),)
    
model = NonhydrostaticModel(; grid,
                              boundary_conditions = (; u = u_boundaries),
                              pressure_solver = ConjugateGradientPoissonSolver(grid, maxiter=10),
                              advection = WENO(; grid, order=5)
                             )
u₀(x, z) = U + 1e-2*rand()
set!(model, u=u₀, enforce_incompressibility=true)

For now this PR also reverts the FFT preconditioner to what it was before https://github.com/CliMA/Oceananigans.jl/pull/4041 because it makes it much easier for simulations with both ImmersedBoundaries and OpenBoundaryConditions to run. Although I guess ideally we wouldn't need to do that since the change is needed for Krylov methods (at least that's what I understood).

Ultimately, I think we should also add a test that catches this in the future.

Related to https://github.com/CliMA/Oceananigans.jl/issues/4451#issuecomment-2903837583

Closes https://github.com/CliMA/Oceananigans.jl/issues/4603

[tomchor]

Also, can someone explain why left open boundary conditions are time-stepped, but not right open BCs? Should we change that so that the implementation of open boundaries matching schemes can be exactly the same for a left boundary and for a right boundary? Then the ad-hoc fix in this PR wouldn't be needed.

[tomchor]

cc @glwagner

[jagoosw]

The left boundary for a face field is part of the "interior", I'm not sure it is actually every necessary to step it so the best solution would probably be to work out how to do that.

I looked through the file changes but I'm not seeing what the fix is here? Also, just to let you know the mwe did work for me with the initialisation in my comment before with ConjugateGradientPoissonSolver.

[tomchor]

The left boundary for a face field is part of the "interior", I'm not sure it is actually every necessary to step it so the best solution would probably be to work out how to do that.

What do you mean by interior here? The right boundary for a face field also appears in interior().

I looked through the file changes but I'm not seeing what the fix is here?

I linked you in a code comment pointing to the fix. Let me know if that's not clear.

Also, just to let you know the mwe did work for me with the initialisation in my comment before with ConjugateGradientPoissonSolver.

Can you try with the MWE here? I started a fresh Manifest with the current main and it produces NaNs.

[tomchor]

I ran many tests today, and I found that the only way that I could somewhat reliably run simulations was to switch the preconditioner to what it was before https://github.com/CliMA/Oceananigans.jl/pull/4041.

So I propose we do that in this PR and add tests to ensure that open boundary conditions + immersed boundaries keep working.

@jagoosw @glwagner @amontoison thoughts?

[jagoosw]

I ran the MWE, and it does not work with the removal of the "secret" boundary storage. If I revert to the commit just before you changed that, then the MWE works as is.

[jagoosw]

Oh wait maybe it works with the latest commit but I definitely do not think it is correct and it will be wrong in some scenarios e.g. when the inflow timescale is not zero.

[tomchor]

Oh wait maybe it works with the latest commit but I definitely do not think it is correct and it will be wrong in some scenarios e.g. when the inflow timescale is not zero.

I'm happy to revert that! But for the record, I tried to break it in many different configurations and every single time there was no visible difference between the original scheme, and the one which ignores the secret indices. For example:

https://github.com/user-attachments/assets/2009a808-1db7-48d6-9093-6db3538b4cd0

https://github.com/user-attachments/assets/d6f19944-30ab-4cd7-b1ee-ba08eb3928e7

(tbf, there is one visible difference which is the first time-step.)

[tomchor]

I think this PR is ready for review. This makes simulations with PerturbationAdvection OBC run as well as they did before https://github.com/CliMA/Oceananigans.jl/pull/4041. The key to make that happen is to revert precondition! to what it was before. In addition to that, this PR also:

1. Makes improvements to PerturbationAdvection, including better initialization.
2. Adds several tests to test_condjugate_gradient_poisson_solver.jl, including one that tests if the simulation runs with IBG + PerturbationAdvection.
3. Cleans up and improves relevant validation scripts

@glwagner @jagoosw I didn't change how the previous timestep is stored for left-side BCs.