Merge pull request #771 from jeremy-lilly/fb-lts

FB-LTS Hurricane Test Case

compass/ocean/tests/hurricane/__init__.py

@@ -19,7 +19,7 @@ def __init__(self, mpas_core):
         storm = 'sandy'
         mesh_name = 'DEQU120at30cr10rr2'
 
-        for use_lts in [False, True]:
+        for use_lts in [False, 'LTS', 'FB_LTS']:
 
             mesh = Mesh(test_group=self,
                         mesh_name=mesh_name,

compass/ocean/tests/hurricane/fblts/forward/namelist.ocean

@@ -0,0 +1,39 @@
+config_start_time = '2012-10-10_00:00:00'
+config_stop_time = '2012-11-03_00:00:00'
+config_run_duration = 'none'
+config_dt = '72'
+config_time_integrator = 'FB_LTS'
+config_dt_scaling_LTS = 4
+config_number_of_time_levels = 5
+config_hmix_scaleWithMesh = .true.
+config_use_mom_del4 = .true.
+config_mom_del4 = 4.0e8
+config_use_cvmix = .false.
+config_use_bulk_wind_stress = .true.
+config_use_time_varying_atmospheric_forcing = .true.
+config_time_varying_atmospheric_forcing_start_time = '2012-10-10_00:00:00'
+config_time_varying_atmospheric_forcing_reference_time = '2012-10-10_00:00:00'
+config_time_varying_atmospheric_forcing_cycle_start = 'none'
+config_time_varying_atmospheric_forcing_cycle_duration = '2-00-00_00:00:00'
+config_use_tidal_potential_forcing = .true.
+config_tidal_potential_reference_time = '2012-10-10_00:00:00'
+config_bottom_drag_mode = 'explicit'
+config_use_topographic_wave_drag = .true.
+config_topographic_wave_drag_coeff = 5.0e-3
+config_pressure_gradient_type = 'ssh_gradient'
+config_vert_coord_movement = 'impermeable_interfaces'
+config_ALE_thickness_proportionality = 'weights_only'
+config_disable_thick_vadv = .true.
+config_disable_thick_sflux = .true.
+config_disable_vel_hmix = .true.
+config_disable_vel_vmix = .true.
+config_disable_vel_vadv = .true.
+config_disable_tr_adv = .true.
+config_disable_tr_hmix = .true.
+config_disable_tr_vmix = .true.
+config_disable_tr_sflux = .true.
+config_disable_tr_nonlocalflux = .true.
+config_disable_redi_k33 = .true.
+config_open_ocean_freezing_temperature_coeff_0 = -1.8
+config_disable_tr_all_tend = .true.
+config_AM_pointwiseStats_enable = .true.

compass/ocean/tests/hurricane/fblts/forward/streams.ocean

@@ -0,0 +1,108 @@
+<streams>
+
+<immutable_stream name="mesh"
+                  type="input"
+                  filename_template="input.nc"
+                  input_interval="initial_only"/>
+
+<immutable_stream name="input"
+                  type="input"
+                  filename_template="input.nc"
+                  input_interval="initial_only"/>
+
+<immutable_stream name="restart"
+                  type="input;output"
+                  filename_template="restarts/restart.$Y-$M-$D_$h.$m.$s.nc"
+                  filename_interval="output_interval"
+                  reference_time="0001-01-01_00:00:00"
+                  clobber_mode="truncate"
+                  input_interval="initial_only"
+                  output_interval="0005_00:00:00"/>
+
+<immutable_stream name="atmospheric_forcing"
+                  type="input"
+                  filename_template="atmospheric_forcing.nc"
+                  filename_interval="none"
+                  packages="timeVaryingAtmosphericForcingPKG"
+                  input_interval="none"/>
+
+<stream name="pointLocationsInput"
+        filename_template="points.nc"
+        type="input"
+        input_interval="initial_only"
+        runtime_format="single_file"
+        mode="forward;analysis">
+
+    <var name="pointCellGlobalID"/>
+</stream>
+
+<stream name="topographic_wave_drag"
+                  type="input"
+                  filename_template="topographic_wave_drag.nc"
+                  filename_interval="none"
+                  packages="topographicWaveDragPKG"
+                  input_interval="initial_only">
+
+    <var name="topographic_wave_drag"/>
+</stream>
+
+<stream name="ltsregion"
+                  type="input"
+                  filename_template="input.nc"
+                  filename_interval="none"
+                  input_interval="initial_only">
+    <var name="LTSRegion"/>
+</stream>
+
+<stream name="output"
+        type="output"
+        filename_template="output.nc"
+        output_interval="12:00:00"
+        precision="double"
+        clobber_mode="truncate">
+
+    <var name="xtime"/>
+    <var name="ssh"/>
+    <var name="kineticEnergyCell"/>
+    <var name="layerThickness"/>
+    <var name="normalVelocity"/>
+    <var name="daysSinceStartOfSim"/>
+    <var name="LTSRegion"/>
+</stream>
+
+<stream name="pointwiseStatsOutput"
+        filename_template="pointwiseStats.nc"
+        type="output"
+        mode="forward;analysis"
+        output_interval="00:30:00"
+        packages="pointwiseStatsAMPKG"
+        clobber_mode="truncate"
+        io_type="netcdf"
+        runtime_format="single_file">
+
+    <var name="pointCellGlobalID"/>
+    <var name="xtime"/>
+    <var name="latCell"/>
+    <var name="lonCell"/>
+    <var name="ssh"/>
+    <var name="layerThickness"/>
+</stream>
+
+<stream name="globalStatsOutput"
+        type="none"
+        filename_template="globalStats.nc"
+        filename_interval="01-00-00_00:00:00"
+        clobber_mode="truncate"
+        packages="globalStatsAMPKG"
+        output_interval="none" >
+
+    <var_array name="minGlobalStats"/>
+    <var_array name="maxGlobalStats"/>
+    <var_array name="avgGlobalStats"/>
+    <var name="xtime"/>
+    <var name="daysSinceStartOfSim"/>
+    <var name="volumeCellGlobal"/>
+    <var name="CFLNumberGlobal"/>
+</stream>
+
+</streams>

compass/ocean/tests/hurricane/forward/__init__.py

@@ -41,8 +41,10 @@ def __init__(self, test_group, mesh, storm, init, use_lts):
         """
         mesh_name = mesh.mesh_name
 
-        if use_lts:
+        if use_lts == 'LTS':
             name = f'{storm}_lts'
+        elif use_lts == 'FB_LTS':
+            name = f'{storm}_fblts'
         else:
             name = storm
         subdir = os.path.join(mesh_name, name)

compass/ocean/tests/hurricane/forward/forward.py

@@ -49,13 +49,22 @@ def __init__(self, test_case, mesh, init, use_lts,
 
         super().__init__(test_case=test_case, name=name)
 
-        if use_lts:
+        if use_lts == 'LTS':
 
             self.add_namelist_file(
                 'compass.ocean.tests.hurricane.lts.forward', 'namelist.ocean')
             self.add_streams_file(
                 'compass.ocean.tests.hurricane.lts.forward', 'streams.ocean')
 
+        elif use_lts == 'FB_LTS':
+
+            self.add_namelist_file(
+                'compass.ocean.tests.hurricane.fblts.forward',
+                'namelist.ocean')
+            self.add_streams_file(
+                'compass.ocean.tests.hurricane.fblts.forward',
+                'streams.ocean')
+
         else:
 
             self.add_namelist_file(

compass/ocean/tests/hurricane/init/__init__.py

@@ -45,8 +45,10 @@ def __init__(self, test_group, mesh, storm, use_lts):
         self.mesh = mesh
         self.use_lts = use_lts
 
-        if use_lts:
+        if use_lts == 'LTS':
             name = 'init_lts'
+        elif use_lts == 'FB_LTS':
+            name = 'init_fblts'
         else:
             name = 'init'
         mesh_name = mesh.mesh_name

compass/ocean/tests/hurricane/lts/mesh/lts_regions.py

@@ -22,6 +22,7 @@ class LTSRegionsStep(Step):
         The culled mesh step containing input files to this step
     """
     def __init__(self, test_case, cull_mesh_step,
+                 num_interface_adjacent,
                  name='lts_regions', subdir='lts_regions'):
         """
         Create a new step
@@ -46,6 +47,7 @@ def __init__(self, test_case, cull_mesh_step,
             self.add_output_file(filename=file)
 
         self.cull_mesh_step = cull_mesh_step
+        self.num_interface_adjacent = num_interface_adjacent
 
     def setup(self):
         """
@@ -90,11 +92,12 @@ def run(self):
         use_progress_bar = self.log_filename is None
         label_mesh(fine_region, mesh='culled_mesh.nc',
                    graph_info='culled_graph.info', num_interface=2,
+                   num_interface_adjacent=self.num_interface_adjacent,
                    logger=self.logger, use_progress_bar=use_progress_bar)
 
 
 def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
-               logger, use_progress_bar):
+               num_interface_adjacent, logger, use_progress_bar):
 
     # read in mesh data
     ds = xr.open_dataset(mesh)
@@ -118,27 +121,28 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
 
     # first layer of cells with label 5
     changed_cells = [[], []]
+    border_cells = []
     for iedge in range(0, n_edges):
         cell1 = cells_on_edge[iedge, 0] - 1
         cell2 = cells_on_edge[iedge, 1] - 1
 
         if (cell1 != -1 and cell2 != -1):
             if (lts_rgn[cell1] == 1 and lts_rgn[cell2] == 2):
 
-                lts_rgn[cell2] = 5
-                changed_cells[0].append(cell2)
+                lts_rgn[cell1] = 5
+                changed_cells[0].append(cell1)
 
             elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 1):
 
-                lts_rgn[cell1] = 5
-                changed_cells[0].append(cell1)
+                lts_rgn[cell2] = 5
+                changed_cells[0].append(cell2)
+
+    border_cells = changed_cells[0]
 
-    # second and third layer of cells with label 5
+    # num_interface_adjacent - 1 layers of cells with label 5
     # only looping over cells changed during loop for previous layer
-    # at the end of this loop, changed_cells[0] will have the list of cells
-    # sharing edegs with the coarse cells
     logger.info('Labeling interface-adjacent fine cells...')
-    for i in range(0, 2):  # this loop creates 2 layers
+    for i in range(0, num_interface_adjacent - 1):
         changed_cells[(i + 1) % 2] = []
 
         for icell in changed_cells[i % 2]:
@@ -149,17 +153,19 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
                     cell2 = cells_on_edge[iedge - 1, 1] - 1
 
                     if (cell1 != -1 and cell2 != -1):
-                        if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 2):
+                        if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 1):
 
                             lts_rgn[cell2] = 5
                             changed_cells[(i + 1) % 2].append(cell2)
 
-                        elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 5):
+                        elif (lts_rgn[cell1] == 1 and lts_rgn[cell2] == 5):
 
                             lts_rgn[cell1] = 5
                             changed_cells[(i + 1) % 2].append(cell1)
 
-    # n layers of interface region with label 4
+    changed_cells[0] = border_cells
+
+    # num_interface layers of interface region with label 4
     logger.info('Labeling interface cells...')
     for i in range(0, num_interface):
         changed_cells[(i + 1) % 2] = []
@@ -200,7 +206,7 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
 
     changed_cells[0] = changed_cells[num_interface % 2]
 
-    # n layers of interface region with label 3
+    # num_interface layers of interface region with label 3
     for i in range(0, num_interface):
         changed_cells[(i + 1) % 2] = []
 
@@ -213,9 +219,9 @@ def label_mesh(fine_region, mesh, graph_info, num_interface,  # noqa: C901
 
                     if (cell1 != -1 and cell2 != -1):
                         # for the first layer, need to check neighbors are
-                        # 5 and 2
-                        # for further layers, need to check neighbors are
                         # 3 and 2
+                        # for further layers, need to check neighbors are
+                        # 4 and 2
                         if (i == 0):
                             if (lts_rgn[cell1] == 3 and lts_rgn[cell2] == 2):
 

compass/ocean/tests/hurricane/mesh/__init__.py

@@ -29,8 +29,10 @@ def __init__(self, test_group, mesh_name, use_lts):
         self.mesh_name = mesh_name
         self.use_lts = use_lts
 
-        if use_lts:
+        if use_lts == 'LTS':
             name = 'mesh_lts'
+        elif use_lts == 'FB_LTS':
+            name = 'mesh_fblts'
         else:
             name = 'mesh'
         subdir = '{}/{}'.format(mesh_name, name)
@@ -60,10 +62,14 @@ def __init__(self, test_group, mesh_name, use_lts):
 
         self.add_step(cull_mesh_step)
 
-        if use_lts:
-
+        if use_lts == 'LTS':
+            self.add_step(LTSRegionsStep(
+                          test_case=self, cull_mesh_step=cull_mesh_step,
+                          num_interface_adjacent=2))
+        elif use_lts == 'FB_LTS':
             self.add_step(LTSRegionsStep(
-                          test_case=self, cull_mesh_step=cull_mesh_step))
+                          test_case=self, cull_mesh_step=cull_mesh_step,
+                          num_interface_adjacent=10))
 
     def configure(self):
         """

docs/users_guide/ocean/test_groups/hurricane.rst

@@ -19,13 +19,16 @@ which can be used to simulate coastal inundation using MPAS-Ocean's
 wetting and drying scheme.
 
 The time stepping options to run the simulations include the fourth
-order Runge-Kutta scheme (RK4), which is the default, and a local
-time-stepping (LTS) scheme, see 
+order Runge-Kutta scheme (RK4), and two local time-stepping schemes.
+The first LTS scheme is based on a strong stability preserving Runge-Kutta
+scheme of order three and is called LTS3, see 
 `Lilly et al. (2023) <https://doi.org/10.1029/2022MS003327>`_
-for details. Each test case in the ``ocean/hurricane`` test group has an 
-LTS counterpart which is identified by the ``_lts`` appendix in the test case name.
-Note that the executable to be used with LTS has to be compiled with
-the following options: ``USE_LTS=true OPENMP=false``.
+for details.
+The second LTS scheme is based on a forward-backward Runge-Kutta scheme
+of order two and is called FB-LTS.
+Each test case in the ``ocean/hurricane`` test group has a counterpart 
+for each LTS scheme which is identified by appending the test case name
+with ``_lts`` for LTS3 and ``_fblts`` for FB-LTS.
 
 Shared config options
 ---------------------
@@ -96,7 +99,7 @@ since this data is not carried over from the cell culling process.
 
 .. _hurricane_mesh_lts:
 
-If the LTS option is selected for the mesh test case, an additional step
+If either LTS option is selected for the mesh test case, an additional step
 is carried out after the mesh culling. This step appropriately flags 
 the cells of the mesh according to a user defined criterion in order to
 use time-steps of different sizes on different regions of the mesh.
@@ -136,7 +139,7 @@ pointWiseStats analysis member for the forward run.
 
 compute topographic wave drag step
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-This step is carried out only if the LTS option is selected for the init test case.
+This step is carried out only if either LTS option is selected for the init test case.
 
 The reciprocal of the e-folding time, ``r_inv``, from the HyCOM model,
 is computed in this step. See 
@@ -165,6 +168,6 @@ is used for the validation.
 
 .. _hurricane_sandy_lts:
 
-If the LTS option is selected for the sandy test case, the LTS scheme
+If either LTS option is selected for the sandy test case, the LTS scheme
 is used to advance the solution in time rather than the default RK4 scheme.