diff --git a/docs/index.rst b/docs/index.rst index a19a622c99..eeac3c13ad 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -73,6 +73,7 @@ components themselves evolves and new functionality is added. :maxdepth: 1 tutorials/dev_add_test_group + tutorials/dev_add_rrm tutorials/dev_add_param_study tutorials/dev_porting_legacy diff --git a/docs/tutorials/dev_add_rrm.rst b/docs/tutorials/dev_add_rrm.rst new file mode 100644 index 0000000000..d14b156b7a --- /dev/null +++ b/docs/tutorials/dev_add_rrm.rst @@ -0,0 +1,1987 @@ +.. _dev_tutorial_add_rrm: + +Developer Tutorial: Adding a new ocean/sea ice regionally refined mesh (RRM) +============================================================================ + +This tutorial presents a step-by-step guide to adding a new mesh to the +``global_ocean`` test group in ``compass`` (see the :ref:`glossary` for +definitions of some relevant terms). In this tutorial, I will add a new +mesh called YAM ("yet another mesh") that is based on techniques used to build +the existing :ref:`dev_ocean_global_ocean_sowisc12to60` as well as +:ref:`dev_ocean_global_ocean_wc14` meshes. + +.. _dev_tutorial_add_rrm_getting_started: + +Getting started +--------------- + +To begin with, you will need to check out the compass repo and create a new +branch from ``main`` for developing the new mesh. For this purpose, we +will stick with the simpler approach in :ref:`dev_compass_repo` here, but feel +free to use the ``worktree`` approach instead if you are comfortable with it. + +.. code-block:: bash + + git clone git@github.com:MPAS-Dev/compass.git add-yet-another-mesh + cd add-yet-another-mesh + git checkout -b add-yet-another-mesh + +Now, you will need to create a conda environment for developing compass, as +described in :ref:`dev_conda_env`. We will assume a simple situation where +you are working on a "supported" machine and using the default compilers and +MPI libraries, but consult the documentation to make an environment to suit +your needs. + +.. code-block:: bash + + # this one will take a while the first time + ./conda/configure_compass_env.py --conda $HOME/mambaforge --env_name compass_yam + +If you don't already have Mambaforge installed in the directory pointed to by +``--conda``, it will be installed automatically for you. If all goes well, you +will have a file named ``load_compass_yam*.sh``, where the details of the +``*`` depend on your specific machine and compilers. For example, on +Chrysalis, you will have ``load_compass_yam_chrysalis_intel_openmpi.sh``, +which will be the example used here: + +.. code-block:: bash + + source load_compass_yam_chrysalis_intel_openmpi.sh + +Now, we're ready to get the MPAS-Ocean source code from the E3SM repository: + +.. code-block:: bash + + # Get the E3SM code -- this one takes a while every time + git submodule update --init --recursive + +Next, we're ready to build the MPAS-Ocean executable: + +.. code-block:: bash + + cd E3SM-Project/components/mpas-ocean/ + make ifort + cd ../../.. + +The make target will be different depending on the machine and compilers, see +:ref:`dev_supported_machines` or :ref:`dev_other_machines` for the right one +for your machine. + +Now, we're ready to start developing! + +.. _dev_tutorial_add_rrm_add_mesh: + +Adding a new mesh +----------------- + +Use any method you like for editing code. If you haven't settled on a method +and are working on your own laptop or desktop, you may want to try an +integrated development environment (`PyCharm `_ +is a really nice one). They have features to make sure your code adheres to +the style required for compass (see :ref:`dev_style`). ``vim``, ``emacs`` or +a similar tool will work fine on supercomputers. Keep in mind that you can +edit code on your laptop or desktop but you will need to use a high-performance +computing (HPC) machine to generate a mesh. One step requires >=360 cores to +remap a high resolution topography dataset to the MPAS mesh. + +Your new mesh will be defined in a directory within the ``global_ocean`` test +group's ``mesh`` package. For this example, we create a new ``yam10to60`` +directory in ``compass/ocean/tests/global_ocean/mesh``. The directory name +should be all lowercase with no underscores. The convention is to give your +mesh a short prefix (``yam`` for "yet another mesh" in our case) and then +the finest and coarsest resolutions in km. In that directory, we +will make a new file called ``__init__.py``. In that file, we will start to +define a ``YAM10to60BaseMesh`` class. + +.. code-block:: bash + + cd compass/ocean/tests/global_ocean/mesh + mkdir yam10to60 + cd yam10to60 + vim __init__.py + +We define the resolution as a function of space using a regular +latitude-longitude map of the distance between cell centers in km. To begin +with , we will start out with something really simple: a quasi-uniform 60 km +mesh. We will make a 2D field called ``cellWidth`` on a 0.1 degree, regular +latitude-longitude grid. + +.. code-block:: python + + import numpy as np + + from compass.mesh import QuasiUniformSphericalMeshStep + + + class YAM10to60BaseMesh(QuasiUniformSphericalMeshStep): + """ + A step for creating YAM10to60 meshes + """ + + def build_cell_width_lat_lon(self): + """ + Create cell width array for this mesh on a regular latitude-longitude + grid + + Returns + ------- + cellWidth : numpy.array + m x n array of cell width in km + + lon : numpy.array + longitude in degrees (length n and between -180 and 180) + + lat : numpy.array + longitude in degrees (length m and between -90 and 90) + """ + + dlon = 0.1 + dlat = dlon + nlon = int(360. / dlon) + 1 + nlat = int(180. / dlat) + 1 + lon = np.linspace(-180., 180., nlon) + lat = np.linspace(-90., 90., nlat) + + cell_width = 60.0 * np.ones((nlat, nlon)) + + return cell_width, lon, lat + +We also need a config file for the mesh. For now, it can just be empty: + +.. code-block:: bash + + touch yam10to60.cfg + +Next, we need to add this mesh to the list of known meshes: + +.. code-block:: bash + + cd .. + vim __init__.py + +.. code-block:: python + :emphasize-lines: 5, 24-25 + + ... + + from compass.ocean.tests.global_ocean.mesh.so12to60 import SO12to60BaseMesh + from compass.ocean.tests.global_ocean.mesh.wc14 import WC14BaseMesh + from compass.ocean.tests.global_ocean.mesh.yam10to60 import YAM10to60BaseMesh + from compass.ocean.tests.global_ocean.metadata import ( + get_author_and_email_from_git, + ) + + ... + + class Mesh(TestCase): + + ... + + def __init__(self, test_group, mesh_name, remap_topography): + + ... + + elif mesh_name.startswith('Kuroshio'): + base_mesh_step = KuroshioBaseMesh(self, name=name, subdir=subdir) + elif mesh_name in ['WC14', 'WCwISC14']: + base_mesh_step = WC14BaseMesh(self, name=name, subdir=subdir) + elif mesh_name in ['YAM10to60', 'YAMwISC10to60']: + base_mesh_step = YAM10to60BaseMesh(self, name=name, subdir=subdir) + else: + raise ValueError(f'Unknown mesh name {mesh_name}') + +In compass, we support 2 versions of nearly every mesh, 1 where everything +south of the Antarctic coast is treated as land and 1 where we include the +ocean cavities below Antarctic ice shelves (``wISC``, meaning "with ice-shelf +cavities") + +Next, we add a test case for making this mesh to the ``global_ocean`` test +group: + +.. code-block:: bash + + cd .. + vim __init__.py + +.. code-block:: python + :emphasize-lines: 18-22 + + ... + + class GlobalOcean(TestGroup): + + ... + + def __init__(self, mpas_core): + + ... + + self._add_tests(mesh_names=['WC14', 'WCwISC14'], + DynamicAdjustment=WC14DynamicAdjustment) + + # Kuroshio meshes without ice-shelf cavities + self._add_tests(mesh_names=['Kuroshio12to60', 'Kuroshio8to60'], + DynamicAdjustment=KuroshioDynamicAdjustment) + + for mesh_name in ['YAM10to60', 'YAMwISC10to60']: + mesh_test = Mesh(test_group=self, mesh_name=mesh_name, + remap_topography=True) + self.add_test_case(mesh_test) + + # A test case for making E3SM support files from an existing mesh + self.add_test_case(FilesForE3SM(test_group=self)) + + +Okay, great! We're ready to test this out! + +.. _dev_tutorial_add_rrm_test_mesh: + +Running the mesh test case +-------------------------- + +First, let's make sure the mesh exists when we list the available test cases +in compass: + +.. code-block:: bash + + compass list + +You should see something like: + +.. code-block:: + :emphasize-lines: 7-8 + + ... + + 250: ocean/global_ocean/Kuroshio8to60/WOA23/init + 251: ocean/global_ocean/Kuroshio8to60/WOA23/performance_test + 252: ocean/global_ocean/Kuroshio8to60/WOA23/dynamic_adjustment + 253: ocean/global_ocean/Kuroshio8to60/WOA23/files_for_e3sm + 254: ocean/global_ocean/YAM10to60/mesh + 255: ocean/global_ocean/YAMwISC10to60/mesh + 256: ocean/global_ocean/files_for_e3sm + 257: ocean/gotm/default + + ... + +On one of the supported HPC machines, let's try setting up and running the mesh +step. We will concentrate on the mesh without ice-shelf cavities ``YAM10to60`` +for now because it is a little simpler. Here's an example that should work on +Chrysalis or Anvil: + +.. code-block:: bash + + compass setup -n 254 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_uniform60km + +The number provided (e.g. ``254``) is the number of the test case in +``compass list`` above. Since these numbers change frequently, you should +check each time you check out or rebase the code with something like +``compass list | grep YAM`` + +The directory you point to with ``-p`` is where you build MPAS-Ocean. We don't +actually need that quite yet but we will soon when we add an initial condition +on the mesh. The base work directory you point to with ``-w`` should be a +location in a scratch space where you can set up temporary tests. We recommend +including a date somewhere in the path just to keep things organized but that +is up to you. + +For simplicity, it is a good idea to open a new terminal for running the test. +In the new terminal window: + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_uniform60km + vim job_script.custom.sh + sbatch job_script.custom.sh + +If you prefer, you can get an interactive node and run the same commands as +in the job script: + +.. code-block:: bash + + source load_compass_env.sh + compass run custom + +Assuming you submitted the job script as above, you can monitor the output +file: + +.. code-block:: + + $ tail -f compass.o* + + Loading conda environment + Done. + + Loading Spack environment... + Done. + + ocean/global_ocean/YAM10to60/mesh + * step: base_mesh + * step: remap_topography + * step: cull_mesh + test execution: SUCCESS + test runtime: 11:31 + Test Runtimes: + 11:31 PASS ocean_global_ocean_YAM10to60_mesh + Total runtime 11:32 + PASS: All passed successfully! + +If things don't go well, you might see something like: + +.. code-block:: + + Loading conda environment + Done. + + Loading Spack environment... + Done. + + ocean/global_ocean/YAM10to60/mesh + * step: base_mesh + Failed + test execution: ERROR + see: case_outputs/ocean_global_ocean_YAM10to60_mesh.log + test runtime: 00:00 + Test Runtimes: + 00:00 FAIL ocean_global_ocean_YAM10to60_mesh + Total runtime 00:01 + FAIL: 1 test failed, see above. + +Hopefully, the contents of the log file, in this case +``case_outputs/ocean_global_ocean_YAM10to60_mesh.log``, can help you debug +what went wrong. + + +Once the ``base_mesh`` step has completed, you should see: + +.. code-block:: bash + + $ ls ocean/global_ocean/YAM10to60/mesh/base_mesh/ + base_mesh.nc geom.msh mesh.cfg opts.log + base_mesh_vtk graph.info mesh.msh spac.msh + cellWidthGlobal.png job_script.sh mesh_triangles.nc step.pickle + cellWidthVsLatLon.nc load_compass_env.sh opts.jig + +The main result is the file ``base_mesh.nc`` + +.. code-block:: bash + + $ cd ocean/global_ocean/YAM10to60/mesh/base_mesh/ + $ source load_compass_env.sh + $ ncdump -h base_mesh.nc + + netcdf base_mesh { + dimensions: + Time = UNLIMITED ; // (0 currently) + nCells = 165049 ; + nEdges = 495141 ; + nVertices = 330094 ; + maxEdges = 8 ; + maxEdges2 = 16 ; + TWO = 2 ; + vertexDegree = 3 ; + variables: + int edgesOnEdge(nEdges, maxEdges2) ; + + ... + +You can get take a look at the map of resolution in ``cellWidthGlobal.png``: + +.. image:: images/qu60.png + :width: 500 px + :align: center + +Not very interesting so far but it will be informative once we start to +vary the resolution later. + +If you want to view the mesh, you can copy the file +``base_mesh_vtk/staticFieldsOnCells.vtp`` over to your laptop or desktop and +open it in ParaView. See +`ParaView VTK Extractor `_ +for more details on the tool used to extract the mesh VTK file. + +.. image:: images/qu60_base_paraview.png + :width: 500 px + :align: center + + +The ``remap_topography`` step will produce: + +.. code-block:: bash + + $ cd ../remap_topography/ + $ ls + base_mesh.nc src_mesh.nc + dst_mesh.nc step.pickle + job_script.sh topography.nc + load_compass_env.sh topography_ncremap.nc + map_0.013x0.013degree_to_YAM10to60_conserve.nc topography_remapped.nc + mesh.cfg + +Here, the main result is ``topography_remapped.nc``, the ocean bathymetry and +Antarctic ice topography remapped to the mesh defined in ``base_mesh.nc``. + +Finally, the ``cull_mesh`` step will remove land cells from the mesh: + +The ``remap_topography`` step will produce: + +.. code-block:: bash + + $ cd ../cull_mesh/ + $ ls + base_mesh.nc job_script.sh + critical_blockages.geojson land_mask.nc + critical_blockages.nc land_mask_with_land_locked_cells.nc + critical_passages.geojson load_compass_env.sh + critical_passages_mask_final.nc mesh.cfg + critical_passages.nc step.pickle + culled_graph.info topography_culled.nc + culled_mesh.nc topography.nc + culled_mesh_vtk + +Here, the main output is ``culled_mesh.nc``. Similarly to the base mesh, you +can look at the the culled mesh in ParaVeiw by copying +``culled_mesh_vtk/staticFieldsOnCells.vtp`` to your laptop or desktop. + +.. image:: images/qu60_culled_paraview.png + :width: 500 px + :align: center + +Here, we have placed a white sphere inside the mesh so the land regions are +easier to see. After culling, the land just appears as holes in the mesh. + +.. _dev_tutorial_add_rrm_ec_base_mesh: + +Switching to an EC30to60 base resolution +---------------------------------------- + +Returning to the terminal where we are developing the code, let's make the mesh +more interesting. + +So far, all E3SM ocean and sea-ice RRMs start with the EC30to60 (eddy-closure +30 to 60 km) mesh as their base resolution. Let's do the same here. Starting +from the base of your development branch: + +.. code-block:: bash + + cd compass/ocean/tests/global_ocean/mesh/yam10to60 + vim __init__.py + +We will replace the constant 60-km mesh resolution with a latitude-dependent +function from the +`mesh_definition_tools `_ +module from MPAS-Tools. The default EC mesh has resolutions of 35 km at the +poles, 60 km at mid-latitudes and 30 km at the equator. + +.. code-block:: python + :emphasize-lines: 1, 17-18 + + import mpas_tools.mesh.creation.mesh_definition_tools as mdt + import numpy as np + + from compass.mesh import QuasiUniformSphericalMeshStep + + + class YAM10to60BaseMesh(QuasiUniformSphericalMeshStep): + + ... + + def build_cell_width_lat_lon(self): + + ... + + lat = np.linspace(-90., 90., nlat) + + cell_width_vs_lat = mdt.EC_CellWidthVsLat(lat) + cell_width = np.outer(cell_width_vs_lat, np.ones([1, lon.size])) + + return cell_width, lon, lat + +At this point, you can set up and test again like you did in +:ref:`dev_tutorial_add_rrm_test_mesh`, but this time you will want to use +a different work directory name, e.g.: + +.. code-block:: bash + + compass setup -n 254 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_ec + +Switch back to your other terminal to submit the job and look at the results. +The map of resolution in ``base_mesh/cellWidthGlobal.png`` should look like: + +.. image:: images/ec30to60.png + :width: 500 px + :align: center + +After culling, the mesh in ``culled_mesh/culled_mesh_vtk/staticFieldsOnCells.vtp`` +should look like: + +.. image:: images/ec30to60_culled_paraview.png + :width: 500 px + :align: center + +.. _dev_tutorial_add_rrm_add_high_res: + +Adding regions of higher resolution +----------------------------------- + +Now, let's add some regions of higher resolution to the mesh. + +We typically define these regions using `geojson `_ +files. The easiest way to create them is to go to `geojson.io `_. +There, you can find your way to the part of the globe you want to refine +and use the polygon tool to make a shape that will act as the boundary for your +high resolution region. + + +.. image:: images/geojson_io_south_atl.png + :width: 800 px + :align: center + +In my case, I have defined a region across the south Atlantic ocean with its +western side centered around the outlet of the Amazon river. My plan is to +define a region of moderately higher resolution across a fairly broad region +first, then define a region of higher resolution close to the Amazon delta +in a subsequent step. + +Let's make an actual ``geojson`` file with this contents. In your terminal for +editing code, from the root of the branch where we're developing: + +.. code-block:: bash + + cd compass/ocean/tests/global_ocean/mesh/yam10to60 + vim northern_south_atlantic.geojson + +Copy the contents of the json code on the right-hand side of the geojson.io +window and paste it into the file: + +.. code-block:: json + + { + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": {}, + "geometry": { + "coordinates": [ + [ + [ + -42.7022201869903, + 28.229943571814303 + ], + [ + -63.8408547092003, + 9.565520467643694 + ], + [ + -54.35184148160458, + -3.0088254981339873 + ], + [ + -37.52116934686214, + -8.341138860925426 + ], + [ + -12.947354056832182, + 10.997433207836309 + ], + [ + -11.493517385995887, + 27.701423680235493 + ], + [ + -42.7022201869903, + 28.229943571814303 + ] + ] + ], + "type": "Polygon" + } + } + ] + } + +Then, modify the ``properties`` dictionary similarly to this example: + +.. code-block:: + :emphasize-lines: 6-11 + + { + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "name": "Northern South Atlantic", + "component": "ocean", + "object": "region", + "author": "Xylar Asay-Davis" + }, + + ... + +These 4 fields are required for compass to be able to use the file. The +``name`` and ``author`` are entirely up to you and are intended to help +document the file in some useful way. The ``component`` must be ``"ocean"`` +and the ``object`` must be ``"region"``. + +Next, let's make the shape available in the code so we can use it later to make +a higher resolution region: + +.. code-block:: bash + + vim __init__.py + +.. code-block:: python + :emphasize-lines: 3, 11-20, 29 + + import mpas_tools.mesh.creation.mesh_definition_tools as mdt + import numpy as np + from geometric_features import read_feature_collection + + from compass.mesh import QuasiUniformSphericalMeshStep + + + class YAM10to60BaseMesh(QuasiUniformSphericalMeshStep): + """ + A step for creating YAM10to60 meshes + """ + def setup(self): + """ + Add some input files + """ + package = 'compass.ocean.tests.global_ocean.mesh.yam10to60' + self.add_input_file(filename='northern_south_atlantic.geojson', + package=package) + super().setup() + + def build_cell_width_lat_lon(self): + + ... + + cell_width_vs_lat = mdt.EC_CellWidthVsLat(lat) + cell_width = np.outer(cell_width_vs_lat, np.ones([1, lon.size])) + + # read the shape + fc = read_feature_collection('northern_south_atlantic.geojson') + + return cell_width, lon, lat + +In the ``setup()`` method above, we add the geojson file as an input to the +step that creates the base mesh. This is how compass finds the geojson file +when it's setting up the work directory where we will build the base mesh. + +In the ``build_cell_width_lat_lon()`` method, we read in a the geojson file +into a "feature collection" (``fc``) object that we will use below to define +the higher resolution region. + +Now, let's make further changes to the same file to use the shape to add a +higher resolution region: + +.. code-block:: bash + + vim __init__.py + +.. code-block:: python + :emphasize-lines: 4-7, 24-46 + + import mpas_tools.mesh.creation.mesh_definition_tools as mdt + import numpy as np + from geometric_features import read_feature_collection + from mpas_tools.cime.constants import constants + from mpas_tools.mesh.creation.signed_distance import ( + signed_distance_from_geojson, + ) + + from compass.mesh import QuasiUniformSphericalMeshStep + + + class YAM10to60BaseMesh(QuasiUniformSphericalMeshStep): + + def build_cell_width_lat_lon(self): + + ... + + cell_width_vs_lat = mdt.EC_CellWidthVsLat(lat) + cell_width = np.outer(cell_width_vs_lat, np.ones([1, lon.size])) + + # read the shape + fc = read_feature_collection('northern_south_atlantic.geojson') + + # How wide in meters the smooth transition between the background + # resolution and the finer resolution regions should be. + # 1200 km is equivalent to about 10 degrees latitude + trans_width = 1200e3 + + # The resolution in km of the finer resolution region + fine_cell_width = 20. + + # the radius of the earth defined in E3SM's shared constants + earth_radius = constants['SHR_CONST_REARTH'] + + # A field defined on the lat-long grid with the signed distance away + # from the boundary of the shape (positive outside and negative inside) + atlantic_signed_distance = signed_distance_from_geojson( + fc, lon, lat, earth_radius, max_length=0.25) + + # A field that goes smoothly from zero inside the shape to one outside + # the shape over the given transition width. + weights = 0.5 * (1 + np.tanh(atlantic_signed_distance / trans_width)) + + # The cell width in km becomes a blend of the background cell width + # and the finer cell width using the weights + cell_width = fine_cell_width * (1 - weights) + cell_width * weights + + return cell_width, lon, lat + +The function ``signed_distance_from_geojson()`` creates a functon on the +lat-lon grid that is the distance from any given point on the globe to the +boundary of the shape defined by the geojson file. The distance is positive +outside the shape and negative inside it. For better accuracy in computing the +distance, we subdivide the shape into segments of ``max_length=0.25`` degrees +latitude or longitude. We use the ``earth_radius`` defined in E3SM's shared +constants. + +Using the signed distance, we create a blending function ``weights`` that goes +from zero inside the shape smoothly to one outside the shape over a distance of +``trans_width`` meters. Then, we use the weights to blend from the fine +resolution inside the shape to the EC30to60 background resolution outside the +shape. + +There are also fancier ways to define gradients in resolution, for example +using the relative distance between the boundaries of 2 shapes. These are not +covered in the tutorial but you can get in touch with the E3SM Ocean Team to +discuss ways to define more complex maps of mesh resolution. + +Once, again, let's set up and run the mesh test case like we did in +:ref:`dev_tutorial_add_rrm_test_mesh`: + +.. code-block:: bash + + compass setup -n 254 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_alt20km + +As before, switch back to your other terminal to submit the job and look at the +results. + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_alt20km + sbatch job_script.custom.sh + tail -f compass.o* + +The map of resolution in ``base_mesh/cellWidthGlobal.png`` should look +like: + +.. image:: images/atl20km.png + :width: 500 px + :align: center + +After culling, the mesh in ``culled_mesh/culled_mesh_vtk/staticFieldsOnCells.vtp`` +should look like: + +.. image:: images/atl20km_culled_paraview.png + :width: 500 px + :align: center + +.. _dev_tutorial_add_rrm_add_very_high_res: + +Adding a very high resolution region +------------------------------------ + +Using the same approach as in the previous section, we can define another +region where we will increase the resolution to 10 km. + +I used geojson.io to create a region around the Amazon River delta: + +.. image:: images/geojson_io_amazon_delta.png + :width: 800 px + :align: center + +Then, I copied the code and pasted it into a file: + +.. code-block:: bash + + cd compass/ocean/tests/global_ocean/mesh/yam10to60 + vim amazon_delta.geojson + +I added the ``properties`` dictionary like in the previous example. + +.. code-block:: json + + { + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "name": "Amazon Delta", + "component": "ocean", + "object": "region", + "author": "Xylar Asay-Davis" + }, + "geometry": { + "coordinates": [ + [ + [ + -33.27493467565196, + 9.398029362516667 + ], + [ + -44.499833304073974, + 11.7502737267192 + ], + [ + -54.422618869265236, + 8.655607226691274 + ], + [ + -60.654712683354944, + 0.9780614705966428 + ], + [ + -54.56296235335806, + -9.767487562476404 + ], + [ + -41.34251704331987, + -9.500764493003032 + ], + [ + -36.85005485733731, + -3.655530642645047 + ], + [ + -33.03465151175149, + 4.644399816423899 + ], + [ + -33.27493467565196, + 9.398029362516667 + ] + ] + ], + "type": "Polygon" + } + } + ] + } + +Using this feature, I added a 10-km region: + +.. code-block:: bash + + vim __init__.py + +.. code-block:: python + :emphasize-lines: 11-12, 27-48 + + ... + + class YAM10to60BaseMesh(QuasiUniformSphericalMeshStep): + def setup(self): + """ + Add some input files + """ + package = 'compass.ocean.tests.global_ocean.mesh.yam10to60' + self.add_input_file(filename='northern_south_atlantic.geojson', + package=package) + self.add_input_file(filename='amazon_delta.geojson', + package=package) + super().setup() + + def build_cell_width_lat_lon(self): + + ... + + # A field that goes smoothly from zero inside the shape to one outside + # the shape over the given transition width. + weights = 0.5 * (1 + np.tanh(atlantic_signed_distance / trans_width)) + + # The cell width in km becomes a blend of the background cell width + # and the finer cell width using the weights + cell_width = fine_cell_width * (1 - weights) + cell_width * weights + + # read the shape + fc = read_feature_collection('amazon_delta.geojson') + + # 400 km is equivalent to about 3 degrees latitude + trans_width = 400e3 + + # The resolution in km of the finer resolution region + fine_cell_width = 10. + + # A field defined on the lat-long grid with the signed distance away + # from the boundary of the shape (positive outside and negative inside) + amazon_delta_signed_distance = signed_distance_from_geojson( + fc, lon, lat, earth_radius, max_length=0.25) + + # A field that goes smoothly from zero inside the shape to one outside + # the shape over the given transition width. + weights = 0.5 * (1 + np.tanh( + amazon_delta_signed_distance / trans_width)) + + # The cell width in km becomes a blend of the background cell width + # and the finer cell width using the weights + cell_width = fine_cell_width * (1 - weights) + cell_width * weights + + return cell_width, lon, lat + +Same procedure as before, set up the test case: + +.. code-block:: bash + + compass setup -n 254 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +Switch back to your other terminal to submit the job and look at the results. + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + sbatch job_script.custom.sh + tail -f compass.o* + +The map of resolution in ``base_mesh/cellWidthGlobal.png`` should look +like: + +.. image:: images/yam10to60.png + :width: 500 px + :align: center + +After culling, the mesh in ``culled_mesh/culled_mesh_vtk/staticFieldsOnCells.vtp`` +should look like: + +.. image:: images/yam10to60_culled_paraview.png + :width: 500 px + :align: center + +.. _dev_tutorial_add_rrm_add_init: + +Adding an initial condition and performance test +------------------------------------------------ + +We now have a horizontal ocean mesh, but there are several more steps in +compass before we can start to integrate the new mesh into E3SM. + +First, we need to add an ``init`` test case that will create the vertical +coordinate and the initial condition. Most of what we need to define for the +initial condition is set with config options: + +.. code-block:: bash + + cd compass/ocean/tests/global_ocean/mesh/yam10to60 + vim yam10to60.cfg + +.. code-block:: ini + + # Options related to the vertical grid + [vertical_grid] + + # the type of vertical grid + grid_type = index_tanh_dz + + # Number of vertical levels + vert_levels = 64 + + # Depth of the bottom of the ocean + bottom_depth = 5500.0 + + # The minimum layer thickness + min_layer_thickness = 10.0 + + # The maximum layer thickness + max_layer_thickness = 250.0 + + # The characteristic number of levels over which the transition between + # the min and max occurs + transition_levels = 28 + +The standard E3SM v3 vertical grid is defined with these parameters. It +transitions from 10 m resolution at the surface to 250 m resolution at a depth +of 5500 m over 64 vertical levels. The transition is centered around the 28th +vertical level. You can modify these parameters or use a different vertical +coordinate (at your own risk). + +Next, we add a section with some required config options including some +metadata: + +.. code-block:: bash + + vim yam10to60.cfg + +.. code-block:: ini + :emphasize-lines: 8-43 + + ... + + # The characteristic number of levels over which the transition between + # the min and max occurs + transition_levels = 28 + + + # options for global ocean testcases + [global_ocean] + + # the approximate number of cells in the mesh + approx_cell_count = 270000 + + ## metadata related to the mesh + # the prefix (e.g. QU, EC, WC, SO) + prefix = YAM + # a description of the mesh and initial condition + mesh_description = Yet Another Mesh (YAM) is regionally refined around the Amazon + River delta. It is used in E3SM version ${e3sm_version} for + studies of blah, blah. It has Enhanced resolution (${min_res} km) + around the Amazon delta, 20-km resolution in the northern South + Atlantic, 35-km resolution at the poles, 60-km at mid + latitudes, and 30-km at the equator. This mesh has <<>> + vertical levels. + # E3SM version that the mesh is intended for + e3sm_version = 3 + # The revision number of the mesh, which should be incremented each time the + # mesh is revised + mesh_revision = 1 + # the minimum (finest) resolution in the mesh + min_res = 10 + # the maximum (coarsest) resolution in the mesh, can be the same as min_res + max_res = 60 + # The URL of the pull request documenting the creation of the mesh + pull_request = https://github.com/MPAS-Dev/compass/pull/XXX + + + # config options related to initial condition and diagnostics support files + # for E3SM + [files_for_e3sm] + + # CMIP6 grid resolution + cmip6_grid_res = 180x360 + +The ``approx_cell_count`` is something you can only determine after you've +made the mesh. In your terminal where you've been running tests: + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + cd ocean/global_ocean/YAM10to60/mesh/cull_mesh + source load_compass_env.sh + ncdump -h culled_mesh.nc | more + +This will tell you the numer of cells (``nCells``). Round this to 2 +significant digits and put this in ``approx_cell_count``. This will be used +to determine an appropriate number of MPI tasks and nodes needed to perform +forward runs with this mesh. + +The ``prefix`` should match the beginning fo the mesh name we have been using +all along, ``YAM`` in this case. + +The ``description`` is used in the metadata of files produced by compass for +this mesh. It should be a fairly detailed description of how resolution is +distributed and what the intended purpose of the mesh is. + +The ``e3sm_version`` is what E3SM version the mesh is intended to be used in. +Presumably, this is 3 for the time being, since no new meshes will be added +to v2 and we don't know much about v4 yet. + +The ``mesh_revision`` should be 1 to begin with but should be incremented to +give each version of the mesh a unique number. Typically, it is time to +increase the revision number if you are altering the mesh and the current +revision number has already been used in a simulation that might be published +or put to some other broader scientific or technical use. If you are still at +the stage of tinkering with the mesh, it may be preferable not to increment the +revision number with each modification. + +``min_res`` *must* be the minimum resolution of the mesh in km, and ``max_res`` +should be the maximum. (``min_res`` is used to determine the time step for +forward runs, so that is why it is required to be correct; ``max_res`` is only +used for metadata and in the mesh name.) + +``pull_request`` points to a pull request (if there is one) where the mesh +was added to compass. This is a useful piece of metadata so folks know where +to look for a discussion of the mesh. + +``cmip6_grid_res`` is the CMIP6 (and presumably 7) resolution to which data +will be remapped for publication. Typically, this is ``180x360`` (i.e. a +1 degree grid) for RRMs because a lot of space is otherwise wasted on +coarse-resolution regions. + +We also need to override some default namelist options with values appropriate +for the mesh. Many of these choices will depend on the details of the mesh +you are making. However, here are some common ones: + +.. code-block:: bash + + vim namelist.split_explicit + +.. code-block:: + + config_time_integrator = 'split_explicit' + config_run_duration = '0000_01:00:00' + config_use_mom_del2 = .true. + config_mom_del2 = 10.0 + config_use_mom_del4 = .true. + config_mom_del4 = 1.5e10 + config_hmix_scaleWithMesh = .true. + config_use_GM = .true. + config_GM_closure = 'constant' + config_GM_constant_kappa = 600.0 + +The ``config_run_duration`` is the length of a performance test run, and should +only be a few time steps. + +The ``config_mom_del2`` and ``config_mom_del4`` are the eddy viscosity and +eddy hyperviscosity, and typically scale in proportion to the cell size and +the cell size cubed, respectively. These are appropriate values for a minimum +resolution of 10 km as in this example. We scale these coefficients with the +cell resolution when ``config_hmix_scaleWithMesh = .true.``. + +The GM coefficients can probably be left as they are here for a typical RRM. + +Next, we will add the ``init`` and ``performance_tests`` test cases for the +new mesh: + +.. code-block:: bash + + cd ../.. + vim __init__.py + +.. code-block:: python + :emphasize-lines: 16-24 + + ... + + class GlobalOcean(TestGroup): + + ... + + def __init__(self, mpas_core): + + ... + + for mesh_name in ['YAM10to60', 'YAMwISC10to60']: + mesh_test = Mesh(test_group=self, mesh_name=mesh_name, + remap_topography=True) + self.add_test_case(mesh_test) + + init_test = Init(test_group=self, mesh=mesh_test, + initial_condition='WOA23', + with_bgc=False) + self.add_test_case(init_test) + + self.add_test_case( + PerformanceTest( + test_group=self, mesh=mesh_test, init=init_test, + time_integrator='split_explicit')) + + # A test case for making E3SM support files from an existing mesh + self.add_test_case(FilesForE3SM(test_group=self)) + +We have indicated that we want an initial condition interpolated from the +World Ocean Atlas 2023 (WOA23) data set, that we do not want to include +BGC tracers, and that we want to use the split-explicit time integrator +(the E3SM default) in our performance test. + +Let's see if the test cases show up: + +.. code-block:: bash + + compass list | grep YAM + +You should see something like: + +.. code-block:: + + 254: ocean/global_ocean/YAM10to60/mesh + 255: ocean/global_ocean/YAM10to60/WOA23/init + 256: ocean/global_ocean/YAM10to60/WOA23/performance_test + 257: ocean/global_ocean/YAMwISC10to60/mesh + 258: ocean/global_ocean/YAMwISC10to60/WOA23/init + 259: ocean/global_ocean/YAMwISC10to60/WOA23/performance_test + +Okay, everything looks good. Let's set up and run the 2 remaining tests: + +.. code-block:: bash + + compass setup -n 255 256 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +We can save the time of rerunning the ``mesh`` test by setting up in the same +work directory as our final mesh run with 10-km finest resolution. + +Switch back to your other terminal to alter the job script and submit the job. + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + sbatch job_script.custom.sh + tail -f compass.o* + +You should see something a lot like this: + +.. code-block:: + + Loading conda environment + Done. + + Loading Spack environment... + Done. + + ocean/global_ocean/YAM10to60/WOA23/init + * step: initial_state + test execution: SUCCESS + test runtime: 01:07 + ocean/global_ocean/YAM10to60/WOA23/performance_test + * step: forward + test execution: SUCCESS + test runtime: 01:35 + Test Runtimes: + 01:07 PASS ocean_global_ocean_YAM10to60_WOA23_init + 01:35 PASS ocean_global_ocean_YAM10to60_WOA23_performance_test + Total runtime 02:42 + PASS: All passed successfully! + +If these tests aren't successful, you'll probably need some expert help from +the E3SM Ocean Team, but you can take a look at the log files and see if you +can diagnose any issues yourself. + +.. _dev_tutorial_add_rrm_add_dyn_adj: + +Adding a dynamic adjustment test +-------------------------------- + +The initial condition we generated in the last section starts with the ocean +at rest. This state is not consistent with the density profile, so there +will be a period of a few months of rapid adjustment involving dissipation of +energetic wave dissipation and acceleration of currents. We call this period +"dynamic adjustment" because the term "spin up" is reserved in Earth system +modeling for reaching a quasi-equilibrium over many centuries. + +Dynamic adjustment is something of an art, and requires some trial and error +in many cases. The basic idea is that we begin with a shorter time step than +we hope to be able to use in production simulations and also with some fairly +strong momentum dissipation. We run forward in time, monitoring the CFL +number and the kinetic energy, which can each tell us if things are going +awry. After several days of simulation, waves have hopefully dissipated to +the point that we can increase the time step and/or decrease the level of +damping. + +We need to create a new subdirectory and add a new class for the dynamic +adjustment test case. From the branch root: + +.. code-block:: bash + + cd compass/ocean/tests/global_ocean/mesh/yam10to60 + mkdir dynamic_adjustment + cd dynamic_adjustment + vim __init__.py + +.. code-block:: python + + from compass.ocean.tests.global_ocean.dynamic_adjustment import ( + DynamicAdjustment, + ) + from compass.ocean.tests.global_ocean.forward import ForwardStep + + + class YAM10to60DynamicAdjustment(DynamicAdjustment): + """ + A test case performing dynamic adjustment (dissipating fast-moving waves) + from an initial condition on the YAM10to60 MPAS-Ocean mesh + + Attributes + ---------- + restart_filenames : list of str + A list of restart files from each dynamic-adjustment step + """ + + def __init__(self, test_group, mesh, init, time_integrator): + """ + Create the test case + + Parameters + ---------- + test_group : compass.ocean.tests.global_ocean.GlobalOcean + The global ocean test group that this test case belongs to + + mesh : compass.ocean.tests.global_ocean.mesh.Mesh + The test case that produces the mesh for this run + + init : compass.ocean.tests.global_ocean.init.Init + The test case that produces the initial condition for this run + + time_integrator : {'split_explicit', 'RK4'} + The time integrator to use for the forward run + """ + if time_integrator != 'split_explicit': + raise ValueError(f'{mesh.mesh_name} dynamic adjustment not ' + f'defined for {time_integrator}') + + restart_times = ['0001-01-03_00:00:00'] + restart_filenames = [ + f'restarts/rst.{restart_time.replace(":", ".")}.nc' + for restart_time in restart_times] + + super().__init__(test_group=test_group, mesh=mesh, init=init, + time_integrator=time_integrator, + restart_filenames=restart_filenames) + + module = self.__module__ + + shared_options = \ + {'config_AM_globalStats_enable': '.true.', + 'config_AM_globalStats_compute_on_startup': '.true.', + 'config_AM_globalStats_write_on_startup': '.true.', + 'config_use_activeTracers_surface_restoring': '.true.'} + + # first step + step_name = 'damped_adjustment_1' + step = ForwardStep(test_case=self, mesh=mesh, init=init, + time_integrator=time_integrator, name=step_name, + subdir=step_name, get_dt_from_min_res=False) + + namelist_options = { + 'config_run_duration': "'00-00-02_00:00:00'", + 'config_dt': "'00:03:00'", + 'config_btr_dt': "'00:00:06'", + 'config_implicit_bottom_drag_type': "'constant_and_rayleigh'", + 'config_Rayleigh_damping_coeff': '1.0e-4'} + namelist_options.update(shared_options) + step.add_namelist_options(namelist_options) + + stream_replacements = { + 'output_interval': '00-00-10_00:00:00', + 'restart_interval': '00-00-02_00:00:00'} + step.add_streams_file(module, 'streams.template', + template_replacements=stream_replacements) + + step.add_output_file(filename=f'../{restart_filenames[0]}') + self.add_step(step) + +This sets up one step called ``damped_adjustment_1`` that runs for 2 days +with 3-minute time steps (we hope to run with a 5 or 6 minute time steps once +we're fully adjusted, given the 10-km minimum resolution), 6-second +subcycling (barotropic or ``btr``) time step, and a strong Rayleigh damping +of 1e-4. Since we're running for 2 days, we set the restart interval to 2 +days. + +You should be able to pick a time step in proprotion to the highest resolution +in your mesh. If a 3-minute time steps works well for a 10 km mesh (and that's +a little on the safe side), the time step for a 1 km mesh would need to be on +the order of 20 seconds. The barotropic time step needs to be 20-30 times +shorter than the full time step. + +We have enabled global stats (``config_AM_globalStats_enable = .true.``) so +we can monitor the progress more easily. + +We have also set up a set of streams for writing out data as we go. The +``streams.template`` file that we will modify looks something like this: + +.. code-block:: bash + + vim streams.template + +.. code-block:: xml + + + + + + + + + + +The ``output_interval`` and ``restart_interval`` will get replaced with +different values in different steps as we add them. + +We need to add the dynamic adjustment test case to the ``global_ocean`` test +group: + +.. code-block:: bash + + cd ../../../ + vim __init__.py + +.. code-block:: python + :emphasize-lines: 6-8, 38-42 + + ... + + from compass.ocean.tests.global_ocean.mesh.wc14.dynamic_adjustment import ( + WC14DynamicAdjustment, + ) + from compass.ocean.tests.global_ocean.mesh.yam10to60.dynamic_adjustment import ( # noqa: E501 + YAM10to60DynamicAdjustment, + ) + from compass.ocean.tests.global_ocean.monthly_output_test import ( + MonthlyOutputTest, + ) + + ... + + class GlobalOcean(TestGroup): + + ... + + def __init__(self, mpas_core): + + ... + + for mesh_name in ['YAM10to60', 'YAMwISC10to60']: + mesh_test = Mesh(test_group=self, mesh_name=mesh_name, + remap_topography=True) + self.add_test_case(mesh_test) + + init_test = Init(test_group=self, mesh=mesh_test, + initial_condition='WOA23', + with_bgc=False) + self.add_test_case(init_test) + + self.add_test_case( + PerformanceTest( + test_group=self, mesh=mesh_test, init=init_test, + time_integrator='split_explicit')) + + dynamic_adjustment_test = YAM10to60DynamicAdjustment( + test_group=self, mesh=mesh_test, init=init_test, + time_integrator='split_explicit') + self.add_test_case(dynamic_adjustment_test) + + # A test case for making E3SM support files from an existing mesh + self.add_test_case(FilesForE3SM(test_group=self)) + + ... + +Let's see if the test cases show up: + +.. code-block:: bash + + compass list | grep YAM + +You should see something like: + +.. code-block:: + + 254: ocean/global_ocean/YAM10to60/mesh + 255: ocean/global_ocean/YAM10to60/WOA23/init + 256: ocean/global_ocean/YAM10to60/WOA23/performance_test + 257: ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment + 258: ocean/global_ocean/YAMwISC10to60/mesh + 259: ocean/global_ocean/YAMwISC10to60/WOA23/init + 260: ocean/global_ocean/YAMwISC10to60/WOA23/performance_test + 261: ocean/global_ocean/YAMwISC10to60/WOA23/dynamic_adjustment + +Okay, everything looks good. Let's set up and run the ``dynamic_adjustment`` test: + +.. code-block:: bash + + compass setup -n 257 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +Switch back to your other terminal to submit the job. + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + sbatch job_script.custom.sh + tail -f compass.o* + +This time, the output should look like: + +.. code-block:: + + Loading conda environment + Done. + + Loading Spack environment... + Done. + + ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment + * step: damped_adjustment_1 + test execution: ERROR + see: case_outputs/ocean_global_ocean_YAM10to60_WOA23_dynamic_adjustment.log + test runtime: 06:03 + Test Runtimes: + 06:03 FAIL ocean_global_ocean_YAM10to60_WOA23_dynamic_adjustment + Total runtime 06:03 + FAIL: 1 test failed, see above. + +This error isn't a big deal. It's related to the fact that we're not done +implementing the test case and it's expecting a step called ``simulation`` that +we haven't added yet: + +.. code-block:: + + $ tail -n 20 case_outputs/ocean_global_ocean_YAM10to60_WOA23_dynamic_adjustment.log + + Running: srun -c 1 -N 6 -n 768 ./ocean_model -n namelist.ocean -s streams.ocean + + compass calling: compass.ocean.tests.global_ocean.mesh.yam10to60.dynamic_adjustment.YAM10to60DynamicAdjustment.validate() + inherited from: compass.ocean.tests.global_ocean.dynamic_adjustment.DynamicAdjustment.validate() + in /gpfs/fs1/home/ac.xylar/compass/add-rrm-tutorial/compass/ocean/tests/global_ocean/dynamic_adjustment.py + + Exception raised in the test case's validate() method + Traceback (most recent call last): + File "/gpfs/fs1/home/ac.xylar/compass/add-rrm-tutorial/compass/run/serial.py", line 335, in _log_and_run_test + test_case.validate() + File "/gpfs/fs1/home/ac.xylar/compass/add-rrm-tutorial/compass/ocean/tests/global_ocean/dynamic_adjustment.py", line 63, in validate + compare_variables(test_case=self, variables=variables, + File "/gpfs/fs1/home/ac.xylar/compass/add-rrm-tutorial/compass/validate.py", line 94, in compare_variables + raise ValueError('{} does not appear to be an output of any step ' + ValueError: simulation/output.nc does not appear to be an output of any step in this test case. + +You can also monitor the result by looking at the global statistics: + +.. code-block:: bash + + cd ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment/damped_adjustment_1 + source load_compass_env.sh + ncdump -v kineticEnergyCellMax analysis_members/globalStats*.nc + ncdump -v CFLNumberGlobal analysis_members/globalStats*.nc + +The kinetic energy should increase gradually (and then likely decrease because +of the damping) but shouldn't spike up during a damped adjustment step. The +CFL number ideally shouldn't exceed about 0.1 during damped adjustment, though +we tend to push it a bit higher during the simulation phase to see how large +we can make it while maintaining stability. + +You can also monitor the MPAS-Ocean progress (e.g. the time stepping) with: + +.. code-block:: bash + + tail log.ocean.0000.out + +If the ``damped_adjustment_1`` step is successful, it's time to add more steps +in which we will ramp down damping and then increase the time step. Let's add a +second step that runs longer (8 days) with less damping: + +.. code-block:: bash + + vim __init__.py + +.. code-block:: python + :emphasize-lines: 11, 18-43 + + ... + + class YAM10to60DynamicAdjustment(DynamicAdjustment): + + ... + + def __init__(self, test_group, mesh, init, time_integrator): + + ... + + restart_times = ['0001-01-03_00:00:00', '0001-01-11_00:00:00'] + + ... + + step.add_output_file(filename=f'../{restart_filenames[0]}') + self.add_step(step) + + # second step + step_name = 'damped_adjustment_2' + step = ForwardStep(test_case=self, mesh=mesh, init=init, + time_integrator=time_integrator, name=step_name, + subdir=step_name, get_dt_from_min_res=False) + + namelist_options = { + 'config_run_duration': "'00-00-08_00:00:00'", + 'config_dt': "'00:03:00'", + 'config_btr_dt': "'00:00:05'", + 'config_implicit_bottom_drag_type': "'constant_and_rayleigh'", + 'config_Rayleigh_damping_coeff': '1.0e-5', + 'config_do_restart': '.true.', + 'config_start_time': f"'{restart_times[0]}'"} + namelist_options.update(shared_options) + step.add_namelist_options(namelist_options) + + stream_replacements = { + 'output_interval': '00-00-10_00:00:00', + 'restart_interval': '00-00-02_00:00:00'} + step.add_streams_file(module, 'streams.template', + template_replacements=stream_replacements) + + step.add_input_file(filename=f'../{restart_filenames[0]}') + step.add_output_file(filename=f'../{restart_filenames[1]}') + self.add_step(step) + +You can set up again and test the second step. In your coding terminal: + +.. code-block:: bash + + compass setup -n 257 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +Back in your terminal in the work directory: + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + cd ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment/damped_adjustment_2 + sbatch job_script.sh + tail -f compass.o* + +Again, you will get errors about the missing ``simulation/output.nc`` file, +but don't worry about those. + +If that goes okay, let's add a third step that runs for 10 days with even less +damping. We can also write out less frequent restarts (every 10 days). In +the coding terminal, which should still be in the ``dynamic_adjustment`` +subdirectory: + +.. code-block:: bash + + vim __init__.py + +.. code-block:: python + :emphasize-lines: 11-12, 19-43 + + ... + + class YAM10to60DynamicAdjustment(DynamicAdjustment): + + ... + + def __init__(self, test_group, mesh, init, time_integrator): + + ... + + restart_times = ['0001-01-03_00:00:00', '0001-01-11_00:00:00', + '0001-01-21_00:00:00'] + + ... + + step.add_input_file(filename=f'../{restart_filenames[0]}') + step.add_output_file(filename=f'../{restart_filenames[1]}') + self.add_step(step) + + # third step + step_name = 'damped_adjustment_3' + step = ForwardStep(test_case=self, mesh=mesh, init=init, + time_integrator=time_integrator, name=step_name, + subdir=step_name, get_dt_from_min_res=False) + + namelist_options = { + 'config_run_duration': "'00-00-10_00:00:00'", + 'config_dt': "'00:03:00'", + 'config_btr_dt': "'00:00:06'", + 'config_implicit_bottom_drag_type': "'constant_and_rayleigh'", + 'config_Rayleigh_damping_coeff': '1.0e-6', + 'config_do_restart': '.true.', + 'config_start_time': f"'{restart_times[1]}'"} + namelist_options.update(shared_options) + step.add_namelist_options(namelist_options) + + stream_replacements = { + 'output_interval': '00-00-10_00:00:00', + 'restart_interval': '00-00-10_00:00:00'} + step.add_streams_file(module, 'streams.template', + template_replacements=stream_replacements) + + step.add_input_file(filename=f'../{restart_filenames[1]}') + step.add_output_file(filename=f'../{restart_filenames[2]}') + self.add_step(step) + +Set up again in the coding terminal: + +.. code-block:: bash + + compass setup -n 257 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +And run this step in the work-directory terminal: + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + cd ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment/damped_adjustment_3 + sbatch job_script.sh + tail -f compass.o* + +Now, we add a fourth that runs for 20 days without any damping, back in the +coding terminal: + +.. code-block:: bash + + vim __init__.py + +.. code-block:: python + :emphasize-lines: 11-12, 19-43 + + ... + + class YAM10to60DynamicAdjustment(DynamicAdjustment): + + ... + + def __init__(self, test_group, mesh, init, time_integrator): + + ... + + restart_times = ['0001-01-03_00:00:00', '0001-01-11_00:00:00', + '0001-01-21_00:00:00', '0001-02-10_00:00:00'] + + ... + + step.add_input_file(filename=f'../{restart_filenames[1]}') + step.add_output_file(filename=f'../{restart_filenames[2]}') + self.add_step(step) + + # fourth step + step_name = 'damped_adjustment_4' + step = ForwardStep(test_case=self, mesh=mesh, init=init, + time_integrator=time_integrator, name=step_name, + subdir=step_name, get_dt_from_min_res=False) + + namelist_options = { + 'config_run_duration': "'00-00-20_00:00:00'", + 'config_dt': "'00:03:00'", + 'config_btr_dt': "'00:00:06'", + 'config_do_restart': '.true.', + 'config_start_time': f"'{restart_times[2]}'"} + namelist_options.update(shared_options) + step.add_namelist_options(namelist_options) + + stream_replacements = { + 'output_interval': '00-00-10_00:00:00', + 'restart_interval': '00-00-10_00:00:00'} + step.add_streams_file(module, 'streams.template', + template_replacements=stream_replacements) + + step.add_input_file(filename=f'../{restart_filenames[2]}') + step.add_output_file(filename=f'../{restart_filenames[3]}') + self.add_step(step) + + +Set up again in the coding terminal: + +.. code-block:: bash + + compass setup -n 257 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +And run this step in the work-directory terminal: + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + cd ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment/damped_adjustment_4 + sbatch job_script.sh + tail -f compass.o* + +Finally, we add one more step where we run for 10 more days with a longer +time step: + +.. code-block:: bash + + vim __init__.py + +.. code-block:: python + :emphasize-lines: 11-13, 20-45 + + ... + + class YAM10to60DynamicAdjustment(DynamicAdjustment): + + ... + + def __init__(self, test_group, mesh, init, time_integrator): + + ... + + restart_times = ['0001-01-03_00:00:00', '0001-01-11_00:00:00', + '0001-01-21_00:00:00', '0001-02-10_00:00:00', + '0001-02-20_00:00:00'] + + ... + + step.add_input_file(filename=f'../{restart_filenames[2]}') + step.add_output_file(filename=f'../{restart_filenames[3]}') + self.add_step(step) + + # final step + step_name = 'simulation' + step = ForwardStep(test_case=self, mesh=mesh, init=init, + time_integrator=time_integrator, name=step_name, + subdir=step_name, get_dt_from_min_res=False) + + namelist_options = { + 'config_run_duration': "'00-00-10_00:00:00'", + 'config_dt': "'00:08:00'", + 'config_btr_dt': "'00:00:15'", + 'config_do_restart': '.true.', + 'config_start_time': f"'{restart_times[3]}'"} + namelist_options.update(shared_options) + step.add_namelist_options(namelist_options) + + stream_replacements = { + 'output_interval': '00-00-10_00:00:00', + 'restart_interval': '00-00-10_00:00:00'} + step.add_streams_file(module, 'streams.template', + template_replacements=stream_replacements) + + step.add_input_file(filename=f'../{restart_filenames[3]}') + step.add_output_file(filename=f'../{restart_filenames[4]}') + step.add_output_file(filename='output.nc') + self.add_step(step) + + self.restart_filenames = restart_filenames + +Set up again in the coding terminal: + +.. code-block:: bash + + compass setup -n 257 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +And run this step in the work-directory terminal: + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + cd ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment/simulation + sbatch job_script.sh + tail -f compass.o* + +The art of this process goes into how you choose to adjust the time step, +duration of each of these runs, and the amount of damping. You may add more +steps or remove some if 5 doesn't work well for your mesh. Make sure that +the restart file that is an output of the previous step is the input to the +next one. + +It's also a good idea to take a look at a few output fields in ParaVeiw. First, +run: + +.. code-block:: bash + + paraview_vtk_field_extractor.py -m init.nc \ + -f output.nc -v normalVelocity,temperature,salinity \ + -d nVertLevels=0 + +Then, transfer the contents of ``vtk_files`` to your laptop or desktop and open +the file ``fieldsOnCells.pvd`` or ``fieldsOnEdges.pvd`` in ParaView. + +.. image:: images/simulation_sst.png + :width: 500 px + :align: center + +Note that ``normalVelocity`` is a pretty noisy field because its orientation +differs by on the order of 60 degrees between adjacent edges. It helps to take +the absolute value as shown below, in which case it it can indicate where +strong currents are present. + +.. image:: images/simulation_norm_vel_mag.png + :width: 500 px + :align: center + +Alternatively, you can add fields like ``kineticEnergyCell`` to the ``output`` +stream in ``streams.ocean`` before running the step, and then you will have +this field available to visualize. + + +.. _dev_tutorial_add_rrm_add_files_for_e3sm: + +Adding a files for E3SM test +---------------------------- + +The final test case to add for a new RRM mesh is ``files_for_e3sm``. This test +case creates a number of important files in the format they are needed for +E3SM or diagnostic software used to analysis E3SM simulations. No additional +customization should be needed for this mesh beyond the config options we +already set up earlier in the tutorial. We just need to add the test case +itself for this mesh to the ``global_ocean`` test group. + +Starting from the root of our development branch: + +.. code-block:: bash + + cd compass/ocean/tests/global_ocean + vim __init__.py + +.. code-block:: python + :emphasize-lines: 15-16 + + ... + + class GlobalOcean(TestGroup): + + ... + + def __init__(self, mpas_core): + + ... + + # Kuroshio meshes without ice-shelf cavities + self._add_tests(mesh_names=['Kuroshio12to60', 'Kuroshio8to60'], + DynamicAdjustment=KuroshioDynamicAdjustment) + + self._add_tests(mesh_names=['YAM10to60', 'YAMwISC10to60'], + DynamicAdjustment=YAM10to60DynamicAdjustment) + + # A test case for making E3SM support files from an existing mesh + self.add_test_case(FilesForE3SM(test_group=self)) + + ... + +We delete the whole ``for`` loop over ``mesh_name`` and instead take advantage +of the fact that the ``_add_tests()`` method of ``GlobalOcean`` will +add the 5 test cases we want by default. (We added them manually, one by one +before so we could test them one or two at a time.) + +Now, when you list the tests, you should see: + +.. code-block:: + + $ compass list | grep YAM + 254: ocean/global_ocean/YAM10to60/mesh + 255: ocean/global_ocean/YAM10to60/WOA23/init + 256: ocean/global_ocean/YAM10to60/WOA23/performance_test + 257: ocean/global_ocean/YAM10to60/WOA23/dynamic_adjustment + 258: ocean/global_ocean/YAM10to60/WOA23/files_for_e3sm + 259: ocean/global_ocean/YAMwISC10to60/mesh + 260: ocean/global_ocean/YAMwISC10to60/WOA23/init + 261: ocean/global_ocean/YAMwISC10to60/WOA23/performance_test + 262: ocean/global_ocean/YAMwISC10to60/WOA23/dynamic_adjustment + 263: ocean/global_ocean/YAMwISC10to60/WOA23/files_for_e3sm + +Once again, you can take advantage of the test cases you've already run, +setting up just the new ``files_for_e3sm`` test from the coding terminal: + +.. code-block:: bash + + compass setup -n 258 \ + -p E3SM-Project/components/mpas-ocean/ \ + -w /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + +And run this test case in the work-directory terminal: + +.. code-block:: bash + + cd /lcrc/group/e3sm/${USER}/compass_tests/tests_20230527/yam10to60_final + cd ocean/global_ocean/YAM10to60/WOA23/files_for_e3sm + sbatch job_script.sh + tail -f compass.o* + +If this goes well, you are ready to ask for help from members of the E3SM +Ocean Team to add support for your mesh to E3SM itself, including uploading +the files produced by the ``files_for_e3sm`` test case to the E3SM +``inputdata`` and ``diagnostics`` directories on our data server. This is +beyond the scope of this tutorial and is not typically something a non-expert +can take on on their own. diff --git a/docs/tutorials/images/atl20km.png b/docs/tutorials/images/atl20km.png new file mode 100644 index 0000000000..f2dd11437b Binary files /dev/null and b/docs/tutorials/images/atl20km.png differ diff --git a/docs/tutorials/images/atl20km_culled_paraview.png b/docs/tutorials/images/atl20km_culled_paraview.png new file mode 100644 index 0000000000..9e9f96a9bd Binary files /dev/null and b/docs/tutorials/images/atl20km_culled_paraview.png differ diff --git a/docs/tutorials/images/ec30to60.png b/docs/tutorials/images/ec30to60.png new file mode 100644 index 0000000000..316c307f55 Binary files /dev/null and b/docs/tutorials/images/ec30to60.png differ diff --git a/docs/tutorials/images/ec30to60_culled_paraview.png b/docs/tutorials/images/ec30to60_culled_paraview.png new file mode 100644 index 0000000000..b20150415d Binary files /dev/null and b/docs/tutorials/images/ec30to60_culled_paraview.png differ diff --git a/docs/tutorials/images/geojson_io_amazon_delta.png b/docs/tutorials/images/geojson_io_amazon_delta.png new file mode 100644 index 0000000000..220641cad0 Binary files /dev/null and b/docs/tutorials/images/geojson_io_amazon_delta.png differ diff --git a/docs/tutorials/images/geojson_io_south_atl.png b/docs/tutorials/images/geojson_io_south_atl.png new file mode 100644 index 0000000000..84150b0cff Binary files /dev/null and b/docs/tutorials/images/geojson_io_south_atl.png differ diff --git a/docs/tutorials/images/qu60.png b/docs/tutorials/images/qu60.png new file mode 100644 index 0000000000..ffe85185da Binary files /dev/null and b/docs/tutorials/images/qu60.png differ diff --git a/docs/tutorials/images/qu60_base_paraview.png b/docs/tutorials/images/qu60_base_paraview.png new file mode 100644 index 0000000000..e8313e9a26 Binary files /dev/null and b/docs/tutorials/images/qu60_base_paraview.png differ diff --git a/docs/tutorials/images/qu60_culled_paraview.png b/docs/tutorials/images/qu60_culled_paraview.png new file mode 100644 index 0000000000..08a4fc6625 Binary files /dev/null and b/docs/tutorials/images/qu60_culled_paraview.png differ diff --git a/docs/tutorials/images/simulation_norm_vel_mag.png b/docs/tutorials/images/simulation_norm_vel_mag.png new file mode 100644 index 0000000000..d10f853a9c Binary files /dev/null and b/docs/tutorials/images/simulation_norm_vel_mag.png differ diff --git a/docs/tutorials/images/simulation_sst.png b/docs/tutorials/images/simulation_sst.png new file mode 100644 index 0000000000..0e23ed0c93 Binary files /dev/null and b/docs/tutorials/images/simulation_sst.png differ diff --git a/docs/tutorials/images/yam10to60.png b/docs/tutorials/images/yam10to60.png new file mode 100644 index 0000000000..f726b4bb78 Binary files /dev/null and b/docs/tutorials/images/yam10to60.png differ diff --git a/docs/tutorials/images/yam10to60_culled_paraview.png b/docs/tutorials/images/yam10to60_culled_paraview.png new file mode 100644 index 0000000000..b63669ca27 Binary files /dev/null and b/docs/tutorials/images/yam10to60_culled_paraview.png differ