Add support for creating mapping files with mbtempest

deploy/bootstrap.py

@@ -291,7 +291,7 @@ def build_conda_env(config, env_type, recreate, mpi, conda_mpi, version,
                             mpi_prefix=mpi_prefix,
                             include_mache=not local_mache)
 
-        for package in ['esmf', 'geometric_features', 'mache', 'metis',
+        for package in ['esmf', 'geometric_features', 'mache', 'metis', 'moab',
                         'mpas_tools', 'netcdf_c', 'netcdf_fortran', 'otps',
                         'parallelio', 'pnetcdf']:
             replacements[package] = config.get('deploy', package)

deploy/conda-dev-spec.template

@@ -23,6 +23,7 @@ mache={{ mache }}
 {% endif %}
 matplotlib-base>=3.9.0
 metis={{ metis }}
+moab={{ moab }}=*_tempest_*
 mpas_tools={{ mpas_tools }}
 nco
 netcdf4=*=nompi_*

docs/developers_guide/framework/remapping.md

@@ -6,30 +6,37 @@ It is frequently useful when working with observational datasets or
 visualizing MPAS data to remap between different global or regional grids and
 meshes.  The [pyremap](https://mpas-dev.github.io/pyremap/stable/) provides
 capabilities for making mapping files (which contain the weights needed to
-interpolate between meshes) and using them to remap files or 
+interpolate between meshes) and using them to remap files or
 {py:class}`xarray.Dataset` objects.  Polaris provides a step for producing
 such a mapping file.  Under the hood, `pyremap` uses the
 [ESMF_RegridWeightGen](https://earthsystemmodeling.org/docs/release/latest/ESMF_refdoc/node3.html#SECTION03020000000000000000)
-tool, which uses MPI parallelism.  To better support task parallelism, it is
+or [mbtempest](https://sigma.mcs.anl.gov/moab/offline-remapping-workflow-with-mbtempest/)
+tools, which use MPI parallelism.  To better support task parallelism, it is
 best to have each MPI task be a separate polaris step.  For this reason, we
 provide {py:class}`polaris.remap.MappingFileStep` to perform remapping.
 
 A remapping step can be added to a task either by creating a
 {py:class}`polaris.remap.MappingFileStep` object directly or by creating a
-step that descends from the class.  Here is an example of using 
+step that descends from the class.  Here is an example of using
 `MappingFileStep` directly to remap data from a WOA 2023 lon-lat grid to an
 MPAS mesh. This could happen in the task's `__init__()` or `configure()`
 method:
 
 ```python
 
 from polaris import Task
+from polaris.config import PolarisConfigParser
 from polaris.remap import MappingFileStep
 
 class MyTestCase(Task):
     def __int__(self, component):
-        step = MappingFileStep(component=component, name='make_map', ntasks=64, 
+        step = MappingFileStep(component=component, name='make_map', ntasks=64,
                                min_tasks=1, method='bilinear')
+        # add required config options related to mapping
+        config = PolarisConfigParser(filepath=filepath)
+        config.add_from_package('polaris.remap', 'mapping.cfg')
+        step.set_shared_config(config, link='my_test_case.cfg')
+
         # indicate the the mesh from another step is an input to this step
         # note: the target is relative to the step, not the task.
         step.add_input_file(filename='woa.nc',
@@ -38,7 +45,7 @@ class MyTestCase(Task):
 
         step.add_input_file(filename='mesh.nc',
                             target='../mesh/mesh.nc')
-        
+
         # you need to specify what type of source and destination mesh you
         # will use
         step.src_from_lon_lat(filename='woa.nc', lon_var='lon', lat_var='lat')
@@ -62,7 +69,7 @@ from polaris.remap import MappingFileStep
 
 
 class VizMap(MappingFileStep):
-    def __init__(self, component, name, subdir, mesh_name):
+    def __init__(self, component, name, subdir, mesh_name, config):
         super().__init__(component=component, name=name, subdir=subdir,
                          ntasks=128, min_tasks=1)
         self.mesh_name = mesh_name
@@ -81,11 +88,26 @@ class VizMap(MappingFileStep):
         super().runtime_setup()
 ```
 
-With either approach, you will need to call one of the `src_*()` methods to
-set up the source mesh or grid and one of the `dst_*()` to configure the 
-destination.  Expect for lon-lat grids, you will need to provide a name for
-the mesh or grid, typically describing its resolution and perhaps its extent
-and the region covered.
+It is important that the task that the step belongs to includes the required
+config options related to mapping.  This could be accomplished either by
+calling:
+```python
+config.add_from_package('polaris.remap', 'mapping.cfg')
+```
+or by including the corresponding config options in the task's config file:
+```cfg
+# config options related to creating mapping files
+[mapping]
+
+# The tool to use for creating mapping files: esmf or moab
+map_tool = moab
+```
+
+Whether you create a `MappingFileStep` object directly or create a subclass,
+you will need to call one of the `src_*()` methods to set up the source mesh or
+grid and one of the `dst_*()` to configure the destination.  Expect for lon-lat
+grids, you will need to provide a name for the mesh or grid, typically
+describing its resolution and perhaps its extent and the region covered.
 
 In nearly all situations, creating the mapping file is only one step in the
 workflow. After that, the mapping file will be used to remap data between

docs/developers_guide/quick_start.md

@@ -182,7 +182,8 @@ this script will also:
   build several libraries with system compilers and MPI library, including:
   [SCORPIO](https://github.com/E3SM-Project/scorpio) (parallel i/o for E3SM
   components) [ESMF](https://earthsystemmodeling.org/) (making mapping files
-  in parallel), [Trilinos](https://trilinos.github.io/),
+  in parallel), [MOAB](https://sigma.mcs.anl.gov/moab-library/),
+  [Trilinos](https://trilinos.github.io/),
   [Albany](https://github.com/sandialabs/Albany),
   [Netlib-LAPACK](http://www.netlib.org/lapack/) and
   [PETSc](https://petsc.org/). **Please uses these flags with caution, as

docs/developers_guide/troubleshooting.md

@@ -117,5 +117,5 @@ of your user config file:
 parallel_executable = mpirun -host localhost
 ```
 
-The [example config file](https://github.com/E3SM-Project/polaris/blob/main/exmaple.cfg)
+The [example config file](https://github.com/E3SM-Project/polaris/blob/main/example.cfg)
 has been updated to include this flag.

docs/users_guide/quick_start.md

@@ -17,7 +17,8 @@ documentation as soon as there is one.  Until then please refer to the
 For each polaris release, we maintain a
 [conda environment](https://docs.conda.io/en/latest/). that includes the
 `polaris` package as well as all of its dependencies and some libraries
-(currently [ESMF](https://earthsystemmodeling.org/) and
+(currently [ESMF](https://earthsystemmodeling.org/),
+[MOAB](https://sigma.mcs.anl.gov/moab-library/) and
 [SCORPIO](https://e3sm.org/scorpio-parallel-io-library/)) built with system
 MPI using [spack](https://spack.io/) on our standard machines (Anvil, Chicoma,
 Chrysalis, Compy, and  Perlmutter).  Once there is a polaris release,

polaris/model_step.py

@@ -177,8 +177,8 @@ def setup(self):
         config = self.config
         component_path = config.get('executables', 'component')
         model_basename = os.path.basename(component_path)
-        self.args = [f'./{model_basename}', '-n', self.namelist,
-                     '-s', self.streams]
+        self.args = [[f'./{model_basename}', '-n', self.namelist,
+                      '-s', self.streams]]
 
     def set_model_resources(self, ntasks=None, min_tasks=None,
                             openmp_threads=None, max_memory=None):

polaris/ocean/tasks/isomip_plus/__init__.py

@@ -34,6 +34,8 @@ def add_isomip_plus_tasks(component, mesh_type):
             config.set('spherical_mesh', 'mpas_mesh_filename',
                        'base_mesh_without_xy.nc')
 
+        config.add_from_package('polaris.remap', 'mapping.cfg')
+
         config.add_from_package('polaris.ocean.tasks.isomip_plus',
                                 'isomip_plus.cfg')
 

polaris/remap/mapping.cfg

@@ -0,0 +1,5 @@
+# config options related to creating mapping files
+[mapping]
+
+# The tool to use for creating mapping files: esmf or moab
+map_tool = moab

polaris/remap/mapping_file_step.py

@@ -366,9 +366,21 @@ def get_remapper(self):
         out_descriptor = _get_descriptor(dst)
 
         if self.map_filename is None:
+            map_tool = self.config.get('mapping', 'map_tool')
+            prefixes = {
+                'esmf': 'esmf',
+                'moab': 'mbtr'
+            }
+            suffixes = {
+                'conserve': 'aave',
+                'bilinear': 'bilin',
+                'neareststod': 'neareststod'
+            }
+            suffix = f'{prefixes[map_tool]}{suffixes[self.method]}'
+
             self.map_filename = \
                 f'map_{in_descriptor.meshName}_to_{out_descriptor.meshName}' \
-                f'_{self.method}.nc'
+                f'_{suffix}.nc'
 
         self.map_filename = os.path.abspath(os.path.join(
             self.work_dir, self.map_filename))
@@ -381,35 +393,62 @@ def runtime_setup(self):
         Create a remapper and set the command-line arguments
         """
         remapper = self.get_remapper()
-        self.args = _build_mapping_file_args(remapper, self.method,
-                                             self.expand_distance,
-                                             self.expand_factor,
+        map_tool = self.config.get('mapping', 'map_tool')
+        _check_remapper(remapper, self.method, map_tool=map_tool)
+
+        src_descriptor = remapper.sourceDescriptor
+        src_descriptor.to_scrip(self.src_mesh_filename)
+
+        dst_descriptor = remapper.destinationDescriptor
+        dst_descriptor.to_scrip(self.dst_mesh_filename,
+                                expandDist=self.expand_distance,
+                                expandFactor=self.expand_factor)
+
+        if map_tool == 'esmf':
+            self.args = _esmf_build_map_args(remapper, self.method,
+                                             src_descriptor,
                                              self.src_mesh_filename,
+                                             dst_descriptor,
+                                             self.dst_mesh_filename)
+        elif map_tool == 'moab':
+            self.args = _moab_build_map_args(remapper, self.method,
+                                             src_descriptor,
+                                             self.src_mesh_filename,
+                                             dst_descriptor,
                                              self.dst_mesh_filename)
 
 
-def _build_mapping_file_args(remapper, method, expand_distance, expand_factor,
-                             src_mesh_filename, dst_mesh_filename):
+def _check_remapper(remapper, method, map_tool):
     """
-    Get command-line arguments for making a mapping file
+    Check for inconsistencies in the remapper
     """
+    if map_tool not in ['moab', 'esmf']:
+        raise ValueError(f'Unexpected map_tool {map_tool}. Valid '
+                         f'values are "esmf" or "moab".')
+
+    if isinstance(remapper.destinationDescriptor,
+                  PointCollectionDescriptor) and \
+            method not in ['bilinear', 'neareststod']:
+        raise ValueError(f'method {method} not supported for destination '
+                         f'grid of type PointCollectionDescriptor.')
 
-    _check_remapper(remapper, method)
+    if map_tool == 'moab' and method == 'neareststod':
+        raise ValueError('method neareststod not supported by mbtempest.')
 
-    src_descriptor = remapper.sourceDescriptor
-    src_descriptor.to_scrip(src_mesh_filename)
 
-    dst_descriptor = remapper.destinationDescriptor
-    dst_descriptor.to_scrip(dst_mesh_filename, expandDist=expand_distance,
-                            expandFactor=expand_factor)
+def _esmf_build_map_args(remapper, method, src_descriptor, src_mesh_filename,
+                         dst_descriptor, dst_mesh_filename):
+    """
+    Get command-line arguments for making a mapping file with
+    ESMF_RegridWeightGen
+    """
 
     args = ['ESMF_RegridWeightGen',
             '--source', src_mesh_filename,
             '--destination', dst_mesh_filename,
             '--weight', remapper.mappingFileName,
             '--method', method,
-            '--netcdf4',
-            '--no_log']
+            '--netcdf4']
 
     if src_descriptor.regional:
         args.append('--src_regional')
@@ -420,22 +459,59 @@ def _build_mapping_file_args(remapper, method, expand_distance, expand_factor,
     if src_descriptor.regional or dst_descriptor.regional:
         args.append('--ignore_unmapped')
 
-    return args
+    return [args]
 
 
-def _check_remapper(remapper, method):
+def _moab_build_map_args(remapper, method, src_descriptor, src_mesh_filename,
+                         dst_descriptor, dst_mesh_filename):
     """
-    Check for inconsistencies in the remapper
+    Get command-line arguments for making a mapping file with mbtempest
     """
-    if isinstance(remapper.destinationDescriptor,
-                  PointCollectionDescriptor) and \
-            method not in ['bilinear', 'neareststod']:
-        raise ValueError(f'method {method} not supported for destination '
-                         'grid of type PointCollectionDescriptor.')
+    fvmethod = {
+        'conserve': 'none',
+        'bilinear': 'bilin'}
+
+    map_filename = remapper.mappingFileName
+    intx_filename = \
+        f'moab_intx_{src_descriptor.meshName}_to_{dst_descriptor.meshName}.h5m'
+
+    intx_args = [
+        'mbtempest',
+        '--type', '5',
+        '--load', src_mesh_filename,
+        '--load', dst_mesh_filename,
+        '--intx', intx_filename
+    ]
+
+    if src_descriptor.regional or dst_descriptor.regional:
+        intx_args.append('--rrmgrids')
+
+    map_args = [
+        'mbtempest',
+        '--type', '5',
+        '--load', src_mesh_filename,
+        '--load', dst_mesh_filename,
+        '--intx', intx_filename,
+        '--weights',
+        '--method', 'fv',
+        '--method', 'fv',
+        '--file', map_filename,
+        '--order', '1',
+        '--order', '1',
+        '--fvmethod', fvmethod[method]
+    ]
+
+    if method == 'conserve' and (src_descriptor.regional or
+                                 dst_descriptor.regional):
+        map_args.append('--rrmgrids')
+
+    return [intx_args, map_args]
 
 
 def _get_descriptor(info):
-    """ Get a mesh descriptor from the mesh info """
+    """
+    Get a mesh descriptor from the mesh info
+    """
     grid_type = info['type']
     if grid_type == 'mpas':
         descriptor = _get_mpas_descriptor(info)
@@ -447,11 +523,16 @@ def _get_descriptor(info):
         descriptor = _get_points_descriptor(info)
     else:
         raise ValueError(f'Unexpected grid type {grid_type}')
+
+    # for compatibility with mbtempest
+    descriptor.format = 'NETCDF3_64BIT_DATA'
     return descriptor
 
 
 def _get_mpas_descriptor(info):
-    """ Get an MPAS mesh descriptor from the given info """
+    """
+    Get an MPAS mesh descriptor from the given info
+    """
     mesh_type = info['mpas_mesh_type']
     filename = info['filename']
     mesh_name = info['name']
@@ -472,7 +553,9 @@ def _get_mpas_descriptor(info):
 
 
 def _get_lon_lat_descriptor(info):
-    """ Get a lon-lat descriptor from the given info """
+    """
+    Get a lon-lat descriptor from the given info
+    """
 
     if 'dlat' in info and 'dlon' in info:
         lon_min = info['lon_min']
@@ -510,7 +593,9 @@ def _get_lon_lat_descriptor(info):
 
 
 def _get_proj_descriptor(info):
-    """ Get a ProjectionGridDescriptor from the given info """
+    """
+    Get a ProjectionGridDescriptor from the given info
+    """
     filename = info['filename']
     grid_name = info['name']
     x = info['x']
@@ -533,7 +618,9 @@ def _get_proj_descriptor(info):
 
 
 def _get_points_descriptor(info):
-    """ Get a PointCollectionDescriptor from the given info """
+    """
+    Get a PointCollectionDescriptor from the given info
+    """
     filename = info['filename']
     collection_name = info['name']
     lon_var = info['lon']

polaris/run/serial.py

@@ -490,10 +490,11 @@ def _run_step(task, step, new_log_file, available_resources,
         if step.args is not None:
             step_logger.info('\nBypassing step\'s run() method and running '
                              'with command line args\n')
-            log_function_call(function=run_command, logger=step_logger)
-            step_logger.info('')
-            run_command(step.args, step.cpus_per_task, step.ntasks,
-                        step.openmp_threads, step.config, step.logger)
+            for args in step.args:
+                log_function_call(function=run_command, logger=step_logger)
+                step_logger.info('')
+                run_command(args, step.cpus_per_task, step.ntasks,
+                            step.openmp_threads, step.config, step.logger)
         else:
             step_logger.info('')
             log_method_call(method=step.run, logger=step_logger)