Feature/cull pairs

src/FSISPH/SolidFSISPHEvaluateDerivatives.cc

@@ -233,6 +233,10 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
       j = pairs[kk].j_node;
       nodeListi = pairs[kk].i_list;
       nodeListj = pairs[kk].j_list;
+      if (pairs[kk].f_couple < 1.0e-16)
+      {
+         continue;
+      }
 
       // Get the state for node i.
       const auto& interfaceFlagsi = interfaceFlags(nodeListi,i);
@@ -556,8 +560,8 @@ secondDerivativesLoop(const typename Dimension::Scalar time,
         if (freeParticle) {
           DvDti += mj*deltaDvDt;
           DvDtj -= mi*deltaDvDt;
-        } 
-      
+        }
+
         // Velocity Gradient
         //-----------------------------------------------------------
         // construct our interface velocity 

src/FSISPH/SolidFSISPHHydroBase.cc

@@ -520,6 +520,8 @@ preStepInitialize(const DataBase<Dimension>& dataBase,
   TIME_BEGIN("SolidFSISPHpreStepInitialize");
   if (mApplySelectDensitySum){
   switch(this->densityUpdate()){
+   case FSIMassDensityMethod::FSIConsistentSumMassDensity:
+     // fallthrough intended
    case FSIMassDensityMethod::FSISumMassDensity:
     {
       const auto& W = this->kernel();
@@ -528,7 +530,8 @@ preStepInitialize(const DataBase<Dimension>& dataBase,
       const auto  mass = state.fields(HydroFieldNames::mass, 0.0);
       const auto  H = state.fields(HydroFieldNames::H, SymTensor::zero);
             auto  massDensity = state.fields(HydroFieldNames::massDensity, 0.0);
-      computeFSISPHSumMassDensity(connectivityMap, W, mSumDensityNodeLists, position, mass, H, massDensity);
+      const auto consistentSum = (this->densityUpdate() == FSIMassDensityMethod::FSIConsistentSumMassDensity);
+      computeFSISPHSumMassDensity(connectivityMap, W, mSumDensityNodeLists, position, mass, H, consistentSum, massDensity);
       for (auto boundaryItr = this->boundaryBegin(); boundaryItr < this->boundaryEnd(); ++boundaryItr) (*boundaryItr)->applyFieldListGhostBoundary(massDensity);
       for (auto boundaryItr = this->boundaryBegin(); boundaryItr < this->boundaryEnd(); ++boundaryItr) (*boundaryItr)->finalizeGhostBoundary();
      break;

src/FSISPH/SolidFSISPHHydroBase.hh

@@ -32,6 +32,7 @@ enum class FSIMassDensityMethod {
   FSISumMassDensity = 0,
   PressureCorrectSumMassDensity = 1,
   HWeightedSumMassDensity = 2,
+  FSIConsistentSumMassDensity = 3,
 };
 
 template<typename Dimension> class State;

src/FSISPH/computeFSISPHSumMassDensity.cc

@@ -14,6 +14,7 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
                             const FieldList<Dimension, typename Dimension::Vector>& position,
                             const FieldList<Dimension, typename Dimension::Scalar>& mass,
                             const FieldList<Dimension, typename Dimension::SymTensor>& H,
+                            const bool consistentSum,
                             FieldList<Dimension, typename Dimension::Scalar>& massDensity) {
 
   // Pre-conditions.
@@ -61,14 +62,15 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
       const auto& ri = position(nodeListi, i);
       const auto& rj = position(nodeListj, j);
       const auto  rij = ri - rj;
+      const auto normalSum = (nodeListi == nodeListj) || consistentSum;
 
       if(sumDensityNodeLists[nodeListi]==1){
         const auto& Hi = H(nodeListi, i);
         const auto  Hdeti = Hi.Determinant();
         const auto  etai = (Hi*rij).magnitude();
         const auto  Wi = W.kernelValue(etai, Hdeti);
 
-        massDensity_thread(nodeListi, i) += (nodeListi == nodeListj ? mj : mi)*Wi; 
+        massDensity_thread(nodeListi, i) += (normalSum ? mj : mi)*Wi; 
       } 
 
       if(sumDensityNodeLists[nodeListj]==1){
@@ -77,7 +79,7 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
         const auto  etaj = (Hj*rij).magnitude();
         const auto  Wj = W.kernelValue(etaj, Hdetj);
 
-        massDensity_thread(nodeListj, j) += (nodeListi == nodeListj ? mi : mj)*Wj;
+        massDensity_thread(nodeListj, j) += (normalSum ? mi : mj)*Wj;
       }
     }
 

src/FSISPH/computeFSISPHSumMassDensity.hh

@@ -22,9 +22,10 @@ computeFSISPHSumMassDensity(const ConnectivityMap<Dimension>& connectivityMap,
                             const FieldList<Dimension, typename Dimension::Vector>& position,
                             const FieldList<Dimension, typename Dimension::Scalar>& mass,
                             const FieldList<Dimension, typename Dimension::SymTensor>& H,
+                            const bool consistentSum,
                             FieldList<Dimension, typename Dimension::Scalar>& massDensity);
 
 }
 
 
- #endif
+ #endif

src/FSISPH/computeFSISPHSumMassDensityInst.cc.py

@@ -14,6 +14,7 @@
                                             const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>&,
                                             const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>&,
                                             const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>&,
+                                            const bool consistentSum,
                                                   FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>&);
 }
 """

src/Neighbor/ConnectivityMap.cc

@@ -134,7 +134,8 @@ ConnectivityMap():
   mNodeTraversalIndices(),
   mKeys(FieldStorageType::CopyFields),
   mCouplingPtr(std::make_shared<NodeCoupling>()),
-  mIntersectionConnectivity() {
+  mIntersectionConnectivity(),
+  mExcludePairs([](int, int, int, int) { return false; }){
 }
 
 //------------------------------------------------------------------------------
@@ -269,10 +270,9 @@ patchConnectivity(const FieldList<Dimension, int>& flags,
       const auto jNodeList = mNodePairList[k].j_list;
       const auto i = mNodePairList[k].i_node;
       const auto j = mNodePairList[k].j_node;
-      if (flags(iNodeList, i) != 0 and flags(jNodeList, j) != 0) {
+      if (flags(iNodeList, i) != 0 and flags(jNodeList, j) != 0)
         culledPairs_thread.push_back(NodePairIdxType(old2new(iNodeList, i), iNodeList,
                                                      old2new(jNodeList, j), jNodeList));
-      }
     }
 #pragma omp critical
     {
@@ -919,9 +919,14 @@ computeConnectivity() {
 
                     // We don't include self-interactions.
                     if ((iNodeList != jNodeList) or (i != j)) {
-                      neighbors[jNodeList].push_back(j);
-                      if (calculatePairInteraction(iNodeList, i, jNodeList, j, firstGhostNodej)) nodePairs_private.push_back(NodePairIdxType(i, iNodeList, j, jNodeList));
-                      if (domainDecompIndependent) keys[jNodeList].push_back(pair<int, Key>(j, mKeys(jNodeList, j)));
+
+                      // Exclude specified pairs.
+                      if (!mExcludePairs(iNodeList, i, jNodeList, j))
+                      {
+                        neighbors[jNodeList].push_back(j);
+                        if (calculatePairInteraction(iNodeList, i, jNodeList, j, firstGhostNodej)) nodePairs_private.push_back(NodePairIdxType(i, iNodeList, j, jNodeList));
+                        if (domainDecompIndependent) keys[jNodeList].push_back(pair<int, Key>(j, mKeys(jNodeList, j)));
+                      }
                     }
                   }
                 }
@@ -1126,6 +1131,8 @@ computeConnectivity() {
   //   }
   // }
 
+  mNodePairList.makeUnique();
+
   // Post conditions.
   BEGIN_CONTRACT_SCOPE
   // Make sure that the correct number of nodes have been completed.
@@ -1145,4 +1152,52 @@ computeConnectivity() {
   TIME_END("ConnectivityMap_computeConnectivity");
 }
 
+//------------------------------------------------------------------------------
+// Remove given pairs from the connectivity in place
+//------------------------------------------------------------------------------
+template<typename Dimension>
+void
+ConnectivityMap<Dimension>::
+removeConnectivity(std::function<bool(int, int, int, int)> excludePairs)
+{
+  const auto domainDecompIndependent = NodeListRegistrar<Dimension>::instance().domainDecompositionIndependent();
+  const auto numNodeLists = mNodeLists.size();
+  for (auto iNodeList = 0u; iNodeList != numNodeLists; ++iNodeList) {
+    const auto numNodes = ((domainDecompIndependent or mBuildGhostConnectivity or mBuildOverlapConnectivity) ? 
+                           mNodeLists[iNodeList]->numNodes() :
+                           mNodeLists[iNodeList]->numInternalNodes());
+    for (auto i = 0u; i < numNodes; ++i) {
+      auto& neighbors = mConnectivity[i];
+      CHECK(neighbors.size() == numNodeLists);
+      for (auto jNodeList = 0u; jNodeList < numNodeLists; ++jNodeList) {
+        auto& neighborsj = neighbors[jNodeList];
+        auto l = 0u;
+        for (auto k = 0u; k < neighborsj.size(); ++k)
+        {
+          const auto j = neighborsj[k];
+          if (!excludePairs(iNodeList, i, jNodeList, j))
+          {
+            neighborsj[l++] = neighborsj[k];
+          }
+        }
+        neighborsj.resize(l);
+      }
+    }
+  }
+
+  const auto npairs = mNodePairList.size();
+  auto l = 0u;
+  for (auto k = 0u; k < npairs; ++k) {
+    const auto iNodeList = mNodePairList[k].i_list;
+    const auto jNodeList = mNodePairList[k].j_list;
+    const auto i = mNodePairList[k].i_node;
+    const auto j = mNodePairList[k].j_node;
+    if (!excludePairs(iNodeList, i, jNodeList, j))
+    {
+      mNodePairList[l++] = mNodePairList[k];
+    }
+  }
+  mNodePairList.resize(l);
+}
+
 }

src/Neighbor/ConnectivityMap.hh

@@ -63,6 +63,9 @@ public:
   // member of a pair (maintaining symmetry).
   void removeConnectivity(const FieldList<Dimension, std::vector<std::vector<int>>>& neighborsToCut);
 
+  // Remove pairs based on a function
+  void removeConnectivity(std::function<bool(int, int, int, int)> excludePairs);
+
   // Are we computing neighbors for ghosts?
   bool buildGhostConnectivity() const;
 
@@ -171,9 +174,11 @@ public:
   // Return which NodeList index in order the given one would be in our connectivity.
   unsigned nodeListIndex(const NodeList<Dimension>* nodeListPtr) const;
 
+  // Exclude (nodeListi, nodei, nodeListj, nodej) pairs if this returns false
+  void setNodePairExclusion(std::function<bool(int, int, int, int)> excludePairs);
+
   // Check that the internal data structure is valid.
   bool valid() const;
-
 private:
   //--------------------------- Private Interface ---------------------------//
   // The set of NodeLists.
@@ -208,6 +213,9 @@ private:
   using IntersectionConnectivityContainer = std::unordered_map<NodePairIdxType, std::vector<std::vector<int>>>;
   IntersectionConnectivityContainer mIntersectionConnectivity;
 
+  // Exclude pairs matching this function
+  std::function<bool(int, int, int, int)> mExcludePairs;
+
   // Internal method to fill in the connectivity, once the set of NodeLists 
   // is determined.
   void computeConnectivity();

src/Neighbor/ConnectivityMapInline.hh

@@ -390,4 +390,15 @@ intersectionConnectivity(const NodePairIdxType& pair) const {
   return itr->second;
 }
 
+//------------------------------------------------------------------------------
+// Set the function for excluding node pairs
+//------------------------------------------------------------------------------
+template<typename Dimension>
+void
+ConnectivityMap<Dimension>::
+setNodePairExclusion(std::function<bool(int, int, int, int)> excludePairs)
+{
+  mExcludePairs = excludePairs;
+}
+
 }

src/Neighbor/NodePairList.cc

@@ -1,5 +1,7 @@
 #include "NodePairList.hh"
 
+#include <algorithm>
+
 namespace Spheral {
 
   NodePairIdxType::NodePairIdxType(int i_n, int i_l, int j_n, int j_l, double f) :
@@ -15,4 +17,31 @@ namespace Spheral {
     mNodePairList.clear();
   }
 
+  void NodePairList::makeUnique() {
+    // Make sure the node pairs are ordered correctly
+    for (auto kk = 0u; kk < mNodePairList.size(); ++kk)
+    {
+      auto& pair = mNodePairList[kk];
+      if (pair.i_list > pair.j_list)
+      {
+        const auto temp_list = pair.j_list;
+        const auto temp_node = pair.j_node;
+        pair.j_list = pair.i_list;
+        pair.j_node = pair.i_node;
+        pair.i_list = temp_list;
+        pair.i_node = temp_node;
+      }
+      if (pair.i_list == pair.j_list && pair.i_node > pair.j_node)
+      {
+        const auto temp = pair.j_node;
+        pair.j_node = pair.i_node;
+        pair.i_node = temp;
+      }
+    }
+
+    // Remove duplicates
+    std::sort(mNodePairList.begin(), mNodePairList.end());
+    mNodePairList.erase(std::unique(mNodePairList.begin(), mNodePairList.end()), mNodePairList.end());
+  }
+
 }

src/Neighbor/NodePairList.hh

@@ -58,8 +58,10 @@ public:
 
   NodePairList();
   void push_back(NodePairIdxType nodePair);
-  void clear(); 
+  void clear();
+  void resize(const size_t n) { mNodePairList.resize(n, NodePairIdxType(-1,-1,-1,-1)); }
   size_t size() const { return mNodePairList.size(); }
+  void makeUnique();
 
   // Iterators
   iterator begin() { return mNodePairList.begin(); }