[cig-commits] r19938 - short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d

brad at geodynamics.org brad at geodynamics.org
Thu Apr 12 11:47:55 PDT 2012


Author: brad
Date: 2012-04-12 11:47:55 -0700 (Thu, 12 Apr 2012)
New Revision: 19938

Added:
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_sliprate.py
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_stressslip.py
Modified:
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/geometry.jou
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/pylithapp.cfg
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_200m.cfg
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv13.cfg
   short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tri3.cfg
Log:
Added some plotting scripts.

Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/geometry.jou
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/geometry.jou	2012-04-11 19:21:15 UTC (rev 19937)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/geometry.jou	2012-04-12 18:47:55 UTC (rev 19938)
@@ -10,7 +10,7 @@
 # ----------------------------------------------------------------------
 # Set units to SI.
 # ----------------------------------------------------------------------
-# {Units('si')}
+#{Units('si')}
 #
 # ----------------------------------------------------------------------
 # Reset geometry.
@@ -36,6 +36,7 @@
 #{faultWidth=15.0*km}
 #
 #{xoffset=0.5*faultWidth*cosd(faultDipAngle)}
+
 brick x {blockWidth} y {blockHeight} z {blockWidth}
 
 # Translate block so the top is at z=0

Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_sliprate.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_sliprate.py	                        (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_sliprate.py	2012-04-12 18:47:55 UTC (rev 19938)
@@ -0,0 +1,123 @@
+#!/usr/bin/env python
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+# Plot initial stress and slip profiles.
+#
+# PREREQUISITES: matplotlib, numpy
+
+sim = "tpv13"
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import matplotlib.pyplot as pyplot
+import sys
+
+sys.path.append("../../../figures")
+import matplotlibext
+
+header = 0.95 
+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.0, 6.25, margins=[[0.45, 0.2, 0.1], [0.35, 0.3, 0.1]], dpi=150)
+
+locs = [0, 3.0, 7.5, 12.0]
+
+nrows = len(locs)
+ncols = 2
+
+for icol in xrange(ncols):
+
+    simdirs = []
+    if icol == 0:
+        labels = ["Quad4, 200m",
+                  "Tri3, 200m",
+                  "Quad4, 100m",
+                  "Tri3 100m",
+                  "Quad4, 50m",
+                  "Tri3, 50m",
+                  ]
+        for dx in [200,100,50]:
+            for cell in ["quad4", "tri3"]:
+                label = "%s, %dm" % (cell.capitalize(), dx)
+                d = "scecfiles/%s_%s_%03dm" % (sim,cell,dx)
+                simdirs.append((label, d))
+    else:
+        cell = "tri3"
+        dx = 100
+        modelers = [('Andrews', "andrews"),
+                    ('Barall', "barall"),
+                    ('Dunham', "dunham"),
+                    ('Ma', "ma"),
+                    ('PyLith (Tri3)', "tri3"),
+                    ]
+        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
+        for (label, simdir) in simdirs:
+            filename = "%s/faultst000dp%03d.dat" % (simdir, int(locs[irow]*10))
+            data = numpy.loadtxt(filename, comments="#", usecols=(0,5),
+                                 converters={0: getval,
+                                             5: 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, 8.0))
+        ax.set_xlabel("Time (s)")
+        ax.set_ylim((0, 16.0))
+        ax.set_ylabel("Slip Rate (m/s)")
+        ax.set_title("%3.1f km Down Dip" % locs[irow],
+                     horizontalalignment="right",
+                     fontweight='bold')
+
+        if irow == 0:
+            ax.legend(loc="lower right",
+                      bbox_to_anchor=(1,1.25), 
+                      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_sliprate" % (sim))
+
+
+# End of file


Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_sliprate.py
___________________________________________________________________
Name: svn:executable
   + *

Added: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_stressslip.py
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_stressslip.py	                        (rev 0)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_stressslip.py	2012-04-12 18:47:55 UTC (rev 19938)
@@ -0,0 +1,92 @@
+#!/usr/bin/env python
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+# Plot initial stress and slip profiles.
+#
+# PREREQUISITES: matplotlib, numpy
+
+sim = "tpv13"
+cell = "tri3"
+dx = 100
+
+infoFilename = "output/%s_%s_%03dm-fault_info.h5" % (sim,cell,dx)
+dataFilename = "output/%s_%s_%03dm-fault.h5" % (sim,cell,dx)
+
+# ----------------------------------------------------------------------
+import tables
+import numpy
+import matplotlib.pyplot as pyplot
+import sys
+
+sys.path.append("../../../figures")
+import matplotlibext
+
+# ----------------------------------------------------------------------
+h5 = tables.openFile(infoFilename, 'r')
+coefdyn = h5.root.vertex_fields.dynamic_coefficient[:].squeeze()
+coefstatic = h5.root.vertex_fields.static_coefficient[:].squeeze()
+cohesion = h5.root.vertex_fields.cohesion[:].squeeze() / 1.0e+6
+h5.close()
+
+h5 = tables.openFile(dataFilename, 'r')
+vertices = h5.root.geometry.vertices[:]
+slip = h5.root.vertex_fields.slip[:].squeeze()
+traction = h5.root.vertex_fields.traction[1,:,:].squeeze() / 1.0e+6
+h5.close()
+
+from math import pi, sin
+dipDist = -vertices[:,1]/sin(60.0*pi/180.0)/1.0e+3
+order = numpy.argsort(dipDist)
+
+figure = matplotlibext.Figure()
+figure.open(3.0, 5.25, margins=[[0.45, 0, 0.15], [0.35, 0.55, 0.1]], dpi=150)
+
+ax = figure.axes(2.0, 1, 1.0, 1)
+ax.plot(traction[order,0], dipDist[order], 
+        color='red', 
+        linewidth=1,
+        dashes=(None,None))
+ax.hold(True)
+ax.plot(traction[order,1], dipDist[order], 
+        color='blue', 
+        linewidth=1,
+        dashes=(3.0,1.5))
+
+ax.plot(cohesion-coefstatic[order]*traction[order,1], dipDist[order], 
+        color='orange', 
+        linewidth=1,
+        dashes=(6.0,1.5))
+ax.plot(cohesion-coefdyn[order]*traction[order,1], dipDist[order], 
+        color='green', 
+        linewidth=1,
+        dashes=(1.5,1.5))
+
+ax.set_xlim((-150,150))
+ax.set_xlabel("Traction (MPa)")
+ax.set_ylim((15.0, 0.0))
+ax.set_ylabel("Dist. Down Dip (km)")
+ax.legend(('$T_\mathit{shear}$', 
+           '$T_\mathit{normal}$',
+           '$T_\mathit{failure}$',
+           '$T_\mathit{sliding}$',
+           ), loc="upper left")
+
+ax = figure.axes(2.0, 1, 2.0, 1)
+itime = slip.shape[0]-1
+ax.plot(slip[itime,order,0], dipDist[order], 
+        color='red', 
+        linewidth=1,
+        dashes=(None,None))
+ax.set_xlabel("Slip (m)")
+ax.set_ylim((15.0, 0.0))
+ax.set_ylabel("Dist. Down Dip (km)")
+
+pyplot.show()
+pyplot.savefig("%s_%s_%03dm_stressslip" % (sim,cell,dx))
+
+
+# End of file


Property changes on: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/plot_stressslip.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-04-11 19:21:15 UTC (rev 19937)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/pylithapp.cfg	2012-04-12 18:47:55 UTC (rev 19938)
@@ -1,5 +1,3 @@
-# -*- Python -*-
-
 [pylithapp]
 
 # ----------------------------------------------------------------------
@@ -120,7 +118,7 @@
 writer = pylith.meshio.DataWriterHDF5Mesh
 
 [pylithapp.timedependent.interfaces.fault.output]
-vertex_info_fields = [strike_dir,normal_dir]
+vertex_info_fields = [strike_dir,normal_dir,static_coefficient,dynamic_coefficient,cohesion]
 vertex_data_fields = [slip,slip_rate,traction]
 output_freq = time_step
 time_step = 0.04999*s

Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg	2012-04-11 19:21:15 UTC (rev 19937)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10.cfg	2012-04-12 18:47:55 UTC (rev 19938)
@@ -1,4 +1,3 @@
-# -*- Python -*-
 [pylithapp]
 
 # ----------------------------------------------------------------------

Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_200m.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_200m.cfg	2012-04-11 19:21:15 UTC (rev 19937)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv10_tri3_200m.cfg	2012-04-12 18:47:55 UTC (rev 19938)
@@ -1,4 +1,3 @@
-# -*- Python -*-
 [pylithapp]
 
 # ----------------------------------------------------------------------

Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv13.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv13.cfg	2012-04-11 19:21:15 UTC (rev 19937)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tpv13.cfg	2012-04-12 18:47:55 UTC (rev 19938)
@@ -20,6 +20,7 @@
 # ----------------------------------------------------------------------
 [pylithapp.timedependent.materials]
 elastic = pylith.materials.DruckerPragerPlaneStrain
+elastic.allow_tensile_yield = True
 
 [pylithapp.timedependent.materials.elastic]
 db_initial_stress = spatialdata.spatialdb.SimpleDB

Modified: short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tri3.cfg
===================================================================
--- short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tri3.cfg	2012-04-11 19:21:15 UTC (rev 19937)
+++ short/3D/PyLith/benchmarks/trunk/dynamic/scecdynrup/tpv210-2d/tri3.cfg	2012-04-12 18:47:55 UTC (rev 19938)
@@ -1,4 +1,3 @@
-# -*- Python -*-
 [pylithapp]
 
 # ----------------------------------------------------------------------



More information about the CIG-COMMITS mailing list