[cig-commits] r18949 - in short/3D/PyLith/trunk: . modulesrc/bc modulesrc/faults modulesrc/utils pylith pylith/problems unittests/pytests/utils

brad at geodynamics.org brad at geodynamics.org
Tue Sep 20 09:36:13 PDT 2011


Author: brad
Date: 2011-09-20 09:36:13 -0700 (Tue, 20 Sep 2011)
New Revision: 18949

Removed:
   short/3D/PyLith/trunk/pylith/problems/ImplicitCUDA.py
Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/modulesrc/bc/AbsorbingDampers.i
   short/3D/PyLith/trunk/modulesrc/bc/BoundaryCondition.i
   short/3D/PyLith/trunk/modulesrc/bc/DirichletBC.i
   short/3D/PyLith/trunk/modulesrc/bc/DirichletBoundary.i
   short/3D/PyLith/trunk/modulesrc/bc/Neumann.i
   short/3D/PyLith/trunk/modulesrc/bc/PointForce.i
   short/3D/PyLith/trunk/modulesrc/faults/Fault.i
   short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i
   short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i
   short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveLagrange.i
   short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i
   short/3D/PyLith/trunk/modulesrc/utils/petsc_general.i
   short/3D/PyLith/trunk/modulesrc/utils/pylith_general.i
   short/3D/PyLith/trunk/pylith/Makefile.am
   short/3D/PyLith/trunk/pylith/problems/Formulation.py
   short/3D/PyLith/trunk/pylith/problems/Implicit.py
   short/3D/PyLith/trunk/pylith/problems/Solver.py
   short/3D/PyLith/trunk/unittests/pytests/utils/TestPetscManager.py
   short/3D/PyLith/trunk/unittests/pytests/utils/TestPylith.py
Log:
Added use_cuda property to Formulation (use CUDA for finite-element integrations) and Solver (use CUDA in solve).

Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/TODO	2011-09-20 16:36:13 UTC (rev 18949)
@@ -2,11 +2,16 @@
 CURRENT ISSUES/PRIORITIES (1.7.0)
 ======================================================================
 
+* configure
+  
+  + Check compatibility of PyLith options with PETSc
+    - HDF5
+    - CUDA
+
 * GPU utilization
 
-  PETSc w/GPU and single precision
-    PyLith use PetscScalar instead of double
-    spatialdata - create single precision versions
+  + Implicit elasticity finite-element integration
+  + Explicit elasticity finite-element integration
 
 * Green's functions
 

Modified: short/3D/PyLith/trunk/modulesrc/bc/AbsorbingDampers.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/AbsorbingDampers.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/bc/AbsorbingDampers.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -52,7 +52,7 @@
        *   direction that is not collinear with surface normal.
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
       
       /** Integrate contributions to residual term (r) for operator.
        *

Modified: short/3D/PyLith/trunk/modulesrc/bc/BoundaryCondition.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/BoundaryCondition.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/bc/BoundaryCondition.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -67,7 +67,7 @@
        */
       virtual
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]) = 0;
+		      const PylithScalar upDir[3]) = 0;
 
     }; // class BoundaryCondition
 

Modified: short/3D/PyLith/trunk/modulesrc/bc/DirichletBC.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/DirichletBC.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/bc/DirichletBC.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -54,7 +54,7 @@
        * @param upDir Vertical direction (somtimes used in 3-D problems).
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
       
       /** Set number of degrees of freedom that are constrained at
        * points in field.

Modified: short/3D/PyLith/trunk/modulesrc/bc/DirichletBoundary.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/DirichletBoundary.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/bc/DirichletBoundary.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -46,7 +46,7 @@
        * @param upDir Vertical direction (somtimes used in 3-D problems).
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
 
       /** Get boundary mesh.
        *

Modified: short/3D/PyLith/trunk/modulesrc/bc/Neumann.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/Neumann.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/bc/Neumann.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -47,7 +47,7 @@
        *   direction that is not collinear with surface normal.
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
 
       /** Integrate contributions to residual term (r) for operator.
        *

Modified: short/3D/PyLith/trunk/modulesrc/bc/PointForce.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/PointForce.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/bc/PointForce.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -46,7 +46,7 @@
        * @param upDir Vertical direction (somtimes used in 3-D problems).
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
       
       /** Integrate contributions to residual term (r) for operator.
        *

Modified: short/3D/PyLith/trunk/modulesrc/faults/Fault.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/Fault.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/faults/Fault.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -97,7 +97,7 @@
        */
       virtual
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]) = 0;
+		      const PylithScalar upDir[3]) = 0;
       
       /** Get mesh associated with fault fields.
        *

Modified: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -62,7 +62,7 @@
        *   be up-dip direction; applies to fault surfaces in 2-D and 3-D).
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
       
       /** Integrate contributions to residual term (r) for operator that
        * do not require assembly across processors.

Modified: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveKin.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -67,7 +67,7 @@
        *   be up-dip direction; applies to fault surfaces in 2-D and 3-D).
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
       
       /** Integrate contributions to residual term (r) for operator that
        * do not require assembly across cells, vertices, or processors.

Modified: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveLagrange.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveLagrange.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveLagrange.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -51,7 +51,7 @@
        */
       virtual
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
       
       /** Split solution field for separate preconditioning.
        *

Modified: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -56,7 +56,7 @@
        *   be up-dip direction; applies to fault surfaces in 2-D and 3-D).
        */
       void initialize(const pylith::topology::Mesh& mesh,
-		      const double upDir[3]);
+		      const PylithScalar upDir[3]);
 
       /** Integrate contribution of cohesive cells to residual term.
        *

Modified: short/3D/PyLith/trunk/modulesrc/utils/petsc_general.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/utils/petsc_general.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/utils/petsc_general.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -41,13 +41,15 @@
 %} // inline
 
 // ----------------------------------------------------------------------
-// sizeofVoidPtr
+// PetscOptionsSetValue
 %inline %{
   int
-  sizeofVoidPtr(void)
-  { // sizeofVoidPtr
-    return sizeof(void*);
-  } // sizeofVoidPtr
+  optionsSetValue(const char* name,
+		  const char* value)
+  { // optionsSetValue
+    PetscErrorCode err = PetscOptionsSetValue(name, value); CHKERRQ(err);
+    return 0;
+  } // optionsSetValue
 %} // inline
 
 

Modified: short/3D/PyLith/trunk/modulesrc/utils/pylith_general.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/utils/pylith_general.i	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/modulesrc/utils/pylith_general.i	2011-09-20 16:36:13 UTC (rev 18949)
@@ -17,6 +17,16 @@
 //
 
 // ----------------------------------------------------------------------
+// sizeofVoidPtr
+%inline %{
+  int
+  sizeofVoidPtr(void)
+  { // sizeofVoidPtr
+    return sizeof(void*);
+  } // sizeofVoidPtr
+%} // inline
+
+// ----------------------------------------------------------------------
 // sizeofPylithScalar
 %inline %{
   int
@@ -27,5 +37,20 @@
 %} // inline
 
 
+// ----------------------------------------------------------------------
+// isCUDAEnabled
+%inline %{
+  bool
+  isCUDAEnabled(void)
+  { // isCUDAEnabled
+#if ENABLE_CUDA
+    return true;
+#else
+    return false;
+#endif
+  } // isCUDAEnabled
+%} // inline
+
+
 // End of file
 

Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/pylith/Makefile.am	2011-09-20 16:36:13 UTC (rev 18949)
@@ -195,8 +195,7 @@
 
 if ENABLE_CUDA
   nobase_pkgpyexec_PYTHON += \
-	feassemble/ElasticityImplicitCUDA.py \
-	problems/ImplicitCUDA.py
+	feassemble/ElasticityImplicitCUDA.py
 endif
 
 if ENABLE_TETGEN

Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py	2011-09-20 16:36:13 UTC (rev 18949)
@@ -29,6 +29,16 @@
 from pylith.utils.profiling import resourceUsageString
 from pyre.units.time import second
 
+# VALIDATORS ///////////////////////////////////////////////////////////
+
+# Validate use of CUDA.
+def validateUseCUDA(value):
+  from pylith.utils.utils import isCUDAEnabled
+  if value and not isCUDAEnabled:
+    raise ValueError("PyLith is not built with CUDA support.")
+  return value
+
+
 # ITEM FACTORIES ///////////////////////////////////////////////////////
 
 def outputFactory(name):
@@ -62,6 +72,7 @@
     ## Python object for managing Formulation facilities and properties.
     ##
     ## \b Properties
+    ## @li \b use_cuda Enable use of CUDA for finite-element integrations.
     ## @li \b matrix_type Type of PETSc sparse matrix.
     ## @li \b split_fields Split solution fields into displacements
     ## @li \b use_custom_constraint_pc Use custom preconditioner for
@@ -77,6 +88,10 @@
 
     import pyre.inventory
 
+    useCUDA = pyre.inventory.bool("use_cuda", default=False,
+                                  validator=validateUseCUDA)
+    useCUDA.meta['tip'] = "Enable use of CUDA for finite-element integrations."
+
     matrixType = pyre.inventory.str("matrix_type", default="unknown")
     matrixType.meta['tip'] = "Type of PETSc sparse matrix."
 
@@ -133,7 +148,6 @@
     self.constraints = None
     self.jacobian = None
     self.fields = None
-    self.solnName = None
     return
 
 
@@ -154,6 +168,7 @@
     self.constraints = []
     self.gravityField = gravityField
 
+    self.solver.preinitialize()
     self._setupMaterials(materials)
     self._setupBC(boundaryConditions)
     self._setupInterfaces(interfaceConditions)
@@ -290,6 +305,7 @@
     Set members based using inventory.
     """
     PetscComponent._configure(self)
+    self.useCUDA = self.inventory.useCUDA
     self.matrixType = self.inventory.matrixType
     self.timeStep = self.inventory.timeStep
     self.solver = self.inventory.solver
@@ -311,6 +327,7 @@
     ModuleFormulation.splitFields(self, self.inventory.useSplitFields)
     ModuleFormulation.useCustomConstraintPC(self,
                                             self.inventory.useCustomConstraintPC)
+
     return
 
 
@@ -334,12 +351,15 @@
               (self.matrixType, matrixMap[self.matrixType])
         self.matrixType = matrixMap[self.matrixType]
     self.blockMatrixOkay = True
+    if self.matrixType == "unknown" and self.solver.useCUDA:
+      self.matrixType = "mpiaijcusp"
     for constraint in self.constraints:
       numDimConstrained = constraint.numDimConstrained()
       if numDimConstrained > 0 and self.mesh.dimension() != numDimConstrained:
         self.blockMatrixOkay = False
     return
 
+
   def _setupMaterials(self, materials):
     """
     Setup materials as integrators.

Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py	2011-09-20 16:36:13 UTC (rev 18949)
@@ -105,8 +105,14 @@
     """
     Get integrator for elastic material.
     """
-    from pylith.feassemble.ElasticityImplicit import ElasticityImplicit
-    return ElasticityImplicit()
+    integrator = None
+    if self.useCUDA:
+      from pylith.feassemble.ElasticityImplicitCUDA import ElasticityImplicitCUDA
+      integrator = ElasticityImplicitCUDA()
+    else:
+      from pylith.feassemble.ElasticityImplicit import ElasticityImplicit
+      integrator = ElasticityImplicit()
+    return integrator
 
 
   def initialize(self, dimension, normalizer):

Deleted: short/3D/PyLith/trunk/pylith/problems/ImplicitCUDA.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/ImplicitCUDA.py	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/pylith/problems/ImplicitCUDA.py	2011-09-20 16:36:13 UTC (rev 18949)
@@ -1,298 +0,0 @@
-#!/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-2011 University of California, Davis
-#
-# See COPYING for license information.
-#
-# ----------------------------------------------------------------------
-#
-
-## @file pylith/problems/Implicit.py
-##
-## @brief Python Implicit object for solving equations using an
-## implicit formulation.
-##
-## Factory: pde_formulation
-
-from Formulation import Formulation
-from problems import Implicit as ModuleImplicit
-from pylith.utils.profiling import resourceUsageString
-
-# Implicit class
-class Implicit(Formulation, ModuleImplicit):
-  """
-  Python Implicit object for solving equations using an implicit
-  formulation.
-
-  The formulation has the general form,
-
-  [A(t+dt)] {u(t+dt)} = {b(t+dt)}. 
-
-  We know the solution at time t, so we write {u(t+dt)} as {u(t)} +
-  {du(t)}, where {du(t)} is the increment in the solution from time t
-  to time t+dt. Thus, we solve
-
-  [A(t+dt)] {du(t)} = {b(t+dt)} - [A(t+dt)]{u(t)}.
-
-  We solve this system by forming the Jacobian, A, and the residual
-
-  {r(t+dt)} = {b(t+dt)} - [A(t+dt)]{u(t)} - [A(t+dt)]{du(t)}
-
-  which we combine into
-
-  {r(t+dt)} = {b(t+dt)} - [A(t+dt)]{u(t)+du(t)}.
-
-  The method reformJacobian() computes [A(t+dt)] and the method
-  reformResidual computes {r(t+dt)}. Note that in forming the residual
-  we compute the action [A(t+dt)]{u(t)+du(t)} and do not perform a
-  matrix-vector multiplication.
-
-  [A(t+dt)] generally depends on {u(t+dt)} as well as the current
-  stresses and additional state variables.  
-
-  For linear elastic or viscoelastic problems with constant time step
-  size, A is a constant (after the elastic solution).  {b(t+dt)}
-  generally depends on the loads applied for time step t+dt (including
-  the contributions to the internal force vector from
-  displacement/velocity BC) as well as the internal force vector
-  computed from the current stresses.
-
-  Factory: pde_formulation.
-  """
-
-  # INVENTORY //////////////////////////////////////////////////////////
-
-  class Inventory(Formulation.Inventory):
-    """
-    Python object for managing Implicit facilities and properties.
-    """
-
-    ## @class Inventory
-    ## Python object for managing Implicit facilities and properties.
-    ##
-    ## \b Properties
-    ## @li None
-    ##
-    ## \b Facilities
-    ## @li None
-
-    import pyre.inventory
-
-
-  # PUBLIC METHODS /////////////////////////////////////////////////////
-
-  def __init__(self, name="implicit"):
-    """
-    Constructor.
-    """
-    Formulation.__init__(self, name)
-    ModuleImplicit.__init__(self)
-    self._loggingPrefix = "TSIm "
-    self._stepCount = None
-    return
-
-
-  def elasticityIntegrator(self):
-    """
-    Get integrator for elastic material.
-    """
-    from pylith.feassemble.ElasticityImplicitCUDA import ElasticityImplicitCUDA
-    return ElasticityImplicitCUDA()
-
-
-  def initialize(self, dimension, normalizer):
-    """
-    Initialize problem for implicit time integration.
-    """
-    logEvent = "%sinit" % self._loggingPrefix
-    self._eventLogger.eventBegin(logEvent)
-
-    self._initialize(dimension, normalizer)
-
-    from pylith.utils.petsc import MemoryLogger
-    memoryLogger = MemoryLogger.singleton()
-    memoryLogger.setDebug(0)
-    memoryLogger.stagePush("Problem")
-
-    # Allocate other fields, reusing layout from dispIncr
-    self._info.log("Creating other fields.")
-    self.fields.add("velocity(t)", "velocity")
-    self.fields.copyLayout("dispIncr(t->t+dt)")
-
-    # Setup fields and set to zero
-    dispT = self.fields.get("disp(t)")
-    dispT.zero()
-    residual = self.fields.get("residual")
-    residual.zero()
-    residual.createScatterMesh(residual.mesh())
-
-    lengthScale = normalizer.lengthScale()
-    timeScale = normalizer.timeScale()
-    velocityScale = lengthScale / timeScale
-    velocityT = self.fields.get("velocity(t)")
-    velocityT.scale(velocityScale.value)
-    velocityT.zero()
-
-    self._debug.log(resourceUsageString())
-    memoryLogger.stagePop()
-
-    # Allocates memory for nonzero pattern and Jacobian
-    self._info.log("Creating Jacobian matrix.")
-    self._setJacobianMatrixType()
-    from pylith.topology.Jacobian import Jacobian
-    self.jacobian = Jacobian(self.fields.solution(),
-                             self.matrixType, self.blockMatrixOkay)
-    self.jacobian.zero() # TEMPORARY, to get correct memory usage
-    self._debug.log(resourceUsageString())
-
-    memoryLogger.stagePush("Problem")
-    self._info.log("Initializing solver.")
-    self.solver.initialize(self.fields, self.jacobian, self)
-    self._debug.log(resourceUsageString())
-
-    # Initial time step solves for total displacement field, not increment
-    self._stepCount = 0
-    for constraint in self.constraints:
-      constraint.useSolnIncr(False)
-    for integrator in self.integratorsMesh + self.integratorsSubMesh:
-      integrator.useSolnIncr(False)
-
-    memoryLogger.stagePop()
-    memoryLogger.setDebug(0)
-    return
-
-
-  def getStartTime(self):
-    """
-    Get time at which time stepping should start.
-    """
-    dt = self.timeStep.timeStep(self.mesh,
-                                self.integratorsMesh + self.integratorsSubMesh)
-    return -dt
-
-
-  def prestep(self, t, dt):
-    """
-    Hook for doing stuff before advancing time step.
-    """
-    
-    # If finishing first time step, then switch from solving for total
-    # displacements to solving for incremental displacements
-    needNewJacobian = False
-    if 1 == self._stepCount:
-      self._info.log("Switching from total field solution to incremental " \
-                     "field solution.")
-      for constraint in self.constraints:
-        constraint.useSolnIncr(True)
-      for integrator in self.integratorsMesh + self.integratorsSubMesh:
-        integrator.useSolnIncr(True)
-      needNewJacobian = True
-
-    self._info.log("Setting constraints.")
-    dispIncr = self.fields.get("dispIncr(t->t+dt)")
-    dispIncr.zero()
-    if 0 == self._stepCount:
-      for constraint in self.constraints:
-        constraint.setField(t+dt, dispIncr)
-    else:
-      for constraint in self.constraints:
-        constraint.setFieldIncr(t, t+dt, dispIncr)
-
-    for integrator in self.integratorsMesh + self.integratorsSubMesh:
-      integrator.timeStep(dt)
-      if integrator.needNewJacobian():
-        needNewJacobian = True
-    if needNewJacobian:
-      self._reformJacobian(t, dt)
-
-    return
-
-
-  def step(self, t, dt):
-    """
-    Advance to next time step.
-    """
-    dispIncr = self.fields.get("dispIncr(t->t+dt)")
-
-    self._reformResidual(t+dt, dt)
-
-    self._info.log("Solving equations.")
-    residual = self.fields.get("residual")
-    self._eventLogger.stagePush("Solve")
-    #self.jacobian.view() # TEMPORARY
-    self.solver.solve(dispIncr, self.jacobian, residual)
-    #dispIncr.view("DISP INCR") # TEMPORARY
-
-    # DEBUGGING Verify solution makes residual 0
-    #self._reformResidual(t+dt, dt)
-    #residual.view("RESIDUAL")
-    
-    self._eventLogger.stagePop()
-
-    return
-
-
-  def poststep(self, t, dt):
-    """
-    Hook for doing stuff after advancing time step.
-    """
-    # Update displacement field from time t to time t+dt.
-    dispIncr = self.fields.get("dispIncr(t->t+dt)")
-    disp = self.fields.get("disp(t)")
-    disp += dispIncr
-    dispIncr.zero()
-
-    # Complete post-step processing, then write data.
-    Formulation.poststep(self, t, dt)
-
-    # Write data. Velocity at time t will be based upon displacement
-    # at time t-dt and t.
-    self._info.log("Writing solution fields.")
-    for output in self.output.components():
-      output.writeData(t+dt, self.fields)
-    self._writeData(t+dt)
-
-    self._stepCount += 1
-    return
-
-
-  def finalize(self):
-    """
-    Cleanup after time stepping.
-    """
-    Formulation.finalize(self)
-    return
-
-
-  # PRIVATE METHODS ////////////////////////////////////////////////////
-
-  def _configure(self):
-    """
-    Set members based using inventory.
-    """
-    Formulation._configure(self)
-
-    import journal
-    self._debug = journal.debug(self.name)
-    return
-
-
-# FACTORIES ////////////////////////////////////////////////////////////
-
-def pde_formulation():
-  """
-  Factory associated with Implicit.
-  """
-  return Implicit()
-
-
-# End of file 

Modified: short/3D/PyLith/trunk/pylith/problems/Solver.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Solver.py	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/pylith/problems/Solver.py	2011-09-20 16:36:13 UTC (rev 18949)
@@ -24,6 +24,16 @@
 
 from pylith.utils.PetscComponent import PetscComponent
 
+# VALIDATORS ///////////////////////////////////////////////////////////
+
+# Validate use of CUDA.
+def validateUseCUDA(value):
+  from pylith.utils.utils import isCUDAEnabled
+  if value and not isCUDAEnabled:
+    raise ValueError("PyLith is not built with CUDA support.")
+  return value
+
+
 # Solver class
 class Solver(PetscComponent):
   """
@@ -43,13 +53,18 @@
     ## Python object for managing Solver facilities and properties.
     ##
     ## \b Properties
-    ## @li None
+    ## @li \b use_cuda Use CUDA in solve if supported by solver.
     ##
     ## \b Facilities
     ## @li None
 
     import pyre.inventory
 
+    useCUDA = pyre.inventory.bool("use_cuda", default=False,
+                                  validator=validateUseCUDA)
+    useCUDA.meta['tip'] = "Enable use of CUDA for finite-element integrations."
+
+
   # PUBLIC METHODS /////////////////////////////////////////////////////
 
   def __init__(self, name="solver"):
@@ -60,6 +75,14 @@
     return
 
 
+  def preinitialize(self):
+    if self.useCUDA:
+      from pylith.utils.petsc import optionsSetValue
+      optionsSetValue("-vec_type", "mpicusp")
+      optionsSetValue("-mat_type", "mpiaijcusp")
+    return
+
+
   # PRIVATE METHODS /////////////////////////////////////////////////////
 
   def _configure(self):
@@ -67,6 +90,8 @@
     Set members based using inventory.
     """
     PetscComponent._configure(self)
+
+    self.useCUDA = self.inventory.useCUDA
     return
 
 

Modified: short/3D/PyLith/trunk/unittests/pytests/utils/TestPetscManager.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/utils/TestPetscManager.py	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/unittests/pytests/utils/TestPetscManager.py	2011-09-20 16:36:13 UTC (rev 18949)
@@ -27,15 +27,16 @@
   """
   Unit testing of PetscManager object.
   """
-  
 
   def test_initfinal(self):
     """
     Test initialize/finalize.
     """
     from pylith.utils.PetscManager import PetscManager
+    from pylith.utils.petsc import optionsSetValue
     manager = PetscManager()
     manager.initialize()
+    optionsSetValue("-vec_type", "seq")
     manager.finalize()
     return
 

Modified: short/3D/PyLith/trunk/unittests/pytests/utils/TestPylith.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/utils/TestPylith.py	2011-09-20 07:03:21 UTC (rev 18948)
+++ short/3D/PyLith/trunk/unittests/pytests/utils/TestPylith.py	2011-09-20 16:36:13 UTC (rev 18949)
@@ -30,9 +30,18 @@
   """
   
 
+  def test_sizeofVoidPtr(self):
+    """
+    Test sizeofVoidPtr().
+    """
+    from pylith.utils.utils import sizeofVoidPtr
+    size = sizeofVoidPtr()
+    return
+
+
   def test_sizeofPylithScalar(self):
     """
-    Test constructor.
+    Test sizeofPylithScalar().
     """
     from pylith.utils.utils import sizeofPylithScalar
     size = sizeofPylithScalar()
@@ -40,4 +49,13 @@
     return
 
 
+  def test_isCUDAEnabled(self):
+    """
+    Test constructor.
+    """
+    from pylith.utils.utils import isCUDAEnabled
+    value = isCUDAEnabled()
+    return
+
+
 # End of file 



More information about the CIG-COMMITS mailing list