[cig-commits] r20589 - in short/3D/PyLith/branches/v1.7-trunk: . applications/utilities pylith pylith/apps tests_auto tests_auto/eqinfo
brad at geodynamics.org
brad at geodynamics.org
Fri Aug 17 17:54:59 PDT 2012
Author: brad
Date: 2012-08-17 17:54:58 -0700 (Fri, 17 Aug 2012)
New Revision: 20589
Added:
short/3D/PyLith/branches/v1.7-trunk/applications/utilities/pylith_eqinfo.in
short/3D/PyLith/branches/v1.7-trunk/pylith/apps/EqInfoApp.py
short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/
short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/eqinfoapp.cfg
short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/generate.py
short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/line_one.h5
Modified:
short/3D/PyLith/branches/v1.7-trunk/applications/utilities/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/configure.ac
short/3D/PyLith/branches/v1.7-trunk/pylith/Makefile.am
Log:
Started work on pylith_eqinfo to compute rupture stats.
Modified: short/3D/PyLith/branches/v1.7-trunk/applications/utilities/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/applications/utilities/Makefile.am 2012-08-17 15:21:44 UTC (rev 20588)
+++ short/3D/PyLith/branches/v1.7-trunk/applications/utilities/Makefile.am 2012-08-18 00:54:58 UTC (rev 20589)
@@ -19,6 +19,7 @@
bin_SCRIPTS = \
pylithinfo \
pylith_genxdmf \
+ pylith_eqinfo \
powerlaw_gendb.py
@@ -36,6 +37,10 @@
$(do_build) < $(srcdir)/pylith_genxdmf.in > $@ || (rm -f $@ && exit 1)
chmod +x $@
+pylith_eqinfo: $(srcdir)/pylith_eqinfo.in Makefile
+ $(do_build) < $(srcdir)/pylith_eqinfo.in > $@ || (rm -f $@ && exit 1)
+ chmod +x $@
+
install-binSCRIPTS: $(bin_SCRIPTS)
@$(NORMAL_INSTALL)
test -z "$(bindir)" || $(mkdir_p) "$(DESTDIR)$(bindir)"
@@ -58,11 +63,13 @@
EXTRA_DIST = \
pylithinfo.in \
pylith_genxdmf.in \
+ pylith_eqinfo.in \
powerlaw_gendb.py
CLEANFILES = \
pylithinfo \
- pylith_genxdmf
+ pylith_genxdmf \
+ pylith_eqinfo
# End of file
Added: short/3D/PyLith/branches/v1.7-trunk/applications/utilities/pylith_eqinfo.in
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/applications/utilities/pylith_eqinfo.in (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/applications/utilities/pylith_eqinfo.in 2012-08-18 00:54:58 UTC (rev 20589)
@@ -0,0 +1,54 @@
+#!@INTERPRETER@
+# -*- Python -*-
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+# This script creates a Python file with earthquake rupture
+# information computed from PyLith output. The rupture information
+# includes:
+#
+# Rupture area
+# Average slip
+# Seismic potency
+# Seismic moment
+# Moment magnitude
+#
+# Usage: pylith_eqinfo [command line arguments]
+#
+# NOTE: Works with HDF5 files, not VTK files.
+
+__requires__ = "PyLith"
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+
+ # re-create the PYTHONPATH at 'configure' time
+ import os.path, sys, site
+ path = '@PYTHONPATH@'.split(':')
+ path.reverse()
+ for directory in path:
+ if directory:
+ directory = os.path.abspath(directory)
+ sys.path.insert(1, directory)
+ site.addsitedir(directory)
+
+ from pylith.apps.EqInfoApp import EqInfoApp
+ from pyre.applications import start
+ start(applicationClass=EqInfoApp)
+
+# End of file
Property changes on: short/3D/PyLith/branches/v1.7-trunk/applications/utilities/pylith_eqinfo.in
___________________________________________________________________
Name: svn:executable
+ *
Modified: short/3D/PyLith/branches/v1.7-trunk/configure.ac
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/configure.ac 2012-08-17 15:21:44 UTC (rev 20588)
+++ short/3D/PyLith/branches/v1.7-trunk/configure.ac 2012-08-18 00:54:58 UTC (rev 20589)
@@ -317,6 +317,7 @@
tests_auto/2d/tri3/Makefile
tests_auto/2d/quad4/Makefile
tests_auto/petsc/Makefile
+ tests_auto/eqinfo/Makefile
tests/Makefile
tests/2d/Makefile
tests/2d/maxwell/Makefile
Modified: short/3D/PyLith/branches/v1.7-trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/pylith/Makefile.am 2012-08-17 15:21:44 UTC (rev 20588)
+++ short/3D/PyLith/branches/v1.7-trunk/pylith/Makefile.am 2012-08-18 00:54:58 UTC (rev 20589)
@@ -21,6 +21,7 @@
apps/__init__.py \
apps/PyLithApp.py \
apps/PetscApplication.py \
+ apps/EqInfoApp.py \
apps/PrepMeshApp.py \
bc/__init__.py \
bc/AbsorbingDampers.py \
Added: short/3D/PyLith/branches/v1.7-trunk/pylith/apps/EqInfoApp.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/pylith/apps/EqInfoApp.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/pylith/apps/EqInfoApp.py 2012-08-18 00:54:58 UTC (rev 20589)
@@ -0,0 +1,313 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/apps/PyLithApp.py
+##
+## @brief Python PyLith application
+
+from pyre.applications.Script import Script as Application
+
+import h5py
+import numpy
+import os
+
+# ======================================================================
+class RuptureStats(object):
+ """
+ Python object to hold rupture stats.
+ """
+
+ def __init__(self, nsnapshots):
+ self.fault = None
+ self.timestamp = numpy.zeros( (nsnapshots,), dtype=numpy.float64)
+ self.ruparea = numpy.zeros( (nsnapshots,), dtype=numpy.float64)
+ self.potency = numpy.zeros( (nsnapshots,), dtype=numpy.float64)
+ self.moment = numpy.zeros( (nsnapshots,), dtype=numpy.float64)
+ return
+
+ def update(self, isnapshot, timestamp, ruparea, potency, moment):
+ self.timestamp[isnapshot] = timestamp
+ self.ruparea[isnapshot] = ruparea
+ self.potency[isnapshot] = potency
+ self.moment[isnapshot] = moment
+ self.recalculate()
+ return
+
+
+ def recalculate(self):
+ self.avgslip = self.potency / self.ruparea
+ self.mommag = 2.0/3.0*numpy.log10(self.moment) - 10.7
+ return
+
+
+ def writeObj(self, fout):
+ fout.write("class RuptureStats(object):\n"
+ " pass\n")
+ return
+
+
+ def write(self, fout):
+ fout.write("%s = RuptureStats()\n" % self.fault)
+
+ self._writeArray("timestamp", fout)
+ self._writeArray("ruparea", fout)
+ self._writeArray("potency", fout)
+ self._writeArray("moment", fout)
+ self._writeArray("avgslip", fout)
+ self._writeArray("mommag", fout)
+
+ return
+
+
+ def _writeArray(self, name, fout):
+
+ g = ("%14.6e" % v for v in self.__getattribute__(name))
+ astr = ", ".join(g)
+ fout.write("%s.%s = [%s]\n" % (self.fault, name, astr))
+ return
+
+
+# ======================================================================
+# EqInfoApp class
+class EqInfoApp(Application):
+ """
+ Python EqInfoApp application.
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(Application.Inventory):
+ """
+ Python object for managing EqInfoApp facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing EqInfoApp facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b faults Array of fault names.
+ ## @li \b filename_pattern Pattern for fault files.
+ ## @li \b snapshots Array of timestamps for slip snapshots.
+ ## @li \b snapshotUnits Units for timestamps in array of snapshots.
+ ## @li \b output_filename Filename for output.
+ ##
+ ## \b Facilities
+ ## @li \b db_properties Spatial database for elastic properties.
+ ## @li \b coordsys Coordinate system associated with mesh.
+
+ import pyre.inventory
+
+ faults = pyre.inventory.list("faults", default=[])
+ faults.meta['tip'] = "Array of fault names."
+
+ filenamePattern = pyre.inventory.str("filename_pattern", default="output/fault_%s.h5")
+ filenamePattern.meta['tip'] = "Pattern for fault files."
+
+ snapshots = pyre.inventory.list("snapshots", default=[-1])
+ snapshots.meta['tip'] = "Array of timestamps for slip snapshots (-1 == last time step)."
+
+ from pyre.units.time import second
+ snapshotUnits = pyre.inventory.dimensional("snapshot_units", default=1*second)
+ snapshotUnits.meta['tip'] = "Units for timestamps in array of snapshots."
+
+ filenameOut = pyre.inventory.str("output_filename", default="eqstats.py")
+ filenameOut.meta['tip'] = "Filename for output."
+
+ from spatialdata.spatialdb.SimpleDB import SimpleDB
+ dbProps = pyre.inventory.facility("db_properties", family="spatial_database", factory=SimpleDB)
+ dbProps.meta['tip'] = "Spatial database for elastic properties."
+
+ from spatialdata.geocoords.CSCart import CSCart
+ cs = pyre.inventory.facility("coordsys", family="coordsys", factory=CSCart)
+ cs.meta['tip'] = "Coordinate system associated with mesh."
+
+ typos = pyre.inventory.str("typos", default="pedantic",
+ validator=pyre.inventory.choice(['relaxed', 'strict', 'pedantic']))
+ typos.meta['tip'] = "Specifies the handling of unknown properties and " \
+ "facilities"
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="eqinfoapp"):
+ """
+ Constructor.
+ """
+ Application.__init__(self, name)
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Run the application.
+ """
+ nfaults = len(self.faults)
+
+ if nfaults == 0:
+ raise ValueError("No faults specified")
+
+ nsnapshots = len(self.snapshots)
+ statsFaults = [RuptureStats(nsnapshots)]*nfaults
+
+ ifault = 0
+ for fault in self.faults:
+ filenameIn = self.filenamePattern % fault
+ if not os.path.isfile(filenameIn):
+ raise IOError("Could not open PyLith fault output file '%s'." % filenameIn)
+
+ h5 = h5py.File(filenameIn, "r", driver='sec2')
+ vertices = h5['geometry/vertices'][:]
+ cells = h5['topology/cells'][:]
+ timestamps = h5['time'][:]
+ cellsArea = self._calcCellArea(cells, vertices)
+ cellsShearMod = self._getShearModulus(cells, vertices)
+
+ stats = statsFaults[ifault]
+ stats.fault = fault
+
+ isnapshot = 0
+ for snapshot in self.snapshots:
+ # Get slip at snapshot
+ istep = self._findTimeStep(snapshot, timestamps)
+ slip = h5['vertex_fields/slip'][istep,:,:]
+
+ slipMag = self._vectorMag(slip)
+ cellsSlipMag = self._ptsToCells(slipMag, cells)
+ mask = cellsSlipMag > 0.0
+
+ ruparea = numpy.sum(cellsArea[mask])
+ potency = numpy.sum(cellsSlipMag*cellsArea)
+ moment = numpy.sum(cellsSlipMag*cellsArea*cellsShearMod)
+
+ stats.update(isnapshot, timestamp=snapshot, ruparea=ruparea, potency=potency, moment=moment)
+
+ isnapshot += 1
+ h5.close()
+ ifault += 1
+
+ statsTotal = RuptureStats(nsnapshots)
+ statsTotal.fault = "all"
+ ruparea = statsTotal.ruparea
+ potency = statsTotal.potency
+ moment = statsTotal.moment
+ for s in statsFaults:
+ ruparea += s.ruparea
+ potency += s.potency
+ moment += s.moment
+
+ statsTotal.recalculate()
+
+ fout = open(self.filenameOut, "w")
+ statsTotal.writeObj(fout)
+ statsTotal.write(fout)
+ for s in statsFaults:
+ s.write(fout)
+ fout.close()
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Setup members using inventory.
+ """
+ Application._configure(self)
+ self.faults = self.inventory.faults
+ self.filenamePattern = self.inventory.filenamePattern
+ self.snapshots = map(float, self.inventory.snapshots)
+ self.snapshotUnits = self.inventory.snapshotUnits
+ self.dbProps = self.inventory.dbProps
+ self.cs = self.inventory.cs
+ self.filenameOut = self.inventory.filenameOut
+ self.typos = self.inventory.typos
+
+ return
+
+
+ def _calcCellArea(self, cells, vertices):
+ (ncells, ncorners) = cells.shape
+ if ncorners == 1:
+ area = numpy.ones( (ncells,), dtype=numpy.float64)
+ elif ncorners == 2:
+ area = self._vectorMag(vertices[cells[:,1]] - vertices[cells[:,0]])
+ elif ncorners == 3:
+ v01 = vertices[cells[:,1]] - vertices[cells[:,0]]
+ v02 = vertices[cells[:,2]] - vertices[cells[:,0]]
+ area = 0.5*self._vectorMag(numpy.cross(v01, v02))
+ elif ncorners == 4:
+ v01 = vertices[cells[:,1]] - vertices[cells[:,0]]
+ v02 = vertices[cells[:,2]] - vertices[cells[:,0]]
+ area = self._vectorMag(numpy.cross(v01, v02))
+ else:
+ raise ValueError("Unknown case for number of cell corners (%d)." % ncorners)
+ return area
+
+
+ def _getShearModulus(self, cells, vertices):
+ coords = self._ptsToCells(vertices, cells)
+ db = self.dbProps
+ db.open()
+ db.queryVals(["density","vs"])
+ (ncells, ndims) = coords.shape
+ data = numpy.zeros( (ncells, 2), dtype=numpy.float64)
+ err = numpy.zeros( (ncells,), dtype=numpy.int32)
+ db.multiquery(data, err, coords, self.cs)
+ db.close()
+ shearMod = data[:,0]*data[:,1]**2
+ return shearMod
+
+
+ def _findTimeStep(self, value, timestamps):
+ if value == -1:
+ i = len(timestamps)-1
+ else:
+ tdiff = numpy.abs(timestamps-value*self.snapshotUnits.value)
+ mindiff = numpy.min(tdiff)
+ i = numpy.where(tdiff < mindiff+1.0e-10)[0]
+ if len(i) > 1:
+ raise ValueError("Could not find snapshot %12.4e s in time stamps." % value)
+ return i
+
+
+ def _vectorMag(self, v):
+ (npts, ndims) = v.shape
+ mag = numpy.zeros( (npts,), dtype=numpy.float64)
+ for i in xrange(ndims):
+ mag += v[:,i]**2
+ mag = mag**0.5
+ return mag
+
+
+ def _ptsToCells(self, valueP, cells):
+ (ncells, ncorners) = cells.shape
+ if len(valueP.shape) > 1:
+ (nvertices, nvals) = valueP.shape
+ valueC = numpy.zeros( (ncells,nvals), dtype=numpy.float64)
+ for i in xrange(ncorners):
+ valueC[:,:] += valueP[cells[:,i],:]
+ else:
+ nvertices = valueP.shape
+ valueC = numpy.zeros( (ncells,), dtype=numpy.float64)
+ for i in xrange(ncorners):
+ valueC[:] += valueP[cells[:,i]]
+ valueC /= ncorners
+ return valueC
+
+
+# End of file
Added: short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/Makefile.am (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/Makefile.am 2012-08-18 00:54:58 UTC (rev 20589)
@@ -0,0 +1,46 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+TESTS = testeqinfo.py
+
+dist_check_SCRIPTS = testeqinfo.py
+
+dist_noinst_PYTHON = \
+ TestLine.py
+
+noinst_DATA = line_one.h5
+
+noinst_TMP = stats_line.py
+
+TESTS_ENVIRONMENT = $(PYTHON)
+
+
+# 'export' the input files by performing a mock install
+export_datadir = $(top_builddir)/tests_auto/petsc
+export-data: $(dist_noinst_PYTHON) $(dist_noinst_DATA)
+ for f in $(dist_noinst_PYTHON) $(dist_noinst_DATA); do $(install_sh_DATA) $(srcdir)/$$f $(export_datadir); done
+
+clean-data:
+ if [ "X$(top_srcdir)" != "X$(top_builddir)" ]; then for f in $(dist_noinst_PYTHON) $(dist_noinst_DATA) $(noinst_TMP); do $(RM) $(RM_FLAGS) $(export_datadir)/$$f; done; fi
+
+
+BUILT_SOURCES = export-data
+clean-local: clean-data
+CLEANFILES = *.pyc
+
+# End of file
Added: short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/eqinfoapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/eqinfoapp.cfg (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/eqinfoapp.cfg 2012-08-18 00:54:58 UTC (rev 20589)
@@ -0,0 +1,15 @@
+[eqinfoapp]
+faults = [one]
+filename_pattern = line_%s.h5
+snapshots = [0.01, 0.99]
+snapshot_units = 1.0*s
+
+output_filename = stats_line.py
+
+db_properties = spatialdata.spatialdb.UniformDB
+db_properties.label = Elastic properties
+db_properties.values = [vs, density]
+db_properties.data = [2.0e+3, 2.5e+3]
+
+coordsys.space_dim = 2
+
Added: short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/generate.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/generate.py (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/generate.py 2012-08-18 00:54:58 UTC (rev 20589)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @brief Generate HDF5 data files for eqinfo tests.
+
+import h5py
+import numpy
+
+# ----------------------------------------------------------------------
+class TestData(object):
+ """
+ Abstract base class for test data.
+ """
+
+ def __init__(self):
+ self.vertices = None
+ self.cells = None
+ self.slip = None
+ self.time = None
+ self.filename = None
+ return
+
+ def write(self):
+ h5 = h5py.File(self.filename, "w", driver='sec2')
+ h5.create_dataset('geometry/vertices', data=self.vertices)
+ h5.create_dataset('topology/cells', data=self.cells)
+ h5.create_dataset('time', data=self.time)
+ slip = h5.create_dataset('vertex_fields/slip', data=self.slip)
+ slip.attrs['vector_field_type'] = 'vector'
+ h5.close()
+ return
+
+
+# ----------------------------------------------------------------------
+class TestDataLineA(TestData):
+ """
+ Test data for fault as line.
+ """
+
+ def __init__(self):
+ TestData.__init__(self)
+ self.vertices = numpy.array([[-1.0, 0.0],
+ [ 0.0, 0.0],
+ [ 0.0, +1.5]],
+ dtype=numpy.float64)
+ self.cells = numpy.array([[0, 1],
+ [1, 2]],
+ dtype=numpy.int32)
+ self.slip = numpy.array([
+ # t=0
+ [[0.4, 0.0],
+ [1.0, 0.0],
+ [0.8, 0.0]],
+ # t=1
+ [[0.0, 0.0],
+ [0.0, 0.0],
+ [0.8, 0.0]]],
+ dtype=numpy.float64)
+ self.time = numpy.array([0.0, 1.0])
+ self.filename = "line_one.h5"
+
+
+# ======================================================================
+data = TestDataLineA()
+data.write()
+
+#data = TestDataTri()
+#data.write()
+
+#data = TestDataQuad()
+#data.write()
+
+# End of file
Property changes on: short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/generate.py
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/line_one.h5
===================================================================
(Binary files differ)
Property changes on: short/3D/PyLith/branches/v1.7-trunk/tests_auto/eqinfo/line_one.h5
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
More information about the CIG-COMMITS
mailing list