From 2bb16489b73638837a703b122d8229f3833c554e Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 26 Sep 2023 15:29:08 -0700 Subject: [PATCH 01/23] Add basic structure for SLM circular icesheet test case This is largely copied from the dome smoke test and has stretches of code from the test left in place as an example. --- compass/landice/__init__.py | 2 + compass/landice/tests/slm/__init__.py | 17 ++ .../tests/slm/circ_icesheet/__init__.py | 44 ++++ .../tests/slm/circ_icesheet/circ_icesheet.cfg | 18 ++ .../tests/slm/circ_icesheet/namelist.landice | 2 + .../tests/slm/circ_icesheet/run_model.py | 72 +++++++ .../tests/slm/circ_icesheet/setup_mesh.py | 181 ++++++++++++++++ .../tests/slm/circ_icesheet/streams.landice | 29 +++ .../tests/slm/circ_icesheet/visualize.py | 200 ++++++++++++++++++ 9 files changed, 565 insertions(+) create mode 100644 compass/landice/tests/slm/__init__.py create mode 100644 compass/landice/tests/slm/circ_icesheet/__init__.py create mode 100644 compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg create mode 100644 compass/landice/tests/slm/circ_icesheet/namelist.landice create mode 100644 compass/landice/tests/slm/circ_icesheet/run_model.py create mode 100644 compass/landice/tests/slm/circ_icesheet/setup_mesh.py create mode 100644 compass/landice/tests/slm/circ_icesheet/streams.landice create mode 100644 compass/landice/tests/slm/circ_icesheet/visualize.py diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index 97b33a010c..338f5237a2 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -16,6 +16,7 @@ from compass.landice.tests.koge_bugt_s import KogeBugtS from compass.landice.tests.mesh_modifications import MeshModifications from compass.landice.tests.mismipplus import MISMIPplus +from compass.landice.tests.slm import Slm from compass.landice.tests.thwaites import Thwaites from compass.mpas_core import MpasCore @@ -49,4 +50,5 @@ def __init__(self): self.add_test_group(KogeBugtS(mpas_core=self)) self.add_test_group(MeshModifications(mpas_core=self)) self.add_test_group(MISMIPplus(mpas_core=self)) + self.add_test_group(Slm(mpas_core=self)) self.add_test_group(Thwaites(mpas_core=self)) diff --git a/compass/landice/tests/slm/__init__.py b/compass/landice/tests/slm/__init__.py new file mode 100644 index 0000000000..0ed3787e7d --- /dev/null +++ b/compass/landice/tests/slm/__init__.py @@ -0,0 +1,17 @@ +from compass.landice.tests.slm.circ_icesheet import CircIcesheetTest +from compass.testgroup import TestGroup + + +class Slm(TestGroup): + """ + A test group for Sea-Level Model test cases + """ + def __init__(self, mpas_core): + """ + mpas_core : compass.landice.Landice + the MPAS core that this test group belongs to + """ + super().__init__(mpas_core=mpas_core, name='slm') + + self.add_test_case( + CircIcesheetTest(test_group=self)) diff --git a/compass/landice/tests/slm/circ_icesheet/__init__.py b/compass/landice/tests/slm/circ_icesheet/__init__.py new file mode 100644 index 0000000000..a40c0f2f94 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/__init__.py @@ -0,0 +1,44 @@ +from compass.landice.tests.slm.circ_icesheet.run_model import RunModel +from compass.landice.tests.slm.circ_icesheet.setup_mesh import SetupMesh +from compass.landice.tests.slm.circ_icesheet.visualize import Visualize +from compass.testcase import TestCase + + +class CircIcesheetTest(TestCase): + """ + This test generates an idealized, circular ice sheet that has a + prescribed thickness evolution for testing coupling between MALI + and the Sea-Level Model. + + Attributes + ---------- + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.dome.Dome + The test group that this test case belongs to + The resolution or type of mesh of the test case + """ + name = 'circular_icesheet_test' + subdir = name + super().__init__(test_group=test_group, name=name, + subdir=subdir) + + self.add_step(SetupMesh(test_case=self)) + + self.add_step(RunModel(test_case=self, ntasks=4, openmp_threads=1, + name='run_step')) + + step = Visualize(test_case=self) + self.add_step(step, run_by_default=True) + + # no configure() method is needed because we will use the default + # config options + + # no run() method is needed because we're doing the default: running all + # steps diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg new file mode 100644 index 0000000000..5f86d01853 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -0,0 +1,18 @@ +# config options for slm circ_icesheet test case +[circ_icesheet] + +# resolution (in m) for the circular ice sheet +dc = 2000.0 + +# Whether to center the dome in the center of the cell that is closest to the +# center of the domain +put_origin_on_a_cell = True + +# config options related to visualization +[circ_icesheet_viz] + +# whether to save image files +save_images = True + +# whether to hide figures (typically when save_images = True) +hide_figs = True diff --git a/compass/landice/tests/slm/circ_icesheet/namelist.landice b/compass/landice/tests/slm/circ_icesheet/namelist.landice new file mode 100644 index 0000000000..13c97da547 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/namelist.landice @@ -0,0 +1,2 @@ +config_dt = '0001-00-00_00:00:00' +config_run_duration = '0002-00-00_00:00:00' diff --git a/compass/landice/tests/slm/circ_icesheet/run_model.py b/compass/landice/tests/slm/circ_icesheet/run_model.py new file mode 100644 index 0000000000..d115d63784 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/run_model.py @@ -0,0 +1,72 @@ +from compass.model import run_model +from compass.step import Step + + +class RunModel(Step): + """ + A step for performing forward MALI runs as part of dome test cases. + + Attributes + ---------- + """ + def __init__(self, test_case, name='run_model', + subdir=None, ntasks=1, min_tasks=None, openmp_threads=1): + """ + Create a new test case + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + name : str, optional + the name of the test case + + subdir : str, optional + the subdirectory for the step. The default is ``name`` + + ntasks : int, optional + the number of tasks the step would ideally use. If fewer tasks + are available on the system, the step will run on all available + tasks as long as this is not below ``min_tasks`` + + min_tasks : int, optional + the number of tasks the step requires. If the system has fewer + than this number of tasks, the step will fail + + openmp_threads : int, optional + the number of OpenMP threads the step will use + """ + if min_tasks is None: + min_tasks = ntasks + super().__init__(test_case=test_case, name=name, subdir=subdir, + ntasks=ntasks, min_tasks=min_tasks, + openmp_threads=openmp_threads) + + self.add_namelist_file( + 'compass.landice.tests.slm.circ_icesheet', 'namelist.landice', + out_name='namelist.landice') + + self.add_streams_file( + 'compass.landice.tests.slm.circ_icesheet', 'streams.landice', + out_name='streams.landice') + + self.add_input_file(filename='landice_grid.nc', + target='../setup_mesh/landice_grid.nc') + self.add_input_file(filename='graph.info', + target='../setup_mesh/graph.info') + self.add_input_file(filename='namelist.sealevel', + package='compass.landice.tests.slm', + copy=True) + + self.add_model_as_input() + + self.add_output_file(filename='output.nc') + + # no setup() is needed + + def run(self): + """ + Run this step of the test case + """ + run_model(step=self) diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py new file mode 100644 index 0000000000..0af42971c4 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -0,0 +1,181 @@ +import numpy +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from mpas_tools.mesh.conversion import convert, cull +from mpas_tools.planar_hex import make_planar_hex_mesh +from netCDF4 import Dataset as NetCDFFile + +from compass.model import make_graph_file +from compass.step import Step + + +class SetupMesh(Step): + """ + A step for creating a mesh with initial conditions plus forcing file for + circular icesheet test + + Attributes + ---------- + """ + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + """ + super().__init__(test_case=test_case, name='setup_mesh') + + self.add_output_file(filename='graph.info') + self.add_output_file(filename='landice_grid.nc') + + # no setup() method is needed + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + section = config['dome'] + + nx = section.getint('nx') + ny = section.getint('ny') + dc = section.getfloat('dc') + + dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, + nonperiodic_x=True, + nonperiodic_y=True) + dsMesh = cull(dsMesh, logger=logger) + dsMesh = convert(dsMesh, logger=logger) + write_netcdf(dsMesh, 'mpas_grid.nc') + + levels = 3 + args = ['create_landice_grid_from_generic_MPAS_grid.py', + '-i', 'mpas_grid.nc', + '-o', 'landice_grid.nc', + '-l', levels] + + check_call(args, logger) + + make_graph_file(mesh_filename='landice_grid.nc', + graph_filename='graph.info') + + _setup_circsheet_initial_conditions(config, logger, + filename='landice_grid.nc') + + +def _setup_circsheet_initial_conditions(config, logger, filename): + """ + Create initial condition for circular ice sheet + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, a combination of the defaults + for the machine, core and configuration + + logger : logging.Logger + A logger for output from the step + + filename : str + file to setup dome + """ + section = config['dome'] + dome_type = section.get('dome_type') + put_origin_on_a_cell = section.getboolean('put_origin_on_a_cell') + shelf = section.getboolean('shelf') + hydro = section.getboolean('hyrdo') + + # Open the file, get needed dimensions + gridfile = NetCDFFile(filename, 'r+') + nVertLevels = len(gridfile.dimensions['nVertLevels']) + # Get variables + xCell = gridfile.variables['xCell'] + yCell = gridfile.variables['yCell'] + xEdge = gridfile.variables['xEdge'] + yEdge = gridfile.variables['yEdge'] + xVertex = gridfile.variables['xVertex'] + yVertex = gridfile.variables['yVertex'] + thickness = gridfile.variables['thickness'] + bedTopography = gridfile.variables['bedTopography'] + layerThicknessFractions = gridfile.variables['layerThicknessFractions'] + SMB = gridfile.variables['sfcMassBal'] + + # Find center of domain + x0 = xCell[:].min() + 0.5 * (xCell[:].max() - xCell[:].min()) + y0 = yCell[:].min() + 0.5 * (yCell[:].max() - yCell[:].min()) + # Calculate distance of each cell center from dome center + r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 + + if put_origin_on_a_cell: + # Center the dome in the center of the cell that is closest to the + # center of the domain. + centerCellIndex = numpy.abs(r[:]).argmin() + xShift = -1.0 * xCell[centerCellIndex] + yShift = -1.0 * yCell[centerCellIndex] + xCell[:] = xCell[:] + xShift + yCell[:] = yCell[:] + yShift + xEdge[:] = xEdge[:] + xShift + yEdge[:] = yEdge[:] + yShift + xVertex[:] = xVertex[:] + xShift + yVertex[:] = yVertex[:] + yShift + # Now update origin location and distance array + x0 = 0.0 + y0 = 0.0 + r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 + + # Assign variable values for dome + # Define dome dimensions - all in meters + r0 = 60000.0 * numpy.sqrt(0.125) + h0 = 2000.0 * numpy.sqrt(0.125) + # Set default value for non-dome cells + thickness[:] = 0.0 + # Calculate the dome thickness for cells within the desired radius + # (thickness will be NaN otherwise) + thickness_field = thickness[0, :] + if dome_type == 'cism': + thickness_field[r < r0] = h0 * (1.0 - (r[r < r0] / r0) ** 2) ** 0.5 + elif dome_type == 'halfar': + thickness_field[r < r0] = h0 * ( + 1.0 - (r[r < r0] / r0) ** (4.0 / 3.0)) ** (3.0 / 7.0) + else: + raise ValueError('Unexpected dome_type: {}'.format(dome_type)) + thickness[0, :] = thickness_field + + # zero velocity everywhere + # normalVelocity[:] = 0.0 + # flat bed at sea level + bedTopography[:] = 0.0 + if shelf: + # this line will make a small shelf: + bedTopography[0, xCell[:] < -10000.0] = -2000.0 + # Setup layerThicknessFractions + layerThicknessFractions[:] = 1.0 / nVertLevels + + # boundary conditions + # Sample values to use, or comment these out for them to be 0. + SMB[:] = 0.0 + # beta[:] = 50000. + # units: m/yr, lapse rate of 1 m/yr with 0 at 500 m + # SMB[:] = 2.0/1000.0 * (thickness[:] + bedTopography[:]) - 1.0 + # Convert from units of m/yr to kg/m2/s using an assumed ice density + SMB[:] = SMB[:] * 910.0 / (3600.0 * 24.0 * 365.0) + + # lapse rate of 5 deg / km + # Tsfc[:, 0] = -5.0/1000.0 * (thickness[0,:] + bedTopography[0,:]) + # G = 0.01 + # BMB[:] = -20.0 # units: m/yr + + if hydro: + gridfile.variables['uReconstructX'][:] = 5.0 / (3600.0 * 24.0 * 365.0) + gridfile.variables['basalMeltInput'][:] = 0.06 / 335000.0 * 50.0 + gridfile.variables['externalWaterInput'][:] = 0.0 + gridfile.variables['waterThickness'][:] = 0.08 + + gridfile.close() + + logger.info('Successfully added dome initial conditions to: {}'.format( + filename)) diff --git a/compass/landice/tests/slm/circ_icesheet/streams.landice b/compass/landice/tests/slm/circ_icesheet/streams.landice new file mode 100644 index 0000000000..3f0648ada8 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/streams.landice @@ -0,0 +1,29 @@ + + + + + + + + + + + + + + + + + diff --git a/compass/landice/tests/slm/circ_icesheet/visualize.py b/compass/landice/tests/slm/circ_icesheet/visualize.py new file mode 100644 index 0000000000..30ea21d857 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/visualize.py @@ -0,0 +1,200 @@ +import matplotlib.pyplot as plt +import netCDF4 +import numpy + +from compass.step import Step + + +class Visualize(Step): + """ + A step for visualizing the output from a dome test case + + Attributes + ---------- + """ + def __init__(self, test_case, name='visualize', subdir=None, + input_dir='run_model'): + """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + name : str, optional + the name of the test case + + subdir : str, optional + the subdirectory for the step. The default is ``name`` + + input_dir : str, optional + The input directory within the test case with a file ``output.nc`` + to visualize + """ + super().__init__(test_case=test_case, name=name, subdir=subdir) + + self.add_input_file(filename='output.nc', + target='../{}/output.nc'.format(input_dir)) + + # depending on settings, this may produce no outputs, so we won't add + # any + + # no setup method is needed + + def run(self): + """ + Run this step of the test case + """ + visualize_slm_circsheet(self.config, self.logger, filename='output.nc') + + +def visualize_slm_circsheet(config, logger, filename): + """ + Plot the output from a dome test case + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, a combination of the defaults + for the machine, core and configuration + + logger : logging.Logger + A logger for output from the step + + filename : str + file to visualize + """ + section = config['dome_viz'] + + time_slice = section.getint('time_slice') + save_images = section.getboolean('save_images') + hide_figs = section.getboolean('hide_figs') + + # Note: this may be slightly wrong for some calendar types! + secInYr = 3600.0 * 24.0 * 365.0 + + f = netCDF4.Dataset(filename, 'r') + + times = f.variables['xtime'] + thickness = f.variables['thickness'] + # dcEdge = f.variables['dcEdge'] + # bedTopography = f.variables['bedTopography'] # not needed + xCell = f.variables['xCell'] + yCell = f.variables['yCell'] + xEdge = f.variables['xEdge'] + yEdge = f.variables['yEdge'] + angleEdge = f.variables['angleEdge'] + temperature = f.variables['temperature'] + lowerSurface = f.variables['lowerSurface'] + upperSurface = f.variables['upperSurface'] + normalVelocity = f.variables['normalVelocity'] + # uReconstructX = f.variables['uReconstructX'] + uReconstructX = f.variables['uReconstructX'] + uReconstructY = f.variables['uReconstructY'] + + vert_levs = len(f.dimensions['nVertLevels']) + + time_length = times.shape[0] + + logger.info("vert_levs = {}; time_length = {}".format(vert_levs, + time_length)) + + var_slice = thickness[time_slice, :] + + fig = plt.figure(1, facecolor='w') + fig.add_subplot(111, aspect='equal') + # C = plt.contourf(xCell, yCell, var_slice ) + plt.scatter(xCell[:], yCell[:], 80, var_slice, marker='h', + edgecolors='none') + plt.colorbar() + plt.title('thickness at time {}'.format(time_slice)) + plt.draw() + if save_images: + logger.info("Saving figures to files.") + plt.savefig('dome_thickness.png') + + fig = plt.figure(2) + fig.add_subplot(121, aspect='equal') + plt.scatter(xCell[:], yCell[:], 80, lowerSurface[time_slice, :], + marker='h', edgecolors='none') + plt.colorbar() + plt.title('lower surface at time {}'.format(time_slice)) + plt.draw() + fig.add_subplot(122, aspect='equal') + plt.scatter(xCell[:], yCell[:], 80, upperSurface[time_slice, :], + marker='h', edgecolors='none') + plt.colorbar() + plt.title('upper surface at time {}'.format(time_slice)) + plt.draw() + if save_images: + plt.savefig('dome_surfaces.png') + + fig = plt.figure(3) + for templevel in range(0, vert_levs): + fig.add_subplot(3, 4, templevel + 1, aspect='equal') + var_slice = temperature[time_slice, :, templevel] + # C = plt.contourf(xCell, yCell, var_slice ) + plt.scatter(xCell[:], yCell[:], 40, var_slice, marker='h', + edgecolors='none') + plt.colorbar() + plt.title('temperature at level {} at time {}'.format(templevel, + time_slice)) + plt.draw() + if save_images: + plt.savefig('dome_temperature.png') + + fig = plt.figure(4) + fig.add_subplot(121, aspect='equal') + plt.scatter(xEdge[:], yEdge[:], 80, + normalVelocity[time_slice, :, vert_levs - 1] * secInYr, + marker='h', edgecolors='none') + plt.colorbar() + normalVel = normalVelocity[time_slice, :, vert_levs - 1] + plt.quiver(xEdge[:], yEdge[:], + numpy.cos(angleEdge[:]) * normalVel * secInYr, + numpy.sin(angleEdge[:]) * normalVel * secInYr) + plt.title('normalVelocity of bottom layer at time {}'.format(time_slice)) + plt.draw() + fig.add_subplot(122, aspect='equal') + plt.scatter(xEdge[:], yEdge[:], 80, + normalVelocity[time_slice, :, 0] * secInYr, marker='h', + edgecolors='none') + plt.colorbar() + normalVel = normalVelocity[time_slice, :, 0] + plt.quiver(xEdge[:], yEdge[:], + numpy.cos(angleEdge[:]) * normalVel * secInYr, + numpy.sin(angleEdge[:]) * normalVel * secInYr) + plt.title('normalVelocity of top layer at time {}'.format(time_slice)) + plt.draw() + if save_images: + plt.savefig('dome_normalVelocity.png') + + fig = plt.figure(5, facecolor='w') + fig.add_subplot(121, aspect='equal') + plt.scatter(xCell[:], yCell[:], 80, + uReconstructX[time_slice, :, 0] * secInYr, marker='h', + edgecolors='none') + plt.colorbar() + plt.quiver(xCell[:], yCell[:], uReconstructX[time_slice, :, 0] * secInYr, + uReconstructY[time_slice, :, 0] * secInYr) + plt.title('uReconstructX of top layer at time {}'.format(time_slice)) + plt.draw() + fig.add_subplot(122, aspect='equal') + plt.scatter(xCell[:], yCell[:], 80, + uReconstructY[time_slice, :, 0] * secInYr, marker='h', + edgecolors='none') + plt.colorbar() + plt.quiver(xCell[:], yCell[:], uReconstructX[time_slice, :, 0] * secInYr, + uReconstructY[time_slice, :, 0] * secInYr) + plt.title('uReconstructY of top layer at time {}'.format(time_slice)) + plt.draw() + if save_images: + plt.savefig('dome_uReconstruct.png') + + if hide_figs: + logger.info("Plot display disabled with hide_plot config option.") + else: + plt.show() + + f.close() From 161eb9c2763a5c34503ad585426017c70fd0a8e9 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Thu, 28 Sep 2023 09:41:44 -0600 Subject: [PATCH 02/23] Edit config options --- .../tests/slm/circ_icesheet/circ_icesheet.cfg | 47 +++++++++++++++++-- 1 file changed, 44 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg index 5f86d01853..52717efdde 100644 --- a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -1,16 +1,57 @@ # config options for slm circ_icesheet test case [circ_icesheet] +# The number of cells in the x,y directions for the planar uniform mesh +# for 3000km mesh, need to have 150 for dc=20km +nx = 150 +ny = 150 + # resolution (in m) for the circular ice sheet -dc = 2000.0 +# i.e. The distance in meters between adjacent cell centers. +dc = 20000.0 + +# the ice shape type ('cylinder', and for dome, 'dome-halfar' or 'dome-cism') +# r0 and h0: initial radius and height of the cylinder ice in meters +ice_type = cylinder +r0 = 15000.0 +h0 = 3000.0 + +# 'True' if manually want to set bedTopography elevation +# bed topography elevation in meters (provide a positive/negative float value +# for grounded above land/marine-based ice) +set_topo_elev = True +topo_elev = -1500.0 -# Whether to center the dome in the center of the cell that is closest to the +# Whether to center the circular ice in the center of the cell that is closest to the # center of the domain -put_origin_on_a_cell = True +put_origin_on_a_cell = False + +# config options for surface mass balance forcing +[smb_forcing] + +# direction in which SMB is applied ('horizontal' or 'vertical') +direction = horizontal + +# config options for the sea-level model +[slm] +# Gauss-Legendre degree +slm_res = 512 + +# mapping method between the MALI and SLM grids +mapping_method_mali_to_slm = conserve +mapping_method_slm_to_mali = bilinear # config options related to visualization [circ_icesheet_viz] +# Area (m^2) of the global ocean for calculating ice sheet contribution to +# sea level. Only used when MALI-SLM coupling is False +Aocn_const = 4.5007E+14 + +# Area (m^2) of the global ocean excluding the marine-based ice region +# only used when MALI-SLM coupling is False +AocnBeta_const = 4.50007E+14 + # whether to save image files save_images = True From c151cd50218b569ad76f7dcb14fab9e06648cff1 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Thu, 28 Sep 2023 14:00:34 -0600 Subject: [PATCH 03/23] set up mesh, smb focring and mapping files --- .../tests/slm/circ_icesheet/setup_mesh.py | 294 +++++++++++++++--- 1 file changed, 248 insertions(+), 46 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index 0af42971c4..c94082cc2c 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -1,8 +1,12 @@ -import numpy +import os +import shutil + +import numpy as np from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call from mpas_tools.mesh.conversion import convert, cull from mpas_tools.planar_hex import make_planar_hex_mesh +from mpas_tools.scrip.from_mpas import scrip_from_mpas from netCDF4 import Dataset as NetCDFFile from compass.model import make_graph_file @@ -17,6 +21,7 @@ class SetupMesh(Step): Attributes ---------- """ + def __init__(self, test_case): """ Create the step @@ -39,12 +44,12 @@ def run(self): """ logger = self.logger config = self.config - section = config['dome'] + section = config['circ_icesheet'] nx = section.getint('nx') ny = section.getint('ny') dc = section.getfloat('dc') - + print('calling the mesh creation function') dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=True, nonperiodic_y=True) @@ -56,7 +61,7 @@ def run(self): args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', 'mpas_grid.nc', '-o', 'landice_grid.nc', - '-l', levels] + '-l', str(levels)] check_call(args, logger) @@ -66,6 +71,13 @@ def run(self): _setup_circsheet_initial_conditions(config, logger, filename='landice_grid.nc') + _create_smb_forcing_file(config, logger, + mali_mesh_file='landice_grid.nc', + filename='smb_forcing.nc') + logger.info('calling building mapping files') + _build_mapping_files(config, logger, + mali_mesh_file='landice_grid.nc') + def _setup_circsheet_initial_conditions(config, logger, filename): """ @@ -81,13 +93,15 @@ def _setup_circsheet_initial_conditions(config, logger, filename): A logger for output from the step filename : str - file to setup dome + file to setup circ_icesheet """ - section = config['dome'] - dome_type = section.get('dome_type') + section = config['circ_icesheet'] + ice_type = section.get('ice_type') + set_topo_elev = section.getboolean('set_topo_elev') + topo_elev = section.get('topo_elev') + r0 = section.get('r0') + h0 = section.get('h0') put_origin_on_a_cell = section.getboolean('put_origin_on_a_cell') - shelf = section.getboolean('shelf') - hydro = section.getboolean('hyrdo') # Open the file, get needed dimensions gridfile = NetCDFFile(filename, 'r+') @@ -102,7 +116,6 @@ def _setup_circsheet_initial_conditions(config, logger, filename): thickness = gridfile.variables['thickness'] bedTopography = gridfile.variables['bedTopography'] layerThicknessFractions = gridfile.variables['layerThicknessFractions'] - SMB = gridfile.variables['sfcMassBal'] # Find center of domain x0 = xCell[:].min() + 0.5 * (xCell[:].max() - xCell[:].min()) @@ -111,9 +124,9 @@ def _setup_circsheet_initial_conditions(config, logger, filename): r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 if put_origin_on_a_cell: - # Center the dome in the center of the cell that is closest to the + # Center the ice in the center of the cell that is closest to the # center of the domain. - centerCellIndex = numpy.abs(r[:]).argmin() + centerCellIndex = np.abs(r[:]).argmin() xShift = -1.0 * xCell[centerCellIndex] yShift = -1.0 * yCell[centerCellIndex] xCell[:] = xCell[:] + xShift @@ -126,56 +139,245 @@ def _setup_circsheet_initial_conditions(config, logger, filename): x0 = 0.0 y0 = 0.0 r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 - - # Assign variable values for dome - # Define dome dimensions - all in meters - r0 = 60000.0 * numpy.sqrt(0.125) - h0 = 2000.0 * numpy.sqrt(0.125) - # Set default value for non-dome cells + print('HH==max r is', max(r)) + print('HH==min r is', min(r)) + # Assign variable values for the circular ice sheet + # Set default value for non-circular cells thickness[:] = 0.0 # Calculate the dome thickness for cells within the desired radius # (thickness will be NaN otherwise) thickness_field = thickness[0, :] - if dome_type == 'cism': + r0 = 6000000.0 * np.sqrt(0.125) + + print('r0', r0) + print('r', r) + if ice_type == 'cylinder': + logger.info('cylinder ice type is chosen') + thickness_field[r < r0] = h0 + elif ice_type == 'dome-cism': thickness_field[r < r0] = h0 * (1.0 - (r[r < r0] / r0) ** 2) ** 0.5 - elif dome_type == 'halfar': + elif ice_type == 'dome-halfar': thickness_field[r < r0] = h0 * ( 1.0 - (r[r < r0] / r0) ** (4.0 / 3.0)) ** (3.0 / 7.0) else: - raise ValueError('Unexpected dome_type: {}'.format(dome_type)) + raise ValueError('Unexpected ice_type: {}'.format(ice_type)) thickness[0, :] = thickness_field - # zero velocity everywhere - # normalVelocity[:] = 0.0 # flat bed at sea level bedTopography[:] = 0.0 - if shelf: + if set_topo_elev: # this line will make a small shelf: - bedTopography[0, xCell[:] < -10000.0] = -2000.0 + bedTopography[:] = topo_elev + # Setup layerThicknessFractions layerThicknessFractions[:] = 1.0 / nVertLevels - # boundary conditions - # Sample values to use, or comment these out for them to be 0. - SMB[:] = 0.0 - # beta[:] = 50000. - # units: m/yr, lapse rate of 1 m/yr with 0 at 500 m - # SMB[:] = 2.0/1000.0 * (thickness[:] + bedTopography[:]) - 1.0 - # Convert from units of m/yr to kg/m2/s using an assumed ice density - SMB[:] = SMB[:] * 910.0 / (3600.0 * 24.0 * 365.0) - - # lapse rate of 5 deg / km - # Tsfc[:, 0] = -5.0/1000.0 * (thickness[0,:] + bedTopography[0,:]) - # G = 0.01 - # BMB[:] = -20.0 # units: m/yr - - if hydro: - gridfile.variables['uReconstructX'][:] = 5.0 / (3600.0 * 24.0 * 365.0) - gridfile.variables['basalMeltInput'][:] = 0.06 / 335000.0 * 50.0 - gridfile.variables['externalWaterInput'][:] = 0.0 - gridfile.variables['waterThickness'][:] = 0.08 - gridfile.close() - logger.info('Successfully added dome initial conditions to: {}'.format( + logger.info('Successfully added circular initial conditions to: {}'.format( filename)) + + +def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): + """ + Create surface mass balance forcing file for circular ice sheet + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, a combination of the defaults + for the machine, core and configuration + + logger : logging.Logger + A logger for output from the step + + mali_mesh_file : str + mesh file created in the current step + + filename : str + file to setup circ_icesheet + """ + section = config['circ_icesheet'] + r0 = section.get('r0') + h0 = section.get('h0') + r0 = 6000000.0 * np.sqrt(0.125) # HH: comment this line out later + + section = config['smb_forcing'] + direction = section.get('direction') + + # other constants used in the fuction + pi = 3.1415926535898 # pi value used in the SLM + rhoi = 910.0 # ice density + + # Open the file, get needed dimensions + shutil.copy(mali_mesh_file, filename) + smbfile = NetCDFFile(filename, 'r+') + + # Get variables + xCell = smbfile.variables['xCell'] + yCell = smbfile.variables['yCell'] + # nCells = smbfile.dimensions['nCells'] + sfcMassBal = smbfile.variables['sfcMassBal'] + + # set xtime variable + # smbfile.variables['xtime'] + + # initialize SMB value everywhere in the mesh + sfcMassBal[:, :] = 0.0 + + # define the time (1 year interval) + dt = 1 + t_array = np.arange(0, 10, dt) + x0 = 0.0 + y0 = 0.0 + r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 + + dVdt = -9.678E13 + dHdt = 20.0 + dRdt = 40000.0 + + if direction == 'vertical': + print('creating a file that prescribes vertical smb') + # dHdt = dVdt / (pi*(r0*1000.0)**2) # change in ice thickness in m + smb = dHdt * rhoi / (365.0 * 24 * 60 * 60) # smb in kg/m^2/s + indx = np.where(r <= r0)[0] + + # for loop through time + for t in t_array: + print((t)) + Ht = h0 + (dHdt * t * dt) # calculate a new height each time + print(f'new height will be: {Ht}km') + sfcMassBal[t, indx] = smb # assign sfcMassBal + smbfile.variables['sfcMassBal'][t, :] = sfcMassBal[t, :] + + elif direction == 'horizontal': + print('creating a file that prescribes horizontal smb') + # dRdt = (abs(Rf - r0)) / len(t_array) + smb = (10000.0) * rhoi / (365.0 * 24.0 * 60.0 * 60.0) + print('-smb is', -smb) + # for loop through time + for t in t_array: + print((t)) + # Rt = abs(np.sqrt(((dVdt * t * dt) / (pi * h0)) + (r0*1000)**2)) + # calculate a new radius at each time; for constant radius change + Rt = (r0 - (dRdt * t * dt)) + if (np.isnan(Rt) or Rt < 0): + Rt = 0.0 + smbfile.variables['sfcMassBal'][t, :] = 0.0 + print(f'new radius is: {Rt}km') + else: + indx = np.where(r >= Rt / 1000)[0] + print('idx', indx) + smbfile.variables['sfcMassBal'][t, indx] = -smb + print(f'new radius is: {Rt / 1000}km') + smbfile.variables['sfcMassBal'][t, :] = sfcMassBal[t, :] + + elif direction == 'dome-halfar': + for t in t_array: + Rt = (((dVdt * dt * t) * 3 / (2 * pi) + + (r0 * 1000) ** 3)) ** (1.0 / 3) + print(f'new radius is: {Rt / 1000}km') + # smb = (dVdt * t * dt) / (2 * pi * (Rt ** 2)) \ + # * rhoi / (365.0 * 24.0 * 60.0 * 60.0)i + indx = np.where(r > r0)[0] + smbfile.variables['sfcMassBal'][t, :] = smb + print(f'smb is: {smb}') + print('ice thickness will change both horizontally and vertically') + + smbfile.close() + + +def _build_mapping_files(self, config, logger, mali_mesh_file): + """ + Build a mapping file if it does not exist. + Mapping file is then used to remap the ismip6 source file in polarstero + coordinate to unstructured mali mesh + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for a ismip6 forcing test case + + cores : int + the number of cores for the ESMF_RegridWeightGen + + logger : logging.Logger + A logger for output from the step + + mali_mesh_file : str, optional + The MALI mesh file if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + + section = config['slm'] + slm_res = section.get['slm_res'] + method_mali_to_slm = section.get['mapping_method_mali_to_slm'] + method_slm_to_mali = section.get['mapping_method_slm_to_mali'] + + section = config['circ_icesheet'] + dc = section.get['dc'] + + mali_scripfile = 'mali_scrip.nc' + slm_scripfile = 'slm_scrip.nc' + + # create slm scripfile + logger.info(f'creating scripfile for the SH-degree {slm_res}' + f'resolution SLM grid') + args = ['ncremap', + '-g', slm_scripfile, + '-G', + f'latlon={slm_res},{int(2*slm_res)}#lat_typ=gss#lat_drc=n2s'] + + check_call(args, logger=self.logger) + + # create scrip files for source and destination grids + logger.info('creating scrip file for the mali mesh') + mali_mesh_copy = f'{mali_mesh_file}_copy' + shutil.copy(mali_mesh_file, mali_mesh_copy) + + args = ['set_lat_lon_fields_in_planar_grid.py', + '--file', mali_mesh_copy, + '--proj', 'ais-bedmap2-sphere'] + + check_call(args, logger=self.logger) + + scrip_from_mpas(mali_mesh_copy, mali_scripfile) + + # create a mapping file using ESMF weight gen + logger.info('Creating a mapping file...') + + mapping_file_mali_to_slm = f'mapping_file_mali{int(dc/1000)}_to' \ + f'_slm{slm_res}_{method_mali_to_slm}.nc' + mapping_file_slm_to_mali = f'mapping_file_slm{slm_res}_to' \ + f'_mali{int(dc/1000)}_{method_slm_to_mali}.nc' + + parallel_executable = config.get("parallel", "parallel_executable") + # split the parallel executable into constituents in case it includes flags + args = parallel_executable.split(' ') + args.extend(['ESMF_RegridWeightGen', + '-s', mali_scripfile, + '-d', slm_scripfile, + '-w', mapping_file_mali_to_slm, + '-m', method_mali_to_slm, + '-i', '-64bit_offset', '--netcdf4', + '--src_regional']) + + check_call(args, logger) + + args.extend(['ESMF_RegridWeightGen', + '-s', slm_scripfile, + '-d', mali_scripfile, + '-w', mapping_file_slm_to_mali, + '-m', method_slm_to_mali, + '-i', '-64bit_offset', '--netcdf4', + '--dst_regional']) + + check_call(args, logger) + + # remove the temporary scripfiles once the mapping file is generated + logger.info('Removing the temporary mesh and scripfiles...') + os.remove(slm_scripfile) + os.remove(mali_scripfile) + os.remove(mali_mesh_copy) From 399fab27e11f65bda0ef6f7bc903cdc32d4d12f7 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Thu, 28 Sep 2023 16:55:34 -0600 Subject: [PATCH 04/23] Clean up the `setup mesh` step --- .../tests/slm/circ_icesheet/setup_mesh.py | 46 ++++++++----------- 1 file changed, 20 insertions(+), 26 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index c94082cc2c..7e33521704 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -74,7 +74,7 @@ def run(self): _create_smb_forcing_file(config, logger, mali_mesh_file='landice_grid.nc', filename='smb_forcing.nc') - logger.info('calling building mapping files') + _build_mapping_files(config, logger, mali_mesh_file='landice_grid.nc') @@ -99,8 +99,8 @@ def _setup_circsheet_initial_conditions(config, logger, filename): ice_type = section.get('ice_type') set_topo_elev = section.getboolean('set_topo_elev') topo_elev = section.get('topo_elev') - r0 = section.get('r0') - h0 = section.get('h0') + r0 = section.getfloat('r0') + h0 = section.getfloat('h0') put_origin_on_a_cell = section.getboolean('put_origin_on_a_cell') # Open the file, get needed dimensions @@ -147,7 +147,6 @@ def _setup_circsheet_initial_conditions(config, logger, filename): # Calculate the dome thickness for cells within the desired radius # (thickness will be NaN otherwise) thickness_field = thickness[0, :] - r0 = 6000000.0 * np.sqrt(0.125) print('r0', r0) print('r', r) @@ -198,9 +197,8 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): file to setup circ_icesheet """ section = config['circ_icesheet'] - r0 = section.get('r0') - h0 = section.get('h0') - r0 = 6000000.0 * np.sqrt(0.125) # HH: comment this line out later + r0 = section.getfloat('r0') + h0 = section.getfloat('h0') section = config['smb_forcing'] direction = section.get('direction') @@ -227,7 +225,7 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): # define the time (1 year interval) dt = 1 - t_array = np.arange(0, 10, dt) + t_array = np.arange(0, 3, dt) x0 = 0.0 y0 = 0.0 r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 @@ -287,7 +285,7 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): smbfile.close() -def _build_mapping_files(self, config, logger, mali_mesh_file): +def _build_mapping_files(config, logger, mali_mesh_file): """ Build a mapping file if it does not exist. Mapping file is then used to remap the ismip6 source file in polarstero @@ -312,25 +310,25 @@ def _build_mapping_files(self, config, logger, mali_mesh_file): """ section = config['slm'] - slm_res = section.get['slm_res'] - method_mali_to_slm = section.get['mapping_method_mali_to_slm'] - method_slm_to_mali = section.get['mapping_method_slm_to_mali'] + slm_res = section.get('slm_res') + method_mali_to_slm = section.get('mapping_method_mali_to_slm') + method_slm_to_mali = section.get('mapping_method_slm_to_mali') section = config['circ_icesheet'] - dc = section.get['dc'] + dc = section.getfloat('dc') - mali_scripfile = 'mali_scrip.nc' - slm_scripfile = 'slm_scrip.nc' + mali_scripfile = f'mali{int(dc/1000)}km_scripfile.nc' + slm_scripfile = f'slm{slm_res}_scripfile.nc' # create slm scripfile - logger.info(f'creating scripfile for the SH-degree {slm_res}' + logger.info(f'creating scripfile for the SH-degree {slm_res} ' f'resolution SLM grid') args = ['ncremap', '-g', slm_scripfile, '-G', - f'latlon={slm_res},{int(2*slm_res)}#lat_typ=gss#lat_drc=n2s'] + f'latlon={slm_res},{2*int(slm_res)}#lat_typ=gss#lat_drc=n2s'] - check_call(args, logger=self.logger) + check_call(args, logger=logger) # create scrip files for source and destination grids logger.info('creating scrip file for the mali mesh') @@ -341,35 +339,31 @@ def _build_mapping_files(self, config, logger, mali_mesh_file): '--file', mali_mesh_copy, '--proj', 'ais-bedmap2-sphere'] - check_call(args, logger=self.logger) + check_call(args, logger=logger) scrip_from_mpas(mali_mesh_copy, mali_scripfile) # create a mapping file using ESMF weight gen logger.info('Creating a mapping file...') - mapping_file_mali_to_slm = f'mapping_file_mali{int(dc/1000)}_to' \ - f'_slm{slm_res}_{method_mali_to_slm}.nc' - mapping_file_slm_to_mali = f'mapping_file_slm{slm_res}_to' \ - f'_mali{int(dc/1000)}_{method_slm_to_mali}.nc' - parallel_executable = config.get("parallel", "parallel_executable") # split the parallel executable into constituents in case it includes flags args = parallel_executable.split(' ') args.extend(['ESMF_RegridWeightGen', '-s', mali_scripfile, '-d', slm_scripfile, - '-w', mapping_file_mali_to_slm, + '-w', 'mapping_file_mali_to_slm.nc', '-m', method_mali_to_slm, '-i', '-64bit_offset', '--netcdf4', '--src_regional']) check_call(args, logger) + args = parallel_executable.split(' ') args.extend(['ESMF_RegridWeightGen', '-s', slm_scripfile, '-d', mali_scripfile, - '-w', mapping_file_slm_to_mali, + '-w', 'mapping_file_slm_to_mali.nc', '-m', method_slm_to_mali, '-i', '-64bit_offset', '--netcdf4', '--dst_regional']) From 101e6ad714428c1baffeb4eb1401bab1d06e74b7 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Fri, 29 Sep 2023 12:02:27 -0600 Subject: [PATCH 05/23] set up SLM-related files and folders --- .../tests/slm/circ_icesheet/run_model.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/run_model.py b/compass/landice/tests/slm/circ_icesheet/run_model.py index d115d63784..77ae5900e7 100644 --- a/compass/landice/tests/slm/circ_icesheet/run_model.py +++ b/compass/landice/tests/slm/circ_icesheet/run_model.py @@ -1,3 +1,5 @@ +import os + from compass.model import run_model from compass.step import Step @@ -58,12 +60,23 @@ def __init__(self, test_case, name='run_model', self.add_input_file(filename='namelist.sealevel', package='compass.landice.tests.slm', copy=True) - + self.add_input_file(filename='smb_forcing.nc', + target='../setup_mesh/smb_forcing.nc') + self.add_input_file(filename='mapping_file_mali_to_slm.nc', + target='../setup_mesh/' + 'mapping_file_mali_to_slm.nc') + self.add_input_file(filename='mapping_file_slm_to_mali.nc', + target='../setup_mesh/' + 'mapping_file_slm_to_mali.nc') self.add_model_as_input() self.add_output_file(filename='output.nc') - # no setup() is needed + def setup(self): + os.makedirs(os.path.join(self.work_dir, 'OUTPUT_SLM/'), + exist_ok='True') + os.makedirs(os.path.join(self.work_dir, 'ICELOAD_SLM/'), + exist_ok='True') def run(self): """ From 5ef2e6491a97add598c788b3955d1d72645ab622 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Tue, 3 Oct 2023 13:37:13 -0600 Subject: [PATCH 06/23] Edit landice namelist and stream files --- .../tests/slm/circ_icesheet/namelist.landice | 55 ++++++++++++++++++- .../tests/slm/circ_icesheet/run_model.py | 4 +- .../tests/slm/circ_icesheet/streams.landice | 46 ++++++++++++++-- 3 files changed, 97 insertions(+), 8 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/namelist.landice b/compass/landice/tests/slm/circ_icesheet/namelist.landice index 13c97da547..e654cc6ee1 100644 --- a/compass/landice/tests/slm/circ_icesheet/namelist.landice +++ b/compass/landice/tests/slm/circ_icesheet/namelist.landice @@ -1,2 +1,53 @@ -config_dt = '0001-00-00_00:00:00' -config_run_duration = '0002-00-00_00:00:00' +&velocity_solver + config_velocity_solver = 'none' +/ +&advection + config_thickness_advection = 'fo' + config_tracer_advection = 'none' +/ +&solidearth + config_uplift_method = 'sealevelmodel' + config_slm_coupling_interval = '0001-00-00_00:00:00' + config_MALI_to_SLM_weights_file = 'mapping_file_mali_to_slm.nc' + config_SLM_to_MALI_weights_file = 'mapping_file_slm_to_mali.nc' +/ +&calving + config_calving = 'none' + config_restore_calving_front = .false. +/ +&thermal_solver + config_thermal_solver = 'none' +/ +&iceshelf_melt + config_basal_mass_bal_float = 'none' + config_front_mass_bal_grounded = 'none' +/ +&physical_parameters + config_ice_density = 910.0 + config_ocean_density = 1028.0 + config_sea_level = 0.0 + config_dynamic_thickness = 10.0 +/ +&time_integration + config_dt = '0001-00-00_00:00:00' + config_time_integration = 'forward_euler' + config_adaptive_timestep = .false. +/ +&time_management + config_do_restart = .false. + config_restart_timestamp_name = 'restart_timestamp' + config_start_time = '2015-01-01_00:00:00' + config_stop_time = '2055-01-01_00:00:00' + config_run_duration = 'none' + config_calendar_type = 'noleap' +/ +&io + config_pio_stride = 128 +/ +&AM_globalStats + config_AM_globalStats_enable = .true. + config_AM_globalStats_compute_interval = 'output_interval' + config_AM_globalStats_stream_name = 'globalStatsOutput' + config_AM_globalStats_compute_on_startup = .true. + config_AM_globalStats_write_on_startup = .true. +/ diff --git a/compass/landice/tests/slm/circ_icesheet/run_model.py b/compass/landice/tests/slm/circ_icesheet/run_model.py index 77ae5900e7..9b0b8bd117 100644 --- a/compass/landice/tests/slm/circ_icesheet/run_model.py +++ b/compass/landice/tests/slm/circ_icesheet/run_model.py @@ -64,10 +64,10 @@ def __init__(self, test_case, name='run_model', target='../setup_mesh/smb_forcing.nc') self.add_input_file(filename='mapping_file_mali_to_slm.nc', target='../setup_mesh/' - 'mapping_file_mali_to_slm.nc') + 'mapping_file_mali_to_slm.nc') self.add_input_file(filename='mapping_file_slm_to_mali.nc', target='../setup_mesh/' - 'mapping_file_slm_to_mali.nc') + 'mapping_file_slm_to_mali.nc') self.add_model_as_input() self.add_output_file(filename='output.nc') diff --git a/compass/landice/tests/slm/circ_icesheet/streams.landice b/compass/landice/tests/slm/circ_icesheet/streams.landice index 3f0648ada8..0a5dc1eaf0 100644 --- a/compass/landice/tests/slm/circ_icesheet/streams.landice +++ b/compass/landice/tests/slm/circ_icesheet/streams.landice @@ -1,29 +1,67 @@ + type="input" + filename_template="landice_grid.nc" + input_interval="initial_only"/> + clobber_mode="overwrite" + precision="double"> + + + + + + + + + + + + + + + + + + + + + From b049af7e2fa539a7e351615ef4896dbb24ae22b9 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Wed, 11 Oct 2023 16:46:18 -0600 Subject: [PATCH 07/23] temp/uncleaned commit: setup mesh and create smb forcing --- .../tests/slm/circ_icesheet/setup_mesh.py | 142 +++++++++--------- 1 file changed, 72 insertions(+), 70 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index 7e33521704..ae88229fe5 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -1,12 +1,14 @@ import os import shutil +import netCDF4 import numpy as np from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call from mpas_tools.mesh.conversion import convert, cull from mpas_tools.planar_hex import make_planar_hex_mesh from mpas_tools.scrip.from_mpas import scrip_from_mpas +from mpas_tools.translate import center, translate from netCDF4 import Dataset as NetCDFFile from compass.model import make_graph_file @@ -22,7 +24,7 @@ class SetupMesh(Step): ---------- """ - def __init__(self, test_case): + def __init__(self, test_case, name='setup_mesh'): """ Create the step @@ -31,7 +33,7 @@ def __init__(self, test_case): test_case : compass.TestCase The test case this step belongs to """ - super().__init__(test_case=test_case, name='setup_mesh') + super().__init__(test_case=test_case, name=name) self.add_output_file(filename='graph.info') self.add_output_file(filename='landice_grid.nc') @@ -49,12 +51,20 @@ def run(self): nx = section.getint('nx') ny = section.getint('ny') dc = section.getfloat('dc') - print('calling the mesh creation function') + # call the mesh creation function dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=True, nonperiodic_y=True) dsMesh = cull(dsMesh, logger=logger) dsMesh = convert(dsMesh, logger=logger) + # translating the mesh center to x=0 & y=0 + center(dsMesh) + # shift the center to a quarter or radius + r0 = section.getfloat('r0') + shift = r0 + print(f'shifting the center by {shift} meters') + translate(dsMesh, shift, shift) + # write to a file write_netcdf(dsMesh, 'mpas_grid.nc') levels = 3 @@ -93,7 +103,7 @@ def _setup_circsheet_initial_conditions(config, logger, filename): A logger for output from the step filename : str - file to setup circ_icesheet + file to setup circular icesheet """ section = config['circ_icesheet'] ice_type = section.get('ice_type') @@ -139,8 +149,7 @@ def _setup_circsheet_initial_conditions(config, logger, filename): x0 = 0.0 y0 = 0.0 r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 - print('HH==max r is', max(r)) - print('HH==min r is', min(r)) + # Assign variable values for the circular ice sheet # Set default value for non-circular cells thickness[:] = 0.0 @@ -148,8 +157,6 @@ def _setup_circsheet_initial_conditions(config, logger, filename): # (thickness will be NaN otherwise) thickness_field = thickness[0, :] - print('r0', r0) - print('r', r) if ice_type == 'cylinder': logger.info('cylinder ice type is chosen') thickness_field[r < r0] = h0 @@ -194,7 +201,7 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): mesh file created in the current step filename : str - file to setup circ_icesheet + file to setup circular icesheet """ section = config['circ_icesheet'] r0 = section.getfloat('r0') @@ -202,9 +209,14 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): section = config['smb_forcing'] direction = section.get('direction') + start_year = int(section.get('start_year')) + end_year = int(section.get('end_year')) + dt_year = int(section.get('dt_year')) + dHdt = section.getfloat('dHdt') + dRdt = section.getfloat('dRdt') # other constants used in the fuction - pi = 3.1415926535898 # pi value used in the SLM + # pi = 3.1415926535898 # pi value used in the SLM rhoi = 910.0 # ice density # Open the file, get needed dimensions @@ -214,73 +226,66 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): # Get variables xCell = smbfile.variables['xCell'] yCell = smbfile.variables['yCell'] - # nCells = smbfile.dimensions['nCells'] - sfcMassBal = smbfile.variables['sfcMassBal'] + nCells = smbfile.dimensions['nCells'].size - # set xtime variable - # smbfile.variables['xtime'] + # time interval for the forcing + t_array = np.arange(start_year, end_year, dt_year) + nt = len(t_array) # initialize SMB value everywhere in the mesh - sfcMassBal[:, :] = 0.0 + sfcMassBal = np.zeros((nt, nCells)) - # define the time (1 year interval) - dt = 1 - t_array = np.arange(0, 3, dt) - x0 = 0.0 - y0 = 0.0 + # Find center of domain + x0 = xCell[:].min() + 0.5 * (xCell[:].max() - xCell[:].min()) + y0 = yCell[:].min() + 0.5 * (yCell[:].max() - yCell[:].min()) + # calculate the radius of circ ice sheet r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 - dVdt = -9.678E13 - dHdt = 20.0 - dRdt = 40000.0 - if direction == 'vertical': - print('creating a file that prescribes vertical smb') - # dHdt = dVdt / (pi*(r0*1000.0)**2) # change in ice thickness in m - smb = dHdt * rhoi / (365.0 * 24 * 60 * 60) # smb in kg/m^2/s + logger.info('creating a file that prescribes vertical smb') + smb = dHdt * rhoi / (365.0 * 24 * 60 * 60) + logger.info(f'prescribed smb is {smb} kg/m2/s') indx = np.where(r <= r0)[0] - # for loop through time - for t in t_array: - print((t)) - Ht = h0 + (dHdt * t * dt) # calculate a new height each time - print(f'new height will be: {Ht}km') + for t in range(len(t_array)): + Ht = h0 + (dHdt * (t + 1)) # calculate a new height each time + logger.info(f'At time {start_year + dt_year}, \ + new height will be {Ht} km') sfcMassBal[t, indx] = smb # assign sfcMassBal - smbfile.variables['sfcMassBal'][t, :] = sfcMassBal[t, :] + smbfile.variables['sfcMassBal'][:, :] = sfcMassBal[:, :] elif direction == 'horizontal': - print('creating a file that prescribes horizontal smb') - # dRdt = (abs(Rf - r0)) / len(t_array) - smb = (10000.0) * rhoi / (365.0 * 24.0 * 60.0 * 60.0) - print('-smb is', -smb) + logger.info('creating a file that prescribes horizontal smb') + smb = -1 * ((10000.0) * rhoi / (365.0 * 24.0 * 60.0 * 60.0)) + print(f'prescribed smb is {smb} kg/m2/s') # for loop through time - for t in t_array: - print((t)) - # Rt = abs(np.sqrt(((dVdt * t * dt) / (pi * h0)) + (r0*1000)**2)) - # calculate a new radius at each time; for constant radius change - Rt = (r0 - (dRdt * t * dt)) - if (np.isnan(Rt) or Rt < 0): + for t in range(len(t_array)): + # calculate a new radius at each time + Rt = r0 + (dRdt * (t + 1)) + if (Rt < 0): Rt = 0.0 - smbfile.variables['sfcMassBal'][t, :] = 0.0 - print(f'new radius is: {Rt}km') + sfcMassBal[t, :] = 0.0 + print(f'At time {start_year+dt_year*t}, \ + new radius will be {Rt} km') else: - indx = np.where(r >= Rt / 1000)[0] - print('idx', indx) - smbfile.variables['sfcMassBal'][t, indx] = -smb - print(f'new radius is: {Rt / 1000}km') - smbfile.variables['sfcMassBal'][t, :] = sfcMassBal[t, :] - - elif direction == 'dome-halfar': - for t in t_array: - Rt = (((dVdt * dt * t) * 3 / (2 * pi) + - (r0 * 1000) ** 3)) ** (1.0 / 3) - print(f'new radius is: {Rt / 1000}km') - # smb = (dVdt * t * dt) / (2 * pi * (Rt ** 2)) \ - # * rhoi / (365.0 * 24.0 * 60.0 * 60.0)i - indx = np.where(r > r0)[0] - smbfile.variables['sfcMassBal'][t, :] = smb - print(f'smb is: {smb}') - print('ice thickness will change both horizontally and vertically') + indx = np.where(r >= Rt)[0] + sfcMassBal[t, indx] = smb + print(f'At time {start_year+dt_year*t}, \ + new radius will be {Rt/1000} km') + smbfile.variables['sfcMassBal'][:, :] = sfcMassBal[:, :] + + # add xtime variable + smbfile.createDimension('StrLen', 64) + xtime = smbfile.createVariable('xtime', 'S1', ('Time', 'StrLen')) + + # initialize SMB value everywhere in the mesh + xtime_str = [] + for t_index in range(len(t_array)): + yr = start_year + (t_index * dt_year) + xtime_str = f'{int(yr)}-01-01_00:00:00'.ljust(64) + xtime_char = netCDF4.stringtochar(np.array([xtime_str], 'S64'), + encoding='utf-8') + xtime[t_index, :] = xtime_char smbfile.close() @@ -330,18 +335,16 @@ def _build_mapping_files(config, logger, mali_mesh_file): check_call(args, logger=logger) - # create scrip files for source and destination grids - logger.info('creating scrip file for the mali mesh') - mali_mesh_copy = f'{mali_mesh_file}_copy' - shutil.copy(mali_mesh_file, mali_mesh_copy) - + # adjust the lat-lon values args = ['set_lat_lon_fields_in_planar_grid.py', - '--file', mali_mesh_copy, + '--file', mali_mesh_file, '--proj', 'ais-bedmap2-sphere'] check_call(args, logger=logger) - scrip_from_mpas(mali_mesh_copy, mali_scripfile) + # create scrip files for source and destination grids + logger.info('creating scrip file for the mali mesh') + scrip_from_mpas(mali_mesh_file, mali_scripfile) # create a mapping file using ESMF weight gen logger.info('Creating a mapping file...') @@ -374,4 +377,3 @@ def _build_mapping_files(config, logger, mali_mesh_file): logger.info('Removing the temporary mesh and scripfiles...') os.remove(slm_scripfile) os.remove(mali_scripfile) - os.remove(mali_mesh_copy) From cf7c2e1642fe6379bd074da81f9ecaceacf7f41f Mon Sep 17 00:00:00 2001 From: hollyhan Date: Wed, 11 Oct 2023 16:57:36 -0600 Subject: [PATCH 08/23] Add config options for the SLM --- .../tests/slm/circ_icesheet/circ_icesheet.cfg | 35 ++++++++++++++----- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg index 52717efdde..8b1191acab 100644 --- a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -1,19 +1,22 @@ -# config options for slm circ_icesheet test case +# Config options for slm circ_icesheet test case [circ_icesheet] # The number of cells in the x,y directions for the planar uniform mesh -# for 3000km mesh, need to have 150 for dc=20km -nx = 150 -ny = 150 +# ny must be an even number +nx = 601 +ny = 602 # resolution (in m) for the circular ice sheet # i.e. The distance in meters between adjacent cell centers. dc = 20000.0 -# the ice shape type ('cylinder', and for dome, 'dome-halfar' or 'dome-cism') +# list of resolutions in km delimited by ',' without space +resolutions = 40,20 + +# Ice shape type ('cylinder', and for dome, 'dome-halfar' or 'dome-cism') # r0 and h0: initial radius and height of the cylinder ice in meters ice_type = cylinder -r0 = 15000.0 +r0 = 2000000.0 h0 = 3000.0 # 'True' if manually want to set bedTopography elevation @@ -26,14 +29,27 @@ topo_elev = -1500.0 # center of the domain put_origin_on_a_cell = False -# config options for surface mass balance forcing +# Config options for surface mass balance forcing [smb_forcing] -# direction in which SMB is applied ('horizontal' or 'vertical') +# Start, end and interval years for SMB forcing +start_year = 2015 +end_year = 2116 +dt_year = 1 + +# Direction in which SMB is applied ('horizontal' or 'vertical') +# and amount of ice to melt direction = horizontal +# Change in radius in meters; used when 'direction == horizontal' +dRdt = -2000.0 +# Change in height in meters;sed when 'direction == vertical' +dHdt = -20.0 # config options for the sea-level model [slm] +# True if MALI-SLM are coupled +coupling = True + # Gauss-Legendre degree slm_res = 512 @@ -41,6 +57,9 @@ slm_res = 512 mapping_method_mali_to_slm = conserve mapping_method_slm_to_mali = bilinear +# ratio of MALI-SLM coupling interval and MALI output interval in integer years +time_stride = 1 + # config options related to visualization [circ_icesheet_viz] From 24bab033717962812b68641aa6ad00a97bb7ba16 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Wed, 11 Oct 2023 17:00:03 -0600 Subject: [PATCH 09/23] Set up steps through the list of MALI grid resolutions --- .../tests/slm/circ_icesheet/__init__.py | 25 +++++++++++++------ 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/__init__.py b/compass/landice/tests/slm/circ_icesheet/__init__.py index a40c0f2f94..d82a209238 100644 --- a/compass/landice/tests/slm/circ_icesheet/__init__.py +++ b/compass/landice/tests/slm/circ_icesheet/__init__.py @@ -23,22 +23,33 @@ def __init__(self, test_group): test_group : compass.landice.tests.dome.Dome The test group that this test case belongs to The resolution or type of mesh of the test case + + name : str, optional + the name of the test case """ name = 'circular_icesheet_test' subdir = name super().__init__(test_group=test_group, name=name, subdir=subdir) - self.add_step(SetupMesh(test_case=self)) - - self.add_step(RunModel(test_case=self, ntasks=4, openmp_threads=1, - name='run_step')) + def configure(self): + """ + Set up the desired mesh-resolution tests + Read the list of resolutions from the config + """ + config = self.config + section = config['circ_icesheet'] + resolutions = section.get('resolutions').split(',') + print('list of resolutione is ', resolutions) + + for res in resolutions: + self.add_step(SetupMesh(test_case=self, + name=f'{res}km_res/setup_mesh')) + self.add_step(RunModel(test_case=self, ntasks=1, openmp_threads=1, + name=f'{res}km_res/run_model')) step = Visualize(test_case=self) self.add_step(step, run_by_default=True) - # no configure() method is needed because we will use the default - # config options - # no run() method is needed because we're doing the default: running all # steps From e49e081e632374cd0e4dcebea557812573c42ecf Mon Sep 17 00:00:00 2001 From: hollyhan Date: Mon, 23 Oct 2023 10:04:13 -0600 Subject: [PATCH 10/23] Change the camel-case named variables to non-camel case for the ice mass change variables --- .../landice/tests/slm/circ_icesheet/circ_icesheet.cfg | 4 ++-- compass/landice/tests/slm/circ_icesheet/setup_mesh.py | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg index 8b1191acab..fbc30ffd64 100644 --- a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -41,9 +41,9 @@ dt_year = 1 # and amount of ice to melt direction = horizontal # Change in radius in meters; used when 'direction == horizontal' -dRdt = -2000.0 +drdt = -2000.0 # Change in height in meters;sed when 'direction == vertical' -dHdt = -20.0 +dhdt = -20.0 # config options for the sea-level model [slm] diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index ae88229fe5..3e3ced37d4 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -212,8 +212,8 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): start_year = int(section.get('start_year')) end_year = int(section.get('end_year')) dt_year = int(section.get('dt_year')) - dHdt = section.getfloat('dHdt') - dRdt = section.getfloat('dRdt') + dhdt = section.getfloat('dhdt') + drdt = section.getfloat('drdt') # other constants used in the fuction # pi = 3.1415926535898 # pi value used in the SLM @@ -243,12 +243,12 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): if direction == 'vertical': logger.info('creating a file that prescribes vertical smb') - smb = dHdt * rhoi / (365.0 * 24 * 60 * 60) + smb = dhdt * rhoi / (365.0 * 24 * 60 * 60) logger.info(f'prescribed smb is {smb} kg/m2/s') indx = np.where(r <= r0)[0] # for loop through time for t in range(len(t_array)): - Ht = h0 + (dHdt * (t + 1)) # calculate a new height each time + Ht = h0 + (dhdt * (t + 1)) # calculate a new height each time logger.info(f'At time {start_year + dt_year}, \ new height will be {Ht} km') sfcMassBal[t, indx] = smb # assign sfcMassBal @@ -261,7 +261,7 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): # for loop through time for t in range(len(t_array)): # calculate a new radius at each time - Rt = r0 + (dRdt * (t + 1)) + Rt = r0 + (drdt * (t + 1)) if (Rt < 0): Rt = 0.0 sfcMassBal[t, :] = 0.0 From 600495fb31802ff43469e9de9f6b5d703c452e5f Mon Sep 17 00:00:00 2001 From: hollyhan Date: Mon, 23 Oct 2023 11:18:32 -0600 Subject: [PATCH 11/23] Add post-processing and plotting functionality to the visualization test case --- .../tests/slm/circ_icesheet/visualize.py | 625 ++++++++++++++---- 1 file changed, 488 insertions(+), 137 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/visualize.py b/compass/landice/tests/slm/circ_icesheet/visualize.py index 30ea21d857..df51b3352c 100644 --- a/compass/landice/tests/slm/circ_icesheet/visualize.py +++ b/compass/landice/tests/slm/circ_icesheet/visualize.py @@ -1,6 +1,8 @@ +import os + import matplotlib.pyplot as plt -import netCDF4 -import numpy +import numpy as np +import xarray as xr from compass.step import Step @@ -34,167 +36,516 @@ def __init__(self, test_case, name='visualize', subdir=None, """ super().__init__(test_case=test_case, name=name, subdir=subdir) - self.add_input_file(filename='output.nc', - target='../{}/output.nc'.format(input_dir)) + def setup(self): + """ - # depending on settings, this may produce no outputs, so we won't add - # any + """ + config = self.config + section = config['circ_icesheet'] + resolutions = section.get('resolutions').split(',') - # no setup method is needed + for res in resolutions: + self.add_input_file(filename=f'output_res{res}km.nc', + target=f'../{res}km_res/run_model/' + 'output/output.nc') def run(self): """ Run this step of the test case """ - visualize_slm_circsheet(self.config, self.logger, filename='output.nc') - - -def visualize_slm_circsheet(config, logger, filename): + logger = self.logger + config = self.config + section = config['circ_icesheet'] + resolutions = section.get('resolutions').split(',') + + section = config['slm'] + coupling = section.getboolean('coupling') + + section = config['circ_icesheet_viz'] + save_images = section.getboolean('save_images') + hide_figs = section.getboolean('hide_figs') + + # visualize run model results + for res in resolutions: + run_path = f'../{res}km_res/run_model/' + logger.info(f'analyzing & visualizing outputs in path: {run_path}') + visualize_slm_circsheet(config, logger, res) + + # calculate and plot rmse in SLC + if coupling and len(resolutions) > 1: + ncells_list = list() + rmse_list = list() + for res in resolutions: + run_data = output_analysis(config, logger, + res, run_path) + ncells_res = run_data.ncells + rmse_res = run_data.rmse_slc + ncells_list.append(ncells_res) + rmse_list.append(rmse_res) + + ncells = np.array(ncells_list) + rmse = np.array(rmse_list) + + # plot rmse errors + p = np.polyfit(np.log10(ncells), np.log10(rmse), 1) + conv = abs(p[0]) * 2.0 + error_fit = ncells**p[0] * 10**p[1] + + plt.figure(1) + plt.loglog(ncells, error_fit, 'k') + plt.loglog(ncells, rmse, 'or') + plt.annotate('Order of Convergence = {}'.format(np.round(conv, 3)), + xycoords='axes fraction', xy=(0.3, 0.95), + fontsize=14) + plt.xlabel('Number of Grid Cells', fontsize=14) + plt.ylabel('L2 Norm', fontsize=14) + plt.title('MALI-SLM RMSE in SLC') + if save_images: + plt.savefig('RMSE.png', bbox_inches='tight') + plt.clf() + if hide_figs: + logger.info("Plot display disabled with" + "hide_plot config option.") + else: + plt.show() + + +def visualize_slm_circsheet(config, logger, res): """ - Plot the output from a dome test case + Plot the output from a circular ice sheet test case Parameters ---------- config : compass.config.CompassConfigParser - Configuration options for this test case, a combination of the defaults - for the machine, core and configuration + Configuration options for this test case, + a combination of the defaults for the machine, + core and configuration logger : logging.Logger A logger for output from the step - filename : str - file to visualize + res : str + resolution of MALI mesh """ - section = config['dome_viz'] - - time_slice = section.getint('time_slice') + section = config['circ_icesheet_viz'] save_images = section.getboolean('save_images') hide_figs = section.getboolean('hide_figs') - # Note: this may be slightly wrong for some calendar types! - secInYr = 3600.0 * 24.0 * 365.0 - - f = netCDF4.Dataset(filename, 'r') - - times = f.variables['xtime'] - thickness = f.variables['thickness'] - # dcEdge = f.variables['dcEdge'] - # bedTopography = f.variables['bedTopography'] # not needed - xCell = f.variables['xCell'] - yCell = f.variables['yCell'] - xEdge = f.variables['xEdge'] - yEdge = f.variables['yEdge'] - angleEdge = f.variables['angleEdge'] - temperature = f.variables['temperature'] - lowerSurface = f.variables['lowerSurface'] - upperSurface = f.variables['upperSurface'] - normalVelocity = f.variables['normalVelocity'] - # uReconstructX = f.variables['uReconstructX'] - uReconstructX = f.variables['uReconstructX'] - uReconstructY = f.variables['uReconstructY'] - - vert_levs = len(f.dimensions['nVertLevels']) - - time_length = times.shape[0] - - logger.info("vert_levs = {}; time_length = {}".format(vert_levs, - time_length)) - - var_slice = thickness[time_slice, :] - - fig = plt.figure(1, facecolor='w') - fig.add_subplot(111, aspect='equal') - # C = plt.contourf(xCell, yCell, var_slice ) - plt.scatter(xCell[:], yCell[:], 80, var_slice, marker='h', - edgecolors='none') - plt.colorbar() - plt.title('thickness at time {}'.format(time_slice)) - plt.draw() - if save_images: - logger.info("Saving figures to files.") - plt.savefig('dome_thickness.png') - - fig = plt.figure(2) - fig.add_subplot(121, aspect='equal') - plt.scatter(xCell[:], yCell[:], 80, lowerSurface[time_slice, :], - marker='h', edgecolors='none') - plt.colorbar() - plt.title('lower surface at time {}'.format(time_slice)) - plt.draw() - fig.add_subplot(122, aspect='equal') - plt.scatter(xCell[:], yCell[:], 80, upperSurface[time_slice, :], - marker='h', edgecolors='none') - plt.colorbar() - plt.title('upper surface at time {}'.format(time_slice)) - plt.draw() - if save_images: - plt.savefig('dome_surfaces.png') - - fig = plt.figure(3) - for templevel in range(0, vert_levs): - fig.add_subplot(3, 4, templevel + 1, aspect='equal') - var_slice = temperature[time_slice, :, templevel] - # C = plt.contourf(xCell, yCell, var_slice ) - plt.scatter(xCell[:], yCell[:], 40, var_slice, marker='h', - edgecolors='none') - plt.colorbar() - plt.title('temperature at level {} at time {}'.format(templevel, - time_slice)) - plt.draw() + section = config['slm'] + coupling = section.getboolean('coupling') + slm_res = int(section.get('slm_res')) + + # get an instance of output analysis class + run_path = f'../{res}km_res/run_model/' + run_data = output_analysis(config, logger, res, run_path) + yrs = run_data.yrs + + # figure 1 + fig1 = plt.figure(1, facecolor='w') + ax1 = fig1.add_subplot(1, 1, 1) + ax1.plot(yrs, run_data.grnd_vol_unscaled, + label='MALI unscaled', linestyle=':', color='k') + ax1.plot(yrs, run_data.grnd_vol, label='MALI scaled', + linestyle='-', color='k') + if coupling: + ax1.plot(yrs, run_data.ice_vol_slm, label='SLM', + linestyle='-.', color='r') + ax1.set_xlabel('Year') + ax1.legend(loc='best', ncol=1, prop={'size': 10}) + ax1.set_title('Grounded Ice Mass') + ax1.set_ylabel('Mass (Gt)') + ax1.grid(True) if save_images: - plt.savefig('dome_temperature.png') - - fig = plt.figure(4) - fig.add_subplot(121, aspect='equal') - plt.scatter(xEdge[:], yEdge[:], 80, - normalVelocity[time_slice, :, vert_levs - 1] * secInYr, - marker='h', edgecolors='none') - plt.colorbar() - normalVel = normalVelocity[time_slice, :, vert_levs - 1] - plt.quiver(xEdge[:], yEdge[:], - numpy.cos(angleEdge[:]) * normalVel * secInYr, - numpy.sin(angleEdge[:]) * normalVel * secInYr) - plt.title('normalVelocity of bottom layer at time {}'.format(time_slice)) - plt.draw() - fig.add_subplot(122, aspect='equal') - plt.scatter(xEdge[:], yEdge[:], 80, - normalVelocity[time_slice, :, 0] * secInYr, marker='h', - edgecolors='none') - plt.colorbar() - normalVel = normalVelocity[time_slice, :, 0] - plt.quiver(xEdge[:], yEdge[:], - numpy.cos(angleEdge[:]) * normalVel * secInYr, - numpy.sin(angleEdge[:]) * normalVel * secInYr) - plt.title('normalVelocity of top layer at time {}'.format(time_slice)) - plt.draw() + plt.savefig(f'ice_mass_mali{res}km_slm{slm_res}.png') + fig1.clf() + + # figure 2 + fig2 = plt.figure(2, facecolor='w') + ax2 = fig2.add_subplot(1, 1, 1) + ax2.plot(yrs, run_data.dgrnd_vol_unscaled, + label='MALI unscaled', linestyle=':', color='k') + ax2.plot(yrs, run_data.dgrnd_vol, label='MALI scaled', + linestyle='-', color='k') + if coupling: + ax2.plot(yrs, run_data.dice_vol_slm, label='SLM', + linestyle='-.', color='r') + ax2.set_xlabel('Year') + ax2.legend(loc='best', ncol=1, prop={'size': 10}) + ax2.set_title('Change in Grounded Ice Mass') + ax2.set_ylabel('Mass (Gt)') + ax2.grid(True) if save_images: - plt.savefig('dome_normalVelocity.png') - - fig = plt.figure(5, facecolor='w') - fig.add_subplot(121, aspect='equal') - plt.scatter(xCell[:], yCell[:], 80, - uReconstructX[time_slice, :, 0] * secInYr, marker='h', - edgecolors='none') - plt.colorbar() - plt.quiver(xCell[:], yCell[:], uReconstructX[time_slice, :, 0] * secInYr, - uReconstructY[time_slice, :, 0] * secInYr) - plt.title('uReconstructX of top layer at time {}'.format(time_slice)) - plt.draw() - fig.add_subplot(122, aspect='equal') - plt.scatter(xCell[:], yCell[:], 80, - uReconstructY[time_slice, :, 0] * secInYr, marker='h', - edgecolors='none') - plt.colorbar() - plt.quiver(xCell[:], yCell[:], uReconstructX[time_slice, :, 0] * secInYr, - uReconstructY[time_slice, :, 0] * secInYr) - plt.title('uReconstructY of top layer at time {}'.format(time_slice)) - plt.draw() + plt.savefig(f'ice_mass_change_mali{res}km_slm{slm_res}.png') + fig2.clf() + + # figure 3 + fig3 = plt.figure(3, facecolor='w') + ax3 = fig3.add_subplot(1, 1, 1) + ax3.plot(yrs, run_data.SLCaf, label='MALI SLC VAF', + linestyle=':', color='k') + ax3.plot(yrs, run_data.SLCcorr_Aocn, + label='MALI SLCcorr (grnded ice flooded)', + linestyle='-.', color='k') + ax3.plot(yrs, run_data.SLCcorr_AocnBeta, + label='MALI SLCcorr (grnded ice NOT flooded)', + linestyle='-', color='b') + ax3.plot(yrs, run_data.SLCcorr_z0_AocnBeta, + label='MALI SLCcorr+z0 (grnded ice NOT flooded)', + linestyle='-', color='y') + if coupling: + ax3.plot(yrs, run_data.SLC_slm_Aocn, + label='SLM (grnded ice flooded)', linestyle='-.', color='r') + ax3.plot(yrs, run_data.SLC_slm_AocnBeta, + label='SLM (grnded ice NOT flooded)', linestyle='-', + color='r') + ax3.set_xlabel('Year') + ax3.legend(loc='best', ncol=1, prop={'size': 10}) + ax3.set_title('Ice-sheet contribution to SLC') + ax3.set_ylabel('Sea-level change (m)') + ax3.grid(True) if save_images: - plt.savefig('dome_uReconstruct.png') + plt.savefig(f'ice_contribution_to_SLC_mali{res}km_slm{slm_res}.png') + fig3.clf() + + # figure 4 & 5 + if coupling: + fig4 = plt.figure(4, facecolor='w') + ax4 = fig4.add_subplot(1, 1, 1) + ax4.plot(yrs, run_data.dgrnd_vol - run_data.dice_vol_slm, + label='MALI minus SLM', linestyle='-', color='k') + ax4.set_xlabel('Year') + ax4.legend(loc='best', ncol=1, prop={'size': 10}) + ax4.set_title('Difference in ice mass change') + ax4.set_ylabel('Mass (Gt)') + ax4.grid(True) + if save_images: + plt.savefig(f'diff_ice_mass_change_mali{res}km_slm{slm_res}.png') + fig4.clf() + + fig5 = plt.figure(5, facecolor='w') + ax5 = fig5.add_subplot(1, 1, 1) + ax5.plot(yrs, run_data.diff_slc, label='MALI minus SLM', + linestyle='-', color='k') + ax5.plot(yrs, run_data.diff_slc_z0, label='MALI minus SLM (z0)', + linestyle='-', color='k') + ax5.set_xlabel('Year') + ax5.legend(loc='best', ncol=1, prop={'size': 10}) + ax5.set_title('Difference in sea-level change') + ax5.set_ylabel('sea-level change (m)') + ax5.grid(True) + if save_images: + plt.savefig(f'diff_sea_level_change_mali{res}km_slm{slm_res}.png') + fig5.clf() if hide_figs: logger.info("Plot display disabled with hide_plot config option.") else: plt.show() - f.close() + +class output_analysis: + """ + Analyze outputs + + Attributes + ---------- + mali_slc : float + ice-sheet contribution to sea-level change calculated + based on MALI outputs + + slm_slc : float + ice-sheet contribution to sea-level change calculated + based on SLM outputs + + rmse : float + root mean square error between mali_slc and slm_slc + """ + def __init__(self, config, logger, res, run_path): + """ + Calculate sea-level change from run outputs + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, + a combination of the defaults for the machine, + core and configuration + + logger : logging.Logger + A logger for output from the step + + res : str + Resolution of MALI mesh + + run_path : str + Path to runs where the output file exists + + coupling : str + Whether MALI is coupled to the SLM + """ + self.config = config + self.logger = logger + self.res = res + self.run_path = run_path + + section = config['slm'] + coupling = section.getboolean('coupling') + + if coupling: + section = config['circ_icesheet_viz'] + Aocn_const = section.getfloat('Aocn_const') + AocnBeta_const = section.getfloat('AocnBeta_const') + + # mali output file name + fname_mali = f'output_res{res}km.nc' + + # read in the MALI outputs + DS = xr.open_mfdataset(fname_mali, combine='nested', + concat_dim='Time', + decode_timedelta=False) + + # default constants + rhoi = 910.0 # ice density in kg/m^3 + rhoo = 1000.0 # ocean density used by the SLM (MALI uses 1028.0) + rhow = 1000.0 # fresh water density + Aocn_const = 4.5007E+14 # area of global ocean in m2 + AocnBeta_const = 4.5007E+14 # ocean area including marine-based area + + self.ncells = DS.dims['nCells'] + cellMask = DS['cellMask'] + bed = DS['bedTopography'] + # bed0 = bed.isel(Time=0).load() + thickness = DS['thickness'] + areaCell = DS['areaCell'][0, :].load() + latCell = DS['latCell'][0, :].load() + lonCell = DS['lonCell'][0, :].load() - np.pi # correct range: [-pi pi] + self.lat_deg = latCell * 180 / np.pi + + # calculate the area distortion (scale) factor for the + # polar stereographic projection + self.k = np.zeros(len(latCell),) + k0 = (1 - np.sin(-71 * np.pi / 180)) / 2 + # center of the latitude + lat0 = -90 * np.pi / 180 + # standard parallel in radians where distortion should be zero (k=1) + lon0 = 0 # center of the longitude + for i in range(len(latCell)): + # expression for 'k' from p. 157 + # eqn. 21-4 in https://pubs.usgs.gov/pp/1395/report.pdf. + # p. 142 provides a table showing 'k' for stereographic projection + self.k[i] = 2 * k0 / (1 + (np.sin(lat0) * np.sin(latCell[i])) + + (np.cos(lat0) * np.cos(latCell[i]) * + np.cos(lonCell[i] - lon0))) + + # default MALI time steps from the MALI outputs + self.yrs_mali = DS['daysSinceStart'].load() / 365.0 + 2015.0 + if coupling: + # path to the SLM output data + fpath_slm = os.path.join(run_path, 'OUTPUT_SLM/') + config = self.config + section = config['slm'] + time_stride = int(section.get('time_stride')) + # reset the time indices to match the # of SLM timesteps + nt = np.arange(0, DS.dims['Time'], time_stride, dtype='i') + self.yrs = self.yrs_mali[nt] + z0 = np.zeros((len(self.yrs), )) + Aocn = np.zeros((len(self.yrs), )) + AocnBeta = np.zeros((len(self.yrs), )) + # get ice mass on the SLM interpolated from the MALI mesh + fname = os.path.join(fpath_slm, 'ice_volume') + self.ice_vol_slm = slm_outputs(fname).data * rhoi / 1.0e12 + self.dice_vol_slm = slm_outputs(fname).change_total * \ + rhoi / 1.0e12 # in Gt + # get slc correction and ocean area from the SLM + fname = os.path.join(fpath_slm, 'gmsle_deltaSL_Ocean_fixed') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + z0 = slm_outputs(fname).change_total + self.SLC_slm_Aocn = slm_outputs(fname).data + + fname = os.path.join(fpath_slm, 'gmsle_deltaSL_OceanBeta_fixed') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + self.SLC_slm_AocnBeta = slm_outputs(fname).data + + fname = os.path.join(fpath_slm, 'ocean_area') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + Aocn = slm_outputs(fname).data + logger.info(f'area of the ocean is: {Aocn}') + fname = os.path.join(fpath_slm, 'oceanBeta_area') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + AocnBeta = slm_outputs(fname).data + logger.info(f'area of the ocean Beta is: {AocnBeta}') + else: + # if not coupled, use default values for MALI outputs + # assuming MALI output interval is 1 year + nt = np.arange(0, DS.dims['Time'], 1, dtype='i') + self.yrs = self.yrs_mali + z0 = np.zeros((len(self.yrs), )) + Aocn = np.zeros((len(self.yrs), )) + AocnBeta = np.zeros(len(self.yrs, )) + + z0[:] = 0.0 + logger.info("'gmsle_change' file doesn't exist. " + "Setting z0 to zeros") + Aocn[:] = Aocn_const + logger.info("'ocean_area' file doesn't exist. Using " + f"constant ocean area defined: {Aocn_const}m2") + AocnBeta[:] = AocnBeta_const + logger.info("'ocean_areaBeta' file doesn't exist. Using" + f"constant oceanBeta area defined: {AocnBeta_const}m2") + + # calculate ice-sheet contribution to sea-level change + # create empty arrays + self.grnd_vol = np.zeros((len(nt), )) + self.grnd_vol_unscaled = np.zeros((len(nt), )) + self.dgrnd_vol = np.zeros((len(nt), )) + self.dgrnd_vol_unscaled = np.zeros((len(nt), )) + self.vaf = np.zeros((len(nt), )) + self.vaf_z0 = np.zeros((len(nt), )) + self.pov = np.zeros((len(nt), )) + self.pov_z0 = np.zeros((len(nt), )) + self.den = np.zeros((len(nt), )) + self.SLEaf = np.zeros((len(nt), )) + self.SLCaf = np.zeros((len(nt), )) + self.SLCpov = np.zeros((len(nt), )) + self.SLCden = np.zeros((len(nt), )) + self.SLCcorr_Aocn = np.zeros((len(nt), )) + self.SLEaf_AocnBeta = np.zeros((len(nt), )) + self.SLEaf_AocnBeta = np.zeros((len(nt), )) + self.SLCaf_AocnBeta = np.zeros((len(nt), )) + self.SLCpov_AocnBeta = np.zeros((len(nt), )) + self.SLCden_AocnBeta = np.zeros((len(nt), )) + self.SLCcorr_AocnBeta = np.zeros((len(nt), )) + self.SLEaf_z0_AocnBeta = np.zeros((len(nt), )) + self.SLCaf_z0_AocnBeta = np.zeros((len(nt), )) + self.SLCpov_z0_AocnBeta = np.zeros((len(nt), )) + self.SLCcorr_z0_AocnBeta = np.zeros((len(nt), )) + + for t in np.arange(0, len(nt)): + # get appropriate time index for the MALI output data + if coupling: + indx = t * time_stride + else: + indx = t + + bedt = bed[indx, :].load() + cellMaskt = cellMask[indx, :].load() + thicknesst = thickness[indx, :].load() + areaCellt = areaCell.sum() + logger.info(f'area of the mali domain is {areaCellt}') + + # calculation of sea-level change following Goelzer et al. 2020 + self.grnd_vol_unscaled[t] = (areaCell * grounded(cellMaskt) * + thicknesst).sum() * rhoi / 1.0e12 + self.grnd_vol[t] = (areaCell * grounded(cellMaskt) * + thicknesst / (self.k**2)).sum() * rhoi / 1.0e12 + self.vaf[t] = (areaCell * grounded(cellMaskt) / (self.k**2) * + np.maximum(thicknesst + (rhoo / rhoi) * + np.minimum(np.zeros((self.ncells,)), + bedt), 0.0)).sum() # eqn. 1 + self.pov[t] = (areaCell / (self.k**2) * + np.maximum(0.0 * bedt, -1.0 * bedt)).sum() # eqn.8 + self.vaf_z0[t] = (areaCell * grounded(cellMaskt) / (self.k**2) * + np.maximum(thicknesst + (rhoo / rhoi) * + np.minimum(np.zeros((self.ncells,)), + bedt + z0[t]), 0.0)).sum() # eqn. 13 + self.pov_z0[t] = (areaCell / (self.k**2) * np.maximum(0.0 * bedt, + -1.0 * bedt + z0[t])).sum() # eqn. 14 + self.den[t] = (areaCell / (self.k**2) * grounded(cellMaskt) * + (rhoi / rhow - rhoi / rhoo)).sum() # eqn. 10 + + # SLC(m) using the ocean area where anywhere below sea level + # is considered as ocean (even where marine-based ice exists) + self.SLEaf[t] = (self.vaf[t] / Aocn[t]) * (rhoi / rhoo) # eqn. 2 + self.SLCaf[t] = -1.0 * (self.SLEaf[t] - self.SLEaf[0]) # eqn. 3 + self.SLCpov[t] = -1.0 * (self.pov[t] / Aocn[t] - self.pov[0] / + Aocn[t]) # eqn. 9 + self.SLCden[t] = -1.0 * (self.den[t] / Aocn[t] - self.den[0] / + Aocn[t]) # eqn. 11 + self.SLCcorr_Aocn[t] = self.SLCaf[t] + self.SLCpov[t] + \ + self.SLCden[t] + + # SLC(m) including the z0 term using the ocean+iced area + # i.e. water can't go where marine-based ice exists + self.SLEaf_AocnBeta[t] = (self.vaf[t] / AocnBeta[t]) * \ + (rhoi / rhoo) + self.SLCaf_AocnBeta[t] = -1.0 * (self.SLEaf_AocnBeta[t] - + self.SLEaf_AocnBeta[0]) + self.SLCpov_AocnBeta[t] = -1.0 * (self.pov[t] / AocnBeta[t] - + self.pov[0] / AocnBeta[t]) + self.SLCden_AocnBeta[t] = -1.0 * (self.den[t] / AocnBeta[t] - + self.den[0] / AocnBeta[t]) + self.SLCcorr_AocnBeta[t] = self.SLCaf_AocnBeta[t] + \ + self.SLCpov_AocnBeta[t] + \ + self.SLCden_AocnBeta[t] + + # same as above but adding the eustatic sea-level change term 'z0' + self.SLEaf_z0_AocnBeta[t] = (self.vaf_z0[t] / AocnBeta[t]) * \ + (rhoi / rhoo) + self.SLCaf_z0_AocnBeta[t] = -1.0 * (self.SLEaf_z0_AocnBeta[t] - + self.SLEaf_z0_AocnBeta[0]) + self.SLCpov_z0_AocnBeta[t] = -1.0 * (self.pov_z0[t] / AocnBeta[t] - + self.pov_z0[0] / AocnBeta[t]) + self.SLCcorr_z0_AocnBeta[t] = self.SLCaf_z0_AocnBeta[t] + \ + self.SLCpov_z0_AocnBeta[t] + \ + self.SLCden_AocnBeta[t] + + # get total ice mass change at each time step + self.dgrnd_vol[t] = (self.grnd_vol[t] - self.grnd_vol[0]) + self.dgrnd_vol_unscaled[t] = (self.grnd_vol_unscaled[t] - + self.grnd_vol_unscaled[0]) + + DS.close() + + # calculate RMSE between MALI and SLM calculation of SLC + if coupling: + self.diff_slc = self.SLCcorr_AocnBeta - self.SLC_slm_AocnBeta + self.diff_slc_z0 = self.SLCcorr_z0_AocnBeta - self.SLC_slm_AocnBeta + self.rmse_slc = np.sqrt((self.diff_slc**2).mean()) + else: + logger.info('MALI is not coupled to the SLM. No RMSE to compute.') + self.rmse_slc = np.nan + + +def grounded(cellMask): + # return ((cellMask&32)//32) & ~((cellMask&4)//4) + return ((cellMask & 32) // 32) * np.logical_not((cellMask & 4) // 4) + + +class slm_outputs: + """ + Read and calculate the SLM outputs + + Attributes + ---------- + data : float + SLM text-formatted outputs + + change_dt : float + Change in data value between each time step + + change_total : float + Change in data value at a time step + with respect to the intial time step + """ + def __init__(self, filename): + """ + Calculate the sea-level model output file + + Parameters + ---------- + filename : str + Filename of SLM output data + """ + self.fname = filename + self.data = np.loadtxt(self.fname) + self.change_dt = np.zeros(len(self.data),) + self.change_total = np.zeros(len(self.data),) + for i in range(len(self.change_dt) - 1): + self.change_dt[i + 1] = (float(self.data[i + 1]) - + float(self.data[i])) + self.change_total[i + 1] = (float(self.data[i + 1]) - + float(self.data[0])) From 9d5d69b574fa86b406189f9613fea14627bfd561 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Fri, 27 Oct 2023 13:46:38 -0600 Subject: [PATCH 12/23] Calculate number of grid cells given the size of the domain --- .../landice/tests/slm/circ_icesheet/circ_icesheet.cfg | 9 ++++----- compass/landice/tests/slm/circ_icesheet/setup_mesh.py | 11 +++++++++-- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg index fbc30ffd64..379f4ff114 100644 --- a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -1,13 +1,12 @@ # Config options for slm circ_icesheet test case [circ_icesheet] -# The number of cells in the x,y directions for the planar uniform mesh -# ny must be an even number -nx = 601 -ny = 602 +# The size of the domain in meters in the x and y directions for the planar mesh +lx = 6000000.0 +ly = 6000000.0 # resolution (in m) for the circular ice sheet -# i.e. The distance in meters between adjacent cell centers. +# # i.e. The distance in meters between adjacent cell centers. dc = 20000.0 # list of resolutions in km delimited by ',' without space diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index 3e3ced37d4..58ac1e1631 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -48,9 +48,16 @@ def run(self): config = self.config section = config['circ_icesheet'] - nx = section.getint('nx') - ny = section.getint('ny') + lx = section.getfloat('lx') + ly = section.getfloat('ly') dc = section.getfloat('dc') + + # calculate the number of grid cells in x and y direction + # for the uniform, hexagonal planar mesh + nx = max(2 * int(0.5 * lx / dc + 0.5), 4) + # factor of 2/sqrt(3) because of hexagonal mesh + ny = max(2 * int(0.5 * ly * (2. / np.sqrt(3)) / dc + 0.5), 4) + # call the mesh creation function dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=True, From 4e1b184fb593a164273667d795fba3601a507239 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Tue, 31 Oct 2023 15:54:41 -0600 Subject: [PATCH 13/23] Calculate k-factors in slice in python --- .../landice/tests/slm/circ_icesheet/visualize.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/visualize.py b/compass/landice/tests/slm/circ_icesheet/visualize.py index df51b3352c..871a95719d 100644 --- a/compass/landice/tests/slm/circ_icesheet/visualize.py +++ b/compass/landice/tests/slm/circ_icesheet/visualize.py @@ -329,13 +329,12 @@ def __init__(self, config, logger, res, run_path): lat0 = -90 * np.pi / 180 # standard parallel in radians where distortion should be zero (k=1) lon0 = 0 # center of the longitude - for i in range(len(latCell)): - # expression for 'k' from p. 157 - # eqn. 21-4 in https://pubs.usgs.gov/pp/1395/report.pdf. - # p. 142 provides a table showing 'k' for stereographic projection - self.k[i] = 2 * k0 / (1 + (np.sin(lat0) * np.sin(latCell[i])) + - (np.cos(lat0) * np.cos(latCell[i]) * - np.cos(lonCell[i] - lon0))) + # expression for 'k' from p. 157 + # eqn. 21-4 in https://pubs.usgs.gov/pp/1395/report.pdf. + # p. 142 provides a table showing 'k' for stereographic projection + self.k[:] = 2 * k0 / (1 + (np.sin(lat0) * np.sin(latCell[:])) + + (np.cos(lat0) * np.cos(latCell[:]) * + np.cos(lonCell[:] - lon0))) # default MALI time steps from the MALI outputs self.yrs_mali = DS['daysSinceStart'].load() / 365.0 + 2015.0 From 3fea3f433c0e326923f4681886fbcbeae03fcd94 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Tue, 31 Oct 2023 19:49:58 -0600 Subject: [PATCH 14/23] Enable steps and analysis through SLM nglv values --- .../tests/slm/circ_icesheet/__init__.py | 22 +++-- .../tests/slm/circ_icesheet/circ_icesheet.cfg | 12 ++- .../tests/slm/circ_icesheet/setup_mesh.py | 7 +- .../tests/slm/circ_icesheet/visualize.py | 82 +++++++++++-------- 4 files changed, 75 insertions(+), 48 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/__init__.py b/compass/landice/tests/slm/circ_icesheet/__init__.py index d82a209238..70e1083324 100644 --- a/compass/landice/tests/slm/circ_icesheet/__init__.py +++ b/compass/landice/tests/slm/circ_icesheet/__init__.py @@ -40,14 +40,20 @@ def configure(self): """ config = self.config section = config['circ_icesheet'] - resolutions = section.get('resolutions').split(',') - print('list of resolutione is ', resolutions) - - for res in resolutions: - self.add_step(SetupMesh(test_case=self, - name=f'{res}km_res/setup_mesh')) - self.add_step(RunModel(test_case=self, ntasks=1, openmp_threads=1, - name=f'{res}km_res/run_model')) + mali_res = section.get('mali_res').split(',') + + section = config['slm'] + slm_nglv = section.get('slm_nglv').split(',') + print(f'list of MALI-mesh resolution is {mali_res} km.') + print(f'list of SLM Gauss-Legendre latitudinal points is {slm_nglv}.') + + for res in mali_res: + for nglv in slm_nglv: + self.add_step(SetupMesh(test_case=self, + name=f'mali{res}km_slm{nglv}/setup_mesh')) + self.add_step(RunModel(test_case=self, ntasks=1, + openmp_threads=1, + name=f'mali{res}km_slm{nglv}/run_model')) step = Visualize(test_case=self) self.add_step(step, run_by_default=True) diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg index 379f4ff114..6ea7ef8dc6 100644 --- a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -9,8 +9,9 @@ ly = 6000000.0 # # i.e. The distance in meters between adjacent cell centers. dc = 20000.0 -# list of resolutions in km delimited by ',' without space -resolutions = 40,20 +# list of MALI mesh resolutions in meters delimited by ',' without space +# i.e. The distance between adjacent cell centers. +mali_res = 20 # Ice shape type ('cylinder', and for dome, 'dome-halfar' or 'dome-cism') # r0 and h0: initial radius and height of the cylinder ice in meters @@ -49,7 +50,12 @@ dhdt = -20.0 # True if MALI-SLM are coupled coupling = True -# Gauss-Legendre degree +# List of the number of Gauss-Legendre points in latitude +# delimited by ',' without space +slm_nglv = 512 + +# Max spherical harmonics degree and order +# i.e. SLM resolution of the SLM slm_res = 512 # mapping method between the MALI and SLM grids diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index 58ac1e1631..d07d2b616f 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -323,6 +323,7 @@ def _build_mapping_files(config, logger, mali_mesh_file): section = config['slm'] slm_res = section.get('slm_res') + slm_nglv = section.get('slm_nglv') method_mali_to_slm = section.get('mapping_method_mali_to_slm') method_slm_to_mali = section.get('mapping_method_slm_to_mali') @@ -330,15 +331,15 @@ def _build_mapping_files(config, logger, mali_mesh_file): dc = section.getfloat('dc') mali_scripfile = f'mali{int(dc/1000)}km_scripfile.nc' - slm_scripfile = f'slm{slm_res}_scripfile.nc' + slm_scripfile = f'slm{slm_res}_nglv{slm_nglv}scripfile.nc' # create slm scripfile logger.info(f'creating scripfile for the SH-degree {slm_res} ' - f'resolution SLM grid') + f'resolution and {slm_nglv} nglv points for the SLM grid') args = ['ncremap', '-g', slm_scripfile, '-G', - f'latlon={slm_res},{2*int(slm_res)}#lat_typ=gss#lat_drc=n2s'] + f'latlon={slm_nglv},{2*int(slm_nglv)}#lat_typ=gss#lat_drc=n2s'] check_call(args, logger=logger) diff --git a/compass/landice/tests/slm/circ_icesheet/visualize.py b/compass/landice/tests/slm/circ_icesheet/visualize.py index 871a95719d..39d1a7aecb 100644 --- a/compass/landice/tests/slm/circ_icesheet/visualize.py +++ b/compass/landice/tests/slm/circ_icesheet/visualize.py @@ -42,12 +42,17 @@ def setup(self): """ config = self.config section = config['circ_icesheet'] - resolutions = section.get('resolutions').split(',') + mali_res = section.get('mali_res').split(',') - for res in resolutions: - self.add_input_file(filename=f'output_res{res}km.nc', - target=f'../{res}km_res/run_model/' - 'output/output.nc') + section = config['slm'] + slm_nglv = section.get('slm_nglv').split(',') + + for res in mali_res: + for nglv in slm_nglv: + self.add_input_file(filename=f'output_mali{res}km_' + 'slm{nglv}.nc', + target=f'../mali{res}km_slm{nglv}/' + 'run_model/output/output.nc') def run(self): """ @@ -56,35 +61,40 @@ def run(self): logger = self.logger config = self.config section = config['circ_icesheet'] - resolutions = section.get('resolutions').split(',') + mali_res = section.get('mali_res').split(',') section = config['slm'] coupling = section.getboolean('coupling') + slm_nglv = section.get('slm_nglv').split(',') section = config['circ_icesheet_viz'] save_images = section.getboolean('save_images') hide_figs = section.getboolean('hide_figs') # visualize run model results - for res in resolutions: - run_path = f'../{res}km_res/run_model/' - logger.info(f'analyzing & visualizing outputs in path: {run_path}') - visualize_slm_circsheet(config, logger, res) + num_files = 0 + for res in mali_res: + for nglv in slm_nglv: + run_path = f'../mali{res}km_slm{nglv}/run_model/' + logger.info(f'analyzing outputs in path: {run_path}') + visualize_slm_circsheet(config, logger, res, nglv) + num_files = num_files + 1 # calculate and plot rmse in SLC - if coupling and len(resolutions) > 1: + if coupling and num_files > 1: ncells_list = list() rmse_list = list() - for res in resolutions: - run_data = output_analysis(config, logger, - res, run_path) - ncells_res = run_data.ncells - rmse_res = run_data.rmse_slc - ncells_list.append(ncells_res) - rmse_list.append(rmse_res) - - ncells = np.array(ncells_list) - rmse = np.array(rmse_list) + for res in mali_res: + for nglv in slm_nglv: + run_data = output_analysis(config, logger, + res, nglv, run_path) + ncells_res = run_data.ncells + rmse_res = run_data.rmse_slc + ncells_list.append(ncells_res) + rmse_list.append(rmse_res) + + ncells = np.array(ncells_list) + rmse = np.array(rmse_list) # plot rmse errors p = np.polyfit(np.log10(ncells), np.log10(rmse), 1) @@ -110,7 +120,7 @@ def run(self): plt.show() -def visualize_slm_circsheet(config, logger, res): +def visualize_slm_circsheet(config, logger, res, nglv): """ Plot the output from a circular ice sheet test case @@ -126,6 +136,9 @@ def visualize_slm_circsheet(config, logger, res): res : str resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre points in latitude in the SLM grid """ section = config['circ_icesheet_viz'] save_images = section.getboolean('save_images') @@ -133,10 +146,10 @@ def visualize_slm_circsheet(config, logger, res): section = config['slm'] coupling = section.getboolean('coupling') - slm_res = int(section.get('slm_res')) + slm_nglv = int(section.get('slm_nglv')) # get an instance of output analysis class - run_path = f'../{res}km_res/run_model/' + run_path = f'../mali{res}km_slm{slm_nglv}/run_model/' run_data = output_analysis(config, logger, res, run_path) yrs = run_data.yrs @@ -156,7 +169,7 @@ def visualize_slm_circsheet(config, logger, res): ax1.set_ylabel('Mass (Gt)') ax1.grid(True) if save_images: - plt.savefig(f'ice_mass_mali{res}km_slm{slm_res}.png') + plt.savefig(f'ice_mass_mali{res}km_slm{slm_nglv}.png') fig1.clf() # figure 2 @@ -175,7 +188,7 @@ def visualize_slm_circsheet(config, logger, res): ax2.set_ylabel('Mass (Gt)') ax2.grid(True) if save_images: - plt.savefig(f'ice_mass_change_mali{res}km_slm{slm_res}.png') + plt.savefig(f'ice_mass_change_mali{res}km_slm{slm_nglv}.png') fig2.clf() # figure 3 @@ -204,7 +217,7 @@ def visualize_slm_circsheet(config, logger, res): ax3.set_ylabel('Sea-level change (m)') ax3.grid(True) if save_images: - plt.savefig(f'ice_contribution_to_SLC_mali{res}km_slm{slm_res}.png') + plt.savefig(f'ice_contribution_to_SLC_mali{res}km_slm{slm_nglv}.png') fig3.clf() # figure 4 & 5 @@ -219,7 +232,7 @@ def visualize_slm_circsheet(config, logger, res): ax4.set_ylabel('Mass (Gt)') ax4.grid(True) if save_images: - plt.savefig(f'diff_ice_mass_change_mali{res}km_slm{slm_res}.png') + plt.savefig(f'diff_ice_mass_change_mali{res}km_slm{slm_nglv}.png') fig4.clf() fig5 = plt.figure(5, facecolor='w') @@ -234,7 +247,7 @@ def visualize_slm_circsheet(config, logger, res): ax5.set_ylabel('sea-level change (m)') ax5.grid(True) if save_images: - plt.savefig(f'diff_sea_level_change_mali{res}km_slm{slm_res}.png') + plt.savefig(f'diff_sea_level_change_mali{res}km_slm{slm_nglv}.png') fig5.clf() if hide_figs: @@ -260,7 +273,7 @@ class output_analysis: rmse : float root mean square error between mali_slc and slm_slc """ - def __init__(self, config, logger, res, run_path): + def __init__(self, config, logger, res, nglv, run_path): """ Calculate sea-level change from run outputs @@ -277,15 +290,16 @@ def __init__(self, config, logger, res, run_path): res : str Resolution of MALI mesh + nglv : str + Number of Gauss-Legendre points in latitude in the SLM + run_path : str Path to runs where the output file exists - - coupling : str - Whether MALI is coupled to the SLM """ self.config = config self.logger = logger self.res = res + self.nglv = nglv self.run_path = run_path section = config['slm'] @@ -297,7 +311,7 @@ def __init__(self, config, logger, res, run_path): AocnBeta_const = section.getfloat('AocnBeta_const') # mali output file name - fname_mali = f'output_res{res}km.nc' + fname_mali = f'output_mali{res}km_slm{nglv}.nc' # read in the MALI outputs DS = xr.open_mfdataset(fname_mali, combine='nested', From 11b6498d222365f0fb936043be284d282e55b779 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Tue, 31 Oct 2023 20:05:00 -0600 Subject: [PATCH 15/23] Hardcode the mesh-center shift to 200 km --- compass/landice/tests/slm/circ_icesheet/setup_mesh.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index d07d2b616f..1c06273cb1 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -67,8 +67,7 @@ def run(self): # translating the mesh center to x=0 & y=0 center(dsMesh) # shift the center to a quarter or radius - r0 = section.getfloat('r0') - shift = r0 + shift = 200000.0 print(f'shifting the center by {shift} meters') translate(dsMesh, shift, shift) # write to a file From 5d013afd312e933d4bbc469ba3b952da5aca0580 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Wed, 8 Nov 2023 13:44:42 -0700 Subject: [PATCH 16/23] Change format of the netcdf MPAS grid file to cdf5 --- .../tests/slm/circ_icesheet/setup_mesh.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index 1c06273cb1..2e44745cff 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -1,11 +1,12 @@ import os import shutil +import mpas_tools.io import netCDF4 import numpy as np from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from mpas_tools.mesh.conversion import convert, cull +from mpas_tools.mesh.conversion import cull from mpas_tools.planar_hex import make_planar_hex_mesh from mpas_tools.scrip.from_mpas import scrip_from_mpas from mpas_tools.translate import center, translate @@ -58,20 +59,27 @@ def run(self): # factor of 2/sqrt(3) because of hexagonal mesh ny = max(2 * int(0.5 * ly * (2. / np.sqrt(3)) / dc + 0.5), 4) + mpas_tools.io.default_format = 'NETCDF4' + mpas_tools.io.default_engine = 'netcdf4' + # call the mesh creation function dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, nonperiodic_x=True, nonperiodic_y=True) dsMesh = cull(dsMesh, logger=logger) - dsMesh = convert(dsMesh, logger=logger) + # adding the time dimension is needed for netcdf4 formatting to work + dsMesh['xtime'] = ('Time', ['2015-01-01_00:00:00'.ljust(64)]) # translating the mesh center to x=0 & y=0 center(dsMesh) # shift the center to a quarter or radius shift = 200000.0 print(f'shifting the center by {shift} meters') translate(dsMesh, shift, shift) - # write to a file - write_netcdf(dsMesh, 'mpas_grid.nc') + + fname_culled = 'culled_mesh_before_cdf5.nc' + write_netcdf(dsMesh, fname_culled) + args = ['ncks', '-O', '-5', fname_culled, 'mpas_grid.nc'] + check_call(args, logger=logger) levels = 3 args = ['create_landice_grid_from_generic_MPAS_grid.py', @@ -94,6 +102,8 @@ def run(self): _build_mapping_files(config, logger, mali_mesh_file='landice_grid.nc') + os.remove(fname_culled) + def _setup_circsheet_initial_conditions(config, logger, filename): """ From c869044c680b4415718b1864c93b9dcd6d8d47b1 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Wed, 8 Nov 2023 19:37:00 -0700 Subject: [PATCH 17/23] Pass down mali resolution and SLM nglv as attributes --- .../tests/slm/circ_icesheet/__init__.py | 3 +- .../tests/slm/circ_icesheet/circ_icesheet.cfg | 6 +-- .../tests/slm/circ_icesheet/run_model.py | 23 +++++++--- .../tests/slm/circ_icesheet/setup_mesh.py | 45 ++++++++++++------- 4 files changed, 48 insertions(+), 29 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/__init__.py b/compass/landice/tests/slm/circ_icesheet/__init__.py index 70e1083324..ab3d8bbb60 100644 --- a/compass/landice/tests/slm/circ_icesheet/__init__.py +++ b/compass/landice/tests/slm/circ_icesheet/__init__.py @@ -50,7 +50,8 @@ def configure(self): for res in mali_res: for nglv in slm_nglv: self.add_step(SetupMesh(test_case=self, - name=f'mali{res}km_slm{nglv}/setup_mesh')) + name=f'mali{res}km_slm{nglv}/' + f'setup_mesh', res=res, nglv=nglv)) self.add_step(RunModel(test_case=self, ntasks=1, openmp_threads=1, name=f'mali{res}km_slm{nglv}/run_model')) diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg index 6ea7ef8dc6..b2b7ebd134 100644 --- a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -5,11 +5,7 @@ lx = 6000000.0 ly = 6000000.0 -# resolution (in m) for the circular ice sheet -# # i.e. The distance in meters between adjacent cell centers. -dc = 20000.0 - -# list of MALI mesh resolutions in meters delimited by ',' without space +# list of MALI mesh resolutions in kilometers delimited by ',' without space # i.e. The distance between adjacent cell centers. mali_res = 20 diff --git a/compass/landice/tests/slm/circ_icesheet/run_model.py b/compass/landice/tests/slm/circ_icesheet/run_model.py index 9b0b8bd117..db3a0b76d4 100644 --- a/compass/landice/tests/slm/circ_icesheet/run_model.py +++ b/compass/landice/tests/slm/circ_icesheet/run_model.py @@ -11,8 +11,8 @@ class RunModel(Step): Attributes ---------- """ - def __init__(self, test_case, name='run_model', - subdir=None, ntasks=1, min_tasks=None, openmp_threads=1): + def __init__(self, test_case, res, nglv, ntasks, name='run_model', + subdir=None, min_tasks=None, openmp_threads=1): """ Create a new test case @@ -21,17 +21,23 @@ def __init__(self, test_case, name='run_model', test_case : compass.TestCase The test case this step belongs to - name : str, optional - the name of the test case + res : str + Resolution of the MALI domain - subdir : str, optional - the subdirectory for the step. The default is ``name`` + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM grid - ntasks : int, optional + ntasks : int the number of tasks the step would ideally use. If fewer tasks are available on the system, the step will run on all available tasks as long as this is not below ``min_tasks`` + name : str, optional + the name of the test case + + subdir : str, optional + the subdirectory for the step. The default is ``name`` + min_tasks : int, optional the number of tasks the step requires. If the system has fewer than this number of tasks, the step will fail @@ -39,6 +45,9 @@ def __init__(self, test_case, name='run_model', openmp_threads : int, optional the number of OpenMP threads the step will use """ + self.res = res + self.nglv = nglv + if min_tasks is None: min_tasks = ntasks super().__init__(test_case=test_case, name=name, subdir=subdir, diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py index 2e44745cff..f6f1f16702 100644 --- a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -23,9 +23,14 @@ class SetupMesh(Step): Attributes ---------- + res : str + Resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM """ - def __init__(self, test_case, name='setup_mesh'): + def __init__(self, test_case, name, res, nglv): """ Create the step @@ -33,12 +38,20 @@ def __init__(self, test_case, name='setup_mesh'): ---------- test_case : compass.TestCase The test case this step belongs to + + res : str + Resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM """ super().__init__(test_case=test_case, name=name) self.add_output_file(filename='graph.info') self.add_output_file(filename='landice_grid.nc') + self.res = res + self.nglv = nglv # no setup() method is needed def run(self): @@ -51,7 +64,7 @@ def run(self): lx = section.getfloat('lx') ly = section.getfloat('ly') - dc = section.getfloat('dc') + dc = float(self.res) * 1000.0 # calculate the number of grid cells in x and y direction # for the uniform, hexagonal planar mesh @@ -99,7 +112,7 @@ def run(self): mali_mesh_file='landice_grid.nc', filename='smb_forcing.nc') - _build_mapping_files(config, logger, + _build_mapping_files(config, logger, self.res, self.nglv, mali_mesh_file='landice_grid.nc') os.remove(fname_culled) @@ -306,7 +319,7 @@ def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): smbfile.close() -def _build_mapping_files(config, logger, mali_mesh_file): +def _build_mapping_files(config, logger, res, nglv, mali_mesh_file): """ Build a mapping file if it does not exist. Mapping file is then used to remap the ismip6 source file in polarstero @@ -317,29 +330,29 @@ def _build_mapping_files(config, logger, mali_mesh_file): config : compass.config.CompassConfigParser Configuration options for a ismip6 forcing test case - cores : int - the number of cores for the ESMF_RegridWeightGen - logger : logging.Logger A logger for output from the step - mali_mesh_file : str, optional - The MALI mesh file if mapping file does not exist + name : str + Name of the step + + res : str + Resolution of MALI mesh - method_remap : str, optional - Remapping method used in building a mapping file + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM + + mali_mesh_file : str + The MALI mesh file """ section = config['slm'] slm_res = section.get('slm_res') - slm_nglv = section.get('slm_nglv') + slm_nglv = int(nglv) method_mali_to_slm = section.get('mapping_method_mali_to_slm') method_slm_to_mali = section.get('mapping_method_slm_to_mali') - section = config['circ_icesheet'] - dc = section.getfloat('dc') - - mali_scripfile = f'mali{int(dc/1000)}km_scripfile.nc' + mali_scripfile = f'mali{int(res)}km_scripfile.nc' slm_scripfile = f'slm{slm_res}_nglv{slm_nglv}scripfile.nc' # create slm scripfile From 82ea070b27e0d1e83bb1b11b2046a9c65057bf02 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Fri, 10 Nov 2023 09:07:53 -0700 Subject: [PATCH 18/23] Calculate and plot the heatmap of MPE in SLC MPE: mean percentage error SLC: sea-level change Also enable the SLM output class to take in an assigned number of time indices of outputs to be analyzied --- .../tests/slm/circ_icesheet/visualize.py | 222 ++++++++++++------ 1 file changed, 154 insertions(+), 68 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/visualize.py b/compass/landice/tests/slm/circ_icesheet/visualize.py index 39d1a7aecb..f9f63f70e0 100644 --- a/compass/landice/tests/slm/circ_icesheet/visualize.py +++ b/compass/landice/tests/slm/circ_icesheet/visualize.py @@ -50,7 +50,7 @@ def setup(self): for res in mali_res: for nglv in slm_nglv: self.add_input_file(filename=f'output_mali{res}km_' - 'slm{nglv}.nc', + f'slm{nglv}.nc', target=f'../mali{res}km_slm{nglv}/' 'run_model/output/output.nc') @@ -69,7 +69,6 @@ def run(self): section = config['circ_icesheet_viz'] save_images = section.getboolean('save_images') - hide_figs = section.getboolean('hide_figs') # visualize run model results num_files = 0 @@ -82,42 +81,112 @@ def run(self): # calculate and plot rmse in SLC if coupling and num_files > 1: - ncells_list = list() - rmse_list = list() + # ncells_list = list() + # rmse_list = list() + pct_err_slc_list = list() + pct_err_slc_z0_list = list() + pct_err_slc_z0_matrix = np.zeros((len(mali_res), len(slm_nglv))) + n = 0 for res in mali_res: for nglv in slm_nglv: + run_path = f'../mali{res}km_slm{nglv}/run_model/' run_data = output_analysis(config, logger, res, nglv, run_path) - ncells_res = run_data.ncells - rmse_res = run_data.rmse_slc - ncells_list.append(ncells_res) - rmse_list.append(rmse_res) - - ncells = np.array(ncells_list) - rmse = np.array(rmse_list) - - # plot rmse errors - p = np.polyfit(np.log10(ncells), np.log10(rmse), 1) - conv = abs(p[0]) * 2.0 - error_fit = ncells**p[0] * 10**p[1] - - plt.figure(1) - plt.loglog(ncells, error_fit, 'k') - plt.loglog(ncells, rmse, 'or') - plt.annotate('Order of Convergence = {}'.format(np.round(conv, 3)), - xycoords='axes fraction', xy=(0.3, 0.95), - fontsize=14) - plt.xlabel('Number of Grid Cells', fontsize=14) - plt.ylabel('L2 Norm', fontsize=14) - plt.title('MALI-SLM RMSE in SLC') - if save_images: - plt.savefig('RMSE.png', bbox_inches='tight') - plt.clf() - if hide_figs: - logger.info("Plot display disabled with" - "hide_plot config option.") - else: - plt.show() + # ncells_res = run_data.ncells + # rmse_res = run_data.rmse_slc + pct_err_slc_res = run_data.pct_err_slc + pct_err_slc_z0_res = run_data.pct_err_slm_z0 + # ncells_list.append(ncells_res) + # rmse_list.append(rmse_res) + pct_err_slc_list.append(pct_err_slc_res) + pct_err_slc_z0_list.append(pct_err_slc_z0_res) + + # ncells = np.array(ncells_list) + # rmse = np.array(rmse_list) + pct_err_slc = np.array(pct_err_slc_list) + pct_err_slc_z0 = np.array(pct_err_slc_z0_list) + + # assign the percentage error array in matrix + pct_err_slc_z0_matrix[n, :] = pct_err_slc_z0 + n += 1 + + # plot mean percentage error for fixed MALI res + plot_pct_err_slc(res, slm_nglv, pct_err_slc, + 'Mean Percent Error in SLC') + plot_pct_err_slc(res, slm_nglv, pct_err_slc_z0, + 'Mean Percent Error in SLC-z0') + # reset the lists for percentage error + del pct_err_slc_list[:], pct_err_slc_z0_list[:] + + # plot the heat map of mean percentage errors + plot_heatmap(pct_err_slc_z0_matrix, mali_res, slm_nglv, + save_images) + + +def plot_heatmap(data, mali_res, slm_nglv, save_images): + """ + Plot a heat map of mean percentage error in SLC. + + Parameters + ---------- + data : float + Matrix of arrays containing mean percentage error values + + mali_res : str + List of MALI-resolution strings + + slm_nglv : str + List of SLM nglv values + + save_images : str + Boolean on whether to save the plot or not + """ + + fig1 = plt.figure(1, facecolot='w') + ax1 = fig1.add_subplot(1, 1, 1) + im = ax1.imshow(data, cmap='Blues') + ax1.figure.colorbar(im, ax=ax1) + ax1.set_title('Mean Percentage Error (%)') + ax1.set_xticks(np.arange(len(slm_nglv)), labels=slm_nglv) + ax1.set_yticks(np.arange(len(mali_res)), labels=mali_res) + ax1.set_xlabel('# of GL nodes in Latitude') + ax1.set_ylabel('MALI domain resolution (km)') + if save_images: + plt.savefig('Heatmap_MPE.png') + fig1.clf() + + +def plot_pct_err_slc(mali_res, slm_nglv, data, figure_title): + """ + Plot mean percentage errors in sea-level change where + MALI-calculation is taken as a true value. + + Parameters + ---------- + mali_res : str + Resolution of the MALI domain + + slm_nglv : int + Array of the number of Gauss-Legendre nodes in latitude + in the SLM grid + + data : float + Data array containing the percentage error values + + figure_title : str + Title of the figure + """ + + fig1 = plt.figure(1, facecolor='w') + ax1 = fig1.add_subplot(1, 1, 1) + ax1.plot(np.array(slm_nglv), data, marker='o') + ax1.set_xlabel('# of Gauss-Legendre nodes in latitude (SLM)', fontsize=14) + ax1.legend(loc='best', ncol=1, prop={'size': 10}) + ax1.set_title(f'{figure_title} for {mali_res}km MALI res') + ax1.set_ylabel('percent error (%)', fontsize=14) + plt.savefig(f'{figure_title} for {mali_res}km MALI res.png', + bbox_inches='tight') + fig1.clf() def visualize_slm_circsheet(config, logger, res, nglv): @@ -146,11 +215,10 @@ def visualize_slm_circsheet(config, logger, res, nglv): section = config['slm'] coupling = section.getboolean('coupling') - slm_nglv = int(section.get('slm_nglv')) # get an instance of output analysis class - run_path = f'../mali{res}km_slm{slm_nglv}/run_model/' - run_data = output_analysis(config, logger, res, run_path) + run_path = f'../mali{res}km_slm{nglv}/run_model/' + run_data = output_analysis(config, logger, res, nglv, run_path) yrs = run_data.yrs # figure 1 @@ -169,7 +237,7 @@ def visualize_slm_circsheet(config, logger, res, nglv): ax1.set_ylabel('Mass (Gt)') ax1.grid(True) if save_images: - plt.savefig(f'ice_mass_mali{res}km_slm{slm_nglv}.png') + plt.savefig(f'ice_mass_mali{res}km_slm{nglv}.png') fig1.clf() # figure 2 @@ -188,7 +256,7 @@ def visualize_slm_circsheet(config, logger, res, nglv): ax2.set_ylabel('Mass (Gt)') ax2.grid(True) if save_images: - plt.savefig(f'ice_mass_change_mali{res}km_slm{slm_nglv}.png') + plt.savefig(f'ice_mass_change_mali{res}km_slm{nglv}.png') fig2.clf() # figure 3 @@ -217,7 +285,7 @@ def visualize_slm_circsheet(config, logger, res, nglv): ax3.set_ylabel('Sea-level change (m)') ax3.grid(True) if save_images: - plt.savefig(f'ice_contribution_to_SLC_mali{res}km_slm{slm_nglv}.png') + plt.savefig(f'ice_contribution_to_SLC_mali{res}km_slm{nglv}.png') fig3.clf() # figure 4 & 5 @@ -232,7 +300,7 @@ def visualize_slm_circsheet(config, logger, res, nglv): ax4.set_ylabel('Mass (Gt)') ax4.grid(True) if save_images: - plt.savefig(f'diff_ice_mass_change_mali{res}km_slm{slm_nglv}.png') + plt.savefig(f'diff_ice_mass_change_mali{res}km_slm{nglv}.png') fig4.clf() fig5 = plt.figure(5, facecolor='w') @@ -240,14 +308,14 @@ def visualize_slm_circsheet(config, logger, res, nglv): ax5.plot(yrs, run_data.diff_slc, label='MALI minus SLM', linestyle='-', color='k') ax5.plot(yrs, run_data.diff_slc_z0, label='MALI minus SLM (z0)', - linestyle='-', color='k') + linestyle='-', color='r') ax5.set_xlabel('Year') ax5.legend(loc='best', ncol=1, prop={'size': 10}) ax5.set_title('Difference in sea-level change') ax5.set_ylabel('sea-level change (m)') ax5.grid(True) if save_images: - plt.savefig(f'diff_sea_level_change_mali{res}km_slm{slm_nglv}.png') + plt.savefig(f'diff_sea_level_change_mali{res}km_slm{nglv}.png') fig5.clf() if hide_figs: @@ -360,45 +428,47 @@ def __init__(self, config, logger, res, nglv, run_path): time_stride = int(section.get('time_stride')) # reset the time indices to match the # of SLM timesteps nt = np.arange(0, DS.dims['Time'], time_stride, dtype='i') + nt = np.arange(0, 50, time_stride, dtype='i') self.yrs = self.yrs_mali[nt] - z0 = np.zeros((len(self.yrs), )) - Aocn = np.zeros((len(self.yrs), )) - AocnBeta = np.zeros((len(self.yrs), )) + nyrs = len(self.yrs) + z0 = np.zeros((nyrs, )) + Aocn = np.zeros((nyrs, )) + AocnBeta = np.zeros((nyrs, )) # get ice mass on the SLM interpolated from the MALI mesh fname = os.path.join(fpath_slm, 'ice_volume') - self.ice_vol_slm = slm_outputs(fname).data * rhoi / 1.0e12 - self.dice_vol_slm = slm_outputs(fname).change_total * \ + self.ice_vol_slm = slm_outputs(fname, nyrs).data * rhoi / 1.0e12 + self.dice_vol_slm = slm_outputs(fname, nyrs).change_total * \ rhoi / 1.0e12 # in Gt # get slc correction and ocean area from the SLM fname = os.path.join(fpath_slm, 'gmsle_deltaSL_Ocean_fixed') if os.path.exists(fname): logger.info(f'reading in file {fname}') - z0 = slm_outputs(fname).change_total - self.SLC_slm_Aocn = slm_outputs(fname).data + z0 = slm_outputs(fname, nyrs).change_total + self.SLC_slm_Aocn = slm_outputs(fname, nyrs).data fname = os.path.join(fpath_slm, 'gmsle_deltaSL_OceanBeta_fixed') if os.path.exists(fname): logger.info(f'reading in file {fname}') - self.SLC_slm_AocnBeta = slm_outputs(fname).data + self.SLC_slm_AocnBeta = slm_outputs(fname, nyrs).data fname = os.path.join(fpath_slm, 'ocean_area') if os.path.exists(fname): logger.info(f'reading in file {fname}') - Aocn = slm_outputs(fname).data - logger.info(f'area of the ocean is: {Aocn}') + Aocn = slm_outputs(fname, nyrs).data + # logger.info(f'area of the ocean is: {Aocn}') fname = os.path.join(fpath_slm, 'oceanBeta_area') if os.path.exists(fname): logger.info(f'reading in file {fname}') - AocnBeta = slm_outputs(fname).data - logger.info(f'area of the ocean Beta is: {AocnBeta}') + AocnBeta = slm_outputs(fname, nyrs).data + # logger.info(f'area of the ocean Beta is: {AocnBeta}') else: # if not coupled, use default values for MALI outputs # assuming MALI output interval is 1 year nt = np.arange(0, DS.dims['Time'], 1, dtype='i') self.yrs = self.yrs_mali - z0 = np.zeros((len(self.yrs), )) - Aocn = np.zeros((len(self.yrs), )) - AocnBeta = np.zeros(len(self.yrs, )) + z0 = np.zeros(nyrs, ) + Aocn = np.zeros(nyrs, ) + AocnBeta = np.zeros(nyrs, ) z0[:] = 0.0 logger.info("'gmsle_change' file doesn't exist. " @@ -447,8 +517,8 @@ def __init__(self, config, logger, res, nglv, run_path): bedt = bed[indx, :].load() cellMaskt = cellMask[indx, :].load() thicknesst = thickness[indx, :].load() - areaCellt = areaCell.sum() - logger.info(f'area of the mali domain is {areaCellt}') + # areaCellt = areaCell.sum() + # logger.info(f'area of the mali domain is {areaCellt}') # calculation of sea-level change following Goelzer et al. 2020 self.grnd_vol_unscaled[t] = (areaCell * grounded(cellMaskt) * @@ -515,12 +585,25 @@ def __init__(self, config, logger, res, nglv, run_path): # calculate RMSE between MALI and SLM calculation of SLC if coupling: - self.diff_slc = self.SLCcorr_AocnBeta - self.SLC_slm_AocnBeta - self.diff_slc_z0 = self.SLCcorr_z0_AocnBeta - self.SLC_slm_AocnBeta + self.diff_slc = self.SLC_slm_AocnBeta - self.SLCcorr_AocnBeta + self.diff_slc_z0 = self.SLC_slm_AocnBeta - \ + self.SLCcorr_z0_AocnBeta self.rmse_slc = np.sqrt((self.diff_slc**2).mean()) + self.pct_err_slc = ((abs(self.diff_slc[1:]) / + abs(self.SLCcorr_AocnBeta[1:])) * + 100).mean() + self.pct_err_slc_z0 = ((abs(self.diff_slc_z0[1:]) / + abs(self.SLCcorr_z0_AocnBeta[1:])) * + 100).mean() + self.pct_err_slm = ((abs(self.diff_slc[1:]) / + abs(self.SLC_slm_AocnBeta[1:])) * + 100).mean() + self.pct_err_slm_z0 = ((abs(self.diff_slc_z0[1:]) / + abs(self.SLC_slm_AocnBeta[1:])) * + 100).mean() + else: - logger.info('MALI is not coupled to the SLM. No RMSE to compute.') - self.rmse_slc = np.nan + logger.info('MALI is not coupled to the SLM. No MPE to compute.') def grounded(cellMask): @@ -544,7 +627,7 @@ class slm_outputs: Change in data value at a time step with respect to the intial time step """ - def __init__(self, filename): + def __init__(self, filename, nyrs): """ Calculate the sea-level model output file @@ -552,11 +635,14 @@ def __init__(self, filename): ---------- filename : str Filename of SLM output data + + nyrs : int + Number of time indices to be analyzed """ self.fname = filename - self.data = np.loadtxt(self.fname) - self.change_dt = np.zeros(len(self.data),) - self.change_total = np.zeros(len(self.data),) + self.data = np.loadtxt(self.fname)[0:nyrs] + self.change_dt = np.zeros(nyrs,) + self.change_total = np.zeros(nyrs,) for i in range(len(self.change_dt) - 1): self.change_dt[i + 1] = (float(self.data[i + 1]) - float(self.data[i])) From 1a48baac47007510196dc02210b824fb39ea0eb6 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Sat, 18 Nov 2023 17:59:22 -0700 Subject: [PATCH 19/23] Change the name of the SLM GMSLC files --- compass/landice/tests/slm/circ_icesheet/visualize.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/visualize.py b/compass/landice/tests/slm/circ_icesheet/visualize.py index f9f63f70e0..b029f4ecde 100644 --- a/compass/landice/tests/slm/circ_icesheet/visualize.py +++ b/compass/landice/tests/slm/circ_icesheet/visualize.py @@ -440,13 +440,13 @@ def __init__(self, config, logger, res, nglv, run_path): self.dice_vol_slm = slm_outputs(fname, nyrs).change_total * \ rhoi / 1.0e12 # in Gt # get slc correction and ocean area from the SLM - fname = os.path.join(fpath_slm, 'gmsle_deltaSL_Ocean_fixed') + fname = os.path.join(fpath_slm, 'gmslc_ocnArea') if os.path.exists(fname): logger.info(f'reading in file {fname}') z0 = slm_outputs(fname, nyrs).change_total self.SLC_slm_Aocn = slm_outputs(fname, nyrs).data - fname = os.path.join(fpath_slm, 'gmsle_deltaSL_OceanBeta_fixed') + fname = os.path.join(fpath_slm, 'gmslc_ocnBetaArea') if os.path.exists(fname): logger.info(f'reading in file {fname}') self.SLC_slm_AocnBeta = slm_outputs(fname, nyrs).data From 65b1061e2f45de4e14d682e5b1acc53cb8b8fde6 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Sat, 18 Nov 2023 18:01:16 -0700 Subject: [PATCH 20/23] Assign variable ntask values depending on MALI res --- .../landice/tests/slm/circ_icesheet/__init__.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/__init__.py b/compass/landice/tests/slm/circ_icesheet/__init__.py index ab3d8bbb60..6315ef6b0f 100644 --- a/compass/landice/tests/slm/circ_icesheet/__init__.py +++ b/compass/landice/tests/slm/circ_icesheet/__init__.py @@ -52,9 +52,18 @@ def configure(self): self.add_step(SetupMesh(test_case=self, name=f'mali{res}km_slm{nglv}/' f'setup_mesh', res=res, nglv=nglv)) - self.add_step(RunModel(test_case=self, ntasks=1, - openmp_threads=1, - name=f'mali{res}km_slm{nglv}/run_model')) + if (int(res) <= 16 and int(res) > 2): + ntasks = 256 + elif (int(res) <= 2): + ntasks = 512 + else: + ntasks = 128 + min_tasks = ntasks + self.add_step(RunModel(test_case=self, res=res, nglv=nglv, + ntasks=ntasks, min_tasks=min_tasks, + openmp_threads=1, + name=f'mali{res}km_slm{nglv}' + '/run_model')) step = Visualize(test_case=self) self.add_step(step, run_by_default=True) From 288c1786dca9bc7b4f5100a3a2d1b3824b21570a Mon Sep 17 00:00:00 2001 From: hollyhan Date: Sat, 18 Nov 2023 21:14:35 -0700 Subject: [PATCH 21/23] Add and edit the SLM namelist template --- .../tests/slm/circ_icesheet/run_model.py | 17 ++++- .../tests/slm/namelist.sealevel.template | 66 +++++++++++++++++++ 2 files changed, 80 insertions(+), 3 deletions(-) create mode 100644 compass/landice/tests/slm/namelist.sealevel.template diff --git a/compass/landice/tests/slm/circ_icesheet/run_model.py b/compass/landice/tests/slm/circ_icesheet/run_model.py index db3a0b76d4..d87888bacf 100644 --- a/compass/landice/tests/slm/circ_icesheet/run_model.py +++ b/compass/landice/tests/slm/circ_icesheet/run_model.py @@ -1,4 +1,7 @@ import os +from importlib import resources + +from jinja2 import Template from compass.model import run_model from compass.step import Step @@ -66,9 +69,6 @@ def __init__(self, test_case, res, nglv, ntasks, name='run_model', target='../setup_mesh/landice_grid.nc') self.add_input_file(filename='graph.info', target='../setup_mesh/graph.info') - self.add_input_file(filename='namelist.sealevel', - package='compass.landice.tests.slm', - copy=True) self.add_input_file(filename='smb_forcing.nc', target='../setup_mesh/smb_forcing.nc') self.add_input_file(filename='mapping_file_mali_to_slm.nc', @@ -87,6 +87,17 @@ def setup(self): os.makedirs(os.path.join(self.work_dir, 'ICELOAD_SLM/'), exist_ok='True') + # change the sealevel namelist + template = Template(resources.read_text + ('compass.landice.tests.slm', + 'namelist.sealevel.template')) + text = template.render(nglv=self.nglv) + + # write out the namelise.sealevel file + file_slm_nl = os.path.join(self.work_dir, 'namelist.sealevel') + with open(file_slm_nl, 'w') as handle: + handle.write(text) + def run(self): """ Run this step of the test case diff --git a/compass/landice/tests/slm/namelist.sealevel.template b/compass/landice/tests/slm/namelist.sealevel.template new file mode 100644 index 0000000000..4676e62bdc --- /dev/null +++ b/compass/landice/tests/slm/namelist.sealevel.template @@ -0,0 +1,66 @@ +&time_config + itersl = 1 + starttime = 2015 + dt1 = 1 + +/ +&io_directory + inputfolder_ice = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/icemodel/GL{{ nglv }}/ice_noGrIS_GL{{ nglv }}/' + inputfolder = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/others/GL{{ nglv }}/' + planetfolder = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/earthmodel/' + gridfolder = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/others/GL{{ nglv }}/' + outputfolder = 'OUTPUT_SLM/' + outputfolder_ice = 'ICELOAD_SLM/' + folder_coupled = '' + +/ +&file_format + ext ='' + fType_in = 'text' + fType_out = 'both' + +/ +&file_name + planetmodel = 'prem_512.l60K2C.sum18p6.dum19p2.tz19p4.lm22' + icemodel = 'iceGlobalDomain_zeroField_GL{{ nglv }}_' + icemodel_out = 'iceload_out_' + timearray = 'times' + topomodel = 'idealized_topo_1500mOcean_SH_GL{{ nglv }}' + topo_initial = 'idealized_topo_1500mOcean_SH_GL{{ nglv }}' + grid_lat = 'GLlat_{{ nglv }}.txt' + grid_lon = 'GLlon_{{ nglv }}.txt' + +/ +&model_config + checkmarine = .false. + tpw = .true. + calcRG = .true. + input_times = .false. + initial_topo = .true. + iceVolume = .true. + coupling = .true. + patch_ice = .true. + +/ +&model_resolution + norder = 512 + nglv = {{ nglv }} + +/ +&timewindow_config + L_sim = 100 + dt1 = 1 + dt2 = 10 + dt3 = 10 + dt4 = 10 + Ldt1 = 100 + Ldt2 = 0 + Ldt3 = 0 + Ldt4 = 0 + +/ +&others + whichplanet = 'earth' + + +/ From 84fc4a6b3a587e41f275299122a6115e54101dd0 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Sat, 31 Aug 2024 11:43:12 -0600 Subject: [PATCH 22/23] clean up some formatting issues --- compass/landice/tests/slm/circ_icesheet/__init__.py | 10 ++-------- .../landice/tests/slm/circ_icesheet/circ_icesheet.cfg | 7 ++++--- compass/landice/tests/slm/circ_icesheet/run_model.py | 6 ++---- 3 files changed, 8 insertions(+), 15 deletions(-) diff --git a/compass/landice/tests/slm/circ_icesheet/__init__.py b/compass/landice/tests/slm/circ_icesheet/__init__.py index 6315ef6b0f..75b373cd7f 100644 --- a/compass/landice/tests/slm/circ_icesheet/__init__.py +++ b/compass/landice/tests/slm/circ_icesheet/__init__.py @@ -9,9 +9,6 @@ class CircIcesheetTest(TestCase): This test generates an idealized, circular ice sheet that has a prescribed thickness evolution for testing coupling between MALI and the Sea-Level Model. - - Attributes - ---------- """ def __init__(self, test_group): @@ -20,12 +17,9 @@ def __init__(self, test_group): Parameters ---------- - test_group : compass.landice.tests.dome.Dome + test_group : compass.landice.tests.slm.Slm The test group that this test case belongs to The resolution or type of mesh of the test case - - name : str, optional - the name of the test case """ name = 'circular_icesheet_test' subdir = name @@ -51,7 +45,7 @@ def configure(self): for nglv in slm_nglv: self.add_step(SetupMesh(test_case=self, name=f'mali{res}km_slm{nglv}/' - f'setup_mesh', res=res, nglv=nglv)) + 'setup_mesh', res=res, nglv=nglv)) if (int(res) <= 16 and int(res) > 2): ntasks = 256 elif (int(res) <= 2): diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg index b2b7ebd134..8525d54101 100644 --- a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -16,9 +16,10 @@ r0 = 2000000.0 h0 = 3000.0 # 'True' if manually want to set bedTopography elevation -# bed topography elevation in meters (provide a positive/negative float value -# for grounded above land/marine-based ice) +# if 'False' topography will be set to 0 everywhere (flat bed at sea level) set_topo_elev = True +# flat bed topography elevation in meters +# (positive value => land, negative => ocean) topo_elev = -1500.0 # Whether to center the circular ice in the center of the cell that is closest to the @@ -47,7 +48,7 @@ dhdt = -20.0 coupling = True # List of the number of Gauss-Legendre points in latitude -# delimited by ',' without space +# list delimited by ',' without space slm_nglv = 512 # Max spherical harmonics degree and order diff --git a/compass/landice/tests/slm/circ_icesheet/run_model.py b/compass/landice/tests/slm/circ_icesheet/run_model.py index d87888bacf..1f9084a658 100644 --- a/compass/landice/tests/slm/circ_icesheet/run_model.py +++ b/compass/landice/tests/slm/circ_icesheet/run_model.py @@ -10,10 +10,8 @@ class RunModel(Step): """ A step for performing forward MALI runs as part of dome test cases. - - Attributes - ---------- """ + def __init__(self, test_case, res, nglv, ntasks, name='run_model', subdir=None, min_tasks=None, openmp_threads=1): """ @@ -93,7 +91,7 @@ def setup(self): 'namelist.sealevel.template')) text = template.render(nglv=self.nglv) - # write out the namelise.sealevel file + # write out the namelist.sealevel file file_slm_nl = os.path.join(self.work_dir, 'namelist.sealevel') with open(file_slm_nl, 'w') as handle: handle.write(text) From 27ad364f3433cb171a98bede3d15d9d21fe95ec3 Mon Sep 17 00:00:00 2001 From: hollyhan Date: Sat, 31 Aug 2024 16:24:03 -0600 Subject: [PATCH 23/23] Fix solidearth nl option dt from string to integer --- compass/landice/tests/slm/circ_icesheet/namelist.landice | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/landice/tests/slm/circ_icesheet/namelist.landice b/compass/landice/tests/slm/circ_icesheet/namelist.landice index e654cc6ee1..37851f0e0d 100644 --- a/compass/landice/tests/slm/circ_icesheet/namelist.landice +++ b/compass/landice/tests/slm/circ_icesheet/namelist.landice @@ -7,7 +7,7 @@ / &solidearth config_uplift_method = 'sealevelmodel' - config_slm_coupling_interval = '0001-00-00_00:00:00' + config_slm_coupling_interval = 5 config_MALI_to_SLM_weights_file = 'mapping_file_mali_to_slm.nc' config_SLM_to_MALI_weights_file = 'mapping_file_slm_to_mali.nc' /