Laplacian operator for `ConjugateGradientPoissonSolver` is nonsymmetric for nonuniform grids

[xkykai]

As seen in this discussion below, the Laplacian operator used is not symmetric. This is because

Oceananigans.jl/src/Solvers/conjugate_gradient_poisson_solver.jl

Line 35 in 9c3c4b8

@inbounds ∇²ϕ[i, j, k] = ∇²ᶜᶜᶜ(i, j, k, grid, ϕ)

uses
https://github.com/CliMA/Oceananigans.jl/blob/9c3c4b8b618b9db17e2bdcbe887f0d55bcfe5e54/src/Operators/laplacian_operators.jl#L36C1-L40C4
which contains the reciprocal volume operator V⁻¹ᶜᶜᶜ(i, j, k, grid) on the LHS, which makes the Laplacian operator nonsymmetric.

This is an issue as strictly speaking a conjugate gradient solver should only be used on symmetric positive semi-definite matrices, thus it will go unstable.

I suggest to build the matrix from the operator and checking that, indeed, the matrix is nonsymmetric in your case.

Indeed, here's a MWE to test compute_laplacian! which is found here

Oceananigans.jl/src/Solvers/conjugate_gradient_poisson_solver.jl

Line 38 in 6383a85

function compute_laplacian!(∇²ϕ, ϕ)

using Oceananigans
using Oceananigans.Solvers: compute_laplacian!
using Oceananigans.BoundaryConditions: fill_halo_regions!
using SparseArrays

xs = collect(range(0, 1, length=6))
xs[2:end-1] .+= rand(length(xs[2:end-1])) * (1 / 6) / 5
@info xs

grid = RectilinearGrid(CPU(), Float64,
                        size = (5), 
                        halo = 4,
                        x = xs,
                        topology = (Bounded, Flat, Flat))

function initialize_matrix(template_field, linear_operator!, args...)
    Nx, Ny, Nz = size(template_field)
    A = spzeros(eltype(template_field.grid), Nx*Ny*Nz, Nx*Ny*Nz)
    make_column(f) = reshape(interior(f), Nx*Ny*Nz)

    eᵢⱼₖ = similar(template_field)
    ∇²eᵢⱼₖ = similar(template_field)

    for k = 1:Nz, j in 1:Ny, i in 1:Nx
        parent(eᵢⱼₖ) .= 0
        parent(∇²eᵢⱼₖ) .= 0
        eᵢⱼₖ[i, j, k] = 1
        fill_halo_regions!(eᵢⱼₖ)
        linear_operator!(∇²eᵢⱼₖ, eᵢⱼₖ, args...)

        A[:, Ny*Nx*(k-1) + Nx*(j-1) + i] .= make_column(∇²eᵢⱼₖ)
    end
    
    return A
end

test_field = CenterField(grid)
A = initialize_matrix(test_field, compute_laplacian!)
display(A)

and the output is

[ Info: [0.0, 0.22960116301095507, 0.4189282613194224, 0.6149861209945792, 0.804943345574523, 1.0]

5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries:
 -20.793    20.793      ⋅         ⋅         ⋅
  25.2161  -52.6269   27.4108     ⋅         ⋅
    ⋅       26.4698  -52.8964   26.4266     ⋅
    ⋅         ⋅       27.2753  -54.6216   27.3463
    ⋅         ⋅         ⋅       26.6313  -26.6313

This shows that compute_laplacian! and hence

Oceananigans.jl/src/Operators/laplacian_operators.jl

Line 36 in 6383a85

@inline function ∇²ᶜᶜᶜ(i, j, k, grid, c)

is nonsymmetric indeed.
Should we fix this? How to fix it? In my mind it seemed to make sense that the Laplacian operator could be nonsymmetric on nonuniform grids, but perhaps I am missing something...

Originally posted by @xkykai in #4563 (comment)