Merge pull request #758 from SeisSol/lts-wiggle-factor

Add LTS Wiggle Factor

Documentation/conf.py

@@ -20,7 +20,7 @@
 # -- Project information -----------------------------------------------------
 
 project = 'SeisSol'
-copyright = '2018, The SeisSol Team'
+copyright = '2023, The SeisSol Team'
 author = 'The SeisSol Team'
 
 # The short X.Y version

Documentation/index.rst

@@ -53,6 +53,7 @@ We gratefully acknowledge the funding of the German Research Foundation (as part
   configuration
   parameter-file
   initial-condition
+  local-timestepping
   left-lateral-right-lateral-normal-reverse.rst
   easi
   fault-tagging

Documentation/local-timestepping.rst

@@ -0,0 +1,90 @@
+Local time-stepping (LTS)
+===================================
+
+You can (and should!) enable local time-stepping for your simulations.
+Generally, this can lead to a better time to solution.
+The following settings are relevant:
+
+.. code-block:: Fortran
+
+    &Discretization
+    ...
+    ClusteredLTS = 2
+
+    LtsWiggleFactorMin = 0.51
+    LtsWiggleFactorStepsize = 0.01
+    LtsWiggleFactorEnforceMaximumDifference = 1
+
+    LtsMaxNumberOfClusters = 20
+    LtsAutoMergeClusters = 1
+    LtsAllowedRelativePerformanceLossAutoMerge = 0.01
+    /
+
+To enable LTS, use the setting :code:`ClusteredLTS = 2`.
+To disable it, use the value :code:`ClusteredLTS = 1`
+Other values are currently not supported.
+
+When using local time-stepping, SeisSol updates elements only as often as needed.
+To do this, SeisSol computes a time step size for each element independently.
+The elements are then grouped into so-called time-clusters, which are updated together.
+The time step size of a cluster is the minimum of the time step sizes of the elements in the cluster.
+Assuming rate-2 LTS (:code:`ClusteredLTS = 2`), the first cluster contains elements of time step sizes in the interval
+
+.. math::
+
+    [\lambda (\Delta t)^\text{min}, 2 \lambda (\Delta t)^\text{min}])
+
+the second cluster elements of size
+
+.. math::
+
+    [2 \lambda (\Delta t)^\text{min}, 4 \lambda (\Delta t)^\text{min}])
+
+and so on.
+The wiggle factor :math:`\lambda` is typically one.
+
+Maximum Difference Property
+----------------------------
+
+SeisSol enforces some constraints on the clustering, for example, neighboring elements are always either in the same cluster,
+or in a cluster which has at most a time step size difference of 2.
+Elements that are connected by a dynamic rupture face have to be in the same time-cluster.
+This is called the maximum difference property.
+
+
+Wiggle factor (experimental)
+----------------------------
+This feature is only supported for rate-2 LTS (:code:`ClusteredLTS = 2`) at the moment.
+The *LtsWiggleFactorMin* parameter sets the minimum allowed value for the wiggle factor :math:`0.5 < \lambda \leq 1`.
+This wiggle factor can be used to reduce the overall number of time-steps in some cases and hence reduce the cost of the simulation.
+Even though it seems counterproductive to reduce the global time step size, it can move the boundaries of the clusters such that
+elements move from a cluster with a smaller time step size to clusters with a larger time step size.
+
+SeisSol tries to find the optimal wiggle factor automatically by performing a grid search in the interval
+
+.. math::
+
+    \text{LtsWiggleFactorMin} \leq \lambda \leq 1
+
+The step size of the grid search is controlled by the parameter *LtsWiggleFactorStepsize*.
+The optimal wiggle factor is the one that minimizes the cost of updates per unit time.
+When the setting *LtsWiggleFactorEnforceMaximumDifference* is set to one, SeisSol enforces the `Maximum Difference Property`_.
+during the grid search. This leads to a better cost estimate, but computing this cost estimate can be costly for large scale simulations with many MPI ranks.
+Hence, this is a trade-off between the time required to initialize the simulation and the time required to perform the actual simulation.
+Typically this setting should be activated and only be deactivated when the initialization time becomes a bottleneck.
+
+The wiggle factor was inspired by the implementation in (Breuer, Heinecke, 2022) [1]_
+
+Enforcing maximum number of clusters (experimental)
+----------------------------------------------------
+You can set a maximum number of clusters by the parameter LtsMaxNumberOfClusters.
+This can lead to better performance in some cases, especially on GPUs.
+You can set a maximum number of clusters by setting :code:`LtsMaxNumberOfClusters=20`.
+SeisSol also supports the automatic merging of clusters.
+For this, you need to set *LtsAutoMergeClusters* to one.
+The parameter *LtsAllowedRelativePerformanceLossAutoMerge* controls the allowed relative performance loss when merging clusters compared
+to the baseline cost without cluster-merging and wiggle factor.
+
+These features should be considered experimental at this point.
+
+.. [1] Breuer, A., & Heinecke, A. (2022). Next-Generation Local Time Stepping for the ADER-DG Finite Element Method. In 2022 IEEE International Parallel and Distributed Processing Symposium (IPDPS) (pp. 402-413). IEEE.

Documentation/parameters.par

@@ -95,6 +95,17 @@ FixTimeStep = 5                      ! Manually chosen maximum time step
 ClusteredLTS = 2                     ! 1 for Global time stepping, 2,3,5,... Local time stepping (advised value 2)
 !ClusteredLTS defines the multi-rate for the time steps of the clusters 2 for Local time stepping
 LtsWeightTypeId = 1                  ! 0=exponential, 1=exponential-balanced, 2=encoded
+
+! Wiggle factor settings:
+! Wiggle factor adjusts time step size by a small factor. This can lead to a slightly better clustering.
+LtsWiggleFactorMin = 1.0 ! Minimal wiggle factor applied to time step size. Should be > 1/rate
+LtsWiggleFactorStepsize = 0.01 ! Stepsize for wiggle factor grid search
+LtsWiggleFactorEnforceMaximumDifference = 1 ! 0 or 1: Enforces the maximum difference between neighboring clusters during wiggle factor search
+LtsMaxNumberOfClusters = 20 ! Enforces a maximal number of clusters
+LtsAutoMergeClusters = 0 !  0 or 1: Activates auto merging of clusters
+LtsAllowedRelativePerformanceLossAutoMerge = 0.1 ! Find minimal max number of clusters such that new computational cost is at most increased by this factor
+
+
 /
 
 &Output

src/Geometry/MeshReaderFBinding.cpp

@@ -344,12 +344,17 @@ void read_mesh_puml_c(const char* meshfile,
 		logError() << error.what();
 	}
 
-	auto ltsWeights = getLtsWeightsImplementation(ltsWeightsType, config);
-	auto meshReader = new seissol::PUMLReader(meshfile, maximumAllowedTimeStep, checkPointFile,
-        ltsWeights.get(), tpwgt, readPartitionFromFile);
-	seissol::SeisSol::main.setMeshReader(meshReader);
-
-	read_mesh(rank, seissol::SeisSol::main.meshReader(), hasFault, displacement, scalingMatrix);
+        const auto* ltsParameters = seissol::SeisSol::main.getMemoryManager().getLtsParameters();
+        auto ltsWeights = getLtsWeightsImplementation(ltsWeightsType, config, ltsParameters);
+        auto meshReader = new seissol::PUMLReader(meshfile,
+                                                  maximumAllowedTimeStep,
+                                                  checkPointFile,
+                                                  ltsWeights.get(),
+                                                  tpwgt,
+                                                  readPartitionFromFile);
+        seissol::SeisSol::main.setMeshReader(meshReader);
+
+        read_mesh(rank, seissol::SeisSol::main.meshReader(), hasFault, displacement, scalingMatrix);
 
 	watch.pause();
 	watch.printTime("Mesh initialized in:");

src/Initializer/MemoryManager.h

@@ -85,6 +85,8 @@
 #include <Initializer/DynamicRupture.h>
 #include <Initializer/Boundary.h>
 #include <Initializer/ParameterDB.h>
+#include <Initializer/time_stepping/LtsParameters.h>
+
 
 #include <DynamicRupture/Factory.h>
 #include <yaml-cpp/yaml.h>
@@ -165,6 +167,7 @@ class seissol::initializers::MemoryManager {
     std::unique_ptr<dr::output::OutputManager> m_faultOutputManager = nullptr;
     std::shared_ptr<dr::DRParameters> m_dynRupParameters = nullptr;
     std::shared_ptr<YAML::Node> m_inputParams = nullptr;
+    std::shared_ptr<time_stepping::LtsParameters> ltsParameters = nullptr;
 
     LTSTree m_boundaryTree;
     Boundary m_boundary;
@@ -351,9 +354,14 @@ class seissol::initializers::MemoryManager {
         return m_dynRupParameters.get();
     }
 
+    inline time_stepping::LtsParameters* getLtsParameters() {
+        return ltsParameters.get();
+    };
+
     void setInputParams(std::shared_ptr<YAML::Node> params) {
       m_inputParams = params;
       m_dynRupParameters = dr::readParametersFromYaml(m_inputParams);
+      ltsParameters = std::make_shared<time_stepping::LtsParameters>(time_stepping::readLtsParametersFromYaml(m_inputParams));
     }
 
 #ifdef ACL_DEVICE

src/Initializer/time_stepping/LtsLayout.cpp

@@ -1191,6 +1191,11 @@ void seissol::initializers::time_stepping::LtsLayout::deriveLayout( enum TimeClu
   double l_perCellSpeedup, l_clusteringSpeedup;
   getTheoreticalSpeedup( l_perCellSpeedup, l_clusteringSpeedup );
 
+  // The speedups above are computed without considering the wiggle factor
+  const auto wiggleFactor = seissol::SeisSol::main.wiggleFactorLts;
+  l_perCellSpeedup *= wiggleFactor;
+  l_clusteringSpeedup *= wiggleFactor;
+
   // get maximum speedup
   logInfo(rank) << "maximum theoretical speedup (compared to GTS):"
                   << l_perCellSpeedup << "per cell LTS," << l_clusteringSpeedup << "with the used clustering.";

src/Initializer/time_stepping/LtsParameters.cpp

@@ -0,0 +1,85 @@
+#include "Initializer/InputAux.hpp"
+#include "LtsParameters.h"
+#include <utils/logger.h>
+
+namespace seissol::initializers::time_stepping {
+
+LtsParameters readLtsParametersFromYaml(std::shared_ptr<YAML::Node>& params) {
+  using namespace seissol::initializers;
+
+  const auto discretizationParams = (*params)["discretization"];
+  const unsigned int rate = getWithDefault(discretizationParams, "clusteredlts", 1);
+  const double wiggleFactorMinimum =
+      getWithDefault(discretizationParams, "ltswigglefactormin", 1.0);
+  const double wiggleFactorStepsize =
+      getWithDefault(discretizationParams, "ltswigglefactorstepsize", 0.01);
+  const bool wiggleFactorEnforceMaximumDifference =
+      getWithDefault(discretizationParams, "ltswigglefactorenforcemaximumdifference", true);
+  const unsigned int maxNumberOfClusters = getWithDefault(
+      discretizationParams, "ltsmaxnumberofclusters", std::numeric_limits<int>::max() - 1);
+  const bool autoMergeClusters =
+      getWithDefault(discretizationParams, "ltsautomergeclusters", false);
+  const double allowedRelativePerformanceLossAutoMerge =
+      getWithDefault(discretizationParams, "ltsallowedrelativeperformancelossautomerge", 0.0);
+  const double allowedPerformanceLossRatioAutoMerge = allowedRelativePerformanceLossAutoMerge + 1.0;
+  return LtsParameters(rate,
+                       wiggleFactorMinimum,
+                       wiggleFactorStepsize,
+                       wiggleFactorEnforceMaximumDifference,
+                       maxNumberOfClusters,
+                       autoMergeClusters,
+                       allowedPerformanceLossRatioAutoMerge);
+}
+
+LtsParameters::LtsParameters(unsigned int rate,
+                             double wiggleFactorMinimum,
+                             double wiggleFactorStepsize,
+                             bool wigleFactorEnforceMaximumDifference,
+                             int maxNumberOfClusters,
+                             bool autoMergeClusters,
+                             double allowedPerformanceLossRatioAutoMerge)
+    : rate(rate), wiggleFactorMinimum(wiggleFactorMinimum),
+      wiggleFactorStepsize(wiggleFactorStepsize),
+      wiggleFactorEnforceMaximumDifference(wigleFactorEnforceMaximumDifference),
+      maxNumberOfClusters(maxNumberOfClusters),
+      autoMergeClusters(autoMergeClusters),
+      allowedPerformanceLossRatioAutoMerge(allowedPerformanceLossRatioAutoMerge) {
+  const bool isWiggleFactorValid =
+      (rate == 1 && wiggleFactorMinimum == 1.0) ||
+      (wiggleFactorMinimum <= 1.0 && wiggleFactorMinimum > (1.0 / rate));
+  if (!isWiggleFactorValid) {
+    logError() << "Minimal wiggle factor of " << wiggleFactorMinimum << "is not valid for rate"
+               << rate;
+  }
+  if (maxNumberOfClusters <= 0) {
+    logError() << "At least one cluster is required. Settings ltsMaxNumberOfClusters is invalid.";
+  }
+  if (allowedPerformanceLossRatioAutoMerge < 1.0) {
+    logError() << "Negative performance loss for auto merge is invalid.";
+  }
+}
+
+bool LtsParameters::isWiggleFactorUsed() const { return wiggleFactorMinimum < 1.0; }
+
+unsigned int LtsParameters::getRate() const { return rate; }
+
+double LtsParameters::getWiggleFactorMinimum() const { return wiggleFactorMinimum; }
+
+double LtsParameters::getWiggleFactorStepsize() const { return wiggleFactorStepsize; }
+
+bool LtsParameters::getWiggleFactorEnforceMaximumDifference() const {
+  return wiggleFactorEnforceMaximumDifference;
+}
+
+int LtsParameters::getMaxNumberOfClusters() const { return maxNumberOfClusters; }
+
+bool LtsParameters::isAutoMergeUsed() const {
+  return autoMergeClusters;
+}
+
+double LtsParameters::getAllowedPerformanceLossRatioAutoMerge() const {
+  return allowedPerformanceLossRatioAutoMerge;
+}
+
+
+} // namespace seissol::initializers::time_stepping

src/Initializer/time_stepping/LtsParameters.h

@@ -0,0 +1,42 @@
+#ifndef SEISSOL_LTSCONFIGURATION_H
+#define SEISSOL_LTSCONFIGURATION_H
+
+#include <memory>
+#include <yaml-cpp/yaml.h>
+
+namespace seissol::initializers::time_stepping {
+
+class LtsParameters {
+  private:
+  unsigned int rate;
+  double wiggleFactorMinimum;
+  double wiggleFactorStepsize;
+  bool wiggleFactorEnforceMaximumDifference;
+  unsigned int maxNumberOfClusters;
+  bool autoMergeClusters;
+  double allowedPerformanceLossRatioAutoMerge;
+
+  public:
+  [[nodiscard]] unsigned int getRate() const;
+  [[nodiscard]] bool isWiggleFactorUsed() const;
+  [[nodiscard]] double getWiggleFactorMinimum() const;
+  [[nodiscard]] double getWiggleFactorStepsize() const;
+  [[nodiscard]] bool getWiggleFactorEnforceMaximumDifference() const;
+  [[nodiscard]] int getMaxNumberOfClusters() const;
+  [[nodiscard]] bool isAutoMergeUsed() const;
+  [[nodiscard]] double getAllowedPerformanceLossRatioAutoMerge() const;
+
+  LtsParameters(unsigned int rate,
+                double wiggleFactorMinimum,
+                double wiggleFactorStepsize,
+                bool wigleFactorEnforceMaximumDifference,
+                int maxNumberOfClusters,
+                bool ltsAutoMergeClusters,
+                double allowedPerformanceLossRatioAutoMerge);
+};
+
+LtsParameters readLtsParametersFromYaml(std::shared_ptr<YAML::Node>& params);
+
+} // namespace seissol::initializers::time_stepping
+
+#endif // SEISSOL_LTSCONFIGURATION_H