[cig-commits] r19579 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d
brad at geodynamics.org
brad at geodynamics.org
Thu Feb 2 10:54:56 PST 2012
Author: brad
Date: 2012-02-02 10:54:55 -0800 (Thu, 02 Feb 2012)
New Revision: 19579
Added:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_faultinfo.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_050m.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_100m.cfg
Modified:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/pylithapp.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_offfaultdata.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_onfaultdata.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tractions.spatialdb
Log:
Updated processing scripts for HDF5 output.
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_faultinfo.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_faultinfo.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_faultinfo.py 2012-02-02 18:54:55 UTC (rev 19579)
@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# ======================================================================
+#
+
+sim = "tpv10"
+cell = "tri3"
+dx = 200
+
+inputRoot = "output/%s_%s_%03dm-fault_info" % (sim,cell,dx)
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import pylab
+
+# ----------------------------------------------------------------------
+h5 = tables.openFile("%s.h5" % inputRoot, 'r')
+vertices = h5.root.geometry.vertices[:]
+traction = h5.root.vertex_fields.initial_traction[:].squeeze()
+h5.close()
+
+from math import pi, sin
+dipDist = vertices[:,1]/sin(60.0*pi/180.0)
+
+# Expected tractions
+mask = numpy.bitwise_and(dipDist <= -10.5e+3, dipDist >= -13.5e+3)
+if sim == "tpv10":
+ shearE = mask*(0.2e+6 + ((0.760+0.0057)*7378*-dipDist)) + \
+ ~mask*0.55*7378*-dipDist
+ normalE = 7378*dipDist
+
+elif sim == "tpv11":
+ shearE = mask*(0.2e+6 + ((0.570+0.0057)*7378*-dipDist)) + \
+ ~mask*0.55*7378*-dipDist
+ normalE = 7378*dipDist
+
+pylab.plot(shearE, dipDist, 'ko',
+ traction[:,0], dipDist, 'rx',
+ normalE, dipDist, 'ko',
+ traction[:,1], dipDist, 'rx')
+pylab.show()
+
+# End of file
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_faultinfo.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/pylithapp.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/pylithapp.cfg 2012-02-02 18:23:53 UTC (rev 19578)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/pylithapp.cfg 2012-02-02 18:54:55 UTC (rev 19579)
@@ -49,8 +49,7 @@
normalizer.shear_wave_speed = 3333*m/s
[pylithapp.timedependent.formulation.time_step]
-#total_time = 15.01*s
-total_time = 0.0*s
+total_time = 15.01*s
dt = 0.01*s
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_offfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_offfaultdata.py 2012-02-02 18:23:53 UTC (rev 19578)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_offfaultdata.py 2012-02-02 18:54:55 UTC (rev 19579)
@@ -5,81 +5,70 @@
# Brad T. Aagaard
# U.S. Geological Survey
#
-# <LicenseText>
-#
# ----------------------------------------------------------------------
#
+sim = "tpv10"
cell = "tri3"
dx = 100
-dt = 0.05
-outputRoot = "output/%s_%3dm_%s" % (cell,dx,"refine")
-outdir = "scecfiles/%s_%3dm_%s/" % (cell,dx,"refine")
+inputRoot = "output/%s_%s_%3dm" % (sim,cell,dx)
+outdir = "scecfiles/%s_%s_%3dm/" % (sim,cell,dx)
# ----------------------------------------------------------------------
-
+import tables
import numpy
import time
-from pylith.utils.VTKDataReader import VTKDataReader
-
# ----------------------------------------------------------------------
-timestamps = numpy.arange(50,15001,50)
-targets = numpy.array([[-3000.0, -0.0, 0.0],
- [-2000.0, -0.0, 0.0],
- [-1000.0, -0.0, 0.0],
- [+1000.0, -0.0, 0.0],
- [+2000.0, -0.0, 0.0],
- [+3000.0, -0.0, 0.0],
- [-1000.0,-300.0, 0.0],
- [ -500.0,-300.0, 0.0],
- [ +500.0,-300.0, 0.0],
- [+1000.0,-300.0, 0.0]])
+targets = numpy.array([[-3000.0, -0.0],
+ [-2000.0, -0.0],
+ [-1000.0, -0.0],
+ [+1000.0, -0.0],
+ [+2000.0, -0.0],
+ [+3000.0, -0.0],
+ [-1000.0, -300.0],
+ [ -500.0, -300.0],
+ [ +500.0, -300.0],
+ [+1000.0, -300.0]])
-reader = VTKDataReader()
tolerance = 1.0e-6
-# Get vertices and find indices of target locations
-filename = "%s_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:
dist = ( (vertices[:,0]-target[0])**2 +
- (vertices[:,1]-target[1])**2 +
- (vertices[:,2]-target[2])**2 )**0.5
+ (vertices[:,1]-target[1])**2 )**0.5
min = numpy.min(dist)
- indices.append(numpy.where(dist <= min+tolerance)[0])
+ 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]
-disp = numpy.zeros((nsteps,ntargets,3)) # 3-D array (time, targets, components)
-vel = numpy.zeros((nsteps,ntargets,3))
-itime = 0
-for timestamp in timestamps:
- filename = "%s_t%05d.vtk" % (outputRoot,timestamp)
- data = reader.read(filename)
- fields = data['vertex_fields']
- disp[itime,0:ntargets,:] = fields['displacement'][indices,:].squeeze()
- vel[itime,0:ntargets,:] = fields['velocity'][indices,:].squeeze()
- itime += 1
+# Get datasets
+disp = h5.root.vertex_fields.displacement[:]
+vel = h5.root.vertex_fields.velocity[:]
+timeStamps = h5.root.time[:].ravel()
+nsteps = timeStamps.shape[0]
+dt = timeStamps[1] - timeStamps[0]
+h5.close()
+# Extract locations
+disp = disp[:,indices,:]
+vel = vel[:,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 7 columns of E14.6:\n" + \
"# Column #1 = time (s)\n" + \
@@ -99,11 +88,10 @@
locName = "%+04dst%03ddp%03d"
lengthScale = 1000.0
-timeScale = 1000.0
+zero = numpy.zeros( (nsteps,), dtype=numpy.float64)
dip = -targets[:,1] / lengthScale
-strike = targets[:,2] / lengthScale
+strike = 0.0*dip
body = -targets[:,0] / lengthScale
-time = timestamps / timeScale
for iloc in xrange(ntargets):
pt = locName % (round(10*body[iloc]),
@@ -118,9 +106,9 @@
strike[iloc],
dip[iloc]))
fout.write(headerB)
- data = numpy.transpose((time,
- disp[:,iloc,2],
- vel[:,iloc,2],
+ data = numpy.transpose((timeStamps,
+ zero,
+ zero,
-disp[:,iloc,1],
-vel[:,iloc,1],
-disp[:,iloc,0],
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_onfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_onfaultdata.py 2012-02-02 18:23:53 UTC (rev 19578)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tabulate_onfaultdata.py 2012-02-02 18:54:55 UTC (rev 19579)
@@ -5,78 +5,69 @@
# Brad T. Aagaard
# U.S. Geological Survey
#
-# <LicenseText>
-#
# ----------------------------------------------------------------------
#
-cell = "quad4"
+sim = "tpv10"
+cell = "tri3"
dx = 100
-dt = 0.05
-outputRoot = "output/%s_%3dm_%s" % (cell,dx,"refine")
-outdir = "scecfiles/%s_%3dm_%s/" % (cell,dx,"refine")
+inputRoot = "output/%s_%s_%3dm-fault" % (sim,cell,dx)
+outdir = "scecfiles/%s_%s_%3dm/" % (sim,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],
- [ -750.0000000, -1299.0381057, 0.0000000],
- [ -1500.0000000, -2598.0762114, 0.0000000],
- [ -2250.0000000, -3897.1143170, 0.0000000],
- [ -3750.0000000, -6495.1905284, 0.0000000],
- [ -6000.0000000, -10392.3048454, 0.0000000]])
+targets = numpy.array([[ 0.0, 0.0, ],
+ [ -750.0, -1299.0381057],
+ [-1500.0, -2598.0762114],
+ [-2250.0, -3897.1143170],
+ [-3750.0, -6495.1905284],
+ [-6000.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:
dist = ( (vertices[:,0]-target[0])**2 +
- (vertices[:,1]-target[1])**2 +
- (vertices[:,2]-target[2])**2 )**0.5
+ (vertices[:,1]-target[1])**2 )**0.5
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 = TPV205\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" + \
"# Column #1 = time (s)\n" + \
@@ -96,12 +87,11 @@
locName = "st%03ddp%03d"
lengthScale = 1000.0
-timeScale = 1000.0
+zero = numpy.zeros( (nsteps,), dtype=numpy.float64)
dip = -2*targets[:,0] / lengthScale
-strike = targets[:,2] / lengthScale
-time = timestamps / timeScale
-print "time", time
+strike = 0.0*dip
+
for iloc in xrange(ntargets):
pt = locName % (round(10*strike[iloc]),
round(10*dip[iloc]))
@@ -112,10 +102,10 @@
fout.write("# num_timesteps = %8d\n" % nsteps)
fout.write(locHeader % (strike[iloc], dip[iloc]))
fout.write(headerB)
- data = numpy.transpose((time,
- slip[:,iloc,2],
- slip_rate[:,iloc,2],
- traction[:,iloc,2]/1e+6,
+ data = numpy.transpose((timeStamps,
+ zero,
+ zero,
+ zero,
slip[:,iloc,0],
slip_rate[:,iloc,0],
traction[:,iloc,0]/1e+6,
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg 2012-02-02 18:23:53 UTC (rev 19578)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg 2012-02-02 18:54:55 UTC (rev 19579)
@@ -4,6 +4,7 @@
# ----------------------------------------------------------------------
# faults
# ----------------------------------------------------------------------
+[pylithapp.timedependent.interfaces.fault]
db_initial_tractions.iohandler.filename = tpv10_tractions.spatialdb
friction.db_properties.iohandler.filename = tpv10_friction.spatialdb
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tractions.spatialdb
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tractions.spatialdb 2012-02-02 18:23:53 UTC (rev 19578)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tractions.spatialdb 2012-02-02 18:54:55 UTC (rev 19579)
@@ -12,7 +12,7 @@
}
}
0.0 0.0 0.0 0.0
--5.2499 0.0 42.60876158 -77.4675244
+-5.2499 0.0 42.60713842 -77.4675244
-5.2500 0.0 59.51801330 -77.4690000
-6.7500 0.0 76.46601710 -99.6030000
-6.7501 0.0 54.78246158 -99.6044756
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_050m.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_050m.cfg (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_050m.cfg 2012-02-02 18:54:55 UTC (rev 19579)
@@ -0,0 +1,29 @@
+# -*- Python -*-
+[pylithapp]
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[pylithapp.mesh_generator.reader]
+filename = tri3_050m.exo
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.timedependent.normalizer]
+wave_period = 0.15*s
+
+[pylithapp.timedependent.formulation.time_step]
+dt = 0.005*s
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem.formulation.output.output]
+writer.filename = output/tpv10_tri3_050m.h5
+
+[pylithapp.timedependent.interfaces.fault.output]
+writer.filename = output/tpv10_tri3_050m-fault.h5
+
+[pylithapp.timedependent.materials.elastic.output]
+writer.filename = output/tpv10_tri3_050m-elastic.h5
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_100m.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_100m.cfg (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_100m.cfg 2012-02-02 18:54:55 UTC (rev 19579)
@@ -0,0 +1,29 @@
+# -*- Python -*-
+[pylithapp]
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[pylithapp.mesh_generator.reader]
+filename = tri3_100m.exo
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.timedependent.normalizer]
+wave_period = 0.3*s
+
+[pylithapp.timedependent.formulation.time_step]
+dt = 0.01*s
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem.formulation.output.output]
+writer.filename = output/tpv10_tri3_100m.h5
+
+[pylithapp.timedependent.interfaces.fault.output]
+writer.filename = output/tpv10_tri3_100m-fault.h5
+
+[pylithapp.timedependent.materials.elastic.output]
+writer.filename = output/tpv10_tri3_100m-elastic.h5
More information about the CIG-COMMITS
mailing list