added: parallel adaptive (multipatch) support

Author: akva2

Apps/HDF5toVTx/HDF5toVTx.C

@@ -80,6 +80,7 @@ bool readBasis (std::vector<ASMbase*>& result, const std::string& name,
     std::string out;
     hdf.readString(geom.str(),out);
     ptype = out.substr(0,10) == "# LRSPLINE" ? ASM::LRSpline : ASM::Spline;
+    isLR = ptype == ASM::LRSpline;
     int gdim = 0;
     if (ptype == ASM::LRSpline)
       gdim = out.substr(11,7) == "SURFACE" ? 2 : 3;

src/ASM/ASMunstruct.h

@@ -47,6 +47,7 @@ namespace LR //! Utilities for LR-splines.
     IntVec    elements; //!< 0-based indices of the elements to refine
     RealArray errors;   //!< List of error indicators for the elements
     std::vector<IntVec> MLGN; //!< MLGN mapping to use for multipatch
+    IntVec pMLGN; //!< Parallel MLGN mapping to use with MPI
 
     //! \brief Default constructor.
     explicit RefineData(bool rs = false) : refShare(rs) {}

src/ASM/GlobalNodes.C

@@ -13,6 +13,8 @@
 
 #include "GlobalNodes.h"
 #include "ASMunstruct.h"
+#include "SIMbase.h"
+#include "DomainDecomposition.h"
 #include "Utilities.h"
 
 
@@ -178,3 +180,80 @@ GlobalNodes::calcGlobalNodes(const GlobalNodes::LRSplineVec& pchs,
 
   return result;
 }
+
+
+GlobalNodes::IntVec
+GlobalNodes::calcDDMapping(const GlobalNodes::LRSplineVec& pchs,
+                           const std::vector<GlobalNodes::IntVec>& MLGN,
+                           const SIMbase& sim, int& nNodes)
+{
+  const ProcessAdm& adm = sim.getProcessAdm();
+
+  int minNode = 0;
+  if (adm.getProcId() > 0)
+    adm.receive(minNode, adm.getProcId()-1);
+
+  IntVec result(*std::max_element(MLGN.back().begin(), MLGN.back().end()) + 1);
+  std::iota(result.begin(), result.end(), minNode);
+  int maxNode = adm.getProcId() == 0 ? 0 : minNode;
+
+  std::map<int,int> old2new;
+  for (const auto& it : adm.dd.ghostConnections) {
+    int sidx = sim.getLocalPatchIndex(it.slave);
+    if (sidx < 1)
+      continue;
+
+    IntVec lNodes = GlobalNodes::getBoundaryNodes(*pchs[sidx-1], it.dim, it.sidx, 0);
+    for (int& i : lNodes)
+      i = result[MLGN[sidx-1][i]];
+
+    int nRecv;
+    adm.receive(nRecv, adm.dd.getPatchOwner(it.master));
+    if (nRecv =! lNodes.size()) {
+      std::cerr <<"\n *** GlobalNodes::calcDDMapping(): "
+                <<" Topology error, boundary size "
+                << nRecv << ", expected " << lNodes.size() << std::endl;
+      return IntVec();
+    }
+    IntVec glbNodes(lNodes.size());
+    adm.receive(glbNodes, adm.dd.getPatchOwner(it.master));
+    for (size_t i = 0; i < lNodes.size(); ++i)
+      old2new[lNodes[i]] = glbNodes[i];
+  }
+
+  // remap ghost nodes
+  for (auto& it : result)
+    utl::renumber(it, old2new, false);
+
+  // remap rest of our nodes
+  for (int i = 0; i < (int)result.size(); ++i)
+    if (old2new.find(i + minNode) == old2new.end()) {
+      std::map<int,int> old2new2;
+      old2new2[i + minNode] = maxNode++;
+      for (auto& it : result)
+        utl::renumber(it, old2new2, false);
+    }
+
+  if (adm.getProcId() < adm.getNoProcs()-1)
+    adm.send(maxNode, adm.getProcId()+1);
+
+  for (const auto& it : adm.dd.ghostConnections) {
+    int midx = sim.getLocalPatchIndex(it.master);
+    if (midx < 1)
+      continue;
+
+    IntVec glbNodes = GlobalNodes::getBoundaryNodes(*pchs[midx-1], it.dim,
+                                                    it.midx, it.orient);
+    for (size_t i = 0; i < glbNodes.size(); ++i)
+      glbNodes[i] = result[MLGN[midx-1][glbNodes[i]]];
+
+    adm.send(int(glbNodes.size()), adm.dd.getPatchOwner(it.slave));
+    adm.send(glbNodes, adm.dd.getPatchOwner(it.slave));
+  }
+
+#ifdef HAVE_MPI
+  nNodes = adm.allReduce(adm.getProcId() == adm.getNoProcs()-1 ? maxNode : 0, MPI_SUM);
+#endif
+
+  return result;
+}

src/ASM/GlobalNodes.h

@@ -19,6 +19,10 @@
 
 #include <vector>
 
+class DomainDecomposition;
+class ProcessAdm;
+class SIMbase;
+
 
 /*!
   \brief Class establishing global node numbers for unstructed FE models.
@@ -39,12 +43,22 @@ class GlobalNodes
   static IntVec getBoundaryNodes(const LR::LRSpline& lr,
                                  int dim, int lidx, int orient);
 
-
   //! \brief Calculate global node numbers for a FE model.
-  //! \param pchs The spline patches in the model
+   //! \param pchs The spline patches in the model
   //! \param interfaces The topological connections for the spline patches
   static std::vector<IntVec> calcGlobalNodes(const LRSplineVec& pchs,
                                              const InterfaceVec& interfaces);
+
+  //! \brief Calculate parallel global node numbers for a FE model.
+  //! \param pchs The spline patches in the model
+  //! \param MLGN Process-global node numbers
+  //! \param sim The simulator holding patch owner information
+  //! \param adm The parallel process administrator
+  //! \param dd The domain decomposition holding ghost connections
+  //! \param[out] nNodes Total number of nodes in the model
+  static IntVec calcDDMapping(const LRSplineVec& pchs,
+                              const std::vector<GlobalNodes::IntVec>& MLGN,
+                              const SIMbase& sim, int& nNodes);
 };
 
 #endif

src/ASM/LR/ASMu2D.C

@@ -1994,10 +1994,17 @@ void ASMu2D::getBoundaryNodes (int lIndex, IntVec& nodes, int basis,
   default: return;
   }
 
-  if (nodes.empty())
+  if (nodes.empty()) {
     nodes = this->getEdgeNodes(edge, basis, 0);
-  else {
+    if (!local)
+      for (int& node : nodes)
+        node = MLGN[node-1];
+  } else {
     IntVec nodes2 = this->getEdgeNodes(edge, basis, 0);
+    if (!local)
+      for (int& node : nodes2)
+        node = MLGN[node-1];
+
     nodes.insert(nodes.end(), nodes2.begin(), nodes2.end());
   }
 

src/SIM/AdaptiveSIM.C

@@ -311,9 +311,9 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec)
   if (eNorm.cols() < 1 || gNorm[adaptor](adNorm) < errTol*refNorm) return false;
 
   // Calculate row index in eNorm of the error norm to adapt based on
-  size_t i, eRow = adNorm;
+  size_t eRow = adNorm;
   NormBase* norm = model.getNormIntegrand();
-  for (i = 0; i < adaptor; i++)
+  for (size_t i = 0; i < adaptor; i++)
     eRow += norm->getNoFields(i+1);
   delete norm;
 
@@ -363,10 +363,12 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec)
 #ifdef HAS_LRSPLINE
       prm.MLGN = GlobalNodes::calcGlobalNodes(refBasis, model.getInterfaces());
       nNodes = *std::max_element(prm.MLGN.back().begin(), prm.MLGN.back().end()) + 1;
+      if (model.getProcessAdm().getNoProcs() > 1)
+        prm.pMLGN = GlobalNodes::calcDDMapping(refBasis, prm.MLGN, model, nNodes);
 #endif
     }
     errors.reserve(nNodes);
-    for (i = 0; i < nNodes; i++)
+    for (int i = 0; i < nNodes; i++)
       errors.push_back(DblIdx(0.0,i));
 
     for (i = 0; i < errors.size(); i++)
@@ -378,18 +380,21 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec)
 
       // extract element norms for this patch
       Vector locNorm(patch->getNoElms());
-      for (i = 1; i <= patch->getNoElms(); ++i)
+      for (size_t i = 1; i <= patch->getNoElms(); ++i)
         locNorm(i) = eNorm(eRow, patch->getElmID(i));
 
       // remap from geometry basis to refinement basis
       Vector locErr(patch->getNoRefineNodes());
       static_cast<ASMunstruct*>(patch)->remapErrors(locErr, locNorm);
 
       // insert into global error array
-      for (i = 0; i < locErr.size(); ++i)
-        if (model.getNoPatches() > 1)
-          errors[prm.MLGN[patch->idx][i]].first += locErr[i];
-        else
+      for (size_t i = 0; i < locErr.size(); ++i)
+        if (model.getNoPatches() > 1) {
+          if (prm.pMLGN.empty())
+            errors[prm.MLGN[patch->idx][i]].first += locErr[i];
+          else
+            errors[prm.pMLGN[prm.MLGN[patch->idx][i]]].first += locErr[i];
+        } else
           errors[i].first += locErr[i];
     }
   }
@@ -400,10 +405,21 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec)
                 <<" is available for isotropic_function only."<< std::endl;
       return false;
     }
-    for (i = 0; i < eNorm.cols(); i++)
+    for (size_t i = 0; i < eNorm.cols(); i++)
       errors.push_back(DblIdx(eNorm(eRow,1+i),i));
   }
 
+#ifdef HAVE_MPI
+  std::vector<double> perr(errors.size(), 0.0);
+  if (model.getProcessAdm().getNoProcs() > 1) {
+    for (size_t i = 0; i < errors.size(); ++i)
+      perr[i] = errors[i].first;
+    model.getProcessAdm().allReduce(perr, MPI_SUM);
+    for (size_t i = 0; i < errors.size(); ++i)
+      errors[i].first = perr[i];
+  }
+#endif
+
   // Sort the elements in the sequence of decreasing errors
   std::sort(errors.begin(),errors.end(),std::greater<DblIdx>());
 
@@ -468,8 +484,15 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec)
     return false;
 
   prm.elements.reserve(refineSize);
-  for (i = 0; i < refineSize; i++)
-    prm.elements.push_back(errors[i].second);
+  for (size_t i = 0; i < refineSize; i++)
+    if (prm.pMLGN.empty())
+      prm.elements.push_back(errors[i].second);
+    else {
+      auto it = std::find(prm.pMLGN.begin(),
+                          prm.pMLGN.end(), errors[i].second);
+      if (it != prm.pMLGN.end())
+        prm.elements.push_back(it-prm.pMLGN.begin());
+    }
 
   // Now refine the mesh
   if (!storeMesh)

src/SIM/SIMbase.C

@@ -99,7 +99,6 @@ void SIMbase::clearProperties ()
   for (auto& i4 : myTracs)
     delete i4.second;
 
-  myPatches.clear();
   myGlb2Loc.clear();
   myScalars.clear();
   myVectors.clear();

src/SIM/SIMinput.C

@@ -1124,7 +1124,7 @@ bool SIMinput::refine (const LR::RefineData& prm,
     }
 
   // Single patch models only pass refinement call to the ASM level
-  if (myModel.size() == 1)
+  if (this->getNoPatches() == 1)
     return (isRefined = pch->refine(prm,sol,fName));
 
   if (!prm.errors.empty()) // refinement by true_beta
@@ -1196,6 +1196,49 @@ bool SIMinput::refine (const LR::RefineData& prm,
     }
   }
 
+#ifdef HAVE_MPI
+  if (this->adm.getNoProcs() > 1) {
+    for (int i = 0; i < this->adm.getNoProcs(); ++i) {
+      int send = 0;
+      if (i == this->adm.getProcId())
+        for (auto& it : conformingIndices)
+          send += it.size();
+
+      int nRecv = adm.allReduce(send, MPI_SUM);
+      std::vector<int> vRecv;
+      if (i == this->adm.getProcId()) {
+        std::vector<int> vSend;
+        vRecv.reserve(nRecv);
+        for (auto& it : conformingIndices)
+          for (const int& node : it)
+            vRecv.push_back(prm.pMLGN[node]);
+      } else
+        vRecv.resize(nRecv);
+
+      MPI_Bcast(vRecv.data(), vRecv.size(), MPI_INT, i, *this->adm.getCommunicator());
+
+      if (i != this->adm.getProcId()) {
+        for (int& node : vRecv) {
+          auto it = std::find(prm.pMLGN.begin(), prm.pMLGN.end(), node);
+          if (it != prm.pMLGN.end()) {
+            int globId = it - prm.pMLGN.begin();
+            int locId;
+            for (size_t j = 0; j < myModel.size(); j++)
+              if (prm.MLGN.empty()) {
+                if ((locId = myModel[j]->getNodeIndex(globId)) > 0)
+                  conformingIndices[j].insert(locId-1);
+              } else {
+                auto it = std::find(prm.MLGN[j].begin(), prm.MLGN[j].end(), globId);
+                if (it != prm.MLGN[j].end())
+                  conformingIndices[j].insert(it-prm.MLGN[j].begin());
+              }
+          }
+        }
+      }
+    }
+  }
+#endif
+
   Vectors lsols;
   lsols.reserve(sol.size()*myModel.size());
   for (size_t i = 0; i < myModel.size(); i++)