`fv_mapz` is substantially different for GEOS

[FlorianDeconinck]

Per Bill Putman, the mapz code of fvdycore has been substantially reworked and will require a re-port.

A quick glance at code diff yields:

• sat_adj2 is entirely different
• Implement the remap_te option (this can clash with existing code). Including the new map1_cubic
• Implement the remap_t option (this can clash with existing code)
• The MULTI_GASES option is out
• Some moist_cv "cases" are removed in GEOS

Strategy should be to diff the original port (https://github.com/ai2cm/fv3gfs-fortran) with the current code to get a sense. Then use serialbox to generate or re-generate the translate tests.

Note: feel free to break this further into G1.7.X

Parent: GEOS-ESM/SMT-Nebulae#55

---

• Remapping in PyFV3 validates under GEOS translate tests (without breaking the original test case)

[FlorianDeconinck]

2.8.1 brought:

• real(kind=8) :: E_Flux = 0. - mixed precision
• real(kind=8), dimension(is:ie,js:je):: te_2d, zsum0, zsum1 - mixed precision

[gmao-ckung]

Note that the Translate test for moist_cv only tests solutions at j == 1 on the grid. This serialization pattern at j == 1 exists at nearly all serialization points within fv_mapz.F90 in fv3gfs_fortran (ex. moist_cv serialization point)

[gmao-ckung]

Initial Sample of Translate Test results for moist_cv

Failure Results for numpy backend : rank 0

cappa error results

Failed count: 943/64800 (1.46%),
Worst failed index (15, 3, 0)
    Computed:0.2857142686843872
    Reference: 0.2857142388820648
    Absolute diff: 2.980e-08
    Metric diff: 5.714e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 2/943
    
Failed count: 960/64800 (1.48%),
Worst failed index (3, 3, 0)
    Computed:0.28571420907974243
    Reference: 0.28571417927742004
    Absolute diff: 2.980e-08
    Metric diff: 5.960e-08
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/960

Failed count: 956/64800 (1.48%),
Worst failed index (3, 3, 0)
    Computed:0.28571420907974243
    Reference: 0.28571417927742004
    Absolute diff: 2.980e-08
    Metric diff: 5.960e-08
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/956

pt error results

Failed count: 1186/64800 (1.83%),
Worst failed index (17, 3, 61)
    Computed:277.7859802246094
    Reference: 277.78564453125
    Absolute diff: 3.357e-04
    Metric diff: 6.714e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 19/1186

Failed count: 1173/64800 (1.81%),
Worst failed index (12, 3, 63)
    Computed:278.1600646972656
    Reference: 278.15972900390625
    Absolute diff: 3.357e-04
    Metric diff: 6.714e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 18/1173

Failed count: 1148/64800 (1.77%),
Worst failed index (12, 3, 58)
    Computed:276.57177734375
    Reference: 276.5714111328125
    Absolute diff: 3.662e-04
    Metric diff: 7.324e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 19/1148

Failure Results for dace:cpu backend : rank 0

cappa error results

Failed count: 1728/64800 (2.67%),
Worst failed index (6, 3, 0)
    Computed:0.28573555116455995
    Reference: 0.2857142388820648
    Absolute diff: 2.131e-05
    Metric diff: 2.857e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1727/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (6, 3, 0)
    Computed:0.2857355515240836
    Reference: 0.2857142388820648
    Absolute diff: 2.131e-05
    Metric diff: 2.857e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (7, 3, 0)
    Computed:0.2857355546886946
    Reference: 0.2857142388820648
    Absolute diff: 2.132e-05
    Metric diff: 2.857e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1727/1728

pt error results

Failed count: 1728/64800 (2.67%),
Worst failed index (24, 3, 71)
    Computed:295.5802535255386
    Reference: 295.4806213378906
    Absolute diff: 9.963e-02
    Metric diff: 2.955e+02
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (24, 3, 71)
    Computed:295.6366142746722
    Reference: 295.5369873046875
    Absolute diff: 9.963e-02
    Metric diff: 2.955e+02
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (24, 3, 71)
    Computed:295.54787338714056
    Reference: 295.44830322265625
    Absolute diff: 9.957e-02
    Metric diff: 2.954e+02
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

q_con error results

Worst failed index (13, 3, 60)
    Computed:0.00021573402045760664
    Reference: 0.00021573402045760304
    Absolute diff: 3.605e-18
    Metric diff: 2.157e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 401/457

Failed count: 315/64800 (0.49%),
Worst failed index (16, 3, 63)
    Computed:0.0004030582103595748
    Reference: 0.00040305822039954364
    Absolute diff: 1.004e-11
    Metric diff: 4.031e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 235/315

Failed count: 444/64800 (0.69%),
Worst failed index (18, 3, 66)
    Computed:0.00037402475288672576
    Reference: 0.0003740247630048543
    Absolute diff: 1.012e-11
    Metric diff: 3.740e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 364/444

[gmao-ckung]

Notes on Fillz

The GEOS and fv3gfs_fortran version of Fillz are the same, so the currently implemented GT4Py stencil of Fillz should work.

One note is that the current version of NDSL/GT4Py accounts only for 8 9 tracers. Within gt4py_utils.py, there's a 8-entry 9-entry tracer list. GEOS has 31 tracers, and if the translate test is executed with q2tracers serialized with all 31 tracers, the translate test will fail to run since translate_fillz.py only accounts for 8 9 tracers. Limiting q2tracers to serialize only 8 9 tracers in a GEOS run will allow the Fillz translate test to pass with both the numpy and dace:cpu backends.

[FlorianDeconinck]

Initial Sample of Translate Test results for moist_cv

Failure Results for numpy backend : rank 0

cappa error results

Failed count: 943/64800 (1.46%),
Worst failed index (15, 3, 0)
    Computed:0.2857142686843872
    Reference: 0.2857142388820648
    Absolute diff: 2.980e-08
    Metric diff: 5.714e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 2/943
    
Failed count: 960/64800 (1.48%),
Worst failed index (3, 3, 0)
    Computed:0.28571420907974243
    Reference: 0.28571417927742004
    Absolute diff: 2.980e-08
    Metric diff: 5.960e-08
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/960

Failed count: 956/64800 (1.48%),
Worst failed index (3, 3, 0)
    Computed:0.28571420907974243
    Reference: 0.28571417927742004
    Absolute diff: 2.980e-08
    Metric diff: 5.960e-08
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/956

pt error results

Failed count: 1186/64800 (1.83%),
Worst failed index (17, 3, 61)
    Computed:277.7859802246094
    Reference: 277.78564453125
    Absolute diff: 3.357e-04
    Metric diff: 6.714e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 19/1186

Failed count: 1173/64800 (1.81%),
Worst failed index (12, 3, 63)
    Computed:278.1600646972656
    Reference: 278.15972900390625
    Absolute diff: 3.357e-04
    Metric diff: 6.714e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 18/1173

Failed count: 1148/64800 (1.77%),
Worst failed index (12, 3, 58)
    Computed:276.57177734375
    Reference: 276.5714111328125
    Absolute diff: 3.662e-04
    Metric diff: 7.324e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 19/1148

Failure Results for dace:cpu backend : rank 0

cappa error results

Failed count: 1728/64800 (2.67%),
Worst failed index (6, 3, 0)
    Computed:0.28573555116455995
    Reference: 0.2857142388820648
    Absolute diff: 2.131e-05
    Metric diff: 2.857e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1727/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (6, 3, 0)
    Computed:0.2857355515240836
    Reference: 0.2857142388820648
    Absolute diff: 2.131e-05
    Metric diff: 2.857e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (7, 3, 0)
    Computed:0.2857355546886946
    Reference: 0.2857142388820648
    Absolute diff: 2.132e-05
    Metric diff: 2.857e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1727/1728

pt error results

Failed count: 1728/64800 (2.67%),
Worst failed index (24, 3, 71)
    Computed:295.5802535255386
    Reference: 295.4806213378906
    Absolute diff: 9.963e-02
    Metric diff: 2.955e+02
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (24, 3, 71)
    Computed:295.6366142746722
    Reference: 295.5369873046875
    Absolute diff: 9.963e-02
    Metric diff: 2.955e+02
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

Failed count: 1728/64800 (2.67%),
Worst failed index (24, 3, 71)
    Computed:295.54787338714056
    Reference: 295.44830322265625
    Absolute diff: 9.957e-02
    Metric diff: 2.954e+02
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1728/1728

q_con error results

Worst failed index (13, 3, 60)
    Computed:0.00021573402045760664
    Reference: 0.00021573402045760304
    Absolute diff: 3.605e-18
    Metric diff: 2.157e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 401/457

Failed count: 315/64800 (0.49%),
Worst failed index (16, 3, 63)
    Computed:0.0004030582103595748
    Reference: 0.00040305822039954364
    Absolute diff: 1.004e-11
    Metric diff: 4.031e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 235/315

Failed count: 444/64800 (0.69%),
Worst failed index (18, 3, 66)
    Computed:0.00037402475288672576
    Reference: 0.0003740247630048543
    Absolute diff: 1.012e-11
    Metric diff: 3.740e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 364/444

This isn't bad, though some are troubling. Looks like a temperature variable that is at e-2 error, that sounds too high

[gmao-ckung]

Updated moist_cv Translate test results on dace:cpu backend

I eventually noticed that I had not set the environment variable PACE_FLOAT_PRECISION=32 when I ran the dace:cpu backend. I reran the dace:cpu translate tests, and the results are similar or better compared to the numpy backend. Note that these results are based on the TinyBC example (see link on where TinyBC resides).

Failure Results for dace:cpu backend : rank 0

pt errors

Failed count: 890/64800 (1.37%),
Worst failed index (19, 3, 70)
    Computed:280.7839660644531
    Reference: 280.78369140625
    Absolute diff: 2.747e-04
    Metric diff: 5.493e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 16/890

Failed count: 890/64800 (1.37%),
Worst failed index (3, 3, 67)
    Computed:288.0103759765625
    Reference: 288.0101013183594
    Absolute diff: 2.747e-04 
    Metric diff: 5.493e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 15/890

Failed count: 879/64800 (1.36%),
Worst failed index (16, 3, 63)
    Computed:276.6903076171875
    Reference: 276.6900329589844
    Absolute diff: 2.747e-04 
    Metric diff: 5.493e-04
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 16/879

cappa errors

Worst failed index (3, 3, 0)
    Computed:0.28571420907974243
    Reference: 0.28571417927742004
    Absolute diff: 2.980e-08 
    Metric diff: 5.960e-08
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/956

Failed count: 960/64800 (1.48%),
Worst failed index (3, 3, 0)
    Computed:0.28571420907974243
    Reference: 0.28571417927742004
    Absolute diff: 2.980e-08 
    Metric diff: 5.960e-08
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/960

Failed count: 943/64800 (1.46%),
Worst failed index (15, 3, 0)
    Computed:0.2857142686843872
    Reference: 0.2857142388820648
    Absolute diff: 2.980e-08 
    Metric diff: 2.857e-01
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 2/943

A final note is that q_con passed the translate tests with both the numpy and dace:cpu backends.

[FlorianDeconinck]

Nice, temperature is much better I think we can move forward

[gmao-ckung]

Additional Note on Fillz

I previously mentioned that gt4py_utils.py currently accounts for 9 tracers in tracer_variables. If I extend tracer_variables to have 31 entries and reran the Fillz translate test with 31 tracers serialized to q2tracers from a GEOS run, the Fillz translate tests pass using both the numpy and dace:cpu backends.

In this case, I wrote tracer_variables as follows.

tracer_variables = [
    "qvapor",
    "qliquid",
    "qrain",
    "qice",
    "qsnow",
    "qgraupel",
    "qo3mr",
    "qsgs_tke",
    "qcld",
    "tracer10",
    "tracer11",
    "tracer12",
    "tracer13",
    "tracer14",
    "tracer15",
    "tracer16",
    "tracer17",
    "tracer18",
    "tracer19",
    "tracer20",
    "tracer21",
    "tracer22",
    "tracer23",
    "tracer24",
    "tracer25",
    "tracer26",
    "tracer27",
    "tracer28",
    "tracer29",
    "tracer30",
    "tracer31",
]

For now, I'll leave tracer_variables with its original 9 entries.

[gmao-ckung]

Updating map_scalar implementation in pyFV3

The next significantly different code difference between the GEOS and GFDL version of fv_mapz.F90 exists in the map_scalar subroutine. In the GEOS version of map_scalar, the subroutine performs an additional field interpolation onto target pressures, an interpolation that does not exist in the GFDL version of map_scalar. This interpolation needs to be incorporated into pyFV3.

The portion of map_scalar that we're interested in modifying exists in pyFV3 within the lagrangian_contributions stencil in map_single.py.

GEO's map_scalar performs the field interpolation by comparing two pressure values, pe1 and pe2, and this comparison results in two k indices LM1 and LP0 from where the interpolation will take place. The tricky part of finding these two indices is that at each kth level, the search for LP0 always starts at k = 1. Given that GT4Py values are indexed relatively (though a feature to have absolute k indexes available to Quantities is in the works), there needs to be a way to get the integer value that stores the "k = 1" level (note this is from the Fortran perspective) at each k. This can be performed by creating an Integer quantity that basically stores the relative k = 1 index for each level of k. With that, I added a computation that takes two Integer quantities, INDEX_LM1 and INDEX_LP0, and makes them store the relative index for k = 1 at each k level. With these two quantities, I can essentially do the same computation based on the GEOS map_scalar code to find the indices needed for the interpolation.

The original AI2 port of lagrangian_contributions does not implement the if statement within map_scalar. I included this if statement in the GEOS implementation of lagrangian_contributions, which updates the lev variable that controls the current k level based on the variable l in the Fortran code. Adding the if statement results in lagrangian_contributions passing the C24 Translate tests with very minimal errors. The following example message shows the error severity when using the numpy backend.

_ test_sequential_savepoint[Lagrangian_Contribution-rank=0-call=2] (varname='q1') _
❌ Numerical failures
Failed count: 1/64800 (0.0%),
Worst failed index (3, 3, 48)
    Computed:260.9957580566406
    Reference: 260.9957275390625
    Absolute diff: 3.052e-05
    Metric diff: 5.220e+02
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/1
All failures:
Index  Computed  Reference  Absloute E  Metric E
(3, 3, 48)  260.9957580566406  260.9957275390625  3.052e-05  5.220e+02

Note that the translate tests do not execute with the dace:cpu or gt:cpu_kfirst backends. Both these backends seems to report errors during the code generation. These backends seem to be tripped up by the use of Integer quantities INDEX_LM1 and INDEX_LP0 in the interpolation calculation.

Update (11/13/2024): It turns out that GEOS will also use lagrangian_contributions calls that do not contain the interpolation mentioned above. Give this situation, I've updated the MapSingle class call with a boolean variable interp (which defaults to False) that gives the developer a choice to use the interpolation or not. Within MapSingle, there are two separate lagrangian_contributions stencils implemented: one with the interpolation (lagrangian_contribtuions_interp) and one without the interpolation.

[gmao-ckung]

NDSL software stack 2024.09.00 presenting translate test runtime issues with map_scalar / lagragian_contribution

When the translate test for lagrangian_contribution uses NDSL version 2024.09.00, the translate test seems to be stuck within a while loop in map_single.py.

Update (11/13/2024): Roman and Tobias fixed the numpy backend that caused the while loop issue. There is an ongoing issue with the while loop in the debug backend that Tobias is investigating.

[gmao-ckung]

Notes comparing GEOS v11.5.2 used for SMT development versus upcoming GEOS v12 regarding fv_mapz

Github can present a visual comparison between the geos/dsl/develop branch and the feature/wmputman/LM_v12rc3_WMP branch of GFDL_atmos_cubed_sphere in the following link: Code Comparison Link. Click the "Files changed" tab to view the file differences between the two branches.

In general, the two code bases structurally look the same. There's consolidation of several routines such as map1_ppm and mapn_tracers into the map_scalar subroutine. I also looked at subroutines remap_z, scalar_profile, and cs_profile, and these subroutines have minor changes but no major algorithmic changes.

[gmao-ckung]

Translate test for pressure updates in latter j loop

I created a stencil pe_pk_delp_peln that computes the pressure updates in the end of the latter j loop in Lagrangian_to_Eulerian. The corresponding translate test has also been written and tested.

The translate tests pass for the dace:cpu and gt:cpu_kfirst backends. For the numpy and debug backends, there are solution errors that are about 7 to 8 orders of magnitude smaller than the solution order, and there are a handful of unique errors.

__ test_sequential_savepoint[PE_pk_delp_peln-rank=0-call=0] (varname='peln_') __
❌ Numerical failures
Failed count: 30/42048 (0.07%),
Worst failed index (4, 0, 42)
    Computed:9.943809509277344
    Reference: 9.943808555603027
    Absolute diff: 9.537e-07
    Metric diff: 9.591e-08
    Metric threshold: 1e-14
  Noise quantification:
    Reference dtype: <class 'numpy.float32'>
    Unique errors: 1/30
All failures:
Index  Computed  Reference  Absloute E  Metric E
(0, 0, 55)  11.167181015014648  11.167180061340332  9.537e-07  8.540e-08
(1, 0, 43)  10.10767936706543  10.107678413391113  9.537e-07  9.435e-08
(1, 0, 50)  10.878594398498535  10.878593444824219  9.537e-07  8.767e-08
(2, 0, 60)  11.326676368713379  11.326675415039062  9.537e-07  8.420e-08
(3, 0, 49)  10.80120849609375  10.801207542419434  9.537e-07  8.829e-08
(3, 0, 69)  11.477877616882324  11.477876663208008  9.537e-07  8.309e-08
...

With the limited repeatable errors that are several orders of magnitude lower than the solution, I declare that the numpy and debug backends have "verified".

[gmao-ckung]

Notes on mpp_global_sum call in fv_mapz.F90

For our purposes, mpp_global_sum is an MPI_ALLREDUCE collective call implemented in the FMS framework. The tracing of this particular call goes through several layers of include files and macro definitions, and the following outline attempts to sketch that call path from fv_mapz.F90.

• mpp_global_sum is called (perhaps multiple times) within the "last step" of the Lagrangian_to_Eulerian subroutine call.
• Implementation-wise, mpp_global_sum itself is a Fortran interface to multiple module procedure names found in mpp_domains.F90.
• The procedures that fall under mpp_global_sum are referenced via a #define macro called MPP_GLOBAL_SUM_ within mpp_domains_reduce.inc. I suspect that the module procedure mpp_global_sum_r8_2d is used, but I did not confirm this.
• MPP_GLOBAL_SUM_ is a "prototype" function implemented in mpp_global_sum.h. For our purposes, the optional parameter flags for MPP_GLOBAL_SUM_ is set to BITWISE_EFP_SUM, a parameter from the mpp_parameter_mod module. This flag leads MPP_GLOBAL_SUM_ down a particular code path to call mpp_reproducing_sum.
• mpp_reproducing_sum is a Fortran interface to multiple module procedure names found in mpp_efp.F90. For our purposes, mpp_reproducing_sum maps into mpp_reproducing_sum_r8_2d.
• mpp_reproducing_sum_r8_2d calls mpp_sum, which is a Fortran interface to multiple module procedure names in mpp.F90. These procedure names are mapped to the macro MPP_SUM_ in mpp_comm_mpi.inc.
• MPP_SUM_ is defined in mpp_sum_mpi.h, and its implementation contains an MPI_ALLREDUCE call.

The code path, while convoluted, is not particularly difficult to implement. The tricky part will be to make sure that we're covering all the mpp- based module parameters and variables referenced from the domain2D type used in fv_mapz.F90