Add partition of unity RBF data mapping

docs/changelog/1483.md

@@ -0,0 +1 @@
+- Added partition of unity RBF mapping as a more efficient alternative to existing RBF mapping for large cases. The mapping can be configured using `<mapping:rbf-pum-direct ... > <basis-function:... /> </mapping:rbf-pum-direct>`.

src/mapping/RadialBasisFctSolver.hpp

@@ -2,6 +2,7 @@
 
 #include <Eigen/Cholesky>
 #include <Eigen/QR>
+#include <Eigen/SVD>
 #include <boost/range/adaptor/indexed.hpp>
 #include <boost/range/irange.hpp>
 #include <numeric>
@@ -27,7 +28,13 @@ class RadialBasisFctSolver {
   /// Default constructor
   RadialBasisFctSolver() = default;
 
-  /// Assembles the system matrices and computes the decomposition of the interpolation matrix
+  /**
+   * assembles the system matrices and computes the decomposition of the interpolation matrix
+   * inputMesh refers to the mesh where the interpolants are built on, i.e., the input mesh
+   * for consistent mappings and the output mesh for conservative mappings
+   * outputMesh refers to the mesh where we evaluate the interpolants, i.e., the output mesh
+   * consistent mappings and the input mesh for conservative mappings
+  */
   template <typename IndexContainer>
   RadialBasisFctSolver(RADIAL_BASIS_FUNCTION_T basisFunction, const mesh::Mesh &inputMesh, const IndexContainer &inputIDs,
                        const mesh::Mesh &outputMesh, const IndexContainer &outputIDs, std::vector<bool> deadAxis, Polynomial polynomial);
@@ -79,6 +86,30 @@ inline double computeSquaredDifference(
   return std::accumulate(v.begin(), v.end(), static_cast<double>(0.), [](auto &res, auto &val) { return res + val * val; });
 }
 
+/// given the active axis, computes sets the axis with the lowest spatial expansion to dead
+template <typename IndexContainer>
+constexpr void reduceActiveAxis(const mesh::Mesh &mesh, const IndexContainer &IDs, std::array<bool, 3> &axis)
+{
+  // make a pair of the axis and the difference
+  std::array<std::pair<int, double>, 3> differences;
+
+  // Compute the difference magnitude per direction
+  for (std::size_t d = 0; d < axis.size(); ++d) {
+    // Ignore dead axis here, i.e., apply the max value such that they are sorted on the last position(s)
+    if (axis[d] == false) {
+      differences[d] = std::make_pair<int, double>(d, std::numeric_limits<double>::max());
+    } else {
+      auto res = std::minmax_element(IDs.begin(), IDs.end(), [&](const auto &a, const auto &b) { return mesh.vertices()[a].rawCoords()[d] < mesh.vertices()[b].rawCoords()[d]; });
+      // Check if we are above or below the threshold
+      differences[d] = std::make_pair<int, double>(d, std::abs(mesh.vertices()[*res.second].rawCoords()[d] - mesh.vertices()[*res.first].rawCoords()[d]));
+    }
+  }
+
+  std::sort(differences.begin(), differences.end(), [](const auto &d1, const auto &d2) { return d1.second < d2.second; });
+  // Disable the axis having the smallest expansion
+  axis[differences[0].first] = false;
+}
+
 // Fill in the polynomial entries
 template <typename IndexContainer>
 inline void fillPolynomialEntries(Eigen::MatrixXd &matrix, const mesh::Mesh &mesh, const IndexContainer &IDs, Eigen::Index startIndex, std::array<bool, 3> activeAxis)
@@ -214,15 +245,34 @@ RadialBasisFctSolver<RADIAL_BASIS_FUNCTION_T>::RadialBasisFctSolver(RADIAL_BASIS
   // In case we deal with separated polynomials, we need dedicated matrices for the polynomial contribution
   if (polynomial == Polynomial::SEPARATE) {
 
-    // 1. Allocate memory for these matrices
     // 4 = 1 + dimensions(3) = maximum number of polynomial parameters
-    const unsigned int polyParams = 4 - std::count(activeAxis.begin(), activeAxis.end(), false);
-    _matrixQ.resize(inputIDs.size(), polyParams);
-    _matrixV.resize(outputIDs.size(), polyParams);
+    auto         localActiveAxis = activeAxis;
+    unsigned int polyParams      = 4 - std::count(localActiveAxis.begin(), localActiveAxis.end(), false);
+
+    do {
+      // First, build matrix Q and check for the condition number
+      _matrixQ.resize(inputIDs.size(), polyParams);
+      fillPolynomialEntries(_matrixQ, inputMesh, inputIDs, 0, localActiveAxis);
+
+      // Compute the condition number
+      Eigen::JacobiSVD<Eigen::MatrixXd> svd(_matrixQ);
+      PRECICE_ASSERT(svd.singularValues().size() > 0);
+      PRECICE_DEBUG("Singular values in polynomial solver: {}", svd.singularValues());
+      const double conditionNumber = svd.singularValues()(0) / std::max(svd.singularValues()(svd.singularValues().size() - 1), math::NUMERICAL_ZERO_DIFFERENCE);
+      PRECICE_DEBUG("Condition number: {}", conditionNumber);
+
+      // Disable one axis
+      if (conditionNumber > 1e5) {
+        reduceActiveAxis(inputMesh, inputIDs, localActiveAxis);
+        polyParams = 4 - std::count(localActiveAxis.begin(), localActiveAxis.end(), false);
+      } else {
+        break;
+      }
+    } while (true);
 
-    // 2. fill the matrices: Q for the inputMesh, V for the outputMesh
-    fillPolynomialEntries(_matrixQ, inputMesh, inputIDs, 0, activeAxis);
-    fillPolynomialEntries(_matrixV, outputMesh, outputIDs, 0, activeAxis);
+    // allocate and fill matrix V for the outputMesh
+    _matrixV.resize(outputIDs.size(), polyParams);
+    fillPolynomialEntries(_matrixV, outputMesh, outputIDs, 0, localActiveAxis);
 
     // 3. compute decomposition
     _qrMatrixQ = _matrixQ.colPivHouseholderQr();

src/mapping/config/MappingConfiguration.cpp

@@ -13,6 +13,7 @@
 #include "mapping/NearestNeighborGradientMapping.hpp"
 #include "mapping/NearestNeighborMapping.hpp"
 #include "mapping/NearestProjectionMapping.hpp"
+#include "mapping/PartitionOfUnityMapping.hpp"
 #include "mapping/PetRadialBasisFctMapping.hpp"
 #include "mapping/RadialBasisFctMapping.hpp"
 #include "mapping/impl/BasisFunctions.hpp"
@@ -41,7 +42,7 @@ void addSubtagsToParents(std::list<xml::XMLTag> &subtags,
 
 // this function uses std::variant in order to add attributes of any type (double, string, bool)
 // to all tags in the list of tags \p storage.
-using variant_t = std::variant<xml::XMLAttribute<double>, xml::XMLAttribute<std::string>, xml::XMLAttribute<bool>>;
+using variant_t = std::variant<xml::XMLAttribute<double>, xml::XMLAttribute<std::string>, xml::XMLAttribute<bool>, xml::XMLAttribute<int>>;
 template <typename TagStorage>
 void addAttributes(TagStorage &storage, const std::vector<variant_t> &attributes)
 {
@@ -54,7 +55,8 @@ void addAttributes(TagStorage &storage, const std::vector<variant_t> &attributes
 // Enum required for the RBF instantiations
 enum struct RBFBackend {
   Eigen,
-  PETSc
+  PETSc,
+  PUM
 };
 
 // Helper in order to resolve the template instantiations.
@@ -78,6 +80,12 @@ struct BackendSelector<RBFBackend::PETSc, RBF> {
 };
 #endif
 
+// Specialization for the RBF PUM backend
+template <typename RBF>
+struct BackendSelector<RBFBackend::PUM, RBF> {
+  typedef mapping::PartitionOfUnityMapping<RBF> type;
+};
+
 // Variant holding all available RBF classes
 using rbf_variant_t = std::variant<CompactPolynomialC0, CompactPolynomialC2, CompactPolynomialC4, CompactPolynomialC6, CompactThinPlateSplinesC2, ThinPlateSplines, VolumeSplines, Multiquadrics, InverseMultiquadrics, Gaussian>;
 
@@ -161,6 +169,8 @@ MappingConfiguration::MappingConfiguration(
       XMLTag{*this, TYPE_RBF_GLOBAL_DIRECT, occ, TAG}.setDocumentation("Radial-basis-function mapping using a direct solver with a gather-scatter parallelism.")};
   std::list<XMLTag> rbfIterativeTags{
       XMLTag{*this, TYPE_RBF_GLOBAL_ITERATIVE, occ, TAG}.setDocumentation("Radial-basis-function mapping using an iterative solver with a distributed parallelism.")};
+  std::list<XMLTag> pumDirectTags{
+      XMLTag{*this, TYPE_RBF_PUM_DIRECT, occ, TAG}.setDocumentation("Radial-basis-function mapping using a partition of unity method, which supports a distributed parallelism.")};
   std::list<XMLTag> rbfAliasTag{
       XMLTag{*this, TYPE_RBF_ALIAS, occ, TAG}.setDocumentation("Alias tag, which auto-selects a radial-basis-function mapping depending on the simulation parameter,")};
 
@@ -188,17 +198,28 @@ MappingConfiguration::MappingConfiguration(
   auto attrPolynomial = makeXMLAttribute(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE)
                             .setDocumentation("Toggles use of the global polynomial")
                             .setOptions({POLYNOMIAL_ON, POLYNOMIAL_OFF, POLYNOMIAL_SEPARATE});
+  auto attrPumPolynomial = makeXMLAttribute(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE)
+                               .setDocumentation("Toggles use a local (per cluster) polynomial")
+                               .setOptions({POLYNOMIAL_OFF, POLYNOMIAL_SEPARATE});
 
   auto attrSolverRtol = makeXMLAttribute(ATTR_SOLVER_RTOL, 1e-9)
                             .setDocumentation("Solver relative tolerance for convergence");
   auto attrPreallocation = makeXMLAttribute(ATTR_PREALLOCATION, PREALLOCATION_TREE)
                                .setDocumentation("Sets kind of preallocation for PETSc RBF implementation")
                                .setOptions({PREALLOCATION_ESTIMATE, PREALLOCATION_COMPUTE, PREALLOCATION_OFF, PREALLOCATION_SAVE, PREALLOCATION_TREE});
 
+  auto verticesPerCluster = XMLAttribute<int>(ATTR_VERTICES_PER_CLUSTER, 100)
+                                .setDocumentation("Average number of vertices per cluster (partition) applied in the rbf partition of unity method.");
+  auto relativeOverlap = makeXMLAttribute(ATTR_RELATIVE_OVERLAP, 0.3)
+                             .setDocumentation("Value between 0 and 1 indicating the relative overlap between clusters. A value of 0.3 is usually a good trade-off between accuracy and efficiency.");
+  auto projectToInput = XMLAttribute<bool>(ATTR_PROJECT_TO_INPUT, true)
+                            .setDocumentation("If enabled, places the cluster centers at the closest vertex of the input mesh. Should be enabled in case of non-uniform point distributions such as for shell structures.");
+
   // Add the relevant attributes to the relevant tags
   addAttributes(projectionTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint});
   addAttributes(rbfDirectTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPolynomial, attrXDead, attrYDead, attrZDead});
   addAttributes(rbfIterativeTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPolynomial, attrXDead, attrYDead, attrZDead, attrSolverRtol, attrPreallocation});
+  addAttributes(pumDirectTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPumPolynomial, attrXDead, attrYDead, attrZDead, verticesPerCluster, relativeOverlap, projectToInput});
   addAttributes(rbfAliasTag, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrXDead, attrYDead, attrZDead});
 
   // Now we take care of the subtag basis function
@@ -217,6 +238,7 @@ MappingConfiguration::MappingConfiguration(
   addAttributes(supportRadiusRBF, {attrSupportRadius});
   addSubtagsToParents(supportRadiusRBF, rbfIterativeTags);
   addSubtagsToParents(supportRadiusRBF, rbfDirectTags);
+  addSubtagsToParents(supportRadiusRBF, pumDirectTags);
   addSubtagsToParents(supportRadiusRBF, rbfAliasTag);
 
   // Now the tags using a shape parameter
@@ -230,6 +252,7 @@ MappingConfiguration::MappingConfiguration(
   addAttributes(shapeParameterRBF, {attrShapeParam});
   addSubtagsToParents(shapeParameterRBF, rbfIterativeTags);
   addSubtagsToParents(shapeParameterRBF, rbfDirectTags);
+  addSubtagsToParents(shapeParameterRBF, pumDirectTags);
   addSubtagsToParents(shapeParameterRBF, rbfAliasTag);
 
   // For the Gaussian, we need default values as the user can pass a support radius or a shape parameter
@@ -240,6 +263,7 @@ MappingConfiguration::MappingConfiguration(
   addAttributes(GaussRBF, {attrShapeParam, attrSupportRadius});
   addSubtagsToParents(GaussRBF, rbfIterativeTags);
   addSubtagsToParents(GaussRBF, rbfDirectTags);
+  addSubtagsToParents(GaussRBF, pumDirectTags);
   addSubtagsToParents(GaussRBF, rbfAliasTag);
 
   // tags without an attribute
@@ -249,12 +273,14 @@ MappingConfiguration::MappingConfiguration(
 
   addSubtagsToParents(attributelessRBFs, rbfIterativeTags);
   addSubtagsToParents(attributelessRBFs, rbfDirectTags);
+  addSubtagsToParents(attributelessRBFs, pumDirectTags);
   addSubtagsToParents(attributelessRBFs, rbfAliasTag);
 
   // Add all tags to the mapping tag
   parent.addSubtags(projectionTags);
   parent.addSubtags(rbfIterativeTags);
   parent.addSubtags(rbfDirectTags);
+  parent.addSubtags(pumDirectTags);
   parent.addSubtags(rbfAliasTag);
 }
 
@@ -279,6 +305,11 @@ void MappingConfiguration::xmlTagCallback(
     std::string strPolynomial = tag.getStringAttributeValue(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE);
     std::string strPrealloc   = tag.getStringAttributeValue(ATTR_PREALLOCATION, PREALLOCATION_TREE);
 
+    // pum related tags
+    int    verticesPerCluster = tag.getIntAttributeValue(ATTR_VERTICES_PER_CLUSTER, 100);
+    double relativeOverlap    = tag.getDoubleAttributeValue(ATTR_RELATIVE_OVERLAP, 0.3);
+    bool   projectToInput     = tag.getBooleanAttributeValue(ATTR_PROJECT_TO_INPUT, true);
+
     // Convert raw string into enum types as the constructors take enums
     if (constraint == CONSTRAINT_CONSERVATIVE) {
       constraintValue = Mapping::CONSERVATIVE;
@@ -294,7 +325,7 @@ void MappingConfiguration::xmlTagCallback(
 
     ConfiguredMapping configuredMapping = createMapping(dir, type, fromMesh, toMesh);
 
-    _rbfConfig = configureRBFMapping(type, context, strPolynomial, strPrealloc, xDead, yDead, zDead, solverRtol);
+    _rbfConfig = configureRBFMapping(type, context, strPolynomial, strPrealloc, xDead, yDead, zDead, solverRtol, verticesPerCluster, relativeOverlap, projectToInput);
 
     checkDuplicates(configuredMapping);
     _mappings.push_back(configuredMapping);
@@ -340,6 +371,8 @@ void MappingConfiguration::xmlTagCallback(
 #else
       PRECICE_CHECK(false, "The global-iterative RBF solver requires a preCICE build with PETSc enabled.");
 #endif
+    } else if (_rbfConfig.solver == RBFConfiguration::SystemSolver::PUMDirect) {
+      mapping.mapping = getRBFMapping<RBFBackend::PUM>(basisFunction, constraintValue, mapping.fromMesh->getDimensions(), supportRadius, shapeParameter, _rbfConfig.deadAxis, _rbfConfig.polynomial, _rbfConfig.verticesPerCluster, _rbfConfig.relativeOverlap, _rbfConfig.projectToInput);
     } else {
       PRECICE_UNREACHABLE("Unknown RBF solver.");
     }
@@ -351,14 +384,19 @@ MappingConfiguration::RBFConfiguration MappingConfiguration::configureRBFMapping
                                                                                  const std::string &              polynomial,
                                                                                  const std::string &              preallocation,
                                                                                  bool xDead, bool yDead, bool zDead,
-                                                                                 double solverRtol) const
+                                                                                 double solverRtol,
+                                                                                 double verticesPerCluster,
+                                                                                 double relativeOverlap,
+                                                                                 bool   projectToInput) const
 {
   RBFConfiguration rbfConfig;
 
   if (type == TYPE_RBF_GLOBAL_ITERATIVE)
     rbfConfig.solver = RBFConfiguration::SystemSolver::GlobalIterative;
   else if (type == TYPE_RBF_GLOBAL_DIRECT)
     rbfConfig.solver = RBFConfiguration::SystemSolver::GlobalDirect;
+  else if (type == TYPE_RBF_PUM_DIRECT)
+    rbfConfig.solver = RBFConfiguration::SystemSolver::PUMDirect;
   else {
     // Rather simple auto-selection (for now)
     // The default is the Eigen backend, as it is always available. Only in certain situations, we will decide for the PETSc backend
@@ -401,6 +439,10 @@ MappingConfiguration::RBFConfiguration MappingConfiguration::configureRBFMapping
   rbfConfig.deadAxis   = {{xDead, yDead, zDead}};
   rbfConfig.solverRtol = solverRtol;
 
+  rbfConfig.verticesPerCluster = verticesPerCluster;
+  rbfConfig.relativeOverlap    = relativeOverlap;
+  rbfConfig.projectToInput     = projectToInput;
+
   return rbfConfig;
 }
 
@@ -488,7 +530,7 @@ const std::vector<MappingConfiguration::ConfiguredMapping> &MappingConfiguration
 
 bool MappingConfiguration::requiresBasisFunction(const std::string &mappingType) const
 {
-  return mappingType == TYPE_RBF_GLOBAL_DIRECT || mappingType == TYPE_RBF_GLOBAL_ITERATIVE || mappingType == TYPE_RBF_ALIAS;
+  return mappingType == TYPE_RBF_PUM_DIRECT || mappingType == TYPE_RBF_GLOBAL_DIRECT || mappingType == TYPE_RBF_GLOBAL_ITERATIVE || mappingType == TYPE_RBF_ALIAS;
 }
 
 BasisFunction MappingConfiguration::parseBasisFunctions(const std::string &basisFctName) const

src/mapping/config/MappingConfiguration.hpp

@@ -65,13 +65,17 @@ class MappingConfiguration : public xml::XMLTag::Listener {
 
     enum struct SystemSolver {
       GlobalDirect,
-      GlobalIterative
+      GlobalIterative,
+      PUMDirect
     };
     SystemSolver        solver{};
     std::array<bool, 3> deadAxis{};
     Polynomial          polynomial{};
     double              solverRtol{};
     Preallocation       preallocation{};
+    int                 verticesPerCluster{};
+    double              relativeOverlap{};
+    bool                projectToInput{};
   };
 
   /// Returns the RBF configuration, which was configured at the latest.
@@ -99,6 +103,7 @@ class MappingConfiguration : public xml::XMLTag::Listener {
   const std::string TYPE_LINEAR_CELL_INTERPOLATION = "linear-cell-interpolation";
   const std::string TYPE_RBF_GLOBAL_DIRECT         = "rbf-global-direct";
   const std::string TYPE_RBF_GLOBAL_ITERATIVE      = "rbf-global-iterative";
+  const std::string TYPE_RBF_PUM_DIRECT            = "rbf-pum-direct";
   const std::string TYPE_RBF_ALIAS                 = "rbf";
 
   const std::string ATTR_DIRECTION  = "direction";
@@ -140,6 +145,11 @@ class MappingConfiguration : public xml::XMLTag::Listener {
   // const std::string PARALLELISM_GATHER_SCATTER = "gather-scatter";
   // const std::string PARALLELISM                = "distributed";
 
+  // For PUM
+  const std::string ATTR_VERTICES_PER_CLUSTER = "vertices-per-cluster";
+  const std::string ATTR_RELATIVE_OVERLAP     = "relative-overlap";
+  const std::string ATTR_PROJECT_TO_INPUT     = "project-to-input";
+
   // We declare the basis function as subtag
   const std::string SUBTAG_BASIS_FUNCTION = "basis-function";
   const std::string RBF_TPS               = "thin-plate-splines";
@@ -195,7 +205,10 @@ class MappingConfiguration : public xml::XMLTag::Listener {
                                        const std::string &              polynomial,
                                        const std::string &              preallocation,
                                        bool xDead, bool yDead, bool zDead,
-                                       double solverRtol) const;
+                                       double solverRtol,
+                                       double verticesPerCluster,
+                                       double relativeOverlap,
+                                       bool   projectToInput) const;
 
   /// Check whether a mapping to and from the same mesh already exists
   void checkDuplicates(const ConfiguredMapping &mapping);