Numerical Instability in Conjugate Gradient Poisson Solver with Immersed Boundary and Non-uniform Grids

[liuchihl]

In the MWE below, I test the ConjugateGradientPoissonSolver with an immersed boundary and non-uniform grids. The model converges successfully when using uniform grids. However, with non-uniform grids (as shown in the MWE), it quickly becomes unstable and blows up. Interestingly, certain non-uniform grid configurations do yield convergent results (e.g., when z_faces_array = [0., 100., 200., 600.] in this particular MWE), but I have yet to pinpoint the exact reason. I have tried addressing this issue by incorporating a preconditioner and regularizer, but these attempts have not resolved the problem.

When switching to an FFT-based Poisson solver, the results converge without any problems.

MWE:

using Printf
using Oceananigans
using Oceananigans.Solvers: ConjugateGradientPoissonSolver, fft_poisson_solver
using Oceananigans.Utils: prettytime
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary

## Simulation parameters
Nx = 4
Ny = 4

tᶠ = 5days # simulation run time
Δtᵒ = 180minutes # interval for saving output
H = 600meters # vertical extent
L = 2*2H # horizontal extent

# uniform grid works without blowing up
# z_faces_array = [0., 200., 400., 600.]
# nonuniform grid blows up
z_faces_array = [0., 50., 200., 600.]

grid = RectilinearGrid(size=(Nx, Ny, length(z_faces_array)-1), 
        x = (0, L),
        y = (0, L), 
        z = z_faces_array,
        halo = (4,4,4),
        topology = (Periodic, Periodic, Bounded)
)
h = 200meters # topographic height
topography = zeros(Nx,Ny)
topography[Nx÷2,Ny÷2] = h
grid_immerse = ImmersedBoundaryGrid(grid, GridFittedBottom(topography))

# IC such that flow is in phase with predicted linear response, but otherwise quiescent
U₀=0.025
uᵢ(x, y, z) = -U₀
vᵢ(x, y, z) = 0.

model = NonhydrostaticModel(
    grid = grid_immerse,
    pressure_solver = ConjugateGradientPoissonSolver(grid_immerse))

set!(model, u=uᵢ, v=vᵢ)

Δt = 30
simulation = Simulation(model, Δt = Δt, stop_time = tᶠ)

## Progress messages
progress_message(s) = @info @sprintf("[%.2f%%], iteration: %d, time: %.3f, max|u|: %.2e, max|w|: %.2e",
                            100 * s.model.clock.time / s.stop_time, s.model.clock.iteration,
                            s.model.clock.time, maximum(abs, model.velocities.u),maximum(abs, model.velocities.w)) 
simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(Δt))
run!(simulation)

Progress message:

[ Info: Initializing simulation...
[ Info: [0.00%], iteration: 0, time: 0.000, max|u|: 9.69e-01, max|w|: 1.78e+00
[ Info:     ... simulation initialization complete (10.942 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (33.779 ms).
[ Info: [0.01%], iteration: 1, time: 30.000, max|u|: 2.04e+05, max|w|: 2.26e+05
[ Info: [0.01%], iteration: 2, time: 60.000, max|u|: 1.83e+33, max|w|: 2.16e+32
[ Info: [0.02%], iteration: 3, time: 90.000, max|u|: 1.60e+252, max|w|: 4.04e+251
[ Info: [0.03%], iteration: 4, time: 120.000, max|u|: NaN, max|w|: NaN

[glwagner]

@liuchihl what happens if you use

pressure_solver = ConjugateGradientPoissonSolver(grid_immerse, maxiter=10)

?

Another thing to look at is #3848 to see if that helps... if it does, we may have finally found motivation to merge that PR.

[liuchihl]

@glwagner
I found that when maxiter=10, the problem persists, but when maxiter is less than 3, the result actually converges.
pressure_solver = ConjugateGradientPoissonSolver(grid_immerse; maxiter=2)

Another thing to look at is #3848 to see if that helps... if it does, we may have finally found motivation to merge that PR.

As for adding preconditioner, the result immediately goes to NaN, and I've tested different values of regularizer, but none of them helps at least in this MWE.

[glwagner]

and I've tested different values of regularizer

What have you tested?

Convergence for smaller maxiter is consistent with a non-convergent solver I think. I would check a few other things. For example preconditioner=DiagnoallyDominantPreconditioner() is probably worth a try.

Also can you check that the preconditioner is indeed FourierTridiagonalPoissonSolver? With the vertically stretched grid the solver has to include the tridiagonal part.

@xkykai may want to see this. It's useful that a vertically stretched grid illustrates the problem.

[liuchihl]

What have you tested?

I have run a loop with regularizer = -1:0.1:10, and these do not have a converging result.

Convergence for smaller maxiter is consistent with a non-convergent solver I think. I would check a few other things. For example preconditioner=DiagnoallyDominantPreconditioner() is probably worth a try.

Could you point me to an example of how I can switch to DiagonallyDominantPreconditioner as the preconditioner? I have trouble making it to work.

Also can you check that the preconditioner is indeed FourierTridiagonalPoissonSolver? With the vertically stretched grid the solver has to include the tridiagonal part.

I think preconditioner is indeed FourierTridiagonalPoissonSolver:

julia> preconditioner = fft_poisson_solver(grid)
FourierTridiagonalPoissonSolver{RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.T

[glwagner]

I think preconditioner is indeed FourierTridiagonalPoissonSolver:

Check your pressure_solver directly to be sure of it.

I have trouble making it to work.

You are the best person to construct this example since I think you are using the conjugate gradient solver the most. If you can't build it, perhaps post the error here? We can also add a PR with a test.

[liuchihl]

Indeed, preconditioner is FourierTridiagonalPoissonSolver

julia> model.pressure_solver
ConjugateGradientPoissonSolver:
├── grid: 4×4×3 ImmersedBoundaryGrid{Float64, Periodic, Periodic, Bounded} on CPU with 4×4×4 halo
└── conjugate_gradient_solver: ConjugateGradientSolver
    ├── maxiter: 1
    ├── reltol: 1.49012e-8
    ├── abstol: 1.49012e-8
    ├── preconditioner: Oceananigans.Solvers.RegularizedPoissonPreconditioner{FourierTridiagonalPoissonSolver{RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, Float64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetV

I realized that the name of preconditioner had been changed to AsymptoticPoissonPreconditioner in this branch, which I've merged to my branch. That's why I have trouble calling DiagonallyDominantPreconditioner.

When using AsymptoticPoissonPreconditioner, the result still does not converge when the grids are stretched with maxiter>2.

[glwagner]

@jm-c see this issue where someone's having trouble with the MITgcm preconditioner AsymptoticPoissonPreconditioner on a stretched grid. Do you have any ideas what might be causing the issue?

[liuchihl]

As previously noted, the numerical instability appears to be sensitive to how the vertical grids are stretched when using the CG Poisson solver.

The goal here is to compare the default FFT Poisson solver with the CG Poisson solver under conditions where the results converge. The simulation setup is straightforward, featuring an immersed boundary at the center with a constant velocity directed to the left. The horizontal boundaries are periodic.

Using the FFT Poisson solver (first video), nonzero velocity divergence is observed near the immersed boundary, motivating me to explore the CG solver as an alternative.

test_CGsolver_with_immersed_nonuniformgrid_FFT.mp4

However, when using CG Poisson solver, I see nonzero velocity divergent everywhere in the domain. Not sure what is wrong, but I am expecting to see divergence-free when using CG.

test_CGsolver_with_immersed_nonuniformgrid_CG.mp4

Here is the code of this simulation:

using Printf
using Oceananigans
using Oceananigans.Solvers: ConjugateGradientPoissonSolver, 
                    fft_poisson_solver, FourierTridiagonalPoissonSolver, AsymptoticPoissonPreconditioner
using Oceananigans.Utils: prettytime
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary
using Oceananigans.Solvers: FFTBasedPoissonSolver
using Oceananigans.Operators: divᶜᶜᶜ

## Simulation parameters
Nx = 32
Ny = 32

tᶠ = 5days # simulation run time
Δtᵒ = 30#180minutes # interval for saving output
H = 600meters # vertical extent
L = 2*2H # horizontal extent

# uniform grid works without blowing up
# z_faces_array = [0., 200., 400., 600.]
# nonuniform grid blows up
# z_faces_array = [range(0,300,13); range(330,600,6)]
z_faces_array = [0., 20., 50., 90., 140., 200., 270., 350., 440., 540., 600.]

grid = RectilinearGrid(size=(Nx, Ny, length(z_faces_array)-1), 
        x = (0, L),
        y = (0, L), 
        z = z_faces_array,
        halo = (4,4,4),
        topology = (Periodic, Periodic, Bounded)
)
h = 200meters # topographic height
topography = zeros(Nx,Ny)
topography[Nx÷2,Ny÷2] = h
grid_immerse = ImmersedBoundaryGrid(grid, GridFittedBottom(topography))

# IC such that flow is in phase with predicted linear response, but otherwise quiescent
U₀=0.01
uᵢ(x, y, z) = -U₀
vᵢ(x, y, z) = 0.
# preconditioner = fft_poisson_solver(grid)
preconditioner = AsymptoticPoissonPreconditioner()
regularizer = 1
model = NonhydrostaticModel(
    grid = grid_immerse,
    # comment out the pressure_solver when using deafult FFT solver
    pressure_solver = ConjugateGradientPoissonSolver(grid_immerse; maxiter=100,preconditioner)
    
)
set!(model, u=uᵢ, v=vᵢ)

Δt = 30
simulation = Simulation(model, Δt = Δt, stop_time = tᶠ)

u, v, w = model.velocities

udiv = KernelFunctionOperation{Center, Center, Center}(divᶜᶜᶜ, model.grid, u, v, w)
fname = "test_CGsolver_with_immersed_nonuniformgrid_CG.nc"
simulation.output_writers[:fields] = NetCDFOutputWriter(model, (udiv=udiv,u=u,v=v,w=w),
                                        schedule = TimeInterval(Δtᵒ),
                                        filename = fname,
                                        overwrite_existing = true)
## Progress messages
progress_message(s) = @info @sprintf("[%.2f%%], iteration: %d, time: %.3f, max|u|: %.2e, max|w|: %.2e",
                            100 * s.model.clock.time / s.stop_time, s.model.clock.iteration,
                            s.model.clock.time, maximum(abs, model.velocities.u),maximum(abs, model.velocities.w)) 
simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(Δt))
run!(simulation)



## plotting animation to see where blows up

using CairoMakie
using NCDatasets
using Printf
fname = "test_CGsolver_with_immersed_nonuniformgrid_FFT.nc"

ds = Dataset(fname,"r")

# grids
zC = ds["zC"]; Nz=length(zC)
zF = ds["zF"]; #Nz=length(zF)
xC = ds["xC"]; Nx=length(xC)
xF = ds["xF"];

yC = ds["yC"]; Ny=length(yC)
t = ds["time"];

u = ds["u"][:,:,:,:];
w = ds["w"][:,:,:,:];
udiv = ds["udiv"][:,:,:,:];

u_center = (u[:,:,:,:].+vcat(u[2:end,:,:,:], u[1:1,:,:,:]))./2
w_center = (w[:,:,1:end-1,:].+w[:,:,2:end,:])./2
u_center[u_center.==0].=NaN
w_center[w_center.==0].=NaN
w[w.==0].=NaN
u[u.==0].=NaN


# plot
n = Observable(1)
uₙ = @lift(u_center[:,16,:,$n])
wₙ = @lift(w_center[:,16,:,$n])
udivₙ = @lift(udiv[:,16,:,$n])

fig = Figure(resolution = (1000, 1000), figure_padding=(10, 40, 10, 10), size=(600,800),fontsize=20)
axis_kwargs = (xlabel = "x (m)",
                  ylabel = "z (m)",
                  limits = ((0, ds["xF"][end]), (0, ds["zF"][end])),
                  )
title = @lift @sprintf("t=%1.2f hrs", t[$n]/3600)
fig[1, :] = Label(fig, title, fontsize=20, tellwidth=false)
                  
                  
ax_u = Axis(fig[2, 1]; title = "u", axis_kwargs...)
ax_w = Axis(fig[3, 1]; title = "w", axis_kwargs...)
ax_udiv = Axis(fig[4, 1]; title = L"∇⋅\vec{u}", axis_kwargs...)



using ColorSchemes
U₀ = 0.01
hm_u = heatmap!(ax_u, xC[:], zC[:], uₙ,
    colorrange = (-U₀, U₀), colormap = :diverging_bwr_20_95_c54_n256,
    lowclip=cgrad(:diverging_bwr_20_95_c54_n256)[1], highclip=cgrad(:diverging_bwr_20_95_c54_n256)[end],
    nan_color = :gray)
hm_w = heatmap!(ax_w, xC[:], zC[:], wₙ,
    colorrange = (-U₀, U₀), colormap = :diverging_bwr_20_95_c54_n256,
    lowclip=cgrad(:diverging_bwr_20_95_c54_n256)[1], highclip=cgrad(:diverging_bwr_20_95_c54_n256)[end],
    nan_color = :gray)
Colorbar(fig[3,2], hm_w; label = "m/s")

hm_udiv = heatmap!(ax_udiv, xC[:], zC[:], udivₙ,
    colorrange = (-1e-8,1e-8), colormap = :diverging_bwr_20_95_c54_n256,
    lowclip=cgrad(:diverging_bwr_20_95_c54_n256)[1], highclip=cgrad(:diverging_bwr_20_95_c54_n256)[end],
    nan_color = :gray)
Colorbar(fig[4,2], hm_udiv; label = "1/s")


frames =  (1:50:length(t))

filename = join(split(fname, ".")[1:end-1], ".")

record(fig, string(filename,".mp4"), frames, framerate=23) do i
    @info "Plotting frame $i of $(frames[end])..."
    n[] = i
end

[glwagner]

However, when using CG Poisson solver, I see nonzero velocity divergent everywhere in the domain.

It doesn't look wrong to me. There is a very small non-zero divergence. How non-divergent the flow is will depend on the tolerance you use. To find expected behavior, show that the divergence decreases when you reduce reltol or abstol for the solver. This is the stopping criteria that the solver uses to figure out when to stop iterating.

[liuchihl]

That's a good point, I'll test that!

[liuchihl]

Thanks, you are exactly right about that! I do see a decrease in divergence when reducing tolerance.