[cig-commits] r4458 - in cs/benchmark/trunk/benchmark: . discretize
short/applications short/examples
baagaard at geodynamics.org
baagaard at geodynamics.org
Thu Aug 31 14:16:09 PDT 2006
Author: baagaard
Date: 2006-08-31 14:16:09 -0700 (Thu, 31 Aug 2006)
New Revision: 4458
Modified:
cs/benchmark/trunk/benchmark/ProjectionQuadPts.py
cs/benchmark/trunk/benchmark/discretize/Tet4.py
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/compare_solns.py
Log:
Cleaned up used of HDF5 files to reduce storage. Switched to 32-bit floats and ints. Compute mappings on the fly; only store jacobians.
Modified: cs/benchmark/trunk/benchmark/ProjectionQuadPts.py
===================================================================
--- cs/benchmark/trunk/benchmark/ProjectionQuadPts.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/ProjectionQuadPts.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -17,6 +17,7 @@
from pyre.components.Component import Component
+import numpy
import tables
import tablewaiter
@@ -92,6 +93,7 @@
"""
self._copySolnInfo(solution)
+ self._setupIntegration()
self._calcProjectPts()
return
@@ -115,27 +117,29 @@
return
- def project(self, filename, label, timestep, field, overwrite=True):
+ def project(self, solnfile, solnmap, label, timestep, field,
+ overwrite=True):
"""
Project solution in to quadrature points.
- @param filename Name of solution file
+ @param solnfile Solution file
+ @param solnmap Mapping to/from reference element for
+ solution discretization
@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
+ from numpy import 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[:]
+ (quadpts, quadwts) = \
+ self.discretization.getQuadPts(self.quadratureOrder)
nquadpts = quadpts.shape[0]
ndims = quadpts.shape[1]
quadcoords = [map(float, tuple(quadpts[q])) \
@@ -145,19 +149,20 @@
basisfns = self.discretization.getBasisFns()
# Tabulate basis fns at quadrature points
- basisfns_quad = array(basisfns.tabulate(quadcoords)).transpose()
+ basisfns_quad = numpy.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)
+ fieldProj = numpy.zeros((nelems, nquadpts, ndims),
+ dtype=numpy.Float32)
- fieldLocal = zeros((connSize, ndims), dtype=float)
for ielem in xrange(nelems):
nodeIndices = simplices[ielem]
- fieldLocal[:] = solnfield[nodeIndices,:]
- fieldProj[ielem] = dot(basisfns_quad, fieldLocal)
+ fieldLocal = solnfield[nodeIndices,:]
+ tmp = dot(basisfns_quad, fieldLocal)
+ fieldProj[ielem] = tmp.astype(numpy.Float32)
group = tablewaiter.createGroup(self.workspace,
self.workspace.root,
@@ -170,15 +175,13 @@
return
- def copy_projection(self, filename, label, timestep, field,
+ def copy_projection(self, filein, 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[:]
@@ -192,8 +195,6 @@
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" % \
@@ -212,7 +213,7 @@
'error'.
"""
- from numpy import zeros, dot, sqrt, reshape
+ from numpy import dot, sqrt, reshape
dataA = tablewaiter.getNode(self.workspace.root,
"projection/data/%s/snapshot%d/%s" % \
@@ -220,16 +221,24 @@
dataB = tablewaiter.getNode(self.workspace.root,
"projection/data/%s/snapshot%d/%s" % \
(labelB, timestep, field))[:]
- quadwts = self.workspace.root.interpolation.quadrature.weights[:]
+ (quadpts, quadwts) = \
+ self.discretization.getQuadPts(self.quadratureOrder)
nquadpts = quadwts.shape[0]
quadwts = reshape(quadwts, (1,nquadpts))
- jacobian = self.workspace.root.interpolation.map.jacobian
+ jacobian = self.workspace.root.integration.jacobian
+ ndims = self.workspace.root._v_attrs.spatial_dim
+
nelems = self.workspace.root._v_attrs.num_simplices
- ndims = self.workspace.root._v_attrs.spatial_dim
- localerr = zeros( (nelems,), dtype=float)
+ localerr = numpy.zeros( (nelems,), dtype=numpy.Float32)
+ globalerr = 0.0
+ globalvol = 0.0
for ielem in xrange(nelems):
diff = reshape(dataA[ielem] - dataB[ielem], (nquadpts, ndims))
- localerr[ielem] = jacobian[ielem] * sum(sum(dot(quadwts, diff**2)))
+ localvol = jacobian[ielem] * sum(quadwts)
+ err = jacobian[ielem] * sum(sum(dot(quadwts, diff**2)))
+ localerr[ielem] = sqrt(err / localvol)
+ globalerr += err
+ globalvol += localvol
group = tablewaiter.createGroup(self.workspace,
self.workspace.root,
"difference/%s_%s/snapshot%d" % \
@@ -239,7 +248,7 @@
dataset = getattr(group, field)
dataset._f_remove()
self.workspace.createArray(group, field, localerr)
- return sqrt(sum(localerr))
+ return sqrt(globalerr/globalvol)
def remove(self, name, recursive=True):
@@ -273,13 +282,25 @@
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 _setupIntegration(self):
+ """
+ Setup quadrature information.
+ """
+ self._info.log("Seeting up integration information...")
+
+ workspace = self.workspace
+ vertices = workspace.root.geometry.vertices[:]
+ simplices = workspace.root.topology.simplices[:]
+ map = self.discretization.calcMap(vertices, simplices)
+ workspace.createGroup("/", "integration")
+ workspace.createArray("/integration", "jacobian", map['jacobian'])
+ return
+
+
def _calcProjectPts(self):
"""
Calculate coordinates of points in projection.
@@ -291,27 +312,26 @@
# 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[:]
+ from numpy import dot
+ vertices = workspace.root.geometry.vertices[:]
+ simplices = workspace.root.topology.simplices[:]
+ mapping = self.discretization.calcMap(vertices, simplices)
+ mapOrigin = mapping['origin'][:]
+ mapFromRef = mapping['fromRef'][:]
nelems = workspace.root._v_attrs.num_simplices
nquadpts = len(quadpts)
ndims = len(quadpts[0])
- points = zeros((nelems, nquadpts, ndims), dtype=float)
+ points = numpy.zeros((nelems, nquadpts, ndims), dtype=numpy.Float32)
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])
+ pts = origin + dot(j, quadpts[iquad])
+ points[ielem, iquad] = pts.astype(numpy.Float32)
group = tablewaiter.createGroup(workspace, workspace.root,
"projection")
Modified: cs/benchmark/trunk/benchmark/discretize/Tet4.py
===================================================================
--- cs/benchmark/trunk/benchmark/discretize/Tet4.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/discretize/Tet4.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -34,15 +34,15 @@
return
- def calcMap(self, coordinates, simplices):
+ def calcMap(self, vertices, 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]
+ nvertices = vertices.shape[0]
+ ndims = vertices.shape[1]
nelems = simplices.shape[0]
jacobian = zeros((nelems,), dtype=float)
toRef = zeros((nelems, ndims, ndims), dtype=float)
@@ -52,10 +52,10 @@
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)
+ x0,y0,z0 = vertices[a] # a = (x0,y0,z0)
+ x1,y1,z1 = vertices[b] # b = (x1,y1,z1)
+ x2,y2,z2 = vertices[c] # c = (x2,y2,z2)
+ x3,y3,z3 = vertices[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)
Modified: cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/short/applications/bm_create_projection.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -95,6 +95,9 @@
# ----------------------------------------------------------------------
if __name__ == "__main__":
+ import journal
+ journal.info("projectionquadpts").activate()
+
app = Creator()
app.run()
Modified: cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/short/applications/femlabtoh5.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -16,7 +16,7 @@
## at projection quadrature points to HDF5 file.
from pyre.applications.Script import Script
-from numpy import zeros
+import numpy
from benchmark.tablewaiter import createGroup
# Translator class
@@ -67,8 +67,8 @@
npts = self._countLines(self.asciiFile)
ndims = 3
- points = zeros( (npts, ndims), dtype=float)
- displacements = zeros( (npts, ndims), dtype=float)
+ points = numpy.zeros( (npts, ndims), dtype=numpy.Float32)
+ displacements = numpy.zeros( (npts, ndims), dtype=numpy.Float32)
filein = open(self.asciiFile, 'r')
ipt = 0
Modified: cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/short/applications/geofesttoh5.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -16,7 +16,7 @@
## format to HDF5.
from pyre.applications.Script import Script
-from numpy import zeros
+import numpy
# Translator class
class Translator(Script):
@@ -73,7 +73,6 @@
self._loadGeometry()
self._loadTopology()
self._loadSolution()
- self._loadInterpolation()
self.h5.close()
return
@@ -115,6 +114,7 @@
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
+ self.h5.root._v_attrs.simplex_type = mesh['simplexType']
return
@@ -124,7 +124,8 @@
"""
self._info.log("Loading geometry...")
- coordinates = zeros((self.numNodes, self.spatialDim), dtype=float)
+ vertices = numpy.zeros((self.numNodes, self.spatialDim),
+ dtype=numpy.Float32)
filename = "%s.coord" % self.scenarioRoot
filein = open(filename, 'r')
@@ -135,7 +136,7 @@
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)
+ vertices[label-1] = (scale*x, scale*y, scale*z)
else:
import pyre.units
parser = pyre.units.parser()
@@ -143,7 +144,7 @@
filein.close()
self.h5.createGroup("/", "geometry")
- self.h5.createArray("/geometry", "coordinates", coordinates)
+ self.h5.createArray("/geometry", "vertices", vertices)
return
@@ -153,7 +154,8 @@
"""
self._info.log("Loading topology...")
- simplices = zeros((self.numElements, self.connSize), dtype=int)
+ simplices = numpy.zeros((self.numElements, self.connSize),
+ dtype=numpy.Int32)
filename = "%s.connect" % self.scenarioRoot
filein = open(filename, 'r')
@@ -186,8 +188,10 @@
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)
+ displacements = numpy.zeros((self.numNodes, self.spatialDim),
+ dtype=numpy.Float32)
+ velocities = numpy.zeros((self.numNodes, self.spatialDim),
+ dtype=numpy.Float32)
filein = open(filename, 'r')
numHeaderLines = 4
for iskip in xrange(numHeaderLines):
@@ -218,25 +222,6 @@
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()
Modified: cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/short/applications/okadaelastictoh5.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -16,7 +16,7 @@
## using Okada at projection quadrature points to HDF5 file.
from pyre.applications.Script import Script
-from numpy import zeros
+import numpy
from benchmark.tablewaiter import createGroup
# Translator class
@@ -67,8 +67,8 @@
npts = self._countLines(self.asciiFile)
ndims = 3
- points = zeros( (npts, ndims), dtype=float)
- displacements = zeros( (npts, ndims), dtype=float)
+ points = numpy.zeros( (npts, ndims), dtype=numpy.Float32)
+ displacements = numpy.zeros( (npts, ndims), dtype=numpy.Float32)
filein = open(self.asciiFile, 'r')
ipt = 0
Modified: cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/short/applications/pylithtoh5.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -16,7 +16,7 @@
## format to HDF5.
from pyre.applications.Script import Script
-from numpy import zeros
+import numpy
# Translator class
class Translator(Script):
@@ -77,7 +77,6 @@
self._loadGeometry()
self._loadTopology()
self._loadSolution()
- self._loadInterpolation()
self.h5.close()
return
@@ -120,6 +119,7 @@
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
+ self.h5.root._v_attrs.simplex_type = mesh['simplexType']
return
@@ -129,7 +129,8 @@
"""
self._info.log("Loading geometry...")
- coordinates = zeros((self.numNodes, self.spatialDim), dtype=float)
+ vertices = numpy.zeros((self.numNodes, self.spatialDim),
+ dtype=numpy.Float32)
for iproc in xrange(self.numProcs):
filename = "%s_%d.%d.coord" % \
@@ -142,7 +143,7 @@
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)
+ vertices[label-1] = (scale*x, scale*y, scale*z)
else:
import pyre.units
parser = pyre.units.parser()
@@ -150,7 +151,7 @@
filein.close()
self.h5.createGroup("/", "geometry")
- self.h5.createArray("/geometry", "coordinates", coordinates)
+ self.h5.createArray("/geometry", "vertices", vertices)
return
@@ -160,7 +161,8 @@
"""
self._info.log("Loading topology...")
- simplices = zeros((self.numElements, self.connSize), dtype=int)
+ simplices = numpy.zeros((self.numElements, self.connSize),
+ dtype=numpy.Int32)
for iproc in xrange(self.numProcs):
filename = "%s_%d.%d.connect" % \
@@ -193,8 +195,10 @@
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)
+ displacements = numpy.zeros((self.numNodes, self.spatialDim),
+ dtype=numpy.Float32)
+ velocities = numpy.zeros((self.numNodes, self.spatialDim),
+ dtype=numpy.Float32)
for istep in self.timeSteps:
displacements[:] = 0.0
velocities[:] = 0.0
@@ -222,25 +226,6 @@
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()
Modified: cs/benchmark/trunk/benchmark/short/examples/compare_solns.py
===================================================================
--- cs/benchmark/trunk/benchmark/short/examples/compare_solns.py 2006-08-31 18:21:45 UTC (rev 4457)
+++ cs/benchmark/trunk/benchmark/short/examples/compare_solns.py 2006-08-31 21:16:09 UTC (rev 4458)
@@ -14,35 +14,100 @@
## @brief Python script to compare benchmark solutions.
-scenario = "tet4_1000m"
-field = "displacements"
+import tables
# ----------------------------------------------------------------------
-def local_main(app):
+def getMap(solnfile):
"""
- Local replacement for application 'main'.
+ Get element interpolation mapping information.
"""
+ vertices = solnfile.root.geometry.vertices[:]
+ simplices = solnfile.root.topology.simplices[:]
+ simplexType = solnfile.root._v_attrs.simplex_type
+
+ element = None
+ if simplexType.lower() == "tet4":
+ from benchmark.discretize.Tet4 import Tet4
+ element = Tet4()
+ else:
+ raise ValueError("Unknown type of finite element '%s'." % \
+ simplexType)
+ return element.calcMap(vertices, simplices)
+
+
+# ----------------------------------------------------------------------
+def main_load_disp(app):
+ """
+ Local replacement for CompareApp::main().
+
+ Load displacement fields.
+ """
+ scenario = "tet4_1000m"
+ field = "displacements"
+
projection = app.projection
projection.open()
- # Time step t=0 --------------------------
- tstep = 0
- print "Time Step: %d" % tstep
-
+ # PyLith ---------------------------------
+ app._info.log("Projecting PyLith solution...")
filename = "output/pylith-0.8/%s.h5" % scenario
- projection.project(filename, "pylith_0_8", tstep, field)
+ solnfile = tables.openFile(filename, 'r')
+ solnmap = getMap(solnfile)
+ for tstep in [0,10,50,100]:
+ projection.project(solnfile, solnmap, "pylith_0_8", tstep, field)
+ solnfile.close()
+ # GeoFEST --------------------------------
+ app._info.log("Projecting GeoFEST solution...")
filename = "output/geofest/%s.h5" % scenario
- projection.project(filename, "geofest", tstep, field)
+ solnfile = tables.openFile(filename, 'r')
+ solnmap = getMap(solnfile)
+ for tstep in [0,10,50,100]:
+ projection.project(solnfile, solnmap, "geofest", tstep, field)
+ solnfile.close()
+ # Femlab ---------------------------------
+ app._info.log("Copying Femlab solution...")
filename = "output/femlab/tet10_2000m.h5"
- projection.copy_projection(filename, "femlab", tstep, field)
+ solnfile = tables.openFile(filename, 'r')
+ for tstep in [0]:
+ projection.copy_projection(solnfile, "femlab", tstep, field)
+ solnfile.close()
+ # Analytic -------------------------------
+ app._info.log("Copying analytic solution...")
filename = "output/analytic/%s.h5" % scenario
- projection.copy_projection(filename, "analytic", tstep, field)
+ solnfile = tables.openFile(filename, 'r')
+ for tstep in [0]:
+ projection.copy_projection(solnfile, "analytic", tstep, field)
+ solnfile.close()
+ # ----------------------------------------
+ projection.close()
+
+ return
+
+# ----------------------------------------------------------------------
+def main_cmp_disp(app):
+ """
+ Local replacement for CompareApp::main().
+
+ Compare displacement fields.
+ """
+
+ scenario = "tet4_1000m"
+ field = "displacements"
+
+ projection = app.projection
+
+ projection.open()
+
+ # Code versus analytic -------------------
+ tstep = 0
+ print "Time Step: %d" % tstep
+
err = projection.difference(tstep, field, "pylith_0_8", "analytic")
print "PyLith-0.8, error wrt analytic: %12.4e" % err
@@ -52,20 +117,10 @@
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 []:
+ # PyLith versus GeoFEST ------------------
+ for tstep in [0,10,50,100]:
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
@@ -74,11 +129,16 @@
return
+
# ----------------------------------------------------------------------
if __name__ == "__main__":
+ import journal
+ journal.info("compareapp").activate()
+
from benchmark.CompareApp import CompareApp
app = CompareApp()
- app.mainfn = local_main
+ app.mainfn = main_cmp_disp
app.run()
+
# End of file
More information about the cig-commits
mailing list