Evaluating volume-averages of functions of x, y, z, t with higher than first-order accuracy

[glwagner]

Our user interface makes extensive use of "continuous" function of x, y, z, t. Throughout our code, we evaluate these functions in control volumes using a first-order accurate method. In other words, the volume-average of a function is approximated by its value at the barycenter of a cell, as in FunctionField:

Oceananigans.jl/src/Fields/function_field.jl

Lines 45 to 46 in 068a7ae

	@inline Base.getindex(f::FunctionField{X, Y, Z, <:Nothing}, i, j, k) where {X, Y, Z} =
	f.func(xnode(X, i, f.grid), ynode(Y, j, f.grid), znode(Z, k, f.grid))

In the above code, getindex is (implicitly) assumed to return the volume average of the function f in the control volume labeled i, j, k. This volume average is approximated by simply evaluating the function x, y, z, t at the barycenter of the control volume in question. This is a valid, but only first-order approximation.

The pitfall of this method is that its accuracy degrades significantly when f varies on a scale comparable to the grid scale. Another issue is that significant differences can arise between an analytically prescribed domain-integrated tracer, for example, and the discrete domain-integrated tracer after using set!.

To fix such problems, higher-order approximation methods for volume averages of functions may be helpful. When using a higher-order approximation method, a function of x, y, z, t would be evaluated multiple times per grid point. This would have an immediate obvious use in set!(field, func::Function).

Such a method may also be useful for forcing functions, boundary conditions, and prescribed background fields (once #960 is resolved). While evaluating a function multiple times per grid point would make a simulation more computationally expensive, the extra cost may be negligible in simulations that are dominated by the cost of memory accesses.

A high-order method we might consider is multi-dimemsional Gaussian quadrature. We could add the necessary quadrature data to FunctionField, and expand the importance of FunctionField by using it within ContinuousForcing, BoundaryFunction, and set!.

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