quad_util_mod rotations about dart HOT 16 OPEN

so this is scale, rotate, transform. Like a perspective transformation?

There's a fairly exhaustive description of the transform in the attached talk (CESM 2013).
The basic steps are to:

1. Define a polar coordinate system (radius, angle) on a plane that's tangent to the sphere at the grid point of interest, which is defined as the origin.
2. Map the 3 other corners of the grid box onto this plane, preserving the distances from the origin and the angles relative to one of the box sides.
3. Tranform these polar locations to cartesian locations. These 2 steps are where the "rotation" happens; the new box always has one side on the new x-axis.
4. Map the x-y coordinates to a unit square for use in the interpolation algorithm

raeder_Homme_DA.pptx

from dart.

this code is an adaptation of the original POP interpolation code. i will try to answer a few questions, but it might be better to have a video call with me and whoever is going to work on this issue when they're ready.

1. rotation: the original POP code interpolated the quads as they are. these quads can be quite deformed and none of the sides may be parallel to either the X or Y axis. when using this code on a regional roms case i saw artifacts in the interpolation. if the namelist item 'do_rotate' is set to true, each quad is rotated so one of the sides is parallel to one of the X or Y axes before doing the interpolation. the interpolation results looked much better at the expense of slightly more computation. it might be that do_rotate should be on all the time if the computation time isn't noticeable.

2. 45 degrees: i am unaware of any problems with quads with sides at a 45 degree angle. if there is a bug there someone needs to construct a test case to show it. i did a lot of testing and never ran into this.

3. the cam-se model_mod was developed before this code was moved into a separate module, as was POP. both should be using it but that had been low priority because the original code was working. for reducing code maintenance they all should use the same interpolation code.

4. i'm not sure what to do about missing values. it's a pain to check for them all over the code but for a model like pop there are land values and quads which have one corner as land should fail to interpolate.

from dart.

I think cam-fv won't be affected by this, because it's a rectangular lat-lon grid.
cam-se would be affected if its interpolation were replaced with quad_utils.

The cam-se rotation is not done at that L2531 point in the code.
It is implicitly embedded in the algorithm which calculates the transformation
from 4 points on Earth's surface to a unit square. That transformation is read from a file into variables
defined at L251; x_ax_bearings,..., and then used in subroutine unit_square_location (L2964).

That formulation may be one reason that quad_utils was developed separately from cam-se.
I believe the philosophy of the quad_utils is to handle each quad on the fly and deal with
rotations (and other problems?) as needed.
The philosophy in the cam-se is to generate the required transforms from the quads on the sphere
to the unit squares (used for the interpolation) before the assimilation. Each transform contains
a rotation so that diagonal points never have lats or lons that are close to each other (except maybe
at the poles, but that's dealt with robustly). The transforms are stored in a file,
which can be reused for multiple assimilations. This is an important reason that the SE assimilations
were no slower than the FV assimilations, despite the cubed sphere grid.
This formulation also deals with quads that have a distorted shape; narrow diamond, or almost triangular.
It does not deal with missing corners, but could. An advantage of this strategy is that
the quads with missing corners could be identified once, before assimilation. But that might not
work if the (ocean) grid depends on the model state (model levels changing depth?).

Another reason the cam-se formulation wasn't recycled may be that it had faded into obscurity
when the quad_utils work started. I was busy with the Reanalysis, so I didn't bring it up.

from dart.

oh ok thanks for the info

algorithm which calculates the transformation
from 4 points on Earth's surface to a unit square

so this is scale, rotate, transform. Like a perspective transformation?

from dart.

i should point out that the quad utils handles 3 cases:

1. completely regular rectangular grids. the current cam-fv code calls this.
2. rectangular grids where all the quads are proper rectangles but the spacing along the X and Y axis can vary. this supports a grid where, for example, the quads are smaller near a boundary and larger away. all the angles of the corners of the quads are still 90 degrees but the spacings of the sides can be larger or smaller.
3. logically a rectangular grid, but the location of all the corners need to be specified and the quads can be completely distorted. this code came from the code in POP written by jeff, and was generalized to have two basic routines - one that returns the field values at the quad corners, and one that interpolates to any interior point in that quad.

i misspoke in my previous post. the rewritten cam-se code was done after the quad utils mod was available and maybe it could have been used. it now computes corner values and weights and does a one-line interpolation based on area weights. there's a reference to the method in the code in the interpolate() routine.

from dart.

historical review (disclaimer: this is as i remember it):

someone (tim, ben?) came to me with a regional roms case where the grid was specified as a regular rectangular grid but the entire grid was rotated about 45 degrees to be parallel to the coastline. (maybe around the carolinas and georgia on the US east coast?) they were seeing interpolation artifacts.

in my test cases to reproduce this i needed to make the field data (the values at the quad corners) a bit noisy, and have the grid rotated relative to the X, Y axes. i saw 2 interpolation artifacts: values in the quad interior could exceed the range of values at the corners, and the interpolated results were not continuous across the quad boundaries.

doing a simple rotation on the quad before the interpolation fixed both problems.

there are test programs in the developer_tests/quad_interpolate directory. if those don't show these problems when you don't rotate the quad, let me know and i'll try to dig up some of my other tests.

from dart.

i misspoke in my previous post. the rewritten cam-se code was done after the quad utils mod was available and maybe it could have been used. it now computes corner values and weights and does a one-line interpolation based on area weights. there's a reference to the method in the code in the interpolate() routine.

I think that both are correct; the original SE code was developed before the quad utils (2012-13)
and the updated SE code was developed/reorganized after the quad_utils (2020).
Jeff opted to keep using the old SE code in the new.

from dart.

from dart.

MOM6 regional restart file example in /glade/work/hkershaw/DART/TestRuns/MOM6/regional/CARIB_012.001.3

from dart.

Related, but may split this off into another issue: Interpolation testing for CESM grids (and others) with quad_utils

Goal is to be confident that the quad_utils has no problems with the various CESM grids.

spec
branch:
https://github.com/NCAR/DART/tree/test_interpolate

from dart.

If we need a highly generalized interpolation package, especially for CESM,
then we might look into the ESMF interpolation package that CESM uses.
Mariana viewed it as a game-changer for them.
But maybe it's too cumbersome to implement at this late stage of DART's development.

from dart.

i'd be cautious. ESMF is implemented in a mix of Fortran and C++ which is a huge pain to compile. they do have a generalized regridding feature to move values from one grid to another. i'm not sure they have a "here is a single random observation point and a field, give me the value" function. even if they did, they'd require our fields to be much more complex than a linearized state vector fortran array.

from dart.

Is this the ESMF regridding? I think about this a lot, is regridding the same as interpolation?

from dart.

it's not the same as our interpolation. they do have to interpolate values from one grid to another, but they know the input grid and the output grid ahead of time. they expect to do a lot of values at once so they can add additional overhead up front by precomputing weights which can be stored. (those are the values to multiply each input point by to get the output values for each output grid point.) they can be stored so it's fast at runtime, but those weights are then set for those 2 input/output grids only. their code also can ensure that the total quantity in the output field is the same as the total quantity in the input field, so no roundoff gains or losses. conserving the overall total value is often important in long running climate applications.

for dart's use, we have a list of scattered observation values which are given, one at a time, to an input field and we need the output value back. the input field changes after each previous observation so the values can't be computed all at once.

from dart.

@nancycollins Good point about the language and compiling.
I would expect that we'd need to provide the locations of the grid corners
and the field values there, just like we do now. 5 years ago there was a mechanism
to interpolate the SE grid values onto an FV grid, for easy plotting of results.
I'm sure that that was just an interpolation, but may have used the spectral element functions
to generate a more faithful interpolation than bilinear would produce.

@hkershaw-brown They're related, but regridding should be a more comprehensive process.
Ideally it conserves a field when it is mapped onto a new grid and back again,
which doesn't happen with many interpolation recipes.

So ESMF may be a bad alley to go into.

from dart.

one other comment - the regridding works great for CESM coupling where the same grids are involved each coupling cycle. the surface values on the atmosphere grid are regridded onto the ocean grid each 15 or 30 minutes of model time, and vice versa. having a fast way to move values from one grid to another is great for them. also, they might only need 2d regridding because the coupler only has to deal with interface values which are 2d, not full 3d fields.

if we had a scenario where the observations never moved (so no GPS, no satellite obs, just a fixed set of ground stations, for example, which reported each assimilation time) then it might be feasible to investigate precomputing some information to speed the interpolation.

but we do all our forward operators in parallel already, and from most of the timings i remember that isn't the limiting factor in our execution time.

from dart.