halos and interpolation

[pjrusello]

I have a couple of questions related to interpolation of results near boundaries and the effect of halo regions.

If I use Oceananigans's interpolate method for Fields, it fails if the specified physical location is too close to a boundary with a bounds error since interpoalte tries to index outside the array bounds. interpoalte also fails if the interpolation location is equal to one of the center or face coordinates (i.e., no interpolation would be needed) since interpolate is still trying to do a linear interpolation rather than falling back to a table lookup (I think).

For the case of interpolating near boundaries though, I was wondering if it is advisable to use the halo region data for interpolation?

From what I understand the halos are there for calculations at the boundaries to maintain higher order accuracy throughout the domain. Is there any reason not to utilize them for lower order linear interpolation to align velocities and scalars in post-processing?

[simone-silvestri]

Hello!

If you interpolate near the boundaries halo points should be used to interpolate. Take care that in case of Bounded domains only one halo point is filled (since we cannot preserve high order in bounded domains).

The interpolation on a Face and Center is, indeed, still doing a linear interpolation that but that should return the correct value regardless.

Can you provide an example for what setting make the interpolation fails?
If I run this script the interpolation I am not able to reproduce the error

grid = RectilinearGrid(size = (5, 1, 1), x = (0, 5), z = (0, 1), y = (0, 1), topology = (Bounded, Bounded, Bounded))
field = CenterField(grid)
set!(field, (x, y, z) -> x)
fill_halo_regions!(field)

julia> interpolate(field, 1.0, 0.5, 0.5)
0.9999999999999998

julia> interpolate(field, 2.0, 0.5, 0.5)
2.0

julia> interpolate(field, 0.0, 0.5, 0.5)
0.0

julia> interpolate(field, 0.1, 0.5, 0.5)
0.10000000000000014

[glwagner]

You can obtain a Field by indexing into FieldTimeSeries at a particular time index:

u_time_series = FieldTimeSeries(langmuir_turbulence_fields.jld2", "u" )
u1 = u_time_series[1]

What are you trying to do?

We should certainly add some documentation for FieldTimeSeries...

[pjrusello]

I wanted to see if the interpolation issue is with how Field and FieldTimeSeries are handled by interpolate. I suspect this is by design, but this behavior isn't what I expected (u_loaded is the FieldTimeSeries above, u is the Field):

interpolate( u_loaded[ 1 ], 124, 38, -33 )
-0.0032558434177190065

julia> interpolate( u, 124, 38, -33 )
-0.0033981014731492613

julia> interpolate( u_loaded, 124, 38, -33 )
ERROR: BoundsError: attempt to access 32×32×32×1 FieldTimeSeries...

I'll chalk this up mostly to me not understanding the difference between a Field and a FieldTimeSeries (or noticing until now I was dealing with two different object types).

Ultimately, I was just trying to figure out the best way to interpolate u and v onto a common grid and was checking to see how Oceananigans or the examples handle this. Which led me down the rabbit hole of seeing what happens at boundaries and how to deal with halo regions and finding the above behavior.

[glwagner]

Ah! interpolate wasn't designed to work with FieldTimeSeries! This is an interesting application though.

interpolate leaves something to be desired, but it depends on what you are trying to use it for. If you would just like u and v on the same grid, then you can probably use interpolate. However, you can also directly output u and v on the tracer grid; for example:

uc = @at (Center, Center, Center) u
vc = @at (Center, Center, Center) v

and then pass uc and vc to an output writer.

But often you don't actually want u and v on the same grid. What is your application?

[pjrusello]

Cool. That is a super handy way to shift data output... Internally what's happening with the @at call? I would guess an interpolation, but is it higher order?

Beyond just a general interest in working with models more (my background is primarily experimental work), in this case we needed some realistic ocean surface velocity data to test a remote sensing model. Nothing that necessarily requires a high-fidelity solution, I'm mostly using it as an excuse to work with numerical models more.

[glwagner]

Ok, that sounds like an appropriate application!

If you inspect uc and vc you'll find that they are UnaryOperation, which is a type of "abstract operation". Abstract operations represent operations (so they are "lazy"; we don't interpolate anything until we evaluate the operation, which typically occurs right before writing the data to disk). It's similar to, for example, something like u2 = u^2.