4th order or higher for center differencing

[francispoulin]

At the moment, derivative_operators.jl defines second order center differencing, which is great.

For ShallowWaterModel, I am splitting up the advection terms and the pressure. For advection I plan to use WENO5 and would similarly like a 4th order center differencing scheme for pressure.

The stencil for 4th order derivatives is [-2, 16, -16, 2]/16.

Based on that, I think a 4th oder derivative could simply be

@inline ∂xᶠᵃᵃ(i, j, k, grid, f::F, args...) where F<:Function = ( -2 f(i+1, j, k, grid, args...) + 16 f(i, j, k, grid, args...) - 16 f(i-1, j, k, grid, args...) + 2 f(i-2, j, k, grid, args...)) / (12 Δx(i, j, k, grid))

My initial ideal was to add a sub or superscript 4 somewhere to denote the fact that this was 4th order. @ali mentioned that maybe we might want to dispatch based on the order.

What do people think is the best way moving forward?

I don't think I would ever go higher than 4th since that makes for a much bigger stencil but if people wanted to, this could be done I suppose.

[ali-ramadhan]

Might be cool to be able to dispatch on the order so the order could be specified as part of the model, e.g.

model = ShallowWaterModel(grid=grid, order=4)

and it would make it easier to use the operators in other models. We could define new types like

struct SecondOrderCenteredDifference end  # one option
struct CenteredDifference{N} end          # another option

then dispatch on ::SecondOrderCenteredDifference or ::CenteredDifference{Val{4}} or we could dispatch on numbers via Val

julia> δ(i, A, ::Val{2}) = A[i] - A[i-1]
δ (generic function with 1 method)

julia> δ(i, A, ::Val{4}) = (-2A[i+1] + 16A[i] - 16A[i-1] +2A[i-2]) / 12
δ (generic function with 2 methods)

julia> δ(10, collect(1:20) .^ 2, Val(2))
19

julia> δ(10, collect(1:20) .^ 2, Val(4))
15.833333333333334

but might have to be careful to avoid performance regressions with Val.

That said it might take a non-trivial amount of refactoring to support and test dispatching on the operator order, at least for the incompressible model.

Maybe it makes sense for ShallowWaterModel to add support for 4th-order operators first (with or without dispatch, probably easier without first) and from there we can investigate how to generalize?

If we go all out and start supporting lots of different operators I wonder if it's worth looking into FiniteDiff.jl or FiniteDifferences.jl. Not sure what role these packages would play. From skimming the FiniteDifferences.jl README it seems that there are no higher-order non-allocating implementations between the two packages.

[glwagner]

I think it's a great idea to use types to control order of approximation. struct CenteredDifference{N} end is more cleanly generalizable. But note that the right concept might be "interpolating the derivative from cell-centered locations to face-centered locations", rather than the order of a finite difference approximation. This distinction will become relevant when we have arbitrarily stretched grids (arbitrary stretching is simple for shallow water models, so this could come sooner than we anticipated previously).

One simple possibility that might involve minimal code modifications is to add an order of approximation annotation to RegularCartesianGrid. We can then define "nth-order" interpolation and differencing operators with special notation (something like δⁿxᶠᵃᵃ) that types can opt into, such as IsotropicDiffusivity or the pressure term for ShallowWaterModel.

The primary application for incompressible models is diffusion operators, I think.

I'm not sure about high-order interpolation for other physics, such as Coriolis forces, etc.

Using grid might avoid the complexity associated with specifying an "order of approximation" for each aspect of the physics separately. I guess we already separate out advection (we could in principle come up with default advection schemes associated with grid's order of approximation...)

[glwagner]

Here's some additional thoughts on a design spec for "general" order-specified differentiation and interpolation operators:

"OperatorOrder" types for dispatch

I think it makes sense to use a generic type for controlling the order, such as

struct Centered{N} <: AbstractOperatorOrder end
struct LeftBiased{N} <: AbstractOperatorOrder end
struct RightBiased{N} <: AbstractOperatorOrder end

We might want to re-use the same types for both differentiation and interpolation so we may not want CenteredDifference. The "biased" operators would presumably only be used for interpolation, and by biased advection schemes. But we still want a unified implementation, since Centered{N} is used by both interpolation and differencing.

The new syntax for operators would be

δxᶠᵃᵃ(i, j, k, grid::AbstractGrid, order::AbstractOperatorOrder, c)

A global default stored in grid might make sense

We may want to add a type parameter to grid so we can write something like

δxᶠᵃᵃ(i, j, k, grid::AbstractGrid, args...) = δxᶠᵃᵃ(i, j, k, grid, grid.order, args...)

We want to allow the operators to be controlled independent of the global grid.order for advection schemes and for differencing pressure.

The above function, the default grid.order == Centered{2}, and specialization of the advection and pressure differencing operations will I think ensure that all existing code works as it did previously, I think.

Boundary conditions need to know about grid.order

Our implementation of boundary conditions assumes second-order operators. We would need to dispatch on grid.order in fill_halo_regions! to support higher-order differentiation. This is eminently doable, just a fair amount of busy work / arithmetic.

A nice bonus is that we would be able to change how high-order advection schemes mutate close to the boundary. For example, if we are using both fourth-order advection and fourth-order diffusion, we don't need to limit to second-order advection at the boundary, since halos will be filled consistent with fourth-order operations. In the implementation I am outlining, this requires using grid.order in our "topologically conditional interpolation":

https://github.com/CliMA/Oceananigans.jl/blob/master/src/Advection/topologically_conditional_interpolation.jl

For WENO5 with fourth-order diffusion, we just need to switch to fourth-order advection close the boundary. This is potentially a powerful increase in global accuracy of our methods.

[francispoulin]

Thanks @glwagner for sharing your thoughts above. I will have to think about this a bit more before before I can reply on the above.

However, I can say that in terms of ShallowWaterModel, we should be able to use the advection schemes for both advection, and also the divergence of mass flux in the continuity (conservation of mass) equation. That means the only place where we need to worry about these high order schemes if for pressure. At the moment, I believe that IncompressibleModel uses 2th order Finite Volume for the gradient of pressure, which is equivalent to 2nd order Center Differencing (I think). Doing 4th order Center Differencing in either model would be easy enough to do but it does mean that the method is not strictly Finite Volume, and therefore does not retain nice conservation properties, in case we have them.

I do have a question though. To use the advection schemes for the mass flux (evolution equation for height), we can use the tracer_flux functions since it's of the same form but c gets replaced by 1. If we want the function below but with c = 1 is it correct to just delete the c in both calls to the difference function?

    1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, mass_flux_x, advection, solution.uh, c) +
                             δyᵃᶜᵃ(i, j, k, grid, mass_flux_y, advection, solution.vh, c))
end

[glwagner]

At the moment, I believe that IncompressibleModel uses 2th order Finite Volume for the gradient of pressure, which is equivalent to 2nd order Center Differencing (I think).

As far as I know this is only true for a regular grid (or maybe this only matters for fluxes, not merely gradients...) ? We would like to support non-regular grids in the future in any case...

However, I can say that in terms of ShallowWaterModel, we should be able to use the advection schemes for both advection, and also the divergence of mass flux in the continuity (conservation of mass) equation.

Are we certain it makes sense to use upwind-biased differencing for the continuity equation?

Does it make sense to add a separate property to ShallowWaterModel (say, divergence or continuity) rather than reusing the property advection to control both advection operators in the momentum equation and the divergence operator in the continuity equation?

To use the advection schemes for the mass flux (evolution equation for height), we can use the tracer_flux functions since it's of the same form but c gets replaced by 1. If we want the function below but with c = 1 is it correct to just delete the c in both calls to the difference function?

You cannot use c=1 because most of the advection stencils attempt to interpolate c assuming that it is an array-like object that can be indexed into with c[i, j, k]. We can either 1) use a special type such as OneField that defines getindex(::OneField{FT}, args..) = one(FT) or 2) define functions _symmetric_interpolate_*(i, j, k, grid, scheme, c::Number) = c. The latter approach is probably preferable since it may be more broadly useful.

Doing 4th order Center Differencing in either model would be easy enough to do but it does mean that the method is not strictly Finite Volume, and therefore does not retain nice conservation properties, in case we have them.

I think we would like Oceananigans.jl to use finite volume methods, not finite difference methods. Perhaps @ali-ramadhan and @christophernhill can comment as well.

We can use a 4th order finite volume methods instead of a 4th order finite different method to obtain 4th order convergence without resorting to a finite difference method. We already have 4th-order interpolation of fields, so the only missing piece for a fully 4th-order method is 4th-order interpolation of field gradients... ? I believe we know how to do 4th order interpolation of gradients via comments in PR #1806.

We probably cannot support a 4th-order treatment of pressure for IncompressibleModel in the near future. The reason is that the FFT-based method we use to solve the pressure Poisson equation for IncompressibleModel explicitly requires 2nd-order discretization of pressure. We do not have an alternative implementation of a pressure Poisson solver for arbitrary-order methods. When we have a hydrostatic pressure solver for the 3D equations, a high-order pressure discretization could be interesting to explore. Alternatively, we could add a more general method for solving the pressure Poisson equation. But this is a non-trivial amount of work.

I think the principal challenge of a high-order implementation is code design. The finite volume stencils themselves are not so hard to write down. But I think we want a clean code design that minimizes bugs and maximizes the utility of the high-order stencil.

[francispoulin]

Can we use advection for mass evolution equation?

Let me try and convince you that using advection for continuity makes sense. If I don't succeed then clearly I need to put more thought into it. The governing equation can be written in two forms $$ \partial_t = -\partial_x (h u) - \partial_y (h v) = - \partial_x U - \partial_y V $$ If we were to solve this in terms of the velocities (u,v), then using advection is clearly a good idea as we are advecting h by the velocity. In the case of upwinding, we pick the direction based on the sign of u (or v). We are using the transports instead of velocities but it is very similar to the advection of a tracer $\nabla \cdot (u c)$ except that we use c = 1. Given that the physics is the same in both cases, I think we should be using the advection scheme for this term.

The evolution equations for h and c have a lot in common, which is helpful here. The one major difference is that h is an active tracer in the sense that it feeds back on the momentum. Even though that's the case, I don't think this should affect how we discretize the advection of h, even though it is disguesed with our use of the transport variable.

If you agree with my reasoning I don't think we need to do anything different for the evolution of height equation, we might just need to have some new functions that account for this, which are essentially equivalent to c = 0. Maybe what you suggest is the way to proceed?

High Order Finite Volume

I agree that keeping things finite volume is a good idea, and that is how I would like to proceed.

Good to know that IncompressibleModel won't generalize to higher order.

I agree that the 4th-order pieces are there and just need to be assembled. The paper you cited above is very helpful but needs to be modified slightly as they present the method for a non-staggered grid, but ours is staggered. That is not problem, I just need to work through the details. Unless they have been done already?

[glwagner]

Good to know that IncompressibleModel won't generalize to higher order.

Just the pressure! We can definitely generalize diffusion schemes to high-order.

[glwagner]

The paper you cited above is very helpful but needs to be modified slightly as they present the method for a non-staggered grid, but ours is staggered. That is not problem, I just need to work through the details. Unless they have been done already?

The staggering shouldn't affect the stencils themselves --- it only affects how the stencils are applied (I think).

[glwagner]

We are using the transports instead of velocities but it is very similar to the advection of a tracer $\nabla \cdot (u c)$ except that we use c = 1. Given that the physics is the same in both cases, I think we should be using the advection scheme for this term.

Ah, I understand the math, I was more asking whether upwinding makes sense as a numerical approximation for this problem in particular. But now I realize that my question is off the mark, since our scheme only upwinds the advected quantity. Symmetric interpolation is used for the advecting velocity. So if you specify advection = WENO5(), you'd get 4th order interpolation in the continuity equation under this assumption. Perhaps this makes sense... my only concern is that it might not be expected behavior for users.

[glwagner]

There could be consistency requirements between the treatment of pressure / height gradients in the momentum equation, and transport gradients in the continuity equation. @francispoulin do you know of any requirements here or is this not a concern?

[francispoulin]

The paper you cited above is very helpful but needs to be modified slightly as they present the method for a non-staggered grid, but ours is staggered. That is not problem, I just need to work through the details. Unless they have been done already?

The staggering shouldn't affect the stencils themselves --- it only affects how the stencils are applied (I think).

I think it affects the coefficients if you are approximating the value at a cell centre or at a face, but I need to think this through more carefully to convince myself.

[francispoulin]

There could be consistency requirements between the treatment of pressure / height gradients in the momentum equation, and transport gradients in the continuity equation. @francispoulin do you know of any requirements here or is this not a concern?

Hmm, I don't believe so. If whatever finite volume scheme we use for the evolution of height, we will ensure mass is conserved. This presumably variable height will drive the motion, and that parts seems very distinct to me. However, if you wanted to ensure global conservation of energy, or something, maybe one needs to be careful about this? I don't believe this is a concern though.

[glwagner]

However, if you wanted to ensure global conservation of energy, or something, maybe one needs to be careful about this?

Yes this, or numerical pitfalls (for example, its not possible to use a biased stencil for the pressure term as @navidcy found, I think). If there's no a priori reason to worry, we might as well test it and find out.

[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.