Validation of threading groups based on SAM data

src/ASM/ASMbase.h

@@ -45,6 +45,7 @@ class RealFunc;
 class VecFunc;
 class Vec3;
 class Tensor;
+class SAM;
 
 namespace ASM { class InterfaceChecker; }
 namespace utl { template<class Arg, class Result> class Function; }
@@ -515,10 +516,13 @@ class ASMbase
   virtual void generateThreadGroupsFromElms(const IntVec&) {}
   //! \brief Generates element groups for multi-threading based on a partition.
   virtual void generateProjThreadGroupsFromElms(const IntVec&) {}
+  //! \brief Validates the threading groups based on the assembly data in %SAM.
+  virtual bool validateThreadGroups(const SAM*) const { return true; }
 
   //! \brief Hook for changing number of threads.
   virtual void changeNumThreads() {}
 
+
   // Methods for integration of finite element quantities.
   // These are the main computational methods of the ASM class hierarchy.
   // ====================================================================

src/ASM/ASMs2D.C

@@ -36,6 +36,8 @@
 #include "Point.h"
 #include "Tensor.h"
 #include "MPC.h"
+#include "SAM.h"
+#include "DiagMatrix.h"
 #include "IFEM.h"
 #include <array>
 #include <utility>
@@ -3356,6 +3358,48 @@ void ASMs2D::generateProjThreadGroupsFromElms (const IntVec& elms)
 }
 
 
+bool ASMs2D::validateThreadGroups (const SAM* sam) const
+{
+  IFEM::cout <<"\nValidating element groups for multi-threaded assembly."
+             << std::endl;
+
+  if (threadGroups[0].size() == 1)
+    return true; // Only one group (no multi-threading)
+
+  size_t nErr = 0;
+  int iTGroup = 0;
+  for (size_t g = 0; g < threadGroups.size(); g++)
+    for (const IntVec& group : threadGroups[g])
+    {
+      ++iTGroup;
+      size_t lErr = nErr;
+      DiagMatrix sysMat(sam->getNoEquations());
+      for (int iel : group)
+        if (IntVec meen; !sam->getElmEqns(meen,MLGE[iel]))
+          ++nErr;
+        else if (!sysMat.assembleStruct(iTGroup,*sam,meen))
+          ++nErr;
+
+      for (size_t ieq = 1; ieq <= sysMat.dim(0); ieq++)
+        if (int count = sysMat(ieq); count > 1 && count/1000 != iTGroup)
+        {
+          std::cerr <<" *** Threading group "<< iTGroup <<" has "<< count
+                    <<" contributors to equation "<< ieq;
+          std::pair<int,int> dof = sam->getNodeAndLocalDof(ieq,true);
+          std::cerr <<" (node "<< dof.first <<" local dof "<< dof.second <<")"
+                    << std::endl;
+          nErr += count;
+        }
+
+      if (lErr == nErr)
+        IFEM::cout <<"   * Thread group "<< iTGroup <<" (size "<< group.size()
+                   <<") is OK"<< std::endl;
+    }
+
+  return nErr == 0;
+}
+
+
 bool ASMs2D::addRigidCpl (int lindx, int ldim, int basis,
                           int& gMaster, const Vec3& Xmaster, bool extraPt)
 {

src/ASM/ASMs2D.h

@@ -698,6 +698,9 @@ class ASMs2D : public ASMstruct, public ASM2D
   //! \brief Generate element groups from a partition.
   virtual void generateProjThreadGroupsFromElms(const IntVec& elms);
 
+  //! \brief Validates the threading groups based on the assembly data in %SAM.
+  virtual bool validateThreadGroups(const SAM* sam) const;
+
   //! \brief Returns 0-based index of first node on integration basis.
   virtual int getFirstItgElmNode() const { return 0; }
   //! \brief Returns 0-based index of last node on integration basis.

src/ASM/ASMu2DLag.C

@@ -463,6 +463,7 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence,
   IntVec status(nel,0);
 
   using IntSet = std::set<int>;
+  IntSet slaveNodes;
   std::vector<IntSet> nodeConn(nnod); // node-to-element connectivity
   std::vector<IntSet> nodeNode(nnod); // node-to-node connectivity via MPCs
 
@@ -471,28 +472,48 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence,
       for (size_t i = 0; i < mpc->getNoMaster(); i++)
         if (int mastr = this->getNodeIndex(mpc->getMaster(i).node); mastr > 0)
         {
-          nodeNode[slave-1].insert(mastr-1);
+          slaveNodes.insert(slave-1);
           nodeNode[mastr-1].insert(slave-1);
         }
 
   size_t fixedElements = 0;
   for (size_t iel = 0; iel < nel; iel++)
     if (separateGroup1noded && MNPC[iel].size() == 1 &&
-        nodeNode[MNPC[iel].front()].empty())
+        slaveNodes.find(MNPC[iel].front()) == slaveNodes.end())
     {
-      status[iel] = 1; // Separate group for all single-noded elements
+      // Separate color for all single-noded elements
+      status[iel] = 1; // whose element node is not a slave
       if (++fixedElements == 1)
         threadGroups[0].resize(1, { static_cast<int>(iel) });
       else
         threadGroups[0].front().push_back(iel);
     }
     else
       for (int node : MNPC[iel])
-      {
         nodeConn[node].insert(iel);
-        for (int master : nodeNode[node])
-          nodeConn[master].insert(iel);
-      }
+
+  // Second pass - account for the MPC couplings
+  for (size_t master = 0; master < nodeNode.size(); master++)
+    if (!nodeNode[master].empty())
+    {
+      IntSet commonElms(nodeConn[master]);
+      for (int slave : nodeNode[master])
+        commonElms.insert(nodeConn[slave].begin(),nodeConn[slave].end());
+      nodeConn[master].insert(commonElms.begin(),commonElms.end());
+      for (int slave : nodeNode[master])
+        nodeConn[slave] = commonElms;
+    }
+
+#if SP_DEBUG > 1
+  std::cout <<"\nNode to element connectivity (incl. MPC couplings):";
+  for (size_t inod = 0; inod < nnod; inod++)
+  {
+    std::cout <<"\n"<< inod <<" ->";
+    for (int e : nodeConn[inod])
+      std::cout <<" "<< e;
+  }
+  std::cout << std::endl;
+#endif
 
   for (size_t nColors = fixedElements > 0; fixedElements < nel; ++nColors)
   {
@@ -501,7 +522,7 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence,
                   [](int& s) { if (s < 0) s = 0; });
 
     // Look for available elements
-    IntVec thisColor;
+    IntVec& thisColor = threadGroups[0].emplace_back();
     for (size_t i = 0; i < nel; ++i)
       if (status[i] == 0)
       {
@@ -514,8 +535,6 @@ void ASMu2DLag::generateThreadGroupsMultiColored (bool silence,
             if (status[j] == 0) // if not assigned a color yet
               status[j] = -1;   // set as unavailable (with current color)
       }
-
-    threadGroups[0].emplace_back(thisColor);
   }
 
   if (!silence)

src/LinAlg/DiagMatrix.C

@@ -62,7 +62,7 @@ void DiagMatrix::dump (std::ostream& os, LinAlg::StorageFormat format,
 
 bool DiagMatrix::assemble (const Matrix& eM, const SAM& sam, int e)
 {
-  if (myMat.size() != (size_t)sam.neq)
+  if (static_cast<int>(myMat.size()) != sam.neq)
     return false;
 
   StdVector B;
@@ -74,22 +74,67 @@ bool DiagMatrix::assemble (const Matrix& eM, const SAM& sam, int e)
 bool DiagMatrix::assemble (const Matrix& eM, const SAM& sam,
                            SystemVector&, const IntVec& meq)
 {
-  if (meq.size() != 1 || eM.size() != 1)
+  const size_t nedof = meq.size();
+  if (eM.rows() < nedof || eM.cols() < nedof)
   {
-    std::cerr <<" *** DiagMatrix::assemble: Only for one-dof elements, nedof = "
-              << meq.size() <<" size(eM) = " << eM.size() << std::endl;
+    std::cerr <<" *** DiagMatrix::assemble: Invalid element matrix, nedof = "
+              << nedof <<" size(eM) = "<< eM.rows() <<","<< eM.cols()
+              << std::endl;
     return false;
   }
 
-  int ieq = meq.front();
-  if (ieq < 1 || ieq > sam.neq)
-  {
-    std::cerr <<" *** DiagMatrix::assemble: ieq="<< ieq <<" is out or range [1,"
-              << sam.neq <<"]"<< std::endl;
-    return false;
-  }
+  // Add elements corresponding to free dofs in eM,
+  // and (appropriately weighted) elements corresponding to constrained dofs,
+  // into the diagonal system matrix
+  for (size_t j = 1; j <= nedof; j++)
+    if (int jeq = meq[j-1]; jeq > 0)
+      myMat(jeq) += eM(j,j);
+    else if (int jceq = -meq[j-1]; jceq > 0)
+      for (int jp = sam.mpmceq[jceq-1]; jp < sam.mpmceq[jceq]-1; jp++)
+        if (int jeq = sam.mmceq[jp] > 0 ? sam.meqn[sam.mmceq[jp]-1]:0; jeq > 0)
+          for (size_t i = 1; i <= nedof; i++)
+            if (meq[i-1] == jeq)
+              myMat(jeq) += sam.ttcc[jp]*(eM(i,j) + eM(j,i));
+            else if (int iceq = -meq[i-1]; iceq > 0)
+              for (int ip = sam.mpmceq[iceq-1]; ip < sam.mpmceq[iceq]-1; ip++)
+                if (sam.mmceq[ip] > 0 && sam.meqn[sam.mmceq[ip]-1] == jeq)
+                  myMat(jeq) += sam.ttcc[ip]*sam.ttcc[jp]*eM(i,j);
+
+  return this->flagNonZeroEqs(meq);
+}
+
+
+/*!
+  This method is doing the same as the above assemble() method,
+  but where all element matrices are assumed identity matrices.
+  The purpose is to verify that each equation only gets a single contribution
+  within each threading group.
+*/
 
-  myMat(ieq) += eM(1,1);
+bool DiagMatrix::assembleStruct (int val, const SAM& sam, const IntVec& meq)
+{
+  // Lambda function marking matrix contributions due to linear couplings.
+  auto&& addEq = [&diag=myMat,val](int ieq)
+  {
+    int current = diag(ieq)/1000.0;
+    if (current != val)
+      diag(ieq) += 1000.0*val;
+  };
+
+  const size_t nedof = meq.size();
+  for (size_t j = 1; j <= nedof; j++)
+    if (int jeq = meq[j-1]; jeq > 0)
+      myMat(jeq) += 1.0;
+    else if (int jceq = -meq[j-1]; jceq > 0)
+      for (int jp = sam.mpmceq[jceq-1]; jp < sam.mpmceq[jceq]-1; jp++)
+        if (int jeq = sam.mmceq[jp] > 0 ? sam.meqn[sam.mmceq[jp]-1]:0; jeq > 0)
+          for (size_t i = 1; i <= nedof; i++)
+            if (meq[i-1] == jeq)
+              addEq(jeq);
+            else if (int iceq = -meq[i-1]; iceq > 0)
+              for (int ip = sam.mpmceq[iceq-1]; ip < sam.mpmceq[iceq]-1; ip++)
+                if (sam.mmceq[ip] > 0 && sam.meqn[sam.mmceq[ip]-1] == jeq)
+                  addEq(jeq);
 
   return this->flagNonZeroEqs(meq);
 }

src/LinAlg/DiagMatrix.h

@@ -101,6 +101,14 @@ class DiagMatrix : public SystemMatrix
   //! \details To be used when there is no underlying SAM
   virtual bool assemble(const Matrix& eM, const IntVec& meq);
 
+  //! \brief Assembles a validation matrix assuming unit element contributions.
+  //! \param[in] val Value identifying current threading group
+  //! \param[in] sam Auxiliary data for FE assembly management
+  //! \param[in] meq Matrix of element equation numbers
+  //!
+  //! \details For validation of threading groups.
+  bool assembleStruct(int val, const SAM& sam, const IntVec& meq);
+
   //! \brief Multiplication with a scalar.
   virtual void mult(Real alpha) { myMat *= alpha; }