[cig-commits] r18790 - in short/3D/PyLith/branches/v1.6-stable: examples/bar_shearwave/tri3/tmp libsrc/pylith/faults libsrc/pylith/feassemble libsrc/pylith/materials libsrc/pylith/meshio libsrc/pylith/topology

brad at geodynamics.org brad at geodynamics.org
Wed Jul 20 20:22:30 PDT 2011


Author: brad
Date: 2011-07-20 20:22:30 -0700 (Wed, 20 Jul 2011)
New Revision: 18790

Modified:
   short/3D/PyLith/branches/v1.6-stable/examples/bar_shearwave/tri3/tmp/pylithapp.cfg
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/IntegratorElasticity.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/materials/Material.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/CellFilterAvg.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriter.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
   short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/topology/Field.cc
Log:
Fixed consistency of creating scatters in parallel (requires synchronization).

Modified: short/3D/PyLith/branches/v1.6-stable/examples/bar_shearwave/tri3/tmp/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/examples/bar_shearwave/tri3/tmp/pylithapp.cfg	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/examples/bar_shearwave/tri3/tmp/pylithapp.cfg	2011-07-21 03:22:30 UTC (rev 18790)
@@ -96,6 +96,7 @@
 id = 1
 
 # Spatial database with physical properties for elastic material
+db_properties.label = Elastic properties
 db_properties.iohandler.filename = ../matprops.spatialdb
 
 # Set the basis functions and quadrature:
@@ -182,12 +183,15 @@
 [pylithapp.timedependent.interfaces.fault.eq_srcs.rupture.slip_function]
 
 # Database specifying the final slip.
+slip.label = Final slip
 slip.iohandler.filename = ../shearwave_slip.spatialdb
 
 # Database specifying rise time.
+rise_time.label = Rise time
 rise_time.iohandler.filename = ../shearwave_risetime.spatialdb
 
 # Database specifying time at which slip begins at each point.
+slip_time.label = Slip time
 slip_time.iohandler.filename = ../shearwave_sliptime.spatialdb
 
 # ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -1035,9 +1035,9 @@
     _allocateBufferVectorField();
     topology::Field<topology::SubMesh>& buffer =
         _fields->get("buffer (vector)");
-    buffer.copy(dirSection);
     buffer.label("strike_dir");
     buffer.scale(1.0);
+    buffer.copy(dirSection);
     return buffer;
 
   } else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
@@ -1049,9 +1049,9 @@
     _allocateBufferVectorField();
     topology::Field<topology::SubMesh>& buffer =
         _fields->get("buffer (vector)");
-    buffer.copy(dirSection);
     buffer.label("dip_dir");
     buffer.scale(1.0);
+    buffer.copy(dirSection);
     return buffer;
 
   } else if (0 == strcasecmp("normal_dir", name)) {
@@ -1065,9 +1065,9 @@
     _allocateBufferVectorField();
     topology::Field<topology::SubMesh>& buffer =
         _fields->get("buffer (vector)");
-    buffer.copy(dirSection);
     buffer.label("normal_dir");
     buffer.scale(1.0);
+    buffer.copy(dirSection);
     return buffer;
 
   } else if (0 == strcasecmp("initial_traction", name)) {
@@ -1097,11 +1097,12 @@
     throw std::runtime_error(msg.str());
   } // else
 
-
   // Satisfy return values
   assert(0 != _fields);
   const topology::Field<topology::SubMesh>& buffer = _fields->get(
     "buffer (vector)");
+  throw std::logic_error("Internal error.");
+
   return buffer;
 } // vertexField
 

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -1976,6 +1976,7 @@
   const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
   buffer.cloneSection(slip);
   buffer.zero();
+  assert(buffer.vectorFieldType() == topology::FieldBase::VECTOR);
 
   logger.stagePop();
 } // _allocateBufferVectorField
@@ -1999,6 +2000,7 @@
   const topology::Field<topology::SubMesh>& area = _fields->get("area");
   buffer.cloneSection(area);
   buffer.zero();
+  assert(buffer.vectorFieldType() == topology::FieldBase::SCALAR);
 
   logger.stagePop();
 } // _allocateBufferScalarField

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/IntegratorElasticity.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/IntegratorElasticity.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -400,6 +400,7 @@
   
   // Return tensor section to satisfy member function definition. Code
   // should never get here.
+  throw std::logic_error("Internal error.");
   topology::Field<topology::Mesh>& buffer = 
     _outputFields->get("buffer (tensor)");    
 

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/materials/Material.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/materials/Material.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -315,7 +315,8 @@
 // ----------------------------------------------------------------------
 // Get physical property or state variable field.
 void
-pylith::materials::Material::getField(topology::Field<topology::Mesh> *field, const char* name) const
+pylith::materials::Material::getField(topology::Field<topology::Mesh> *field,
+				      const char* name) const
 { // getField
   // Logging of allocation is handled by getField() caller since it
   // manages the memory for the field argument.

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/CellFilterAvg.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/CellFilterAvg.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -128,11 +128,6 @@
     _fieldAvg->newSection(fieldIn, fiberDim);
     _fieldAvg->allocate();
   } else if (_fieldAvg->sectionSize() != cells->size()*fiberDim) {
-#if 1 // :BUG: Avoid memory leak in section->clear()
-    _fieldAvg->clear();
-#else
-    delete _fieldAvg; _fieldAvg = new field_type(fieldIn.mesh());
-#endif
     _fieldAvg->newSection(fieldIn, fiberDim);
     _fieldAvg->allocate();
   } // else

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriter.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriter.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriter.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -56,8 +56,13 @@
 
   const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
   assert(!sieveMesh.isNull());
-  _context = std::string("output_") + sieveMesh->getName() + 
-    ((label) ? std::string("_") + std::string(label) : std::string(""));
+
+  ostringstream s;
+  s << "output_"
+    << sieveMesh->getName();
+  if (label)
+    s << "_" << label << labelId;
+  _context = s.str();
 } // open
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -201,14 +201,12 @@
     err = VecDestroy(&elemVec); CHECK_PETSC_ERROR(err);
     delete[] tmpVertices; tmpVertices = 0;
 
-    if (!rank) {
-      hid_t h5 = -1;
-      err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
-      assert(h5 >= 0);
-      const int cellDim = mesh.dimension();
-      HDF5::writeAttribute(h5, "/topology/cells", "cell_dim", (void*)&cellDim,
-			   H5T_NATIVE_INT);
-    } // if
+    hid_t h5 = -1;
+    err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
+    assert(h5 >= 0);
+    const int cellDim = mesh.dimension();
+    HDF5::writeAttribute(h5, "/topology/cells", "cell_dim", (void*)&cellDim,
+			 H5T_NATIVE_INT);
   } catch (const std::exception& err) {
     std::ostringstream msg;
     msg << "Error while opening HDF5 file "
@@ -235,12 +233,16 @@
   _timesteps.clear();
   _tstampIndex = 0;
 
-  Xdmf metafile;
-  const std::string& hdf5filename = _hdf5Filename();
-  const int indexExt = hdf5filename.find(".h5");
-  std::string xdmfFilename = 
-    std::string(hdf5filename, 0, indexExt) + ".xmf";
-  metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  int rank = 0;
+  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
+  if (!rank) {
+    Xdmf metafile;
+    const std::string& hdf5filename = _hdf5Filename();
+    const int indexExt = hdf5filename.find(".h5");
+    std::string xdmfFilename = 
+      std::string(hdf5filename, 0, indexExt) + ".xmf";
+    metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  } // if
 } // close
 
 // ----------------------------------------------------------------------
@@ -312,7 +314,7 @@
     err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
     err = VecSetBlockSize(vector, blockSize); CHECK_PETSC_ERROR(err);
 
-    if (!rank && 0 == istep) {
+    if (0 == istep) {
       hid_t h5 = -1;
       err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
       assert(h5 >= 0);
@@ -366,7 +368,6 @@
     field.scatterSectionToVector(context);
     PetscVec vector = field.vector(context);
     assert(vector);
-
     const ALE::Obj<typename mesh_type::RealSection>& section = field.section();
     assert(!section.isNull());      
     assert(!sieveMesh->getLabelStratum(labelName, depth).isNull());
@@ -401,7 +402,7 @@
     err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
     err = VecSetBlockSize(vector, blockSize); CHECK_PETSC_ERROR(err);
 
-    if (!rank && 0 == istep) {
+    if (0 == istep) {
       hid_t h5 = -1;
       err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
       assert(h5 >= 0);
@@ -451,7 +452,7 @@
 { // _writeTimeStamp
   PetscErrorCode err = 0;
 
-  if (!rank) {
+  if (0 == rank) {
     err = VecSetValue(_tstamp, 0, t, INSERT_VALUES); CHECK_PETSC_ERROR(err);
   } // if
   err = VecAssemblyBegin(_tstamp); CHECK_PETSC_ERROR(err);

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -261,11 +261,15 @@
   _tstampIndex = 0;
   deallocate();
 
-  Xdmf metafile;
-  const std::string& hdf5filename = _hdf5Filename();
-  const int indexExt = hdf5filename.find(".h5");
-  std::string xdmfFilename = std::string(hdf5filename, 0, indexExt) + ".xmf";
-  metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  int rank = 0;
+  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
+  if (!rank) {
+    Xdmf metafile;
+    const std::string& hdf5filename = _hdf5Filename();
+    const int indexExt = hdf5filename.find(".h5");
+    std::string xdmfFilename = std::string(hdf5filename, 0, indexExt) + ".xmf";
+    metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  } // if
 } // close
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/topology/Field.cc	2011-07-20 23:04:58 UTC (rev 18789)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/topology/Field.cc	2011-07-21 03:22:30 UTC (rev 18790)
@@ -87,6 +87,7 @@
       err = VecDestroy(&s_iter->second.scatterVec);CHECK_PETSC_ERROR(err);
     } // if
   } // for
+  _scatters.clear();
 } // deallocate
 
 // ----------------------------------------------------------------------
@@ -146,6 +147,9 @@
 void
 pylith::topology::Field<mesh_type, section_type>::newSection(void)
 { // newSection
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
 
@@ -167,8 +171,12 @@
 { // newSection
   typedef typename mesh_type::SieveMesh::point_type point_type;
 
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
+
   if (fiberDim < 0) {
     std::ostringstream msg;
     msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
@@ -203,6 +211,9 @@
 { // newSection
   typedef typename mesh_type::SieveMesh::point_type point_type;
 
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
   if (fiberDim < 0) {
@@ -261,6 +272,9 @@
 pylith::topology::Field<mesh_type, section_type>::newSection(const Field& src,
 					       const int fiberDim)
 { // newSection
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
 
@@ -299,21 +313,22 @@
 void
 pylith::topology::Field<mesh_type, section_type>::cloneSection(const Field& src)
 { // cloneSection
+  std::string origLabel = _metadata.label;
+
+  // Clear memory
+  clear();
+
+  const ALE::Obj<section_type>& srcSection = src.section();
+  if (!srcSection.isNull() && _section.isNull()) {
+    newSection();
+  } // if
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
 
-  deallocate();
-  std::string origLabel = _metadata.label;
   _metadata = src._metadata;
   label(origLabel.c_str());
 
-  const ALE::Obj<section_type>& srcSection = src.section();
-  if (!srcSection.isNull() && _section.isNull()) {
-    logger.stagePop();
-    newSection();
-    logger.stagePush("Field");
-  }
-
   if (!_section.isNull()) {
     if (!srcSection->sharedStorage()) {
       _section->setAtlas(srcSection->getAtlas());
@@ -1006,14 +1021,17 @@
   err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);  
 
 #if 0
-  std::cout << "CONTEXT: " << context 
+  std::cout << "["<<sieveMesh->commRank()<<"] CONTEXT: " << context 
 	    << ", orderLabel: " << orderLabel
 	    << ", section size w/BC: " << _section->sizeWithBC()
 	    << ", section size: " << _section->size()
 	    << ", section storage size: " << _section->getStorageSize()
 	    << ", global numbering size: " << numbering->getGlobalSize()
-	    << ", global size: " << order->getGlobalSize()
+	    << ", global order size: " << order->getGlobalSize()
+	    << ", local numbering size: " << numbering->getLocalSize()
+	    << ", local order size: " << order->getLocalSize()
 	    << ", scatter from size: " << sinfo.scatter->from_n
+	    << ", scatter: " << sinfo.scatter
 	    << std::endl;
 #endif
   
@@ -1139,8 +1157,33 @@
 { // _getScatter
   assert(context);
 
-  const bool isNewScatter = _scatters.find(context) == _scatters.end();
+  bool isNewScatter = _scatters.find(context) == _scatters.end();
 
+  // Synchronize creation of scatter (empty sections may have
+  // leftover, reusable scatters that need to be cleared out).
+  int numNewScatterLocal = (isNewScatter) ? 1 : 0;
+  int numNewScatter = 0;
+  MPI_Allreduce(&numNewScatterLocal, &numNewScatter, 1, MPI_INT, MPI_MAX,
+		_mesh.comm());
+  if (numNewScatter && !isNewScatter) {
+    // remove old scatter
+    ScatterInfo& sinfo = _scatters[context];
+    PetscErrorCode err = 0;
+    if (sinfo.vector) {
+      err = VecDestroy(&sinfo.vector);CHECK_PETSC_ERROR(err);
+    } // if
+    if (sinfo.scatter) {
+      err = VecScatterDestroy(&sinfo.scatter);CHECK_PETSC_ERROR(err);
+    } // if
+
+    if (sinfo.scatterVec) {
+      err = VecDestroy(&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } // if
+
+    _scatters.erase(context);
+    isNewScatter = true;
+  } // if
+
   if (isNewScatter && !createOk) {
     std::ostringstream msg;
     msg << "Scatter for context '" << context << "' does not exist.";



More information about the CIG-COMMITS mailing list