Feature/variable resolution damage

scripts/lc-builds/toss3_gcc8.1.0-release-mpi-cxx.sh

@@ -0,0 +1,25 @@
+#!/usr/bin/env bash
+
+SCRIPT_PATH=${0%/*}
+. "$SCRIPT_PATH/utils/parse-args.sh"
+
+# Inherit build directory name from script name
+BUILD_SUFFIX="lc_$(TMP=${BASH_SOURCE##*/}; echo ${TMP%.*})"
+
+rm -rf ${BUILD_SUFFIX} 2>/dev/null
+mkdir -p ${BUILD_SUFFIX}/install
+mkdir -p ${BUILD_SUFFIX}/build && cd ${BUILD_SUFFIX}/build
+
+module load cmake/3.14.5
+module load gcc/8.1.0
+
+cmake \
+  ${SRC_DIR} \
+  -DCMAKE_BUILD_TYPE=Release \
+  -DCMAKE_CXX_COMPILER=/usr/tce/packages/gcc/gcc-8.1.0/bin/g++ \
+  -C ${HOST_CONFIGS_DIR}/lc-builds/toss3/gcc8.1.0_tpl.cmake \
+  -DENABLE_OPENMP=On \
+  -DENABLE_MPI=On \
+  -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR \
+  -DENABLE_STATIC_CXXONLY=On \
+  $CMAKE_ARGS \

src/Boundary/ConstantBoundary.cc

@@ -58,7 +58,7 @@ ConstantBoundary(DataBase<Dimension>& dataBase,
     mNodeFlags[*itr] = 1;
   }
 
-  // // Issue a big old warning!
+  // Issue a big old warning!
   // if (Process::getRank() == 0) cerr << "WARNING: ConstantBoundary is currently not compatible with redistributing nodes!\nMake sure you don't allow redistribution with this Boundary condition." << endl;
 }
 

src/Damage/CMakeLists.txt

@@ -2,6 +2,8 @@ include_directories(.)
 set(Damage_inst
     DamageModel
     TensorDamageModel
+    ProbabilisticDamageModel
+    ProbabilisticDamagePolicy
     StrainPolicy
     TensorStrainPolicy
     TensorDamagePolicy
@@ -39,7 +41,10 @@ set(Damage_headers
     StrainPolicy.hh
     TensorDamageModel.hh
     TensorDamageModelInline.hh
+    ProbabilisticDamageModel.hh
+    ProbabilisticDamageModelInline.hh
     TensorDamagePolicy.hh
+    ProbabilisticDamagePolicy.hh
     TensorStrainPolicy.hh
     YoungsModulusPolicy.hh
     computeFragmentField.hh
@@ -53,6 +58,7 @@ spheral_install_python_files(
   WeibullDamage.py
   GradyKippTensorDamage.py
   JohnsonCookDamageFactories.py
+  ShadowProbabilisticDamageModel.py
   )
 
 spheral_add_cxx_library(Damage)

src/Damage/DamageModel.cc

@@ -7,9 +7,6 @@
 //----------------------------------------------------------------------------//
 #include "FileIO/FileIO.hh"
 #include "DamageModel.hh"
-#include "YoungsModulusPolicy.hh"
-#include "LongitudinalSoundSpeedPolicy.hh"
-#include "DamagedSoundSpeedPolicy.hh"
 #include "Strength/SolidFieldNames.hh"
 #include "NodeList/SolidNodeList.hh"
 #include "DataBase/DataBase.hh"
@@ -57,17 +54,12 @@ DamageModel<Dimension>::
 DamageModel(SolidNodeList<Dimension>& nodeList,
             const TableKernel<Dimension>& W,
             const double crackGrowthMultiplier,
-            const DamageCouplingAlgorithm damageCouplingAlgorithm,
-            const FlawStorageType& flaws):
+            const DamageCouplingAlgorithm damageCouplingAlgorithm):
   Physics<Dimension>(),
-  mFlaws(SolidFieldNames::flaws, flaws),
   mNodeList(nodeList),
   mW(W),
   mCrackGrowthMultiplier(crackGrowthMultiplier),
-  mCriticalNodesPerSmoothingScale(0.99),
   mDamageCouplingAlgorithm(damageCouplingAlgorithm),
-  mYoungsModulus(SolidFieldNames::YoungsModulus, nodeList),
-  mLongitudinalSoundSpeed(SolidFieldNames::longitudinalSoundSpeed, nodeList),
   mExcludeNode("Nodes excluded from damage", nodeList, 0),
   mNodeCouplingPtr(new NodeCoupling()),
   mComputeIntersectConnectivity(false),
@@ -97,7 +89,6 @@ computeScalarDDDt(const DataBase<Dimension>& /*dataBase*/,
 
   // Pre-conditions.
   REQUIRE(DDDt.nodeListPtr() == &mNodeList);
-  REQUIRE(mFlaws.nodeListPtr() == &mNodeList);
 
   // Get the state fields.
   const auto  clKey = State<Dimension>::buildFieldKey(SolidFieldNames::longitudinalSoundSpeed, mNodeList.name());
@@ -135,29 +126,6 @@ computeScalarDDDt(const DataBase<Dimension>& /*dataBase*/,
   END_CONTRACT_SCOPE
 }
 
-//------------------------------------------------------------------------------
-// Register our state.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-void
-DamageModel<Dimension>::
-registerState(DataBase<Dimension>& dataBase,
-              State<Dimension>& state) {
-
-  // Register Youngs modulus and the longitudinal sound speed.
-  typedef typename State<Dimension>::PolicyPointer PolicyPointer;
-  PolicyPointer EPolicy(new YoungsModulusPolicy<Dimension>());
-  PolicyPointer clPolicy(new LongitudinalSoundSpeedPolicy<Dimension>());
-  state.enroll(mYoungsModulus, EPolicy);
-  state.enroll(mLongitudinalSoundSpeed, clPolicy);
-
-  // Set the initial values for the Youngs modulus, sound speed, and pressure.
-  typename StateDerivatives<Dimension>::PackageList dummyPackages;
-  StateDerivatives<Dimension> derivs(dataBase, dummyPackages);
-  EPolicy->update(state.key(mYoungsModulus), state, derivs, 1.0, 0.0, 0.0);
-  clPolicy->update(state.key(mLongitudinalSoundSpeed), state, derivs, 1.0, 0.0, 0.0);
-}
-
 //------------------------------------------------------------------------------
 // initialize (before evaluateDerivatives)
 // This is where we update the damage coupling, which is stored in the
@@ -238,30 +206,13 @@ finalize(const Scalar /*time*/,
     const auto ntot = std::max(1, dataBase.globalNumInternalNodes());
     const auto dfrac = double(nD)/double(ntot);
     mComputeIntersectConnectivity = (dfrac > 0.2);  // Should tune this number...
-    if (Process::getRank() == 0) std::cout << "DamageModel dfrac = " << nD << "/" << ntot << " = " << dfrac << " : " << mComputeIntersectConnectivity << std::endl;
+    // if (Process::getRank() == 0) std::cout << "DamageModel dfrac = " << nD << "/" << ntot << " = " << dfrac << " : " << mComputeIntersectConnectivity << std::endl;
   }
 
   TIME_DamageModel_finalize.stop();
   TIME_Damage.stop();
 }
 
-//------------------------------------------------------------------------------
-// Cull the flaws on each node to the single weakest one.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-void
-DamageModel<Dimension>::
-cullToWeakestFlaws() {
-#pragma omp parallel for
-  for (auto i = 0u; i < mNodeList.numInternalNodes(); ++i) {
-    auto& flaws = mFlaws[i];
-    if (flaws.size() > 0) {
-      const auto maxVal = *max_element(flaws.begin(), flaws.end());
-      flaws = vector<double>(maxVal);
-    }
-  }
-}
-
 //------------------------------------------------------------------------------
 // Access the set of nodes to be excluded from damage.
 //------------------------------------------------------------------------------
@@ -289,39 +240,6 @@ excludeNodes(vector<int> ids) {
   }
 }
 
-//------------------------------------------------------------------------------
-// Compute a Field with the sum of the activation energies per node.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-Field<Dimension, typename Dimension::Scalar>
-DamageModel<Dimension>::
-sumActivationEnergiesPerNode() const {
-  Field<Dimension, Scalar> result("Sum activation energies", mNodeList);
-  for (auto i = 0u; i != mNodeList.numInternalNodes(); ++i) {
-    const vector<double>& flaws = mFlaws(i);
-    for (vector<double>::const_iterator flawItr = flaws.begin();
-         flawItr != flaws.end();
-         ++flawItr) {
-      result(i) += *flawItr;
-    }
-  }
-  return result;
-}
-
-//------------------------------------------------------------------------------
-// Compute a Field with the number of flaws per node.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-Field<Dimension, typename Dimension::Scalar>
-DamageModel<Dimension>::
-numFlawsPerNode() const {
-  Field<Dimension, Scalar> result("num flaws", mNodeList);
-  for (auto i = 0u; i != mNodeList.numInternalNodes(); ++i) {
-    result(i) = flawsForNode(i).size();
-  }
-  return result;
-}
-
 //------------------------------------------------------------------------------
 // Dump the current state to the given file.
 //------------------------------------------------------------------------------
@@ -330,7 +248,6 @@ void
 DamageModel<Dimension>::
 dumpState(FileIO& file, const string& pathName) const {
   file.write(mCrackGrowthMultiplier, pathName + "/crackGrowthMultiplier");
-  file.write(mFlaws, pathName + "/flaws");
   file.write(mExcludeNode, pathName + "/excludeNode");
   file.write(mComputeIntersectConnectivity, pathName + "/computeIntersectConnectivity");
 }
@@ -343,7 +260,6 @@ void
 DamageModel<Dimension>::
 restoreState(const FileIO& file, const string& pathName) {
   file.read(mCrackGrowthMultiplier, pathName + "/crackGrowthMultiplier");
-  file.read(mFlaws, pathName + "/flaws");
   file.read(mExcludeNode, pathName + "/excludeNode");
   file.read(mComputeIntersectConnectivity, pathName + "/computeIntersectConnectivity");
 }

src/Damage/DamageModel.hh

@@ -64,8 +64,7 @@ public:
   DamageModel(SolidNodeList<Dimension>& nodeList,
               const TableKernel<Dimension>& W,
               const double crackGrowthMultiplier,
-              const DamageCouplingAlgorithm damageCouplingAlgorithm,
-              const FlawStorageType& flaws);
+              const DamageCouplingAlgorithm damageCouplingAlgorithm);
   virtual ~DamageModel();
 
   // Compute the generic Grady-Kipp (ala Benz-Asphaug) scalar damage time 
@@ -79,15 +78,11 @@ public:
 
   //...........................................................................
   // Provide a subset of the required physics package interface.
-  virtual void registerState(DataBase<Dimension>& dataBase,
-                             State<Dimension>& state) override;
-
-  virtual 
-  void initialize(const Scalar time,
-                  const Scalar dt,
-                  const DataBase<Dimension>& dataBase,
-                  State<Dimension>& state,
-                  StateDerivatives<Dimension>& derivatives) override;
+  virtual void initialize(const Scalar time,
+                          const Scalar dt,
+                          const DataBase<Dimension>& dataBase,
+                          State<Dimension>& state,
+                          StateDerivatives<Dimension>& derivatives) override;
 
   virtual void finalize(const Scalar time, 
                         const Scalar dt,
@@ -96,17 +91,9 @@ public:
                         StateDerivatives<Dimension>& derivs);
 
   virtual bool requireGhostConnectivity() const override;
-
   virtual bool requireIntersectionConnectivity() const override;
 
   //...........................................................................
-  // Optional method to cull the set of flaws to the single weakest one on
-  // each point.
-  void cullToWeakestFlaws();
-
-  // Get the set of flaw activation energies for the given node index.
-  const std::vector<double> flawsForNode(const size_t index) const;
-
   // Access the SolidNodeList we're damaging.
   SolidNodeList<Dimension>& nodeList();
   const SolidNodeList<Dimension>& nodeList() const;
@@ -123,44 +110,19 @@ public:
   std::vector<int> excludeNodes() const;
   void excludeNodes(std::vector<int> ids);
 
-  // Provide access to the state fields we maintain.
-  const Field<Dimension, Scalar>& youngsModulus() const;
-  const Field<Dimension, Scalar>& longitudinalSoundSpeed() const;
-
-  // Access the flaw field.
-  const FlawStorageType& flaws() const;
-  FlawStorageType& flaws();
-
-  // Compute a Field with the sum of the activation energies per node.
-  Field<Dimension, Scalar> sumActivationEnergiesPerNode() const;
-
-  // Compute a Field with the number of flaws per node.
-  Field<Dimension, Scalar> numFlawsPerNode() const;
-
-  // The effective critical number of nodes per smoothing scale, below which we
-  // assume all flaws are active on a node.
-  double criticalNodesPerSmoothingScale() const;
-  void criticalNodesPerSmoothingScale(double x);
-
   //**************************************************************************
   // Restart methods.
   virtual std::string label() const override { return "DamageModel"; }
   virtual void dumpState(FileIO& file, const std::string& pathName) const;
   virtual void restoreState(const FileIO& file, const std::string& pathName);
   //**************************************************************************
 
-protected:
-  //-------------------------- Protected Interface --------------------------//
-  FlawStorageType mFlaws;
-
 private:
   //--------------------------- Private Interface ---------------------------//
   SolidNodeList<Dimension>& mNodeList;
   const TableKernel<Dimension>& mW;
-  double mCrackGrowthMultiplier, mCriticalNodesPerSmoothingScale;
+  double mCrackGrowthMultiplier;
   DamageCouplingAlgorithm mDamageCouplingAlgorithm;
-  Field<Dimension, Scalar> mYoungsModulus;
-  Field<Dimension, Scalar> mLongitudinalSoundSpeed;
   Field<Dimension, int> mExcludeNode;
   NodeCouplingPtr mNodeCouplingPtr;
   bool mComputeIntersectConnectivity;

src/Damage/DamageModelInline.hh

@@ -21,40 +21,6 @@ requireIntersectionConnectivity() const {
   return mComputeIntersectConnectivity;
 }
 
-//------------------------------------------------------------------------------
-// The effective critical number of nodes per smoothing scale, below which we
-// assume all flaws are active on a node.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-inline
-double
-DamageModel<Dimension>::
-criticalNodesPerSmoothingScale() const {
-  return mCriticalNodesPerSmoothingScale;
-}
-
-template<typename Dimension>
-inline
-void
-DamageModel<Dimension>::
-criticalNodesPerSmoothingScale(double x) {
-  VERIFY(x >= 0.0);
-  mCriticalNodesPerSmoothingScale = x;
-}
-
-//------------------------------------------------------------------------------
-// Return the set of flaw activation energies for the given node.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-inline
-const std::vector<double>
-DamageModel<Dimension>::
-flawsForNode(const size_t index) const {
-  REQUIRE(index < mNodeList.numInternalNodes());
-  REQUIRE(mFlaws.nodeListPtr() == &mNodeList);
-  return mFlaws(index);
-}
-
 //------------------------------------------------------------------------------
 // Access the internal parameters of the model.
 //------------------------------------------------------------------------------
@@ -106,42 +72,4 @@ nodeCoupling() const {
   return *mNodeCouplingPtr;
 }
 
-//------------------------------------------------------------------------------
-// Access the state fields.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-inline
-const Field<Dimension, typename Dimension::Scalar>&
-DamageModel<Dimension>::
-youngsModulus() const {
-  return mYoungsModulus;
-}
-
-template<typename Dimension>
-inline
-const Field<Dimension, typename Dimension::Scalar>&
-DamageModel<Dimension>::
-longitudinalSoundSpeed() const {
-  return mLongitudinalSoundSpeed;
-}
-
-//------------------------------------------------------------------------------
-// Access the flaw activation strains.
-//------------------------------------------------------------------------------
-template<typename Dimension>
-inline
-const typename DamageModel<Dimension>::FlawStorageType&
-DamageModel<Dimension>::
-flaws() const {
-  return mFlaws;
-}
-
-template<typename Dimension>
-inline
-typename DamageModel<Dimension>::FlawStorageType&
-DamageModel<Dimension>::
-flaws() {
-  return mFlaws;
-}
-
 }