[cig-commits] r16667 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d
surendra at geodynamics.org
surendra at geodynamics.org
Fri May 7 18:18:50 PDT 2010
Author: surendra
Date: 2010-05-07 18:18:50 -0700 (Fri, 07 May 2010)
New Revision: 16667
Modified:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py
Log:
Clean up of python scripts to extract data from VTK files for tpv205-2d benchmark
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py 2010-05-07 20:52:54 UTC (rev 16666)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py 2010-05-08 01:18:50 UTC (rev 16667)
@@ -24,7 +24,7 @@
from pylith.utils.VTKDataReader import VTKDataReader
# ----------------------------------------------------------------------
-timestamps = numpy.arange(0,12000,50)
+timestamps = numpy.arange(50,12001,50)
targets = numpy.array([[3000.0, -12000.0, 0.0],
[3000.0, +12000.0, 0.0]])
@@ -66,7 +66,7 @@
# Write data
headerA = \
"# problem = TPV205\n" + \
- "# author = SurendraSomala\n" + \
+ "# author = Surendra N. Somala\n" + \
"# date = %s\n" % (time.asctime()) + \
"# code = PyLith\n" + \
"# code_version = 1.5.0\n" + \
@@ -90,7 +90,7 @@
locName = "%+04dst%+04ddp%03d"
lengthScale = 1000.0
-timeScale = 100.0
+timeScale = 1000.0
dip = 7.5
body = targets[:,0] / lengthScale
strike = targets[:,1] / lengthScale
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py 2010-05-07 20:52:54 UTC (rev 16666)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py 2010-05-08 01:18:50 UTC (rev 16667)
@@ -10,86 +10,75 @@
# ----------------------------------------------------------------------
#
-outputRoot = "output/quad4_100m"
-outdir = "scecfiles/quad4_100m/"
+cell = "quad4"
dx = 100
+dt = 0.05
+outputRoot = "output_pacha/%s_%3dm" % (cell,dx)
+outdir = "scecfiles/%s_%3dm/" % (cell,dx)
+
import numpy
import time
from pylith.utils.VTKDataReader import VTKDataReader
-from numpy import *
# ----------------------------------------------------------------------
-def feq(a,b):
- if numpy.abs(a-b)<0.0001:
- return 1
- else:
- return 0
+timestamps = numpy.arange(50,12001,50)
+if dx == 200:
+ targets = numpy.array([[0.0, -12000.0, 0.0],
+ [0.0, -7600.0, 0.0],
+ [0.0, -4400.0, 0.0],
+ [0.0, 0.0, 0.0],
+ [0.0, +4400.0, 0.0],
+ [0.0, +7600.0, 0.0],
+ [0.0, +12000.0, 0.0]])
+elif dx == 100:
+ targets = numpy.array([[0.0, -12000.0, 0.0],
+ [0.0, -7500.0, 0.0],
+ [0.0, -4500.0, 0.0],
+ [0.0, 0.0, 0.0],
+ [0.0, +4500.0, 0.0],
+ [0.0, +7500.0, 0.0],
+ [0.0, +12000.0, 0.0]])
+
-# ----------------------------------------------------------------------
-# Extract On Fault stations data ( Slip, Slip Rate, Tractions )
-# ----------------------------------------------------------------------
-nfileNames = numpy.arange(1,1202,5)
-filename = "%s-fault_t%04d.vtk" % (outputRoot,nfileNames[0])
-
-dt = 0.05
-nsteps = len(nfileNames)
-print nsteps
-
-fieldNames = ["slip","slip_rate","traction"]
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)
-values = numpy.array(data['vertices'])
-(nvertices, spaceDim) = values.shape
-x_OnFltStat = numpy.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
-#y_OnFltStat = numpy.array([-12000.0, -7600.0, -4400.0, 0.0, 4400.0, 7600.0, 12000.0]) # CHECK THESE VALUES
-y_OnFltStat = numpy.array([-12000.0, -7500.0, -4500.0, 0.0, 4500.0, 7500.0, 12000.0]) # CHECK THESE VALUES
-Index_OnFltStat = numpy.zeros(7)
+vertices = numpy.array(data['vertices'])
+ntargets = targets.shape[0]
+indices = []
+for target in targets:
+ dist = ( (vertices[:,0]-target[0])**2 +
+ (vertices[:,1]-target[1])**2 +
+ (vertices[:,2]-target[2])**2 )**0.5
+ min = numpy.min(dist)
+ indices.append(numpy.where(dist <= min+tolerance)[0])
-# Find the indices of On Fault Stations
-nFltStat = -1
-for i in xrange(nvertices):
- if feq(values[i,1],y_OnFltStat[0]) or \
- feq(values[i,1],y_OnFltStat[1]) or \
- feq(values[i,1],y_OnFltStat[2]) or \
- feq(values[i,1],y_OnFltStat[3]) or \
- feq(values[i,1],y_OnFltStat[4]) or \
- feq(values[i,1],y_OnFltStat[5]) or \
- feq(values[i,1],y_OnFltStat[6]):
- nFltStat += 1
- Index_OnFltStat[nFltStat] = i
+print "Indices", indices
+print "Coordinates of selected points:",vertices[indices,:]
-Index_OnFltStat = numpy.int_(Index_OnFltStat)
-print "On-Fault indices", Index_OnFltStat
-
# Extract values
-slip = numpy.zeros((nsteps,7,3)) # 3-D array (ntime, nstats, ncomponents)
-slip_rate = numpy.zeros((nsteps,7,3))
-traction = numpy.zeros((nsteps,7,3))
-nFile = -1
-for i in nfileNames:
- nFile += 1
- filename = "%s-fault_t%04d.vtk" % (outputRoot,i)
+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)
data = reader.read(filename)
- nField = 0
- for name in fieldNames:
- nField += 1
- values = data['vertex_fields'][name]
- if nField == 1:
- slip[nFile,:,:] = values[Index_OnFltStat,:]
- elif nField == 2:
- slip_rate[nFile,:,:] = values[Index_OnFltStat,:]
- elif nField == 3:
- traction[nFile,:,:] = values[Index_OnFltStat,:]
+ 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
-
-# ----------------------------------------------------------------------
-# Write On Fault stations data to file ( Slip, Slip Rate, Tractions )
-# ----------------------------------------------------------------------
-
+# Write data
headerA = \
"# problem = TPV205\n" + \
"# author = Surendra N. Somala\n" + \
@@ -115,15 +104,13 @@
locName = "st%+04ddp%03d"
lengthScale = 1000.0
-timeScale = 100.0
+timeScale = 1000.0
dip = 7.5
-strike = y_OnFltStat[:] / lengthScale
-time = nfileNames[:] / timeScale
+strike = targets[:,1] / lengthScale
+time = timestamps / timeScale
print "time", time
-# Write data
-nlocs = slip.shape[1]
-for iloc in xrange(nlocs):
+for iloc in xrange(ntargets):
pt = locName % (round(10*strike[iloc]),
round(10*dip))
filename = "%sfault%s.dat" % (outdir,pt)
More information about the CIG-COMMITS
mailing list