Variable grid spacing in horizontal dimension (x or y) causes NonhydrostaticModel instantiation to fail

[pjrusello]

Using a variable grid spacing like the Ocean wind mixing and convection example in a direction other than z raises a MethodError when creating a NonhydrostaticModel.

Modified from OWM&C:

using Oceananigans
Nz = 24          # number of points in the vertical direction
Lz = 32          # (m) domain depth

refinement = 1.2 # controls spacing near surface (higher means finer spaced)
stretching = 12  # controls rate of stretching at bottom

# Normalized height ranging from 0 to 1
h(k) = (k - 1) / Nz

# Linear near-surface generator
ζ₀(k) = 1 + (h(k) - 1) / refinement

# Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))

# Generating function
z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

grid = RectilinearGrid(size = (32, Nz, 32),
                          x = (0, 64),
                          y = z_faces,
                          z = ( 0, 64 ))

model = NonhydrostaticModel( grid = grid )

I get the following error if I try to use the above code or put the variable grid spacing in the x-direction

ERROR: MethodError: no method matching PressureSolver(::CPU, ::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, CPU})

[tomchor]

I think this is indeed a limitation of Oceananigans at the moment. We have implemented a pressure solver for when the z coordinate is stretched, but not the other ones...

[pjrusello]

I think I remember reading that someplace in an issue or discussion thread... that would certainly explain the PressureSolver call failing.

[navidcy]

Could you comment on what version of Oceananigans you are using?

Perhaps we need to add a validation that spits out a meaningful warning to the user when this happens.

[simone-silvestri]

Indeed the pressure solver at the moment allows only a variable z-direction. This is because for a regular grid we can use a full FFT solve while with a singular stretched direction we can solve a tridiagonal system.

For more than one stretched direction we should use an iterative solver to solve the Poisson equation (the solvers are there but the infrastructure not quite yet). We should implement the iterative solver as a default for two or more stretched directions. And set up the FFT-tridiagonal solver as default when we have one stretched direction (either x, y or z)

In the meantime, as a hotfix, If you only need one stretched direction I would suggest you to rotate your problem to stretch it in z-direction

[glwagner]

I think the error is surprisingly clear actually! There's no pressure solver for grids that are stretched in x or y.

But it won't be much effort to add a more specific warning in the constructor for NonhydrostaticModel.

@simone-silvestri good point that the FFT/tridiagonal solver should work in other directions! Unfortunately though I think the tridiagonal solver is hard-coded to batch in i, j and solve in k.

[navidcy]

The error is not as clear, @glwagner.

ERROR: MethodError: no method matching PressureSolver(::CPU, ::RectilinearGrid{Float64, Periodic, Periodic, Bounded, Float64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, CPU})

does not immediately reads "no pressure solver for grids that are stretched in x or y". One needs to infer that from the long types of the function arguments.

[glwagner]

I guess I just read the first part "no method matching PressureSolver" and that made sense to me: there's no pressure solver for my input!

But you're right that you don't know which input is the problem (only I know that only the grid matters).

So yeah, I think some input validation for grid would help in the NonhydrostaticModel constructor. We only support RegRectilinearGrid and HRegRectilinearGrid. All other grid types are invalid.

[navidcy]

I guess I just read the first part "no method matching PressureSolver" and that made sense to me: there's no pressure solver for my input!

Fair enough! But perhaps not what would come in mind to everyone ;)

[tomchor]

I guess I just read the first part "no method matching PressureSolver" and that made sense to me: there's no pressure solver for my input!

Fair enough! But perhaps not what would come in mind to everyone ;)

I tend to agree with @navidcy. For that error to be clear a user would have to be familiar with Julia's multiple dispatch feature and what "method" means. Clear for us, but not super clear for a new Julia user.

@simone-silvestri good point that the FFT/tridiagonal solver should work in other directions! Unfortunately though I think the tridiagonal solver is hard-coded to batch in i, j and solve in k.

Is it worth to rotate x- and y-stretched grids so that the stretching is in the z direction and then use the tridiagonal solver behind the scenes? I have no idea how slow something like that would be on a GPU.

[simone-silvestri]

Is it worth to rotate x- and y-stretched grids so that the stretching is in the z direction and then use the tridiagonal solver behind the scenes? I have no idea how slow something like that would be on a GPU.

I think transposing at each time step would be quite expensive and memory-intensive.

[pjrusello]

I'm using v0.68.6, previously I was on v0.67.1 and saw the same error (not surprising).

I was hoping to use a stretched x-grid with a stretched z-grid while I play around with some 2D internal wave simulations. I'm still new to Julia and Oceananigans so I've been trying to simulate the classic vertically oscillating sphere or cylinder as my source. I was playing around with refining the mesh near the oscillation to see if the wave beams were a little cleaner. I can make due with stretched z-grids for now!

--
Regarding the MethodError, and while I am new to Julia, but maybe in this case a NotImplementedError would work better? I find Julia's stacktraces a bit long winded because of multi-dispatch and I obviously didn't figure out that variable grid spacing is only implemented for the vertical coordinate for the pressure solver. But getting a NotImplementedError when I tried would have been painfully obvious even to me. :-)

[glwagner]

I was hoping to use a stretched x-grid with a stretched z-grid while I play around with some 2D internal wave simulations. I'm still new to Julia and Oceananigans so I've been trying to simulate the classic vertically oscillating sphere or cylinder as my source. I was playing around with refining the mesh near the oscillation to see if the wave beams were a little cleaner. I can make due with stretched z-grids for now!

True!

You likely don't want a grid that's stretched in two directions for this application, because even if we had a valid solver for such a case, it would be much slower than the regular-x case. So even though you might be able to save a few grid points, you might not end up with shorter time-to-solution.

[glwagner]

But getting a NotImplementedError when I tried would have been painfully obvious even to me. :-)

NotImplementedError indicates that we plan to implement an appropriate pressure solver in the future (rather than ArgumentError, which is agnostic on that point). NotImplementedError might be appropriate...

[glwagner]

Is it worth to rotate x- and y-stretched grids so that the stretching is in the z direction and then use the tridiagonal solver behind the scenes? I have no idea how slow something like that would be on a GPU.

I think it might be more fruitful to generalize BatchedTridiagonalSolver. It's a little annoying (our fault for using i, j, k rather than CartesianIndex...) --- but we can probably figure out how to rotate indices in

Oceananigans.jl/src/Solvers/batched_tridiagonal_solver.jl

Line 91 in a314f1f

@kernel function solve_batched_tridiagonal_system_kernel!(ϕ, a, b, c, f, t, grid, p, args...)

Not a small refactor, so I'd argue not priority 1 right now.

[glwagner]

We added a more helpful error for this