[cig-commits] r17889 - in short/3D/PyLith/branches/pylith-scecdynrup: . libsrc/materials libsrc/topology modulesrc/mpi pylith/apps
brad at geodynamics.org
brad at geodynamics.org
Thu Feb 17 07:58:19 PST 2011
Author: brad
Date: 2011-02-17 07:58:18 -0800 (Thu, 17 Feb 2011)
New Revision: 17889
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/TODO
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/materials/DruckerPrager3D.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/topology/RefineVol8Face4Edges2.cc
short/3D/PyLith/branches/pylith-scecdynrup/merge.sh
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/mpi.i
short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PetscApplication.py
Log:
Merge from trunk.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/TODO 2011-02-17 02:23:25 UTC (rev 17888)
+++ short/3D/PyLith/branches/pylith-scecdynrup/TODO 2011-02-17 15:58:18 UTC (rev 17889)
@@ -26,26 +26,40 @@
Add time dataset to vertex_fields and cell_fields.
- Add creation of .xmf metadata file for VTK.
+ Add creation of .xmf metadata file for ParaView/Visit.
Will this work for [#timesteps, #points, fiberDim]?
-* nondimensional viscosity in explicit (0.1 is what is in SPECFEM)
-
* preprocess mesh
+ Adjust topology, reorder cells, partition
+ Write partitioned mesh.
+
* Reimplement Fields to use a single section.
Fields
FieldsNew -> CombinedFields
- Use CombinedFields for acc, vel, disp(t+dt), disp(t), etc?
- Use Field for dispIncr(t->t+dt)?
+ SolutionFields
+ Use CombinedFields for acc, vel, disp(t+dt), disp(t), etc?
+ Use Field for dispIncr(t->t+dt), residual(t)
* Cleanup
+ memory model
+ full-scale testing
+ Code cleanup
+----------------------------------------------------------------------
+SECONDARY PRIORITIES
+----------------------------------------------------------------------
+
+* Paper
+
+ General paper - focus on fault implementation
+
+ Custom preconditioner
+
+* Cleanup
+
Add elasticPrestep() to Formulation
Remove solnIncr, keep setField()
@@ -58,17 +72,6 @@
+ Refactor friction sensitivity solve, fault preconditioner, and
adjustSolnLumped()
-
-----------------------------------------------------------------------
-SECONDARY PRIORITIES
-----------------------------------------------------------------------
-
-* Paper
-
- General paper - focus on fault implementation
-
- Custom preconditioner
-
* Optimization
+ inline methods (what isn't getting inlined) -Winline
+ Specialized elasticity integrator objects for Quad4, Hex8
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/materials/DruckerPrager3D.cc 2011-02-17 02:23:25 UTC (rev 17888)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/materials/DruckerPrager3D.cc 2011-02-17 15:58:18 UTC (rev 17889)
@@ -32,6 +32,7 @@
#include <cassert> // USES assert()
#include <cstring> // USES memcpy()
#include <sstream> // USES std::ostringstream
+#include <iostream> // USES std::cout
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
@@ -545,11 +546,20 @@
strainPPTpdt[4]/ae + devStressInitial[4],
strainPPTpdt[5]/ae + devStressInitial[5]};
const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
- const double yieldFunction = 3.0* alphaYield * trialMeanStress +
+ const double stressInvar2 =
sqrt(0.5 *
- pylith::materials::ElasticMaterial::scalarProduct3D( trialDevStress,
- trialDevStress)) -
- beta;
+ pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
+ trialDevStress));
+ const double yieldFunction = 3.0 * alphaYield * trialMeanStress +
+ stressInvar2 - beta;
+#if 0 // DEBUGGING
+ std::cout << "Function _calcStressElastoPlastic: elastic" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " trialMeanStress: " << trialMeanStress << std::endl;
+ std::cout << " stressInvar2: " << stressInvar2 << std::endl;
+ std::cout << " yieldFunction: " << yieldFunction << std::endl;
+#endif
PetscLogFlops(76);
// If yield function is greater than zero, compute elastoplastic stress.
@@ -558,11 +568,13 @@
pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
devStressInitial);
const double strainPPTpdtProd =
- pylith::materials::ElasticMaterial::scalarProduct3D(strainPPTpdt, strainPPTpdt);
- const double d = sqrt(ae * ae * devStressInitialProd +
- 2.0 * ae *
- pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial, strainPPTpdt) +
- strainPPTpdtProd);
+ pylith::materials::ElasticMaterial::scalarProduct3D(strainPPTpdt,
+ strainPPTpdt);
+ const double d =
+ sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
+ pylith::materials::ElasticMaterial::scalarProduct3D(devStressInitial,
+ strainPPTpdt) +
+ strainPPTpdtProd);
const double plasticMult = 2.0 * ae * am *
(3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
(6.0 * alphaYield * alphaFlow * ae + am);
@@ -623,6 +635,31 @@
PetscLogFlops(31);
} // else
+#if 0 // DEBUGGING
+ const double alphaYield = properties[p_alphaYield];
+ const double beta = properties[p_beta];
+ const double alphaFlow = properties[p_alphaFlow];
+ const double meanStressTest = (stress[0] + stress[1] + stress[2])/3.0;
+ const double devStressTest[] = { stress[0] - meanStressTest,
+ stress[1] - meanStressTest,
+ stress[2] - meanStressTest,
+ stress[3],
+ stress[4],
+ stress[5]};
+ const double stressInvar2Test =
+ sqrt(0.5 *
+ pylith::materials::ElasticMaterial::scalarProduct3D(devStressTest,
+ devStressTest));
+
+ const double yieldFunctionTest = 3.0 * alphaYield * meanStressTest +
+ stressInvar2Test - beta;
+ std::cout << "Function _calcStressElastoPlastic: end" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " meanStressTest: " << meanStressTest << std::endl;
+ std::cout << " stressInvar2Test: " << stressInvar2Test << std::endl;
+ std::cout << " yieldFunctionTest: " << yieldFunctionTest << std::endl;
+#endif
} // _calcStressElastoplastic
@@ -808,11 +845,20 @@
strainPPTpdt[4]/ae + devStressInitial[4],
strainPPTpdt[5]/ae + devStressInitial[5]};
const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
- const double yieldFunction = 3.0* alphaYield * trialMeanStress +
+ const double stressInvar2 =
sqrt(0.5 *
pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
- trialDevStress)) -
- beta;
+ trialDevStress));
+ const double yieldFunction = 3.0 * alphaYield * trialMeanStress +
+ stressInvar2 - beta;
+#if 0 // DEBUGGING
+ std::cout << "Function _calcElasticConstsElastoPlastic:" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " trialMeanStress: " << trialMeanStress << std::endl;
+ std::cout << " stressInvar2: " << stressInvar2 << std::endl;
+ std::cout << " yieldFunction: " << yieldFunction << std::endl;
+#endif
PetscLogFlops(76);
// If yield function is greater than zero, compute elastoplastic stress and
@@ -1067,11 +1113,20 @@
strainPPTpdt[4]/ae + devStressInitial[4],
strainPPTpdt[5]/ae + devStressInitial[5]};
const double trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
- const double yieldFunction = 3.0* alphaYield * trialMeanStress +
+ const double stressInvar2 =
sqrt(0.5 *
pylith::materials::ElasticMaterial::scalarProduct3D(trialDevStress,
- trialDevStress)) -
- beta;
+ trialDevStress));
+ const double yieldFunction = 3.0 * alphaYield * trialMeanStress +
+ stressInvar2 - beta;
+#if 0 // DEBUGGING
+ std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
+ std::cout << " alphaYield: " << alphaYield << std::endl;
+ std::cout << " beta: " << beta << std::endl;
+ std::cout << " trialMeanStress: " << trialMeanStress << std::endl;
+ std::cout << " stressInvar2: " << stressInvar2 << std::endl;
+ std::cout << " yieldFunction: " << yieldFunction << std::endl;
+#endif
PetscLogFlops(76);
// If yield function is greater than zero, compute plastic strains.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/topology/RefineVol8Face4Edges2.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/topology/RefineVol8Face4Edges2.cc 2011-02-17 02:23:25 UTC (rev 17888)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/topology/RefineVol8Face4Edges2.cc 2011-02-17 15:58:18 UTC (rev 17889)
@@ -294,7 +294,6 @@
const Obj<mesh_type>& oldMesh,
const MeshOrder& orderOldMesh)
{ // overlapAddNewVertices
-#if 0
assert(!newMesh.isNull());
assert(!oldMesh.isNull());
@@ -334,10 +333,10 @@
std::insert_iterator<std::set<int> >(ranks, ranks.begin()));
if(ranks.size()) {
- newVerticesSection->addFiberDimension(std::min(e_iter->first.first, e_iter->first.second)+localOffset, 1);
- for(std::set<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
- bndryEdgeToRank[e_iter->first].push_back(*r_iter);
- } // for
+ newVerticesSection->addFiberDimension(std::min(e_iter->first.first, e_iter->first.second)+localOffset, 1);
+ for(std::set<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
+ bndryEdgeToRank[e_iter->first].push_back(*r_iter);
+ } // for
} // if
} // if
} // for
@@ -414,22 +413,42 @@
int v = 0;
value_type* values = (dim > 0) ? new value_type[dim] : 0;
-#if 0
for(std::map<FaceType, std::vector<int> >::const_iterator f_iter = bndryFaceToRank.begin(); f_iter != bndryFaceToRank.end() && v < dim; ++f_iter) {
- if (std::min(f_iter->first.first, f_iter->first.second)+localOffset == p) {
- values[v++] = FaceType(std::max(f_iter->first.first, f_iter->first.second)+localOffset, _faceToVertex[f_iter->first]);
+ const point_type first = f_iter->first.points[0];
+ point_type minVertex = first;
+
+ for(int i = 1; i < 4; ++i) {
+ const point_type nextVertex = f_iter->first.points[i];
+ minVertex = std::min(nextVertex, minVertex);
+ }
+
+ if (minVertex+localOffset == p) {
+ FaceType face;
+ int k = 0;
+
+ for(int i = 0; i < 4; ++i) {
+ if (f_iter->first.points[i] != minVertex) {
+ face.points[k++] = f_iter->first.points[i];
+ }
+ }
+ assert(k == 3);
+ std::sort(&face.points[0], &face.points[3]);
+ face.points[3] = _faceToVertex[f_iter->first];
+ values[v++] = face;
} // if
} // for
-#endif
+ assert(v == dim);
newFaceVerticesSection->updatePoint(p, values);
delete [] values;
} // for
// Copy across overlap
typedef ALE::Pair<int, point_type> overlap_point_type;
- Obj<ALE::Section<overlap_point_type, EdgeType> > overlapVertices = new ALE::Section<overlap_point_type, EdgeType>(oldMesh->comm());
+ Obj<ALE::Section<overlap_point_type, EdgeType> > overlapVertices = new ALE::Section<overlap_point_type, EdgeType>(oldMesh->comm());
+ Obj<ALE::Section<overlap_point_type, FaceType> > overlapFaceVertices = new ALE::Section<overlap_point_type, FaceType>(oldMesh->comm());
- ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
+ ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
+ ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newFaceVerticesSection, overlapFaceVertices);
// Merge by translating edge to local points, finding edge in _edgeToVertex, and adding (local new vetex, remote new vertex) to overlap
for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end(); ++e_iter) {
const point_type localPoint = _edgeToVertex[e_iter->first];
@@ -440,43 +459,84 @@
const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanks = newSendOverlap->support(e_iter->first.first+localOffset);
for(mesh_type::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks->begin(); lr_iter != leftRanks->end(); ++lr_iter) {
- if (rank == *lr_iter) {
- remoteLeft = lr_iter.color();
- break;
- } // if
+ if (rank == *lr_iter) {
+ remoteLeft = lr_iter.color();
+ break;
+ } // if
} // for
const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = newSendOverlap->support(e_iter->first.second+localOffset);
for(mesh_type::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks->begin(); rr_iter != rightRanks->end(); ++rr_iter) {
- if (rank == *rr_iter) {
- remoteRight = rr_iter.color();
- break;
- } // if
+ if (rank == *rr_iter) {
+ remoteRight = rr_iter.color();
+ break;
+ } // if
} // for
const point_type remoteMin = std::min(remoteLeft, remoteRight);
const point_type remoteMax = std::max(remoteLeft, remoteRight);
const int remoteSize = overlapVertices->getFiberDimension(overlap_point_type(rank, remoteMin));
- const EdgeType *remoteVals = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
+ const EdgeType *remoteVals = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
point_type remotePoint = -1;
for(int d = 0; d < remoteSize; ++d) {
- if (remoteVals[d].first == remoteMax) {
- remotePoint = remoteVals[d].second;
- break;
- } // if
+ if (remoteVals[d].first == remoteMax) {
+ remotePoint = remoteVals[d].second;
+ break;
+ } // if
} // for
newSendOverlap->addArrow(localPoint, rank, remotePoint);
newRecvOverlap->addArrow(rank, localPoint, remotePoint);
} // for
} // for
+ // Merge by translating face to local points, finding face in _faceToVertex, and adding (local new vetex, remote new vertex) to overlap
+ for(std::map<FaceType, std::vector<int> >::const_iterator f_iter = bndryFaceToRank.begin(); f_iter != bndryFaceToRank.end(); ++f_iter) {
+ const point_type localPoint = _faceToVertex[f_iter->first];
+
+ for(std::vector<int>::const_iterator r_iter = f_iter->second.begin(); r_iter != f_iter->second.end(); ++r_iter) {
+ FaceType remoteVertices(-1);
+ const int rank = *r_iter;
-#if 1
+ for(int i = 0; i < 4; ++i) {
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& faceRanks = newSendOverlap->support(f_iter->first.points[i]+localOffset);
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator fr_iter = faceRanks->begin(); fr_iter != faceRanks->end(); ++fr_iter) {
+ if (rank == *fr_iter) {
+ remoteVertices.points[i] = fr_iter.color();
+ break;
+ } // if
+ } // for
+ }
+ const point_type remoteMin = std::min(std::min(std::min(remoteVertices.points[0], remoteVertices.points[1]), remoteVertices.points[2]), remoteVertices.points[3]);
+ const int remoteSize = overlapFaceVertices->getFiberDimension(overlap_point_type(rank, remoteMin));
+ const FaceType *remoteVals = overlapFaceVertices->restrictPoint(overlap_point_type(rank, remoteMin));
+ point_type remotePoint = -1;
+ int k = 0;
+ FaceType remoteMax;
+
+ for(int i = 0; i < 4; ++i) {
+ if (remoteVertices.points[i] == remoteMin) continue;
+ remoteMax.points[k++] = remoteVertices.points[i];
+ }
+ assert(k == 3);
+ std::sort(&remoteMax.points[0], &remoteMax.points[3]);
+ for(int d = 0; d < remoteSize; ++d) {
+ int i = 0;
+
+ for(i = 0; i < 3; ++i) {
+ if (remoteVals[d].points[i] != remoteMax.points[i]) break;
+ }
+ if (i == 3) {
+ remotePoint = remoteVals[d].points[3];
+ break;
+ } // if
+ } // for
+ newSendOverlap->addArrow(localPoint, rank, remotePoint);
+ newRecvOverlap->addArrow(rank, localPoint, remotePoint);
+ } // for
+ } // for
+
oldSendOverlap->view("OLD SEND OVERLAP");
oldRecvOverlap->view("OLD RECV OVERLAP");
newSendOverlap->view("NEW SEND OVERLAP");
newRecvOverlap->view("NEW RECV OVERLAP");
-#endif
-
-#endif
} // overlapAddNewVertces
Modified: short/3D/PyLith/branches/pylith-scecdynrup/merge.sh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/merge.sh 2011-02-17 02:23:25 UTC (rev 17888)
+++ short/3D/PyLith/branches/pylith-scecdynrup/merge.sh 2011-02-17 15:58:18 UTC (rev 17889)
@@ -1 +1,2 @@
+svn up
svn merge -r 16682:HEAD svn+ssh://svn@geodynamics.org/cig/short/3D/PyLith/trunk .
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/Makefile.am 2011-02-17 02:23:25 UTC (rev 17888)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/Makefile.am 2011-02-17 15:58:18 UTC (rev 17889)
@@ -25,7 +25,8 @@
swig_sources = \
mpi.i \
- mpi_comm.i
+ mpi_comm.i \
+ mpi_error.i
swig_generated = \
mpi_wrap.cxx \
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/mpi.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/mpi.i 2011-02-17 02:23:25 UTC (rev 17888)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/mpi/mpi.i 2011-02-17 15:58:18 UTC (rev 17889)
@@ -28,6 +28,8 @@
// Interfaces
%include "mpi_comm.i"
+%include "mpi_error.i"
+
// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PetscApplication.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PetscApplication.py 2011-02-17 02:23:25 UTC (rev 17888)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PetscApplication.py 2011-02-17 15:58:18 UTC (rev 17889)
@@ -59,13 +59,30 @@
Run the application in parallel on the compute nodes.
"""
self.petsc.initialize()
+
+ # Test for CUDA
try:
import pycuda
import pycuda.autoinit
self._info.log('Initialized CUDA')
except ImportError:
self._info.log('Could not initialize CUDA')
- self.main(*args, **kwds)
+ # :TODO: Set some flag here to disable CUDA. Optionally also
+ # allow user to turn off use of CUDA.
+
+ try:
+
+ self.main(*args, **kwds)
+
+ except Exception, err:
+ self.cleanup() # Attempt to clean up memory.
+ print "Fatal error while running PyLith:"
+ print err
+ print "Calling MPI_Abort() to abort PyLith application."
+ from pylith.mpi import mpi
+ errorCode = -1
+ mpi.mpi_abort(mpi.petsc_comm_world(), errorCode)
+
self.cleanup()
self.petsc.finalize()
return
More information about the CIG-COMMITS
mailing list