Add baroclinic gyre tests #547

alicebarthel · 2023-02-28T19:45:12Z

This PR aims to add the MITgcm baroclinic test case to compass. Started as a Work In Progress [WIP] PR to allow discussion and debugging before merging.

Checklist

User's Guide has been updated
Developer's Guide has been updated
API documentation in the Developer's Guide (api.rst) has any new or modified class, method and/or functions listed
Documentation has been built locally and changes look as expected
Document (in a comment titled Testing in this PR) any testing that was used to verify the changes

xylar · 2023-03-07T18:14:06Z

Nice! All tests passing!

compass/ocean/__init__.py

cbegeman · 2023-03-07T19:36:40Z

compass/ocean/tests/mitgcm_baroclinic_gyre/__init__.py

+        """
+        super().__init__(mpas_core=mpas_core, name='mitgcm_baroclinic_gyre')
+
+        for resolution in ['80km']:


Consider adding a comment that resolutions here must be added with the 'km' suffix (since later you chop off the last 2 characters to get the float)

An alternative would be to provide the resolution as a float and add the km suffix to the subdirectory. I think that might be cleaner. (I know we were following ziso, so that one should probably be fixed in a similar way.)

By the way, I have something along these lines here, where the resolution is always given in meters

compass/compass/ocean/tests/turbulence_closure/default/__init__.py

Lines 38 to 41 in f8af494

if resolution >= 1e3:

res_name = f'{int(resolution/1e3)}km'

else:

res_name = f'{int(resolution)}m'

@cbegeman Working on this resolution change now. A few questions/suggestions:

Are you satisfied with resolution in meters by default? It makes sense for turbulence closure since we go to fine resolution. Is it best to be consistent across cases? Happy to switch mine, just wanted to check.

the comment blocks in your forward (L12, L25) and initial_state (L22, L34) still carry comments where resolution is expected to be string, though it is the float that is passed (as far as I understand). It may be good to propagate the update.

The run in forward still assumes a string resolution with units. I couldn't find a call for it so maybe it's just stale code, but worth updating in case it ever get exercised.

If the resolution variable is only present in your case-specific python code (e.g., not in the config file), I think you can use whatever units are convenient for your cases. We care more about standardizing dc (in config files) and the resolution directory name across cases.

Thanks for pointing this out. I'm not actively developing that branch (and case) any more because we didn't have success validating mpas_les. If I return to it, I'll fix this up.

Same as (2)

cbegeman · 2023-03-07T19:40:38Z

compass/ocean/tests/mitgcm_baroclinic_gyre/initial_state.py

+    bottT = temperature[0, 200, -1].values
+    surfT = temperature[0, 200, 0].values
+    print(f'bottom T: {bottT} and surface : {surfT}')
+    salinity = 34.0 * xr.ones_like(temperature)


Consider making certain initial temperature and salinity parameters a part of the mitgcm_baroclinic_gyre section in your config file

The initial temperature is a function of depth. @alicebarthel has constructed a function that's a pretty good approximation of the table of numbers used in the MITgcm version of the test case. It's full of "magic" numbers so they wouldn't be super easy to add to the config file. But maybe there would be a way to give the top (z=0) and bottom (z=bottom_depth) temperature and have them vary between those two with your prescribed functional form? That would seem cleaner to me.

I wouldn't hardcode sampling bottT and surfT at the 200th grid cell here. Either you should take grid cell 0 or you'll want to have a number somehow based on nCells. We don't necessarily know that the mesh will have a 200th grid cell (though it's obviously highly likely) since we're flexible about the resolution.

Regarding the initial salinity, I agree for sure that that needs to be a config option rather than a magic number here.

Good point @xylar re. cell 200. @cbegeman I like the idea of it being more general (setting T in cfg) but I have not come up with a satisfying way to do that at this stage. For the original purpose of matching MITgcm, hard-coding the structure was easiest. It can extend to a deeper set-up, but approximates their values over 0-1800m. I'll look to see about setting it from top and bottom T keeping the existing functional form. We could do a prettier functional form with, say, a thermocline depth parameter setting the change in T=f(z) shape, but I have not looked into it that closely.

@alicebarthel, couldn't you have config options that are the temperatures in the middle of the top and bottom layers respectively? Then, couldn't you take your existing functional form but stretch it linearly so it hits the desired top and bottom temperatures?

Alternatively, if you want to be very particular about having the same temperature profile independent of resolution, couldn't you do the same but instead of the middle of the top and bottom layers, you would use z = 0 and z=1800 m. The default values would be chosen such that the temperature in the middle of the top and bottom layers is quite close to the value you currently get with your functional form.

Moved the constant salinity value to the cfg file. The initial temperature profile is now set (still with a "magical" functional form and a fitted parameter) from the top and bottom temperature set at the surface (z=0m) and bottom (z = bottom_depth). The default values approximate the mitgcm stratification, but we could add a table to reach to be more exact.

compass/ocean/tests/mitgcm_baroclinic_gyre/initial_state.py

compass/ocean/tests/mitgcm_baroclinic_gyre/gyre_test_case.py

cbegeman · 2023-03-07T19:52:52Z

@alicebarthel I think this is coming together really nicely! I left only minor comments at this point.

alicebarthel · 2023-03-11T00:23:30Z

Thanks @cbegeman and @xylar for your comments so far!

compass/ocean/tests/mitgcm_baroclinic_gyre/initial_state.py

alicebarthel · 2024-03-28T21:10:57Z

I made some offline plots for sanity checks on the initial and forcing files. The more useful ones are these 3:

They contain numbers (mean, std) on the other variables to flag if something is modified in other fields.

alicebarthel · 2024-03-29T01:08:31Z

I made some (offline) plots to check the simulation state -- as a sanity check. They still need minor adjustments (e.g. I am unsure about the surface heat flux calculation) but are promising. At this stage, the "mean" surface state was calculated from 30-day snapshots over year 3 (the last year of the 3_year_test).

I plotted a spin-up plot up to year 33. As expected, 3 years is far from enough for the spin-up (the MITgcm example show year 100).

I'll update when I have later results.

xylar

I had some issues mainly with the docs. A few things to clean up there.

compass/ocean/tests/baroclinic_gyre/forward.py

docs/developers_guide/ocean/api.rst

docs/developers_guide/ocean/test_groups/baroclinic_gyre.rst

docs/users_guide/ocean/test_groups/baroclinic_gyre.rst

docs/developers_guide/ocean/api.rst

docs/users_guide/ocean/test_groups/baroclinic_gyre.rst

xylar · 2024-03-29T12:10:22Z

I successfully ran the 80km/performance_test, 80km/3_year_test and 20km/performance_test. I didn't try to run the 20km/3_year_test.

xylar · 2024-03-29T12:12:06Z

@alicebarthel, the outputs you posted above look great to me! I think we're almost there...

alicebarthel · 2024-04-01T15:03:22Z

Thanks @xylar, that's super helpful. Attached are updated mean state plots for year 100 of the 80km. It looks nicely spun-up to me, at least dynamically:

cbegeman · 2024-04-01T15:59:10Z

@alicebarthel This is awesome progress! Is this ready for review? Are you planning on adding a visualization or analysis step to aid comparison with MITgcm https://mitgcm.readthedocs.io/en/latest/examples/baroclinic_gyre/baroclinic_gyre.html#model-solution? The heat flux figure you posted looks like a good comparison with Fig 2.7. Do you think we want the analogs to Figs 4.8 or 4.9?

alicebarthel · 2024-04-01T17:52:48Z

Thanks @cbegeman! For now, all the analysis is offline. Yes, there are a few more plots I want to make, including Figure 4.9b. Figure 4.8 may be good too - although at this stage the SSH contour looks good enough for me. Adding the barotropic streamfunction could be nice but feels lower priority for me.
A question for you two:
3 years is not enough for a full spin-up. The plots for MITgcm were made for year 100.

Does it make sense to keep the 3_year_test even if it does not bring things to equilibrium? The expectation is for users to run it further as needed. For 80km, it looks like ~50 years is needed.
Does it make sense to add an analysis/visualization step if the default run is not usually spun-up?

At this stage, I am focusing on:

Testing/debugging the restart capabilities. I have had issues but I suspect this could be a simple storage issue at this stage.
enhancing the offline analysis to increase our confidence of the spun-up state (including at 20-km resolution) and include a readily available AMOC estimate.

I welcome comments and suggestions as I work on this.

cbegeman · 2024-04-01T19:20:01Z

Does it make sense to keep the 3_year_test even if it does not bring things to equilibrium? The expectation is for users to run it further as needed. For 80km, it looks like ~50 years is needed.
Does it make sense to add an analysis/visualization step if the default run is not usually spun-up?

How about keeping the long test and having analysis/viz steps included in that test (and plot the last time slice in the output) but not include these steps in the performance test? I think it's reasonable to have the users change the duration. I didn't read through the docs yet, but let's make sure that ~50 year duration is recommended there.

xylar · 2024-04-03T20:43:32Z

compass/ocean/tests/baroclinic_gyre/forward.py

+        if self.name == '3_year_test':
+            replacements = {'config_do_restart': '.true.',
+                            'config_start_time': "'file'"}
+            self.update_namelist_at_runtime(replacements)


This isn't working because the name of the step is always forward, not 3_year_test. One fix is:

Suggested change

if self.name == '3_year_test':

replacements = {'config_do_restart': '.true.',

'config_start_time': "'file'"}

self.update_namelist_at_runtime(replacements)

if self.test_case.name == '3_year_test':

replacements = {'config_do_restart': '.true.',

'config_start_time': "'file'"}

self.update_namelist_at_runtime(replacements)

Or I think it might be more logical to store the long argument as an attribute and use it here:

Suggested change

if self.name == '3_year_test':

replacements = {'config_do_restart': '.true.',

'config_start_time': "'file'"}

self.update_namelist_at_runtime(replacements)

if self.long:

replacements = {'config_do_restart': '.true.',

'config_start_time': "'file'"}

self.update_namelist_at_runtime(replacements)

Merge branch 'add-mitgcm-moc-analysis' into add-mitgcm-baroclinic-gyre

xylar · 2024-08-12T17:11:39Z

compass/ocean/tests/baroclinic_gyre/initial_state.py

+    print(f'analytical bottom T: {val_bot} at \
+            depth : {bottom_depth}')


Here and in similar places, this is the preferred way to continue an f-string on a subsequent line:

Suggested change

print(f'analytical bottom T: {val_bot} at \

depth : {bottom_depth}')

print(f'analytical bottom T: {val_bot} at '

f'depth : {bottom_depth}')

Thanks @xylar , I modified these. Let me know if I missed others.

alicebarthel · 2024-08-13T21:08:27Z

@xylar @cbegeman This is ready for a review. I ran the performance and 3_year tests on Chrysalis without issue. I also ran the 3_year test for longer (50years, which takes ~2hours for 80km set-up).
Thanks @xylar for working with me on the heat fluxes, this looks much better!
Results after 3 years:

Results after 50 years (average over last 10years):

xylar · 2024-08-14T16:10:57Z

@alicebarthel, that looks awesome! I'm glad we got the surface restoring worked out.

I'll review it tomorrow.

xylar

@alicebarthel, this is some amazing work! It's very nearly there!

I was able to run the 80km tests without a hitch! I'm running the 20 km tests now. In the meantime, I have several small suggestions that you could make, or we could discuss further.

compass/ocean/tests/baroclinic_gyre/initial_state.py

xylar · 2024-08-15T12:51:56Z

docs/developers_guide/ocean/api.rst

+   BaroclinicGyre
+
+   GyreTestCase
+   GyreTestCase.configure
+
+   forward.Forward
+   forward.Forward.run
+
+   initial_state.InitialState
+   initial_state.InitialState.run


The other classes should be added here as well (CullMesh, Moc and Viz).

Thanks! Is this supposed to be alphabetical or does the order matter?

Yes, alphabetical is easiest.

xylar · 2024-08-15T12:53:38Z

docs/developers_guide/ocean/test_groups/baroclinic_gyre.rst

+At this stage, the test case is available at 80-km and 20-km horizontal
+resolution.  By default, the 15 vertical layers vary linearly in thickness with depth, from 50m at the surface to 190m at depth (full depth: 1800m).
+
+The test group includes 2 test cases, called ``performance_test`` for a short (10-day) run, and ``3_year_test`` for a 3-year simulation. Note that 3 years are insufficient to bring the standard test case to full equilibrium.  Both test cases have 2 steps,


Suggested change

The test group includes 2 test cases, called ``performance_test`` for a short (10-day) run, and ``3_year_test`` for a 3-year simulation. Note that 3 years are insufficient to bring the standard test case to full equilibrium. Both test cases have 2 steps,

The test group includes 2 test cases, called ``performance_test`` for a short (3-time-step) run, and ``3_year_test`` for a 3-year simulation. Note that 3 years are insufficient to bring the standard test case to full equilibrium. Both test cases have 2 steps,

xylar · 2024-08-15T12:57:31Z

docs/developers_guide/ocean/test_groups/baroclinic_gyre.rst

+Framework
+---------


It seems like cull_mesh, moc and viz need to be added to this.

xylar · 2024-08-15T13:39:19Z

compass/ocean/tests/baroclinic_gyre/viz.py

+        maxloc = 'at lat = {} and z = {}m'.format(
+            latbins[idx[-1]].values, int(dsMesh.refInterfaces[idx[0]].values))
+        maxval = 'max MOC = {:.1e} at lat={}-{}'.format(
+            round(np.max(moc[:, 175]), 1),
+            latbins[175 - 1].values, latbins[175].values)


Ideally, break this into several lines where you store these various inputs in variables and then use f-strings rather than the .format() methods to put them in place.

compass/ocean/tests/baroclinic_gyre/viz.py

cbegeman

I think this looks great! I just visually inspected the code. Would you like me to run tests on any particular machine or build docs?

cbegeman · 2024-08-15T14:37:59Z

compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg

+# Number of vertical levels
+vert_levels = 15
+
+# Total water column depth


Suggested change

# Total water column depth

# Total water column depth in m

cbegeman · 2024-08-15T14:38:34Z

compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg

+# the type of vertical grid
+grid_type = linear_dz
+
+# the linear rate of thickness increase for linear_dz


cbegeman · 2024-08-15T14:39:21Z

compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg

+restoring_temp_timescale = 30.
+
+# config options for the mean state visualization
+[mean_state_viz]


Suggested change

[mean_state_viz]

[baroclinic_gyre_viz]

Thanks! I modified this to include both viz and moc calculation parameters.

compass/ocean/tests/baroclinic_gyre/cull_mesh.py

compass/ocean/tests/baroclinic_gyre/initial_state.py

cbegeman · 2024-08-15T16:02:32Z

docs/developers_guide/ocean/test_groups/baroclinic_gyre.rst

+:py:class:`compass.mesh.QuasiUniformSphericalMeshStep`.  Then, the mesh is
+culled to keep a single ocean basin (with lat, lon bounds set in `.cfg` file).  A vertical grid is generated,
+with 15 layers of thickness that increases linearly with depth (10m increase by default), from 50m at the surface to 190m at depth (full depth: 1800m).
+Finally, the initial temperature field is initialized with a vertical profile to approximate the discrete values set in the `MITgcm test case <https://mitgcm.readthedocs.io/en/latest/examples/baroclinic_gyre/baroclinic_gyre.html>`_, uniform in horizontal space. The surface and bottom values are set in the `.cfg` file. Salinity is set to a constant value (34 psu, set in the `.cfg` file)  and initial


Break up into multiple lines here and elsewhere. I think max length should be 89 char.

Yes, 79 characters actually. The width of a default terminal is 80, so that's the reasoning.

alicebarthel · 2024-08-16T21:07:56Z

docs/users_guide/ocean/test_groups/index.rst

@@ -12,6 +12,7 @@ coming months.
   :titlesonly:

   baroclinic_channel
+   baroclinic_gyre


@xylar Is this what you meant in your review of the users_guide above?

xylar

Looks great! Thanks for making the requested changes.

alicebarthel added enhancement New feature or request ocean in progress This PR is not ready for review or merging labels Feb 28, 2023

alicebarthel force-pushed the add-mitgcm-baroclinic-gyre branch 2 times, most recently from 4ccb04f to a3baa9c Compare March 7, 2023 17:35

cbegeman reviewed Mar 7, 2023

View reviewed changes

compass/ocean/__init__.py Outdated Show resolved Hide resolved

cbegeman reviewed Mar 7, 2023

View reviewed changes

compass/ocean/tests/mitgcm_baroclinic_gyre/initial_state.py Outdated Show resolved Hide resolved

cbegeman reviewed Mar 7, 2023

View reviewed changes

compass/ocean/tests/mitgcm_baroclinic_gyre/gyre_test_case.py Outdated Show resolved Hide resolved

xylar self-assigned this Mar 29, 2023

xylar reviewed Oct 30, 2023

View reviewed changes

compass/ocean/tests/mitgcm_baroclinic_gyre/initial_state.py Outdated Show resolved Hide resolved

xylar mentioned this pull request Nov 1, 2023

added linear dz in vertical grid #724

Draft

8 tasks

alicebarthel force-pushed the add-mitgcm-baroclinic-gyre branch from 0c341eb to 3bac95a Compare November 6, 2023 16:46

alicebarthel force-pushed the add-mitgcm-baroclinic-gyre branch from 4145dce to b976f7a Compare March 27, 2024 16:48

alicebarthel force-pushed the add-mitgcm-baroclinic-gyre branch from 85775b3 to 6a48d21 Compare March 29, 2024 00:05

xylar requested changes Mar 29, 2024

View reviewed changes

xylar reviewed Apr 3, 2024

View reviewed changes

xylar mentioned this pull request May 21, 2024

Switch wind stress to use isotropic interpolation from cells to edge normals E3SM-Ocean-Discussion/E3SM#98

Open

alicebarthel added 13 commits July 31, 2024 19:38

WIP - debugging moc analysis

d61771a

debugging moc calc

dd7db80

add time-varying moc

517a81c

added moc viz step

51e7422

added spinup timeseries plot

7e37003

removed some static fields from output and monthly streams

4f8ea0c

adding the moc analysis and visualization steps (moc and spinup)

a7776ba

Merge branch 'add-mitgcm-moc-analysis' into add-mitgcm-baroclinic-gyre

added averaging period for moc

1e74310

added time mean state plots

05f2b79

WIP - correction restoring velocity

1964db4

modified doc for restoring correction and estimated spin-up time

3d5fff7

correction to heat flux plot

e976366

WIP - added surface flux variable to streams as test

77b6e27

xylar reviewed Aug 12, 2024

View reviewed changes

alicebarthel added 3 commits August 12, 2024 12:18

correction to plot title

a85d655

corrected the heat flux viz

b7139d3

best practice modif for strings

5d571a9

alicebarthel force-pushed the add-mitgcm-baroclinic-gyre branch from 6aeafb9 to 5d571a9 Compare August 13, 2024 20:15

alicebarthel requested review from xylar and cbegeman August 13, 2024 21:04

xylar changed the title ~~Add mitgcm baroclinic gyre - WIP~~ Add baroclinic gyre tests Aug 15, 2024

xylar removed the in progress This PR is not ready for review or merging label Aug 15, 2024

xylar reviewed Aug 15, 2024

View reviewed changes

cbegeman approved these changes Aug 15, 2024

View reviewed changes

alicebarthel commented Aug 16, 2024

View reviewed changes

PR review changes

36f3f91

xylar approved these changes Aug 17, 2024

View reviewed changes

xylar merged commit 12f956f into MPAS-Dev:main Aug 17, 2024
4 checks passed

	if resolution >= 1e3:
	res_name = f'{int(resolution/1e3)}km'
	else:
	res_name = f'{int(resolution)}m'

		print(f'analytical bottom T: {val_bot} at \
		depth : {bottom_depth}')

	The test group includes 2 test cases, called ``performance_test`` for a short (10-day) run, and ``3_year_test`` for a 3-year simulation. Note that 3 years are insufficient to bring the standard test case to full equilibrium. Both test cases have 2 steps,
	The test group includes 2 test cases, called ``performance_test`` for a short (3-time-step) run, and ``3_year_test`` for a 3-year simulation. Note that 3 years are insufficient to bring the standard test case to full equilibrium. Both test cases have 2 steps,

		Framework
		---------

Add baroclinic gyre tests #547

Add baroclinic gyre tests #547

Conversation

alicebarthel commented Feb 28, 2023 • edited by xylar Loading

xylar commented Mar 7, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cbegeman commented Mar 7, 2023

alicebarthel commented Mar 11, 2023

alicebarthel commented Mar 28, 2024 • edited Loading

alicebarthel commented Mar 29, 2024 • edited Loading

xylar left a comment

Choose a reason for hiding this comment

xylar commented Mar 29, 2024

xylar commented Mar 29, 2024

alicebarthel commented Apr 1, 2024

cbegeman commented Apr 1, 2024

alicebarthel commented Apr 1, 2024

cbegeman commented Apr 1, 2024

Choose a reason for hiding this comment

xylar Aug 12, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alicebarthel commented Aug 13, 2024

xylar commented Aug 14, 2024

xylar left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

cbegeman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

xylar left a comment

Choose a reason for hiding this comment

alicebarthel commented Feb 28, 2023 •

edited by xylar

Loading

alicebarthel commented Mar 28, 2024 •

edited

Loading

alicebarthel commented Mar 29, 2024 •

edited

Loading

xylar Aug 12, 2024 •

edited

Loading