Skip to content

ConjugateGradientPoissonSolver's FFTBasedPoissonSolver doesn't need to flip sign since it's already semi-positive definite #4581

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Open
xkykai opened this issue Jun 5, 2025 · 4 comments · May be fixed by #4582
Open

ConjugateGradientPoissonSolver's FFTBasedPoissonSolver doesn't need to flip sign since it's already semi-positive definite #4581

xkykai opened this issue Jun 5, 2025 · 4 comments · May be fixed by #4582
Labels
immersed boundaries ⛰️ Less Ocean, more anigans numerics 🧮 So things don't blow up and boil the lobsters alive

Comments

@xkykai
Copy link
Collaborator

xkykai commented Jun 5, 2025 •
edited
Loading

As discussed with @glwagner and @tomchor and outlined in #4539, reverting

Oceananigans.jl/src/Solvers/conjugate_gradient_poisson_solver.jl

Line 132 in 51a77cc

p .*= -1
into its pre- #4041 produces improvement when open boundary conditions are involved. Furthermore, the FFTBasedPoissonSolver solves $\nabla^2 \phi = b$ rather than $-\nabla^2 \phi = b$. Therefore the operator of FFTBasedPoissonSolver is (semi-) positive definite. Here's a script to illustrate this point:

using Oceananigans
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Solvers: FourierTridiagonalPoissonSolver, FFTBasedPoissonSolver, solve!
using Oceananigans.Operators
using Oceananigans.Utils: launch!
using Statistics
using Random
using KernelAbstractions: @kernel, @index
using Oceananigans.Architectures: architecture

arch = CPU()
Nx = Nz = 64
topology = (Periodic, Flat, Bounded)

rng = Xoshiro(123)
x = z = (0, 1)

grid  = RectilinearGrid(arch; x, z, topology, size=(Nx, Nz))

lhs = CenterField(grid)

solver = FFTBasedPoissonSolver(grid)

rhs = CenterField(grid)
set!(rhs, (x, z) -> rand())
rhs .-= mean(rhs)

solver.storage .= interior(rhs) # Copy rhs to storage

solve!(lhs, solver)

@kernel function laplacian!(∇²ϕ, grid, ϕ)
    i, j, k = @index(Global, NTuple)
    @inbounds ∇²ϕ[i, j, k] = ∇²ᶜᶜᶜ(i, j, k, grid, ϕ)
end

function compute_laplacian!(∇²ϕ, ϕ)
    grid = ϕ.grid
    arch = architecture(grid)
    fill_halo_regions!(ϕ)
    launch!(arch, grid, :xyz, laplacian!, ∇²ϕ, grid, ϕ)
    return nothing
end

rhs_∇² = CenterField(grid)
compute_laplacian!(rhs_∇², lhs)

@info sum((interior(rhs_∇²) .≈ interior(rhs))[:, 1, :]) == Nx * Nz

which outputs true, which illustrate that FFTBasedPoissonSolver does indeed invert ∇²ᶜᶜᶜ(i, j, k, grid, ϕ).

Beyond this we should also remove the shift in the preconditioner

Oceananigans.jl/src/Solvers/conjugate_gradient_poisson_solver.jl

Line 131 in 51a77cc

solve!(p, preconditioner, preconditioner.storage, shift)
since it improves performance for open boundary conditions according to @tomchor, and the (immersed-aware) Laplacian operator is already (semi-) positive definite anyway.

@amontoison
@xkykai xkykai added numerics 🧮 So things don't blow up and boil the lobsters alive immersed boundaries ⛰️ Less Ocean, more anigans labels Jun 5, 2025
@xkykai xkykai linked a pull request Jun 5, 2025 that will close this issue
Open
@glwagner
Copy link
Member

glwagner commented Jun 5, 2025

(or make the "shift" an optional argument when constructing the preconditioner perhaps)

@xkykai
Copy link
Collaborator Author

xkykai commented Jun 5, 2025 •
edited
Loading

(or make the "shift" an optional argument when constructing the preconditioner perhaps)

Do you think this implementation should be (yet another - oh man!) PR that updates the API for both FFTBasedPoissonSolver and FourierTridiagonalPoissonSolver to move the shift argument to the solver construction level? In my opinion this should be done after #4580 is sorted out

@glwagner
Copy link
Member

glwagner commented Jun 5, 2025

yes that is fine

@amontoison
Copy link
Collaborator

It is important to ensure that the preconditioner is strictly SPD numerically.
A very small shift is enough but it is just to protect against a null-curvature direction.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
immersed boundaries ⛰️ Less Ocean, more anigans numerics 🧮 So things don't blow up and boil the lobsters alive
Projects
None yet
Milestone
No milestone
3 participants
@glwagner @amontoison @xkykai