[cig-commits] r19828 - in short/3D/PyLith/branches/pylith-scecdynrup: . examples/3d/hex8 libsrc/pylith/feassemble libsrc/pylith/meshio modulesrc/meshio pylith/meshio pylith/problems unittests/libtests/meshio
brad at geodynamics.org
brad at geodynamics.org
Tue Mar 20 09:41:07 PDT 2012
Author: brad
Date: 2012-03-20 09:41:07 -0700 (Tue, 20 Mar 2012)
New Revision: 19828
Added:
short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/GreensFns.py
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/TODO
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/feassemble/CellGeometry.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.hh
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i
short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestDataWriterVTKPoints.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc
Log:
Merge from trunk.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-03-20 16:41:07 UTC (rev 19828)
@@ -1,4 +1,28 @@
======================================================================
+CURRENT BUGS
+======================================================================
+
+* Deallocation with field split
+
+ Not all PETSc memory is being freed when using field split.
+
+ Test cases (examine True_0):
+ cd examples/3d/tet4
+ pylith step02.cfg ../../../share/debug_malloc.cfg
+ pylith step04.cfg ../../../share/debug_malloc.cfg
+
+* Point interpolation
+
+ Fails in parallel.
+
+ Test case:
+ cd examples/3d/hex8
+ pylith step19.cfg --nodes=1 # OK
+ pylith step19.cfg --nodes=2 # ERROR
+ [0]PETSC ERROR: Petsc has generated inconsistent data!
+ [0]PETSC ERROR: Invalid number of points located 9 should be 8!
+
+======================================================================
CURRENT ISSUES/PRIORITIES
======================================================================
@@ -6,12 +30,11 @@
- Order of tensor components for Xdmf files
- Drucker Prager fit to yield surface
+ - Drucker Prage allow tensile yield
- Slip-weakening healing parameter
- Green's functions
+ Interpolation
- * Approximate for quad4 and hex8 cells
- (assume Jacobian is uniform throughout cell =? rhomboid)
- (should be much better than nearest)
+ * Accuracy of interpolation?
- Examples/Tutorials
+ Table with concepts/features
@@ -31,11 +54,8 @@
* 2-D materials
- + DruckerPragerPlaneStrain (Drucker-Prager plane strain ) [Charles]
+ PowerLawPlaneStrain (power law plane strain ) [Charles]
-* GenMaxwellQpQs [BRAD]
-
* Cleanup
Add elasticPrestep() to Formulation (called from Problem)
@@ -49,6 +69,8 @@
spatial and temporal variation in shear/normal stress for nucleation
+* GenMaxwellQpQs [BRAD]
+
* GPU utilization
Finite-element integrations
@@ -75,6 +97,9 @@
* Coupling
+ Use slip rate threshold for quasi-static -> dynamic transition.
+ Use velocity threshold for dynamic -> quasi-static transition.
+
* Field split.
Add flag to material [default is false] for creating null vector
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt 2012-03-20 16:41:07 UTC (rev 19828)
@@ -7,7 +7,7 @@
# coordsys component of the OutputSolnPoints object. The points will
# be transformed into the coordinate system of the mesh before
# interpolation.
-10.0 10.0 -500.0
-10.0 10.0 -1500.0
-10.0 10.0 -2500.0
-10.0 10.0 -3500.0
+0.0001 0.0001 -500.0
+0.0001 0.0001 -1500.0
+0.0001 0.0001 -2500.0
+0.0001 0.0001 -3500.0
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg 2012-03-20 16:41:07 UTC (rev 19828)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# This is not a self-contained simulation configuration file. This
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/feassemble/CellGeometry.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/feassemble/CellGeometry.cc 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/feassemble/CellGeometry.cc 2012-03-20 16:41:07 UTC (rev 19828)
@@ -249,7 +249,15 @@
PylithScalar p2 = upDir[0]*r1 - upDir[1]*r0;
// Make unit vector
mag = sqrt(p0*p0 + p1*p1 + p2*p2);
- assert(mag > 0.0);
+ if (0.0 == mag) {
+ std::ostringstream msg;
+ msg << "Error computing orientation of cell face.\nUp direction ("
+ << upDir[0] << ", " << upDir[1] << ", " << upDir[2] << ") "
+ << " cannot be parallel to the face normal ("
+ << r0 << ", " << r1 << ", " << r2 << ").\n"
+ << "Adjust the up_dir parameter.";
+ throw std::runtime_error(msg.str());
+ } // if
p0 /= mag;
p1 /= mag;
p2 /= mag;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-20 16:41:07 UTC (rev 19828)
@@ -84,7 +84,8 @@
pylith::meshio::OutputSolnPoints::setupInterpolator(topology::Mesh* mesh,
const PylithScalar* points,
const int numPoints,
- const int spaceDim)
+ const int spaceDim,
+ const spatialdata::units::Nondimensional& normalizer)
{ // setupInterpolator
assert(mesh);
assert(points);
@@ -123,6 +124,10 @@
cells, numCells, numCorners, meshDim,
interpolate);
_pointsMesh->coordsys(_mesh->coordsys());
+ _pointsMesh->nondimensionalize(normalizer);
+#if 1 // DEBUGGING
+ _pointsMesh->view("POINTS MESH");
+#endif
// Setup interpolator object
DM dm;
@@ -139,7 +144,13 @@
err = DMMeshInterpolationSetDim(dm, spaceDim,
_interpolator);CHECK_PETSC_ERROR(err);
- err = DMMeshInterpolationAddPoints(dm, numPoints, (PetscReal*)points,
+
+ assert(!_pointsMesh->sieveMesh().isNull());
+ assert(_pointsMesh->sieveMesh()->hasRealSection("coordinates"));
+ const ALE::Obj<topology::Mesh::RealSection>& coordinatesSection = _pointsMesh->sieveMesh()->getRealSection("coordinates");
+ assert(!coordinatesSection.isNull());
+ const PylithScalar* coordinates = coordinatesSection->restrictSpace();
+ err = DMMeshInterpolationAddPoints(dm, numPoints, (PetscReal*)coordinates,
_interpolator);CHECK_PETSC_ERROR(err);
err = DMMeshInterpolationSetUp(dm, _interpolator);CHECK_PETSC_ERROR(err);
err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.hh 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.hh 2012-03-20 16:41:07 UTC (rev 19828)
@@ -64,14 +64,16 @@
/** Setup interpolator.
*
* @param mesh Domain mesh.
- * @param points Array of coordinates for points [numPoints*spaceDim].
+ * @param points Array of dimensioned coordinates for points [numPoints*spaceDim].
* @param numPoints Number of points.
* @param spaceDim Spatial dimension for coordinates.
+ * @param normalizer Nondimensionalizer.
*/
void setupInterpolator(topology::Mesh* mesh,
const PylithScalar* points,
const int numPoints,
- const int spaceDim);
+ const int spaceDim,
+ const spatialdata::units::Nondimensional& normalizer);
/** Prepare for output.
*
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i 2012-03-20 16:41:07 UTC (rev 19828)
@@ -52,7 +52,7 @@
/** Setup interpolator.
*
* @param mesh Domain mesh.
- * @param points Array of coordinates for points [numPoints*spaceDim].
+ * @param points Array of dimensioned coordinates for points [numPoints*spaceDim].
* @param numPoints Number of points.
* @param spaceDim Spatial dimension for coordinates.
*/
@@ -64,7 +64,8 @@
void setupInterpolator(pylith::topology::Mesh* mesh,
const PylithScalar* points,
const int numPoints,
- const int spaceDim);
+ const int spaceDim,
+ const spatialdata::units::Nondimensional& normalizer);
%clear(const PylithScalar* points, const int numPoints, const int spaceDim);
/** Prepare for output.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py 2012-03-20 16:41:07 UTC (rev 19828)
@@ -113,10 +113,8 @@
# Convert to mesh coordinate system
from spatialdata.geocoords.Converter import convert
convert(points, mesh.coordsys(), self.coordsys)
- print "LENGTH SCALE",normalizer.lengthScale()
- points /= normalizer.lengthScale().value
- ModuleOutputSolnPoints.setupInterpolator(self, mesh, points)
+ ModuleOutputSolnPoints.setupInterpolator(self, mesh, points, normalizer)
self.mesh = ModuleOutputSolnPoints.pointsMesh(self)
self._eventLogger.eventEnd(logEvent)
Copied: short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/GreensFns.py (from rev 19827, short/3D/PyLith/trunk/pylith/problems/GreensFns.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/GreensFns.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/GreensFns.py 2012-03-20 16:41:07 UTC (rev 19828)
@@ -0,0 +1,209 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/problems/GreensFns.py
+##
+## @brief Python abstract base class for time dependent crustal
+## dynamics problems.
+##
+## Factory: problem.
+
+from Problem import Problem
+
+# GreensFns class
+class GreensFns(Problem):
+ """
+ Python abstract base class for time dependent crustal dynamics problems.
+
+ Factory: problem.
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(Problem.Inventory):
+ """
+ Python object for managing GreensFns facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing GreensFns facilities and properties.
+ ##
+ ## \b Properties
+ ## None
+ ##
+ ## \b Facilities
+ ## @li \b formulation Formulation for solving PDE.
+ ## @li \b checkpoint Checkpoint manager.
+
+ import pyre.inventory
+
+ from Implicit import Implicit
+ formulation = pyre.inventory.facility("formulation",
+ family="pde_formulation",
+ factory=Implicit)
+ formulation.meta['tip'] = "Formulation for solving PDE."
+
+ from pylith.utils.CheckpointTimer import CheckpointTimer
+ checkpointTimer = pyre.inventory.facility("checkpoint",
+ family="checkpointer",
+ factory=CheckpointTimer)
+ checkpointTimer.meta['tip'] = "Checkpoint manager."
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="greensfns"):
+ """
+ Constructor.
+ """
+ Problem.__init__(self, name)
+ self._loggingPrefix = "PrGF "
+ return
+
+
+ def preinitialize(self, mesh):
+ """
+ Setup integrators for each element family (material/quadrature,
+ bc/quadrature, etc.).
+ """
+ self._setupLogging()
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Pre-initializing problem.")
+ self.mesh = mesh
+ self.formulation.preinitialize(mesh, self.materials, self.bc,
+ self.interfaces, self.gravityField)
+ return
+
+
+ def verifyConfiguration(self):
+ """
+ Verify compatibility of configuration.
+ """
+ Problem.verifyConfiguration(self)
+ self.formulation.verifyConfiguration()
+ return
+
+
+ def initialize(self):
+ """
+ Setup integrators for each element family (material/quadrature,
+ bc/quadrature, etc.).
+ """
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Initializing problem.")
+ self.checkpointTimer.initialize(self.normalizer)
+ self.formulation.initialize(self.dimension, self.normalizer)
+ return
+
+
+ def run(self, app):
+ """
+ Compute Green's functions associated with fault slip.
+ """
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Computing Green's functions.")
+ self.checkpointTimer.toplevel = app # Set handle for saving state
+
+ nimpulses = 0
+ ipulse = 0;
+ while ipulse < nimpulses:
+ self._eventLogger.stagePush("Prestep")
+ if 0 == comm.rank:
+ self._info.log("Main loop, impulse %d of %d." % (ipulse+1, nimpulses))
+
+ # Checkpoint if necessary
+ self.checkpointTimer.update(t)
+
+ if 0 == comm.rank:
+ self._info.log("Preparing impulse %d of %d." % \
+ (ipulse+1, nimpulses))
+ self.formulation.prestep(t, dt)
+ self._eventLogger.stagePop()
+
+ if 0 == comm.rank:
+ self._info.log("Computing response to impulse %d of %d." %
+ (ipulse+1, nimpulses))
+ self._eventLogger.stagePush("Step")
+ self.formulation.step(t, dt)
+ self._eventLogger.stagePop()
+
+ if 0 == comm.rank:
+ self._info.log("Finishing impulse %d of %d." % \
+ (ipulse+1, nimpulses))
+ self._eventLogger.stagePush("Poststep")
+ self.formulation.poststep(t, dt)
+ self._eventLogger.stagePop()
+
+ # Update time
+ t += dt
+ return
+
+
+ def finalize(self):
+ """
+ Cleanup after running problem.
+ """
+ self.formulation.finalize()
+ return
+
+
+ def checkpoint(self):
+ """
+ Save problem state for restart.
+ """
+ Problem.checkpoint()
+
+ # Save state of this object
+ raise NotImplementedError, "GreensFns::checkpoint() not implemented."
+
+ # Save state of children
+ self.formulation.checkpoint()
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based using inventory.
+ """
+ Problem._configure(self)
+ self.formulation = self.inventory.formulation
+ self.checkpointTimer = self.inventory.checkpointTimer
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def problem():
+ """
+ Factory associated with GreensFns.
+ """
+ return GreensFns()
+
+
+# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestDataWriterVTKPoints.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestDataWriterVTKPoints.cc 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestDataWriterVTKPoints.cc 2012-03-20 16:41:07 UTC (rev 19828)
@@ -74,12 +74,13 @@
OutputSolnPoints output;
DataWriterVTK<topology::Mesh, MeshField> writer;
+ spatialdata::units::Nondimensional normalizer;
writer.filename(_data->timestepFilename);
writer.timeFormat(_data->timeFormat);
output.writer(&writer);
output.setupInterpolator(_mesh, _data->points,
- _data->numPoints, _data->spaceDim);
+ _data->numPoints, _data->spaceDim, normalizer);
const PylithScalar t = _data->time;
const int numTimeSteps = 1;
@@ -109,6 +110,7 @@
OutputSolnPoints output;
DataWriterVTK<topology::Mesh, MeshField> writer;
+ spatialdata::units::Nondimensional normalizer;
topology::Fields<MeshField> vertexFields(*_mesh);
_createVertexFields(&vertexFields);
@@ -117,7 +119,7 @@
writer.timeFormat(_data->timeFormat);
output.writer(&writer);
output.setupInterpolator(_mesh, _data->points,
- _data->numPoints, _data->spaceDim);
+ _data->numPoints, _data->spaceDim, normalizer);
const int nfields = _data->numVertexFields;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc 2012-03-20 16:31:13 UTC (rev 19827)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc 2012-03-20 16:41:07 UTC (rev 19828)
@@ -113,6 +113,8 @@
topology::Mesh mesh;
spatialdata::geocoords::CSCart cs;
+ spatialdata::units::Nondimensional normalizer;
+
cs.setSpaceDim(spaceDim);
cs.initialize();
mesh.coordsys(&cs);
@@ -122,7 +124,7 @@
OutputSolnPoints output;
CPPUNIT_ASSERT(data.points);
- output.setupInterpolator(&mesh, data.points, numPoints, spaceDim);
+ output.setupInterpolator(&mesh, data.points, numPoints, spaceDim, normalizer);
const topology::Mesh& pointsMesh = output.pointsMesh();
const ALE::Obj<SieveMesh>& sievePointsMesh = pointsMesh.sieveMesh();
More information about the CIG-COMMITS
mailing list