From b546dc06b3fdfe1062e437a733acbc3608d13463 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 6 Jan 2025 07:39:38 -0600 Subject: [PATCH] Update code to follow remap_topography approach MOAB commands are updated to be more like recent development in remap_topography and are organized into more helper methods. --- .../mesh/remap_mali_topography/__init__.py | 168 +++++++++++------- 1 file changed, 99 insertions(+), 69 deletions(-) diff --git a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py index 80d1b3c2a..cf189be3f 100644 --- a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py @@ -1,4 +1,5 @@ import os +import pathlib import numpy as np import xarray as xr @@ -63,7 +64,7 @@ def setup(self): mali_filename = self.config.get('remap_mali_topography', 'mali_filename') self.add_input_file( - filename='mali_topography.nc', + filename='mali_topography_orig.nc', target=mali_filename, database='mali_topo') @@ -76,33 +77,38 @@ def run(self): self._combine_topo() def _remap_mali_topo(self): - config = self.config - logger = self.logger - in_mesh_name = self.mali_ais_topo - in_descriptor = MpasCellMeshDescriptor(fileName='mali_topography.nc', - meshName=in_mesh_name) + in_descriptor = MpasCellMeshDescriptor( + fileName='mali_topography_orig.nc', + meshName=in_mesh_name) in_descriptor.format = 'NETCDF3_64BIT' - in_descriptor.to_scrip('mali_scrip.nc') + in_descriptor.to_scrip('mali.scrip.64bit.nc') + self._partition_scrip_file('mali.scrip.64bit.nc') out_mesh_name = self.mesh_name - out_descriptor = MpasCellMeshDescriptor(fileName='base_mesh.nc', - meshName=out_mesh_name) + out_descriptor = MpasCellMeshDescriptor( + fileName='base_mesh.nc', + meshName=out_mesh_name) out_descriptor.format = 'NETCDF3_64BIT' - out_descriptor.to_scrip('mpaso_scrip.nc') + out_descriptor.to_scrip('mpaso.scrip.64bit.nc') + self._partition_scrip_file('mpaso.scrip.64bit.nc') + + map_filename = \ + f'map_{in_mesh_name}_to_{out_mesh_name}_mbtraave.nc' - args = ['mbtempest', - '--type', '5', - '--load', 'mali_scrip.nc', - '--load', 'mpaso_scrip.nc', - '--intx', 'moab_intx_mali_mpaso.h5m', - '--rrmgrids'] + self._create_weights_tempest(map_filename) - run_command(args=args, cpus_per_task=self.cpus_per_task, - ntasks=self.ntasks, openmp_threads=self.openmp_threads, - config=config, logger=logger) + self._remap_to_target(map_filename) - ds_mali = xr.open_dataset('mali_topography.nc') + def _modify_mali_topo(self): + """ + Modify MALI topography to have desired fields and names + """ + logger = self.logger + config = self.config + logger.info('Modifying MALI topography fields and names') + + ds_mali = xr.open_dataset('mali_topography_orig.nc') if 'Time' in ds_mali.dims: ds_mali = ds_mali.isel(Time=0) bed = ds_mali.bedTopography @@ -142,57 +148,79 @@ def _remap_mali_topo(self): ds_in['oceanFrac'] = ocean_frac ds_in['maliFrac'] = mali_frac - mbtempest_args = {'conserve': ['--order', '1', - '--order', '1', - '--fvmethod', 'none', - '--rrmgrids'], - 'bilinear': ['--order', '1', - '--order', '1', - '--fvmethod', 'bilin']} - - suffix = {'conserve': 'mbtraave', - 'bilinear': 'mbtrbilin'} - - method = 'conserve' - mapping_file_name = \ - f'map_{in_mesh_name}_to_{out_mesh_name}_{suffix[method]}.nc' - - # split the parallel executable into constituents in case it - # includes flags - args = ['mbtempest', - '--type', '5', - '--load', 'mali_scrip.nc', - '--load', 'mpaso_scrip.nc', - '--intx', 'moab_intx_mali_mpaso.h5m', - '--weights', - '--method', 'fv', - '--method', 'fv', - '--file', mapping_file_name] + mbtempest_args[method] - - if method == 'bilinear': - # unhappy in parallel for now - check_call(args, logger) - else: - - run_command(args=args, cpus_per_task=self.cpus_per_task, - ntasks=self.ntasks, - openmp_threads=self.openmp_threads, - config=config, logger=logger) - - in_filename = f'mali_topography_{method}.nc' - out_filename = f'mali_topography_ncremap_{method}.nc' - write_netcdf(ds_in, in_filename) - - # remapping with the -P mpas flag leads to undesired - # renormalization - args = ['ncremap', - '-m', mapping_file_name, - '--vrb=1', - in_filename, - out_filename] + write_netcdf(ds_in, 'mali_topography_mod.nc') + + logger.info(' Done.') + + def _partition_scrip_file(self, in_filename): + """ + Partition SCRIP file for parallel mbtempest use + """ + logger = self.logger + logger.info('Partition SCRIP file') + + stem = pathlib.Path(in_filename).stem + h5m_filename = f'{stem}.h5m' + part_filename = f'{stem}.p{self.ntasks}.h5m' + + # Convert source SCRIP to mbtempest + args = [ + 'mbconvert', '-B', + in_filename, + h5m_filename, + ] + check_call(args, logger) + + # Partition source SCRIP + args = [ + 'mbpart', f'{self.ntasks}', + '-z', 'RCB', + h5m_filename, + part_filename, + ] check_call(args, logger) - ds_remapped = xr.open_dataset(out_filename) + logger.info(' Done.') + + def _create_weights_tempest(self, map_filename): + """ + Create mapping weights file using TempestRemap + """ + logger = self.logger + logger.info('Create weights file') + + args = [ + 'mbtempest', '--type', '5', + '--load', f'mali.scrip.64bit.p{self.ntasks}.h5m', + '--load', f'mpaso.scrip.64bit.p{self.ntasks}.h5m', + '--file', map_filename, + '--weights', '--gnomonic', + '--boxeps', '1e-9', + ] + + run_command( + args, self.cpus_per_task, self.ntasks, + self.openmp_threads, self.config, self.logger, + ) + + logger.info(' Done.') + + def _remap_to_target(self, map_filename): + """ + Remap combined bathymetry onto MPAS target mesh + """ + logger = self.logger + logger.info('Remap to target') + + args = [ + 'ncremap', + '-m', map_filename, + '--vrb=1', + 'mali_topography_mod.nc', 'mali_topography_ncremap.nc', + ] + check_call(args, logger) + + ds_remapped = xr.open_dataset('mali_topography_ncremap.nc') ds_out = xr.Dataset() for var_name in ds_remapped: var = ds_remapped[var_name] @@ -204,6 +232,8 @@ def _remap_mali_topo(self): write_netcdf(ds_out, 'mali_topography_remapped.nc') + logger.info(' Done.') + def _combine_topo(self): os.rename('topography_remapped.nc', 'bedmachine_topography_remapped.nc')