Remove shift in FFTBasedPoissonSolver and multiply by negative one in preconditioning for Conjugate GradientPoissonSolver

[xkykai]

Closes #4581, and details of the origin of this can also be found #4539.

@glwagner @tomchor @amontoison

[tomchor]

For the record, it seems that running the code below (which fails on main in the sense that the velocities start blowing up) using this branch plus #4554 works:

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),
                             )
u₀(x, z) = U + 1e-2*rand()
set!(model, u=u₀, enforce_incompressibility=true)

Produces:

julia> model.velocities.u
9×1×8 Field{Face, Center, Center} on ImmersedBoundaryGrid on CPU
├── grid: 8×1×8 ImmersedBoundaryGrid{Float64, Bounded, Flat, Bounded} on CPU with 4×0×4 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Open{Oceananigans.BoundaryConditions.PerturbationAdvection{Float64}}, east: Open{Oceananigans.BoundaryConditions.PerturbationAdvection{Float64}}, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 17×1×16 OffsetArray(::Array{Float64, 3}, -3:13, 1:1, -3:12) with eltype Float64 with indices -3:13×1:1×-3:12
    └── max=1.00883, min=0.996508, mean=1.0023

and it converges in only 7 iterations. So I think we're on the right track!

[xkykai]

Yes I was going to report on this as well! @amontoison weirdly enough the tiny shift in the preconditioner causes a significant blowup of velocities in this MWE. I have little familiarity of open boundary conditions and its implementation, so I can't comment on why that would be the case.

#4581 (comment)

[amontoison]

@xkykai let's drop the shift in that case. If we want to ensure positive-definiteness of the preconditioner, we probably need to add a SPD matrix related to the physic application.
But we should investigate it only if we don't have good numerical results.
It seems to work well now.