Can we rewrite the pressure equation and calculate each term in it?

[zwei961120]

Hello everyone,

I'm trying to calculate the pressure terms shown in the following figure (Eq.(2.5) from Pearson et al., 2019):

The equation in the above figure I think is another form of writing the pressure equation to maintain incompressibility, but different from the common way ($p = p_{HY'}+ p_{NHS}$), it may provide new insights into how different processes contribute to pressures and help to improve second-moment vertical mixing parameterizations, which is the motivation that I want to calculate them.

What I'm trying to do is write the code that repeats the pressure solver algorithm for non-hydrostatic pressure (cause this is also an elliptic problem I think) to calculate the terms in the above equation, and re-run the pressure solver separately for each of these new terms. Lastly, I want to save these new pressure terms in a way that will not overwrite the original pressures.

I tried to start with the second term ($p_r'$) in the RHS of the equation. This term is the component that results from rapid distortion by the current shear. For this term, the vertical boundary condition is $\partial p_r' / \partial x_3 = 0$ (see more details about boundary conditions in Appendix A in Pearson et al., 2019). So the goal is to solve an equation $\nabla^2p_r' = -2\frac{\partial U_i}{\partial x_j}\frac{\partial u_j'}{\partial x_i}$ with applying the boundary condition $\partial p_r' / \partial x_3 = 0$.

I have written a part of the code (I start from the calculation of $p_r'$ because there is no need to consider the boundary conditions, and the calculation I think can directly follow the calculation of $p_{NHS}$). Now I have tested it on the CPU, and output the pressure (I call it “prapid” in my code) like model.pressures.prapid , but there are also some problems. I add part of the code in https://github.com/CliMA/Oceananigans.jl/blob/main/src/Models/NonhydrostaticModels/solve_for_pressure.jl following the cacluation of $p_{NHS}$ like:

@kernel function calculate_pressure_source_term_fft_based_solver_rapid!(rhs, grid, Δt, U★)
   i, j, k = @index(Global, NTuple)
   U = Field(Average(U★.u, dims=(1, 2)))
   V = Field(Average(U★.v, dims=(1, 2)))
   W = Field(Average(U★.w, dims=(1, 2)))
   u_prime = Field(U★.u - U)
   v_prime = Field(U★.v - V)
   w_prime = Field(U★.w - W)

   @inbounds rhs[i, j, k] = -2*(  ∂xᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, w_prime) \
                                + ∂xᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, w_prime) \
                                + ∂xᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, w_prime))
end

@kernel function calculate_permuted_pressure_source_term_fft_based_solver_rapid!(rhs, grid, Δt, U★, ::ZXYPermutation)
   i, j, k = @index(Global, NTuple)
   U = Field(Average(U★.u, dims=(1, 2)))
   V = Field(Average(U★.v, dims=(1, 2)))
   W = Field(Average(U★.w, dims=(1, 2)))
   u_prime = Field(U★.u - U)
   v_prime = Field(U★.v - V)
   w_prime = Field(U★.w - W)

   @inbounds rhs[k, i, j] = -2*(  ∂xᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, w_prime) \
                                + ∂xᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, w_prime) \
                                + ∂xᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, w_prime))
end


@kernel function calculate_permuted_pressure_source_term_fft_based_solver_rapid!(rhs, grid, Δt, U★, ::ZYXPermutation)
   i, j, k = @index(Global, NTuple)
   U = Field(Average(U★.u, dims=(1, 2)))
   V = Field(Average(U★.v, dims=(1, 2)))
   W = Field(Average(U★.w, dims=(1, 2)))
   u_prime = Field(U★.u - U)
   v_prime = Field(U★.v - V)
   w_prime = Field(U★.w - W)

   @inbounds rhs[k, j, i] = -2*(  ∂xᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, w_prime) \
                                + ∂xᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, w_prime) \
                                + ∂xᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, W)*∂zᶜᶜᶜ(i, j, k, grid, w_prime))
end 


function solve_for_pressure_rapid!(pressure, solver::DistributedFFTBasedPoissonSolver, Δt, U★)
   rhs = parent(first(solver.storage))
   arch = architecture(solver)
   grid = solver.local_grid

   launch!(arch, grid, :xyz, calculate_permuted_pressure_source_term_fft_based_solver_rapid!,
           rhs, grid, Δt, U★, solver.input_permutation,)

   # Solve pressure Poisson equation for pressure, given rhs
   solve!(pressure, solver)

   return pressure
end

I tried to run it on CPU (Oceananigans v0.84.1 & Julia v1.8.4) like (A simple case of Langmuir turbulence case in a shallow water):

using Oceananigans
using Printf
using Oceananigans.Units

# Grid
const H = 30
grid = RectilinearGrid(GPU();
                      size = (100, 100, 48),
                      halo = (3, 3, 3),
                         x = (0, 2π*H),
                         y = (0, 2π*H),
                         z = (-H, 0))
# Buoyancy
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-4, haline_contraction=8e-4))

# Surface boundary conditions
u₁₀ = 10    
cᴰ = 1.25e-3 
ρₐ = 1.225  
ρₒ = 1026.0 
Qᵘ = - ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀) 


# Bottom boundary conditions
z₀ = 0.01
κ = 0.4 
z₁ = grid.Δzᵃᵃᶠ/2 
cᴰ = (κ / log(z₁ / z₀))^2 

@inline drag_u(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + v^2) * u
@inline drag_v(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + v^2) * v

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(; cᴰ))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(; cᴰ))

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)


# Coriolis
f₀ = 1e-4
coriolis = FPlane(f=f₀)

# Background current                                                                                                                                                                             
U∞ = 0.25 # m s⁻¹

# Stokes drift
using Oceananigans.BuoyancyModels: g_Earth

const amplitude = 0.8 
const wavelength = 60  
const wavenumber = 2π / wavelength 
const frequency = sqrt(g_Earth * wavenumber * tanh(wavenumber * H)) 

## Stokes drift velocity at the surface
const Uˢ = amplitude^2 * wavenumber * frequency # m s⁻¹ # 0.0676  

# The vertical derivative of the Stokes drift is
∂z_uˢ(z, t) = Uˢ * wavenumber * sinh(2 * wavenumber * (z + H)) / (sinh(wavenumber * H))^2

# ## Forcing
#
# pressure gradient force due Coriolis in geostrophic balance
#

const pgrad = f₀ * U∞
@inline pressure_gradient(x, y, z, t) = pgrad
v_forcing = Forcing(pressure_gradient)

# ## Model instantiation
model = NonhydrostaticModel(; grid,
               advection = WENO(grid; order=5),
             timestepper = :RungeKutta3,
                 tracers = (:T, :S),
                 coriolis = coriolis,
                buoyancy = buoyancy,
                 closure = SmagorinskyLilly(),
              stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
             boundary_conditions = (;u = u_bcs, v = v_bcs),
                 forcing = (;v=v_forcing))
S₀ = 35.0
T₀ = 20.0
dTdz = 0.1 # K m⁻¹
T0(z) = T₀ + dTdz * z
Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # noise

## Horizontal velocity initial condition: geostrophic flow with random noise superposed
uᵢ(x, y, z) = U∞ + 1e-5 * Ξ(z)

## Vertical velocity initial condition: random noise scaled by the friction velocity.
wᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-3 * Ξ(z)

Tᵢ(x, y, z) = T0(z) + dTdz * model.grid.Lz * 1e-6 * Ξ(z)

## `set!` the `model` fields using functions or constants:
set!(model, u=uᵢ, w=wᵢ, T=Tᵢ, S=S₀)

When I tried to output model.pressures.prapid, I got:

100×100×48 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 100×100×48 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 106×106×54 OffsetArray(::Array{Float64, 3}, -2:103, -2:103, -2:51) with eltype Float64 with indices -2:103×-2:103×-2:51
   └── max=NaN, min=NaN, mean=NaN

The NaN value appears. But for $p_{NHS}$ and $p_{HY'}$, the results are reasonable. So any comments or ideas about this error?

Another issue is about the calculation of the pressure terms for which the boundary conditions are inhomogeneous Neumann boundary conditions ($p_b'$, $p_{st}'$ and $p_{sg}'$, see Appendix A for details in Pearson et al., 2019). For this issue, I have checked the related issue #1232. I wanna ask do I need to subtly update the boundary values of $p_b'$, $p_{st}'$ and $p_{sg}'$ to convert the inhomogeneous Neumann problem to the homogeneous Neumann problem? If so, how do I determine the values the boundary should update?

I hope I've described this issue clearly, but please let me know if anything confuse you. In addition, this is my first attempt to make a possible contribution to an open source project, and I would like to ask what is the most efficient way to conduct the discussion? Should we talk about this here firstly or can I start with a PR to better let you know what specific changes I've made to the code?

Thank you very much for any comments, suggestions or ideas!

[glwagner]

This is very interesting, thank you for starting this! Do you mind if I convert this to a discussion and we can continue the conversation over there? We're trying to use discussions for science discussions and illustrations of how to use the code, whereas issues are for bugs and problems that require source code changes.

[glwagner]

Something I see that needs to be fixed right away is that the "mean" and "prime" terms have to be precomputed outside a kernel. For example,

@kernel function calculate_pressure_source_term_fft_based_solver_rapid!(rhs, grid, Δt, U★)
   i, j, k = @index(Global, NTuple)
   U = Field(Average(U★.u, dims=(1, 2)))
   V = Field(Average(U★.v, dims=(1, 2)))
   W = Field(Average(U★.w, dims=(1, 2)))

probably won't work. The other issue is that we have to call compute! to precompute these fields, which is not being done.

Also, just as a matter of approach, I would suggest not changing the source code to do this problem but instead using the solvers directly to compute the solution to this new Poisson problem. That will allow you to more quickly compare the pressure computed by the native algorithm to the "new" pressure that is computed offline.

The first step towards solving this problem is to show that you can solve Poisson's equation correctly using FFTBasedPoissonSolver! (This is the solver that will work for grids that are regularly spaced in all directions.) Once you have code that can solve the Poisson equation for a simple RHS, then I suggest building up the RHS terms that you have piece by piece. When you add a new term to the RHS, plot the term to make sure it looks right. Then solve Poisson's equation with the new term, and make sure the solution looks good. Keep adding terms ensuring that each time you build up the complexity the result is as expected and you don't get NaNs.

[glwagner]

To be specific, here are the steps that have to be taken:

1. Define the mean velocities and perturbations:

u, v, w = model.velocities

# These Fields will have to be computed every time step
U = Field(Average(u, dims=(1, 2)))
V = Field(Average(v, dims=(1, 2)))

# These are operations, not Fields
u′ = u - U
v′ = v - V

Note that it is best to leave the perturbations as AbstractOperations rather than wrapping them in Fields. This will simplify your computation and also save computational expense. Furthermore, it's unlikely you have a mean vertical velocity. You don't need to compute the horizontally-averaged vertical velocity. That also simplifies things.

2. Every time-step, you must compute the means.

Thus somewhere in your code, before you calculate the RHS of the Poisson equation, you must call

compute!(U)
compute!(V)

This will compute the horizontal averages. Because horizontal averages are reductions, they cannot be computed "on the fly". They have to be computed in advance. However, the perturbation velocities can be computed on the fly, so it is best to leave them as abstract operations.

Finally, a difficulty I can forsee is that the equations you have written are not all located at cell centers on the staggered grid. It's unclear to me the best way to "average" the terms so that they end up at cell centers. Either way, you will have to take care of the staggered grid to make these computations I think. Alternatively, you can try to express everything as an AbstractOperation. In that case, locating the terms will be done automatically (though not necessarily in a way that is consistent with a discrete derivation of this pressure Poisson equation...)

[BrodiePearson]

Thanks for putting your thoughts about this together Zheng! I support Greg's suggestion of moving this to the discussion section.

In response to your BC question: since I don't yet understand the nuances of the Oceananigans pressure solver and time stepping, I would follow Greg's guidance with the aim of writing code that calculates all of the individual pressure components using the existing /default Oceananigans homogeneous Neumann BC implementation for each term. Once you have confidence that your code works and is using the correct RHS for each pressure component, I would check the residual (how much of the total non-hydrostatic pressure is explained by the sum of the pressure components). If the residual is very small I do not think it would be necessary to delve into the BC details. If the residual is large, you probably need to look at the BC implementation in more detail.

[glwagner]

The way the pressure boundary conditions enter is subtle. To solve the Poisson equation we use an eigenfunction expansion of the discrete second-order Poisson operator (combined with a vertical tridiagonal solve, if the grid is vertically-stretched). In this algorithm, any contribution from the boundary conditions must be added to the right hand side. Otherwise, the eigenfunction boundary conditions apply, which are homogeneous Neumann for Bounded directions and periodic for Periodic directions.

In a fractional step algorithm, the "extra" RHS terms that enforce boundary conditions are added in a subtle way by enforcing no-penetration on the predictor velocity field. The predictor velocity field will exhibit flow through the boundary in the event that a boundary pressure gradient exists. By eliminating this penetration velocity, the proper terms corresponding to a non-zero pressure gradient boundary condition are inserted into the Poisson equation RHS.

It should be possible to compute the boundary contributions empirically (by computing tendencies on the boundaries) but this will take a little bit of work.

[BrodiePearson]

Thanks @glwagner!

Another question: Is the velocity being used in @zwei961120's above examples representing the Lagrangian momentum or the Eulerian momentum?

If it is the Lagrangian momentum the Stokes drift will need to be subtracted out before calculating the rapid term, otherwise you will actually be calculating the sum of the rapid term and part of the Stokes term (specifically the part that is not associated with the curl calculated in StokesDrifts.jl).

[glwagner]

Oceananigans prognoses the Lagrangian-mean velocity. So model.velocities are Lagrangian-mean. More information here: https://clima.github.io/OceananigansDocumentation/stable/physics/surface_gravity_waves/

[glwagner]

Note, the "effective rotation term" that StokesDrifts.jl implements appears only in the Lagrangian-mean formulation of the wave-averaged Boussinesq equations.

[glwagner]

To get started with the pressure solver, there is a script here that illustrates how to use it:

https://github.com/glwagner/Krylovable.jl/blob/main/fft_versus_pcg_solver.jl

here's the gist:

using Oceananigans
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Solvers: FFTBasedPoissonSolver, solve!

using Statistics
using Random
Random.seed!(123)

#####
##### Grid setup
#####

arch = CPU() # try changing to GPU()
Nx = Ny = Nz = 64
topology = (Periodic, Periodic, Periodic)

# Regularly-spaced grid for now:
x = y = z = (0, 2π)

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

#####
##### Right-hand side generation: sum of randomly-place Gaussians
#####

δ = 0.1 # Gaussian width
gaussian(ξ, η, ζ) = exp(-(ξ^2 + η^2 + ζ^2) / 2δ^2)
Ngaussians = 17
Ξs = [2π * rand(3) for _ = 1:Ngaussians]

function many_gaussians(ξ, η, ζ)
    val = zero(ξ)
    for Ξ₀ in Ξs
        ξ₀, η₀, ζ₀ = Ξ₀
        val += gaussian(ξ - ξ₀, η - η₀, ζ - ζ₀)
    end
    return val
end

#####
##### FFT-based Poisson solver
#####

# Solve lap(p) = r
fft_solver = FFTBasedPoissonSolver(grid)

# Note: we use in-place transforms, so the RHS has to be AbstractArray{Complex{T}}.
# So, we first fill up "r" and then copy it into "rc = fft_solver.storage",
# which has the correct type.
r = CenterField(grid)
set!(r, many_gaussians)
parent(r) .-= mean(interior(r))
rc = fft_solver.storage 
rc .= interior(r)

pfft = CenterField(grid)
solve!(pfft, fft_solver, rc)

This will solve the Poisson equation with the right hand side r.

[zwei961120]

Thank you so much for very useful information and code, I will have a try :-)

[zwei961120]

I'm trying to calculate the RHS and the first thing I did is calculating the Stokes drift velocity profile like:

### Calculate the constant Stokes drift velocity profile (only applied in x-direction in my case)

uˢ(x, y, z) = Uˢ * cosh(2 * wavenumber * (z + H)) / (2*(sinh(wavenumber * H))^2)
us = CenterField(grid)
set!(us, uˢ)

And I got:

196×196×96 Field{Center, Center, Center} on RectilinearGrid on GPU
├── grid: 196×196×96 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on GPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 202×202×102 OffsetArray(::CUDA.CuArray{Float64, 3, CUDA.Mem.DeviceBuffer}, -2:199, -2:199, -2:99) with eltype Float64 with indices -2:199×-2:199×-2:99
    └── max=0.0658545, min=0.000254276, mean=0.0108277

And then I tried to write a function to calculate the RHS like:

### Calculate the Eulerian x-velocity by substracting Stokes drift in x-direction

u = model.velocities.u - us
v = model.velocities.v
w = model.velocities.w
U = Field(Average(u, dims=(1, 2)))
V = Field(Average(v, dims=(1, 2)))
compute!(U)
compute!(V)
u_prime =  u - U
v_prime =  v - V

### Define a function to calculate the RHS and try to make sure that each variable is located at cell center
function calculate_RHS_test(i, j, k, grid, U, V, u_prime, v_prime, w)
    return -2*(  ∂xᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, U)*∂xᶜᶜᶜ(i, j, k, grid, w) \
                 ∂xᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, u_prime) + ∂yᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, v_prime) + ∂zᶜᶜᶜ(i, j, k, grid, V)*∂yᶜᶜᶜ(i, j, k, grid, w))
end

RHS = KernelFunctionOperation{Center, Center, Center}(calculate_RHS_test, grid, U, V, u_prime, v_prime, w)
RHS_F = Field(RHS)
compute!(RHS_F)

Then I got:

196×196×96 Field{Center, Center, Center} on RectilinearGrid on GPU
├── grid: 196×196×96 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on GPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
├── operand: KernelFunctionOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 202×202×102 OffsetArray(::CUDA.CuArray{Float64, 3, CUDA.Mem.DeviceBuffer}, -2:199, -2:199, -2:99) with eltype Float64 with indices -2:199×-2:199×-2:99
    └── max=2.22636e-28, min=-2.22237e-28, mean=NaN

The value is very small and the NaN value appears. Any comments or suggestions about this? Thanks! (Julia version: 1.9.3, Oceananigans version: 0.90.6)

[glwagner]

Your code contains invalid operations. For example ∂xᶜᶜᶜ(i, j, k, grid, V), which attemps to compute the x derivative of V, ending up at cell centers. But V is at faces in y, and differentiating in x will not change it's location. To use those low-level operators, you need to know the precise location of every quantity on the staggered grid and utilize the convention embedded in three-letter superscripts (as in ᶜᶜᶜ) correctly. Another option is to write everything with abstract operations, which is guaranteed to be correct wrt to the staggered grid. That might work too.

[BrodiePearson]

@zwei961120: The easiest way I see to fix the problem @glwagner pointed out is to emulate the Oceananigans budget calculations, but only considering the budget terms relevant to a specific pressure component. For example:

1. Calculate the budget term of interest for each velocity component by using the same function as the tendency equation calculation. For the rapid term, the contribution from the u-budget would look like this (perhaps with the Stokes drift subtracted from the input velocities):

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl

Line 71 in 6b22a79

return ( - div_𝐯u(i, j, k, grid, advection, total_velocities, velocities.u)

2. You should get three budget terms (one for each velocity component) that are co-located with their respective velocity on the grid. You could store these in a data structure/collection(?) analogous to U★, for example U★_tendency.u co-located with U★.u, U★_tendency.v colocated with U★.v, U★_tendency.w colocated with U★.w

3. Feed these three velocity tendency components into the FFT calculation, which is already set up to calculate the divergence with this grid configuration, by replacing U★ $\rightarrow$ U★_tendency when calling this kernel

Oceananigans.jl/src/Models/NonhydrostaticModels/solve_for_pressure.jl

Lines 15 to 18 in 6b22a79

	@kernel function calculate_pressure_source_term_fft_based_solver!(rhs, grid, Δt, U★)
	i, j, k = @index(Global, NTuple)
	@inbounds rhs[i, j, k] = divᶜᶜᶜ(i, j, k, grid, U★.u, U★.v, U★.w) / Δt
	end

[zwei961120]

Thanks for your suggestions! I tried to figure out the precise location of each quantity and modified my script like:

uˢ(x, y, z) = Uˢ * cosh(2 * wavenumber * (z + H)) / (2*(sinh(wavenumber * H))^2)
us = CenterField(grid)
set!(us, uˢ)


u = model.velocities.u - us
v = model.velocities.v
w = model.velocities.w
U = Field(Average(u, dims=(1, 2)))
V = Field(Average(v, dims=(1, 2)))
compute!(U)
compute!(V)
u_prime =  u - U
v_prime =  v - V


@inline function calculate_RHS_test(i, j, k, grid, U, V, u_prime, v_prime, w)
    
    Ux = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, U) # F, F, F -> C, F, F -> C, C, C
    Uy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, U)
    Uz = ℑxyᶜᶜᵃ(i, j, k, grid, ∂zᶠᶠᶜ, U)

    Vx = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, V)
    Vy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, V)
    Vz = ℑxyᶜᶜᵃ(i, j, k, grid, ∂zᶠᶠᶜ, V)

    ux = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, u_prime)
    vx = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, v_prime)
    wx = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, w)

    uy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, u_prime)
    vy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, v_prime)
    wy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, w)

    return -2*(Ux*ux + Uy*vx + Uz*wx + Vx*uy + Vy*vy + Vz*wy) 
end

RHS = KernelFunctionOperation{Center, Center, Center}(calculate_RHS_test, grid, U, V, u_prime, v_prime, w)
RHS_F = Field(RHS)
compute!(RHS_F)

I output the maximum of the absolute value of RHS_F at each time step shown as:

[ Info: Initializing simulation...
Iteration: 0000, time: 0 seconds, Δt: 5 seconds, max(|RHSF|) = 1.50592e-07 ms⁻¹, wall time: 0 seconds, advective CFL: 1.30e+00, diffusive CFL: 0.00e+00
[ Info:     ... simulation initialization complete (17.109 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.633 seconds).
Iteration: 0100, time: 1.896 minutes, Δt: 329.868 ms, max(|RHSF|) = 1.62724e-08 ms⁻¹, wall time: 27.924 seconds, advective CFL: 1.00e-01, diffusive CFL: 5.66e-03
Iteration: 0200, time: 2.439 minutes, Δt: 321.483 ms, max(|RHSF|) = 1.74939e-08 ms⁻¹, wall time: 33.353 seconds, advective CFL: 1.00e-01, diffusive CFL: 5.37e-03
Iteration: 0300, time: 2.969 minutes, Δt: 314.523 ms, max(|RHSF|) = 1.96655e-08 ms⁻¹, wall time: 38.799 seconds, advective CFL: 1.00e-01, diffusive CFL: 5.11e-03
Iteration: 0400, time: 3.489 minutes, Δt: 308.528 ms, max(|RHSF|) = 2.23742e-08 ms⁻¹, wall time: 44.250 seconds, advective CFL: 1.00e-01, diffusive CFL: 4.89e-03
Iteration: 0500, time: 3.999 minutes, Δt: 303.231 ms, max(|RHSF|) = 2.82702e-08 ms⁻¹, wall time: 49.725 seconds, advective CFL: 1.00e-01, diffusive CFL: 4.69e-03
Iteration: 0600, time: 4.501 minutes, Δt: 298.479 ms, max(|RHSF|) = 3.25192e-08 ms⁻¹, wall time: 55.122 seconds, advective CFL: 1.00e-01, diffusive CFL: 4.51e-03
Iteration: 0700, time: 4.995 minutes, Δt: 294.170 ms, max(|RHSF|) = 3.89417e-08 ms⁻¹, wall time: 1.009 minutes, advective CFL: 1.00e-01, diffusive CFL: 4.36e-03
Iteration: 0800, time: 5.482 minutes, Δt: 290.227 ms, max(|RHSF|) = 4.78641e-08 ms⁻¹, wall time: 1.098 minutes, advective CFL: 1.00e-01, diffusive CFL: 4.22e-03
Iteration: 0900, time: 5.963 minutes, Δt: 286.590 ms, max(|RHSF|) = 5.78322e-08 ms⁻¹, wall time: 1.188 minutes, advective CFL: 1.00e-01, diffusive CFL: 4.09e-03
Iteration: 1000, time: 6.438 minutes, Δt: 283.212 ms, max(|RHSF|) = 6.86305e-08 ms⁻¹, wall time: 1.277 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.97e-03
Iteration: 1100, time: 6.908 minutes, Δt: 280.054 ms, max(|RHSF|) = 8.08776e-08 ms⁻¹, wall time: 1.367 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.87e-03
Iteration: 1200, time: 7.373 minutes, Δt: 277.090 ms, max(|RHSF|) = 1.02543e-07 ms⁻¹, wall time: 1.457 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.77e-03
Iteration: 1300, time: 7.832 minutes, Δt: 274.297 ms, max(|RHSF|) = 1.30010e-07 ms⁻¹, wall time: 1.548 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.68e-03
Iteration: 1400, time: 8.287 minutes, Δt: 271.656 ms, max(|RHSF|) = 1.68695e-07 ms⁻¹, wall time: 1.639 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.60e-03
Iteration: 1500, time: 8.738 minutes, Δt: 269.153 ms, max(|RHSF|) = 2.13026e-07 ms⁻¹, wall time: 1.730 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.52e-03
Iteration: 1600, time: 9.185 minutes, Δt: 266.772 ms, max(|RHSF|) = 2.70865e-07 ms⁻¹, wall time: 1.820 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.44e-03
Iteration: 1700, time: 9.628 minutes, Δt: 264.501 ms, max(|RHSF|) = 3.62140e-07 ms⁻¹, wall time: 1.910 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.37e-03
Iteration: 1800, time: 10.066 minutes, Δt: 262.337 ms, max(|RHSF|) = 4.76646e-07 ms⁻¹, wall time: 2.000 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.31e-03
Iteration: 1900, time: 10.501 minutes, Δt: 260.256 ms, max(|RHSF|) = 6.22043e-07 ms⁻¹, wall time: 2.090 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.25e-03
Iteration: 2000, time: 10.934 minutes, Δt: 258.257 ms, max(|RHSF|) = 8.34393e-07 ms⁻¹, wall time: 2.181 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.19e-03
Iteration: 2100, time: 11.363 minutes, Δt: 256.334 ms, max(|RHSF|) = 1.05704e-06 ms⁻¹, wall time: 2.271 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.14e-03
Iteration: 2200, time: 11.788 minutes, Δt: 254.481 ms, max(|RHSF|) = 1.39786e-06 ms⁻¹, wall time: 2.361 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.08e-03
Iteration: 2300, time: 12.211 minutes, Δt: 252.691 ms, max(|RHSF|) = 1.78541e-06 ms⁻¹, wall time: 2.450 minutes, advective CFL: 1.00e-01, diffusive CFL: 3.03e-03
Iteration: 2400, time: 12.631 minutes, Δt: 250.960 ms, max(|RHSF|) = 2.31677e-06 ms⁻¹, wall time: 2.541 minutes, advective CFL: 1.00e-01, diffusive CFL: 2.99e-03
Iteration: 2500, time: 13.048 minutes, Δt: 249.283 ms, max(|RHSF|) = 2.94794e-06 ms⁻¹, wall time: 2.631 minutes, advective CFL: 1.00e-01, diffusive CFL: 2.94e-03
...
Iteration: 11800, time: 43.349 minutes, Δt: 174.345 ms, max(|RHSF|) = 5.43150e-03 ms⁻¹, wall time: 10.970 minutes, advective CFL: 1.00e-01, diffusive CFL: 2.85e-03
Iteration: 11900, time: 43.639 minutes, Δt: 174.971 ms, max(|RHSF|) = 5.33748e-03 ms⁻¹, wall time: 11.059 minutes, advective CFL: 1.00e-01, diffusive CFL: 2.84e-03
Iteration: 12000, time: 43.931 minutes, Δt: 176.169 ms, max(|RHSF|) = 5.07422e-03 ms⁻¹, wall time: 11.148 minutes, advective CFL: 1.00e-01, diffusive CFL: 2.85e-03
Iteration: 12100, time: 44.226 minutes, Δt: 177.017 ms, max(|RHSF|) = 4.64305e-03 ms⁻¹, wall time: 11.237 minutes, advective CFL: 1.00e-01, diffusive CFL: 2.91e-03

The result I think is reasonable. Then I added the part to calculate the pressure term like:

fft_solver = FFTBasedPoissonSolver(grid)
r = CenterField(grid)
set!(r, RHS_F)
rc = fft_solver.storage 
rc .= interior(r)

pfft = CenterField(grid)
solve!(pfft, fft_solver, rc)

I tried to output the maximum of the absolute value of pfft at each time step shown as:

Iteration: 0100, time: 3.422 minutes, Δt: 1.546 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 27.178 seconds, advective CFL: 5.00e-01, diffusive CFL: 2.46e-02
Iteration: 0200, time: 5.908 minutes, Δt: 1.435 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 31.907 seconds, advective CFL: 5.00e-01, diffusive CFL: 2.05e-02
Iteration: 0300, time: 8.240 minutes, Δt: 1.360 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 36.803 seconds, advective CFL: 5.00e-01, diffusive CFL: 1.80e-02
Iteration: 0400, time: 10.459 minutes, Δt: 1.302 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 41.649 seconds, advective CFL: 5.00e-01, diffusive CFL: 1.63e-02
Iteration: 0500, time: 12.593 minutes, Δt: 1.256 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 46.461 seconds, advective CFL: 5.00e-01, diffusive CFL: 1.50e-02
Iteration: 0600, time: 14.655 minutes, Δt: 1.215 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 51.190 seconds, advective CFL: 5.00e-01, diffusive CFL: 1.45e-02
Iteration: 0700, time: 16.651 minutes, Δt: 1.178 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 55.906 seconds, advective CFL: 5.00e-01, diffusive CFL: 1.44e-02
Iteration: 0800, time: 18.584 minutes, Δt: 1.134 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.011 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.46e-02
Iteration: 0900, time: 20.433 minutes, Δt: 1.065 seconds, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.089 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.62e-02
Iteration: 1000, time: 22.160 minutes, Δt: 994.143 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.169 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.90e-02
Iteration: 1100, time: 23.765 minutes, Δt: 928.111 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.248 minutes, advective CFL: 5.00e-01, diffusive CFL: 2.08e-02
Iteration: 1200, time: 25.291 minutes, Δt: 918.704 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.327 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.90e-02
Iteration: 1300, time: 26.841 minutes, Δt: 944.256 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.407 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.73e-02
Iteration: 1400, time: 28.429 minutes, Δt: 955.870 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.490 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.57e-02
Iteration: 1500, time: 30.032 minutes, Δt: 961.020 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.577 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.49e-02
Iteration: 1600, time: 31.579 minutes, Δt: 921.443 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.664 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.43e-02
Iteration: 1700, time: 33.147 minutes, Δt: 948.858 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.750 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.42e-02
Iteration: 1800, time: 34.683 minutes, Δt: 907.696 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.837 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.46e-02
Iteration: 1900, time: 36.170 minutes, Δt: 894.066 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 1.923 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.42e-02
Iteration: 2000, time: 37.682 minutes, Δt: 882.857 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.009 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.46e-02
Iteration: 2100, time: 39.182 minutes, Δt: 902.168 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.095 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.56e-02
Iteration: 2200, time: 40.661 minutes, Δt: 907.427 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.183 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.44e-02
Iteration: 2300, time: 42.149 minutes, Δt: 908.291 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.269 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.41e-02
Iteration: 2400, time: 43.639 minutes, Δt: 891.806 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.356 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.41e-02
Iteration: 2500, time: 45.101 minutes, Δt: 876.494 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.442 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.38e-02
Iteration: 2600, time: 46.563 minutes, Δt: 849.315 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.529 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.40e-02
Iteration: 2700, time: 47.999 minutes, Δt: 857.372 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.615 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.43e-02
Iteration: 2800, time: 49.458 minutes, Δt: 878.759 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.701 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.36e-02
Iteration: 2900, time: 50.895 minutes, Δt: 877.943 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.788 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.46e-02
Iteration: 3000, time: 52.341 minutes, Δt: 851.701 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.875 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.36e-02
Iteration: 3100, time: 53.778 minutes, Δt: 857.786 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 2.961 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.46e-02
Iteration: 3200, time: 55.192 minutes, Δt: 833.145 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 3.047 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.37e-02
Iteration: 3300, time: 56.582 minutes, Δt: 835.902 ms, max(|pfft|) = 1.93537e-08 ms⁻¹, wall time: 3.133 minutes, advective CFL: 5.00e-01, diffusive CFL: 1.43e-02

But it seems like pfft does not update with time? Any ideas about this problem? Thanks!

[glwagner]

But it seems like pfft does not update with time? Any ideas about this problem? Thanks!

Are you recomputing it every time step?

[zwei961120]

I think I forgot this step. I tried to use Callback to recompute the pressure at each time step like:

pfft_aux = CenterField(grid)

function calculate_pfft!(sim)
    fft_solver = FFTBasedPoissonSolver(sim.model.grid)
    r = CenterField(sim.model.grid)
    set!(r, RHS_F)
    rc = fft_solver.storage
    rc .= interior(RHS_F)
    pfft = CenterField(sim.model.grid)
    solve!(pfft, fft_solver, rc)
    pfft_aux .= pfft
    return nothing
end

simulation.callbacks[:calculate_pfft] = Callback(calculate_pfft!)

But it also doesn't work...

[glwagner]

I'm not sure why that doesn't work but note that you don't need to create a new solver every time-step. You just need to set the RHS correctly and then invoke solve!:

function calculate_pfft!(sim)
    set!(r, RHS_F)
    rc = fft_solver.storage
    rc .= interior(RHS_F)
    solve!(pfft, fft_solver, rc)
    pfft_aux .= pfft
    return nothing
end

I would test this outside of time-stepping first. Run a simple simulation and then make a computation at the end. (Or, save data then load it back). Once that works, it's just bookkeeping to compute it online while the simulation is running.

[zwei961120]

I think I figured out the reason though I'm not entirely sure how it works. I just saved the RHS data and pressure data after calculating them, and leave the progress_message step for last like:

@inline function calculate_RHS_test(i, j, k, grid, U, V, u_prime, v_prime, w)
    
    Ux = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, U) # F, F, F -> C, F, F -> C, C, C
    Uy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, U)
    Uz = ℑxyᶜᶜᵃ(i, j, k, grid, ∂zᶠᶠᶜ, U)

    Vx = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, V)
    Vy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, V)
    Vz = ℑxyᶜᶜᵃ(i, j, k, grid, ∂zᶠᶠᶜ, V)

    ux = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, u_prime)
    vx = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, v_prime)
    wx = ℑyzᵃᶜᶜ(i, j, k, grid, ∂xᶜᶠᶠ, w)

    uy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, u_prime)
    vy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, v_prime)
    wy = ℑxzᶜᵃᶜ(i, j, k, grid, ∂yᶠᶜᶠ, w)

    return -2*(Ux*ux + Uy*vx + Uz*wx + Vx*uy + Vy*vy + Vz*wy) 
end

RHS = KernelFunctionOperation{Center, Center, Center}(calculate_RHS_test, grid, U, V, u_prime, v_prime, w)
RHS_F = Field(RHS)
compute!(RHS_F)

simulation.output_writers[:fieldsRHS] =
    NetCDFOutputWriter(model, (; RHS_F),
                     schedule = AveragedTimeInterval(10minute, window=10minute),
                     filename = "RHS.nc",
                     overwrite_existing = true)

fft_solver = FFTBasedPoissonSolver(grid)
pfft_aux = CenterField(grid)
pfft = CenterField(grid)

function calculate_pfft!(sim)
    rc = fft_solver.storage
    rc .= interior(RHS_F)
    solve!(pfft, fft_solver, rc)
    parent(pfft_aux) .= parent(pfft)

    return nothing
end

simulation.callbacks[:calculate_pfft] = Callback(calculate_pfft!)

simulation.output_writers[:fieldsPressure] =
    NetCDFOutputWriter(model, (; pfft_aux),
                     schedule = AveragedTimeInterval(10minute, window=10minute),
                     filename = "pfft_aux.nc",
                     overwrite_existing = true)

progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|pfft_aux|) = %.5e ms⁻¹, wall time: %s, advective CFL: %.2e, diffusive CFL: %.2e\n",
                                iteration(sim),
                                prettytime(sim),
                                prettytime(sim.Δt),
                                maximum(interior(pfft_aux)),
                                prettytime(sim.run_wall_time),
                                AdvectiveCFL(sim.Δt)(sim.model),
                                DiffusiveCFL(sim.Δt)(sim.model))
simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))

Then it works.

[glwagner]

Probably there is some scoping issue (eg overwritting pfft_aux) that means the pfft entering in your callback is not what you thought it was. Or something like. Also, note that callbacks are executed in the order that you add them to simulation.

[glwagner]

You need to call compute!(RHS_F) inside calculate_pfft!.

Think of it this way: nothing happens unless you explicitly ask for it. So in order for

rc .= interior(RHS_F)

to make sense, then you have to compute RHS_F first.