quad_util_mod rotations

[hkershaw-brown]

Reported in the standup that the quad_utils_mod may have problems when the quad is rotated at 45 degrees.

model_mods that would be affected by this:

• cam-fv
• mom6 (unreleased)
• quad_utils is used by people with their own DART development for TIEGCM

Not confirmed (yet).

There is no record of the error (or fix) in the log: https://github.com/NCAR/DART/commits/main/models/utilities/quad_utils_mod.f90

There is the namelist option do_rotate (default .false.) https://github.com/NCAR/DART/blob/16595cad595b6349605c9bc8e021c7ecc2c2391a/models/utilities/quad_utils_mod.f90#L2179

The notes talk about the degenerate case where the grid is rotated 90 degrees: https://github.com/NCAR/DART/blob/16595cad595b6349605c9bc8e021c7ecc2c2391a/models/utilities/quad_utils_mod.f90#L2201-L2202

Mentioned at the standup (note I might not be recalling this correctly): cam-se has a fix for this rotation problem. Is this it? https://github.com/NCAR/DART/blob/16595cad595b6349605c9bc8e021c7ecc2c2391a/models/cam-se/model_mod.f90#L2531

Q. Why does cam-se not use quad_utils_mod? Q. What is the goal of rotating? Q. Why is the default do_rotate = .false.?

I think quad_interp_mod is a generalization (and replacement?) for the original POP model_mod interpolation code: https://github.com/NCAR/DART/blob/16595cad595b6349605c9bc8e021c7ecc2c2391a/models/utilities/quad_utils_mod.f90#L894-L897 Q. Can we switch out the POP model_mod code?

There are 21 fixmes is the quad_util_mod module.

Other quad_util_mod issues:

• see https://github.com/NCAR/DART/issues/249 for a note only only one of the quad_utils routines checking for missing values in the state.

• Note from pull #247

  quad_utils_mod.f90 uses  invals(4, nitems)
  cam-fv/model_mod.f90 uses quad_vals(ens_size,4) !< array of interpolated values
  

  which seems ripe for bugs such as https://github.com/NCAR/DART/issues/246. I'm not sure why the change in order was made, maybe should be reverted.

[nancycollins]

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.

[kdraeder]

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.

[hkershaw-brown]

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?

[nancycollins]

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.

[nancycollins]

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.

[kdraeder]

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

[kdraeder]

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.

[jlaucar]

Kevin is correct. I originally planned to use quad_utils in the SE rewrite, but I ran into difficulties reproducing the existing code and gave up given my time constraints. Jeff

On Wed, Feb 8, 2023 at 11:40 AM kdraeder @.***> wrote:

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.

— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/453#issuecomment-1423075517, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANDHUITBHUIIAJJWSA5RVCLWWPSCFANCNFSM6AAAAAAUUMLQC4 . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[hkershaw-brown]

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

[hkershaw-brown]

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

[kdraeder]

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.

[nancycollins]

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.

[hkershaw-brown]

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

[nancycollins]

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.

[kdraeder]

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

[nancycollins]

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.