[cig-commits] r16134 - short/3D/PyLith/trunk/playpen/powerlaw
brad at geodynamics.org
brad at geodynamics.org
Fri Jan 15 13:10:53 PST 2010
Author: brad
Date: 2010-01-15 13:10:53 -0800 (Fri, 15 Jan 2010)
New Revision: 16134
Added:
short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.cfg
short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.py
short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_params.spatialdb
short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_points.txt
Removed:
short/3D/PyLith/trunk/playpen/powerlaw/mat_elastic.spatialdb
short/3D/PyLith/trunk/playpen/powerlaw/mat_powerlaw_in.spatialdb
short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.cfg
short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.py
Log:
Initial implementation of power-law parameters utility. Need info from Charles to complete.
Deleted: short/3D/PyLith/trunk/playpen/powerlaw/mat_elastic.spatialdb
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/mat_elastic.spatialdb 2010-01-15 18:47:34 UTC (rev 16133)
+++ short/3D/PyLith/trunk/playpen/powerlaw/mat_elastic.spatialdb 2010-01-15 21:10:53 UTC (rev 16134)
@@ -1,26 +0,0 @@
-// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting)
-//
-// This spatial database specifies the distribution of material
-// properties. In this case, the material properties are uniform.
-//
-#SPATIAL.ascii 1
-SimpleDB {
- num-values = 3 // number of material property values
- value-names = density vs vp // names of the material property values
- value-units = kg/m**3 m/s m/s // units
- num-locs = 1 // number of locations
- data-dim = 0
- space-dim = 3
- cs-data = cartesian {
- to-meters = 1.0
- space-dim = 3
- }
-}
-// Columns are
-// (1) x coordinate (m)
-// (2) y coordinate (m)
-// (3) z coordinate (m)
-// (4) density (kg/m^3)
-// (5) vs (m/s)
-// (6) vp (m/s)
-0.0 0.0 0.0 2500.0 3000.0 5291.502622129181
Deleted: short/3D/PyLith/trunk/playpen/powerlaw/mat_powerlaw_in.spatialdb
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/mat_powerlaw_in.spatialdb 2010-01-15 18:47:34 UTC (rev 16133)
+++ short/3D/PyLith/trunk/playpen/powerlaw/mat_powerlaw_in.spatialdb 2010-01-15 21:10:53 UTC (rev 16134)
@@ -1,37 +0,0 @@
-// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting)
-//
-// This spatial database defines the power-law parameters for two materials:
-// Wet granite and dry olivine. The properties are from Strehlau and
-// Meissner (1987).
-// Note that the flow constant is expressed as log10(flow-constant), where
-// the units are GPa^(-n)/s. The activation energy (Q) is expressed in
-// kJ/mol. To accommodate the strange units of the flow constant, the
-// given units of are 'None', and we specify a new parameter
-// (flow-constant-multiplier), with a value of 9, to indicate that the
-// underlying units are GPa. A value of 6 would indicate MPa, etc.
-//
-#SPATIAL.ascii 1
-SimpleDB {
- num-values = 4 // number of material property values
- value-names = flow-constant activation-energy power-law-exponent flow-constant-multiplier // names of the material property values
- value-units = None kJ/mol None None // units
- num-locs = 4 // number of locations
- data-dim = 1
- space-dim = 3
- cs-data = cartesian {
- to-meters = 1000.0
- space-dim = 3
- }
-}
-// Columns are
-// (1) x coordinate (m)
-// (2) y coordinate (m)
-// (3) z coordinate (m)
-// (4) flow-constant (None)
-// (5) activation-energy (kJ/mol)
-// (6) power-law-exponent (None)
-// (7) activation-energy-multiplier (None)
-0.0 0.0 0.0 2.0 137.0 1.5 9.0
-0.0 0.0 -24.9 2.0 137.0 1.5 9.0
-0.0 0.0 -25.1 15.5 535.0 3.5 9.0
-0.0 0.0 -400. 15.5 535.0 3.5 9.0
Deleted: short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.cfg 2010-01-15 18:47:34 UTC (rev 16133)
+++ short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.cfg 2010-01-15 21:10:53 UTC (rev 16134)
@@ -1,13 +0,0 @@
-# -*- Python -*-
-[pl3dcomputeparams]
-
-# Top-level parameters
-use_reference_stress = False
-reference_stress = 1.0*MPa
-reference_strain_rate = 1.0e-6/s
-output_db_name = mat_powerlaw.spatialdb
-db_template = db_input_power_law_properties
-
-# Top-level facilities
-# Need to specify all of these along with their facilities and properties.
-db_input_elastic_properties =
Deleted: short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.py
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.py 2010-01-15 18:47:34 UTC (rev 16133)
+++ short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.py 2010-01-15 21:10:53 UTC (rev 16134)
@@ -1,207 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file powerlaw/pl3dcomputeparams
-
-## @brief Python application to compute power-law parameters used by PyLith,
-## given a spatial database describing the temperature and another spatial
-## database describing the laboratory-derived properties for the various
-## materials. The output is another spatial database containing the power-law
-## parameters for PyLith. The generated spatial database will have the same
-## dimensions and coordinate specifications as the input properties database.
-
-import numpy
-import math
-
-from pyre.applications.Script import Script as Application
-
-class Pl3dComputeParams(Application):
- """
- Python application to compute power-law parameters used by PyLith,
- given input spatial databases describing the temperature and the
- laboratory-derived properties for the various materials. The output is
- another spatial database containing the power-law parameters for PyLith.
- The generated spatial database will have the same dimensions and coordinate
- specifications as the input properties database.
- """
-
- class Inventory(Application.Inventory):
- """
- Python object for managing Pl3dComputeParams facilities and properties.
- """
-
- ## @class Inventory
- ## Python object for managing Pl3dComputeParams facilities and properties.
- ##
- ## \b Properties
- ## @li \b use_reference_stress Use reference stress to compute reference strain rate.
- ## @li \b reference_stress Value of reference stress.
- ## @li \b reference_strain_rate Value of reference strain rate.
- ## @li \b db_template Which database to use as template for output.
- ## @li \b point_locations_file File containing points to output.
- ## @li \b output_db_name File name for output spatial database.
- ##
- ## \b Facilities
- ## @li \b db_input_power_law_properties Database of input power law properties.
- ## @li \b db_input_temperature Database of input temperature values.
- ## @li \b db_output_properties Database of output material properties.
-
- import pyre.inventory
- from pyre.units.pressure import Pa
- from pyre.units.time import s
-
- useReferenceStress = pyre.inventory.bool("use_reference_stress",
- default=False)
- useReferenceStress.meta['tip'] = "Use reference stress to compute reference strain rate."
-
- referenceStress = pyre.inventory.dimensional("reference_stress",
- default=1.0*MPa)
- referenceStress.meta['tip'] = "Reference stress value."
-
- referenceStrainRate = pyre.inventory.dimensional("reference_strain_rate",
- default=1.0e-6/s)
- referenceStrainRate.meta['tip'] = "Reference strain rate value."
-
- pointLocationsFile = pyre.inventory.str("point_locations_file",
- default="points.txt")
- pointLocationsFile.meta['tip'] = "File containing points to output."
-
- outputDbName = pyre.inventory.str("output_db_name",
- default="mat_powerlaw.spatialdb")
- outputDbName.meta['tip'] = "File name of output spatial database."
-
- from spatialdata.spatialdb.SimpleDB import SimpleDB
-
- dbInputPowerLawProperties = pyre.inventory.facility(
- "db_input_power_law_properties",
- family="spatial_database",
- factory=SimpleDB)
- dbInputPowerLawProperties.meta['tip'] = "DB for input power law properties."
-
- dbInputTemperature = pyre.inventory.facility(
- "db_input_temperature",
- family="spatial_database",
- factory=SimpleDB)
- dbInputTemperature.meta['tip'] = "DB for input temperature values."
-
- dbOutputProperties = pyre.inventory.facility(
- "db_output_properties",
- family="spatial_database",
- factory=SimpleDB)
- dbOutputProperties.meta['tip'] = "DB for output material properties."
-
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="pl3dcomputeparams"):
- Application.__init__(self, name)
-
- self.numPoints = 0
- self.points = None
-
- return
-
-
- def main(self):
- # import pdb
- # pdb.set_trace()
- self._readPoints()
- self._computeParams()
- return
-
-
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _configure(self):
- """
- Setup members using inventory.
- """
- Application._configure(self)
-
- # Parameters
- self.useReferenceStress = self.inventory.useReferenceStress
- self.referenceStress = self.inventory.referenceStress.value
- self.referenceStrainRate = self.inventory.referenceStrainRate.value
- self.pointLocationsFile = self.inventory.pointLocationsFile
- self.outputDbName = self.inventory.outputDbName
-
- # This script is presently restricted to 3D.
- self.spaceDim = 3
-
- # Facilities
- self.dbInputPowerLawProperties = self.inventory.dbInputPowerLawProperties
- self.dbInputTemperature = self.inventory.dbInputTemperature
- self.dbOutputProperties = self.inventory.dbOutputProperties
-
- return
-
-
- def _readPoints(self):
- """
- Read point coordinates to use in generating output.
- """
- inFile = file(self.pointLocationsFile, 'r')
- listCoords = []
-
- for line in inFile:
- if not line.strip():
- continue
- if re.compile('^#').search(line) is not None:
- continue
- else:
- self.numPoints += 1
- data = line.split()
- for dim in range(self.spaceDim):
- listCoords.append(float(data[dim]))
-
- self.points = numpy.array(
- listCoords, dtype=numpy.float64).reshape(self.numPoints, self.spaceDim)
-
- return
-
-
- def _computeParams(self):
- """
- Get parameters from different databases and output power-law properties
- at specified points.
- """
-
- from pyre.handbook.constants.fundamental import R
-
- # Open input databases.
- self.dbInputPowerLawProperties.open()
- self.dbInputTemperature.open()
-
- temp_queryVals = ["temperature"]
- pl_queryVals = ["flow-constant", "activation-energy",
- "power-law-exponent", "flow-constant-multiplier"]
-
- # Actually, I can't find a way of determining the number of locations
- # a priori. We may need a different type of loop that continues until
- # the data in the template database is exhausted.
- # Alternatively, we can define a set of points for the output database
- # and just loop over those. One method might be:
- # specify min and max values for each coordinate direction along with
- # either an interval or number of increments. If the min and max values are
- # equal for a direction, that direction is just held constant.
-
- for loc in range(nlocs):
-
- return
-
-
-# ----------------------------------------------------------------------
-if __name__ == '__main__':
- app = Pl3dComputeParams()
- app.run()
-
-# End of file
Copied: short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.cfg (from rev 16133, short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.cfg)
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.cfg (rev 0)
+++ short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.cfg 2010-01-15 21:10:53 UTC (rev 16134)
@@ -0,0 +1,28 @@
+# -*- Python -*-
+[powerlaw_gendb]
+
+reference_value = strain_rate
+
+[powerlaw_gendb.db_exponent]
+label = Power-law exponent
+iohandler.filename = powerlaw_params.spatialdb
+
+[powerlaw_gendb.db_activation_energy]
+label = Activation energy
+iohandler.filename = powerlaw_params.spatialdb
+
+[powerlaw_gendb.db_powerlaw_coefficient]
+label = Experimentally derived power-law coefficient
+iohandler.filename = powerlaw_params.spatialdb
+
+[powerlaw_gendb.db_temperature]
+label = Temperature
+iohandler.filename = powerlaw_params.spatialdb
+
+[powerlaw_gendb.geometry]
+data_dim = 1
+reader = spatialdata.utils.PointsStream
+reader.filename = powerlaw_points.txt
+
+[powerlaw_gendb.iohandler]
+filename = powerlaw_properties.spatialdb
Copied: short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.py (from rev 16132, short/3D/PyLith/trunk/playpen/powerlaw/pl3dcomputeparams.py)
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.py (rev 0)
+++ short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_gendb.py 2010-01-15 21:10:53 UTC (rev 16134)
@@ -0,0 +1,189 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file playpen/powerlaw/powerlaw_gendb.py
+
+## @brief Python application to compute power-law parameters used by
+## PyLith, given a spatial database describing the temperature and
+## another spatial database describing the laboratory-derived
+## properties for the various materials. The output is another spatial
+## database containing the power-law parameters for PyLith.
+
+import numpy
+import math
+
+from pyre.applications.Script import Script as Application
+
+class PowerLawApp(Application):
+ """
+ Python application to compute power-law parameters used by PyLith,
+ given input spatial databases describing the temperature and the
+ laboratory-derived properties for the various materials. The output is
+ another spatial database containing the power-law parameters for PyLith.
+ """
+
+ ## \b Properties
+ ## @li \b reference_value Indicates whether reference stress or
+ ## reference strain rate is provided as input.
+ ## @li \b reference_stress Value for reference stress.
+ ## @li \b reference_strain_rate Value for reference strain rate.
+ ##
+ ## \b Facilities
+ ## @li \b db_exponent Spatial db for power-law exponent, n.
+ ## @li \b db_activation_energy Spatial db for activation energy, Q.
+ ## @li \b db_temperature Spatial db for temperature, T.
+ ## @li \b db_powerlaw_coefficient Spatial db for power-law coefficient, Ae.
+ ## @li \b geometry Geometry for output database.
+ ## @li \b iohandler Object for writing database.
+
+ import pyre.inventory
+ from pyre.units.pressure import MPa
+ from pyre.units.time import s
+
+ refSelection = pyre.inventory.str("reference_value",
+ default="strain_rate",
+ validator=pyre.inventory.choice(['stress',
+ 'strain_rate']))
+ refSelection.meta['tip'] = "Indicates whether reference stress or " \
+ "reference strain rate is provided as input."
+
+ refStress = pyre.inventory.dimensional("reference_stress", default=1.0*MPa)
+ refStress.meta['tip'] = "Reference stress value."
+
+ refStrainRate = pyre.inventory.dimensional("reference_strain_rate",
+ default=1.0e-6/s)
+ refStrainRate.meta['tip'] = "Reference strain rate value."
+
+ from spatialdata.spatialdb.SimpleDB import SimpleDB
+
+ dbExponent = pyre.inventory.facility("db_exponent",
+ family="spatial_database",
+ factory=SimpleDB)
+ dbExponent.meta['tip'] = "Spatial db for power-law exponent, n."
+
+ dbActivationE = pyre.inventory.facility("db_activation_energy",
+ family="spatial_database",
+ factory=SimpleDB)
+ dbActivationE.meta['tip'] = "Spatial db for activation energy, Q."
+
+ dbTemperature = pyre.inventory.facility("db_tempoerature",
+ family="spatial_database",
+ factory=SimpleDB)
+ dbTemperature.meta['tip'] = "Spatial db for temperature, T."
+
+ dbAe = pyre.inventory.facility("db_powerlaw_coefficient",
+ family="spatial_database",
+ factory=SimpleDB)
+ dbAe.meta['tip'] = "Spatial db for power-law coefficient, Ae."
+
+ from spatialdata.spatialdb.generator.Geometry import Geometry
+ geometry = pyre.inventory.facility("geometry", family="geometry",
+ factory=Geometry)
+ geometry.meta['tip'] = "Geometry for output database."
+
+ from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+ iohandler = pyre.inventory.facility("iohandler", family="simpledb_io",
+ factory=SimpleIOAscii)
+ iohandler.meta['tip'] = "Object for writing database."
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="powerlaw_gendb"):
+ Application.__init__(self, name)
+
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Application driver.
+ """
+ # Get output points
+ self._info.log("Reading geometry.")
+ self.geometry.read()
+ points = self.geometry.vertices
+ coordsys = self.geometry.coordsys
+
+ (npoints, spaceDim) = points.shape
+ refStrainRate = numpy.zeros( (npoints,), dtype=numpy.float64)
+ refStress = numpy.zeros( (npoints,), dtype=numpy.float64)
+
+ # Query databases to get inputs at output points
+ self._info.log("Querying for parameters at output points.")
+ n = self._queryDB(self.dbExponent, "power-law-exponent", points, coordsys)
+ Q = self._queryDB(self.dbActivationE, "activation-energy", points, coordsys)
+ Ae = self._queryDB(self.dbAe, "power-law-coefficient", points, coordsys)
+ T = self._queryDB(self.dbTemperature, "temperature", points, coordsys)
+
+ # Compute power-law parameters
+ self._info("Computing parameters at output points.")
+ from pyre.handbook.constants.fundamental import R
+ At = 3**(0.5*(n+1))/2.0 * Ae * numpy.exp(-Q/(R*T))
+
+ if self.refSelection == "stress":
+ refStrainRate = self.refStress**n / At
+ elif self.refSelection == "strain_rate":
+ refStress = (self.refStrainRate / At)**(1.0/n)
+ else:
+ raise ValueError("Invalid value (%s) for reference value." % \
+ self.refSelection)
+
+ refStressInfo = {'name': "reference-stress",
+ 'units': "Pa",
+ 'data': refStress}
+ refStrainRateInfo = {'name': "reference-strain-rate",
+ 'units': "1/s",
+ 'data': refStrainRate}
+ exponentInfo = {'name': "powerlaw-exponent",
+ 'units': "none",
+ 'data': n}
+
+ # Write database
+ self._info.log("Writing database.")
+ data = {'points': points,
+ 'coordsys': coordsys,
+ 'data_dim': self.geometry.dataDim,
+ 'values': [refStressInfo, refStrainInfo, exponentInfo]}
+ self.iohandler.write(data)
+ return
+
+
+ def _queryDB(self, db, valueName, points, cs):
+ """
+ Query spatial database
+ """
+
+ (npoints, spaceDim) = points.shape
+ data = numpy.zeros( (npoints,1), dtype=numpy.float64)
+ err = numpy.zeros( (npoints,), dtype=numpy.int32)
+
+ db.open()
+ db.queryVals([valueName])
+ db.multiquery(data, err, points, cs)
+ db.close()
+ errSum = numpy.sum(err)
+ if errSum > 0:
+ msg = "Query for %s failed at %d points.\n" \
+ "Coordinates of points:\n" % (valueName, errSum)
+ msg += "%s" % points[err,:]
+ raise ValueError(msg)
+
+ return data
+
+
+# ----------------------------------------------------------------------
+if __name__ == '__main__':
+ app = PowerLawApp()
+ app.run()
+
+# End of file
Copied: short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_params.spatialdb (from rev 16133, short/3D/PyLith/trunk/playpen/powerlaw/mat_powerlaw_in.spatialdb)
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_params.spatialdb (rev 0)
+++ short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_params.spatialdb 2010-01-15 21:10:53 UTC (rev 16134)
@@ -0,0 +1,37 @@
+// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting)
+//
+// This spatial database defines the power-law parameters for two materials:
+// Wet granite and dry olivine. The properties are from Strehlau and
+// Meissner (1987).
+// Note that the flow constant is expressed as log10(flow-constant), where
+// the units are GPa^(-n)/s. The activation energy (Q) is expressed in
+// kJ/mol. To accommodate the strange units of the flow constant, the
+// given units of are 'None', and we specify a new parameter
+// (flow-constant-multiplier), with a value of 9, to indicate that the
+// underlying units are GPa. A value of 6 would indicate MPa, etc.
+//
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 4 // number of material property values
+ value-names = flow-constant activation-energy power-law-exponent flow-constant-multiplier // names of the material property values
+ value-units = None kJ/mol None None // units
+ num-locs = 4 // number of locations
+ data-dim = 1
+ space-dim = 3
+ cs-data = cartesian {
+ to-meters = 1000.0
+ space-dim = 3
+ }
+}
+// Columns are
+// (1) x coordinate (m)
+// (2) y coordinate (m)
+// (3) z coordinate (m)
+// (4) flow-constant (None)
+// (5) activation-energy (kJ/mol)
+// (6) power-law-exponent (None)
+// (7) activation-energy-multiplier (None)
+0.0 0.0 0.0 2.0 137.0 1.5 9.0
+0.0 0.0 -24.9 2.0 137.0 1.5 9.0
+0.0 0.0 -25.1 15.5 535.0 3.5 9.0
+0.0 0.0 -400. 15.5 535.0 3.5 9.0
Added: short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_points.txt
===================================================================
--- short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_points.txt (rev 0)
+++ short/3D/PyLith/trunk/playpen/powerlaw/powerlaw_points.txt 2010-01-15 21:10:53 UTC (rev 16134)
@@ -0,0 +1,2 @@
+0.0 0.0 0.0
+0.0 0.0 -1.0e+3
More information about the CIG-COMMITS
mailing list