Run ISSM over multiple input files (multiple ice sheets or glaciers)

examples/GEOSInput/Antarctica_input_control.py

@@ -0,0 +1,103 @@
+import devpath
+import os
+ISSM_DIR = os.getenv('ISSM_DIR') # for binaries
+
+import numpy as np
+from triangle import triangle
+from model import *
+from netCDF4 import Dataset
+from InterpFromGridToMesh import InterpFromGridToMesh
+from bamg import bamg
+from loadmodel import loadmodel
+from setmask import setmask
+from parameterize import parameterize
+from clusters.discover_geos import export_discover
+from marshall import marshall
+
+if not os.path.exists('./Models'):
+    os.mkdir('./Models')
+
+# Step 1: Mesh generation
+# Generate initial uniform mesh (resolution = 60000 m)
+# project mesh onto new coordinate system
+md = triangle(model(), '../Data/Ais.exp', 60000)
+
+print('   Loading velocities data from NetCDF')
+nsidc_vel = Dataset('../Data/Antarctica_ice_velocity.nc')
+xmin = nsidc_vel.xmin
+xmin = float(xmin.lstrip()[0:10])
+ymax = nsidc_vel.ymax
+ymax = float(ymax.lstrip()[0:10])
+spacing = nsidc_vel.spacing
+spacing = float((spacing.lstrip())[0:4])
+nx = nsidc_vel.nx
+ny = nsidc_vel.ny
+velx = nsidc_vel['vx'][:].data
+vely = nsidc_vel['vy'][:].data
+# Build coordinates
+x2 = xmin + np.arange(nx + 1) * spacing
+y2 = (ymax - ny * spacing) + np.arange(ny + 1) * spacing
+
+# print('   Set observed velocities')
+md.initialization.vx = InterpFromGridToMesh(x2, y2, np.flipud(velx), md.mesh.x, md.mesh.y, 0)
+md.initialization.vy = InterpFromGridToMesh(x2, y2, np.flipud(vely), md.mesh.x, md.mesh.y, 0)
+md.initialization.vz = np.zeros(md.mesh.numberofvertices)
+md.initialization.vel = np.sqrt(md.initialization.vx**2 + md.initialization.vy**2)
+del velx, vely
+
+md = bamg(md, 'hmax', 75000, 'hmin', 5000, 'gradation', 1.4, 'field', md.initialization.vel, 'err', 8)
+
+# export mesh
+export_discover(md, './Models/Antarctica_mesh.nc')
+
+# Step 2: parameterize model
+md = loadmodel('./Models/Antarctica_mesh.nc')
+
+md = setmask(md, '', '')
+md = parameterize(md, './Antarctica_parameterize.py')
+
+
+# export parameterization
+export_discover(md, "./Models/Antarctica_parameterization.nc",delete_rundir=True)
+
+# Step 3: basal friction inversion 
+	# Control general
+md.inversion.nsteps = 400
+md.inversion.iscontrol=1
+md.inversion.maxsteps=400
+md.inversion.maxiter=400
+md.inversion.dxmin=0.01
+md.inversion.gttol=1.0e-8
+
+md.inversion.step_threshold = 0.99 * np.ones((md.inversion.nsteps))
+md.inversion.maxiter_per_step = 40 * np.ones((md.inversion.nsteps))
+
+md.inversion.gradient_scaling = 50 * np.ones((md.inversion.nsteps, 1))
+md.inversion.min_parameters = 1 * np.ones((md.mesh.numberofvertices, 1))
+md.inversion.max_parameters = 200 * np.ones((md.mesh.numberofvertices, 1))
+
+#Cost functions
+md.inversion.cost_functions = [101, 103, 501]
+md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, 3))
+md.inversion.cost_functions_coefficients[:, 0] = 1
+md.inversion.cost_functions_coefficients[:, 1] = 1
+md.inversion.cost_functions_coefficients[:, 2] = 2e-10
+
+# Controls
+md.inversion.control_parameters = ['FrictionCoefficient']
+md.inversion.min_parameters=1*np.ones(np.shape(md.mesh.x))
+md.inversion.max_parameters=200*np.ones(np.shape(md.mesh.x))
+
+# Additional parameters
+md.stressbalance.restol=0.0000001
+md.stressbalance.reltol=0.0000001
+md.stressbalance.abstol=np.nan
+
+# Solve
+md.private.solution = 'Stressbalance'
+md.settings.waitonlock = 0
+md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits')
+marshall(md) # create .bin file
+
+# export configuration for loading solution in next step
+export_discover(md,'./Models/Antarctica_inversion.nc',delete_rundir=True)

examples/GEOSInput/Antarctica_input_transient.py

@@ -0,0 +1,29 @@
+import devpath
+import os,sys
+ISSM_DIR = os.getenv('ISSM_DIR') # for binaries
+
+from model import *
+from loadmodel import loadmodel
+from clusters.discover_geos import export_discover
+from loadresultsfromdisk import loadresultsfromdisk
+from marshall import marshall
+from verbose import verbose 
+
+md = loadmodel('./Models/Antarctica_inversion.nc')
+md = loadresultsfromdisk(md,'AntarcticaGEOS.outbin')
+md.friction.coefficient = md.results.StressbalanceSolution.FrictionCoefficient
+
+
+# Write the binary input file
+# Additional options
+md.inversion.iscontrol = 0
+md.transient.requested_outputs = ['default']
+md.transient.isthermal=0
+md.settings.waitonlock = 0
+md.private.solution = 'Transient'
+md.verbose = verbose('000000000')
+md.toolkits = toolkits()
+marshall(md) # create .bin file
+md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits')
+export_discover(md,'./Models/Antarctica_initialization.nc',delete_rundir=True)
+

examples/GEOSInput/Antarctica_parameterize.py

@@ -0,0 +1,145 @@
+import numpy as np
+from paterson import paterson
+from netCDF4 import Dataset
+from xy2ll import xy2ll
+from InterpFromGridToMesh import InterpFromGridToMesh
+from SetMarineIceSheetBC import SetMarineIceSheetBC
+from m1qn3inversion import m1qn3inversion 
+from setflowequation import setflowequation
+import os
+ISSM_DIR = os.getenv('ISSM_DIR')
+
+
+#Name and Coordinate system
+md.miscellaneous.name='AntarcticaGEOS'
+md.mesh.epsg=3031
+
+nsidc_vel = Dataset('../Data/Antarctica_ice_velocity.nc')
+xmin = nsidc_vel.xmin
+xmin = float(xmin.lstrip()[0:10])
+ymax = nsidc_vel.ymax
+ymax = float(ymax.lstrip()[0:10])
+spacing = nsidc_vel.spacing
+spacing = float((spacing.lstrip())[0:4])
+nx = nsidc_vel.nx
+ny = nsidc_vel.ny
+velx = nsidc_vel['vx'][:].data
+vely = nsidc_vel['vy'][:].data
+# Build coordinates
+x2 = xmin + np.arange(nx + 1) * spacing
+y2 = (ymax - ny * spacing) + np.arange(ny + 1) * spacing
+
+# print('   Set observed velocities')
+md.initialization.vx = InterpFromGridToMesh(x2, y2, np.flipud(velx), md.mesh.x, md.mesh.y, 0)
+md.initialization.vy = InterpFromGridToMesh(x2, y2, np.flipud(vely), md.mesh.x, md.mesh.y, 0)
+md.initialization.vz = np.zeros(md.mesh.numberofvertices)
+md.initialization.vel = np.sqrt(md.initialization.vx**2 + md.initialization.vy**2)
+del velx, vely
+
+# Parameters to change/Try
+friction_coefficient = 10  # default [10]
+Temp_change = 0           # default [0 K]
+
+# NetCDF Loading
+print('   Loading SeaRISE data from NetCDF')
+ncdata = Dataset('../Data/Antarctica_5km_withshelves_v0.75.nc')
+x1 = ncdata['x1'][:].data
+y1 = ncdata['y1'][:].data
+usrf = ncdata['usrf'][:].data[0]
+topg = ncdata['topg'][:].data[0]
+temp = ncdata['presartm'][:].data[0]
+smb = ncdata['presprcp'][:].data[0]
+gflux = ncdata['bheatflx_fox'][:].data[0]
+
+# Geometry
+print('   Interpolating surface and ice base')
+md.geometry.base = InterpFromGridToMesh(x1, y1, topg, md.mesh.x, md.mesh.y, 0)
+md.geometry.surface = InterpFromGridToMesh(x1, y1, usrf, md.mesh.x, md.mesh.y, 0)
+del usrf, topg
+
+thkmask=ncdata['thkmask'][:].data[0]
+
+##interpolate onto our mesh vertices
+groundedice=   InterpFromGridToMesh(x1,y1,thkmask,md.mesh.x,md.mesh.y,0)
+groundedice[groundedice<=0]=-1
+del thkmask
+
+#fill in the md.mask structure
+md.mask.ocean_levelset = groundedice  #ice is grounded for mask equal one
+md.mask.ice_levelset = -1*np.ones(np.shape(md.mesh.x))  #ice is present when negatvie
+
+print('   Constructing thickness')
+md.geometry.thickness = md.geometry.surface - md.geometry.base
+
+# Ensure hydrostatic equilibrium on ice shelf
+di = md.materials.rho_ice / md.materials.rho_water
+
+# Get the node numbers of floating nodes
+pos = np.where(md.mask.ocean_levelset < 0)
+
+# Apply flotation criterion
+md.geometry.thickness[pos] = 1 / (1 - di) * md.geometry.surface[pos]
+md.geometry.base[pos] = md.geometry.surface[pos] - md.geometry.thickness[pos]
+md.geometry.hydrostatic_ratio = np.ones(md.mesh.numberofvertices)  
+
+# Set min thickness to 1 meter
+pos0 = np.where(md.geometry.thickness <= 1)
+md.geometry.thickness[pos0] = 1
+md.geometry.surface = md.geometry.thickness + md.geometry.base
+md.geometry.bed = md.geometry.base.copy()
+md.geometry.bed[pos] = md.geometry.base[pos] - 1000
+
+
+# Initialization parameters
+print('   Interpolating temperatures')
+md.initialization.temperature = InterpFromGridToMesh(
+    x1, y1, temp, md.mesh.x, md.mesh.y, 0
+) + 273.15 + Temp_change
+
+[md.mesh.lat, md.mesh.long] = xy2ll(md.mesh.x, md.mesh.y, -1)
+
+
+print('   Set Pressure')
+md.initialization.pressure = md.materials.rho_ice * md.constants.g * md.geometry.thickness
+
+print('   Construct ice rheological properties')
+md.materials.rheology_n = 3 * np.ones(md.mesh.numberofelements)
+md.materials.rheology_B = paterson(md.initialization.temperature)
+
+# Forcings
+print('   Interpolating surface mass balance')
+mass_balance = InterpFromGridToMesh(x1, y1, smb, md.mesh.x, md.mesh.y, 0)
+md.smb.mass_balance = mass_balance * md.materials.rho_water / md.materials.rho_ice
+
+print('   Set geothermal heat flux')
+md.basalforcings.geothermalflux = InterpFromGridToMesh(x1, y1, gflux, md.mesh.x, md.mesh.y, 0)
+
+# Friction and inversion set up
+print('   Construct basal friction parameters')
+md.friction.coefficient = friction_coefficient * np.ones(md.mesh.numberofvertices)
+md.friction.p = np.ones(md.mesh.numberofelements)
+md.friction.q = np.ones(md.mesh.numberofelements)
+
+# No friction applied on floating ice
+pos = np.where(md.mask.ocean_levelset < 0)[0]
+md.friction.coefficient[pos] = 0
+md.groundingline.migration = 'SubelementMigration'
+
+md.inversion = m1qn3inversion()
+md.inversion.vx_obs = md.initialization.vx
+md.inversion.vy_obs = md.initialization.vy
+md.inversion.vel_obs = md.initialization.vel
+
+print('   Set flow equations')
+md = setflowequation(md,'SSA','all')
+
+print('   Set boundary conditions')
+md = SetMarineIceSheetBC(md)
+md.basalforcings.floatingice_melting_rate = np.zeros(md.mesh.numberofvertices)
+md.basalforcings.groundedice_melting_rate = np.zeros(md.mesh.numberofvertices)
+md.thermal.spctemperature = md.initialization.temperature
+md.masstransport.spcthickness = np.full(md.mesh.numberofvertices, np.nan)
+
+
+
+

examples/GEOSInput/Greenland_input_control.py

@@ -9,12 +9,10 @@
 from InterpFromGridToMesh import InterpFromGridToMesh
 from bamg import bamg
 from xy2ll import xy2ll
-from plotmodel import plotmodel
 from loadmodel import loadmodel
 from setmask import setmask
 from parameterize import parameterize
 from setflowequation import setflowequation
-from solve import solve
 from ll2xy import ll2xy
 from clusters.discover_geos import export_discover
 from marshall import marshall

examples/GEOSInput/Greenland_input_transient.py

@@ -2,7 +2,6 @@
 import os,sys
 ISSM_DIR = os.getenv('ISSM_DIR') # for binaries
 
-import numpy as np
 from model import *
 from loadmodel import loadmodel
 from clusters.discover_geos import export_discover
@@ -18,7 +17,8 @@
 # Write the binary input file
 # Additional options
 md.inversion.iscontrol = 0
-md.transient.requested_outputs = ['IceVolume', 'TotalSmb', 'SmbMassBalance']
+md.transient.requested_outputs = ['default']
+md.transient.isthermal=0
 md.settings.waitonlock = 0
 md.private.solution = 'Transient'
 md.verbose = verbose('000000000')

examples/GEOSInput/generate_antarctica_input.sh

@@ -0,0 +1,5 @@
+source issm_env_withpy
+rm -f AntarcticaGEOS.bin AntarcticaGEOS.outbin AntarcticaGEOS.outlog AntarcticaGEOS.errlog
+(export LD_LIBRARY_PATH=$PYTHON_LIB:$LD_LIBRARY_PATH; python Antarctica_input_control.py)
+${ISSM_DIR}/bin/issm.exe StressbalanceSolution $(pwd) AntarcticaGEOS 2>> AntarcticaGEOS.errlog
+(export LD_LIBRARY_PATH=$PYTHON_LIB:$LD_LIBRARY_PATH; python Antarctica_input_transient.py)

examples/GEOSInput/generate_geos_input.sh → examples/GEOSInput/generate_greenland_input.sh

@@ -1,5 +1,4 @@
 source issm_env_withpy
-rm -rf Models
 rm -f GreenlandGEOS.bin GreenlandGEOS.outbin GreenlandGEOS.outlog GreenlandGEOS.errlog
 (export LD_LIBRARY_PATH=$PYTHON_LIB:$LD_LIBRARY_PATH; python Greenland_input_control.py)
 ${ISSM_DIR}/bin/issm.exe StressbalanceSolution $(pwd) GreenlandGEOS 2>> GreenlandGEOS.errlog