Add advection operator convergence test

[xylar]

This merge contains the COMPASS changes from: MPAS-Dev/MPAS-Model#583

Adds the Skamarock and Gassmann cosine bell test case and a regression
suite across resolutions to test operator convergence

This PR requires the addition of the cosine bell test case in MPAS-Model init mode, which is in MPAS-Dev/MPAS-Model#583, or changing the initial_state step to a python script.

Reference: Section 2a in Skamarock, W.C. and A. Gassmann, 2011: Conservative Transport Schemes for Spherical Geodesic Grids: High-Order Flux Operators for ODE-Based Time Integration. Mon. Wea. Rev., 139, 2962–2975, https://doi.org/10.1175/MWR-D-10-05056.1

[xylar]

We will want to take out the COMPASS changes from MPAS-Dev/MPAS-Model#583. Then, that PR and this one will be our first test of making changes to both MPAS-Ocean init mode and COMPASS at the same time. Let's see how it goes...

[vanroekel]

approved by visual inspection and verifying test set up correctly on grizzly

[mark-petersen]

Well, I tricked myself the very first time. I tried to run this with the head of ocean/develop, forgetting that this PR depends on changes from MPAS-Model init mode. Unfortunately, there is no MPAS-Ocean init mode error if you pick a test case that doesn't exist. Instead, you get an indecipherable error during simulation start-up:

App launch reported: 1 (out of 1) daemons - 2 (out of 2) procs
[gr1281:28737] *** An error occurred in MPI_Isend
[gr1281:28737] *** reported by process [2296512513,1]
[gr1281:28737] *** on communicator MPI COMMUNICATOR 3 DUP FROM 0
[gr1281:28737] *** MPI_ERR_COUNT: invalid count argument
[gr1281:28737] *** MPI_ERRORS_ARE_FATAL (processes in this communicator will now abort,
[gr1281:28737] ***    and potentially your MPI job)

Once I pointed the compass config to the executable and namelists from MPAS-Dev/MPAS-Model#583, it works fine.

This is a good chance to switch a test case from init mode to a python script. I have a working example with xarray here:
https://github.com/mark-petersen/MPAS-Model/blob/internal_tide_test_case/testing_and_setup/compass/ocean/internal_tide/3km/RK4/add_initial_state.py
that could be altered for this test case.

[xylar]

@mark-petersen, in your example, you use things like:

    vars2D = ['ssh', 'bottomDepth', 'bottomDepthObserved',
        'surfaceStress', 'atmosphericPressure', 'boundaryLayerDepth']
    for var in vars2D:
        globals()[var] = np.nan * np.ones(nCells)

This is fine for a one-off but I don't think this is acceptable as a general practice for python coding. anytime you find yourself using globals(), there's got to be a better, more "pythonic" way to do it. I'm only pointing this out because I don't want this example to become a template for future test cases initialized in python.

That being said, I'm working on initializing a test case in python right now as well, and I really appreciate this example for comparison, since I'm running into a bit of trouble.

Given the work @vanroekel has already put into MPAS-Dev/MPAS-Model#583, I think the right thing to do would be to rebase that branch, removing changes in testing_and_setup/compass and to merge that PR and this one at roughly the same time, rather than trying to remove the init-mode changes in this case. We can encourage that approach in the future, though.

[xylar]

@mark-petersen, Regarding the confusing error, can you track down why this isn't giving a more useful error message and fix it somehow? Or at least create an issue on MPAS-Model?

[mark-petersen]

I reran with the rebased code and executable from MPAS-Dev/MPAS-Model#583. Everything works great. Here is that wonderful convergence plot:

see

/lustre/scratch3/turquoise/mpeterse/runs/nightly/ocean_model_201114_cosine_bell2/ocean/convergence_global

[mark-petersen]

I used these lines:

    vars2D = ['ssh', 'bottomDepth', 'bottomDepthObserved',
        'surfaceStress', 'atmosphericPressure', 'boundaryLayerDepth']
    for var in vars2D:
        globals()[var] = np.nan * np.ones(nCells)

because otherwise for xarray you need one line per variable at the top, and another one at the bottom, and you have to ensure that the dimensions match. That seemed like a lot of bookkeeping, so I switched it to a few lists instead.

[mark-petersen]

See MPAS-Dev/MPAS-Model#758 for issue.

[mark-petersen]

Ran convergence test with rebased PR here: MPAS-Dev/MPAS-Model#583 with gnu debug an optimized, produces same plot as above. Takes 12 minutes on grizzly with 4 nodes.

[pep8speaks]

Hello @xylar! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file ocean/convergence_global/analysis/cosine_bell/compute_rmse.py:

Line 7:1: E402 module level import not at top of file

Comment last updated at 2021-01-20 12:11:47 UTC

[xylar]

I rebased and update the submodule for MPAS-Model/ocean/develop. I reran the test suite on Anvil. The results look similar to before:

But I wanted to note that the order of convergence is noticeably lower and the error in each test is noticeably higher than in September when @mark-petersen last posted a plot: MPAS-Dev/MPAS-Model#583 (review)

Did parameters change that would account for this? Longer integration time?

[xylar]

@vanroekel and @mark-petersen, thank you for your reviews earlier. I'm going to go ahead and merge because the order of convergence isn't related to this PR. But it would have been good if this had been noticed and discussed before MPAS-Dev/MPAS-Model#583 went in. It also would have been good if we had had this regression suite since September as a tool to find out when things changed. Part of our slow ocean/develop to E3SM process that I think is hindering us more than it's helping.