Skip to content

Commit

Permalink
Update code to follow remap_topography approach
Browse files Browse the repository at this point in the history
MOAB commands are updated to be more like recent development in
remap_topography and are organized into more helper methods.
  • Loading branch information
xylar committed Jan 6, 2025
1 parent d1683b1 commit b546dc0
Showing 1 changed file with 99 additions and 69 deletions.
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
import os
import pathlib

import numpy as np
import xarray as xr
Expand Down Expand Up @@ -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')

Expand All @@ -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
Expand Down Expand Up @@ -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]
Expand All @@ -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')
Expand Down

0 comments on commit b546dc0

Please sign in to comment.