<p><b>mhoffman@lanl.gov</b> 2013-03-07 15:25:14 -0700 (Thu, 07 Mar 2013)</p><p>BRANCH COMMIT - land ice<br>
<br>
Minor updates to a few scripts and the addition of a script that visualizes the layout of blocks in a domain.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice_projects/grid_tools/calibrate_beta.py
===================================================================
--- branches/land_ice_projects/grid_tools/calibrate_beta.py 2013-03-07 21:43:23 UTC (rev 2560)
+++ branches/land_ice_projects/grid_tools/calibrate_beta.py 2013-03-07 22:25:14 UTC (rev 2561)
@@ -1,5 +1,6 @@
#!/usr/bin/python
# Script for calibrating beta using MPAS
+# Note that the script currenly assumes that the target velocity is an observed surface velocity, as opposed to a balance velocity - but that distinction may not matter too much.
# Matt Hoffman, March 2013
import os, sys
@@ -23,9 +24,9 @@
################ 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
+vobsfile='observedSpeed.nc' # This is the file that has the observed surface velocity field - it can have anything else, too.
+vobsvariable='observedSpeed' # this is the name of the observed velocity field variable in the above file. It should have dimensions (Time, nCells)
+numiters = 5 # 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.
Modified: branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py
===================================================================
--- branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py 2013-03-07 21:43:23 UTC (rev 2560)
+++ branches/land_ice_projects/grid_tools/copy_CISM_grid_to_MPAS_grid.py 2013-03-07 22:25:14 UTC (rev 2561)
@@ -130,10 +130,12 @@
except:
print 'Output file is missing the dimension nVertLevels. Might not be a problem.'
+ try:
+ # 1d vertical fields
+ layerThicknessFractions = outfile.variables['layerThicknessFractions'][:]
+ except:
+ print 'Output file is missing the variable layerThicknessFractions. Might not be a problem.'
- # 1d vertical fields
- layerThicknessFractions = outfile.variables['layerThicknessFractions'][:]
-
# '2d' spatial fields on cell centers
xCell = outfile.variables['xCell'][:]
print 'xCell min/max:', xCell.min(), xCell.max()
@@ -273,6 +275,25 @@
+try:
+ balvel = infile.variables['balvel']
+ try:
+ sf = balvel.scale_factor
+ except:
+ sf = 1.0
+ if netCDF_module == 'netCDF4':
+ sf = 1.0 # netCDF4 should autoscale!
+ observedSpeed = outfile.variables['observedSpeed']
+ print '</font>
<font color="blue">balvel min/max', balvel[:].min()*sf, balvel[:].max()*sf
+ observedSpeed[timelevout,:] = BilinearInterp(x0, y0, balvel[timelev,:,:]*sf, xCell, yCell)
+ print 'new observedSpeed min/max', observedSpeed[:].min(), observedSpeed[:].max()
+ del balvel, observedSpeed
+except:
+ print '</font>
<font color="black">problem with balvel field (e.g. not found in input file), skipping...</font>
<font color="gray">'
+
+
+
+
# Variables on edges...?
#normalVelocity = outfile.variables['normalVelocity'][:]
Modified: branches/land_ice_projects/grid_tools/plot_mpas_field.py
===================================================================
--- branches/land_ice_projects/grid_tools/plot_mpas_field.py 2013-03-07 21:43:23 UTC (rev 2560)
+++ branches/land_ice_projects/grid_tools/plot_mpas_field.py 2013-03-07 22:25:14 UTC (rev 2561)
@@ -57,7 +57,10 @@
f = NetCDFFile(options.filename,'r')
# Get grid stuff
-xtime = f.variables['xtime'][:] # Not needed unless trying to print actual time stamp
+try:
+ xtime = f.variables['xtime'][:] # Not needed unless trying to print actual time stamp
+except:
+ xtime = numpy.zeros((2,2))
if not options.time:
print "No time provided. Using time 0, if time is needed."
@@ -69,12 +72,20 @@
else:
time_slices = int(options.time)
print "Using time level " + options.time
-
-xCell = f.variables['xCell']
-yCell = f.variables['yCell']
-xEdge = f.variables['xEdge']
-yEdge = f.variables['yEdge']
-angleEdge = f.variables['angleEdge']
+try:
+ xCell = f.variables['xCell']
+ yCell = f.variables['yCell']
+except:
+ print 'xCell and/or yCell are missing or have problems. Might not be a problem.'
+try:
+ xEdge = f.variables['xEdge']
+ yEdge = f.variables['yEdge']
+except:
+ print 'xEdge and/or yEdge are missing or have problems. Might not be a problem.'
+try:
+ angleEdge = f.variables['angleEdge']
+except:
+ print 'angleEdge is missing or has problems. Might not be a problem.'
nCells = xCell.shape[0] # This is an easier way to get nCells than querying the nCells dimension because the way to do that is module-specific
Added: branches/land_ice_projects/grid_tools/visualize_blocks.py
===================================================================
--- branches/land_ice_projects/grid_tools/visualize_blocks.py (rev 0)
+++ branches/land_ice_projects/grid_tools/visualize_blocks.py 2013-03-07 22:25:14 UTC (rev 2561)
@@ -0,0 +1,98 @@
+#!/usr/bin/python
+# Script capable of plotting most MPAS fields. Invoke with 'python plot_mpas_field.py --help for details about how to use.
+# Matt Hoffman, June 14, 2012
+
+import sys
+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')
+
+
+print "** Gathering information."
+parser = OptionParser()
+parser.add_option("-g", "--gridfile", dest="gridfile", help="grid file to visualize; default: land_ice_grid.nc", metavar="FILE")
+parser.add_option("-b", "--blockfile", dest="blockfile", help="block file to visualize; default: graph.info.part.2", metavar="FILE")
+
+
+options, args = parser.parse_args()
+
+if not options.gridfile:
+ print "No grid filename provided. Using land_ice_grid.nc."
+ options.gridfile = "land_ice_grid.nc"
+
+if not options.blockfile:
+ print "No block filename provided. Using graph.info.part.2."
+ options.blockfile = "graph.info.part.2"
+
+try:
+ f = NetCDFFile(options.gridfile,'r')
+except:
+ sys.exit('Error loading grid file.')
+
+
+# Get grid stuff
+#times = f.variables['xtime'] # Not needed unless trying to print actual time stamp
+try:
+ xCell = f.variables['xCell']
+ yCell = f.variables['yCell']
+except:
+ sys.exit('xCell and/or yCell is/are missing.')
+
+try:
+ blocks = np.genfromtxt(options.blockfile)
+except:
+ sys.exit('block file could not be read properly.')
+
+# Get the needed slice and determine the plot title. Make some assumptions about how dimensions are arranged:
+plottitle = 'processor decomposition for ' + options.blockfile
+
+
+# Determine color axis max & min values to use.
+#if not options.maximum:
+# maxval = var_slice[:].max()
+#else:
+# maxval = float(options.maximum)
+#
+#if not options.minimum:
+# minval = var_slice[:].min()
+#else:
+# minval = float(options.minimum)
+
+
+# MAKE THE PLOT
+print '** Beginning to create plot.'
+fig = plt.figure(1, facecolor='w')
+ax = fig.add_subplot(111, aspect='equal')
+plt.scatter(xCell[:], yCell[:], 60, blocks, edgecolors='none') #, vmin=minval, vmax=maxval)
+plt.colorbar()
+plt.title( plottitle )
+plt.draw()
+
+
+#if options.saveimages:
+# plotname = varname + '.' + options.filename + '.png'
+# plt.savefig(plotname)
+# print 'Saved plot as ' + plotname
+
+#if options.hidefigs:
+# print "Plot display disabled with -n argument."
+#else:
+# print 'Showing plot... Close plot window to exit.'
+# plt.show()
+
+plt.show()
+
+f.close()
+
Property changes on: branches/land_ice_projects/grid_tools/visualize_blocks.py
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
</font>
</pre>