diff --git a/compass/ocean/suites/lts.txt b/compass/ocean/suites/lts.txt
new file mode 100644
index 0000000000..53d10afead
--- /dev/null
+++ b/compass/ocean/suites/lts.txt
@@ -0,0 +1,9 @@
+ocean/dam_break/40cm/default_lts
+ocean/dam_break/40cm/ramp_lts
+ocean/dam_break/120cm/default_lts
+ocean/dam_break/120cm/ramp_lts
+ocean/hurricane/DEQU120at30cr10rr2/mesh_lts
+ocean/hurricane/DEQU120at30cr10rr2/init_lts
+ocean/hurricane/DEQU120at30cr10rr2/sandy_lts
+ocean/parabolic_bowl/standard/ramp_lts
+ocean/parabolic_bowl/standard/noramp_lts
diff --git a/compass/ocean/suites/wetdry.txt b/compass/ocean/suites/wetdry.txt
index 0a1307c10c..2577929041 100644
--- a/compass/ocean/suites/wetdry.txt
+++ b/compass/ocean/suites/wetdry.txt
@@ -20,6 +20,14 @@ ocean/dam_break/40cm/default
ocean/dam_break/40cm/ramp
ocean/dam_break/120cm/default
ocean/dam_break/120cm/ramp
+ocean/dam_break/40cm/default_lts
+ocean/dam_break/40cm/ramp_lts
+ocean/dam_break/120cm/default_lts
+ocean/dam_break/120cm/ramp_lts
+ocean/parabolic_bowl/standard/ramp
+ocean/parabolic_bowl/standard/noramp
+ocean/parabolic_bowl/standard/ramp_lts
+ocean/parabolic_bowl/standard/noramp_lts
ocean/isomip_plus/planar/5km/sigma/thin_film_drying_Ocean0
ocean/isomip_plus/planar/5km/sigma/thin_film_wetting_Ocean0
ocean/isomip_plus/planar/5km/single_layer/thin_film_tidal_forcing_Ocean0
diff --git a/compass/ocean/tests/dam_break/__init__.py b/compass/ocean/tests/dam_break/__init__.py
index c7b365f7fa..c063637b34 100644
--- a/compass/ocean/tests/dam_break/__init__.py
+++ b/compass/ocean/tests/dam_break/__init__.py
@@ -16,7 +16,10 @@ def __init__(self, mpas_core):
super().__init__(mpas_core=mpas_core, name='dam_break')
for resolution in [0.04, 0.12]:
- self.add_test_case(
- Default(test_group=self, resolution=resolution))
- self.add_test_case(
- Ramp(test_group=self, resolution=resolution))
+ for use_lts in [True, False]:
+ self.add_test_case(
+ Default(test_group=self, resolution=resolution,
+ use_lts=use_lts))
+ self.add_test_case(
+ Ramp(test_group=self, resolution=resolution,
+ use_lts=use_lts))
diff --git a/compass/ocean/tests/dam_break/dam_break.cfg b/compass/ocean/tests/dam_break/dam_break.cfg
index 6849d56ba3..e69de29bb2 100644
--- a/compass/ocean/tests/dam_break/dam_break.cfg
+++ b/compass/ocean/tests/dam_break/dam_break.cfg
@@ -1,5 +0,0 @@
-# Options related to the vertical grid
-[vertical_grid]
-
-# Number of vertical levels # redundant with config_dam_break_vert_levels
-vert_levels = 3
diff --git a/compass/ocean/tests/dam_break/default/__init__.py b/compass/ocean/tests/dam_break/default/__init__.py
index 3f400a82f5..95a558ac1e 100644
--- a/compass/ocean/tests/dam_break/default/__init__.py
+++ b/compass/ocean/tests/dam_break/default/__init__.py
@@ -1,8 +1,10 @@
from math import floor
-from compass.testcase import TestCase
-from compass.ocean.tests.dam_break.initial_state import InitialState
+
from compass.ocean.tests.dam_break.forward import Forward
+from compass.ocean.tests.dam_break.initial_state import InitialState
+from compass.ocean.tests.dam_break.lts.lts_regions import LTSRegions
from compass.ocean.tests.dam_break.viz import Viz
+from compass.testcase import TestCase
from compass.validate import compare_variables
@@ -17,7 +19,7 @@ class Default(TestCase):
"""
- def __init__(self, test_group, resolution):
+ def __init__(self, test_group, resolution, use_lts):
"""
Create the test case
@@ -29,22 +31,36 @@ def __init__(self, test_group, resolution):
resolution : float
The resolution of the test case in m
+ use_lts : bool
+ Whether local time-stepping is used
+
"""
- name = 'default'
+ self.use_lts = use_lts
+
+ if use_lts:
+ name = 'default_lts'
+ else:
+ name = 'default'
self.resolution = resolution
if resolution < 1.:
res_name = f'{int(resolution*1e3)}cm'
else:
res_name = f'{int(resolution)}m'
- min_tasks = int(40/(resolution/0.04)**2)
- ntasks = 10*min_tasks
+ min_tasks = int(40 / (resolution / 0.04) ** 2)
+ ntasks = 10 * min_tasks
subdir = f'{res_name}/{name}'
super().__init__(test_group=test_group, name=name,
subdir=subdir)
- self.add_step(InitialState(test_case=self))
+ init_step = InitialState(test_case=self, use_lts=use_lts)
+ self.add_step(init_step)
+
+ if use_lts:
+ self.add_step(LTSRegions(test_case=self, init_step=init_step))
+
self.add_step(Forward(test_case=self, resolution=resolution,
+ use_lts=use_lts,
ntasks=ntasks, min_tasks=min_tasks,
openmp_threads=1))
self.add_step(Viz(test_case=self))
@@ -56,11 +72,12 @@ def configure(self):
resolution = self.resolution
config = self.config
+ use_lts = self.use_lts
dc = resolution # cell width in m
dx = 13 # width of the domain in m
dy = 28 # length of the domain in m
- nx = round(dx/dc)
- ny = int(2*floor(dy/(2*dc))) # guarantee that ny is even
+ nx = round(dx / dc)
+ ny = int(2 * floor(dy / (2 * dc))) # guarantee that ny is even
config.set('dam_break', 'nx', f'{nx}', comment='the number of '
'mesh cells in the x direction')
@@ -69,10 +86,19 @@ def configure(self):
config.set('dam_break', 'dc', f'{dc}', comment='the distance '
'between adjacent cell centers')
+ if use_lts:
+ vert_levels = 1
+ config.set('vertical_grid', 'vert_levels', f'{vert_levels}',
+ comment='number of vertical levels')
+ else:
+ vert_levels = 3
+ config.set('vertical_grid', 'vert_levels', f'{vert_levels}',
+ comment='number of vertical levels')
+
def validate(self):
"""
Validate variables against a baseline
"""
variables = ['layerThickness', 'normalVelocity', 'ssh']
compare_variables(test_case=self, variables=variables,
- filename1=f'forward/output.nc')
+ filename1=f'{"forward/output.nc"}')
diff --git a/compass/ocean/tests/dam_break/forward.py b/compass/ocean/tests/dam_break/forward.py
index e7480630e5..ae6d4ba5e9 100644
--- a/compass/ocean/tests/dam_break/forward.py
+++ b/compass/ocean/tests/dam_break/forward.py
@@ -7,7 +7,8 @@ class Forward(Step):
A step for performing forward MPAS-Ocean runs as part of dam break
test cases.
"""
- def __init__(self, test_case, resolution, name='forward', subdir=None,
+ def __init__(self, test_case, resolution, use_lts,
+ name='forward', subdir=None,
ntasks=1, min_tasks=None, openmp_threads=1):
"""
Create a new test case
@@ -20,6 +21,9 @@ def __init__(self, test_case, resolution, name='forward', subdir=None,
resolution : str
The resolution of the test case
+ use_lts : bool
+ Whether local time-stepping is used
+
name : str
the name of the test case
@@ -50,18 +54,36 @@ def __init__(self, test_case, resolution, name='forward', subdir=None,
self.add_namelist_file('compass.ocean.tests.dam_break',
'namelist.forward')
- self.add_streams_file('compass.ocean.tests.dam_break',
- 'streams.forward')
-
- input_path = '../initial_state'
- self.add_input_file(filename='mesh.nc',
- target=f'{input_path}/culled_mesh.nc')
-
- self.add_input_file(filename='init.nc',
- target=f'{input_path}/ocean.nc')
- self.add_input_file(filename='graph.info',
- target=f'{input_path}/culled_graph.info')
+ if use_lts:
+ self.add_namelist_options(
+ {'config_time_integrator': "'LTS'",
+ 'config_dt_scaling_LTS': "4",
+ 'config_number_of_time_levels': "4",
+ 'config_pressure_gradient_type': "'ssh_gradient'"})
+
+ self.add_streams_file('compass.ocean.tests.dam_break.lts',
+ 'streams.forward')
+ input_path = '../lts_regions'
+ self.add_input_file(filename='mesh.nc',
+ target=f'{input_path}/lts_mesh.nc')
+ self.add_input_file(filename='graph.info',
+ target=f'{input_path}/lts_graph.info')
+ self.add_input_file(filename='init.nc',
+ target=f'{input_path}/lts_ocean.nc')
+
+ else:
+ self.add_streams_file('compass.ocean.tests.dam_break',
+ 'streams.forward')
+ input_path = '../initial_state'
+ self.add_input_file(filename='mesh.nc',
+ target=f'{input_path}/culled_mesh.nc')
+
+ self.add_input_file(filename='init.nc',
+ target=f'{input_path}/ocean.nc')
+
+ self.add_input_file(filename='graph.info',
+ target=f'{input_path}/culled_graph.info')
self.add_model_as_input()
diff --git a/compass/ocean/tests/dam_break/initial_state.py b/compass/ocean/tests/dam_break/initial_state.py
index 94e065f5af..e6bf4cfb38 100644
--- a/compass/ocean/tests/dam_break/initial_state.py
+++ b/compass/ocean/tests/dam_break/initial_state.py
@@ -1,13 +1,10 @@
import xarray
-import numpy
-import matplotlib.pyplot as plt
-
-from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.io import write_netcdf
from mpas_tools.mesh.conversion import convert, cull
+from mpas_tools.planar_hex import make_planar_hex_mesh
-from compass.step import Step
from compass.model import run_model
+from compass.step import Step
class InitialState(Step):
@@ -15,7 +12,7 @@ class InitialState(Step):
A step for creating a mesh and initial condition for dam break test
cases
"""
- def __init__(self, test_case):
+ def __init__(self, test_case, use_lts):
"""
Create the step
@@ -27,8 +24,14 @@ def __init__(self, test_case):
super().__init__(test_case=test_case, name='initial_state', ntasks=1,
min_tasks=1, openmp_threads=1)
- self.add_namelist_file('compass.ocean.tests.dam_break',
- 'namelist.init', mode='init')
+ self.use_lts = use_lts
+
+ if use_lts:
+ self.add_namelist_file('compass.ocean.tests.dam_break.lts',
+ 'namelist.init', mode='init')
+ else:
+ self.add_namelist_file('compass.ocean.tests.dam_break',
+ 'namelist.init', mode='init')
self.add_streams_file('compass.ocean.tests.dam_break',
'streams.init', mode='init')
@@ -52,7 +55,7 @@ def run(self):
dc = section.getfloat('dc')
self.update_namelist_at_runtime(
- {'config_dam_break_dc': f'{dc}'})
+ {'config_dam_break_dc': f'{dc}'})
logger.info(' * Make planar hex mesh')
dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=True,
diff --git a/compass/ocean/tests/dam_break/lts/__init__.py b/compass/ocean/tests/dam_break/lts/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/compass/ocean/tests/dam_break/lts/lts_regions.py b/compass/ocean/tests/dam_break/lts/lts_regions.py
new file mode 100644
index 0000000000..e3a6101834
--- /dev/null
+++ b/compass/ocean/tests/dam_break/lts/lts_regions.py
@@ -0,0 +1,335 @@
+import os
+
+import netCDF4 as nc
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+from mpas_tools.viz.paraview_extractor import extract_vtk
+
+from compass.step import Step
+
+
+class LTSRegions(Step):
+ """
+ A step for adding LTS regions to a global MPAS-Ocean mesh
+
+ Attributes
+ ----------
+ initial_state_step :
+ compass.ocean.tests.dam_break.initial_state.InitialState
+ The initial step containing input files to this step
+ """
+ def __init__(self, test_case, init_step,
+ name='lts_regions', subdir='lts_regions'):
+ """
+ Create a new step
+
+ Parameters
+ ----------
+ test_case : compass.Testcase
+ The test case this step belongs to
+
+ init_step :
+ compass.ocean.tests.dam_break.initial_state.InitialState
+ The initial state step containing input files to this step
+
+ name : str, optional
+ the name of the step
+
+ subdir : str, optional
+ the subdirectory for the step. The default is ``name``
+ """
+ super().__init__(test_case, name=name, subdir=subdir)
+
+ for file in ['lts_mesh.nc', 'lts_graph.info', 'lts_ocean.nc']:
+ self.add_output_file(filename=file)
+
+ self.init_step = init_step
+
+ def setup(self):
+ """
+ Set up the test case in the work directory, including downloading any
+ dependencies.
+ """
+ super().setup()
+
+ init_path = self.init_step.path
+ tgt1 = os.path.join(init_path, 'culled_mesh.nc')
+ self.add_input_file(filename='culled_mesh.nc', work_dir_target=tgt1)
+
+ tgt2 = os.path.join(init_path, 'culled_graph.info')
+ self.add_input_file(filename='culled_graph.info', work_dir_target=tgt2)
+
+ tgt3 = os.path.join(init_path, 'ocean.nc')
+ self.add_input_file(filename='ocean.nc', work_dir_target=tgt3)
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+
+ use_progress_bar = self.log_filename is None
+ label_mesh(init='ocean.nc',
+ mesh='culled_mesh.nc',
+ graph_info='culled_graph.info', num_interface=2,
+ logger=self.logger, use_progress_bar=use_progress_bar)
+
+
+def label_mesh(init, mesh, graph_info, num_interface, # noqa: C901
+ logger, use_progress_bar):
+
+ # read in mesh data
+ ds = xr.open_dataset(mesh)
+ n_cells = ds['nCells'].size
+ n_edges = ds['nEdges'].size
+ area_cell = ds['areaCell'].values
+ cells_on_edge = ds['cellsOnEdge'].values
+ edges_on_cell = ds['edgesOnCell'].values
+ x_cell = ds['xCell']
+
+ # start by setting all cells to coarse
+ lts_rgn = [2] * n_cells
+
+ # check each cell, if in the fine region, label as fine
+ logger.info('Labeling fine cells...')
+ for icell in range(0, n_cells):
+ if x_cell[icell] > 8.0:
+ lts_rgn[icell] = 1
+
+ # first layer of cells with label 5
+ changed_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)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 1):
+
+ lts_rgn[cell1] = 5
+ changed_cells[0].append(cell1)
+
+ # second and third layer 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
+ changed_cells[(i + 1) % 2] = []
+
+ for icell in changed_cells[i % 2]:
+ edges = edges_on_cell[icell]
+ for iedge in edges:
+ if iedge != 0:
+ cell1 = cells_on_edge[iedge - 1, 0] - 1
+ cell2 = cells_on_edge[iedge - 1, 1] - 1
+
+ if (cell1 != -1 and cell2 != -1):
+ if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 5
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 5):
+
+ lts_rgn[cell1] = 5
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ # n layers of interface region with label 4
+ logger.info('Labeling interface cells...')
+ for i in range(0, num_interface):
+ changed_cells[(i + 1) % 2] = []
+
+ for icell in changed_cells[i % 2]:
+ edges = edges_on_cell[icell]
+ for iedge in edges:
+ if iedge != 0:
+ cell1 = cells_on_edge[iedge - 1, 0] - 1
+ cell2 = cells_on_edge[iedge - 1, 1] - 1
+
+ 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
+ if (i == 0):
+ if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 3
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 5):
+
+ lts_rgn[cell1] = 3
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ else:
+ if (lts_rgn[cell1] == 3 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 3
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 3):
+
+ lts_rgn[cell1] = 3
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ changed_cells[0] = changed_cells[num_interface % 2]
+
+ # n layers of interface region with label 3
+ for i in range(0, num_interface):
+ changed_cells[(i + 1) % 2] = []
+
+ for icell in changed_cells[i % 2]:
+ edges = edges_on_cell[icell]
+ for iedge in edges:
+ if iedge != 0:
+ cell1 = cells_on_edge[iedge - 1, 0] - 1
+ cell2 = cells_on_edge[iedge - 1, 1] - 1
+
+ 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
+ if (i == 0):
+ if (lts_rgn[cell1] == 3 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 4
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 3):
+
+ lts_rgn[cell1] = 4
+ changed_cells[(i + 1) % 2].append(cell1)
+ else:
+ if (lts_rgn[cell1] == 4 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 4
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 4):
+
+ lts_rgn[cell1] = 4
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ # create lts_mesh.nc
+
+ logger.info('Creating lts_mesh...')
+
+ # open mesh nc file to be copied
+
+ ds_msh = xr.open_dataset(mesh)
+ ds_ltsmsh = ds_msh.copy(deep=True)
+ ltsmsh_name = 'lts_mesh.nc'
+ write_netcdf(ds_ltsmsh, ltsmsh_name)
+ mshnc = nc.Dataset(ltsmsh_name, 'a', format='NETCDF4_64BIT_OFFSET')
+
+ try:
+ # try to get LTSRegion and assign new value
+ lts_rgn_NC = mshnc.variables['LTSRegion']
+ lts_rgn_NC[:] = lts_rgn[:]
+ except KeyError:
+ # create new variable
+ ncells_NC = mshnc.dimensions['nCells'].name
+ lts_rgn_NC = mshnc.createVariable('LTSRegion', np.int32, (ncells_NC,))
+
+ # set new variable
+ lts_rgn_NC[:] = lts_rgn[:]
+
+ mshnc.close()
+
+ # open init nc file to be copied
+
+ ds_init = xr.open_dataset(init)
+ ds_ltsinit = ds_init.copy(deep=True)
+ ltsinit_name = 'lts_ocean.nc'
+ write_netcdf(ds_ltsinit, ltsinit_name)
+ initnc = nc.Dataset(ltsinit_name, 'a', format='NETCDF4_64BIT_OFFSET')
+
+ try:
+ # try to get LTSRegion and assign new value
+ lts_rgn_NC = initnc.variables['LTSRegion']
+ lts_rgn_NC[:] = lts_rgn[:]
+ except KeyError:
+ # create new variable
+ ncells_NC = initnc.dimensions['nCells'].name
+ lts_rgn_NC = initnc.createVariable('LTSRegion', np.int32, (ncells_NC,))
+
+ # set new variable
+ lts_rgn_NC[:] = lts_rgn[:]
+
+ initnc.close()
+
+ extract_vtk(ignore_time=True, lonlat=0,
+ dimension_list=['maxEdges='],
+ variable_list=['allOnCells'],
+ filename_pattern=ltsmsh_name,
+ out_dir='lts_mesh_vtk', use_progress_bar=use_progress_bar)
+
+ # label cells in graph.info
+
+ logger.info('Weighting ' + graph_info + '...')
+
+ fine_cells = 0
+ coarse_cells = 0
+
+ newf = ""
+ with open(graph_info, 'r') as f:
+ lines = f.readlines()
+ # OLD METIS SETUP: fine and interface are separate for METIS
+ # CURRENT METIS SETUP: fine and interface are together for METIS
+
+ # OLD METIS SETUP:
+ # newf += lines[0].strip() + " 010 3 \n"
+ # CURRENT METIS SETUP:
+ newf += lines[0].strip() + " 010 2 \n"
+
+ for icell in range(1, len(lines)):
+ if (lts_rgn[icell - 1] == 1 or lts_rgn[icell - 1] == 5): # fine
+ # OLD METIS SETUP:
+ # newf+= "0 1 0 " + lines[icell].strip() + "\n"
+ # CURRENT METIS SETUP:
+ newf += "0 1 " + lines[icell].strip() + "\n"
+ fine_cells = fine_cells + 1
+
+ elif (lts_rgn[icell - 1] == 2): # coarse
+ # OLD METIS SETUP:
+ # newf+= "1 0 0 " + lines[icell].strip() + "\n"
+ # CURRENT METIS SETUP:
+ newf += "1 0 " + lines[icell].strip() + "\n"
+ coarse_cells = coarse_cells + 1
+
+ else: # interface 1 and 2
+ # OLD METIS SETUP:
+ # newf+= "0 0 1 " + lines[icell].strip() + "\n"
+ # CURRENT METIS SETUP:
+ newf += "0 1 " + lines[icell].strip() + "\n"
+ coarse_cells = coarse_cells + 1
+
+ with open('lts_graph.info', 'w') as f:
+ f.write(newf)
+
+ max_area = max(area_cell)
+ min_area = min(area_cell)
+ max_width = 2 * np.sqrt(max_area / np.pi) / 1000
+ min_width = 2 * np.sqrt(min_area / np.pi) / 1000
+ area_ratio = max_area / min_area
+ width_ratio = max_width / min_width
+ number_ratio = coarse_cells / fine_cells
+
+ txt = 'number of fine cells = {}\n'.format(fine_cells)
+ txt += 'number of coarse cells = {}\n'.format(coarse_cells)
+ txt += 'ratio of coarse to fine cells = {}\n'.format(number_ratio)
+ txt += 'ratio of largest to smallest cell area = {}\n'.format(area_ratio)
+ txt += 'ratio of largest to smallest cell width = {}\n'.format(width_ratio)
+ txt += 'number of interface layers = {}\n'.format(num_interface)
+
+ logger.info(txt)
+
+ with open('lts_mesh_info.txt', 'w') as f:
+ f.write(txt)
diff --git a/compass/ocean/tests/dam_break/lts/namelist.init b/compass/ocean/tests/dam_break/lts/namelist.init
new file mode 100644
index 0000000000..e1c69416a7
--- /dev/null
+++ b/compass/ocean/tests/dam_break/lts/namelist.init
@@ -0,0 +1,5 @@
+config_ocean_run_mode = 'init'
+config_init_configuration='dam_break'
+config_dam_break_vert_levels=1
+config_dam_break_R0=12.0
+config_drying_min_cell_height=1.0e-2
diff --git a/compass/ocean/tests/dam_break/lts/streams.forward b/compass/ocean/tests/dam_break/lts/streams.forward
new file mode 100644
index 0000000000..9ab90e7783
--- /dev/null
+++ b/compass/ocean/tests/dam_break/lts/streams.forward
@@ -0,0 +1,32 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/dam_break/namelist.forward b/compass/ocean/tests/dam_break/namelist.forward
index 48bb1d3b8b..4ce329af42 100644
--- a/compass/ocean/tests/dam_break/namelist.forward
+++ b/compass/ocean/tests/dam_break/namelist.forward
@@ -13,4 +13,4 @@ config_drying_min_cell_height=1e-6
config_thickness_flux_type='upwind'
config_vert_coord_movement='impermeable_interfaces'
config_use_debugTracers=.false.
-config_disable_tr_all_tend = .true.
+config_disable_tr_all_tend=.true.
diff --git a/compass/ocean/tests/dam_break/namelist.init b/compass/ocean/tests/dam_break/namelist.init
index e6933120ca..5818bff3b1 100644
--- a/compass/ocean/tests/dam_break/namelist.init
+++ b/compass/ocean/tests/dam_break/namelist.init
@@ -1,4 +1,4 @@
-config_ocean_run_mode = 'init'
+config_ocean_run_mode='init'
config_init_configuration='dam_break'
config_dam_break_vert_levels=3
config_dam_break_R0=12.0
diff --git a/compass/ocean/tests/dam_break/ramp/__init__.py b/compass/ocean/tests/dam_break/ramp/__init__.py
index d72b9234fc..91553be037 100644
--- a/compass/ocean/tests/dam_break/ramp/__init__.py
+++ b/compass/ocean/tests/dam_break/ramp/__init__.py
@@ -2,6 +2,7 @@
from compass.ocean.tests.dam_break.forward import Forward
from compass.ocean.tests.dam_break.initial_state import InitialState
+from compass.ocean.tests.dam_break.lts.lts_regions import LTSRegions
from compass.ocean.tests.dam_break.viz import Viz
from compass.testcase import TestCase
from compass.validate import compare_variables
@@ -18,7 +19,7 @@ class Ramp(TestCase):
"""
- def __init__(self, test_group, resolution):
+ def __init__(self, test_group, resolution, use_lts):
"""
Create the test case
@@ -30,8 +31,14 @@ def __init__(self, test_group, resolution):
resolution : float
The resolution of the test case in m
+ use_lts : bool
+ Whether local time-stepping is used
+
"""
- name = 'ramp'
+ if use_lts:
+ name = 'ramp_lts'
+ else:
+ name = 'ramp'
self.resolution = resolution
if resolution < 1.:
@@ -44,8 +51,14 @@ def __init__(self, test_group, resolution):
super().__init__(test_group=test_group, name=name,
subdir=subdir)
- self.add_step(InitialState(test_case=self))
+ init_step = InitialState(test_case=self, use_lts=use_lts)
+ self.add_step(init_step)
+
+ if use_lts:
+ self.add_step(LTSRegions(test_case=self, init_step=init_step))
+
forward_step = Forward(test_case=self, resolution=resolution,
+ use_lts=use_lts,
ntasks=ntasks, min_tasks=min_tasks,
openmp_threads=1)
forward_step.add_namelist_options({'config_zero_drying_velocity_ramp':
diff --git a/compass/ocean/tests/dam_break/viz/__init__.py b/compass/ocean/tests/dam_break/viz/__init__.py
index 765c6746a9..ab01f8db8a 100644
--- a/compass/ocean/tests/dam_break/viz/__init__.py
+++ b/compass/ocean/tests/dam_break/viz/__init__.py
@@ -1,10 +1,11 @@
-import xarray
-import numpy
+from collections import OrderedDict
+
import matplotlib.pyplot as plt
+import numpy
import pandas as pd
-from scipy import spatial
+import xarray
from PIL import Image
-from collections import OrderedDict
+from scipy import spatial
from compass.step import Step
@@ -87,22 +88,22 @@ def run(self):
ii = 0
- fig = plt.figure()
for cell in station:
ii += 1
- ax = plt.subplot(3, 2, ii+1)
+ ax = plt.subplot(3, 2, ii + 1)
# MPAS-O simulation results
- mpaso = plt.plot(numpy.arange(0, dt*nt, dt), wct[:, station[cell]],
+ mpaso = plt.plot(numpy.arange(0, dt * nt, dt),
+ wct[:, station[cell]],
color='#228B22', linewidth=2, alpha=0.6)
# Measured data
- data = pd.read_csv('Station'+cell+'.csv', header=None)
+ data = pd.read_csv('Station' + cell + '.csv', header=None)
measured = plt.scatter(data[0], data[1], 4,
marker='o', color='k')
# ROMS simulation results (Warner et al., 2013)
- roms_data = pd.read_csv(cell+'-sim.csv', header=None)
+ roms_data = pd.read_csv(cell + '-sim.csv', header=None)
roms = plt.scatter(roms_data[0], roms_data[1], 4,
marker='v', color='b')
@@ -110,7 +111,7 @@ def run(self):
plt.xticks(numpy.arange(0, 11, 2))
plt.ylim(0, 0.7)
plt.yticks(numpy.arange(0, 0.7, 0.2))
- plt.text(3.5, 0.55, 'Station '+cell)
+ plt.text(3.5, 0.55, 'Station ' + cell)
if ii % 2 == 0:
plt.ylabel('h (m)')
diff --git a/compass/ocean/tests/parabolic_bowl/__init__.py b/compass/ocean/tests/parabolic_bowl/__init__.py
index 167b3842bb..23b24fed93 100644
--- a/compass/ocean/tests/parabolic_bowl/__init__.py
+++ b/compass/ocean/tests/parabolic_bowl/__init__.py
@@ -15,6 +15,10 @@ def __init__(self, mpas_core):
super().__init__(mpas_core=mpas_core, name='parabolic_bowl')
for wetdry in ['standard']: # , 'subgrid']:
for ramp_type in ['ramp', 'noramp']:
- self.add_test_case(
- Default(test_group=self, ramp_type=ramp_type,
- wetdry=wetdry))
+ # note: LTS has only standard W/D
+ for use_lts in [True, False]:
+ self.add_test_case(
+ Default(test_group=self,
+ ramp_type=ramp_type,
+ wetdry=wetdry,
+ use_lts=use_lts))
diff --git a/compass/ocean/tests/parabolic_bowl/default/__init__.py b/compass/ocean/tests/parabolic_bowl/default/__init__.py
index ccb7fd437f..c1b48cf8fc 100644
--- a/compass/ocean/tests/parabolic_bowl/default/__init__.py
+++ b/compass/ocean/tests/parabolic_bowl/default/__init__.py
@@ -3,6 +3,7 @@
from compass.config import CompassConfigParser
from compass.ocean.tests.parabolic_bowl.forward import Forward
from compass.ocean.tests.parabolic_bowl.initial_state import InitialState
+from compass.ocean.tests.parabolic_bowl.lts.lts_regions import LTSRegions
from compass.ocean.tests.parabolic_bowl.viz import Viz
from compass.testcase import TestCase
from compass.validate import compare_variables
@@ -18,7 +19,7 @@ class Default(TestCase):
The type of vertical coordinate (``ramp``, ``noramp``, etc.)
"""
- def __init__(self, test_group, ramp_type, wetdry):
+ def __init__(self, test_group, ramp_type, wetdry, use_lts):
"""
Create the test case
@@ -32,30 +33,41 @@ def __init__(self, test_group, ramp_type, wetdry):
wetdry : str
The type of wetting and drying used (``standard``, ``subgrid``)
- """
- name = f'{wetdry}_{ramp_type}'
- subdir = f'{wetdry}/{ramp_type}'
+ use_lts : bool
+ Whether local time-stepping is used
+ """
+ if use_lts:
+ name = f'{wetdry}_{ramp_type}_lts'
+ else:
+ name = f'{wetdry}_{ramp_type}'
+
+ if use_lts:
+ subdir = f'{wetdry}/{ramp_type}_lts'
+ else:
+ subdir = f'{wetdry}/{ramp_type}'
super().__init__(test_group=test_group, name=name,
subdir=subdir)
self.resolutions = None
self.wetdry = wetdry
self.ramp_type = ramp_type
+ self.use_lts = use_lts
# add the steps with default resolutions so they can be listed
config = CompassConfigParser()
config.add_from_package('compass.ocean.tests.parabolic_bowl',
'parabolic_bowl.cfg')
- self._setup_steps(config)
+ self._setup_steps(config, use_lts)
def configure(self):
"""
Set config options for the test case
"""
config = self.config
+ use_lts = self.use_lts
# set up the steps again in case a user has provided new resolutions
- self._setup_steps(config)
+ self._setup_steps(config, use_lts)
self.update_cores()
@@ -97,7 +109,7 @@ def update_cores(self):
str(min_tasks),
comment=f'Minimum core count for {res_name} mesh')
- def _setup_steps(self, config):
+ def _setup_steps(self, config, use_lts):
""" setup steps given resolutions """
default_resolutions = '20, 10, 5'
@@ -124,16 +136,24 @@ def _setup_steps(self, config):
res_name = f'{resolution}km'
- self.add_step(InitialState(test_case=self,
- name=f'initial_state_{res_name}',
- resolution=resolution,
- wetdry=self.wetdry))
+ init_step = InitialState(test_case=self,
+ name=f'initial_state_{res_name}',
+ resolution=resolution,
+ wetdry=self.wetdry)
+ self.add_step(init_step)
+ if use_lts:
+ self.add_step(LTSRegions(test_case=self,
+ init_step=init_step,
+ name=f'lts_regions_{res_name}',
+ subdir=f'lts_regions_{res_name}'))
self.add_step(Forward(test_case=self,
name=f'forward_{res_name}',
+ use_lts=use_lts,
resolution=resolution,
ramp_type=self.ramp_type,
wetdry=self.wetdry))
- self.add_step(Viz(test_case=self, resolutions=resolutions))
+ self.add_step(Viz(test_case=self, resolutions=resolutions,
+ use_lts=use_lts))
def validate(self):
"""
diff --git a/compass/ocean/tests/parabolic_bowl/forward.py b/compass/ocean/tests/parabolic_bowl/forward.py
index 62fdd4d49d..5863bd9de2 100644
--- a/compass/ocean/tests/parabolic_bowl/forward.py
+++ b/compass/ocean/tests/parabolic_bowl/forward.py
@@ -9,7 +9,8 @@ class Forward(Step):
A step for performing forward MPAS-Ocean runs as part of parabolic bowl
test cases.
"""
- def __init__(self, test_case, resolution, name,
+ def __init__(self, test_case, resolution,
+ name, use_lts,
ramp_type='ramp', coord_type='single_layer',
wetdry='standard'):
"""
@@ -24,16 +25,19 @@ def __init__(self, test_case, resolution, name,
The resolution of the test case
name : str
- the name of the test case
+ The name of the test case
+
+ use_lts : bool
+ Whether local time-stepping is used
subdir : str, optional
- the subdirectory for the step. The default is ``name``
+ The subdirectory for the step. The default is ``name``
coord_type : str, optional
- vertical coordinate configuration
+ Vertical coordinate configuration
ramp_type : str, optional
- vertical coordinate configuration
+ Vertical coordinate configuration
"""
self.resolution = resolution
@@ -51,18 +55,32 @@ def __init__(self, test_case, resolution, name,
self.add_namelist_file('compass.ocean.tests.parabolic_bowl',
'namelist.ramp.forward')
- self.add_streams_file('compass.ocean.tests.parabolic_bowl',
- 'streams.forward')
-
- input_path = f'../initial_state_{res_name}'
- self.add_input_file(filename='mesh.nc',
- target=f'{input_path}/culled_mesh.nc')
-
- self.add_input_file(filename='init.nc',
- target=f'{input_path}/ocean.nc')
-
- self.add_input_file(filename='graph.info',
- target=f'{input_path}/culled_graph.info')
+ if use_lts:
+ self.add_namelist_options(
+ {'config_time_integrator': "'LTS'",
+ 'config_dt_scaling_LTS': "4",
+ 'config_number_of_time_levels': "4",
+ 'config_pressure_gradient_type': "'ssh_gradient'"})
+ self.add_streams_file('compass.ocean.tests.parabolic_bowl.lts',
+ 'streams.forward')
+ input_path = f'../lts_regions_{res_name}'
+ self.add_input_file(filename='mesh.nc',
+ target=f'{input_path}/lts_mesh.nc')
+ self.add_input_file(filename='graph.info',
+ target=f'{input_path}/lts_graph.info')
+ self.add_input_file(filename='init.nc',
+ target=f'{input_path}/lts_ocean.nc')
+
+ else:
+ self.add_streams_file('compass.ocean.tests.parabolic_bowl',
+ 'streams.forward')
+ input_path = f'../initial_state_{res_name}'
+ self.add_input_file(filename='mesh.nc',
+ target=f'{input_path}/culled_mesh.nc')
+ self.add_input_file(filename='init.nc',
+ target=f'{input_path}/ocean.nc')
+ self.add_input_file(filename='graph.info',
+ target=f'{input_path}/culled_graph.info')
self.add_model_as_input()
diff --git a/compass/ocean/tests/parabolic_bowl/lts/__init__.py b/compass/ocean/tests/parabolic_bowl/lts/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/compass/ocean/tests/parabolic_bowl/lts/lts_regions.py b/compass/ocean/tests/parabolic_bowl/lts/lts_regions.py
new file mode 100644
index 0000000000..4885459736
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/lts/lts_regions.py
@@ -0,0 +1,340 @@
+import math
+import os
+
+import netCDF4 as nc
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+from mpas_tools.viz.paraview_extractor import extract_vtk
+
+from compass.step import Step
+
+
+class LTSRegions(Step):
+ """
+ A step for adding LTS regions to a global MPAS-Ocean mesh
+
+ Attributes
+ ----------
+ initial_state_step :
+ compass.ocean.tests.dam_break.initial_state.InitialState
+ The initial step containing input files to this step
+ """
+ def __init__(self, test_case, init_step,
+ name='lts_regions', subdir='lts_regions'):
+ """
+ Create a new step
+
+ Parameters
+ ----------
+ test_case : compass.Testcase
+ The test case this step belongs to
+
+ init_step :
+ compass.ocean.tests.dam_break.initial_state.InitialState
+ The initial state step containing input files to this step
+
+ name : str, optional
+ the name of the step
+
+ subdir : str, optional
+ the subdirectory for the step. The default is ``name``
+ """
+ super().__init__(test_case, name=name, subdir=subdir)
+
+ for file in ['lts_mesh.nc', 'lts_graph.info', 'lts_ocean.nc']:
+ self.add_output_file(filename=file)
+
+ self.init_step = init_step
+
+ def setup(self):
+ """
+ Set up the test case in the work directory, including downloading any
+ dependencies.
+ """
+ super().setup()
+
+ init_path = self.init_step.path
+ tgt1 = os.path.join(init_path, 'culled_mesh.nc')
+ self.add_input_file(filename='culled_mesh.nc', work_dir_target=tgt1)
+
+ tgt2 = os.path.join(init_path, 'culled_graph.info')
+ self.add_input_file(filename='culled_graph.info', work_dir_target=tgt2)
+
+ tgt3 = os.path.join(init_path, 'ocean.nc')
+ self.add_input_file(filename='ocean.nc', work_dir_target=tgt3)
+
+ def run(self):
+ """
+ Run this step of the test case
+ """
+
+ use_progress_bar = self.log_filename is None
+ label_mesh(init='ocean.nc',
+ mesh='culled_mesh.nc',
+ graph_info='culled_graph.info', num_interface=2,
+ logger=self.logger, use_progress_bar=use_progress_bar)
+
+
+def label_mesh(init, mesh, graph_info, num_interface, # noqa: C901
+ logger, use_progress_bar):
+
+ # read in mesh data
+ ds = xr.open_dataset(mesh)
+ n_cells = ds['nCells'].size
+ n_edges = ds['nEdges'].size
+ area_cell = ds['areaCell'].values
+ cells_on_edge = ds['cellsOnEdge'].values
+ edges_on_cell = ds['edgesOnCell'].values
+ x_cell = ds['xCell']
+ y_cell = ds['yCell']
+
+ # start by setting all cells to coarse
+ lts_rgn = [2] * n_cells
+
+ # check each cell, if in the fine region, label as fine
+ logger.info('Labeling fine cells...')
+ for icell in range(0, n_cells):
+ xC = 745000.0
+ yC = 701480.577065
+ radius = np.sqrt((x_cell[icell] - xC) ** 2 + (y_cell[icell] - yC) ** 2)
+ if radius <= 192500.0:
+ lts_rgn[icell] = 1
+
+ # first layer of cells with label 5
+ changed_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)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 1):
+
+ lts_rgn[cell1] = 5
+ changed_cells[0].append(cell1)
+
+ # second and third layer 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
+ changed_cells[(i + 1) % 2] = []
+
+ for icell in changed_cells[i % 2]:
+ edges = edges_on_cell[icell]
+ for iedge in edges:
+ if iedge != 0:
+ cell1 = cells_on_edge[iedge - 1, 0] - 1
+ cell2 = cells_on_edge[iedge - 1, 1] - 1
+
+ if (cell1 != -1 and cell2 != -1):
+ if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 5
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 5):
+
+ lts_rgn[cell1] = 5
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ # n layers of interface region with label 4
+ logger.info('Labeling interface cells...')
+ for i in range(0, num_interface):
+ changed_cells[(i + 1) % 2] = []
+
+ for icell in changed_cells[i % 2]:
+ edges = edges_on_cell[icell]
+ for iedge in edges:
+ if iedge != 0:
+ cell1 = cells_on_edge[iedge - 1, 0] - 1
+ cell2 = cells_on_edge[iedge - 1, 1] - 1
+
+ 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
+ if (i == 0):
+ if (lts_rgn[cell1] == 5 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 3
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 5):
+
+ lts_rgn[cell1] = 3
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ else:
+ if (lts_rgn[cell1] == 3 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 3
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 3):
+
+ lts_rgn[cell1] = 3
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ changed_cells[0] = changed_cells[num_interface % 2]
+
+ # n layers of interface region with label 3
+ for i in range(0, num_interface):
+ changed_cells[(i + 1) % 2] = []
+
+ for icell in changed_cells[i % 2]:
+ edges = edges_on_cell[icell]
+ for iedge in edges:
+ if iedge != 0:
+ cell1 = cells_on_edge[iedge - 1, 0] - 1
+ cell2 = cells_on_edge[iedge - 1, 1] - 1
+
+ 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
+ if (i == 0):
+ if (lts_rgn[cell1] == 3 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 4
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 3):
+
+ lts_rgn[cell1] = 4
+ changed_cells[(i + 1) % 2].append(cell1)
+ else:
+ if (lts_rgn[cell1] == 4 and lts_rgn[cell2] == 2):
+
+ lts_rgn[cell2] = 4
+ changed_cells[(i + 1) % 2].append(cell2)
+
+ elif (lts_rgn[cell1] == 2 and lts_rgn[cell2] == 4):
+
+ lts_rgn[cell1] = 4
+ changed_cells[(i + 1) % 2].append(cell1)
+
+ # create lts_mesh.nc
+
+ logger.info('Creating lts_mesh...')
+
+ # open mesh nc file to be copied
+
+ ds_msh = xr.open_dataset(mesh)
+ ds_ltsmsh = ds_msh.copy(deep=True)
+ ltsmsh_name = 'lts_mesh.nc'
+ write_netcdf(ds_ltsmsh, ltsmsh_name)
+ mshnc = nc.Dataset(ltsmsh_name, 'a', format='NETCDF4_64BIT_OFFSET')
+
+ try:
+ # try to get LTSRegion and assign new value
+ lts_rgn_NC = mshnc.variables['LTSRegion']
+ lts_rgn_NC[:] = lts_rgn[:]
+ except KeyError:
+ # create new variable
+ ncells_NC = mshnc.dimensions['nCells'].name
+ lts_rgn_NC = mshnc.createVariable('LTSRegion', np.int32, (ncells_NC,))
+
+ # set new variable
+ lts_rgn_NC[:] = lts_rgn[:]
+
+ mshnc.close()
+
+ # open init nc file to be copied
+
+ ds_init = xr.open_dataset(init)
+ ds_ltsinit = ds_init.copy(deep=True)
+ ltsinit_name = 'lts_ocean.nc'
+ write_netcdf(ds_ltsinit, ltsinit_name)
+ initnc = nc.Dataset(ltsinit_name, 'a', format='NETCDF4_64BIT_OFFSET')
+
+ try:
+ # try to get LTSRegion and assign new value
+ lts_rgn_NC = initnc.variables['LTSRegion']
+ lts_rgn_NC[:] = lts_rgn[:]
+ except KeyError:
+ # create new variable
+ ncells_NC = initnc.dimensions['nCells'].name
+ lts_rgn_NC = initnc.createVariable('LTSRegion', np.int32, (ncells_NC,))
+
+ # set new variable
+ lts_rgn_NC[:] = lts_rgn[:]
+
+ initnc.close()
+
+ extract_vtk(ignore_time=True, lonlat=0,
+ dimension_list=['maxEdges='],
+ variable_list=['allOnCells'],
+ filename_pattern=ltsmsh_name,
+ out_dir='lts_mesh_vtk', use_progress_bar=use_progress_bar)
+
+ # label cells in graph.info
+
+ logger.info('Weighting ' + graph_info + '...')
+
+ fine_cells = 0
+ coarse_cells = 0
+
+ newf = ""
+ with open(graph_info, 'r') as f:
+ lines = f.readlines()
+ # OLD METIS SETUP: fine and interface are separate for METIS
+ # CURRENT METIS SETUP: fine and interface are together for METIS
+
+ # OLD METIS SETUP:
+ # newf += lines[0].strip() + " 010 3 \n"
+ # CURRENT METIS SETUP:
+ newf += lines[0].strip() + " 010 2 \n"
+
+ for icell in range(1, len(lines)):
+ if (lts_rgn[icell - 1] == 1 or lts_rgn[icell - 1] == 5): # fine
+ # OLD METIS SETUP:
+ # newf+= "0 1 0 " + lines[icell].strip() + "\n"
+ # CURRENT METIS SETUP:
+ newf += "0 1 " + lines[icell].strip() + "\n"
+ fine_cells = fine_cells + 1
+
+ elif (lts_rgn[icell - 1] == 2): # coarse
+ # OLD METIS SETUP:
+ # newf+= "1 0 0 " + lines[icell].strip() + "\n"
+ # CURRENT METIS SETUP:
+ newf += "1 0 " + lines[icell].strip() + "\n"
+ coarse_cells = coarse_cells + 1
+
+ else: # interface 1 and 2
+ # OLD METIS SETUP:
+ # newf+= "0 0 1 " + lines[icell].strip() + "\n"
+ # CURRENT METIS SETUP:
+ newf += "0 1 " + lines[icell].strip() + "\n"
+ coarse_cells = coarse_cells + 1
+
+ with open('lts_graph.info', 'w') as f:
+ f.write(newf)
+
+ max_area = max(area_cell)
+ min_area = min(area_cell)
+ max_width = 2 * np.sqrt(max_area / math.pi) / 1000
+ min_width = 2 * np.sqrt(min_area / math.pi) / 1000
+ area_ratio = max_area / min_area
+ width_ratio = max_width / min_width
+ number_ratio = coarse_cells / fine_cells
+
+ txt = 'number of fine cells = {}\n'.format(fine_cells)
+ txt += 'number of coarse cells = {}\n'.format(coarse_cells)
+ txt += 'ratio of coarse to fine cells = {}\n'.format(number_ratio)
+ txt += 'ratio of largest to smallest cell area = {}\n'.format(area_ratio)
+ txt += 'ratio of largest to smallest cell width = {}\n'.format(width_ratio)
+ txt += 'number of interface layers = {}\n'.format(num_interface)
+
+ logger.info(txt)
+
+ with open('lts_mesh_info.txt', 'w') as f:
+ f.write(txt)
diff --git a/compass/ocean/tests/parabolic_bowl/lts/streams.forward b/compass/ocean/tests/parabolic_bowl/lts/streams.forward
new file mode 100644
index 0000000000..92f71fac06
--- /dev/null
+++ b/compass/ocean/tests/parabolic_bowl/lts/streams.forward
@@ -0,0 +1,43 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/compass/ocean/tests/parabolic_bowl/viz/__init__.py b/compass/ocean/tests/parabolic_bowl/viz/__init__.py
index 92d8bcad30..26bc2b9de4 100644
--- a/compass/ocean/tests/parabolic_bowl/viz/__init__.py
+++ b/compass/ocean/tests/parabolic_bowl/viz/__init__.py
@@ -18,7 +18,7 @@ class Viz(Step):
Attributes
----------
"""
- def __init__(self, test_case, resolutions):
+ def __init__(self, test_case, resolutions, use_lts):
"""
Create the step
@@ -30,6 +30,7 @@ def __init__(self, test_case, resolutions):
super().__init__(test_case=test_case, name='viz')
self.resolutions = resolutions
+ self.use_lts = use_lts
for res in resolutions:
self.add_input_file(filename=f'output_{res}km.nc',
@@ -193,9 +194,15 @@ def rmse_plots(self):
Plot convergence curves
"""
+ ramp_name = 'ramp'
+ noramp_name = 'noramp'
+ if self.use_lts:
+ ramp_name = 'ramp_lts'
+ noramp_name = 'noramp_lts'
+
comparisons = []
- cases = {'standard_ramp': '../../../standard/ramp/viz',
- 'standard_noramp': '../../../standard/noramp/viz'}
+ cases = {'standard_ramp': f'../../../standard/{ramp_name}/viz',
+ 'standard_noramp': f'../../../standard/{noramp_name}/viz'}
for case in cases:
include = True
for res in self.resolutions:
diff --git a/docs/users_guide/ocean/test_groups/dam_break.rst b/docs/users_guide/ocean/test_groups/dam_break.rst
index 5c9ace7152..485b6cb818 100644
--- a/docs/users_guide/ocean/test_groups/dam_break.rst
+++ b/docs/users_guide/ocean/test_groups/dam_break.rst
@@ -51,3 +51,11 @@ ramp
the factor that scales velocities and velocity tendencies is ramped over a
given layer thickness range rather than a binary switch at the minimum
thickness. ``RES`` is either 40cm or 120cm.
+
+lts
+---
+
+Both the ``default`` and ``ramp`` test cases can be run with the ``lts`` variant
+which uses local time-stepping (LTS) as time integrator. Note that the tests
+verify the ability of the LTS scheme to run correctly with wetting and drying
+and are not designed to leverage the LTS capability of producing faster runs.
diff --git a/docs/users_guide/ocean/test_groups/parabolic_bowl.rst b/docs/users_guide/ocean/test_groups/parabolic_bowl.rst
index 8030a395ef..2093f4d677 100644
--- a/docs/users_guide/ocean/test_groups/parabolic_bowl.rst
+++ b/docs/users_guide/ocean/test_groups/parabolic_bowl.rst
@@ -169,3 +169,11 @@ viz
The visualization step can be configured to plot the timeseries for an
arbitrary set of coordinates by setting ``points``. Also, the interval
between contour plot time snaps can be controlled with ``plot_interval``.
+
+lts
+~~~
+
+Both the ``ramp`` and ``noramp`` test cases can be run with the ``lts`` variant
+which uses local time-stepping (LTS) as time integrator. Note that the tests
+verify the ability of the LTS scheme to run correctly with wetting and drying
+and are not designed to leverage the LTS capability of producing faster runs.