Skip to content

Commit

Permalink
Modified ocean/mesh/remap_topography to use mbtempest weights file ge…
Browse files Browse the repository at this point in the history
…neration.
  • Loading branch information
bmooremaley committed Oct 7, 2024
1 parent c4082e4 commit 0475c64
Show file tree
Hide file tree
Showing 2 changed files with 162 additions and 54 deletions.
10 changes: 9 additions & 1 deletion compass/ocean/mesh/remap_topography.cfg
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@
[remap_topography]

# the name of the topography file in the bathymetry database
topo_filename = BedMachineAntarctica_v3_and_GEBCO_2023_0.0125_degree_20240828.nc
topo_filename = BedMachineAntarctica_v3_and_GEBCO_2023_ne3000_20241004.nc
src_scrip_filename = ne3000_20241004.scrip.nc

# variable names in topo_filename
lon_var = lon
Expand Down Expand Up @@ -31,3 +32,10 @@ renorm_threshold = 0.01

# the density of land ice from MALI (kg/m^3)
ice_density = 910.0

# smoothing parameters
# no smoothing:
# expandDist = 0 [m]
# expandFactor = 1 [cell fraction]
expandDist = 0
expandFactor = 1
206 changes: 153 additions & 53 deletions compass/ocean/mesh/remap_topography.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,10 @@
import xarray as xr
from mpas_tools.cime.constants import constants
from mpas_tools.io import write_netcdf
from pyremap import LatLonGridDescriptor, MpasCellMeshDescriptor, Remapper
from mpas_tools.logging import check_call
from pyremap import MpasCellMeshDescriptor

from compass.parallel import run_command
from compass.step import Step


Expand All @@ -23,8 +25,10 @@ class RemapTopography(Step):
The name of the MPAS mesh to include in the mapping file
"""

def __init__(self, test_case, base_mesh_step, name='remap_topography',
subdir=None, mesh_name='MPAS_mesh'):
def __init__(
self, test_case, base_mesh_step, name='remap_topography', subdir=None,
mesh_name='MPAS_mesh',
):
"""
Create a new step
Expand Down Expand Up @@ -59,20 +63,30 @@ def setup(self):
"""
super().setup()
config = self.config
topo_filename = config.get('remap_topography', 'topo_filename')
section = config['remap_topography']
topo_filename = section.get('topo_filename')
src_scrip_filename = section.get('src_scrip_filename')

self.add_input_file(
filename='topography.nc',
target=topo_filename,
database='bathymetry_database')
database='bathymetry_database',
)
self.add_input_file(
filename='source.scrip.nc',
target=src_scrip_filename,
database='bathymetry_database',
)

base_path = self.base_mesh_step.path
base_filename = self.base_mesh_step.config.get(
'spherical_mesh', 'mpas_mesh_filename')
'spherical_mesh', 'mpas_mesh_filename',
)
target = os.path.join(base_path, base_filename)
self.add_input_file(filename='base_mesh.nc', work_dir_target=target)

self.ntasks = config.getint('remap_topography', 'ntasks')
self.min_tasks = config.getint('remap_topography', 'min_tasks')
self.ntasks = section.getint('ntasks')
self.min_tasks = section.getint('min_tasks')

def constrain_resources(self, available_resources):
"""
Expand All @@ -93,65 +107,151 @@ def run(self):
"""
Run this step of the test case
"""
self._create_target_scrip_file()
self._partition_scrip_file('source.scrip.nc')
self._partition_scrip_file('target.scrip.nc')
self._create_weights()
self._remap_to_target()
self._modify_remapped_bathymetry()

def _create_target_scrip_file(self):
"""
Create target SCRIP file from MPAS mesh file.
"""
logger = self.logger
logger.info('Create source SCRIP file')

config = self.config
section = config['remap_topography']
expandDist = section.getfloat('expandDist')
expandFactor = section.getfloat('expandFactor')

descriptor = MpasCellMeshDescriptor(
fileName='base_mesh.nc',
meshName=self.mesh_name,
)
descriptor.to_scrip(
'target.scrip.nc',
expandDist=expandDist,
expandFactor=expandFactor,
)

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')

# Convert source SCRIP to mbtempest
args = [
'mbconvert', '-B',
in_filename,
in_filename.replace('.nc', '.h5m'),
]
check_call(args, logger)

# Partition source SCRIP
args = [
'mbpart', f'{self.ntasks}',
'-z', 'RCB',
in_filename.replace('.nc', '.h5m'),
in_filename.replace('.nc', f'.p{self.ntasks}.h5m'),
]
check_call(args, logger)

logger.info(' Done.')

def _create_weights(self):
"""
Create mapping weights file using mbtempest
"""
logger = self.logger
logger.info('Create weights file')

args = [
'mbtempest', '--type', '5',
'--load', f'source.scrip.p{self.ntasks}.h5m',
'--load', f'target.scrip.p{self.ntasks}.h5m',
'--file', f'mapfv_source_to_target.nomask_{self.ntasks}_gnom.nc',
'--intx', 'moab_intx_source_target.h5m',
'--weights', '--verbose', '--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):
"""
Remap combined bathymetry onto MPAS target mesh
"""
logger = self.logger
parallel_executable = config.get('parallel', 'parallel_executable')

lon_var = config.get('remap_topography', 'lon_var')
lat_var = config.get('remap_topography', 'lat_var')
method = config.get('remap_topography', 'method')
renorm_threshold = config.getfloat('remap_topography',
'renorm_threshold')
ice_density = config.getfloat('remap_topography', 'ice_density')
logger.info('Remap to target')

# Build command args
# Unused options:
# -P mpas, handles some MPAS-specific index ordering, CF, etc...
# -C climatology, basically bypasses fill values
args = [
'ncremap',
'-m', f'mapfv_source_to_target.nomask_{self.ntasks}_gnom.nc',
'--vrb=1',
'topography.nc', 'topography_ncremap.nc',
]
check_call(args, logger)

logger.info(' Done.')

def _modify_remapped_bathymetry(self):
"""
Modify remapped bathymetry
"""
logger = self.logger
logger.info('Modify remapped bathymetry')

config = self.config
section = config['remap_topography']
renorm_threshold = section.getfloat('renorm_threshold')
ice_density = section.getfloat('ice_density')
ocean_density = constants['SHR_CONST_RHOSW']
g = constants['SHR_CONST_G']

in_descriptor = LatLonGridDescriptor.read(fileName='topography.nc',
lonVarName=lon_var,
latVarName=lat_var)

in_mesh_name = in_descriptor.meshName

out_mesh_name = self.mesh_name
out_descriptor = MpasCellMeshDescriptor(fileName='base_mesh.nc',
meshName=self.mesh_name)

mapping_file_name = \
f'map_{in_mesh_name}_to_{out_mesh_name}_{method}.nc'
remapper = Remapper(in_descriptor, out_descriptor, mapping_file_name)

remapper.build_mapping_file(method=method, mpiTasks=self.ntasks,
tempdir='.', logger=logger,
esmf_parallel_exec=parallel_executable)

remapper.remap_file(inFileName='topography.nc',
outFileName='topography_ncremap.nc',
logger=logger)

ds_in = xr.open_dataset('topography_ncremap.nc')
ds_in = ds_in.rename({'ncol': 'nCells'})
ds_out = xr.Dataset()
rename = {'bathymetry_var': 'bed_elevation',
'ice_thickness_var': 'landIceThkObserved',
'ice_frac_var': 'landIceFracObserved',
'grounded_ice_frac_var': 'landIceGroundedFracObserved',
'ocean_frac_var': 'oceanFracObserved',
'bathy_frac_var': 'bathyFracObserved'}

ds_out = xr.Dataset()
rename = {
'bathymetry_var': 'bed_elevation',
'ice_thickness_var': 'landIceThkObserved',
'ice_frac_var': 'landIceFracObserved',
'grounded_ice_frac_var': 'landIceGroundedFracObserved',
'ocean_frac_var': 'oceanFracObserved',
'bathy_frac_var': 'bathyFracObserved',
}
for option, out_var in rename.items():
in_var = config.get('remap_topography', option)
in_var = section.get(option)
ds_out[out_var] = ds_in[in_var]

ds_out['landIceFloatingFracObserved'] = \
ds_out.landIceFracObserved - ds_out.landIceGroundedFracObserved

# make sure fractions don't exceed 1
for var in ['landIceFracObserved', 'landIceGroundedFracObserved',
'landIceFloatingFracObserved', 'oceanFracObserved',
'bathyFracObserved']:
# Make sure fractions don't exceed 1
varNames = [
'landIceFracObserved',
'landIceGroundedFracObserved',
'landIceFloatingFracObserved',
'oceanFracObserved',
'bathyFracObserved',
]
for var in varNames:
ds_out[var] = np.minimum(ds_out[var], 1.)

# renormalize elevation variables
# Renormalize elevation variables
norm = ds_out.bathyFracObserved
valid = norm > renorm_threshold
for var in ['bed_elevation', 'landIceThkObserved']:
Expand All @@ -163,13 +263,13 @@ def run(self):
# not allowed to be thicker than the flotation thickness
thickness = np.minimum(thickness, flotation_thickness)
ds_out['landIceThkObserved'] = thickness

ds_out['landIcePressureObserved'] = ice_density * g * thickness

# compute the ice draft to be consistent with the land ice pressure
# and using E3SM's density of seawater

ds_out['landIceDraftObserved'] = \
- (ice_density / ocean_density) * thickness

write_netcdf(ds_out, 'topography_remapped.nc')

logger.info(' Done.')

0 comments on commit 0475c64

Please sign in to comment.