[cig-commits] r4451 - in cs/benchmark/trunk: . benchmark
benchmark/discretize benchmark/short
benchmark/short/applications benchmark/short/examples
baagaard at geodynamics.org
baagaard at geodynamics.org
Tue Aug 29 16:54:11 PDT 2006
Author: baagaard
Date: 2006-08-29 16:54:11 -0700 (Tue, 29 Aug 2006)
New Revision: 4451
Added:
cs/benchmark/trunk/Makefile.am
cs/benchmark/trunk/README
cs/benchmark/trunk/TODO
cs/benchmark/trunk/benchmark/
cs/benchmark/trunk/benchmark/CompareApp.py
cs/benchmark/trunk/benchmark/Makefile.am
cs/benchmark/trunk/benchmark/ProjectionQuadPts.py
cs/benchmark/trunk/benchmark/__init__.py
cs/benchmark/trunk/benchmark/discretize/
cs/benchmark/trunk/benchmark/discretize/Tet4.py
cs/benchmark/trunk/benchmark/discretize/__init__.py
cs/benchmark/trunk/benchmark/short/
cs/benchmark/trunk/benchmark/short/applications/
cs/benchmark/trunk/benchmark/short/applications/Makefile.am
cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py
cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py
cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py
cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py
cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py
cs/benchmark/trunk/benchmark/short/examples/
cs/benchmark/trunk/benchmark/short/examples/compare_solns.py
cs/benchmark/trunk/benchmark/short/examples/mesh_catalog.py
cs/benchmark/trunk/benchmark/tablewaiter.py
cs/benchmark/trunk/configure.ac
cs/benchmark/trunk/tidy.sh
Log:
Factored tasks into different scripts and HDF files. Pyrized code.
Added: cs/benchmark/trunk/Makefile.am
===================================================================
--- cs/benchmark/trunk/Makefile.am 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/Makefile.am 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,16 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+SUBDIRS = \
+ benchmark
+
+# End of file
Added: cs/benchmark/trunk/README
===================================================================
--- cs/benchmark/trunk/README 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/README 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,22 @@
+0. Translate simulation output into HDF5 files.
+
+ For example:
+
+ pylithtoh5.py \
+ --mesh_name=tet4_1000m \
+ --scenario_root=output/pylith-0.8/tet4_1000m
+
+1. Construct workspace for playing with projected data.
+
+ For example:
+
+ create_projection.py \
+ --solution=output/pylith-0.8/tet4_1000m.h5 \
+ --dump_filename=output/analytic/tet4_1000m_points.in
+
+2. Compute analytical solution at projection quadrature points and
+ translate solution into HDF5 file.
+
+3. Load solutions into projection workspace and compute differences.
+
+4. Plot results
Added: cs/benchmark/trunk/TODO
===================================================================
--- cs/benchmark/trunk/TODO 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/TODO 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,8 @@
+Vectorize as many loops over elements or nodes as possible to improve
+speed.
+
+Create projection object that will actually project solution
+to points using efficient search algorithm.
+
+Pyrize discretize objects so it is easy to specify reference element
+basis functions and quadrature points.
Added: cs/benchmark/trunk/benchmark/CompareApp.py
===================================================================
--- cs/benchmark/trunk/benchmark/CompareApp.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/CompareApp.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,83 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file benchmark/CompareApp.py
+
+## @brief Python application to compare projected solutions.
+
+from pyre.applications.Script import Script
+
+# CompareApp class
+class CompareApp(Script):
+ """
+ Python application to compare projected solutions.
+ """
+
+ # INVENTORY ////////////////////////////////////////////////////////
+
+ class Inventory(Script.Inventory):
+ """
+ Python object for managing CompareApp facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing CompareApp facilities and properties.
+ ##
+ ## \b Properties
+ ## @li none
+ ##
+ ## \b Facilities
+ ## @li \b projection Projection to use in projecting solutions
+
+ import pyre.inventory
+
+ from ProjectionQuadPts import ProjectionQuadPts
+ projection = pyre.inventory.facility("projection",
+ factory=ProjectionQuadPts)
+ projection.meta['tip'] = "Projection to use in projecting solutions."
+
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="compareapp"):
+ """
+ Constructor.
+ """
+ Script.__init__(self, name)
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Entry point for application.
+
+ Don't know what user wants to do, so leave this empty and let
+ them do whatever they want by adding stuff to a driver after a
+ call to run(). By calling run(), Pythia handles ingestion of
+ all of the parameters.
+ """
+ return
+
+
+ # PRIVATE METHODS //////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based on inventory.
+ """
+ Script._configure(self)
+
+ self.projection = self.inventory.projection
+ return
+
+
+# End of file
Added: cs/benchmark/trunk/benchmark/Makefile.am
===================================================================
--- cs/benchmark/trunk/benchmark/Makefile.am 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/Makefile.am 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,31 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+# Applications
+dist_bin_SCRIPTS = \
+ short/applications/bm_create_projection.py \
+ short/applications/femlabtoh5.py \
+ short/applications/geofesttoh5.py \
+ short/applications/okadaelastictoh5.py \
+ short/applications/pylithtoh5.py
+
+# Components
+nobase_pkgpyexec_PYTHON = \
+ CompareApp.py \
+ ProjectionQuadPts.py \
+ tablewaiter.py \
+ __init__.py \
+ discretize/Tet4.py \
+ discretize/__init__.py
+
+
+# End of file
Added: cs/benchmark/trunk/benchmark/ProjectionQuadPts.py
===================================================================
--- cs/benchmark/trunk/benchmark/ProjectionQuadPts.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/ProjectionQuadPts.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,334 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file ProjectionQuadPts.py
+
+## @brief Python object for projecting solutions to element quadrature
+## points. Projected solutions are stored in an HDF5 file.
+
+from pyre.components.Component import Component
+
+import tables
+import tablewaiter
+
+# ProjectionQuadPts class
+class ProjectionQuadPts(Component):
+ """
+ Python object for projecting solutions to element quadrature
+ points. Projected solutions are stored in an HDF5 file.
+ """
+
+ class Inventory(Component.Inventory):
+ """
+ Python object for managing ProjectionQuadPts facilities and
+ properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing ProjectionQuadPts facilities and
+ ## properties.
+ ##
+ ## \b Properties
+ ## @li \b filename Name of filename for workspace
+ ## @li \b quad_order Order of quadrature to use
+ ##
+ ## \b Facilities
+ ## @li \b discretization Interpolation and quadrature information
+
+ import pyre.inventory
+
+ filename = pyre.inventory.str("filename", default="quadpts_error.h5")
+ filename.meta['tip'] = "Name of filename for workspace."
+
+ quadratureOrder = pyre.inventory.int("quad_order", default=1)
+ quadratureOrder.meta['tip'] = "Order of quadrature to use."
+
+ from discretize.Tet4 import Tet4
+ discretization = pyre.inventory.facility("discretization",
+ factory=Tet4)
+ discretization.meta['tip'] = "Interpolation and quadrature " \
+ "information."
+
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="projectionquadpts"):
+ """
+ Constructor.
+ """
+ Component.__init__(self, name, facility="projection")
+ return
+
+
+ def open(self, mode='a'):
+ """
+ Open workspace for use.
+ """
+ self.workspace = tables.openFile(self.filename, mode)
+ return
+
+
+ def close(self):
+ """
+ Close workspace.
+ """
+ self.workspace.close()
+ return
+
+
+ def create(self, solution):
+ """
+ Create workspace (HDF5 file) to hold projection
+ information. Workspace must be already open.
+ """
+
+ self._copySolnInfo(solution)
+ self._calcProjectPts()
+ return
+
+
+ def dumpProjectPts(self, filename):
+ """
+ Dump projection quadrature points to ASCII file. Workspace
+ must already be open.
+ """
+
+ points = self.workspace.root.projection.points[:]
+
+ nelems = points.shape[0]
+ nquadpts = points.shape[1]
+ fileout = open(filename, 'w')
+ for ielem in xrange(nelems):
+ for iquad in xrange(nquadpts):
+ pt = tuple(points[ielem, iquad])
+ fileout.write("%14.6e %14.6e %14.6e\n" % pt)
+ fileout.close()
+ return
+
+
+ def project(self, filename, label, timestep, field, overwrite=True):
+ """
+ Project solution in to quadrature points.
+
+ @param filename Name of solution file
+ @param label Label for solution
+ @param timestep Index of timestep to project
+ @param field Name of field to project
+ """
+
+ from numpy import array, zeros, dot
+
+ # Get field from solution file
+ solnfile = tables.openFile(filename, 'r')
+ solnfield = tablewaiter.getNode(solnfile.root,
+ "solution/snapshot%d/%s" % \
+ (timestep, field))[:]
+ solnfile.close()
+
+ # Get coordinates of quadrature points in reference element
+ quadpts = self.workspace.root.interpolation.quadrature.points[:]
+ nquadpts = quadpts.shape[0]
+ ndims = quadpts.shape[1]
+ quadcoords = [map(float, tuple(quadpts[q])) \
+ for q in xrange(nquadpts)]
+
+ # Get basis functions
+ basisfns = self.discretization.getBasisFns()
+
+ # Tabulate basis fns at quadrature points
+ basisfns_quad = array(basisfns.tabulate(quadcoords)).transpose()
+
+ simplices = self.workspace.root.topology.simplices[:]
+ nelems = simplices.shape[0]
+ connSize = simplices.shape[1]
+
+ fieldProj = zeros((nelems, nquadpts, ndims), dtype=float)
+
+ fieldLocal = zeros((connSize, ndims), dtype=float)
+ for ielem in xrange(nelems):
+ nodeIndices = simplices[ielem]
+ fieldLocal[:] = solnfield[nodeIndices,:]
+ fieldProj[ielem] = dot(basisfns_quad, fieldLocal)
+
+ group = tablewaiter.createGroup(self.workspace,
+ self.workspace.root,
+ "projection/data/%s/snapshot%d" % \
+ (label, timestep))
+ if field in group and overwrite:
+ dataset = getattr(group, field)
+ dataset._f_remove()
+ self.workspace.createArray(group, field, fieldProj)
+ return
+
+
+ def copy_projection(self, filename, label, timestep, field,
+ overwrite=True):
+ """
+ Copy precomputed projection from file.
+ """
+ from numpy import dot
+
+ filein = tables.openFile(filename, 'r')
+
+ # Check to make sure projection is same as local one
+ pointsExt = filein.root.projection.points[:]
+ points = self.workspace.root.projection.points[:]
+ dd = pointsExt.flatten() - points.flatten()
+ tolerance = 1.0e-4 * len(dd)
+ if dot(dd, dd) > tolerance:
+ raise ValueError("Projection in file '%s' does not match "
+ "projection in workspace." % filename)
+
+ # Copy field evaluated at projection points to workspace
+ fieldExt = tablewaiter.getNode(filein.root,
+ "projection/data/snapshot%d/%s" % \
+ (timestep, field))[:]
+ filein.close()
+
+ group = tablewaiter.createGroup(self.workspace,
+ self.workspace.root,
+ "projection/data/%s/snapshot%d" % \
+ (label, timestep))
+ if field in group and overwrite:
+ dataset = getattr(group, field)
+ dataset._f_remove()
+ self.workspace.createArray(group, field, fieldExt)
+ return
+
+
+ def difference(self, timestep, field, labelA, labelB, overwrite=True):
+ """
+ Compute difference in field 'field' at time step 'timestep'
+ for data labelA and labelB (labelA-labelB). Integrate difference to get
+ 'error'.
+ """
+
+ from numpy import zeros, dot, sqrt, reshape
+
+ dataA = tablewaiter.getNode(self.workspace.root,
+ "projection/data/%s/snapshot%d/%s" % \
+ (labelA, timestep, field))[:]
+ dataB = tablewaiter.getNode(self.workspace.root,
+ "projection/data/%s/snapshot%d/%s" % \
+ (labelB, timestep, field))[:]
+ quadwts = self.workspace.root.interpolation.quadrature.weights[:]
+ nquadpts = quadwts.shape[0]
+ quadwts = reshape(quadwts, (1,nquadpts))
+ jacobian = self.workspace.root.interpolation.map.jacobian
+ nelems = self.workspace.root._v_attrs.num_simplices
+ ndims = self.workspace.root._v_attrs.spatial_dim
+ localerr = zeros( (nelems,), dtype=float)
+ for ielem in xrange(nelems):
+ diff = reshape(dataA[ielem] - dataB[ielem], (nquadpts, ndims))
+ localerr[ielem] = jacobian[ielem] * sum(sum(dot(quadwts, diff**2)))
+ group = tablewaiter.createGroup(self.workspace,
+ self.workspace.root,
+ "difference/%s_%s/snapshot%d" % \
+ (labelA, labelB, timestep))
+
+ if field in group and overwrite:
+ dataset = getattr(group, field)
+ dataset._f_remove()
+ self.workspace.createArray(group, field, localerr)
+ return sqrt(sum(localerr))
+
+
+ def remove(self, name, recursive=True):
+ """
+ Remove object 'name' from workspace.
+ """
+
+ node = tablewaiter.getNode(self.workspace.root, name)
+ node._f_remove(recursive=recursive)
+
+ return
+
+
+ # PRIVATE METHODS //////////////////////////////////////////////////
+
+ def _copySolnInfo(self, solution):
+ """
+ Copy topology from solution file to workspace.
+ """
+ self._info.log("Copying solution info...")
+
+ solnfile = tables.openFile(solution, 'r')
+ workspace = self.workspace
+
+ solnattrs = solnfile.root._v_attrs
+ workspace.root._v_attrs.spatial_dim = solnattrs.spatial_dim
+ workspace.root._v_attrs.num_vertices = solnattrs.num_vertices
+ workspace.root._v_attrs.num_simplices = solnattrs.num_simplices
+ workspace.root._v_attrs.simplex_size = solnattrs.simplex_size
+
+ solnfile.root.topology._f_copy(workspace.root, recursive=True)
+ solnfile.root.geometry._f_copy(workspace.root, recursive=True)
+
+ tablewaiter.createGroup(workspace, workspace.root, "interpolation")
+ solnfile.root.interpolation.map._f_copy(workspace.root.interpolation,
+ recursive=True)
+ solnfile.close()
+ return
+
+
+ def _calcProjectPts(self):
+ """
+ Calculate coordinates of points in projection.
+ """
+ self._info.log("Calculating quadrature points for projection...")
+
+ workspace = self.workspace
+
+ # Get integration quadrature points.
+ (quadpts, quadwts) = \
+ self.discretization.getQuadPts(self.quadratureOrder)
+ group = tablewaiter.createGroup(workspace, workspace.root,
+ "interpolation/quadrature")
+ workspace.createArray(group, "points", quadpts)
+ workspace.createArray(group, "weights", quadwts)
+
+ # Compute coordinates of quadrature points for each element
+ from numpy import dot, zeros
+
+ mapOrigin = workspace.root.interpolation.map.origin[:]
+ mapFromRef = workspace.root.interpolation.map.from_ref[:]
+
+ nelems = workspace.root._v_attrs.num_simplices
+ nquadpts = len(quadpts)
+ ndims = len(quadpts[0])
+
+ points = zeros((nelems, nquadpts, ndims), dtype=float)
+ for ielem in xrange(nelems):
+ j = mapFromRef[ielem]
+ origin = mapOrigin[ielem]
+ for iquad in xrange(nquadpts): # TODO: vectorize this loop
+ points[ielem, iquad] = origin + dot(j, quadpts[iquad])
+
+ group = tablewaiter.createGroup(workspace, workspace.root,
+ "projection")
+ workspace.createArray(group, "points", points)
+ return
+
+
+ def _configure(self):
+ """
+ Set members based on inventory.
+ """
+ Component._configure(self)
+
+ self.filename = self.inventory.filename
+ self.quadratureOrder = self.inventory.quadratureOrder
+ self.discretization = self.inventory.discretization
+ return
+
+
+# End of file
Added: cs/benchmark/trunk/benchmark/__init__.py
===================================================================
--- cs/benchmark/trunk/benchmark/__init__.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/__init__.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,28 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file benchmark/__init__.py
+
+## @brief Python top-level benchmark module initialization
+
+
+def copyright():
+ return "benchmark pyre module: Copyright (c) 2006 Computational " \
+ "Infrastructure for Geodynamics"
+
+
+__all__ = ['CompareApp',
+ 'ProjectionQuadPts',
+ 'tablewaiter']
+
+
+# End of file
Added: cs/benchmark/trunk/benchmark/discretize/Tet4.py
===================================================================
--- cs/benchmark/trunk/benchmark/discretize/Tet4.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/discretize/Tet4.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,100 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file benchmark/discretize/Tet4.py
+
+## @brief Python object describing interpolation and quadrature
+## information for linear tetrahedral finite-element.
+
+from pyre.components.Component import Component
+
+# Tet4 class
+class Tet4(Component):
+ """
+ Python object describing interpolation and quadrature information
+ for linear tetrahedral finite-element.
+ """
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="elemtet4"):
+ """
+ Constructor.
+ """
+ Component.__init__(self, name, facility="discretization")
+ return
+
+
+ def calcMap(self, coordinates, simplices):
+ """
+ Calculate mapping to/from reference element.
+ """
+ from numpy.linalg import inv, det
+ from numpy import zeros, array
+
+ nvertices = coordinates.shape[0]
+ ndims = coordinates.shape[1]
+ nelems = simplices.shape[0]
+ jacobian = zeros((nelems,), dtype=float)
+ toRef = zeros((nelems, ndims, ndims), dtype=float)
+ fromRef = zeros((nelems, ndims, ndims), dtype=float)
+ origin = zeros((nelems, 3), dtype=float)
+
+ for e in xrange(nelems):
+ # print "processing element", e
+ a,b,c,d = simplices[e]
+ x0,y0,z0 = coordinates[a] # a = (x0,y0,z0)
+ x1,y1,z1 = coordinates[b] # b = (x1,y1,z1)
+ x2,y2,z2 = coordinates[c] # c = (x2,y2,z2)
+ x3,y3,z3 = coordinates[d] # d = (x3,y3,z3)
+ j = array([[(x1 - x0)/2, (x2 - x0)/2, (x3 - x0)/2],
+ [(y1 - y0)/2, (y2 - y0)/2, (y3 - y0)/2],
+ [(z1 - z0)/2, (z2 - z0)/2, (z3 - z0)/2]], dtype=float)
+ jacobian[e] = det(j)
+ fromRef[e] = j
+ toRef[e] = inv(j)
+ origin[e] = ((x1+x2+x3-x0)/2, (y1+y2+y3-y0)/2, (z1+z2+z3-z0)/2)
+
+ map = {'jacobian': jacobian,
+ 'toRef': toRef,
+ 'fromRef': fromRef,
+ 'origin': origin}
+ return map
+
+
+ def getQuadPts(self, order):
+ """
+ Get quadrature points for element given quadrature order.
+ """
+ from FIAT.shapes import TETRAHEDRON
+ from FIAT.quadrature import make_quadrature
+ from numpy import array
+
+ quadrature = make_quadrature(TETRAHEDRON, order)
+ pts = array(quadrature.x)
+ wts = array(quadrature.w)
+ return (pts, wts)
+
+
+ def getBasisFns(self):
+ """
+ Get basis functions.
+ """
+ from FIAT.Lagrange import Lagrange
+ from FIAT.shapes import TETRAHEDRON
+
+ degree = 1
+ basis = Lagrange(TETRAHEDRON, degree)
+ return basis.function_space()
+
+
+# End of file
Added: cs/benchmark/trunk/benchmark/discretize/__init__.py
===================================================================
--- cs/benchmark/trunk/benchmark/discretize/__init__.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/discretize/__init__.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file benchmark/discreteize/__init__.py
+
+## @brief Python benchmark discretization module initialization
+
+
+def copyright():
+ return "benchmark discretization pyre module: " \
+ "Copyright (c) 2006 Computational Infrastructure for Geodynamics"
+
+
+__all__ = ['Tet4']
+
+
+# End of file
Added: cs/benchmark/trunk/benchmark/short/applications/Makefile.am
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/Makefile.am 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/applications/Makefile.am 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,23 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+dist_bin_SCRIPTS = \
+ bm_create_projection.py \
+ femlabtoh5.py \
+ geofesttoh5.py \
+ okadaelastictoh5.py \
+ pylithtoh5.py
+
+# version
+# $Id$
+
+# End of file
Added: cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file benchmark/short/applications/bm_create_projection.py
+
+## @brief Python application to create projection and dump projection
+## quadrature points to file.
+
+from pyre.applications.Script import Script
+
+# Creator class
+class Creator(Script):
+ """
+ Python application to create projection and dump projection
+ quadrature points to file.
+ """
+
+ # INVENTORY ////////////////////////////////////////////////////////
+
+ class Inventory(Script.Inventory):
+ """
+ Python object for managing Creator facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing Creator facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b solution Name of solution file to use for projection
+ ## @li \b dumpFilename Filename for quadrature points dump.
+ ##
+ ## \b Facilities
+ ## @li \b projection Projection to use
+
+ import pyre.inventory
+
+ solution = pyre.inventory.str("solution", default="")
+ solution.meta['tip'] = "Name of solution file to use for projection."
+
+ dumpFilename = pyre.inventory.str("dump_filename",
+ default="qpoints.in")
+ dumpFilename.meta['tip'] = "Filename for quadrature points dump."
+
+ from benchmark.ProjectionQuadPts import ProjectionQuadPts
+ projection = pyre.inventory.facility("projection",
+ factory=ProjectionQuadPts)
+ projection.meta['tip'] = "Projection to use."
+
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="bm_create_projection"):
+ """
+ Constructor.
+ """
+ Script.__init__(self, name)
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Entry point for application.
+ """
+ self.projection.open('w')
+ self.projection.create(self.solution)
+ self.projection.dumpProjectPts(self.dumpFilename)
+ self.projection.close()
+ return
+
+
+ # PRIVATE METHODS //////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based on inventory.
+ """
+ Script._configure(self)
+
+ if self.inventory.solution == "":
+ raise ValueError("Property 'solution' must be supplied.")
+ self.solution = self.inventory.solution
+ self.dumpFilename = self.inventory.dumpFilename
+ self.projection = self.inventory.projection
+ return
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+ app = Creator()
+ app.run()
+
+# End of file
Property changes on: cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,126 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file benchmark/short/applications/femlabtoh5.py
+
+## @brief Python application to translate Femlab solution computed
+## at projection quadrature points to HDF5 file.
+
+from pyre.applications.Script import Script
+from numpy import zeros
+from benchmark.tablewaiter import createGroup
+
+# Translator class
+class Translator(Script):
+ """
+ Python application to translate Femlab solution computed at
+ projection quadrature points to HDF5 file.
+ """
+
+ # INVENTORY ////////////////////////////////////////////////////////
+
+ class Inventory(Script.Inventory):
+ """
+ Python object for managing Translator facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing Translator facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b ascii_file Name of file containing analytical solution.
+ ## @li \b h5_file Name of HDF5 file to create.
+ ##
+ ## \b Facilities
+ ## @li None
+
+ import pyre.inventory
+
+ asciiFile = pyre.inventory.str("ascii_file", default="points.disp")
+ asciiFile.meta['tip'] = "Name of file containing analytical solution."
+
+ h5File = pyre.inventory.str("h5_file", default="points.h5")
+ h5File.meta['tip'] = "Name of HDF5 file to create."
+
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="femlabtoh5"):
+ Script.__init__(self, name)
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Entry point for application.
+ """
+
+ npts = self._countLines(self.asciiFile)
+ ndims = 3
+
+ points = zeros( (npts, ndims), dtype=float)
+ displacements = zeros( (npts, ndims), dtype=float)
+
+ filein = open(self.asciiFile, 'r')
+ ipt = 0
+ for line in filein:
+ values = map(float, line.split())
+ points[ipt] = tuple(values[0:3])
+ displacements[ipt] = tuple(values[3:6])
+ ipt += 1
+ filein.close()
+
+ import tables
+ fileOut = tables.openFile(self.h5File, 'w')
+
+ group = createGroup(fileOut, fileOut.root, "projection")
+ fileOut.createArray(group, "points", points)
+
+ group = createGroup(fileOut, fileOut.root,
+ "projection/data/snapshot0")
+ fileOut.createArray(group, "displacements", displacements)
+ fileOut.close()
+ return
+
+
+ # PRIVATE METHODS //////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based on inventory.
+ """
+ Script._configure(self)
+
+ self.asciiFile = self.inventory.asciiFile
+ self.h5File = self.inventory.h5File
+ return
+
+
+ def _countLines(self, filename):
+ """
+ Count number of lines in file.
+ """
+
+ filein = open(filename, 'r')
+ count = 0
+ for line in filein:
+ count += 1
+ filein.close()
+ return count
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+ app = Translator()
+ app.run()
+
+# End of file
Property changes on: cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,245 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file benchmark/short/applications/geofesttoh5.py
+
+## @brief Python application to translate GeoFEST output from UCD
+## format to HDF5.
+
+from pyre.applications.Script import Script
+from numpy import zeros
+
+# Translator class
+class Translator(Script):
+ """
+ Python application to translate GeoFEST output from UCD format to HDF5.
+ """
+
+ # INVENTORY ////////////////////////////////////////////////////////
+
+ class Inventory(Script.Inventory):
+ """
+ Python object for managing Translator facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing Translator facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b mesh_name Name of finite-element mesh
+ ## @li \b scenarioRoot Filename root for simulation scenario
+ ## @li \b timeSteps Indices of time steps
+ ##
+ ## \b Facilities
+ ## @li \b discretization Interpolation and quadrature information
+
+ import pyre.inventory
+
+ meshName = pyre.inventory.str("mesh_name", default="")
+ meshName.meta['tip'] = "Name of finite-element mesh."
+
+ scenarioRoot = pyre.inventory.str("scenario_root",
+ default="output/geofest")
+ scenarioRoot.meta['tip'] = "Filename root for simulation scenario."
+
+ timeSteps = pyre.inventory.list("time_steps", default=[0,10,50,100])
+ timeSteps.meta['tip'] = "Indices of time steps."
+
+ from benchmark.discretize.Tet4 import Tet4
+ discretization = pyre.inventory.facility("discretization",
+ factory=Tet4)
+ discretization.meta['tip'] = \
+ "Interpolation and quadrature information."
+
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="geofesttoh5"):
+ Script.__init__(self, name)
+ return
+
+
+ def main(self, *args, **kwds):
+ self._initialize()
+ self._loadGeometry()
+ self._loadTopology()
+ self._loadSolution()
+ self._loadInterpolation()
+ self.h5.close()
+ return
+
+
+ # PRIVATE METHODS //////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based on inventory.
+ """
+ Script._configure(self)
+
+ if self.inventory.meshName == "":
+ raise ValueError("Property mesh_name must be supplied.")
+ self.meshName = self.inventory.meshName
+ self.scenarioRoot = self.inventory.scenarioRoot
+ self.timeSteps = self.inventory.timeSteps
+ self.discretization = self.inventory.discretization
+ return
+
+
+ def _initialize(self):
+ """
+ Prepare for translation.
+ """
+ self._info.log("Initializing translator...")
+
+ import mesh_catalog
+ mesh = mesh_catalog.catalog[self.meshName]
+ self.spatialDim = mesh['spatialDim']
+ self.numNodes = mesh['numNodes']
+ self.numElements = mesh['numElements']
+ self.connSize = mesh['connSize']
+
+ import tables
+ filename = "%s.h5" % self.scenarioRoot
+ self.h5 = tables.openFile(filename, 'w')
+ self.h5.root._v_attrs.spatial_dim = self.spatialDim
+ self.h5.root._v_attrs.num_vertices = self.numNodes
+ self.h5.root._v_attrs.num_simplices = self.numElements
+ self.h5.root._v_attrs.simplex_size = self.connSize
+ return
+
+
+ def _loadGeometry(self):
+ """
+ Load mesh geometry.
+ """
+ self._info.log("Loading geometry...")
+
+ coordinates = zeros((self.numNodes, self.spatialDim), dtype=float)
+
+ filename = "%s.coord" % self.scenarioRoot
+ filein = open(filename, 'r')
+ scale = 1.0
+ for line in filein:
+ if not "#" == line[0]:
+ fields = line.split()
+ if not "coord_units" == fields[0]:
+ label = int(fields[0])
+ x,y,z = map(float, (fields[1], fields[2], fields[3]))
+ coordinates[label-1] = (scale*x, scale*y, scale*z)
+ else:
+ import pyre.units
+ parser = pyre.units.parser()
+ scale = parser.parse(fields[2]).value
+ filein.close()
+
+ self.h5.createGroup("/", "geometry")
+ self.h5.createArray("/geometry", "coordinates", coordinates)
+ return
+
+
+ def _loadTopology(self):
+ """
+ Load mesh topology.
+ """
+ self._info.log("Loading topology...")
+
+ simplices = zeros((self.numElements, self.connSize), dtype=int)
+
+ filename = "%s.connect" % self.scenarioRoot
+ filein = open(filename, 'r')
+ for line in filein:
+ if not "#" == line[0]:
+ fields = line.split()
+ label = int(fields[0])
+ simplex = map(int, fields[4:5+self.connSize])
+ # Change simplices from node labels to node indices
+ for i in xrange(self.connSize):
+ simplex[i] -= 1
+ simplices[label-1] = tuple(simplex)
+ filein.close()
+ self.h5.createGroup("/", "topology")
+ self.h5.createArray("/topology", "simplices", simplices)
+ return
+
+
+ def _loadSolution(self):
+ """
+ Load displacements and velocities.
+ """
+ self._info.log("Loading solution...")
+
+ self.h5.createGroup("/", "solution")
+
+ dispStart = 5
+ dispEnd = dispStart + self.spatialDim
+ velStart = 8
+ velEnd = velStart + self.spatialDim
+ filename = "%sA.out" % self.scenarioRoot
+
+ displacements = zeros((self.numNodes, self.spatialDim), dtype=float)
+ velocities = zeros((self.numNodes, self.spatialDim), dtype=float)
+ filein = open(filename, 'r')
+ numHeaderLines = 4
+ for iskip in xrange(numHeaderLines):
+ line = filein.readline()
+ for istep in self.timeSteps:
+ # Read node information
+ numHeaderLines = 5
+ for iskip in xrange(numHeaderLines):
+ line = filein.readline()
+ for i in xrange(self.numNodes):
+ fields = filein.readline().split()
+ label = int(fields[1])
+ displacements[label-1] = tuple(map(float,
+ fields[dispStart:dispEnd]))
+ velocities[label-1] = tuple(map(float,
+ fields[velStart:velEnd]))
+ self.h5.createGroup("/solution", "snapshot%d" % istep)
+ self.h5.createArray("/solution/snapshot%d" % istep,
+ "displacements", displacements)
+ self.h5.createArray("/solution/snapshot%d" % istep,
+ "velocities", velocities)
+
+ # Read element information (header only)
+ numHeaderLines = 8
+ for iskip in xrange(numHeaderLines):
+ line = filein.readline()
+ filein.close()
+ return
+
+
+ def _loadInterpolation(self):
+ """
+ Load interpolation and quadrature information.
+ """
+ self._info.log("Loading interpolation mapping...")
+
+ from numpy.linalg import inv, det
+ coordinates = self.h5.root.geometry.coordinates[:]
+ simplices = self.h5.root.topology.simplices[:]
+ map = self.discretization.calcMap(coordinates, simplices)
+ self.h5.createGroup("/", "interpolation")
+ self.h5.createGroup("/interpolation", "map")
+ self.h5.createArray("/interpolation/map", "jacobian", map['jacobian'])
+ self.h5.createArray("/interpolation/map", "to_ref", map['toRef'])
+ self.h5.createArray("/interpolation/map", "from_ref", map['fromRef'])
+ self.h5.createArray("/interpolation/map", "origin", map['origin'])
+ return
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+ app = Translator()
+ app.run()
+
+# End of file
Property changes on: cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,126 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file benchmark/short/applications/okadaelastictoh5.py
+
+## @brief Python application to translate analytical elastic solution computed
+## using Okada at projection quadrature points to HDF5 file.
+
+from pyre.applications.Script import Script
+from numpy import zeros
+from benchmark.tablewaiter import createGroup
+
+# Translator class
+class Translator(Script):
+ """
+ Python application to translate analytical elastic solution
+ computed using Okada at projection quadrature points to HDF5 file.
+ """
+
+ # INVENTORY ////////////////////////////////////////////////////////
+
+ class Inventory(Script.Inventory):
+ """
+ Python object for managing Translator facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing Translator facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b ascii_file Name of file containing analytical solution.
+ ## @li \b h5_file Name of HDF5 file to create.
+ ##
+ ## \b Facilities
+ ## @li None
+
+ import pyre.inventory
+
+ asciiFile = pyre.inventory.str("ascii_file", default="points.disp")
+ asciiFile.meta['tip'] = "Name of file containing analytical solution."
+
+ h5File = pyre.inventory.str("h5_file", default="points.h5")
+ h5File.meta['tip'] = "Name of HDF5 file to create."
+
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="okadaelastictoh5"):
+ Script.__init__(self, name)
+ return
+
+
+ def main(self, *args, **kwds):
+ """
+ Entry point for application.
+ """
+
+ npts = self._countLines(self.asciiFile)
+ ndims = 3
+
+ points = zeros( (npts, ndims), dtype=float)
+ displacements = zeros( (npts, ndims), dtype=float)
+
+ filein = open(self.asciiFile, 'r')
+ ipt = 0
+ for line in filein:
+ values = map(float, line.split())
+ points[ipt] = tuple(values[0:3])
+ displacements[ipt] = tuple(values[3:6])
+ ipt += 1
+ filein.close()
+
+ import tables
+ fileOut = tables.openFile(self.h5File, 'w')
+
+ group = createGroup(fileOut, fileOut.root, "projection")
+ fileOut.createArray(group, "points", points)
+
+ group = createGroup(fileOut, fileOut.root,
+ "projection/data/snapshot0")
+ fileOut.createArray(group, "displacements", displacements)
+ fileOut.close()
+ return
+
+
+ # PRIVATE METHODS //////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based on inventory.
+ """
+ Script._configure(self)
+
+ self.asciiFile = self.inventory.asciiFile
+ self.h5File = self.inventory.h5File
+ return
+
+
+ def _countLines(self, filename):
+ """
+ Count number of lines in file.
+ """
+
+ filein = open(filename, 'r')
+ count = 0
+ for line in filein:
+ count += 1
+ filein.close()
+ return count
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+ app = Translator()
+ app.run()
+
+# End of file
Property changes on: cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,249 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file benchmark/short/applications/pylithtoh5.py
+
+## @brief Python application to translate PyLith output from UCD
+## format to HDF5.
+
+from pyre.applications.Script import Script
+from numpy import zeros
+
+# Translator class
+class Translator(Script):
+ """
+ Python application to translate PyLith output from UCD format to HDF5.
+ """
+
+ # INVENTORY ////////////////////////////////////////////////////////
+
+ class Inventory(Script.Inventory):
+ """
+ Python object for managing Translator facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing Translator facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b mesh_name Name of finite-element mesh
+ ## @li \b scenarioRoot Filename root for simulation scenario
+ ## @li \b numProcs Number of procesors used in simulation
+ ## @li \b timeSteps Indices of time steps
+ ##
+ ## \b Facilities
+ ## @li \b discretization Interpolation and quadrature information
+
+ import pyre.inventory
+
+ meshName = pyre.inventory.str("mesh_name", default="")
+ meshName.meta['tip'] = "Name of finite-element mesh."
+
+ scenarioRoot = pyre.inventory.str("scenario_root",
+ default="output/pylith-0.8")
+ scenarioRoot.meta['tip'] = "Filename root for simulation scenario."
+
+ numProcs = pyre.inventory.int("num_procs", default=1)
+ numProcs.meta['tip'] = "Number of processors used in simulation."
+
+ timeSteps = pyre.inventory.list("time_steps", default=[0,10,50,100])
+ timeSteps.meta['tip'] = "Indices of time steps."
+
+ from benchmark.discretize.Tet4 import Tet4
+ discretization = pyre.inventory.facility("discretization",
+ factory=Tet4)
+ discretization.meta['tip'] = \
+ "Interpolation and quadrature information."
+
+
+ # PUBLIC METHODS ///////////////////////////////////////////////////
+
+ def __init__(self, name="pylithtoh5"):
+ Script.__init__(self, name)
+ return
+
+
+ def main(self, *args, **kwds):
+ self._initialize()
+ self._loadGeometry()
+ self._loadTopology()
+ self._loadSolution()
+ self._loadInterpolation()
+ self.h5.close()
+ return
+
+
+ # PRIVATE METHODS //////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based on inventory.
+ """
+ Script._configure(self)
+
+ if self.inventory.meshName == "":
+ raise ValueError("Property mesh_name must be supplied.")
+ self.meshName = self.inventory.meshName
+ self.scenarioRoot = self.inventory.scenarioRoot
+ self.numProcs = self.inventory.numProcs
+ self.timeSteps = self.inventory.timeSteps
+ self.discretization = self.inventory.discretization
+ return
+
+
+ def _initialize(self):
+ """
+ Prepare for translation.
+ """
+ self._info.log("Initializing translator...")
+
+ import mesh_catalog
+ mesh = mesh_catalog.catalog[self.meshName]
+ self.spatialDim = mesh['spatialDim']
+ self.numNodes = mesh['numNodes']
+ self.numElements = mesh['numElements']
+ self.connSize = mesh['connSize']
+
+ import tables
+ filename = "%s.h5" % self.scenarioRoot
+ self.h5 = tables.openFile(filename, 'w')
+ self.h5.root._v_attrs.spatial_dim = self.spatialDim
+ self.h5.root._v_attrs.num_vertices = self.numNodes
+ self.h5.root._v_attrs.num_simplices = self.numElements
+ self.h5.root._v_attrs.simplex_size = self.connSize
+ return
+
+
+ def _loadGeometry(self):
+ """
+ Load mesh geometry.
+ """
+ self._info.log("Loading geometry...")
+
+ coordinates = zeros((self.numNodes, self.spatialDim), dtype=float)
+
+ for iproc in xrange(self.numProcs):
+ filename = "%s_%d.%d.coord" % \
+ (self.scenarioRoot, self.numProcs, iproc)
+ filein = open(filename, 'r')
+ scale = 1.0
+ for line in filein:
+ if not "#" == line[0]:
+ fields = line.split()
+ if not "coord_units" == fields[0]:
+ label = int(fields[0])
+ x,y,z = map(float, (fields[1], fields[2], fields[3]))
+ coordinates[label-1] = (scale*x, scale*y, scale*z)
+ else:
+ import pyre.units
+ parser = pyre.units.parser()
+ scale = parser.parse(fields[2]).value
+ filein.close()
+
+ self.h5.createGroup("/", "geometry")
+ self.h5.createArray("/geometry", "coordinates", coordinates)
+ return
+
+
+ def _loadTopology(self):
+ """
+ Load mesh topology.
+ """
+ self._info.log("Loading topology...")
+
+ simplices = zeros((self.numElements, self.connSize), dtype=int)
+
+ for iproc in xrange(self.numProcs):
+ filename = "%s_%d.%d.connect" % \
+ (self.scenarioRoot, self.numProcs, iproc)
+ filein = open(filename, 'r')
+ for line in filein:
+ if not "#" == line[0]:
+ fields = line.split()
+ label = int(fields[0])
+ simplex = map(int, fields[4:5+self.connSize])
+ # Change simplices from node labels to node indices
+ for i in xrange(self.connSize):
+ simplex[i] -= 1
+ simplices[label-1] = tuple(simplex)
+ filein.close()
+ self.h5.createGroup("/", "topology")
+ self.h5.createArray("/topology", "simplices", simplices)
+ return
+
+
+ def _loadSolution(self):
+ """
+ Load displacements and velocities.
+ """
+ self._info.log("Loading solution...")
+
+ self.h5.createGroup("/", "solution")
+
+ dispStart = 1
+ dispEnd = dispStart + self.spatialDim
+ velStart = 4
+ velEnd = velStart + self.spatialDim
+ displacements = zeros((self.numNodes, self.spatialDim), dtype=float)
+ velocities = zeros((self.numNodes, self.spatialDim), dtype=float)
+ for istep in self.timeSteps:
+ displacements[:] = 0.0
+ velocities[:] = 0.0
+ for iproc in xrange(self.numProcs):
+ filename = "%s_%d.%d.mesh.time.%05d.inp" % \
+ (self.scenarioRoot, self.numProcs, iproc, istep)
+ filein = open(filename, 'r')
+ numHeaderLines = 7
+ for iskip in xrange(numHeaderLines):
+ line = filein.readline()
+ for i in xrange(self.numNodes):
+ fields = filein.readline().split()
+ label = int(fields[0])
+ displacements[label-1] = tuple(map(float,
+ fields[dispStart:dispEnd]))
+ velocities[label-1] = tuple(map(float,
+ fields[velStart:velEnd]))
+ filein.close()
+
+ self.h5.createGroup("/solution", "snapshot%d" % istep)
+ self.h5.createArray("/solution/snapshot%d" % istep,
+ "displacements", displacements)
+ self.h5.createArray("/solution/snapshot%d" % istep,
+ "velocities", velocities)
+ return
+
+
+ def _loadInterpolation(self):
+ """
+ Load interpolation and quadrature information.
+ """
+ self._info.log("Loading interpolation mapping...")
+
+ from numpy.linalg import inv, det
+ coordinates = self.h5.root.geometry.coordinates[:]
+ simplices = self.h5.root.topology.simplices[:]
+ map = self.discretization.calcMap(coordinates, simplices)
+ self.h5.createGroup("/", "interpolation")
+ self.h5.createGroup("/interpolation", "map")
+ self.h5.createArray("/interpolation/map", "jacobian", map['jacobian'])
+ self.h5.createArray("/interpolation/map", "to_ref", map['toRef'])
+ self.h5.createArray("/interpolation/map", "from_ref", map['fromRef'])
+ self.h5.createArray("/interpolation/map", "origin", map['origin'])
+ return
+
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+ app = Translator()
+ app.run()
+
+# End of file
Property changes on: cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/benchmark/trunk/benchmark/short/examples/compare_solns.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/examples/compare_solns.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/examples/compare_solns.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file benchmark/short/applications/bm_compare_solns.py
+
+## @brief Python script to compare benchmark solutions.
+
+scenario = "tet4_1000m"
+field = "displacements"
+
+# ----------------------------------------------------------------------
+if __name__ == "__main__":
+ from benchmark.CompareApp import CompareApp
+ app = CompareApp()
+ app.run()
+
+ projection = app.projection
+
+ projection.open()
+
+ # Time step t=0 --------------------------
+ tstep = 0
+ print "Time Step: %d" % tstep
+
+ filename = "output/pylith-0.8/%s.h5" % scenario
+ projection.project(filename, "pylith_0_8", tstep, field)
+
+ filename = "output/geofest/%s.h5" % scenario
+ projection.project(filename, "geofest", tstep, field)
+
+ filename = "output/femlab/tet10_2000m.h5"
+ projection.copy_projection(filename, "femlab", tstep, field)
+
+ filename = "output/analytic/%s.h5" % scenario
+ projection.copy_projection(filename, "analytic", tstep, field)
+
+ err = projection.difference(tstep, field, "pylith_0_8", "analytic")
+ print "PyLith-0.8, error wrt analytic: %12.4e" % err
+
+ err = projection.difference(tstep, field, "geofest", "analytic")
+ print "GeoFEST, error wrt analytic: %12.4e" % err
+
+ err = projection.difference(tstep, field, "femlab", "analytic")
+ print "Femlab, error wrt analytic: %12.4e" % err
+
+ err = projection.difference(tstep, field, "pylith_0_8", "geofest")
+ print "PyLith-GeoFEST, difference: %12.4e" % err
+
+ # Time step 10,100 --------------------------
+ #for tstep in [10,100]:
+ for tstep in []:
+ print "Time Step: %d" % tstep
+
+ filename = "output/pylith-0.8/%s.h5" % scenario
+ projection.project(filename, "pylith_0_8", tstep, field)
+
+ filename = "output/geofest/%s.h5" % scenario
+ projection.project(filename, "geofest", tstep, field)
+
+ err = projection.difference(tstep, field, "pylith_0_8", "geofest")
+ print "PyLith-GeoFEST, difference: %12.4e" % err
+
+ # ----------------------------------------
+ projection.close()
+
+# End of file
Property changes on: cs/benchmark/trunk/benchmark/short/examples/compare_solns.py
___________________________________________________________________
Name: svn:executable
+ *
Added: cs/benchmark/trunk/benchmark/short/examples/mesh_catalog.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/examples/mesh_catalog.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/short/examples/mesh_catalog.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,27 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+catalog = {
+ 'tet4_1000m': {'name': "tet4_1000m",
+ 'spatialDim': 3,
+ 'connSize': 4,
+ 'numNodes': 16071,
+ 'numElements': 85368},
+
+ 'tet4_0500m': {'name': "tet4_0500m",
+ 'spatialDim': 3,
+ 'connSize': 4,
+ 'numNodes': 119451,
+ 'numElements': 673396}
+ }
+
+# End of file
Added: cs/benchmark/trunk/benchmark/tablewaiter.py
===================================================================
--- cs/benchmark/trunk/benchmark/tablewaiter.py 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/benchmark/tablewaiter.py 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file tablewaiter.py
+
+## @brief Python functions for doing simple operations on
+## PyTables. These are high level functions.
+
+import tables
+
+def getNode(root, name):
+ """
+ Get node in table where name specifies path relative to root.
+ """
+ parent = root
+ lineage = name.split('/')
+ for child in lineage:
+ if not child == "":
+ parent = getattr(parent, child)
+ return parent
+
+
+def createGroup(h5file, root, name):
+ """
+ Create group where name is path relative to root. Any ancestors are
+ created as necessary (ala mkdir -p) and the desired group is returned.
+ """
+
+ parent = root
+ lineage = name.split('/')
+ for child in lineage:
+ if not child == "":
+ if not child in parent:
+ h5file.createGroup(parent, child)
+ parent = getattr(parent, child)
+ return parent
+
+
+# End of file
Added: cs/benchmark/trunk/configure.ac
===================================================================
--- cs/benchmark/trunk/configure.ac 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/configure.ac 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,51 @@
+# -*- autoconf -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+AC_PREREQ(2.59)
+AC_INIT([benchmark], [0.3.0], [baagaard at usgs.gov])
+AC_CONFIG_HEADER([portinfo])
+AC_CONFIG_AUX_DIR([./aux-config])
+AM_INIT_AUTOMAKE([foreign])
+
+# ----------------------------------------------------------------------
+# PYTHON
+AC_ARG_VAR(PYTHON, [Python interpreter])
+AC_ARG_VAR(PYTHONPATH, [Python module search path])
+
+# ----------------------------------------------------------------------
+
+AC_PROG_INSTALL
+AM_PATH_PYTHON([2.3])
+
+# ----------------------------------------------------------------------
+# PYTHON
+AC_CACHE_CHECK([for $am_display_PYTHON include directory],
+ [PYTHON_INCDIR],
+ [PYTHON_INCDIR=`$PYTHON -c "from distutils import sysconfig; print sysconfig.get_python_inc()" 2>/dev/null ||
+ echo "$PYTHON_PREFIX/include/python$PYTHON_VERSION"`])
+AC_SUBST([PYTHON_INCDIR], [$PYTHON_INCDIR])
+
+# PYTHIA
+# TODO: Add test for Pythia components. We only use Python Pyre stuff,
+# so perhaps just try to import journal?
+
+# ----------------------------------------------------------------------
+AC_CONFIG_FILES([Makefile
+ benchmark/Makefile])
+
+AC_OUTPUT
+
+
+# version
+# $Id$
+
+# End of file
Added: cs/benchmark/trunk/tidy.sh
===================================================================
--- cs/benchmark/trunk/tidy.sh 2006-08-29 23:44:53 UTC (rev 4450)
+++ cs/benchmark/trunk/tidy.sh 2006-08-29 23:54:11 UTC (rev 4451)
@@ -0,0 +1,23 @@
+#!/bin/bash
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+# Remove autoconf stuff
+autoconf_files="aclocal.m4 configure portinfo.in"
+autoconf_dirs="autom4te.cache aux-config"
+rm -r $autoconf_dirs
+rm $autoconf_files
+rm `find . -name Makefile.in`
+
+# Remove emacs backup stuff
+rm `find . -name "*~"`
+
+# End of file
Property changes on: cs/benchmark/trunk/tidy.sh
___________________________________________________________________
Name: svn:executable
+ *
More information about the cig-commits
mailing list