From 64f63aeddf9bd30019d8d775e85d1028a53e1373 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Mon, 20 Jul 2026 13:39:03 -0700 Subject: [PATCH] Add FastIsostasy coupling setup for ISMIP7 AIS test case This commit adds functionality to set up coupled MALI-FastIsostasy simulations in the ISMIP7 AIS test case, following the same pattern as the existing Sea Level Model (SLM) integration. New files: - create_fastiso_mapping_files.py: Step class for generating mapping files * Includes create_ismip7_grid_file() function copied from MPAS-Tools * Creates MALI -> FastIsostasy mapping (conserve method) * Creates FastIsostasy -> MALI mapping (bilinear method) Modified files: - __init__.py: Added CreateFastIsoMappingFiles step when fastisostasy=true - set_up_experiment.py: Added FastIsostasy configuration and validation - ismip7_ais.cfg: Added fastisostasy, fastiso_res_km, and icesheet options - ismip7_ais_test.cfg: Added same configuration options Configuration options: - fastisostasy (bool): Enable FastIsostasy coupling - fastiso_res_km (int): Grid resolution in km (2, 4, 8, or 16 for AIS) - icesheet (str): Ice sheet identifier ('AIS' or 'GrIS') --- .../tests/ismip7_run/ismip7_ais/__init__.py | 21 +- .../create_fastiso_mapping_files.py | 184 ++++++++++++++++++ .../ismip7_run/ismip7_ais/ismip7_ais.cfg | 11 ++ .../ismip7_run/ismip7_ais/ismip7_ais_test.cfg | 11 ++ .../ismip7_ais/namelist.fastisostasy.template | 83 ++++++++ .../ismip7_ais/set_up_experiment.py | 40 ++++ 6 files changed, 349 insertions(+), 1 deletion(-) create mode 100644 compass/landice/tests/ismip7_run/ismip7_ais/create_fastiso_mapping_files.py create mode 100644 compass/landice/tests/ismip7_run/ismip7_ais/namelist.fastisostasy.template diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/__init__.py b/compass/landice/tests/ismip7_run/ismip7_ais/__init__.py index c6e3a41c94..ad707eb789 100644 --- a/compass/landice/tests/ismip7_run/ismip7_ais/__init__.py +++ b/compass/landice/tests/ismip7_run/ismip7_ais/__init__.py @@ -1,5 +1,8 @@ import os +from compass.landice.tests.ismip7_run.ismip7_ais.create_fastiso_mapping_files import ( # noqa + CreateFastIsoMappingFiles, +) from compass.landice.tests.ismip7_run.ismip7_ais.create_slm_mapping_files import ( # noqa CreateSlmMappingFiles, ) @@ -145,6 +148,19 @@ def configure(self): subdir=subdir)) self.steps_to_run.append('mapping_files') + # Optionally set up FastIsostasy mapping files + fastisostasy = config.getboolean('ismip7_run_ais', 'fastisostasy') + if fastisostasy: + subdir = 'fastiso_mapping_files' + if os.path.exists(os.path.join(self.work_dir, subdir)): + print(f"WARNING: {subdir} path already exists; skipping.") + else: + self.add_step( + CreateFastIsoMappingFiles(test_case=self, + name='fastiso_mapping_files', + subdir=subdir)) + self.steps_to_run.append('fastiso_mapping_files') + def run(self): """ A dummy run method @@ -154,4 +170,7 @@ def run(self): "level for this test. Please submit the job script in each " "experiment's subdirectory manually instead. " "To create Sea-Level Model mapping files, submit job script " - "or execute 'compass run' from the 'mapping_files' subdirectory.") + "or execute 'compass run' from the 'mapping_files' subdirectory. " + "To create FastIsostasy mapping files, submit job script " + "or execute 'compass run' from the 'fastiso_mapping_files' " + "subdirectory.") diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/create_fastiso_mapping_files.py b/compass/landice/tests/ismip7_run/ismip7_ais/create_fastiso_mapping_files.py new file mode 100644 index 0000000000..61e424cab9 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/create_fastiso_mapping_files.py @@ -0,0 +1,184 @@ +import shutil +import sys + +import netCDF4 +import numpy as np +from mpas_tools.logging import check_call +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from mpas_tools.scrip.from_planar import main as scrip_from_planar + +from compass.step import Step + +# AIS: polar stereographic EPSG:3031, standard parallel 71S, central meridian 0 +# Cell centres span -3,040,000 m to 3,040,000 m in both x and y. +AIS_X_MIN = -3040000 +AIS_X_MAX = 3040000 +AIS_Y_MIN = -3040000 +AIS_Y_MAX = 3040000 + + +def create_ismip7_grid_file(icesheet, res_km, output_file): + """ + Create a minimal ISMIP7 standard grid file containing only x and y + projected coordinate variables (metres). + + This function is copied from MPAS-Tools to create the regular grid + for FastIsostasy, which uses the same standard grid as ISMIP7. + + Parameters + ---------- + icesheet : str + Ice sheet domain: 'AIS' or 'GrIS'. + res_km : int + Grid resolution in kilometres. + output_file : str + Path for the output NetCDF file. + """ + res_m = int(res_km) * 1000 + if icesheet == 'AIS': + x = np.arange(AIS_X_MIN, AIS_X_MAX + res_m, res_m, dtype=float) + y = np.arange(AIS_Y_MIN, AIS_Y_MAX + res_m, res_m, dtype=float) + elif icesheet == 'GrIS': + # GrIS: polar stereographic EPSG:3413, standard parallel 70N, + # central meridian 315E + # Cell centres span x: -720,000 to 960,000 m; + # y: -3,450,000 to -570,000 m + GRIS_X_MIN = -720000 + GRIS_X_MAX = 960000 + GRIS_Y_MIN = -3450000 + GRIS_Y_MAX = -570000 + x = np.arange(GRIS_X_MIN, GRIS_X_MAX + res_m, res_m, dtype=float) + y = np.arange(GRIS_Y_MIN, GRIS_Y_MAX + res_m, res_m, dtype=float) + else: + raise ValueError(f"Unknown icesheet '{icesheet}'.") + + ds = netCDF4.Dataset(output_file, 'w') + ds.createDimension('x', len(x)) + ds.createDimension('y', len(y)) + xv = ds.createVariable('x', 'f8', ('x',)) + yv = ds.createVariable('y', 'f8', ('y',)) + xv[:] = x + yv[:] = y + xv.units = 'm' + xv.standard_name = 'projection_x_coordinate' + xv.long_name = 'x' + yv.units = 'm' + yv.standard_name = 'projection_y_coordinate' + yv.long_name = 'y' + ds.close() + print(f"Created ISMIP7 grid file: {output_file} " + f"({len(x)} x {len(y)} cells at {res_km} km)") + + +class CreateFastIsoMappingFiles(Step): + """ + A step for creating mapping files for FastIsostasy bedrock model + """ + + def __init__(self, test_case, name, subdir): + """ + Initialize step + """ + super().__init__(test_case=test_case, name=name, subdir=subdir) + + def setup(self): + print(" Setting up fastiso_mapping_files subdirectory") + + def run(self): + """ + Run this step of the test case + """ + config = self.config + logger = self.logger + section = config['ismip7_run_ais'] + fastisostasy = section.getboolean('fastisostasy') + if fastisostasy: + self._build_mapping_files(config, logger) + + def _build_mapping_files(self, config, logger): + """ + Build mapping files between the MALI mesh and the FastIsostasy grid. + + FastIsostasy uses a regular ISMIP7 grid (same as used for output), + so we create that grid and build bidirectional mapping files. + """ + section = config['ismip7_run_ais'] + init_cond_path = section.get('init_cond_path') + fastiso_res_km = section.getint('fastiso_res_km') + icesheet = section.get('icesheet') + section = config['parallel'] + ntasks = section.getint('cores_per_node') + + mali_scripfile = 'mali_scripfile.nc' + fastiso_scripfile = f'fastiso_{fastiso_res_km}km_scripfile.nc' + mali_meshfile = 'mali_meshfile.nc' + fastiso_gridfile = f'fastiso_grid_{fastiso_res_km}km.nc' + + # Create FastIsostasy ISMIP7 grid file + logger.info(f'Creating FastIsostasy grid file with ' + f'{fastiso_res_km} km resolution') + create_ismip7_grid_file(icesheet, fastiso_res_km, fastiso_gridfile) + + # Create FastIsostasy scripfile + logger.info('Creating scripfile for the FastIsostasy grid') + + # Determine projection based on ice sheet + if icesheet == 'AIS': + projection = 'ais-bedmap2' + elif icesheet == 'GrIS': + projection = 'gis-gimp' + else: + raise ValueError(f"Unknown icesheet '{icesheet}'") + + # Use mpas_tools.scrip.from_planar module directly + # Set up sys.argv to mimic command-line arguments + old_argv = sys.argv + sys.argv = ['scrip_from_planar', + '--input', fastiso_gridfile, + '--scrip', fastiso_scripfile, + '--proj', projection, + '--rank', '2'] + + try: + scrip_from_planar() + finally: + sys.argv = old_argv + + # Create MALI scripfile + logger.info('Creating scripfile for the MALI mesh') + shutil.copy(init_cond_path, mali_meshfile) + scrip_from_mpas(mali_meshfile, mali_scripfile) + + # MALI -> FastIsostasy mapping file (conserve) + logger.info('Creating MALI -> FastIsostasy grid mapfile ' + 'with conserve method') + + parallel_executable = config.get("parallel", "parallel_executable") + args = parallel_executable.split(' ') + args.extend(['-n', f'{ntasks}', + 'ESMF_RegridWeightGen', + '-s', mali_scripfile, + '-d', fastiso_scripfile, + '-w', 'mapfile_mali_to_fastiso.nc', + '-m', 'conserve', + '-i', '-64bit_offset', '--netcdf4', + '--src_regional', '--dst_regional']) + + check_call(args, logger) + + # FastIsostasy -> MALI mapping file (bilinear) + logger.info('Creating FastIsostasy -> MALI mesh mapfile ' + 'with bilinear method') + args = parallel_executable.split(' ') + args.extend(['-n', f'{ntasks}', + 'ESMF_RegridWeightGen', + '-s', fastiso_scripfile, + '-d', mali_scripfile, + '-w', 'mapfile_fastiso_to_mali.nc', + '-m', 'bilinear', + '-i', '-64bit_offset', '--netcdf4', + '--src_regional', '--dst_regional']) + + check_call(args, logger) + + logger.info('FastIsostasy mapping file creation complete') diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg index 745e4c89f7..d4c9ff620e 100644 --- a/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg +++ b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg @@ -78,5 +78,16 @@ slm_input_others = NotAvailable # Number of gauss-legendre nodes in latitude (typically multiple of 512) nglv = 2048 +# True if running coupled MALI-FastIsostasy simulation +fastisostasy = false + +# Resolution of FastIsostasy grid in kilometers (should match ISMIP7 grid) +# Valid resolutions for AIS: 2, 4, 8, 16 +# Note: for now, resolution finer than 8 km yield numerically unstable results +fastiso_res_km = 8 + +# Ice sheet identifier for grid creation ('AIS' or 'GrIS') +icesheet = AIS + [parallel] parallel_executable = srun --label --cpu-bind=cores diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg index 4d9f05318b..cc892395d3 100644 --- a/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg +++ b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg @@ -78,5 +78,16 @@ slm_input_others = NotAvailable # Number of gauss-legendre nodes in latitude (typically multiple of 512) nglv = 2048 +# True if running coupled MALI-FastIsostasy simulation +fastisostasy = false + +# Resolution of FastIsostasy grid in kilometers (should match ISMIP7 grid) +# Valid resolutions for AIS: 2, 4, 8, 16 +# Note: for now, resolution finer than 8 km yield numerically unstable results +fastiso_res_km = 8 + +# Ice sheet identifier for grid creation ('AIS' or 'GrIS') +icesheet = AIS + [parallel] parallel_executable = srun --label --cpu-bind=cores diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/namelist.fastisostasy.template b/compass/landice/tests/ismip7_run/ismip7_ais/namelist.fastisostasy.template new file mode 100644 index 0000000000..1afd229bfd --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/namelist.fastisostasy.template @@ -0,0 +1,83 @@ +&isostasy + ! FastIsostasy namelist for Antarctic Ice Sheet ISMIP7 simulations + ! Coupled with MALI via compass test case + + ! ======================================================================== + ! Main Options + ! ======================================================================== + heterogeneous_ssh = .true. ! Sea level varies spatially + interactive_sealevel = .true. ! Sea level interacts with solid Earth + include_elastic = .true. ! Include elastic deformation + correct_compression = .true. ! Account for compression effects + correct_distortion = .false. ! Account for distortion (false for stability) + method = 3 ! 3: LV-ELVA (most realistic) + dt_diagnostics = 10. ! [yr] Diagnostics timestep + pad = 0. ! [km] No padding needed + + ! ======================================================================== + ! Rheology + ! ======================================================================== + mantle = "uniform" ! Uniform mantle viscosity + lithosphere = "uniform" ! Uniform lithospheric thickness + layering = "equalized" ! Equalized layering + lumping = "freqdomain" ! Frequency domain lumping + viscosity_scaling_method = "none" + viscosity_scaling = 1.0 ! Only used if method = "scale" + + ! ======================================================================== + ! Rheology Parameters (Antarctic values) + ! ======================================================================== + nl = 2 ! Number of viscous layers + zl = 120.0, 600.0 ! [km] Layer depths + viscosities = 0.5e21, 2.0e21 ! [Pa s] Layer viscosities + ref_viscosity = 1.0e21 ! [Pa s] Reference viscosity + tau = 3000.0 ! [yr] Relaxation time + rheo_smoothing_radius = 40.0 ! [km] Smoothing radius + litho_thickness_scaling = 1.0 ! Lithosphere scaling + + ! ======================================================================== + ! Physical Constants + ! ======================================================================== + E = 66.0 ! [GPa] Young modulus + nu = 0.28 ! [-] Poisson ratiols + + rho_ice = 910.0 ! [kg/m^3] Ice density + rho_seawater = 1028.0 ! [kg/m^3] Seawater density + rho_water = 1000.0 ! [kg/m^3] Freshwater density + rho_uppermantle = 3400.0 ! [kg/m^3] Upper mantle density + rho_litho = 3200.0 ! [kg/m^3] Lithosphere density + g = 9.81 ! [m/s^2] Gravity + r_earth = 6.378e3 ! [km] Earth radius + m_earth = 5.972e24 ! [kg] Earth mass + + ! ======================================================================== + ! Input Files + ! ======================================================================== + restart = "None" ! Restart file (not used) + mask_file = "None" ! Mask file (MALI provides) + rheology_file = "None" ! Rheology file + + ! ======================================================================== + ! Time Stepping + ! ======================================================================== + dt_method = "tsit54" ! Adaptive time stepping + dt_init = 1.0 ! [yr] Initial timestep + dt_min = 1e-10 ! [yr] Minimum timestep + dt_max = 10.0 ! [yr] Maximum timestep + rtol = 1e-4 ! Relative tolerance + atol = 1e-2 ! Absolute tolerance + q = 1.0 ! Timestep adaptation exponent +/ + +&barysealevel + ! ======================================================================== + ! Barystatic Sea Level + ! ======================================================================== + method = "fastiso" ! Compute from ice mass changes + bsl_init = 0.0 ! [m] Initial sea level + sl_path = "input/sl_deglaciation.dat" + sl_name = "none" + A_ocean_pd = 3.625e14 ! [m^2] Ocean area (Goelzer et al. 2020) + A_ocean_path = "None" + restart = "None" +/ diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py b/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py index 98d59d4c40..58f88d1aab 100644 --- a/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py +++ b/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py @@ -64,6 +64,7 @@ def setup(self): # noqa: C901 reference_surface_fname = os.path.split(reference_surface_path)[-1] calving_method = section.get('calving_method') sea_level_model = section.getboolean('sea_level_model') + fastisostasy = section.getboolean('fastisostasy') exp_info = self.exp_info scenario = exp_info['scenario'] @@ -347,6 +348,34 @@ def setup(self): # noqa: C901 os.symlink(os.path.join(map_dir, map_file), os.path.join(self.work_dir, map_file)) + # FastIsostasy options + if fastisostasy: + options = { + 'config_uplift_method': "'fastisostasy'", + 'config_MALI_to_FASTISOSTASY_weights_file': + "'mapfile_mali_to_fastiso.nc'", + 'config_FASTISOSTASY_to_MALI_weights_file': + "'mapfile_fastiso_to_mali.nc'", + 'config_fastisostasy_parameter_file': + "'fastisostasy.nml'" + } + self.add_namelist_options(options=options, + out_name='namelist.landice') + + # Copy FastIsostasy namelist file + self.add_input_file( + filename='fastisostasy.nml', + target='namelist.fastisostasy.template', + package=resource_location, + copy=True) + + # Symlink mapping files + map_dir = os.path.join('..', 'fastiso_mapping_files') + for map_file in ('mapfile_mali_to_fastiso.nc', + 'mapfile_fastiso_to_mali.nc'): + os.symlink(os.path.join(map_dir, map_file), + os.path.join(self.work_dir, map_file)) + # --- Symlink restart for projections/ctrl --- if not is_historical: hist_exp = f"historical_{model}" @@ -398,5 +427,16 @@ def run(self): "Please run the 'mapping_files' step " "before proceeding.") + fastisostasy = section.getboolean('fastisostasy') + if fastisostasy: + map_dir = os.path.join('..', 'fastiso_mapping_files') + for map_file in ('mapfile_mali_to_fastiso.nc', + 'mapfile_fastiso_to_mali.nc'): + if not os.path.isfile(os.path.join(map_dir, map_file)): + sys.exit(f"ERROR: 'fastiso_mapping_files/{map_file}' " + "does not exist in workdir. " + "Please run the 'fastiso_mapping_files' step " + "before proceeding.") + run_model(step=self, namelist='namelist.landice', streams='streams.landice')