Source code for compass.landice.tests.hydro_radial.setup_mesh

import numpy as np
from netCDF4 import Dataset as NetCDFFile

from mpas_tools.planar_hex import make_planar_hex_mesh
from mpas_tools.io import write_netcdf
from mpas_tools.mesh.conversion import convert, cull
from mpas_tools.logging import check_call

from compass.model import make_graph_file
from compass.step import Step


[docs] class SetupMesh(Step): """ A step for creating a mesh and initial condition for dome test cases Attributes ---------- initial_condition : {'zero', 'exact'} The type of initial condition to set up. 'zero' means nearly zero ice thickness. 'exact' uses a precomputed near exact solution from a file. """
[docs] def __init__(self, test_case, initial_condition): """ Create the step Parameters ---------- test_case : compass.TestCase The test case this step belongs to initial_condition : {'zero', 'exact'} The type of initial condition to set up. 'zero' means nearly zero ice thickness. 'exact' uses a precomputed near exact solution from a file. """ super().__init__(test_case=test_case, name='setup_mesh') self.initial_condition = initial_condition if initial_condition == 'exact': self.add_input_file(filename='near_exact_solution_r_P_W.txt', package='compass.landice.tests.hydro_radial') elif initial_condition != 'zero': raise ValueError("Unknown initial condition type specified " "{}.".format(initial_condition)) self.add_output_file(filename='graph.info') self.add_output_file(filename='landice_grid.nc')
# no setup() method is needed
[docs] def run(self): """ Run this step of the test case """ initial_condition = self.initial_condition logger = self.logger section = self.config['hydro_radial'] 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) write_netcdf(dsMesh, 'grid.nc') dsMesh = cull(dsMesh, logger=logger) dsMesh = convert(dsMesh, logger=logger) write_netcdf(dsMesh, 'mpas_grid.nc') levels = section.get('levels') args = ['create_landice_grid_from_generic_MPAS_grid.py', '-i', 'mpas_grid.nc', '-o', 'landice_grid.nc', '-l', levels, '--hydro', '--diri'] check_call(args, logger) make_graph_file(mesh_filename='landice_grid.nc', graph_filename='graph.info') _setup_hydro_radial_initial_conditions( logger, filename='landice_grid.nc', initial_condition=initial_condition)
def _setup_hydro_radial_initial_conditions(logger, filename, initial_condition): """ Add the initial condition to the given MPAS mesh file Parameters ---------- logger : logging.Logger A logger for output from the step filename : str file to setup hydro_radial initial_condition : {'zero', 'exact'} the type of initial condition """ # 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'] # 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 # Center the dome in the center of the cell that is closest to the center # of the domain. # NOTE: for some meshes, maybe we don't want to do this - could add # command-line argument controlling this later. putOriginOnACell = True if putOriginOnACell: centerCellIndex = np.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 # center thickness (m) h0 = 500.0 # sliding velocity at margin (m/s) v0 = 100.0 / (3600.0 * 24.0 * 365.0) # ideal ice cap radius (m) R0 = 25.0e3 # onset of sliding (m) R1 = 5.0e3 # actual margin location (m) L = 0.9 * R0 thickness[0, r < R0] = h0 * (1.0 - (r[r < R0] / R0)**2) thickness[0, r > L] = 0.0 # flat bed bedTopography[:] = 0.0 # Setup layerThicknessFractions layerThicknessFractions[:] = 1.0 / nVertLevels # melt gridfile.variables['basalMeltInput'][:] = 0.0 # 20 cm/yr as SI mass rate gridfile.variables['basalMeltInput'][:] = \ 0.2 / (365.0 * 24.0 * 3600.0) * 1000.0 # Use this line to only add a source term to the center cell - useful for # debugging divergence # value from ramp # gridfile.variables['basalMeltInput'][0,r==0.0] = 4.0e-10 * 1000.0 *100 # velocity gridfile.variables['uReconstructX'][:] = 0.0 velo = v0 * (r - R1)**5 / (L - R1)**5 velo[r < R1] = 0.0 gridfile.variables['uReconstructX'][0, :, -1] = velo gridfile.variables['uReconstructX'][0, thickness[0, :] == 0.0, :] = 0.0 if initial_condition == 'zero': logger.info("Using 'zero' option for initial condition.") # set some small initial value to keep adaptive time stepper from # taking a huge time step initially gridfile.variables['waterThickness'][0, :] = 0.01 gridfile.variables['waterPressure'][0, :] = 0.0 elif initial_condition == 'exact': logger.info("Using 'exact' option for initial condition.") # IC on thickness # import exact solution fnameSoln = 'near_exact_solution_r_P_W.txt' soln = np.loadtxt(fnameSoln, delimiter=',') rsoln = soln[:, 0] Psoln = soln[:, 1] Wsoln = soln[:, 2] Wmpas = np.interp(r, rsoln, Wsoln) # apply exact solution Wmpas[np.isnan(Wmpas)] = 0.0 gridfile.variables['waterThickness'][0, :] = Wmpas # IC on water pressure # apply exact solution Pmpas = np.interp(r, rsoln, Psoln) Pmpas[np.isnan(Pmpas)] = 0.0 gridfile.variables['waterPressure'][0, :] = Pmpas else: raise ValueError("Unknown initial condition type specified " "{}.".format(initial_condition)) gridfile.close() logger.info('Successfully added hydro_radial initial conditions to: ' '{}'.format(filename))