<p><b>mhoffman@lanl.gov</b> 2013-03-05 20:14:35 -0700 (Tue, 05 Mar 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
A script for calibrating the basal traction parameter field, beta, using observed (or balance) surface velocities. This is roughly the method described in Price et al 2011 (PNAS) with CISM. In the long term we want to use a control method, but this method provides a reasonable alternative in the meantime. <br>
</p><hr noshade><pre><font color="gray">Added: branches/land_ice_projects/grid_tools/calibrate_beta.py
===================================================================
--- branches/land_ice_projects/grid_tools/calibrate_beta.py (rev 0)
+++ branches/land_ice_projects/grid_tools/calibrate_beta.py 2013-03-06 03:14:35 UTC (rev 2549)
@@ -0,0 +1,131 @@
+#!/usr/bin/python
+# Script for calibrating beta using MPAS
+# Matt Hoffman, March 2013
+
+import os, sys
+import subprocess
+import numpy as np
+#from optparse import OptionParser
+#import matplotlib.pyplot as plt
+try:
+ from Scientific.IO.NetCDF import NetCDFFile
+ netCDF_module = 'Scientific.IO.NetCDF'
+except ImportError:
+ try:
+ from netCDF4 import Dataset as NetCDFFile
+ netCDF_module = 'netCDF4'
+ except ImportError:
+ print 'Unable to import any of the following python modules:'
+ print ' Scientific.IO.NetCDF </font>
<font color="blue"> netcdf4 '
+ print 'One of them must be installed.'
+ raise ImportError('No netCDF module found')
+
+
+################ SET THESE FOR YOUR SETUP ###########
+runcmd = 'time openmpirun -n 4 land_ice_model.exe' # command to run MPAS on your system
+vobsfile='land_ice_grid.nc.balvel' # This is the file that has the observed surface velocity field - it can have anything else, too.
+vobsvariable='vsurfOb' # this is the name of the observed velocity field variable in the above file. It should have dimensions (Time, nCells)
+numiters = 3 # Number of iterations to run the algorithm. the first time through uses a constant no-sliding beta to get the first guess of the deformation velocity
+#####################################################
+
+# These are probably fine.
+infile='land_ice_grid.nc'
+outfile='output.nc'
+modelfiles='namelist.input graph.info.part* dataIce *xml ' + infile # land_ice_model.exe gets symlinked
+
+# remove previous attempts
+os.system('rm -rf iter*')
+
+
+# Loop and iterate these steps -
+# (add relaxation or smarter guess later)
+for i in range(numiters):
+ print 'starting iteration ' + str(i)
+
+ ########### Calculate a guess for beta ########
+ if i == 0:
+ # Setup the directory for the first run
+ os.system('mkdir iter0')
+ os.system('cp ' + modelfiles + ' iter0')
+ os.chdir('iter0')
+ os.system('ln -s ../land_ice_model.exe .')
+ print 'Working in ' + os.getcwd()
+ # Make an initial guess for beta
+ fi=NetCDFFile(infile,'r+')
+ fi.variables['betaTimeSeries'][:]=1.0e8
+ fi.close()
+
+ else:
+ # Setup the directory for the first run
+ os.system('mkdir iter'+str(i))
+ os.system('cp ' + modelfiles + ' iter'+str(i) )
+ os.chdir('iter'+str(i) )
+ os.system('ln -s ../land_ice_model.exe .')
+ print 'Working in ' + os.getcwd()
+
+ ########## Get Taub ###############
+ # Read the previous input, output files
+ fo_old=NetCDFFile('../iter'+str(i-1)+'/'+outfile,'r')
+ fi_old=NetCDFFile('../iter'+str(i-1)+'/'+infile,'r')
+ # Get ub and beta
+ beta_old=fi_old.variables['betaTimeSeries'][:,0] # just use the first time level from the input file since beta isn't in the output
+ #cellMask=fi.variables['cellMask'][0,:] # need to treat floating ice separately? i guess not since it gets handled independently...
+ # use the reconstructed velocities because beta is on cell centers
+ ux=fo_old.variables['uReconstructX'][0,:,:] # use the first time level
+ uy=fo_old.variables['uReconstructY'][0,:,:] # use the first time level
+ ub_old = (ux[:,-1]**2 + uy[:,-1]**2) ** 0.5 # sliding velocity - bottom vertical layer
+ ud_old = (ux[:,0]**2 + uy[:,-0]**2) ** 0.5 - ub_old # deformation velocity = surface velocity - sliding velocity
+ # Calculate previous model's Taub=beta*ub since we aren't outputting Taub
+ Taub = beta_old * ub_old
+
+ # calculate what the sliding velocity target is now that we have deformation velocity
+ fv = NetCDFFile('../' + vobsfile,'r')
+ vsurfOb = fv.variables[vobsvariable][0,:]
+ ubTarget = vsurfOb - ud_old # target sliding = observed surface velocity - modeled deformation
+ # check 0 and negative targets to something small so we always have some sliding
+ ubTarget[ubTarget<0.1] = 0.1
+
+ # Calculate new beta guess: beta = Taub_M / ub_T M=modeled, T=target
+ betaguess = Taub / ubTarget
+ betaguess[betaguess>1.0e8] = 1.0e8 # Cap how big beta can be
+ fi=NetCDFFile(infile,'r+')
+ fi.variables['betaTimeSeries'][:,0] = betaguess[:]
+
+ # use the old velo soln as the guess for the next run
+ fi.variables['normalVelocity'][:] = fo_old.variables['normalVelocity'][:] # TODO make sure this works properly
+
+ fi.close()
+ fi_old.close()
+ fo_old.close()
+ fv.close()
+
+
+ ########### Run the model - diagnostic solve
+ print 'Starting model run'
+ try:
+ subprocess.check_call(runcmd, shell=True)
+ except:
+ sys.exit('model run failed. exiting...')
+
+ ###### check the results
+ fo=NetCDFFile(outfile, 'r')
+ thickness=fo.variables['thickness'][0,:]
+ ux=fo.variables['uReconstructX'][0,:,:] # use the first time level
+ uy=fo.variables['uReconstructY'][0,:,:] # use the first time level
+ us = (ux[:,0]**2 + uy[:,0]**2) ** 0.5 # surface velocity - top vertical layer
+
+ fv = NetCDFFile('../' + vobsfile,'r')
+ vsurfOb = fv.variables[vobsvariable][0,:]
+
+ rmse = ( ((us[thickness>0.0]-vsurfOb[thickness>0.0])**2).sum() / len(us[thickness>0.0]) )**0.5 # TODO check if this is calculating right
+ print "RMSE="+str(rmse)
+ fo.close()
+ fv.close()
+
+ os.chdir('..')
+
+
+print 'Completed all iterations.'
+
+
+
Property changes on: branches/land_ice_projects/grid_tools/calibrate_beta.py
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
</font>
</pre>