Merge pull request #292 from ls1mardyn/static-irreg-dd

added functionality for static irregular rectilinear grids

examples/Mamico-couette/basic4.csv

@@ -0,0 +1,3 @@
+1,10
+1,1,1
+1

examples/Mamico-couette/ls1configNoCP.xml

@@ -56,13 +56,15 @@
 							<lower> <x>0</x> <y>0</y> <z>0</z> </lower>
 							<upper> <x>30</x> <y>30</y> <z>30</z> </upper>
 						</object>
-						<velocityAssigner type="EqualVelocityDistribution" enableRandomSeed="1"></velocityAssigner>
+						<velocityAssigner type="EqualVelocityDistribution" enableRandomSeed="true"></velocityAssigner>
 					</objectgenerator>
 				</generator>
 			</phasespacepoint>
 		</ensemble>
 		<algorithm>
-			<parallelisation type="DomainDecomposition"/>
+			<parallelisation type="StaticIrregDomainDecomposition">
+            <subdomainWeightsCSV>basic4.csv</subdomainWeightsCSV>
+			</parallelisation>
 			<datastructure type="LinkedCells">
 			<cellsInCutoffRadius>1</cellsInCutoffRadius>
 			</datastructure>

examples/all-options.xml

@@ -197,6 +197,27 @@
           <overlappingStartAtStep>5</overlappingStartAtStep> <!-- Start overlapping at given step (default: 5), only relevant if overlappingCollectives==True -->
           <overlappingP2P>False</overlappingP2P> <!-- Defines whether to use overlapping p2p communication or not. Default: False -->
       </parallelisation>
+      <parallelisation type="StaticIrregDomainDecomposition">
+          <!-- Provide CSV file with weights that determine the relative widths of each subdomain
+          EX: A file with its first line being 1,2,1 defines an x axis with subdomain lengths in the 1:2:1 ratio
+          The number of digits per column are the number of ranks in that column, and the number of columns need not be equal per row
+          The number of ranks per row are therefore calculated from the column length, and cannot be given explicitly
+          The program must be run with ranks equal to the ranks in the CSV (examples/Mamico-couette/basic4.csv must be run with 6 ranks, for example)
+          The actual subdomain lengths are calculated from the total length and given ratios
+          If not given, default behaviour is an equally spaced grid, same as DomainDecomposition
+            -->
+          <subdomainWeightsCSV>filename.csv</subdomainWeightsCSV>
+      </parallelisation>
+      <parallelisation type="GeneralDomainDecomposition">
+          <updateFrequency>INTEGER</updateFrequency>
+          <initialPhaseTime>INTEGER</initialPhaseTime><!--time for initial rebalancing phase-->
+          <initialPhaseFrequency>INTEGER</initialPhaseFrequency><!--frequency for initial rebalancing phase-->
+          <gridSize>STRING</gridSize><!--default: 0; if non-zero, the process boundaries are fixed to multiples of
+              gridSize. Comma separated string to define three different grid sizes for the different dimensions is
+              possible.-->
+          <loadBalancer type="STRING"><!--STRING...type of the load balancer, currently supported: ALL-->
+          </loadBalancer>
+      </parallelisation>
       <datastructure type="LinkedCells">
           <cellsInCutoffRadius>1</cellsInCutoffRadius>
           <!-- select traversal algorithm

src/Simulation.cpp

@@ -31,6 +31,7 @@
 #include "parallel/DomainDecomposition.h"
 #include "parallel/KDDecomposition.h"
 #include "parallel/GeneralDomainDecomposition.h"
+#include "parallel/StaticIrregDomainDecomposition.h"
 #endif
 
 #include "particleContainer/adapter/ParticlePairs2PotForceAdapter.h"
@@ -326,6 +327,17 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 #ifdef MAMICO_COUPLING	
 				coupling::interface::LS1StaticCommData::getInstance().getLocalCommunicator(),
 				coupling::interface::LS1StaticCommData::getInstance().getDomainGridDecomp()
+#endif
+				);
+			}
+			else if(parallelisationtype == "StaticIrregDomainDecomposition") {
+				delete _domainDecomposition;
+				_domainDecomposition = new StaticIrregDomainDecomposition(
+				_domain
+#ifdef MAMICO_COUPLING
+				,
+				coupling::interface::LS1StaticCommData::getInstance().getLocalCommunicator(),
+				coupling::interface::LS1StaticCommData::getInstance().getSubdomainWeights()
 #endif
 				);
 			}

src/parallel/StaticIrregDomainDecomposition.cpp

@@ -0,0 +1,153 @@
+/**
+ * @file StaticIrregDomainDecompostion.cpp
+ * @author amartyads
+ * @date 21.02.24
+ */
+
+#include "StaticIrregDomainDecomposition.h"
+#include "Domain.h"
+#include "utils/Logger.h"
+#include <fstream>
+#include <numeric>
+#include <sstream>
+
+StaticIrregDomainDecomposition::StaticIrregDomainDecomposition(Domain *domain)
+    : StaticIrregDomainDecomposition(domain, MPI_COMM_WORLD,
+                                     {std::vector<unsigned int>{}, {}, {}}) {}
+
+StaticIrregDomainDecomposition::StaticIrregDomainDecomposition(
+    Domain *domain, MPI_Comm comm,
+    const std::array<std::vector<unsigned int>, DIMgeom> &subdomainWeights)
+    : DomainDecomposition(comm, {(int)subdomainWeights[0].size(),
+                                 (int)subdomainWeights[1].size(),
+                                 (int)subdomainWeights[2].size()}),
+      _subdomainWeights(subdomainWeights), _domainLength{
+                                               domain->getGlobalLength(0),
+                                               domain->getGlobalLength(1),
+                                               domain->getGlobalLength(2)} {
+  if (_subdomainWeights[0].size() == 0 && _subdomainWeights[1].size() == 0 &&
+      _subdomainWeights[2].size() == 0) { // default behaviour, regular grid
+    /* If we have empty vectors for subdomain weights, we successfully
+    initialize a regular grid by leveraging DomainDecomposition's default
+    behaviour. However, we now need to use the updated _coords and _gridSize to
+    write these weights back to _subdomainWeights, since otherwise the box
+    dimension calculation will be broken. So for each axis, we insert a number
+    of ones equal to the number of subdomains, to denote that we have equal
+    subdomains on this axis, amounting to _gridSize.
+    */
+    for (int i = 0; i < 3; i++) {
+      _subdomainWeights[i].reserve(_gridSize[i]);
+      for (int j = 0; j < _gridSize[i]; j++) {
+        _subdomainWeights[i].push_back(1); // equal subdomains, so weight = 1
+      }
+    }
+  }
+  updateSubdomainDimensions();
+}
+
+void StaticIrregDomainDecomposition::readXML(XMLfileUnits &xmlconfig) {
+  DomainDecompMPIBase::readXML(xmlconfig);
+  // bypass DomainDecomposition readXML to avoid reading MPIGridDims
+
+  std::string filename = xmlconfig.getNodeValue_string("subdomainWeightsCSV");
+  if (!filename.empty()) {
+    updateSubdomainWeightsFromFile(filename);
+    for (int i = 0; i < _subdomainWeights.size(); ++i) {
+      _gridSize[i] = static_cast<int>(
+          _subdomainWeights[i]
+              .size()); //_gridSize still contains the number of ranks per
+                        // dimension, just not the actual "size" of subdomains
+    }
+    initMPIGridDims();           // to recalculate _coords
+    updateSubdomainDimensions(); // recalculate sizes from _coords
+  }
+}
+
+double StaticIrregDomainDecomposition::getBoundingBoxMin(int dimension,
+                                                         Domain *domain) {
+  return _boxMin[dimension];
+}
+
+double StaticIrregDomainDecomposition::getBoundingBoxMax(int dimension,
+                                                         Domain *domain) {
+  return _boxMax[dimension];
+}
+
+void StaticIrregDomainDecomposition::updateSubdomainDimensions() {
+  for (int i = 0; i < 3; i++) {
+    const auto backWeight =
+        std::reduce(_subdomainWeights[i].begin(),
+                    _subdomainWeights[i].begin() + _coords[i], 0u);
+    const auto totalWeight =
+        std::reduce(_subdomainWeights[i].begin() + _coords[i],
+                    _subdomainWeights[i].end(), backWeight);
+    Log::global_log->debug()
+        << "Dim: " << i << " totalWeight: " << totalWeight
+        << " backWeight: " << backWeight << " coords: " << _coords[0] << ", "
+        << _coords[1] << ", " << _coords[2] << std::endl;
+
+    // calculate box bounds from cumulative weights of previous ranks, and the
+    // weight of the current rank
+    _boxMin[i] =
+        static_cast<double>(backWeight) * _domainLength[i] / totalWeight;
+    _boxMax[i] =
+        _boxMin[i] + (static_cast<double>(_subdomainWeights[i][_coords[i]]) *
+                      _domainLength[i] / totalWeight);
+  }
+}
+
+void StaticIrregDomainDecomposition::updateSubdomainWeightsFromFile(
+    const std::string &filename) {
+  std::ifstream file(filename.c_str());
+  if (!file.good()) {
+    Log::global_log->fatal() << "CSV file to read domain decomposition from "
+                                "does not exist! Please check config file!";
+    Simulation::exit(5001);
+  }
+
+  std::string line;
+  for (int i = 0; i < 3; i++) { // only reads the first 3 lines, theoretically
+                                // the rest of the file can contain whatever
+    getline(file, line);
+    if (line.empty()) {
+      Log::global_log->fatal()
+          << "CSV has less than 3 lines! Please check CSV file!";
+      Simulation::exit(5002);
+    }
+    _subdomainWeights[i].clear();
+    std::stringstream ss(line);
+    if (!ss.good()) {
+      Log::global_log->fatal() << "EOF or I/O error occured on line " << i
+                               << " of CSV. Please check CSV file!";
+      Simulation::exit(5004);
+    }
+    while (ss.good()) {
+      int temp;
+      ss >> temp;
+      if (temp <= 0) {
+        Log::global_log->fatal() << "CSV has non-natural number! Only weights "
+                                    "> 0 allowed, please check CSV file!";
+        Simulation::exit(5003);
+      }
+      _subdomainWeights[i].push_back(temp);
+      if (ss.peek() == ',' || ss.peek() == ' ') // skip commas and spaces
+        ss.ignore();
+    }
+    if (_subdomainWeights[i].empty()) {
+      Log::global_log->fatal()
+          << "Weights empty, failed reading operation, please check CSV file!";
+      Simulation::exit(5005);
+    }
+  }
+  Log::global_log->info() << "Weights for subdomains for "
+                             "StaticIrregDomainDecomposition have been read"
+                          << std::endl;
+  for (int i = 0; i < 3; i++) {
+    std::stringstream ss;
+    for (auto w : _subdomainWeights[i]) {
+      ss << w << " ";
+    }
+    Log::global_log->info()
+        << "Weights for axis " << i << ": " << ss.str() << std::endl;
+  }
+}

src/parallel/StaticIrregDomainDecomposition.h

@@ -0,0 +1,137 @@
+/**
+ * @file StaticIrregDomainDecompostion.h
+ * @author amartyads
+ * @date 21.02.24
+ */
+
+#pragma once
+
+#include "DomainDecompMPIBase.h"
+#include "DomainDecomposition.h"
+
+/**
+ * Extends DomainDecomposition to implement static irregular grids
+ */
+class StaticIrregDomainDecomposition : public DomainDecomposition {
+public:
+  /**
+   * Default constructor. Passes default values to the main constructor.
+   *
+   * The constructor passes the default value of _subDomainWeights to the main
+   * constructor, so that the size() of individual vectors is zero, to trigger
+   * the default MPI coord setup from DomainDecomposition.
+   *
+   * @param *domain The domain object defined in Simulation, needed to extract
+   * the global simulation bounds.
+   */
+  StaticIrregDomainDecomposition(Domain *domain);
+
+  /**
+   * Main constructor. Passes values to DomainDecomposition, and inits
+   * class members.
+   *
+   * In case _subdomainWeights is blank (to trigger initial coord breakdown from
+   * DomainDecomposition), the generated default _gridSize is taken and used to
+   * populate the weights. The behaviour becomes identical to
+   * DomainDecomposition, with a regular equally-spaced grid.
+   *
+   * @note This constructor is directly called when MaMiCo coupling is used. In
+   * this case the communicator and the weights are provided by MaMiCo.
+   *
+   * @param *domain The domain object defined in Simulation, needed to extract
+   * the global simulation bounds.
+   * @param comm The local communicator for the simulation
+   * @param _subdomainWeights An array containing 3 vectors, each containing the
+   * numeric weights of each subdomain, ordered by axes.
+   */
+  StaticIrregDomainDecomposition(
+      Domain *domain, MPI_Comm comm,
+      const std::array<std::vector<unsigned int>, DIMgeom> &subdomainWeights);
+
+  /**
+   * Reads in XML configuration for StaticIrregDomainDecomposition.
+   *
+   * The only configuration allowed right now is a CSV file, which contains the
+   actual domain breakdown.
+   * Even though this class subclasses DomainDecomposition, it bypasses the
+   readXML() mehod of DomainDecomposition
+   * because MPIGridDims is supposed to be calculated from the CSV
+   configuration.
+   *
+   * The following xml object structure is handled by this method:
+   * \code{.xml}
+     <parallelisation type="StaticIrregDomainDecomposition">
+       <!-- structure handled by DomainDecompMPIBase -->
+       <subdomainWeightsCSV> STRING.csv </subdomainWeightsCSV>
+     </parallelisation>
+     \endcode
+   *
+   * A file with its first line being 1,2,1 defines an x axis with subdomain
+   lengths in the 1:2:1 ratio
+   * If file not given, default behaviour is an equally spaced grid, same as
+   DomainDecomposition
+   *
+   * @param &xmlconfig The xml node from which to read the CSV filename with
+   weights
+   */
+  void readXML(XMLfileUnits &xmlconfig) override;
+
+  double getBoundingBoxMin(int dimension, Domain *domain) override;
+
+  double getBoundingBoxMax(int dimension, Domain *domain) override;
+
+  /**
+   * Assuming _subdomainWeights is up-to-date, calculates bounds of
+   * current subdomain and updates _boxMin and _boxMax.
+   *
+   * For an explanation on what the weights signify, please see the
+   * documentation for the member _subdomainWeights.
+   */
+  void updateSubdomainDimensions();
+
+  /**
+   * Reads in the CSV file given by the XML config, and updates
+   * _subdomainWeights.
+   *
+   * The CSV file is expected to have 3 lines of comma-separated integers, with
+   * the integer signifying the "weight" (relative width) of the subdomain. The
+   * lines, in order, are expected to be the weights for the x, y and z
+   * dimension. Consequently, the number of integers for a dimension signify the
+   * number of ranks in that dimension, and is calculated as such.
+   *
+   * @param &filename The CSV file from which to read weights. Obtained from the
+   * XML config.
+   */
+  void updateSubdomainWeightsFromFile(const std::string &filename);
+
+private:
+  /**
+   * Stores the weights from the given CSV file.
+   *
+   * The weights denote the relative width of that subdomain relative to the
+   * others in the same dimension. Ex: if the weights in the x dimension are
+   * 1,3,1, and the size of the domain in x dimension is 100, the total weight
+   * is 1+3+1=5 and the divisions are 1/5, 3/5 and 1/5. Hence, all subdomains
+   * with coords (0,y,z) are 20 units in x direction, (1,y,z) are 60 units and
+   * (2,y,z) are 20 units. Weights are only relative for the dimension, and
+   * weights for different dimensions are independent.
+   */
+  std::array<std::vector<unsigned int>, 3> _subdomainWeights{{{}, {}, {}}};
+
+  /**
+   * Stores the start of the subdomain. Calculated by
+   * updateSubdomainDimensions().
+   */
+  std::array<double, 3> _boxMin{0, 0, 0};
+
+  /**
+   * Stores the end of the subdomain. Calculated by
+   * updateSubdomainDimensions().
+   */
+  std::array<double, 3> _boxMax{0, 0, 0};
+
+  /**
+   * Stores the domain lengths. Set in the constructor.
+   */
+  std::array<double, 3> _domainLength;
+};