Add in-line self attraction and loading for global tidal simulation

[sbrus89]

Self attraction and loading (SAL) forcing involves a spherical harmonic transform (SHT) of the SSH (Ray, 1998). The spherical harmonic coefficients are scaled by the load Love numbers, then an inverse transformation is performed to compute the SAL field. The SAL forcing is incorporated in the pressure gradient term along with the tidal potential.

Two different approaches are implemented for performing the SHT:

1. SSH is gathered onto a single MPI rank and interpolated to a Gaussian mesh. A very fast library, SHTns, is used to compute the forward/inverse SHT. The SAL field is then interpolated back to the unstructured mesh and scattered back to the other MPI ranks. This work was done by @knbarton.
2. The SHT integration is performed on each subdomain and the SH coefficients are globally summed with an all-reduce.

Approach 1 ) performs better at low core-counts, but 2) scales better for larger runs.

This PR includes an init mode option for verification testing of the forward/inverse SHT.

This work was done as a part of the tides task of the ICoM project.

[NML]
[BFB] (for all standard E3SM tests, and for tests with SAL on when config_parallel_self_attraction_loading_bfb = .true.)

[sbrus89]

Here are some verification results for both SHT algorithms:

Approximation errors with different SH orders are shown for two different analytical functions. The first is on a quasi uniform resolution mesh with resolutions: QU60 (165k cells), QU30 (659k cells), and QU15 (2.6M cells). In all cases the SH order for the serial approach is N = Nlat/2 - 1. Where Nlat is the number of latitudes in the Gaussian grid.

The second analytical function is on an "equatorial refined" mesh with resolutions: ER120to60 (60k cells), ER12to30 (118k cells), and ER12to15 (286k cells). The oscillations vary according to mesh resolution such that there are approximately 30 grid cells per period.

Convergence results for the QU meshes are:

The convergence of the serial approach is slower due to interpolation errors. Convergence for the parallel approach is limited by discretization error due to mid-point rule integration.

And convergence results for the ER meshes are:

Looking across equivalent accuracy levels, the performance of the two SHT algorithms for the QU meshes are:

QU60

QU30

QU15

The performance comparisons for the ER meshes are:

ER120to60

ER120to30

ER120to15

All runs were done on Compy with Intel compilers.

[knbarton]

I performed harmonic analysis comparison to benchmark TPXO8 tidal data from a single-layer barotropic run on Icosahedron 9 (16km) uniform mesh. This shows the serial version of SAL in this pull request gives the same tidal results as the serial version built previously. I have also attached a comparison with the scalar SAL approximation to show the improvement due to full SAL calculation.

This was build on Cori with these settings : make gnu-nersc USE_PIO2=true OPENMP=false DEBUG=false GEN_F90=true

---

M2 Constituent

M2 Complex RMSE Values in Deep water (<1000m) between 66N and 66S :
Serial before : 0.06399980070877255m
Serial currently: 0.06399973949098071m
Scalar approximation: 0.072925093765029m

K1 Constituent

K1 Complex RMSE Values in Deep water (<1000m) between 66N and 66S :
Serial before : 0.02360560164062904m
Serial currently: 0.023605630922693552m
Scalar approximation: 0.02605653002810445m

N2 Constituent

N2 Complex RMSE Values in Deep water (<1000m) between 66N and 66S :
Serial before : 0.01316337401656041m
Serial currently: 0.013163374004691027m
Scalar approximation: 0.014364836226592457m

O1 Constituent

O1 Complex RMSE Values in Deep water (<1000m) between 66N and 66S :
Serial before : 0.014989959013047288m
Serial currently: 0.014989994509038505m
Scalar approximation: 0.021686930435006916m

S2 Constituent

S2 Complex RMSE Values in Deep water (<1000m) between 66N and 66S :
Serial before : 0.02929243376361691m
Serial currently: 0.02929240084208756m
Scalar approximation: 0.034033753955931854m

---

We do expect a slight difference in results between the previous serial SAL calculation and the one in this PR due to changing from N = Nlat - 1 to N = Nlat/2 - 1.

Currently, there is an error when using the parallel SAL. When that is resolved, I will share those harmonic analysis results as well.

[knbarton]

Here are the results for the parallel SAL implementation on Icosahedron 9 mesh using the same settings as before.

RMSE Values in Deep water (<1000m) between 66N and 66S in meters :
M2 Constituent: 0.06416788614061514
K1 Constituent: 0.023685223640391814
N2 Constituent: 0.013199779825957164
O1 Constituent: 0.015113880782811768
S2 Constituent: 0.02942663706041666

M2 Constituent plots:

[mark-petersen]

@sbrus89 thanks. Please check my last commit d6daa07. Stand-alone was not compiling with SHTNS=false because it referenced some shtns variables. I just put those lines within the cpp constructs, so should be fine for the parallel version.

[mark-petersen]

Note, this parenthesis correction was added to this PR because config_use_implicit_bottom_roughness is only used in the tidal simulations used in this PR.

[mark-petersen]

@sbrus89 we typically place OMP pragmas on the outer loop, e.g.

               !$omp parallel
               !$omp do schedule(runtime)
               do iCell = 1, nCellsAll
                   ...
               end do
               !$omp end do
               !$omp end parallel

with private clauses for internal indices as needed. That said, I think the $OMP would go on the outermost loop m here, if m iterations are independent, which they appear to be.

OpenMP is particularly important for the parallel implementation, but could also be used for serial now that we found that bash error.

[mark-petersen]

BTW, OMP is not needed for init mode. It is much less important for init subroutines in forward mode.

[mark-petersen]

@sbrus89 and/or @knbarton could you, with parallel and serial with full SAL on, conduct these three tests after the OpenMP commit, in optimized mode:

1. restart test: run (say) a day, restart an extra day, and separately two days, compare restarts at end
2. partition test: run identical short test with two partitions (say 16 vs 32 cores for a day)
3. thread test: run with OPENMP=false, run for a day, recompile with OPENMP=true, set OMP_NUM_THREADS=2, run for day.

for all three, use a bfb comparison of last restart file, with python or ncdiff. Simple ncdiff commands:

ncdiff file1.nc file2.nc diff.nc
ncwa -O -y mabs -v normalVelocity,layerThickness diff.nc diffmabs.n

max absolute difference should be exactly zero. You can use compass for this, but if you already have SAL runs set up, might be easier to just copy your settings from there.

[mark-petersen]

Merged into master locally and retested. Passes

SMS_D_Ln9.T62_oQU120_ais20.MPAS_LISIO_TEST.chrysalis_gnu
SMS_D_Ln9.T62_oQU120_ais20.MPAS_LISIO_TEST.chrysalis_intel
PET_Ln3.T62_oEC60to30v3wLI.GMPAS-DIB-IAF-ISMF.chrysalis_gnu

Also compiled with gnu and intel on grizzly with SHTNS off (default). Passes pr test suite for both, with bfb match to compass reference. Compiled gnu with SHTNS=true but off using config flags. That also passes pr suite with both debug and optimized gnu.

[sbrus89]

Here are the results I got for BFB testing for the spherical harmonic transform (SHT) only (via init mode, config_init_configuration = 'test_sht') :

Compy: Intel 19.0.5, DEBUG=false

	serial max diff	parallel max diff
1 - 2 nodes	BFB	1.55431223447522e-15
2 - 4 nodes	BFB	1.77635683940025e-15
unthreaded - 2 threads (2 nodes)	(compile error for openMP)	3.5527136788005e-15
unthreaded - 4 threads (2 nodes)	(compile error for openMP)	3.99680288865056e-15

The non-BFB behavior for the parallel SHT is due to the reproducibility of the MPI_AllReduce call and !$omp atomic statements.

For a forward mode run on the icos9 mesh, the parallel version is BFB for restarts, but not across partitions due to the differences in the SHT above. The max difference after two days between a 32 node and 64 node run was 8.83346729096957e-11 for layerThickness and 1.7656587997239e-11 for normalVelocity. The serial version is BFB for restarts and across partitions for this case.

[mark-petersen]

Approving, based on lengthy review and testing by @sbrus89 @knbarton and myself described above. Only remaining issue is that the parallel version of the SHT is non-bfb across partitions due to operation order of the global sums (see previous comment), but it is machine precision. I think that result is sufficient for merging in the PR.

[jonbob]

@philipwjones is working on implementing a reproducible sum in MPAS, which would solve the problem of non-BFB results on different processor counts.

[jonbob]

@mark-petersen - I re-requested your review on this PR because it would be more helpful once the reproducible sum has been added. As-is, it really isn't ready to be merged

[sbrus89]

@jonbob and @mark-petersen: Now that #4700 has been submitted, I'll work on testing this with the reproducible sum implemented by @philipwjones .

[mark-petersen]

I retested with today's commits using gnu debug on badger with SAL and tidal potential on, and bfb flag first on. The full nightly suite passes.

@jonbob this is ready to merge when the repo is open.

[jonbob]

@philipwjones and @matthewhoffman - can you please review after the changes have been made? Thanks

[hollyhan]

This comment is to confirm that the issue I experienced in my previous comment (#4472 (comment)) has been addressed by @knbarton's commit (badff01).

Previously, I experienced a bound-check error during coupled MALI - Sea Level Model simulations in which the subroutine interpolation adopted from the code components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F is utilized. The error only appeared in the DEBUG mode. The coupled MALI-SLM simulations were terminated with a runtime error with the following message: "Fortran runtime error: Index '1627' of dimension 1 of array 'rowvalues' above upper bound of 1626" (The array size of 'rowValues', 'colValues' and 'sValues (or n_S)' in the mapfile I used is 1626.)

After incorporating this fix, coupled MALI - Sea Level Model simulations compiled in DEBUG mode do not experience any errors, and the simulation results (i.e., MALI thickness data interpolated between the MALI mesh and the global Gaussian grid) are correct.

Thanks again, @knbarton, @mark-petersen, @sbrus89!

[sbrus89]

Thanks for pointing this out and helping to test the fix @hollyhan!

[mark-petersen]

@philipwjones can you check the acc changes I just made in 03a5551?

[philipwjones]

Changes for acc look fine. Approving based on all the other discussion/testing.

[jonbob]

@sbrus89 - the scripts that update the bld files to be consistent with Registry changes only came up with a missing blank line in namelist_definition_mpaso.xml. Is it OK to push to your branch? Or I can tell you where the line is if you prefer to handle it

[sbrus89]

@jonbob, thanks for catching that. Sure, it's totally fine to push to my branch.

[rljacob]

Is this ready?

[sbrus89]

As far as I know, it is.

[jonbob]

Passes:

• ERS.ne11_oQU240.WCYCL1850NS.chrysalis_intel
• SMS_D_Ld3.T62_oQU120.CMPASO-IAF.chrysalis_intel
• SMS_D_Ld1.ne30pg2_r05_EC30to60E2r2.WCYCL1850.chrysalis_intel
• PEM_Ln9.ne30pg2_EC30to60E2r2.WCYCL1850.chrysalis_intel

with expected NML DIFFs.

[jonbob]

merged to next

[jonbob]

Missed changes that needed to go in the mpas-ocean build-namelist file. Retesting passes:

• ERS.ne11_oQU240.WCYCL1850NS.chrysalis_intel
• SMS_D_Ld3.T62_oQU120.CMPASO-IAF.chrysalis_intel
• SMS_D_Ld1.ne30pg2_r05_EC30to60E2r2.WCYCL1850.chrysalis_intel
• PEM_Ln9.ne30pg2_EC30to60E2r2.WCYCL1850.chrysalis_intel

Still with expected NML DIFFs. Re-merged to next

[jonbob]

merged to master and expected NML DIFFs blessed