Switch wind stress to use isotropic interpolation from cells to edge normals

[xylar]

This merge adds a namelist option that defaults to the current anisotropic approach to interpolating wind stresses in MPAS-Ocean but which can be set to use a new isotropic interpolation approach.

This merge also adds 2 new subroutines to the framework that interpolate a 1d field of 2d (zonal and meridional) vectors from cell centers to edges.

[xylar]

This passes the SMS_D_Ld5.T62_oQU240wLI.GMPAS-IAF.chrysalis_intel test. The next thing to do would be to add a stealth test that uses the isotropic interpolation.

[proteanplanet]

I have started testing and will report back here on results from a G-case comparison.

[xylar]

@proteanplanet, I haven't tested the isotropic interpolation at all yet so maybe hold off until I have a chance to run a similar test to the one above but with config_bulk_wind_stress_interp_isotropic = .true.. Right now, E3SM won't even know about that config option so even setting a user namelist option won't work, I don't think.

[xylar]

@proteanplanet, sorry for the confusion last night. I now have this tested in both "anisotropic" and "isotropic" modes with a simple smoke test:

SMS_D_Ld5.T62_oQU240wLI.GMPAS-IAF.chrysalis_intel.mpaso-isotropic_wind_stress
SMS_D_Ld5.T62_oQU240wLI.GMPAS-IAF.chrysalis_intel

I believe it is ready for you to try out. We should plot the resulting normal component of wind stress at edges as early in the simulation as possible to avoid wasting a lot of compute time if there are obvious bugs.

[proteanplanet]

Testing so far, using this branch merged with the current E3SM master in https://github.com/proteanplanet/E3SM/tree/proteanplanet/ocn/isowindstress:

This passes the e3sm_cryo_developer suite on Chrysalis:

     "SMS_D_Ld1.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF.mpaso-jra_1958",
     "ERS_Ld5.T62_oQU240wLI.GMPAS-DIB-IAF-PISMF",
     "PEM_Ln5.T62_oQU240wLI.GMPAS-DIB-IAF-PISMF",
     "PET_Ln5.T62_oQU240wLI.GMPAS-DIB-IAF-PISMF",
     "ERS_Ld5.T62_oQU240wLI.GMPAS-DIB-IAF-DISMF",
     "PEM_Ln5.T62_oQU240wLI.GMPAS-DIB-IAF-DISMF",
     "PET_Ln5.T62_oQU240wLI.GMPAS-DIB-IAF-DISMF",

This hangs on Chrysalis upon initialization, but runs on Anvil:

COMPSET[G]="GMPAS-JRA1p5"
RESOLUTION[G]="TL319_EC30to60E2r2"
PELAYOUT[G]="L"

Analysis of Anvil results will continue while the hanging problem on Chrysalis is investigated.

[jonbob]

I can verify that a G-case with the default layout hangs. The chrysalis default layout is threaded but has the ocean running on its own set of processors. As @proteanplanet reported, the run simply hangs and the resulting traceback doesn't offer much information. But it does not appear to be a threading issue, since both of these tests pass:
SMS_D_P480x1.TL319_EC30to60E2r2.GMPAS-JRA1p5
SMS_D_P480x2.TL319_EC30to60E2r2.GMPAS-JRA1p5

[proteanplanet]

Here are some early results from G-case testing with GMPAS-JRA1p5 and TL319_EC30to60E2r2:

[cbegeman]

@proteanplanet Interesting. I would have expected new-old kineticEnergyCellSum to be positive if more momentum was being transferred into the ocean with the more accurate interpolation, but it looks negative in your plot. Maybe the flow is more oscillatory in the old case due to more interpolation error? The idealized tests we've planned seem useful for understanding this.

[jonbob]

An update on the strange testing hang. It appears to be a seaice partition file that's causing the problem, specifically:

/lcrc/group/e3sm/data/inputdata/ice/mpas-seaice/EC30to60E2r2/partitions/mpas-seaice.graph.info.230313.part.160

The layout that has illustrated the problem is the default one, with all components using 2 threads and atm, ice, cpl, lnd and rof on 160 pes and the ocn on its own 320 pes. Changing the atm/ice/cpl/lnd/rof count to 128 pes runs fine, as does 192. When I tested with 1 thread, everything worked until the ice ran on 160 pes as well. It also ran fine using the older partition file

mpas-seaice.graph.info.200908.part.160

I checked a semi-recent version of master and it behaved similarly, so it does not appear to be anything introduced in this PR.

[jonbob]

I checked and that graph file does send points to each processor. And it does have the correct number of points

[xylar]

@proteanplanet, the main thing I'm finding in idealized testing is that the new wind-stress reconstruction results in about twice as much stress at edges at land boundaries than the previous version:

anisotropic:

isotropic:

This is because the anisotropic method normalizes the stress at an edge by the total kite area even if there is only 1 adjacent cell, whereas the new, isotropic approach normalizes by the total area associated with cells that exist. I can't see any way to normalize the isotropic approach by the total area associated with both ocean and land cells because we have removed the land cells and their associated kite areas.

[xylar]

By default, the windstress in the idealized test case I'm using (MPAS-Dev/compass#547) varies quite smoothly. If I increase the frequency by a factor of 128 (so it varies at close to the grid scale), I see:

zonal wind stress at cell centers (either approach):

anisotropic interpolation to edges:

isotropic interpolation to edges:

It is interesting that the isotropic case produces stronger interpolated values for zonally-oriented edges in this case. I had expected smoother results but it wasn't clear that that would mean higher zonal values.

[xylar]

I'm going to try a case that's a little less extreme.

[xylar]

the main thing I'm finding in idealized testing is that the new wind-stress reconstruction results in about twice as much stress at edges at land boundaries than the previous version

I believe this is a moot point. The normal velocity at edges adjacent to land is constrained to be zero (no flow into or out of land) so whatever wind stress gets computed at these edges won't be applied anyway.

[xylar]

Same as above but for a factor of 32 higher frequency than the usual wind stress (rather than 128):

The zonal component at cell centers:

The normal component with anisotropic interpolation:

The normal component with isotropic interpolation:

Analytic solution coming soon...

[xylar]

Here's the comparison between the respective interpolation methods and the analytic solution computed directly at edges for the same wind-stress field plotted in my previous message.

anisotropic - analytic:

isotropic - analytic:

[xylar]

Ignoring land-adjacent edges, where the surface stress will be ignored, the error is much smaller with the anisotropic method, presumably because of the extra smoothing introduced by the isotropic method. But it seems like the error is on the order of a few percent for the isotropic case so I don't have any concerns that it isn't implemented correctly.

[cbegeman]

@xylar This is fantastic! Thanks so much for doing that verification test.