`FourierTridiagonalPoissonSolver` for Flat dimensions

[tomchor]

I'm working on a tilted 2D bottom boundary layer example for the docs for which I'm trying to use a Flat dimension. I then get this warning:

┌ Warning: FourierTridiagonalPoissonSolver is probably wrong for topologies that contain Flat dimensions.
└ @ Oceananigans.Solvers ~/repos/Oceananigans.jl/src/Solvers/fourier_tridiagonal_poisson_solver.jl:32

and then the simulation fails with a NaN in the first iteration. This doesn't happen if I don't use Flat and just manually set Ny=1.

A few questions:

• So, what exactly is missing for FourierTridiagonalPoissonSolver to deal with Flat dimensions.
• Are there other limitations that Flat dimensions have? We should probably mention those in the docs.
• Should we even allow a model to be built in these situations if we know it's gonna lead to NaNs?

[francispoulin]

Hmnm, I thought there existed a Poisson solver that worked for each topology. Maybe I'm wrong?

I looked at your code and see you want to use (Periodic, Flat, Bounded). I also checked in test_poisson_solver.jl and see this does test that particular topology but the test is for instantiating.

Is it possible that the model can instantiate but not solve the poisson problem for this topology? If yes then do we want a stronger test?

[glwagner]

I think Flat sets some grid spacings to 0 and other grid metrics to 1 (so that areas and volumes are correct).

We'll have to read through the code to figure out where the issue is:

https://github.com/CliMA/Oceananigans.jl/blob/master/src/Solvers/fourier_tridiagonal_poisson_solver.jl

[glwagner]

As for issues with other topologies, we could possibly use the "forced flow, fixed-slip" convergence test with 2D slices oriented in various directions to uncover issues:

https://github.com/CliMA/Oceananigans.jl/blob/master/validation/convergence_tests/src/ForcedFlowFixedSlip.jl

[glwagner]

This problem may be an issue with the vertically stretched grid rather than the FourierTridiagonalPoissonSolver.

The mixed FFT/Tridiagonal Poisson solver uses similar constructs as the FFTBased Poisson solver. The two ingredients are 1) a Fourier transform in the Flat direction, and 2) computation of the eigenvalues:

https://github.com/CliMA/Oceananigans.jl/blob/5fbd8cd20c5db8e9b11b6175984e7592a08fc874/src/Solvers/poisson_eigenvalues.jl#L8-L11

https://github.com/CliMA/Oceananigans.jl/blob/5fbd8cd20c5db8e9b11b6175984e7592a08fc874/src/Solvers/poisson_eigenvalues.jl#L31

The first eigenvalue is 0, so the Flat case is correct. It's also used for the fully FFT-based solver and we know it works in that context. The transforms must be fine too.

The grids differ though in how grid metrics are computed, so I suspect the problem is there. I think we've also had issues with HydrostaticFreeSurfaceModel and VerticallyStretchedRectilinearGrid with Flat topologies, but I haven't looked into it systematically.

We might be able to convert the internal wave setup dynamics test to a vertically bounded domain and put it on a stretched grid to test these issues.

[tomchor]

I think I get the idea, but I still can't imagine what could be wrong with the vertically stretched grids themselves. They seems pretty straightforward. Could you clarify what specific metrics you're talking about that are different?

We might be able to convert the internal wave setup dynamics test to a vertically bounded domain and put it on a stretched grid to test these issues.

This seems like a good idea.

[glwagner]

I think I get the idea, but I still can't imagine what could be wrong with the vertically stretched grids themselves. They seems pretty straightforward. Could you clarify what specific metrics you're talking about that are different?

We might be able to convert the internal wave setup dynamics test to a vertically bounded domain and put it on a stretched grid to test these issues.

This seems like a good idea.

I don't know which ones are problematic --- once we find those, the problem is solved.

The answer might be obvious; it looks like we only redefine Flat metrics for RegularRectilinearGrid:

https://github.com/CliMA/Oceananigans.jl/blob/326f22aff244ad1f9d0778b9d06184f348db211b/src/Operators/spacings_and_areas_and_volumes.jl#L57-L62

We may just have to use AbstractGrid rather than RegularRectilinearGrid. I'm not sure why the above functions are specific to RegularRectilinearGrid in the first place. @francispoulin you helped with this, right? Do you know?

We're missing many of the horizontal metrics too so I think we should add those, otherwise Flat won't work with curvilinear grids, either. EDIT: definitions for flat area metrics are probably not very useful / represent edge cases that need special consideration. So we might just focus on flat linear metrics for now (dx, dy at various locations). Right now additional horizontal linear metrics are irrelevant because there are no rectilinear grids that are stretched in horizontal directions. However, we will need these once we have general stretched rectilinear grids and it might be good to avoid confusion like what's happening in the present issue...

I'm putting together a test for rectilinear grids that catches this. I suggest we fix the problem after we have the test. If anyone has more information about Flat please chime in.

[glwagner]

On https://github.com/CliMA/Oceananigans.jl/pull/1865 we find that both the incompressible and hydrostatic model are instant NaN-factories with vertically stretched grid and flat dimensions. Therefore I think the issue is not the Poisson solver. But we'll see.