[cig-commits] r19401 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16
brad at geodynamics.org
brad at geodynamics.org
Fri Jan 20 12:02:45 PST 2012
Author: brad
Date: 2012-01-20 12:02:44 -0800 (Fri, 20 Jan 2012)
New Revision: 19401
Added:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_offfaultdata.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_onfaultdata.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_ruptime.py
Modified:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/createbc.jou
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/geometry.jou
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/pylithapp.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4_150m.jou
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_075m.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_150m.cfg
Log:
Addd tabulating scripts. Fixed some typos. Small tweaks to some settings.
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/createbc.jou
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/createbc.jou 2012-01-20 19:15:44 UTC (rev 19400)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/createbc.jou 2012-01-20 20:02:44 UTC (rev 19401)
@@ -61,12 +61,12 @@
# Create nodeset for subset of +z face
# ----------------------------------------------------------------------
group "face_zpos_subset" add node in group face_zpos
-group "face_zpos_subset" remove node with y_coord > +24.001e+3
-group "face_zpos_subset" remove node with y_coord < -24.001e+3
-group "face_zpos_subset" remove node with x_coord > +10.001e+3
-group "face_zpos_subset" remove node with x_coord < -10.001e+3
+group "face_zpos_subset" remove node with y_coord > {+24.001*km}
+group "face_zpos_subset" remove node with y_coord < {-24.001*km}
+group "face_zpos_subset" remove node with x_coord > {+10.001*km}
+group "face_zpos_subset" remove node with x_coord < {-10.001*km}
nodeset 17 group face_zpos_subset
nodeset 17 name "face zpos_subset"
-
+# End of file
\ No newline at end of file
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/geometry.jou
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/geometry.jou 2012-01-20 19:15:44 UTC (rev 19400)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/geometry.jou 2012-01-20 20:02:44 UTC (rev 19401)
@@ -3,7 +3,7 @@
# ----------------------------------------------------------------------
# Block is 48.0 km x 72.0 km x 32.0 km
# -24.0 km <= x <= 24.0 km
-# -32.0 km <= y <= 32.0 km
+# -36.0 km <= y <= 36.0 km
# -32.0 km <= z <= 0.0 km
#
#{Units('si')}
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/pylithapp.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/pylithapp.cfg 2012-01-20 19:15:44 UTC (rev 19400)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/pylithapp.cfg 2012-01-20 20:02:44 UTC (rev 19401)
@@ -133,7 +133,7 @@
[pylithapp.problem.formulation.output.domain]
vertex_data_fields=[displacement,velocity]
output_freq = time_step
-time_step = 12.0*s
+time_step = 15.0*s
writer = pylith.meshio.DataWriterHDF5ExtMesh
[pylithapp.problem.formulation.output.subdomain]
@@ -154,7 +154,7 @@
cell_data_fields = []
cell_filter = pylith.meshio.CellFilterAvgMesh
output_freq = time_step
-time_step = 12.0*s
+time_step = 15.0*s
writer = pylith.meshio.DataWriterHDF5ExtMesh
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_offfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_offfaultdata.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_offfaultdata.py 2012-01-20 20:02:44 UTC (rev 19401)
@@ -0,0 +1,121 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+sim = "tpv16"
+cell = "tet4"
+dx = 150
+dt = 0.05
+
+inputRoot = "output/%s_%s_%03dm-groundsurf" % (sim,cell,dx)
+outputRoot = "scecfiles/%s_%s_%03dm/" % (sim,cell,dx)
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import time
+
+# ----------------------------------------------------------------------
+targets = numpy.array([[-6000.0, -9000.0, 0.0],
+ [ -200.0, -9000.0, 0.0],
+ [ +200.0, -9000.0, 0.0],
+ [+6000.0, -9000.0, 0.0],
+ [-6000.0, 0.0, 0.0],
+ [ -200.0, 0.0, 0.0],
+ [ +200.0, 0.0, 0.0],
+ [+6000.0, 0.0, 0.0],
+ [-6000.0, +9000.0, 0.0],
+ [ -200.0, +9000.0, 0.0],
+ [ +200.0, +9000.0, 0.0],
+ [+6000.0, +9000.0, 0.0]])
+
+
+h5 = tables.openFile("%s.h5" % (inputRoot), 'r')
+vertices = h5.root.geometry.vertices[:]
+ntargets = targets.shape[0]
+indices = []
+tolerance = 1.0e-6
+for target in targets:
+ dist = ( (vertices[:,0]-target[0])**2 +
+ (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:\n",vertices[indices,:]
+
+# 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 = %s-2D\n" % sim.upper() + \
+ "# author = Brad Aagaard\n" + \
+ "# date = %s\n" % (time.asctime()) + \
+ "# code = PyLith\n" + \
+ "# code_version = 1.7.0a (scecdynrup branch)\n" + \
+ "# element_size = %s\n" % dx
+headerB = \
+ "# Time series in 7 columns of E14.6:\n" + \
+ "# Column #1 = time (s)\n" + \
+ "# Column #2 = horizontal fault-parallel displacement (m)\n" + \
+ "# Column #3 = horizontal fault-parallel velocity (m/s)\n" + \
+ "# Column #4 = vertical displacement (m)\n" + \
+ "# Column #5 = vertical velocity (m/s)\n" + \
+ "# Column #6 = horizontal fault-normal displacement (m)\n" + \
+ "# Column #7 = horizontal fault-normal velocity (m/s)\n" + \
+ "#\n" + \
+ "# Data fields\n" + \
+ "t h-disp h-vel v-disp v-vel n-disp n-vel\n" + \
+ "#\n"
+
+locHeader = "# location = %3.1f km off fault, %3.1f km along strike " \
+ "and %3.1f km depth\n"
+locName = "%+04dst%+04ddp%03d"
+
+lengthScale = 1000.0
+timeScale = 1000.0
+dip = 0.0
+body = -targets[:,0] / lengthScale
+strike = (targets[:,1]) / lengthScale
+
+for iloc in xrange(ntargets):
+ pt = locName % (round(10*body[iloc]),
+ round(10*strike[iloc]),
+ round(10*dip))
+ filename = "%sbody%s.dat" % (outputRoot, pt)
+ fout = open(filename, 'w')
+ fout.write(headerA)
+ fout.write("# time_step = %14.6E\n" % dt)
+ fout.write("# num_timesteps = %8d\n" % nsteps)
+ fout.write(locHeader % (body[iloc],
+ strike[iloc],
+ dip))
+ fout.write(headerB)
+ data = numpy.transpose((timeStamps,
+ +disp[:,iloc,1],
+ +vel[:,iloc,1],
+ -disp[:,iloc,2],
+ -vel[:,iloc,2],
+ -disp[:,iloc,0],
+ -vel[:,iloc,0]))
+ numpy.savetxt(fout, data, fmt='%14.6e')
+ fout.close()
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_offfaultdata.py
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_onfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_onfaultdata.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_onfaultdata.py 2012-01-20 20:02:44 UTC (rev 19401)
@@ -0,0 +1,111 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+
+sim = "tpv16"
+cell = "tet4"
+dx = 150
+
+inputRoot = "output/%s_%s_%3dm-fault" % (sim,cell,dx)
+outdir = "scecfiles/%s_%s_%3dm/" % (sim,cell,dx)
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import time
+
+# ----------------------------------------------------------------------
+targets = numpy.array([[ 0.0, -9000.0, 0.0],
+ [ 0.0, 0.0, 0.0],
+ [ 0.0, +9000.0, 0.0],
+ [ 0.0, -9000.0, -9000.0],
+ [ 0.0, 0.0, -9000.0],
+ [ 0.0, +9000.0, -9000.0]])
+tolerance = 1.0e-6
+
+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
+ min = numpy.min(dist)
+ indices.append(numpy.where(dist <= min+tolerance)[0][0])
+
+print "Indices: ", indices
+print "Coordinates of selected points:\n",vertices[indices,:]
+
+
+# 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 = TPV16\n" + \
+ "# author = Brad Aagaard\n" + \
+ "# date = %s\n" % (time.asctime()) + \
+ "# code = PyLith\n" + \
+ "# code_version = 1.7.0a (scecdynrup branch)\n" + \
+ "# element_size = %s\n" % dx
+headerB = \
+ "# 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" + \
+ "# Column #4 = horizontal right-lateral shear stress (MPa)\n" + \
+ "# Column #5 = vertical up-dip slip (m)\n" + \
+ "# Column #6 = vertical up-dip slip-rate (m/s)\n" + \
+ "# Column #7 = vertical up-dip shear stress (MPa)\n" + \
+ "# Column #8 = normal stress 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" + \
+ "#\n"
+
+locHeader = "# location = on fault, %3.1f km along strike and %3.1f km depth\n"
+locName = "st%+04ddp%03d"
+
+lengthScale = 1.0e+3
+strike = targets[:,1] / lengthScale
+dip = -targets[:,2] / lengthScale
+
+for iloc in xrange(ntargets):
+ pt = locName % (round(10*strike[iloc]),
+ round(10*dip[iloc]))
+ filename = "%sfault%s.dat" % (outdir,pt)
+ fout = open(filename, 'w');
+ fout.write(headerA)
+ fout.write("# time_step = %14.6E\n" % dt)
+ fout.write("# num_timesteps = %8d\n" % nsteps)
+ fout.write(locHeader % (strike[iloc], dip[iloc]))
+ fout.write(headerB)
+ 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,
+ +traction[:,iloc,2]/1e+6))
+ numpy.savetxt(fout, data, fmt='%14.6e')
+ fout.close()
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_onfaultdata.py
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_ruptime.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_ruptime.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_ruptime.py 2012-01-20 20:02:44 UTC (rev 19401)
@@ -0,0 +1,87 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+
+sim = "tpv16"
+cell = "tet4"
+dx = 150
+
+inputRoot = "output/%s_%s_%3dm-fault" % (sim,cell,dx)
+outdir = "scecfiles/%s_%s_%3dm/" % (sim,cell,dx)
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import time
+
+# ----------------------------------------------------------------------
+threshold = 0.001 # threshold for detecting slip has started
+maxTime = 1.0e+10 # Large value for default rupture time
+
+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]
+
+h5.close()
+
+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 = timeStamps[itime]
+
+ # Compute magnitude of slip rate
+ slipRateMag = (slipRate[itime,:,0]**2 + slipRate[itime,:,1]**2)**0.5
+
+ # Set rupture time at locations where threshold is exceeded
+ mask = slipRateMag > threshold
+ tmpTime[:] = maxTime
+ tmpTime[mask] = t
+
+ indices = numpy.where(tmpTime < rupTime)[0]
+ rupTime[indices] = t
+
+ itime += 1
+
+# Write data
+headerA = \
+ "# problem = TPV16\n" + \
+ "# author = Brad Aagaard\n" + \
+ "# date = %s\n" % (time.asctime()) + \
+ "# code = PyLith\n" + \
+ "# code_version = 1.7.0a (scecdynrup branch)\n" + \
+ "# element_size = %s\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" + \
+ "# Column #3 = Rupture time (s)\n" + \
+ "#\n" + \
+ "# Data fields\n" + \
+ "j k t\n" + \
+ "#\n"
+
+
+distDip = -vertices[:,2]
+distStrike = vertices[:,1]
+
+filename = "%sruptime.dat" % (outdir)
+fout = open(filename, "w")
+fout.write(headerA)
+data = numpy.transpose((distStrike, distDip, rupTime))
+numpy.savetxt(fout, data, fmt='%14.6e')
+fout.close()
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tabulate_ruptime.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4.cfg 2012-01-20 19:15:44 UTC (rev 19400)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4.cfg 2012-01-20 20:02:44 UTC (rev 19401)
@@ -51,4 +51,3 @@
[pylithapp.timedependent.interfaces.fault]
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.shape = triangle
-quadrature.cell.quad_order = 2
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4_150m.jou
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4_150m.jou 2012-01-20 19:15:44 UTC (rev 19400)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tet4_150m.jou 2012-01-20 20:02:44 UTC (rev 19401)
@@ -1,5 +1,5 @@
# ----------------------------------------------------------------------
-# Create tet4 mesh at 100m resolution.
+# Create tet4 mesh at 150m resolution.
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
@@ -44,3 +44,6 @@
# Export exodus file
# ----------------------------------------------------------------------
export mesh "tet4_150m.exo" dimension 3 overwrite
+
+
+# End of file
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_075m.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_075m.cfg 2012-01-20 19:15:44 UTC (rev 19400)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_075m.cfg 2012-01-20 20:02:44 UTC (rev 19401)
@@ -12,10 +12,10 @@
# problem
# ----------------------------------------------------------------------
[pylithapp.timedependent.normalizer]
-wave_period = 0.3*s
+wave_period = 0.2*s
[pylithapp.timedependent.formulation.time_step]
-dt = 0.005*s
+dt = 0.003125*s
# ----------------------------------------------------------------------
# output
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_150m.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_150m.cfg 2012-01-20 19:15:44 UTC (rev 19400)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/tpv16_tet4_150m.cfg 2012-01-20 20:02:44 UTC (rev 19401)
@@ -14,7 +14,7 @@
wave_period = 0.4*s
[pylithapp.timedependent.formulation.time_step]
-dt = 0.005*s
+dt = 0.00625*s
# ----------------------------------------------------------------------
# output
More information about the CIG-COMMITS
mailing list