Fixes bug in FourierTridiagonalPoissonSolver that swapped ΔzF and ΔzC

[glwagner]

This PR fixes a bug in the FourierTridiagonalPoissonSolver that almost certainly meant it produced incorrect answers. There were two concurrent bugs that cancelled each other out, which allowed tests to pass:

1. The Laplacian operator (previously called ∇², now called ∇²ᶜᶜᶜ) was incorrect, because "Face" and "Center" superscripts were swapped. In other words, the Laplacian operator was correct for an object located at (Face, Face, Face), rather than (Center, Center, Center) as it was being applied.
2. ΔzF and ΔzC were also swapped in derivation of the FourierTridiagonalPoissonSolver. This means that both the docs and the code were incorrect.

The tests passed because these two bugs effectively cancel each other out. This PR fixes both bugs. There still may be a lingering issue however, because the docs multiply the entire Poisson equation by the vertical grid spacing (which in the docs is written ΔzC, but should be ΔzF), including the source term. I didn't see immediately where to replace ΔzC with ΔzF in solve_poisson_equation!. The tests may pass anyways because there is no test that incompressibility is maintained on a stretched grid.

I also cleaned up the Laplacian operators a bit, and the grid spacing operators, and added a convenience function set_source_term! so that users don't have to know about the special formulation that FourierTridiagonalPoissonSolver uses.

TODO:

• Update the docs
• Test incompressibility on a stretched grid

[tomchor]

Nice, catch! What simulations use this solver? Just the ones with VerticallyStretchedGrid?

[glwagner]

Nice, catch! What simulations use this solver? Just the ones with VerticallyStretchedGrid?

That's right!

[francispoulin]

Thanks for catching this. The changes look good to me but I admit that I am a bit confuses as to when we have superscripts a vs c or f.

[glwagner]

The changes look good to me but I admit that I am a bit confuses as to when we have superscripts a vs c or f.

The superscript a refers to "Any" location. For example, the difference operator δxᶜᵃᵃ is independent of the y and z directions:

Oceananigans.jl/src/Operators/difference_operators.jl

Line 7 in 6bfc155

@inline δxᶜᵃᵃ(i, j, k, grid, u) = @inbounds u[i+1, j, k] - u[i, j, k]

Thus this operator is valid if y and z are respectively Center or Face.

When this is not true (for example, on a horizontally-curvilinear grid, the grid spacing in the x-direction depends on both x and y location), we need to explicitly specify the location of the operator / object.

[ali-ramadhan]

Thanks for catching this! Shows the importance of good notation I suppose 😅

Should we tag a new version with this PR? Vertically stretched grid users (I only know of @tomchor and @mukund-gupta) would probably want to upgrade.

[navidcy]

One question: Should we modify the tests so that they would be able to catch the bug? Or is that too much?

[glwagner]

One question: Should we modify the tests so that they would be able to catch the bug? Or is that too much?

We need to add a new test. I already changed the test that demonstrates the Poisson equation is solved correctly (with the correct Laplace operator as a diagnostic). However, something missing is a straightforward test that mass is conserved when using this solver on a stretched grid (not a regular grid, where this issue doesn't crop up). I'm working on that now.

[glwagner]

@ali-ramadhan @christophernhill it looks like the incompressible-in-time tests are failing for a hyperbolically VerticallyStretchedRectilinearGrid. However, a curious aspect of the test is that we compare maximum(|δ|) to 0, where δ is the divergence. Models that use VerticallyStretchedRectilinearGrid fail this test because maximum(|δ|) is greater than 5e-8 (its around 1e-4 for AdamsBashforth, and 3e-5 for RK3).

However, it looks like sum(δ) is smaller than machine precision; thus mass does not accumulate in time.

The larger-than-normal error only shows up when the grid is stretched; when the coordinate spacing in VerticallyStretchedRectilinearGrid is constant, the tests pass.

Should we modify this test? What do we want to test for and what do we expect? Presumably both non-accumulation of mass and a small divergence are both important...

[glwagner]

@ali-ramadhan @christophernhill disregard this previous comment. I found a 3rd bug in the divergence operator. I believe the tests pass now.

[glwagner]

Are we ok to merge this PR? It's important and might not be worth waiting for docs.