Add advection operator convergence test #583

vanroekel · 2020-06-01T20:41:26Z

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

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

vanroekel · 2020-06-01T20:42:16Z

@mark-petersen and @xylar here is my WIP advection operator convergence for the global advective test case. If you have some time, I'd greatly appreciate your comments.

vanroekel · 2020-06-01T20:42:36Z

also tagging @jamilgafur

xylar · 2020-06-01T20:46:14Z

testing_and_setup/compass/ocean/jigsaw_to_MPAS/build_mesh.py

-    write_netcdf(convert(xarray.open_dataset('mesh_triangles.nc')),
-                 'base_mesh.nc')
+    write_netcdf(convert(xarray.open_dataset('mesh_triangles.nc'),
+                    graphInfoFileName='graph.info'),'base_mesh.nc')


Perfect! My only comment here would be to indent another 5 spaces so the formatting follows the PEP8 convention.

Just so I'm clear, indent the second line 5 more spaces? Is this so it is easier to see it is part of the convert function?

If you load e3sm_unified you can run

autopep8 --in-place --aggressive --aggressive foo.py

on your python files for automatic pep8 spacing.

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/config_driver.xml

xylar · 2020-06-01T20:51:04Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+            ii = line.find('=')+1
+            pd = float(line[ii:])
+
+    init = xr.open_dataset('../../'+resTag+'/default/init_step/initial_state.nc')


I really prefer people alter strings with .format() instead of string "addition". This may be may own prejudice but I don't see string addition used by seasoned python programmers.

Suggested change

init = xr.open_dataset('../../'+resTag+'/default/init_step/initial_state.nc')

init = xr.open_dataset('../../{}/default/init_step/initial_state.nc'.format(resTag))

xylar · 2020-06-01T20:51:22Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+
+    init = xr.open_dataset('../../'+resTag+'/default/init_step/initial_state.nc')
+    #find time since the beginning of run
+    ds = xr.open_dataset('../../'+resTag+'/default/forward/output/output.0001-01-01_00.00.00.nc')


Similar to above, maybe use .format().

xylar · 2020-06-01T20:53:56Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+    for i in range(len(latC)):
+        temp = R*np.arccos(np.sin(latCent)*np.sin(latC[i]) + np.cos(latCent)*np.cos(latC[i])*np.cos(lonC[i] - newLon))
+        if temp < radius:
+            tracer[i] = psi0 / 2.0 * (1.0 + np.cos(3.1415926*temp/radius))


Never loop if you don't have to in python.

Suggested change

for i in range(len(latC)):

temp = R*np.arccos(np.sin(latCent)*np.sin(latC[i]) + np.cos(latCent)*np.cos(latC[i])*np.cos(lonC[i] - newLon))

if temp < radius:

tracer[i] = psi0 / 2.0 * (1.0 + np.cos(3.1415926*temp/radius))

temp = R*np.arccos(np.sin(latCent)*np.sin(latC) + np.cos(latCent)*np.cos(latC)*np.cos(lonC - newLon))

mask = temp < radius

tracer[mask] = psi0 / 2.0 * (1.0 + np.cos(numpy.py*temp[mask]/radius))

that's nice! thanks for the suggestion

xylar · 2020-06-01T20:57:20Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+import matplotlib.pyplot as plt
+import os 
+
+%matplotlib inline


I haven't seen this one before.

this was an error in importing from jupyter notebook to python script. It's ipython magic. I've removed this.

xylar · 2020-06-01T20:58:54Z

@vanroekel, most of my comments are just little things that you're free to take or leave. Looks like a decent approach to me. Happy to do some testing once the time comes.

vanroekel

Thanks @xylar I will make the changes. As I test this, I'm having a hard time controlling the plot size and position of the convergence label (annotate) command. Do you have any suggestions? Here is what a plot looks like right now

xylar · 2020-06-01T21:04:08Z

Thanks @xylar I will make the changes. As I test this, I'm having a hard time controlling the plot size and position of the convergence label (annotate) command. Do you have any suggestions? Here is what a plot looks like right now

Oh, yeah, that doesn't look great. Let me see what I can figure out.

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

xylar · 2020-06-01T21:07:00Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+matplotlib.rcParams['figure.figsize'] = (16, 16) # Large figures
+dpi=200;
+
+### axis_font = {'fontname':'Arial', 'size':'18'}    
+title_font = {'fontname':'Arial', 'size':'32', 'color':'black', 'weight':'normal'}
+matplotlib.rc('xtick', labelsize=28)
+matplotlib.rc('ytick', labelsize=28)


I think this is the stuff that's causing your figure to look bad. Try taking it out. I'll make a few recommendations below once you do that.

xylar · 2020-06-01T21:08:24Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+p = np.polyfit(np.log10(xdata),np.log10(ydata),1)
+conv = abs(p[0])
+
+plt.loglog(xdata,ydata,'k')


If the figure isn't the size you want, set the size here:

plt.figure(figsize=[<width>, <height>])

where width and height are in inches.

xylar · 2020-06-01T21:09:09Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+plt.loglog(xdata,ydata,'k')
+plt.scatter(xdata,ydata,s=40,c='r')
+plt.annotate('Order of Convergence = {}'.format(np.round(conv,2)),xy=((xdata[-1]+xdata[0])*0.4,ydata[0]),fontsize=28)
+plt.savefig('convergence.png')


If the png is low resolution, set dpi=200 here. I wouldn't bother doing that earlier unless you're going to use plt.show().

xylar · 2020-06-01T21:13:30Z

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py

+
+plt.loglog(xdata,ydata,'k')
+plt.scatter(xdata,ydata,s=40,c='r')
+plt.annotate('Order of Convergence = {}'.format(np.round(conv,2)),xy=((xdata[-1]+xdata[0])*0.4,ydata[0]),fontsize=28)


28 is a big font size. And it might work better to just put it in the upper left corner or something. There are probably other tricks for putting it "where it fits bets" like is done for legends. If you want that, I can look into it.

vanroekel · 2020-06-01T21:53:15Z

thanks @xylar! Your comments have been extremely helpful. Here is an updated plot

vanroekel · 2020-06-02T05:05:40Z

@mark-petersen I've run through the advection test and here is the result

The convergence order is much lower than expected. When you have a chance could you look through this test and see if I have an error?

mark-petersen · 2020-06-02T12:26:06Z

@vanroekel thanks for your work on this. The low order of convergence is certainly surprising. The direct comparison is Skamarock and Gassmann figure 4:

They don't give numbers for the computed order on this plot, but everything is at least 2. MPAS-Ocean should be using the fourth order scheme. I would expect that limiting is not turning on, but we could check that.

vanroekel · 2020-06-02T13:29:28Z

@mark-petersen monotone is on here. We also are using the blended 3rd/4th order scheme (coeff3rdOrder =0.25). I've tried for various values of that coefficient and it doesn't help convergence order. I've also tried with non monotone and the order is roughly the same.

vanroekel · 2020-06-02T15:20:43Z

One other thing, there is a pretty strong time step sensitivity. Cutting dt to 1 min changes convergence strongly (not for the better). I'll add some plots of the tracer vs exact solution shortly

vanroekel · 2020-06-03T20:31:15Z

a quick update, using RK4 improves convergence, but only slightly

jamilgafur · 2020-06-08T18:57:52Z

src/core_ocean/mode_init/Makefile

             mpas_ocn_init_hurricane.o \
-             mpas_ocn_init_tidal_boundary.o
+             mpas_ocn_init_tidal_boundary.o \
+						 mpas_ocn_init_cosine_bell.o


Small edit, but spacing here is weird.

mark-petersen · 2020-09-16T14:31:56Z

Rebased, compiled and tested with gnu and intel debug for initial condition. Ran QU240 test through forward.

mark-petersen

This looks great! Thanks for your work @vanroekel. I cleaned up the COMPASS files so they are linked to templates rather than a bunch of individual files. Then it was very easy to add more resolutions in between. This still spans QU60 to QU240:

This is beautiful! The extra test cases give us more confidence in the slope.

I also was able to clean up the simulations. It had been writing output and restart files very frequently. It now can run all seven cases in 12 minutes on 4 nodes on grizzly with gnu.

mark-petersen · 2020-09-16T23:52:06Z

@jonbob This one is ready to merge when we have an open E3SM slot. It changes both COMPASS and init mode Fortran code.

xylar · 2020-11-06T12:41:02Z

I have a draft PR for moving the COMPASS changes from this PR to the new COMPASS repo: MPAS-Dev/compass#8

Once MPAS-Dev/compass#2 has been merged, we need to:

rebase Add advection operator convergence test compass#8
remove COMPASS changes from this PR
Re-test within the new COMPASS branch, pointing to a built of MPAS-Ocean from this PR
Remove the COMPASS label

xylar · 2020-11-14T19:52:28Z

remove COMPASS changes from this PR

@vanroekel, would you like to remove the COMPASS changes from this PR, or would you prefer I do it? Here is a link that I think describes the approach we would use: https://stackoverflow.com/a/17824718/7728169

xylar · 2020-11-14T21:56:58Z

@vanroekel, I went ahead and tried the rebase. The link I posted above wasn't much help, I just needed to keep deleting modified contents in testing_and_setup/compass. Here's my version of the branch, which just has 2 commits:
https://github.com/xylar/MPAS-Model/tree/ocean/addCosineBellTest
I think the second commit could be squashed with the first.

vanroekel · 2020-11-15T05:24:40Z

@xylar thanks for taking the time to pull the COMPASS changes out of that branch. It looks good to me.

xylar · 2020-11-15T08:38:18Z

@vanroekel, I'm going to take that as tacit permission to alter your branch to update this PR.

xylar · 2020-11-15T08:41:33Z

@mark-petersen, please merge this when it is convenient for E3SM. Then, we can update the ocean/develop submodule in MPAS-Dev/compass#8 and merge that PR.

mark-petersen · 2020-11-15T16:19:35Z

@xylar, thanks for rebasing. Hopefully we can continue to merge these peripheral PRs during the mushy freeze.

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

mark-petersen · 2021-01-15T14:18:16Z

Rebased. Corrected author and test case names. Tested successfully using COMPASS case, which is here: MPAS-Dev/compass#8, with gnu optimized and debug.

Merge PR #583 'vanroekel/ocean/addCosineBellTest' into ocean/develop

mark-petersen · 2021-01-15T21:42:13Z

Rebased version passes nightly regression suite with gnu debug. Ran convergence test here: MPAS-Dev/compass#8 with gnu debug an optimized, produces same plot as above. Takes 12 minutes on grizzly with 4 nodes.

Add cosine bell advection test case This PR brings in a new mpas-source submodule with changes only to the ocean core. It adds initial condition for cosine bell test case, used for converge test of tracer advection. See MPAS-Dev/MPAS-Model#583 and MPAS-Dev/compass#8. [BFB]

Merge PR MPAS-Dev#583 'vanroekel/ocean/addCosineBellTest' into ocean/develop

vanroekel added Ocean COMPASS in progress labels Jun 1, 2020

vanroekel requested review from mark-petersen and xylar June 1, 2020 20:41

vanroekel marked this pull request as draft June 1, 2020 20:41

xylar reviewed Jun 1, 2020

View reviewed changes

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/config_driver.xml Outdated Show resolved Hide resolved

xylar reviewed Jun 1, 2020

View reviewed changes

vanroekel commented Jun 1, 2020

View reviewed changes

xylar reviewed Jun 1, 2020

View reviewed changes

testing_and_setup/compass/ocean/cosine_bell/RMSE/default/compute_rmse.py Outdated Show resolved Hide resolved

xylar reviewed Jun 1, 2020

View reviewed changes

mark-petersen self-assigned this Jun 2, 2020

jamilgafur reviewed Jun 8, 2020

View reviewed changes

mark-petersen force-pushed the ocean/addCosineBellTest branch from 850201c to 45bda53 Compare September 16, 2020 14:30

mark-petersen approved these changes Sep 16, 2020

View reviewed changes

xylar mentioned this pull request Nov 6, 2020

Add advection operator convergence test MPAS-Dev/compass#8

Merged

xylar removed the COMPASS label Nov 14, 2020

xylar force-pushed the ocean/addCosineBellTest branch from 7db1d00 to 905dcd4 Compare November 15, 2020 08:39

Adds cosine bell case and advection test case

a01d197

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

mark-petersen force-pushed the ocean/addCosineBellTest branch from 905dcd4 to 0f1b805 Compare January 15, 2021 14:13

Correct comments and spacing

f3f049a

mark-petersen force-pushed the ocean/addCosineBellTest branch from 0f1b805 to f3f049a Compare January 15, 2021 14:16

mark-petersen mentioned this pull request Jan 15, 2021

Add cosine bell advection test case E3SM-Project/E3SM#4037

Merged

Add to cmake

f17c5d0

mark-petersen added a commit that referenced this pull request Jan 15, 2021

Merge branch 'ocean/develop' into e3sm/develop

bcebf2a

Merge PR #583 'vanroekel/ocean/addCosineBellTest' into ocean/develop

mark-petersen merged commit b3d870d into MPAS-Dev:ocean/develop Jan 19, 2021

mark-petersen added a commit to mark-petersen/MPAS-Model that referenced this pull request Mar 16, 2021

Merge branch 'ocean/develop' into e3sm/develop

80b705e

Merge PR MPAS-Dev#583 'vanroekel/ocean/addCosineBellTest' into ocean/develop

matthewhoffman added COMPASS in progress Ocean labels Mar 17, 2021

	init = xr.open_dataset('../../'+resTag+'/default/init_step/initial_state.nc')
	init = xr.open_dataset('../../{}/default/init_step/initial_state.nc'.format(resTag))

Add advection operator convergence test #583

Add advection operator convergence test #583

Uh oh!

Conversation

vanroekel commented Jun 1, 2020 • edited by mark-petersen Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

vanroekel commented Jun 1, 2020

Uh oh!

vanroekel commented Jun 1, 2020

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

xylar commented Jun 1, 2020

Uh oh!

vanroekel left a comment

Choose a reason for hiding this comment

Uh oh!

xylar commented Jun 1, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

xylar Jun 1, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

vanroekel commented Jun 1, 2020

Uh oh!

vanroekel commented Jun 2, 2020

Uh oh!

mark-petersen commented Jun 2, 2020

Uh oh!

vanroekel commented Jun 2, 2020

Uh oh!

vanroekel commented Jun 2, 2020

Uh oh!

vanroekel commented Jun 3, 2020

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mark-petersen commented Sep 16, 2020

Uh oh!

mark-petersen left a comment

Choose a reason for hiding this comment

Uh oh!

mark-petersen commented Sep 16, 2020

Uh oh!

xylar commented Nov 6, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

xylar commented Nov 14, 2020

Uh oh!

xylar commented Nov 14, 2020

Uh oh!

vanroekel commented Nov 15, 2020

Uh oh!

xylar commented Nov 15, 2020

vanroekel commented Jun 1, 2020 •

edited by mark-petersen

Loading

xylar commented Jun 1, 2020 •

edited

Loading

xylar Jun 1, 2020 •

edited

Loading

xylar commented Nov 6, 2020 •

edited

Loading