[cig-commits] r16607 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d
surendra at geodynamics.org
surendra at geodynamics.org
Fri Apr 30 17:56:26 PDT 2010
Author: surendra
Date: 2010-04-30 17:56:26 -0700 (Fri, 30 Apr 2010)
New Revision: 16607
Added:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.exo.gz
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.jou
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py
Modified:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/pylithapp.cfg
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_200m.cfg
Log:
Added mesh of 100m resolution for tpv205-2d and post-processing python scripts to extract data needed by SCEC benchmarkwebsite from VTK files
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/pylithapp.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/pylithapp.cfg 2010-05-01 00:51:04 UTC (rev 16606)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/pylithapp.cfg 2010-05-01 00:56:26 UTC (rev 16607)
@@ -61,7 +61,7 @@
bc.y_pos = pylith.bc.AbsorbingDampers
[pylithapp.timedependent.formulation.time_step]
-total_time = 1.00*s
+total_time = 12.00*s
dt = 0.01*s
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.cfg (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.cfg 2010-05-01 00:56:26 UTC (rev 16607)
@@ -0,0 +1,37 @@
+# -*- Python -*-
+
+[pylithapp]
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[pylithapp.mesh_generator.reader]
+filename = quad4_100m.exo
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+# Give basename for VTK domain output of solution over domain.
+[pylithapp.problem.formulation.output.output]
+vertex_data_fields=[displacement,velocity]
+output_freq = time_step
+time_step = 0.04999*s
+writer.filename = output/quad4_100m.vtk
+writer.time_format = %05.2f
+
+# Give basename for VTK fault output.
+[pylithapp.timedependent.interfaces.fault.output]
+vertex_info_fields = [strike_dir,normal_dir,initial_traction]
+vertex_data_fields = [slip,slip_rate,traction]
+output_freq = time_step
+time_step = 0.04999*s
+writer.filename = output/quad4_100m-fault.vtk
+writer.time_format = %05.2f
+
+# Give basename for VTK output of state variables.
+[pylithapp.timedependent.materials.elastic.output]
+cell_filter = pylith.meshio.CellFilterAvgMesh
+output_freq = time_step
+time_step = 12.0*s
+writer.filename = output/quad4_100m-elastic.vtk
+writer.time_format = %05.2f
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.exo.gz
===================================================================
(Binary files differ)
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.exo.gz
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.jou
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.jou (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_100m.jou 2010-05-01 00:56:26 UTC (rev 16607)
@@ -0,0 +1,30 @@
+# ----------------------------------------------------------------------
+# Create quad4 mesh at 100m resolution.
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# Generate geometry
+# ----------------------------------------------------------------------
+playback 'geometry.jou'
+
+# ----------------------------------------------------------------------
+# Set discretization size
+# ----------------------------------------------------------------------
+surface all size 100
+
+# ----------------------------------------------------------------------
+# Generate the mesh
+# ----------------------------------------------------------------------
+mesh surface all
+
+# ----------------------------------------------------------------------
+# Setup boundary conditions.
+# ----------------------------------------------------------------------
+playback 'createbc.jou'
+
+# ----------------------------------------------------------------------
+# Export exodus file
+# ----------------------------------------------------------------------
+export mesh "quad4_100m.exo" dimension 2 overwrite
+
+
Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_200m.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_200m.cfg 2010-05-01 00:51:04 UTC (rev 16606)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/quad4_200m.cfg 2010-05-01 00:56:26 UTC (rev 16607)
@@ -13,16 +13,18 @@
# ----------------------------------------------------------------------
# Give basename for VTK domain output of solution over domain.
[pylithapp.problem.formulation.output.output]
+vertex_data_fields=[displacement,velocity]
output_freq = time_step
-time_step = 0.05*s
+time_step = 0.04999*s
writer.filename = output/quad4_200m.vtk
writer.time_format = %05.2f
# Give basename for VTK fault output.
[pylithapp.timedependent.interfaces.fault.output]
vertex_info_fields = [strike_dir,normal_dir,initial_traction]
+vertex_data_fields = [slip,slip_rate,traction]
output_freq = time_step
-time_step = 0.05*s
+time_step = 0.04999*s
writer.filename = output/quad4_200m-fault.vtk
writer.time_format = %05.2f
@@ -30,6 +32,6 @@
[pylithapp.timedependent.materials.elastic.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
output_freq = time_step
-time_step = 1.0*s
+time_step = 12.0*s
writer.filename = output/quad4_200m-elastic.vtk
writer.time_format = %05.2f
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_offfaultdata.py 2010-05-01 00:56:26 UTC (rev 16607)
@@ -0,0 +1,138 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+outputRoot = "output/quad4_100m"
+outdir = "scecfiles/quad4_100m/"
+dx = 100
+
+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
+
+# ----------------------------------------------------------------------
+# Extract Off Fault stations data ( Displacements and Velocities )
+# ----------------------------------------------------------------------
+nfileNames = numpy.arange(1,1202,5)
+filename = "%s_t%04d.vtk" % (outputRoot,nfileNames[0])
+
+dt = 0.05
+nsteps = len(nfileNames)
+print nsteps
+
+fieldNames = ["displacement","velocity"]
+reader = VTKDataReader()
+data = reader.read(filename)
+values = numpy.array(data['vertices'])
+(nvertices, spaceDim) = values.shape
+
+x_OffFltStat = numpy.array([3000.0, 3000.0])
+y_OffFltStat = numpy.array([-12000.0, 12000.0]) # CHECK THESE VALUES
+Index_OffFltStat = numpy.zeros(2)
+
+# Find the indices of On Fault Stations
+nFltStat = -1
+for i in xrange(nvertices):
+ if (feq(values[i,0],x_OffFltStat[0]) and feq(values[i,1],y_OffFltStat[0])) or \
+ (feq(values[i,0],x_OffFltStat[1]) and feq(values[i,1],y_OffFltStat[1])):
+ nFltStat += 1
+ Index_OffFltStat[nFltStat] = i
+
+Index_OffFltStat = numpy.int_(Index_OffFltStat)
+print "Off-fault indices", Index_OffFltStat
+
+# Extract values
+disp = numpy.zeros((nsteps,2,3)) # 3-D array (ntime, nstats, ncomponents)
+vel = numpy.zeros((nsteps,2,3))
+nFile = -1
+for i in nfileNames:
+ nFile += 1
+ filename = "%s_t%04d.vtk" % (outputRoot,i)
+ data = reader.read(filename)
+ nField = 0
+ for name in fieldNames:
+ nField += 1
+ values = data['vertex_fields'][name]
+ if nField == 1:
+ disp[nFile,:,:] = values[Index_OffFltStat,:]
+ elif nField == 2:
+ vel[nFile,:,:] = values[Index_OffFltStat,:]
+
+
+# ----------------------------------------------------------------------
+# Extract Off Fault stations data ( Displacements and Velocities )
+# ----------------------------------------------------------------------
+
+headerA = \
+ "# problem = TPV205\n" + \
+ "# author = BradAagaard\n" + \
+ "# date = %s\n" % (time.asctime()) + \
+ "# code = PyLith\n" + \
+ "# code_version = 1.4.3\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 = 100.0
+dip = 7.5
+body = x_OffFltStat[:] / lengthScale
+strike = y_OffFltStat[:] / lengthScale
+time = (nfileNames[:] - 1) / timeScale
+print "time", time
+
+nlocs = disp.shape[1]
+for iloc in xrange(nlocs):
+ pt = locName % (round(10*body[iloc]),
+ round(10*strike[iloc]),
+ round(10*dip))
+ filename = "%sbody%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 % (body[iloc],
+ strike[iloc],
+ dip))
+ fout.write(headerB)
+ data = numpy.transpose((time,
+ +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()
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv205-2d/tabulate_onfaultdata.py 2010-05-01 00:56:26 UTC (rev 16607)
@@ -0,0 +1,144 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+outputRoot = "output/quad4_100m"
+outdir = "scecfiles/quad4_100m/"
+dx = 100
+
+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
+
+# ----------------------------------------------------------------------
+# 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()
+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)
+
+# 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
+
+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)
+ 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,:]
+
+
+
+# ----------------------------------------------------------------------
+# Write On Fault stations data to file ( Slip, Slip Rate, Tractions )
+# ----------------------------------------------------------------------
+
+headerA = \
+ "# problem = TPV205\n" + \
+ "# author = BradAagaard\n" + \
+ "# date = %s\n" % (time.asctime()) + \
+ "# code = PyLith\n" + \
+ "# code_version = 1.4.3\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" + \
+ "#\n" + \
+ "# Data fields\n" + \
+ "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress\n" + \
+ "#\n"
+
+locHeader = "# location = on fault, %3.1f km along strike and %3.1f km depth\n"
+locName = "st%+04ddp%03d"
+
+lengthScale = 1000.0
+timeScale = 100.0
+dip = 7.5
+strike = y_OnFltStat[:] / lengthScale
+time = (nfileNames[:] - 1) / timeScale
+print "time", time
+
+# Write data
+nlocs = slip.shape[1]
+for iloc in xrange(nlocs):
+ pt = locName % (round(10*strike[iloc]),
+ round(10*dip))
+ 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))
+ fout.write(headerB)
+ data = numpy.transpose((time,
+ -slip[:,iloc,0],
+ -slip_rate[:,iloc,0],
+ -traction[:,iloc,0]/1e+6,
+ +slip[:,iloc,1],
+ +slip_rate[:,iloc,1],
+ +traction[:,iloc,1]/1e+6))
+ numpy.savetxt(fout, data, fmt='%14.6e')
+ fout.close()
More information about the CIG-COMMITS
mailing list