Pressure solver algorithm involves two superfluous fields

[glwagner]

It appears that we allocate memory for a few more 3D arrays than needed in order to calculate pressure in NonhydrostaticModel. We currently allocate 3 arrays: one for hydrostatic pressure, one for non-hydrostatic pressure, and a third array with Complex{Float64} elements that's used for in-place FFTs while solving the pressure Poisson equation.

For our algorithm to be correct, however, we only need one array with Complex{Float64} elements. This array can than be used as scratch space when calculating the hydrostatic pressure and computing hydrostatic pressure gradients. It will then be overwritten when calculating the nonhydrostatic pressure component.

The distinction between the solver field with Complex{Float64} elements and the nonhydrostatic pressure is just that the nonhydrostatic pressure is real, which we enforce when copying it over:

Oceananigans.jl/src/Solvers/fft_based_poisson_solver.jl

Lines 102 to 105 in 4f80c38

	@kernel function copy_real_component!(ϕ, ϕc)
	i, j, k = @index(Global, NTuple)
	@inbounds ϕ[i, j, k] = real(ϕc[i, j, k])
	end

However, we could just as easily extract the real component when computing the pressure gradient:

Oceananigans.jl/src/Models/NonhydrostaticModels/pressure_correction.jl

Lines 34 to 40 in 4f80c38

	@kernel function _pressure_correct_velocities!(U, grid, Δt, pNHS)
	i, j, k = @index(Global, NTuple)
	@inbounds U.u[i, j, k] -= ∂xᶠᶜᵃ(i, j, k, grid, pNHS) * Δt
	@inbounds U.v[i, j, k] -= ∂yᶜᶠᵃ(i, j, k, grid, pNHS) * Δt
	@inbounds U.w[i, j, k] -= ∂zᵃᵃᶠ(i, j, k, grid, pNHS) * Δt
	end

[glwagner]

cc @tomchor

[tomchor]

Nice catch. Does this only happen for the fft-based solver?

[glwagner]

All of the pressure solvers require storage arrays whose memory space can, in principle, be reused to hold either the hydrostatic or non-hydrostatic pressure.

The main detriment is that one of the pressure fields has to be overwritten. In our algorithm I think that means it would no longer be possible to output the nonhydrostatic pressure component.

[glwagner]

I'm closing this issue because I'm judging that it's not of current, timely relevance to Oceananigans development. If you would like to make it a higher priority or if you think the issue was closed in error please feel free to re-open.