[cig-commits] r13028 - in short/3D/PyLith/trunk: libsrc/feassemble modulesrc/utils pylith/problems

[brad at geodynamics.org]

Author: brad
Date: 2008-10-14 14:49:38 -0700 (Tue, 14 Oct 2008)
New Revision: 13028

Modified:
   short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
   short/3D/PyLith/trunk/modulesrc/utils/petsc.pyxe.src
   short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
Fixed assembly of Jacobian. Use MAT_FLUSH_ASSEMBLY mode for assembled portion.

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc	2008-10-14 21:00:53 UTC (rev 13027)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc	2008-10-14 21:49:38 UTC (rev 13028)
@@ -300,6 +300,20 @@
   _basisDerivPre->setFiberDimension(cells, _numQuadPts*_numBasis*_spaceDim);
   _basisDerivPre->allocatePoint();
 
+  const int ncells = cells->size();
+  const int nbytes = 
+    (
+     _numQuadPts*_spaceDim + // quadPts
+    _numQuadPts*_cellDim*_spaceDim + // jacobian
+    //_numQuadPts*_cellDim*_spaceDim + // jacobianInv
+    _numQuadPts + // jacobianDet
+     _numQuadPts*_numBasis*_spaceDim // basisDeriv
+     ) * ncells * sizeof(double);
+    
+  std::cout << "Quadrature::precomputeGeometry() allocating "
+	    << nbytes/(1024*1024) << " MB."
+	    << std::endl;
+
   for(Mesh::label_sequence::iterator c_iter = cells->begin();
       c_iter != end;
       ++c_iter) {

Modified: short/3D/PyLith/trunk/modulesrc/utils/petsc.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/utils/petsc.pyxe.src	2008-10-14 21:00:53 UTC (rev 13027)
+++ short/3D/PyLith/trunk/modulesrc/utils/petsc.pyxe.src	2008-10-14 21:49:38 UTC (rev 13028)
@@ -13,6 +13,8 @@
 #header{
 #include <petsc.h>
 #include <petscmat.h>
+#include <string.h>
+#include <stdexcept>
 #}header
 
 # ----------------------------------------------------------------------
@@ -82,21 +84,30 @@
   return
 
 
-def mat_assemble(mat):
+def mat_assemble(mat, mode):
   """
   Assemble matrix.
   """
   # create shim for 'MatAssemblyBegin/MatAssemblyEnd'
-  #embed{ int Mat_assemble(void* matVptr)
+  #embed{ int Mat_assemble(void* matVptr, char* mode)
   Mat* mat = (Mat*) matVptr;
   PetscErrorCode err = 0;
-  err = MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY); CHKERRQ(err);
-  err = MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY); CHKERRQ(err);
+  if (0 == strcmp(mode, "final_assembly")) {
+    err = MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY); CHKERRQ(err);
+    err = MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY); CHKERRQ(err);
+  } else if (0 == strcmp(mode, "flush_assembly")) {
+    err = MatAssemblyBegin(*mat, MAT_FLUSH_ASSEMBLY); CHKERRQ(err);
+    err = MatAssemblyEnd(*mat, MAT_FLUSH_ASSEMBLY); CHKERRQ(err);
+  } else {
+    throw std::runtime_error("Unknown mode");
+  }  
   return 0;
   #}embed
   cdef void* matVptr
   matVptr = PyCObject_AsVoidPtr(mat)
-  err = Mat_assemble(matVptr)
+  if mode != "final_assembly" and mode != "flush_assembly":
+    raise ValueError("Unknown mode '%s' for mat_assemble." % mode)
+  err = Mat_assemble(matVptr, mode)
   if err:
     raise RuntimError("Error assembling matrix.")
 

Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py	2008-10-14 21:00:53 UTC (rev 13027)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py	2008-10-14 21:49:38 UTC (rev 13028)
@@ -411,7 +411,6 @@
     """
     Reform Jacobian matrix for operator.
     """
-    self._info.log("Reforming Jacobian of operator.")
     self._debug.log(resourceUsageString())
     import pylith.utils.petsc as petsc
     petsc.mat_setzero(self.jacobian)
@@ -423,16 +422,17 @@
       integrator.integrateJacobianAssembled(self.jacobian, t+dt,
                                             self.fields)
       self._debug.log(resourceUsageString()) # TEMPORARY
-    self._info.log("Assembling Jacobian of operator.")
-    petsc.mat_assemble(self.jacobian)
+    self._info.log("Flushing assembly of Jacobian of operator.")
+    petsc.mat_assemble(self.jacobian, "flush_assembly")
 
     # Add in contributions that require assembly
+    self._info.log("Reforming unassembled portion of Jacobian of operator.")
     for integrator in self.integrators:
       integrator.timeStep(dt)
       integrator.integrateJacobian(self.jacobian, t+dt, self.fields)
       self._debug.log(resourceUsageString()) # TEMPORARY
-    self._info.log("Assembling Jacobian of operator.")
-    petsc.mat_assemble(self.jacobian)
+    self._info.log("Doing final assembly of Jacobian of operator.")
+    petsc.mat_assemble(self.jacobian, "final_assembly")
 
     if self.viewJacobian:
       filename = self._createJacobianFilename(t+dt)
@@ -470,8 +470,6 @@
     for integrator in self.integrators:
       integrator.timeStep(dt)
       integrator.integrateResidual(residual, t, self.fields)
-    self._info.log("Completing residual.")
-    bindings.completeSection(self.mesh.cppHandle, residual)
 
     # Add in contributions that do not require assembly
     self._info.log("Integrating assembled residual term in operator.")
@@ -479,6 +477,9 @@
       integrator.timeStep(dt)
       integrator.integrateResidualAssembled(residual, t, self.fields)
 
+    self._info.log("Completing residual.")
+    bindings.completeSection(self.mesh.cppHandle, residual)
+
     return