[cig-commits] r19836 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 21 10:25:02 PDT 2012
Author: brad
Date: 2012-03-21 10:25:02 -0700 (Wed, 21 Mar 2012)
New Revision: 19836
Removed:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/output/
Modified:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_onfaultdata.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_ruptime.py
Log:
Remove output dir from repository.
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_onfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_onfaultdata.py 2012-03-21 17:24:24 UTC (rev 19835)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_onfaultdata.py 2012-03-21 17:25:02 UTC (rev 19836)
@@ -10,43 +10,35 @@
# ----------------------------------------------------------------------
#
-cell = "hex8"
+sim = "tpv10"
+cell = "tet4"
dx = 200
-dt = 0.05
-# outputRoot = "output/%s_%3dm_%s" % (cell,dx,"refine")
-# outdir = "scecfiles/%s_%3dm_%s/" % (cell,dx,"refine")
+inputRoot = "output/%s_%s_%03dm-fault" % (sim,cell,dx)
+outdir = "scecfiles/%s_%s_%03dm/" % (sim,cell,dx)
-outputRoot = "output/%s_%3dm" % (cell,dx)
-outdir = "scecfiles/%s_%3dm/" % (cell,dx)
-
+# ----------------------------------------------------------------------
+import tables
import numpy
import time
-from pylith.utils.VTKDataReader import VTKDataReader
-
# ----------------------------------------------------------------------
-timestamps = numpy.arange(50,15001,50)
-targets = numpy.array([[ 0.0000000, 0.0000000, 0.0000000],
- [ 0.0000000, 4500.0000000, 0.0000000],
- [ 0.0000000, 12000.0000000, 0.0000000],
- [ -750.0000000, 0.0000000, -1299.0381057],
- [ -1500.0000000, 0.0000000, -2598.0762114],
- [ -2250.0000000, 0.0000000, -3897.1143170],
- [ -3750.0000000, 0.0000000, -6495.1905284],
- [ -3750.0000000, 4500.0000000, -6495.1905284],
- [ -3750.0000000, 12000.0000000, -6495.1905284],
- [ -6000.0000000, 0.0000000, -10392.3048454]])
+targets = numpy.array([[ 0.0, 0.0, 0.0 ],
+ [ 0.0, 4500.0, 0.0 ],
+ [ 0.0, 12000.0, 0.0 ],
+ [ -750.0, 0.0, -1299.0381057],
+ [ -1500.0, 0.0, -2598.0762114],
+ [ -2250.0, 0.0, -3897.1143170],
+ [ -3750.0, 0.0, -6495.1905284],
+ [ -3750.0, 4500.0, -6495.1905284],
+ [ -3750.0, 12000.0, -6495.1905284],
+ [ -6000.0, 0.0, -10392.3048454]])
-reader = VTKDataReader()
tolerance = 1.0e-6
-# Get vertices and find indices of target locations
-filename = "%s-fault_t%05d.vtk" % (outputRoot,timestamps[0])
-data = reader.read(filename)
-
-vertices = numpy.array(data['vertices'])
+h5 = tables.openFile("%s.h5" % inputRoot, 'r')
+vertices = h5.root.geometry.vertices[:]
ntargets = targets.shape[0]
indices = []
for target in targets:
@@ -56,36 +48,35 @@
min = numpy.min(dist)
indices.append(numpy.where(dist <= min+tolerance)[0][0])
-print "Indices", indices
-print "Coordinates of selected points:",vertices[indices,:]
+print "Indices: ", indices
+print "Coordinates of selected points:\n",vertices[indices,:]
-# Extract values
-nsteps = timestamps.shape[0]
-slip = numpy.zeros((nsteps,ntargets,3)) # 3-D array (time, targets, components)
-slip_rate = numpy.zeros((nsteps,ntargets,3))
-traction = numpy.zeros((nsteps,ntargets,3))
-itime = 0
-for timestamp in timestamps:
- filename = "%s-fault_t%05d.vtk" % (outputRoot,timestamp)
-# print "filename", filename
- data = reader.read(filename)
- fields = data['vertex_fields']
- slip[itime,0:ntargets,:] = fields['slip'][indices,:].squeeze()
- slip_rate[itime,0:ntargets,:] = fields['slip_rate'][indices,:].squeeze()
- traction[itime,0:ntargets,:] = fields['traction'][indices,:].squeeze()
- itime += 1
+# Get datasets
+slip = h5.root.vertex_fields.slip[:]
+slip_rate = h5.root.vertex_fields.slip_rate[:]
+traction = h5.root.vertex_fields.traction[:]
+timeStamps = h5.root.time[:].ravel()
+nsteps = timeStamps.shape[0]
+dt = timeStamps[1] - timeStamps[0]
+h5.close()
+
+# Extract locations
+slip = slip[:,indices,:]
+slip_rate = slip_rate[:,indices,:]
+traction = traction[:,indices,:]
+
# Write data
headerA = \
- "# problem = TPV210\n" + \
- "# author = Surendra N. Somala\n" + \
+ "# problem = %s\n" % sim.upper() + \
+ "# author = Brad Aagaard\n" + \
"# date = %s\n" % (time.asctime()) + \
"# code = PyLith\n" + \
- "# code_version = 1.5.0a\n" + \
- "# element_size = %s\n" % dx
+ "# code_version = 1.7.0a (scecdynrup branch)\n" + \
+ "# element_size = %s m\n" % dx
headerB = \
- "# Time series in 8 columns of E14.6:\n" + \
+ "# Time series in 7 columns of E14.6:\n" + \
"# Column #1 = time (s)\n" + \
"# Column #2 = horizontal right-lateral slip (m)\n" + \
"# Column #3 = horizontal right-lateral slip rate (m/s)\n" + \
@@ -96,19 +87,17 @@
"# Column #8 = normal stress (MPa)\n" + \
"#\n" + \
"# Data fields\n" + \
- "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress \n" + \
+ "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress\n" + \
"#\n"
locHeader = "# location = on fault, %3.1f km along strike and %3.1f km depth\n"
-locName = "st%03ddp%03d"
+locName = "st%+04ddp%03d"
-lengthScale = 1000.0
-timeScale = 1000.0
+lengthScale = 1.0e+3
dip = -2*targets[:,0] / lengthScale
strike = targets[:,1] / lengthScale
-time = timestamps / timeScale
-print "time", time
+
for iloc in xrange(ntargets):
pt = locName % (round(10*strike[iloc]),
round(10*dip[iloc]))
@@ -119,13 +108,13 @@
fout.write("# num_timesteps = %8d\n" % nsteps)
fout.write(locHeader % (strike[iloc], dip[iloc]))
fout.write(headerB)
- data = numpy.transpose((time,
+ data = numpy.transpose((timeStamps,
-slip[:,iloc,0],
-slip_rate[:,iloc,0],
-traction[:,iloc,0]/1e+6,
- slip[:,iloc,1],
- slip_rate[:,iloc,1],
- traction[:,iloc,1]/1e+6,
+ -slip[:,iloc,1],
+ -slip_rate[:,iloc,1],
+ -traction[:,iloc,1]/1e+6,
traction[:,iloc,2]/1e+6))
numpy.savetxt(fout, data, fmt='%14.6e')
fout.close()
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_ruptime.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_ruptime.py 2012-03-21 17:24:24 UTC (rev 19835)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210/tabulate_ruptime.py 2012-03-21 17:25:02 UTC (rev 19836)
@@ -5,93 +5,67 @@
# Brad T. Aagaard
# U.S. Geological Survey
#
-# <LicenseText>
-#
# ----------------------------------------------------------------------
#
-cell = "hex8"
+sim = "tpv10"
+cell = "tet4"
dx = 200
-dt = 0.05
-outputRoot = "output/%s_%3dm" % (cell,dx)
-outdir = "scecfiles/%s_%3dm/" % (cell,dx)
+inputRoot = "output/%s_%s_%03dm-fault" % (sim,cell,dx)
+outdir = "scecfiles/%s_%s_%03dm/" % (sim,cell,dx)
+# ----------------------------------------------------------------------
+import tables
import numpy
import time
-from pylith.utils.VTKDataReader import VTKDataReader
-
# ----------------------------------------------------------------------
-timestamps = numpy.arange(50,12001,50)
-
-
-reader = VTKDataReader()
threshold = 0.001 # threshold for detecting slip has started
+maxTime = 1.0e+10 # Large value for default rupture time
-filename = "%s-fault_t%05d.vtk" % (outputRoot,timestamps[0])
-data = reader.read(filename)
-vertices = numpy.array(data['vertices'])
-npts = vertices.shape[0]
+h5 = tables.openFile("%s.h5" % inputRoot, 'r')
+vertices = h5.root.geometry.vertices[:]
+slip = h5.root.vertex_fields.slip[:]
+slipRate = h5.root.vertex_fields.slip_rate[:]
+timeStamps = h5.root.time[:].ravel()
+dt = timeStamps[1] - timeStamps[0]
-# Extract values
-nsteps = timestamps.shape[0]
+h5.close()
-# Set default rupture time to a large value (1.0e+9)
-rupTime = 1.0e+9 * numpy.ones( (npts,), dtype=numpy.float64)
-#rupTime = numpy.zeros( (npts,), dtype=numpy.float64)
+nsteps = timeStamps.shape[0]
+npts = vertices.shape[0]
+rupTime = maxTime * numpy.ones( (npts,), dtype=numpy.float64)
+
# Create buffer for current rupture time
tmpTime = numpy.zeros( (npts,), dtype=numpy.float64)
itime = 0
-for timestamp in timestamps:
- t = (itime+1)*dt # itime=0, t=dt
+for timestamp in timeStamps:
+ t = timeStamps[itime]
- # Get slip rate field
- filename = "%s-fault_t%05d.vtk" % (outputRoot,timestamp)
- data = reader.read(filename)
- fields = data['vertex_fields']
- slipRate = fields['slip'][:,:].squeeze()
-
# Compute magnitude of slip rate
- slipRateMag = (slipRate[:,0]**2 + slipRate[:,1]**2)**0.5
+ slipRateMag = (slipRate[itime,:,0]**2 + slipRate[itime,:,1]**2)**0.5
# Set rupture time at locations where threshold is exceeded
- tmpTime[slipRateMag > threshold] = t
+ mask = slipRateMag > threshold
+ tmpTime[:] = maxTime
+ tmpTime[mask] = t
- #print "slipRateMag \n", slipRateMag
- #print "slipRateMag > threshold \n ", slipRateMag > threshold
-
- #print "tmpTime \n", tmpTime[1:50]
- #print "time", t
-
- # Get indicates where current time is less than current rupture
- # time (this is only the locations that just started slipping)
- #indices = numpy.where((tmpTime < rupTime) and (rupTime < 0.00001))[0]
- #indices = ((tmpTime < rupTime) & (rupTime < testarr))
- indices = numpy.where((tmpTime < rupTime) & (tmpTime > 0.00001))[0]
-
- #print "indices \n", indices[1:200]
-
-
+ indices = numpy.where(tmpTime < rupTime)[0]
rupTime[indices] = t
- #print "rupTime \n", rupTime[1:200]
-
itime += 1
-# print "rupTime \n", rupTime
-# rupTime[rupTime < 0.00001] = 1.0e+9
-
# Write data
headerA = \
- "# problem = TPV210\n" + \
- "# author = Surendra N. Somala\n" + \
+ "# problem = %s\n" % sim.upper() + \
+ "# author = Brad Aagaard\n" + \
"# date = %s\n" % (time.asctime()) + \
"# code = PyLith\n" + \
- "# code_version = 1.5.0a\n" + \
- "# element_size = %s\n" % dx + \
+ "# code_version = 1.7.0a (scecdynrup branch)\n" + \
+ "# element_size = %s m\n" % dx + \
"# Contour data in 3 columns of E14.6:\n" + \
"# Column #1 = Distance along strike from hypocenter (m)\n" + \
"# Column #2 = Distance down-dip from surface (m)\n" + \
@@ -102,12 +76,8 @@
"#\n"
-lengthScale = 1000.0
-timeScale = 1000.0
distDip = -2*vertices[:,0]
distStrike = vertices[:,1]
-time = timestamps / timeScale
-print "time", time
filename = "%sruptime.dat" % (outdir)
fout = open(filename, "w")
More information about the CIG-COMMITS
mailing list