Implement do_cube_transform in make_hgrid

[lharris4]

Is your feature request related to a problem? Please describe. The current 'do_schmidt' option in make_hgrid (taken from FV3) rotates the south pole (center of tile 6 on the cubed sphere) to the target_lat and target_lon. This has the effect that tile 6 is 'upside down' in the new orientation, with south at the top. Since tile 6 is the the highest-resolution grid when doing grid stretching (stretch_fac > 1.0), is where nested grids are usually placed, and is the foundation for the regional domain, the native grids are all upside down. This is not only aesthetically unpleasing, but makes placement of nested grids extremely difficult since [ij]offset are given in terms of local grid coordinates, and so counterintuitively i points westward and j points southward. Plus this has led to significant confusion for users of the nest and the regional domain.

Describe the solution you'd like A revision to do_schmidt, called do_cube_transform, has been implemented within FV3. This uses the routine fv_grid_utils_mod::cube_transform(). For idealized tests in either SHiELD or solo_core this works very well to create a tile 6 which is oriented as expected, with north at the top. I would like this option to be implemented within make_hgrid as well for generating the supergrid. With the supergrid and grid mosaic, the input datasets can be generated to run full-physics SHiELD or UFS. The new make_hgrid can then easily be inserted within preprocessing workflows.

The name is also much more indicative than do_schmidt---who is Schmidt and what does he have to do with grids? (Someone who appears to have written a single article in 1977 in a German journal.)

Describe alternatives you've considered I had thought of having FV3 write out supergrid tiles itself, replicating the capabilities of make_hgrid within the dycore itself. I realized this may be quite a bit of work since I would have to define a supergrid within the solver, including defining a domain structure on the supergrid, and then create the I/O routines to write out to disk. This may be a powerful way of rapid testing of new grids though.

[ngs333]

@lharris4 Hi! I have noticed that your request may already be in make_hgrid. To me it appears that the FV3 do_cube_transform you link above is equivalent to the NCTools cube_transorm() function (in file tools/make_hgrid/create_gnomonic_cubic_grid.c (lines 1396ff) or https://github.com/NOAA-GFDL/FRE-NCtools/blob/master/tools/make_hgrid/create_gnomonic_cubic_grid.c ).

When create_hgrid is called, the use of the cube_transform may be selected with the option "do_cube_transform" instead of the option "do_schmidt". It is possible, however, that it is has not been used very much (it at all) and it appears that none of the unit tests in directory FRE-NCTools/t use the do_cube_transform option.

[lharris4]

Hi, Miguel. Thanks for pointing this out. Indeed I do see that do_cube_transform is implemented and does work, for the most part.

I ran this command

make_hgrid --grid_type gnomonic_ed --nlon 48 --do_cube_transform --stretch_factor 5 --target_lat 40 --target_lon 20

and produced a stretched cubed-sphere grid, which has the smallest grid cells on tile 6, and indeed has north at the 'top'. So this is very convenient, and much easier to consider than the existing do_schmidt. Compare, for instance:

make_hgrid --grid_type gnomonic_ed --nlon 48 --do_schmidt --stretch_factor 5 --target_lat 40 --target_lon 20

One problem with the --do_cube_transform is that the 'area' array is incorrect. Instead of showing areas that match what I would expect from dx and dy, we instead see the same (incorrect) area array on all six tiles, as shown in the attached files. (They show dx and area, in order, and show tiles 1 to 6 with 1 at the bottom). This isn't a problem for FV3 and may not be for chgres either, although I am not sure if some other program uses the area program.

Thanks, Lucas

[image: image.png]

[image: image.png]

On Fri, Sep 23, 2022 at 12:10 PM Miguel R Zuniga @.***> wrote:

@lharris4 https://github.com/lharris4 Hi! I have noticed that your request may already be in make_hgrid. To me it appears that the FV3 do_cube_transform you link above is equivalent to the NCTools cube_transorm() function (in file tools/make_hgrid/create_gnomonic_cubic_grid.c (lines 1396ff) or https://github.com/NOAA-GFDL/FRE-NCtools/blob/master/tools/make_hgrid/create_gnomonic_cubic_grid.c ).

When create_hgrid is called, the use of the cube_transform may be selected with the option "do_cube_transform" instead of the option "do_schmidt". It is possible, however, that it is has not been used very much (it at all) and it appears that none of the unit tests in directory FRE-NCTools/t use the do_cube_transform option.

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/FRE-NCtools/issues/179#issuecomment-1256402349, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVC7UDQM6URBJ2637F3V7XI5DANCNFSM6AAAAAAQHF3ATA . You are receiving this because you were mentioned.Message ID: @.***>

[lharris4]

Images shown here

Images are different sizes since dx is staggered and area is cell-mean.

[lharris4]

One last thing: the dx values use the 'low' range in ncview, so they emphasize low values (purple). The area values use a linear color scale.

Lucas

[ngs333]

@lharris4 Could there be a problem with ncview when multiple files are read together. I see images similar (not exact) to what you see when the six nc tile files are read together, but if ncview is used with just one file at a time, ten the images are different. Here is the result, using do_cube_transform options, for tile6, dx:

tile6 area: tile1 dx tile1 area

[lharris4]

Hi, Miguel. Your dx looks right, but area appears to be wrong. On tile 6, area should be smallest in the center, not on the edge, and should look different than tile 1.

Thanks, Lucas

On Fri, Sep 23, 2022 at 4:51 PM Miguel R Zuniga @.***> wrote:

@lharris4 https://github.com/lharris4 Could there might be a problem with ncview when multiple files are read together. I see images similar (not exact) to what you see when the six tile nc files are read together, but ncview is used with just one file at a time, the images are different. Here is the result, using do_cube_transform options, for tile6, dx: [image: tile6_dx] https://user-images.githubusercontent.com/42479054/192054352-a47abb55-507e-47f9-8aa0-e89385cb2f6b.png

tile6 area: [image: tile6_area] https://user-images.githubusercontent.com/42479054/192054202-8e36250d-4f84-427a-a72f-885a2c66c017.png tile1 dx [image: tile1_dx] https://user-images.githubusercontent.com/42479054/192054258-d3255caa-f692-4c00-9897-fe879348f517.png tile1 area [image: tile1_area] https://user-images.githubusercontent.com/42479054/192054290-99e59294-25dc-4212-87d0-f8166222e59b.png

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/FRE-NCtools/issues/179#issuecomment-1256660961, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVA52SWYSTX5YUWTL7TV7YJ3NANCNFSM6AAAAAAQHF3ATA . You are receiving this because you were mentioned.Message ID: @.***>

[ngs333]

@lharris4 Thanks Lucas, I'll see if I can find the bug.

[ngs333]

@lharris4 Looks like make_hgrid had a simple bug in code to calculate the areas when the cube transform is applied. I am not totally sure this was the only bug.

In comparison to do_schmidt, the cube_transform x and y fields for tile 6 look perfectly flipped (or rotated 180 degrees). Here they are for the cube transform (y values are the second image):

Could there still be an issue with the area calculation?. For tile 6 it does now finally look different than tile 1, but it looks the same as the area produced by the do_schmidt transform, and the smallest values are in the center:

[lharris4]

Hi, Miguel. I believe this is the correct answer. The area pattern looks right, with largest areas in the center and smallest in the corners. Also the latitudes on tile 6 in do_cube_transform should be flipped 180deg from those in do_schmidt.

Thanks, Lucas

On Tue, Oct 11, 2022 at 3:21 PM Miguel R Zuniga @.***> wrote:

@lharris4 https://github.com/lharris4 Looks like make_hgrid had a simple bug in code to calculate the areas when the cube transform is applied. I am not totally sure this was the only bug.

In comparison to do_schmidt, the cube_transform x and y fields for tile 6 look perfectly flipped (or rotated 180 degrees). Here they are for the cube transform: [image: tile_6_x_ct] https://user-images.githubusercontent.com/42479054/195178922-4177d31e-8f8b-401f-b235-a77dbcae9baa.png [image: tile_6_y_ct] https://user-images.githubusercontent.com/42479054/195178936-9b18ece4-0650-4fef-9350-6c8e2664f03d.png

Could there still be an issue with the area calculation?. For tile 6 it does now finally look different than tile 1, but it looks the same as the area produced by the do_schmidt transform, and the smallest values are in the center:

[image: tile_6_area_ct] https://user-images.githubusercontent.com/42479054/195180081-4fbcd441-beec-4897-8a5a-80e8c58136b2.png

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/FRE-NCtools/issues/179#issuecomment-1275165955, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVGUFG6RZ6AZOSGH273WCW44ZANCNFSM6AAAAAAQHF3ATA . You are receiving this because you were mentioned.Message ID: @.***>

[ngs333]

@lharris4 Hi Lucas, I’ve been able to do some comparisons of the SHIELD grid output vs. files generated with the do_transform option it the new version of make_hgrid. With the help of @laurenchilutti we ran on Gaea the SHIELD experiment C128r3.solo.TC. That experiment was compiled with gcc and its input.nml had these parameters”

       do_cube_transform = .true. t
       target_lat = 17.5
       target_lon = 172.5
       stretch_fac = 3.0

The experiment produced files grid_spec_tile[1 - 6].nc files. These files were then converted (using the NCTools work we did last year for making supergrid/mosaics from internal FV3 grid files,).

Finally for the sake of comparison the equivalent files were directly generated from scratch by the new make_hgrid:

make_hgrid --grid_type gnomonic_ed --nlon 256 --do_cube_transform  --stretch_factor 
   3 --target_lat 17.5 --target_lon 172.5

A final conversion was required in order to compare with nccmp, and that was a NetCDF type conversion type via nco’s ncks. The nccmp comparison results were that all data values in each of the six grid tile files match with a tolerance of 1.0e-7% for fields x and y, and a tolerance 1.5e-4% for fields dx and dy. I though that it would be a little better than that, but given the differences in language (C vs Fortran), in Hardware (Gaea vs Analysis) and possible compiler options, plus all the conversions, maybe the differences are acceptable? Thanks MZ

[lharris4]

Hi, Miguel. Thanks for working on this. I think a tolerance of one part in 10^9 is acceptable, although it is a bit concerning that it is higher for dx and dy. As long as the pattern is fairly uniform worldwide (could you show a plot of differences over a tile?) it shouldn't be a problem.

Is make_hgrid using single or double precision, and is it possible to compile make_hgrid with higher precision still? If so, then this would be a good way to test whether it is roundoff/arithmetic errors, or an actual implementation bug. If the errors go down with higher precision, then it is just the arithmetic.

Lucas

On Thu, Nov 17, 2022 at 1:53 PM Miguel R Zuniga @.***> wrote:

@lharris4 https://github.com/lharris4 Hi Lucas, I’ve been able to do some comparisons of the SHIELD grid output vs. files generated with the do_transform option it the new version of make_hgrid. With the help of @laurenchilutti https://github.com/laurenchilutti we ran on Gaea the SHIELD experiment C128r3.solo.TC. That experiment was compiled with gcc and its input.nml had these parameters”

   do_cube_transform = .true. t

   target_lat = 17.5

   target_lon = 172.5

   stretch_fac = 3.0

The experiment produced files grid_spec_tile[1 - 6].nc files. These files were then converted (using the NCTools work we did last year for making supergrid/mosaics from internal FV3 grid files,).

Finally for the sake of comparison the equivalent files were directly generated from scratch by the new make_hgrid:

make_hgrid --grid_type gnomonic_ed --nlon 256 --do_cube_transform --stretch_factor

3 --target_lat 17.5 --target_lon 172.5

A final conversion was required in order to compare with nccmp, and that was a NetCDF type conversion type via nco’s ncks. The nccmp comparison results were that all data values in each of the six grid tile files match with a tolerance of 1.0e-7% for fields x and y, and a tolerance 1.5e-4% for fields dx and dy. I though that it would be a little better than that, but given the differences in language (C vs Fortran), in Hardware (Gaea vs Analysis) and possible compiler options, plus all the conversions, maybe the differences are acceptable? Thanks MZ

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/FRE-NCtools/issues/179#issuecomment-1319064828, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVFXDSFXDNP2YAO6GZTWIZ5LFANCNFSM6AAAAAAQHF3ATA . You are receiving this because you were mentioned.Message ID: @.***>

[ngs333]

@lharris4 I just realized that the SHiELD experiment was writing the FV3 grid lats and longs as degrees of seven digits (and the decimal point)! I'll see if I can find a parameter to change that (to 15 digits, hopefully).

[ngs333]

@lharris4 I think the 32 bit grid data output issue is resolved, but in the course of testing I ran into two different issues in other parts of the code, which I describe way below.

First, in experiment C128r3.solo.TC I changed the grid spec output format from 2 to 1. This changes the resultant diag_table, and the grid spec is now printed with 15 digits (plus decimal). With this all the x,y field values are now in tolerance of 10^-9 %. I.e. this check on every file shows no discrepancies:

nccmp -df --Tolerance=1e-9 --variable=x,y grid.tile1.nc ../fv3_sgrids/grid_spec.tile1.nc

And for fields dx,dy, the Tolerance ls 2.0*10^-9%.

However, for grid cell area there is an issue, most values are in Tolerance of 1.0e-4, but for tiles 1 and 3 there is a small number of points with a very large (up to 27%) error. Here is the output for one comparison:

gaea14$ nccmp -df --Tolerance=1e0 --variable=area grid.tile3.nc ../fv3_sgrids/grid_spec.tile3.nc
DIFFER : VARIABLE : area : POSITION : [125,201] : VALUES : 6.35747e+09 <> 6.42859e+09 : PERCENT : 1.10627
DIFFER : VARIABLE : area : POSITION : [126,199] : VALUES : 6.63873e+09 <> 6.55313e+09 : PERCENT : 1.28942
DIFFER : VARIABLE : area : POSITION : [126,201] : VALUES : 6.36169e+09 <> 6.50816e+09 : PERCENT : 2.25049
DIFFER : VARIABLE : area : POSITION : [126,202] : VALUES : 6.16146e+09 <> 6.03806e+09 : PERCENT : 2.00288
DIFFER : VARIABLE : area : POSITION : [127,199] : VALUES : 6.63976e+09 <> 6.72003e+09 : PERCENT : 1.19451
DIFFER : VARIABLE : area : POSITION : [127,200] : VALUES : 6.43222e+09 <> 6.18753e+09 : PERCENT : 3.80405
DIFFER : VARIABLE : area : POSITION : [127,201] : VALUES : 6.36265e+09 <> 8.61327e+09 : PERCENT : 26.1296
DIFFER : VARIABLE : area : POSITION : [128,199] : VALUES : 6.63976e+09 <> 6.72003e+09 : PERCENT : 1.19451
DIFFER : VARIABLE : area : POSITION : [128,200] : VALUES : 6.43222e+09 <> 6.18753e+09 : PERCENT : 3.80405
DIFFER : VARIABLE : area : POSITION : [128,201] : VALUES : 6.36265e+09 <> 8.61327e+09 : PERCENT : 26.1296
DIFFER : VARIABLE : area : POSITION : [129,199] : VALUES : 6.63873e+09 <> 6.55313e+09 : PERCENT : 1.28942
DIFFER : VARIABLE : area : POSITION : [129,201] : VALUES : 6.36169e+09 <> 6.50816e+09 : PERCENT : 2.25049
DIFFER : VARIABLE : area : POSITION : [129,202] : VALUES : 6.16146e+09 <> 6.03806e+09 : PERCENT : 2.00288
DIFFER : VARIABLE : area : POSITION : [130,201] : VALUES : 6.35747e+09 <> 6.42859e+09 : PERCENT : 1.10627

This is almost certainly the area discontinuity problem near the pole (https://github.com/NOAA-GFDL/FRE-NCtools/issues/44) . There is an existing github issue on this and Niki may have recently addressed through another work that may need to be incorporated in issue 44.

The other issue is that the generation of the supergrids from FV3 internal files are setting the angle_dx and angle_dy to 0 - which seems like a bug (and probably a simple one) and I should write a ticket for it.

[lharris4]

Hi, Miguel. Thanks for running the alternate-precision test, it looks like things are pretty well set other than the two problems you have identified. It is odd that the large percentage error is on tile 1, which is an equatorial tile.

Lucas

On Sat, Nov 19, 2022 at 2:02 PM Miguel R Zuniga @.***> wrote:

@lharris4 https://github.com/lharris4 I think the 32 bit grid data output issue is resolved, but in the course of testing I ran into two different issues in other parts of the code, which I describe way below.

First, in experiment C128r3.solo.TC I changed the grid spec output format from 2 to 1. This changes the resultant diag_table, and the grid spec is printed with 15 digits (plus decimal). With this now all the x,y field values are in tolerance of 10^-9 %. I.e. this check on every file shows no discrepancies:

nccmp -df --Tolerance=1e-9 --variable=x,y grid.tile1.nc ../fv3_sgrids/grid_spec.tile1.nc

And for fields dx,dy, the Tolerance ls 2.0*10^-9%.

However, for grid cell area there is an issue, most values are in Tolerance of 1.0e-4, but for tiles 1 and 3 there is a small number of points with a very large (up to 27%) error. Here is the output for one comparison:

gaea14$ nccmp -df --Tolerance=1e0 --variable=area grid.tile3.nc ../fv3_sgrids/grid_spec.tile3.nc DIFFER : VARIABLE : area : POSITION : [125,201] : VALUES : 6.35747e+09 <> 6.42859e+09 : PERCENT : 1.10627 DIFFER : VARIABLE : area : POSITION : [126,199] : VALUES : 6.63873e+09 <> 6.55313e+09 : PERCENT : 1.28942 DIFFER : VARIABLE : area : POSITION : [126,201] : VALUES : 6.36169e+09 <> 6.50816e+09 : PERCENT : 2.25049 DIFFER : VARIABLE : area : POSITION : [126,202] : VALUES : 6.16146e+09 <> 6.03806e+09 : PERCENT : 2.00288 DIFFER : VARIABLE : area : POSITION : [127,199] : VALUES : 6.63976e+09 <> 6.72003e+09 : PERCENT : 1.19451 DIFFER : VARIABLE : area : POSITION : [127,200] : VALUES : 6.43222e+09 <> 6.18753e+09 : PERCENT : 3.80405 DIFFER : VARIABLE : area : POSITION : [127,201] : VALUES : 6.36265e+09 <> 8.61327e+09 : PERCENT : 26.1296 DIFFER : VARIABLE : area : POSITION : [128,199] : VALUES : 6.63976e+09 <> 6.72003e+09 : PERCENT : 1.19451 DIFFER : VARIABLE : area : POSITION : [128,200] : VALUES : 6.43222e+09 <> 6.18753e+09 : PERCENT : 3.80405 DIFFER : VARIABLE : area : POSITION : [128,201] : VALUES : 6.36265e+09 <> 8.61327e+09 : PERCENT : 26.1296 DIFFER : VARIABLE : area : POSITION : [129,199] : VALUES : 6.63873e+09 <> 6.55313e+09 : PERCENT : 1.28942 DIFFER : VARIABLE : area : POSITION : [129,201] : VALUES : 6.36169e+09 <> 6.50816e+09 : PERCENT : 2.25049 DIFFER : VARIABLE : area : POSITION : [129,202] : VALUES : 6.16146e+09 <> 6.03806e+09 : PERCENT : 2.00288 DIFFER : VARIABLE : area : POSITION : [130,201] : VALUES : 6.35747e+09 <> 6.42859e+09 : PERCENT : 1.10627

This is almost certainly the area discontinuity problem near the pole (#44 https://github.com/NOAA-GFDL/FRE-NCtools/issues/44) . There is an existing github issue on this and Niki may have recently addressed through another work that may need to be incorporated in issue 44.

The other issue is that the generation of the supergrids from FV3 internal files are setting the angle_dx and angle_dy to 0 - which seems like a bug and I should write a ticket for it.

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/FRE-NCtools/issues/179#issuecomment-1320947950, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVALXEQVVYVPXPZVZPLWJEP3JANCNFSM6AAAAAAQHF3ATA . You are receiving this because you were mentioned.Message ID: @.***>

[ngs333]

@lharris4 The PR (https://github.com/NOAA-GFDL/FRE-NCtools/pull/195) for this issue also adds a new unit test, which shows that grid files generated by make_hgrid with the do_transform option are within an nccmp Tolerance of 1.0e-8 (i.e. less than 10^-8 % error) of the FV3 grids for variables x,y,dx,dy and area. This is across a wide range of hardware, with different hardware and compilers for the FV3 runs vs the NCTools runs.