[cig-commits] r19947 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16
brad at geodynamics.org
brad at geodynamics.org
Mon Apr 16 10:59:38 PDT 2012
Author: brad
Date: 2012-04-16 10:59:37 -0700 (Mon, 16 Apr 2012)
New Revision: 19947
Added:
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_ruptime.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_sliprate.py
short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py
Log:
Started working on stress/slip contour plot.
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_ruptime.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_ruptime.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_ruptime.py 2012-04-16 17:59:37 UTC (rev 19947)
@@ -0,0 +1,126 @@
+#!/usr/bin/env python
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+# Plot initial stress and slip profiles.
+#
+# PREREQUISITES: matplotlib, numpy, gmt
+
+sim = "tpv16"
+griddata = False
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import matplotlib.pyplot as pyplot
+import sys
+import subprocess
+
+sys.path.append("../../../figures")
+import matplotlibext
+
+header = 0.45
+lineStyle = [("red", 'dashed'),
+ ("blue", 'dashdot'),
+ ("orange", 'dotted'),
+ ("black", 'solid'),
+ ]
+
+# ----------------------------------------------------------------------
+def getval(v):
+ try:
+ d = float(v)
+ except ValueError:
+ d = None
+ return d
+
+# ----------------------------------------------------------------------
+figure = matplotlibext.Figure()
+figure.open(6.0, 2.25, margins=[[0.4, 0.3, 0.1], [0.35, 0, 0.1]], dpi=150)
+
+nrows = 1
+ncols = 2
+
+for icol in xrange(ncols):
+
+ simdirs = []
+ if icol == 0:
+ labels = ["Tet4, 150m",
+ "Tet4 75m",
+ ]
+ for dx in [150,75]:
+ for cell in ["tet4"]:
+ label = "%s, %dm" % (cell.capitalize(), dx)
+ d = "scecfiles/%s_%s_%03dm" % (sim,cell,dx)
+ simdirs.append((label, d))
+ else:
+ cell = "tet4"
+ dx = 75
+ modelers = [('Barall', "barall"),
+ ('Kaneko', "kaneko"),
+ ('PyLith', "tet4"),
+ ]
+ for (label,modeler) in modelers:
+ d = "scecfiles/%s_%s_%03dm" % (sim,modeler,dx)
+ simdirs.append((label, d))
+
+ for irow in xrange(nrows):
+ ax = figure.axes(nrows+header, ncols, irow+1+header, icol+1)
+
+ isim = 0
+ labels = []
+ lines = []
+ for (label, simdir) in simdirs:
+ if griddata:
+ cmd = "awk '/^j/ {next}; { print $0 }' %s/ruptime.dat | triangulate -G%s/ruptime.grd -I75.0/75.0 -R-24.0e+3/+24.0e+3/0.0/19.5e+3 -Jx1.0 > /dev/null " % (simdir, simdir)
+ subprocess.call(cmd, shell=True)
+ cmd = "grd2xyz %s/ruptime.grd > %s/ruptime.txtgrd" % \
+ (simdir, simdir)
+ subprocess.call(cmd, shell=True)
+ filename = "%s/ruptime.txtgrd" % (simdir)
+ data = numpy.loadtxt(filename)
+
+ data = numpy.reshape(data, (261, 641, 3))
+
+ if isim != len(simdirs)-1:
+ style = lineStyle[isim]
+ else:
+ style = lineStyle[len(lineStyle)-1]
+ contours = ax.contour(data[:,:,0]/1.0e+3, data[:,:,1]/1e+3, data[:,:,2],
+ levels=numpy.arange(0.0, 8.01, 0.5),
+ colors=style[0],
+ linewidths=1.0,
+ linestyles=style[1])
+ lines.append(contours.collections[0])
+ ax.hold(True)
+ labels.append(label)
+ isim += 1
+
+ ax.set_xlim((-15.0, +15.0))
+ ax.set_xlabel("Dist. Along Strike (km)")
+ ax.set_ylim((15.0, 0.0))
+ ax.set_yticks(numpy.arange(0, 15.01, 5.0))
+ ax.set_ylabel("Dist. Down Dip (km)")
+ ax.set_aspect('equal')
+
+ if irow == 0:
+ ax.legend(lines, labels, loc="lower right",
+ bbox_to_anchor=(1,1.1),
+ borderaxespad=0)
+
+ if irow+1 < nrows:
+ ax.set_xticklabels([])
+ ax.set_xlabel("")
+ if icol > 0:
+ ax.set_title("")
+ ax.set_yticklabels([])
+ ax.set_ylabel("")
+
+pyplot.show()
+pyplot.savefig("%s_ruptime" % (sim))
+
+
+# End of file
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_ruptime.py
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_sliprate.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_sliprate.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_sliprate.py 2012-04-16 17:59:37 UTC (rev 19947)
@@ -0,0 +1,124 @@
+#!/usr/bin/env python
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+# Plot initial stress and slip profiles.
+#
+# PREREQUISITES: matplotlib, numpy
+
+sim = "tpv16"
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import matplotlib.pyplot as pyplot
+import sys
+
+sys.path.append("../../../figures")
+import matplotlibext
+
+header = 0.45
+lineStyle = [("red", (2.0, 1.0)),
+ ("blue", (4.0, 1.0)),
+ ("purple", (6.0, 1.0)),
+ ("green", (3.0, 1.0, 1.5, 1.0)),
+ ("orange", (6.0, 1.0, 1.5, 1.0)),
+ ("black", (None, None)),
+ ]
+
+# ----------------------------------------------------------------------
+def getval(v):
+ try:
+ d = float(v)
+ except ValueError:
+ d = None
+ return d
+
+# ----------------------------------------------------------------------
+figure = matplotlibext.Figure()
+figure.open(6.75, 4.5, margins=[[0.45, 0.2, 0.1], [0.35, 0.6, 0.1]], dpi=150)
+
+locs = [(-9.0, 0.0), (0.0, 0.0), (9.0, 0.0),
+ (-9.0, 9.0), (0.0, 9.0), (9.0, 9.0)]
+
+nrows = 2
+ncols = 3
+
+cell = "tet4"
+dx = 75
+simdirs = []
+modelers = [('Barall', "barall"),
+ ('Kaneko', "kaneko"),
+ ('PyLith', "tet4"),
+ ]
+for (label,modeler) in modelers:
+ d = "scecfiles/%s_%s_%03dm" % (sim,modeler,dx)
+ simdirs.append((label, d))
+
+iloc = 0
+for irow in xrange(nrows):
+ for icol in xrange(ncols):
+ ax = figure.axes(nrows+header, ncols, irow+1+header, icol+1)
+
+ isim = 0
+ for (label, simdir) in simdirs:
+ filename = "%s/faultst%+04ddp%03d.dat" % \
+ (simdir, int(locs[iloc][0]*10), int(locs[iloc][1]*10))
+ data = numpy.loadtxt(filename, comments="#", usecols=(0,2),
+ converters={0: getval,
+ 2: getval})
+ if isim != len(simdirs)-1:
+ style = lineStyle[isim]
+ else:
+ style = lineStyle[len(lineStyle)-1]
+ ax.plot(data[:,0], data[:,1], color=style[0],
+ linewidth=1,
+ dashes=style[1],
+ label=label)
+ ax.hold(True)
+ isim += 1
+
+ ax.set_xlim((0.0, 10.5))
+ ax.set_xlabel("Time (s)")
+ if irow == 0:
+ ax.set_ylim((0.0, 10.0))
+ ymin0, ymax0 = ax.get_ylim()
+ else:
+ ax.set_ylim((0.0, 5.0))
+
+ if icol == 0:
+ ymin, ymax = ax.get_ylim()
+ ax.text(-1.5, (ymax0+2.2)/ymax0*ymax,
+ "%3.1f km Down Dip" % locs[iloc][1],
+ fontweight='bold',
+ horizontalalignment='left')
+
+
+ ax.set_ylabel("Slip Rate (m/s)")
+
+ ax.set_title("%3.1f km Along Strike" % locs[iloc][0],
+ fontweight='bold',
+ horizontalalignment='center')
+
+ if irow == 0 and icol == ncols-1:
+ ax.legend(loc="lower right",
+ bbox_to_anchor=(1,1.3),
+ borderaxespad=0)
+
+ if irow+1 < nrows:
+ ax.set_xticklabels([])
+ ax.set_xlabel("")
+ if icol > 0:
+ ax.set_yticklabels([])
+ ax.set_ylabel("")
+
+ iloc += 1
+
+pyplot.show()
+pyplot.savefig("%s_sliprate" % (sim))
+
+
+# End of file
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_sliprate.py
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py 2012-04-16 17:59:37 UTC (rev 19947)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+# Plot initial stress and slip profiles.
+#
+# PREREQUISITES: matplotlib, numpy, gmt, tables
+
+sim = "tpv16"
+cell = "tet4"
+dx = 75
+griddata = True
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import matplotlib.pyplot as pyplot
+import sys
+import subprocess
+
+sys.path.append("../../../figures")
+import matplotlibext
+
+lineStyle = ("black", 'solid')
+
+# ----------------------------------------------------------------------
+figure = matplotlibext.Figure()
+figure.open(3.0, 2.5, margins=[[0.4, 0.3, 0.1], [0.35, 0.4, 0.1]], dpi=150)
+
+nrows = 2
+ncols = 1
+icol = 0
+
+for irow in xrange(nrows):
+ if griddata:
+ inputRoot = "output/%s_%s_%03dm-fault" % (sim,cell,dx)
+ h5 = tables.openFile("%s.h5" % inputRoot, 'r')
+ vertices = h5.root.geometry.vertices[:]
+ distDip = -vertices[:,2]
+ distStrike = vertices[:,1]
+ if irow == 0:
+ tractionShear = h5.root.vertex_fields.traction[1,:,0].squeeze()
+ data = numpy.transpose((distStrike, dispDip, tractionShear))
+ else:
+ slip = h5.root.vertex_fields.slip[-1,:,:].squeeze()
+ slipMag = (slip[:,0]**2 + slip[:,1]**2)**0.5
+ data = numpy.transpose((distStrike, distDip, slipMag))
+ h5.close()
+ numpy.savetxt("tmp.dat", data, fmt='%14.63')
+ cmd = "awk '/^j/ {next}; { print $0 }' tmp.dat | triangulate -Gtmp.grd -I75.0/75.0 -R-24.0e+3/+24.0e+3/0.0/19.5e+3 -Jx1.0 > /dev/null "
+ subprocess.call(cmd, shell=True)
+ cmd = "grd2xyz tmp.grd > tmp.txtgrd"
+ subprocess.call(cmd, shell=True)
+
+
+ ax = figure.axes(nrows, ncols, irow+1, icol+1)
+ filename = "tmp.txtgrd"
+ data = numpy.loadtxt(filename)
+ data = numpy.reshape(data, (261, 641, 3))
+ contours = ax.contour(data[:,:,0]/1.0e+3, data[:,:,1]/1e+3, data[:,:,2],
+ levels=numpy.arange(0.0, 10.01, 0.5),
+ colors=lineStyle[0],
+ linewidths=1.0,
+ linestyles=lineStyle[1])
+
+ ax.set_xlim((-24.0, +24.0))
+ ax.set_xlabel("Dist. Along Strike (km)")
+ ax.set_ylim((19.5, 0.0))
+ ax.set_yticks(numpy.arange(0, 20.01, 5.0))
+ ax.set_ylabel("Dist. Down Dip (km)")
+ ax.set_aspect('equal')
+
+ if irow+1 < nrows:
+ ax.set_xticklabels([])
+ ax.set_xlabel("")
+ if icol > 0:
+ ax.set_title("")
+ ax.set_yticklabels([])
+ ax.set_ylabel("")
+
+pyplot.show()
+pyplot.savefig("%s_stressslip" % (sim))
+
+
+# End of file
Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv16/plot_stressslip.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list