Fix bugs in ESMF interpolation method #544

matthewhoffman · 2023-12-17T06:34:24Z

This merge fixes a number of issues with the 'esmf' method for interpolation in the MALI interpolation script:

indexing was off by one, which made interpolation with unstructured source meshes garbage. For structured source meshes, it made interpolation shifted by one grid cell.
support MPAS source fields with a vertical dimension when using the 'esmf' method
refactor to use sparse matrix multiply, which speeds up interpolation a few hundred times
add destination mesh area normalization support. This is necessary when using the ESMF 'conserve' method for destination cells that are only partly overlapped by the source cells.

Two issues remain for interpolating between two MPAS meshes with the 'esmf' method:

If the destination mesh is larger than the source mesh, those locations are filled with zeros. The ESMF 'conserve' method does not support extrapolation, and there is no obvious solution to this issue and it would need to be handled manually on a case by case basis.
Some fields are only defined on subdomains (e.g. temperature is only defined where ice thickness is nonzero) and the script currently has no mechanism for masking them. This results in garbage values getting interpolated in, e.g. temperature values near the margin will have values around 100 K because of interpolating realistic values around 250K with garbage values of 0K. This issue applies to the barycentric method as well. However, this situation can be worked around by performing extrapolation of the temperature field on the source mesh before doing interpolation between meshes.

This PR also includes a refactoring of create_SCRIP_file_from_planar_rectangular_grid.py that speeds it up by orders of magnitude for large meshes, as well as a minor update to define_cullMask.py

This bug yielded garbage results when interpolating between two MPAS meshes. When interpolating from a regular grid to an MPAS grid, it would cause a shift of one cell on the regular grid.

For interpolating from 4km AIS (385379 cells) to 8km AIS (98341 cells) this reduced interpolation time for one variable from 13 seconds to 0.04 seconds (~300x speedup).

This was never attempted before and needed additional handling.

See section 12.5 of ESMF_RegridWeightGen documentation: https://earthsystemmodeling.org/docs/release/ESMF_8_0_1/ESMF_refdoc/node3.html#SECTION03025000000000000000

matthewhoffman · 2023-12-18T22:48:27Z

@trhille , I decided to leave this branch as is and ask for review now (no rush). The issue about the conserve method filling 0s in regions of extrapolation of the source mesh is a larger problem that it is not obvious how to resolve, so I left it as a known issue for now. But if you have ideas about appropriate functionality for that situation, I'm happy to do further work in this PR.

matthewhoffman · 2023-12-21T07:04:10Z

@xylar , there is a commit in this PR that makes an update to mesh_tools/create_SCRIP_files/create_SCRIP_file_from_planar_rectangular_grid.py. However, it is unclear to me if I should also be updating the copy of this file at conda_package/mesh_tools/create_SCRIP_files/create_SCRIP_file_from_planar_rectangular_grid.py. It seems that I need to do that for compass to be able to see my changes, but when I try to commit that change, I get this error:

The following paths are ignored by one of your .gitignore files:
conda_package/mesh_tools

Can you explain to me how this is supposed to be handled? I've never fully understood how the conda package interacts (or doesn't) with the legacy scripts.

xylar · 2023-12-21T20:38:40Z

@matthewhoffman, when you do a local install of the conda_package into a conda environment, e.g. with:

python -m pip install conda_package

it makes a copy of the mesh_tools, ocean, land_ice and other directories that contain scripts because it can't include anything in the python package that lives outside of conda_package. So if you're modifying scripts outside of conda_package, you should run the command above again to first copy the script(s) over again and then install them into the conda environment.

Ideally, if we want to have mpas_tools be a "proper" python package, we would stop making scripts and switch to entry points that live within conda_package/mpas_tools. But given that that's not the workflow that you all are used to, this copying over is necessary.

matthewhoffman · 2023-12-21T21:53:41Z

@xylar , thanks for clarifying that. I had run python -m pip install conda_package but I hadn't appreciated that because of the legacy scripts, that that command would need to be rerun every time one of those scripts was modified. Now I understand the workflow and have more incentive to move things into the conda package. :)

For large grids, this script was very slow. For the 500m BedMachine-AIS dataset, it's been running for over half an hour and is not complete. After refactoring, it took about 3 minutes.

matthewhoffman · 2024-01-04T20:08:51Z

Testing

These updated scripts were used for the 2km AIS mesh generation in MPAS-Dev/compass#750 and it worked as expected end to end. The refactoring sped up those individual operations dramatically.

trhille

@matthewhoffman, this is fantastic! Thanks for making these changes, which collectively reduce the pain of mesh generation by a huge margin. I have one measly complaint about an inscrutable variable name, but otherwise this is ready to go.

trhille · 2024-01-04T23:13:42Z

landice/mesh_tools_li/interpolate_to_mpasli_grid.py

    wfile.close()
    #----------------------------

+    # convert to sparse matrix format
+    wt_csr = scipy.sparse.coo_array((S, (row - 1, col - 1)), shape=(n_b, n_a)).tocsr()


I dislike the variable name wt_csr, because it requires googling to figure out what it means. Can we use a longer, more descriptive name? Also, is wt short for weight?

Updated in b0e6147.

trhille · 2024-01-04T23:14:32Z

landice/mesh_tools_li/interpolate_to_mpasli_grid.py

-    sourceFieldFlat = sourceField.flatten()  # Flatten source field
-    for i in range(len(row)):
-      destinationField[row[i]-1] = destinationField[row[i]-1] + S[i] * sourceFieldFlat[col[i]]
+    source_csr = scipy.sparse.csr_matrix(sourceField.flatten()[:, np.newaxis])


Same complaint as above about the csr in the variable name, although here it's adjacent to the scipy method that helps explain it some.

Updated in b0e6147.

trhille · 2024-01-04T23:16:52Z

landice/mesh_tools_li/interpolate_to_mpasli_grid.py

@@ -328,7 +337,7 @@ def interpolate_field_with_layers(MPASfieldName):
        if filetype=='cism':
           print('  Input layer {}, layer {} min/max: {} {}'.format(z, InputFieldName, InputField[z,:,:].min(), InputField[z,:,:].max()))
        elif filetype=='mpas':
-           print('  Input layer {}, layer {} min/max: {} {}'.format(z, InputFieldName, InputField[:,z].min(), InputField[z,:].max()))


Was this just a typo?

Yeah, this was a typo that has been present from whenever this print statement was introduced. It wasn't obviously incorrect from the output because the typo version still yielded reasonable values.

xylar · 2024-01-08T15:08:38Z

I'm about to make a new release to fix #547. If this can be merged today, I'll wait for it. Otherwise, it will be in the next release.

xylar · 2024-01-08T15:09:53Z

If this can be merged today

On second thought, I'll be testing in E3SM-Unified for a few more days most likely so if this gets merged in the next 2 or 3 days, that should be soon enough.

matthewhoffman · 2024-01-18T03:17:22Z

@trhille , thanks for the feedback on these changes. I've added a little more detail to the comments and improved the new variable name. If you prefer more detail, let me know.

trhille · 2024-01-18T15:43:22Z

@matthewhoffman, thanks! That looks great. This is ready to merge, as far as I'm concerned, although I see that two of the CI checks were not successful.

matthewhoffman · 2024-01-18T18:45:35Z

@trhille , thanks for catching that. Not sure why they didn't run, but I manually retriggered them and they all pass now. @xylar , did you have any other feedback on this given it touches one non-MALI specific tool?

xylar

The SCRIP changes look good to me!

xylar · 2024-01-19T10:05:35Z

I need to make an MPAS-Tools release ASAP so let me know if this can get merged today.

matthewhoffman added 4 commits December 16, 2023 20:34

Fix bug in indexing when using esmf interplation

880e35f

This bug yielded garbage results when interpolating between two MPAS meshes. When interpolating from a regular grid to an MPAS grid, it would cause a shift of one cell on the regular grid.

Refactor loop to use scipy sparse matrix multiply

4cbe511

For interpolating from 4km AIS (385379 cells) to 8km AIS (98341 cells) this reduced interpolation time for one variable from 13 seconds to 0.04 seconds (~300x speedup).

Add support for ESMF interp of MPAS fields with a vertical dimension

625e505

This was never attempted before and needed additional handling.

Add dst_area normalization

3cd10a2

See section 12.5 of ESMF_RegridWeightGen documentation: https://earthsystemmodeling.org/docs/release/ESMF_8_0_1/ESMF_refdoc/node3.html#SECTION03025000000000000000

matthewhoffman added bug enhancement don't merge clean up landice and removed don't merge labels Dec 17, 2023

matthewhoffman requested a review from trhille December 18, 2023 22:46

matthewhoffman added the don't merge label Dec 21, 2023

matthewhoffman added 2 commits December 21, 2023 21:11

Refactor create_SCRIP_file_from_planar_rectangular_grid.py

abdc035

For large grids, this script was very slow. For the 500m BedMachine-AIS dataset, it's been running for over half an hour and is not complete. After refactoring, it took about 3 minutes.

Fix print statements in landice cull mask script

f758259

matthewhoffman force-pushed the landice/esmf_interp_fix branch from 5b9e298 to f758259 Compare December 22, 2023 04:12

matthewhoffman removed the don't merge label Dec 22, 2023

matthewhoffman mentioned this pull request Dec 22, 2023

Updated data interpolation for Antarctica mesh generation case MPAS-Dev/compass#750

Merged

5 tasks

trhille approved these changes Jan 4, 2024

View reviewed changes

Update comments and variable name to be more descriptive

b0e6147

matthewhoffman self-assigned this Jan 18, 2024

xylar approved these changes Jan 19, 2024

View reviewed changes

matthewhoffman merged commit 6c89e44 into MPAS-Dev:master Jan 19, 2024
7 checks passed

matthewhoffman deleted the landice/esmf_interp_fix branch January 19, 2024 16:23

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix bugs in ESMF interpolation method #544

Fix bugs in ESMF interpolation method #544

matthewhoffman commented Dec 17, 2023 •

edited

Loading

matthewhoffman commented Dec 18, 2023

matthewhoffman commented Dec 21, 2023

xylar commented Dec 21, 2023 •

edited

Loading

matthewhoffman commented Dec 21, 2023

matthewhoffman commented Jan 4, 2024

trhille left a comment

trhille Jan 4, 2024

matthewhoffman Jan 18, 2024

trhille Jan 4, 2024

matthewhoffman Jan 18, 2024

trhille Jan 4, 2024

matthewhoffman Jan 18, 2024

xylar commented Jan 8, 2024

xylar commented Jan 8, 2024

matthewhoffman commented Jan 18, 2024

trhille commented Jan 18, 2024

matthewhoffman commented Jan 18, 2024

xylar left a comment

xylar commented Jan 19, 2024

Fix bugs in ESMF interpolation method #544

Fix bugs in ESMF interpolation method #544

Conversation

matthewhoffman commented Dec 17, 2023 • edited Loading

matthewhoffman commented Dec 18, 2023

matthewhoffman commented Dec 21, 2023

xylar commented Dec 21, 2023 • edited Loading

matthewhoffman commented Dec 21, 2023

matthewhoffman commented Jan 4, 2024

Testing

trhille left a comment

Choose a reason for hiding this comment

trhille Jan 4, 2024

Choose a reason for hiding this comment

matthewhoffman Jan 18, 2024

Choose a reason for hiding this comment

trhille Jan 4, 2024

Choose a reason for hiding this comment

matthewhoffman Jan 18, 2024

Choose a reason for hiding this comment

trhille Jan 4, 2024

Choose a reason for hiding this comment

matthewhoffman Jan 18, 2024

Choose a reason for hiding this comment

xylar commented Jan 8, 2024

xylar commented Jan 8, 2024

matthewhoffman commented Jan 18, 2024

trhille commented Jan 18, 2024

matthewhoffman commented Jan 18, 2024

xylar left a comment

Choose a reason for hiding this comment

xylar commented Jan 19, 2024

matthewhoffman commented Dec 17, 2023 •

edited

Loading

xylar commented Dec 21, 2023 •

edited

Loading