Merge pull request #5781 from nilsvu/random_amr

Allow random hp AMR for testing

src/ParallelAlgorithms/Amr/Criteria/Random.cpp

@@ -3,42 +3,60 @@
 
 #include "ParallelAlgorithms/Amr/Criteria/Random.hpp"
 
+#include <cstddef>
 #include <random>
+#include <unordered_map>
 
 #include "Domain/Amr/Flag.hpp"
+#include "Utilities/Algorithm.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
 
 namespace amr::Criteria {
-Random::Random(const double do_something_fraction,
+Random::Random(std::unordered_map<amr::Flag, size_t> probability_weights,
                const size_t maximum_refinement_level)
-    : do_something_fraction_(do_something_fraction),
+    : probability_weights_(std::move(probability_weights)),
       maximum_refinement_level_(maximum_refinement_level) {}
 
 Random::Random(CkMigrateMessage* msg) : Criterion(msg) {}
 
 // NOLINTNEXTLINE(google-runtime-references)
 void Random::pup(PUP::er& p) {
   Criterion::pup(p);
-  p | do_something_fraction_;
+  p | probability_weights_;
   p | maximum_refinement_level_;
 }
 
-amr::Flag Random::random_flag(const size_t current_refinement_level) const {
+namespace detail {
+amr::Flag random_flag(
+    const std::unordered_map<amr::Flag, size_t>& probability_weights) {
+  const size_t total_weight =
+      alg::accumulate(probability_weights, 0_st,
+                      [](const size_t total, const auto& flag_and_weight) {
+                        return total + flag_and_weight.second;
+                      });
+  if (total_weight == 0) {
+    return amr::Flag::DoNothing;
+  }
+  if (probability_weights.size() == 1) {
+    return probability_weights.begin()->first;
+  }
+
   static std::random_device r;
   static const auto seed = r();
   static std::mt19937 generator(seed);
-  static std::uniform_real_distribution<> distribution(0.0, 1.0);
-
-  const double random_number = distribution(generator);
-  if (random_number > do_something_fraction_) {
-    return amr::Flag::DoNothing;
-  }
-  const double join_fraction =
-      current_refinement_level / static_cast<double>(maximum_refinement_level_);
-  if (random_number < join_fraction * do_something_fraction_) {
-    return amr::Flag::Join;
+  std::uniform_int_distribution<size_t> distribution{0, total_weight - 1};
+
+  const size_t random_number = distribution(generator);
+  size_t cumulative_weight = 0;
+  for (const auto& [flag, probability_weight] : probability_weights) {
+    cumulative_weight += probability_weight;
+    if (random_number < cumulative_weight) {
+      return flag;
+    }
   }
-  return amr::Flag::Split;
+  ERROR("Should never reach here");
 }
+}  // namespace detail
 
 PUP::able::PUP_ID Random::my_PUP_ID = 0;  // NOLINT
 }  // namespace amr::Criteria

src/ParallelAlgorithms/Amr/Criteria/Random.hpp

@@ -4,7 +4,9 @@
 #pragma once
 
 #include <cstddef>
+#include <limits>
 #include <pup.h>
+#include <unordered_map>
 
 #include "Domain/Amr/Flag.hpp"
 #include "Domain/Structure/ElementId.hpp"
@@ -16,49 +18,54 @@
 #include "Utilities/TMPL.hpp"
 
 namespace amr::Criteria {
+namespace detail {
+amr::Flag random_flag(
+    const std::unordered_map<amr::Flag, size_t>& probability_weights);
+}  // namespace detail
+
 /*!
- * \brief Randomly h-refine (or coarsen) an Element in each dimension.
+ * \brief Randomly refine (or coarsen) an Element in each dimension.
  *
- * \details Let \f$f\f$ be `ChangeRefinementFraction`, \f$L_{max}\f$ be
- * `MaximumRefinementLevel`, and \f$L_d\f$ be the current refinement level
- * of an Element in a particular dimension.  In each dimension, a random
- * number \f$r_d \in [0, 1]\f$ is generated.  If \f$r_d > f\f$ the refinement
- * flag is set to amr::Flags::DoNothing.  If \f$r_d < f L_d / L_{max}\f$
- * the refinement flag is set to amr::Flags::Join.  Otherwise the
- * refinement flag is set to amr::Flag::Split.
+ * You can specify a probability for each possible `amr::Flag`. It is evaluated
+ * in each dimension separately. Details:
  *
- * \note This criterion is primarily useful for testing the mechanics of
- * h-refinement
+ * - Probabilities are specified as integer weights. The probability for an
+ *   `amr::Flag` is its weight over the sum of all weights.
+ * - Flags with weight zero do not need to be specified.
+ * - If all weights are zero, `amr::Flag::DoNothing` is always chosen.
  */
 class Random : public Criterion {
  public:
-  /// The fraction of the time random refinement does changes the grid
-  struct ChangeRefinementFraction {
-    using type = double;
+  struct ProbabilityWeights {
+    using type = std::unordered_map<amr::Flag, size_t>;
     static constexpr Options::String help = {
-        "The fraction of the time that random refinement will change the "
-        "grid."};
-    static double lower_bound() { return 0.0; }
-    static double upper_bound() { return 1.0; }
+        "Possible refinement types and their probability, specified as integer "
+        "weights. The probability for a refinement type is its weight over the "
+        "sum of all weights. For example, set 'Split: 1' and 'DoNothing: 4' to "
+        "split each element with 20% probability. The refinement is evaluated "
+        "in each dimension separately."};
   };
 
-  /// The maximum allowed refinement level
+  /// The maximum allowed refinement level.
+  /// Can be deleted once the max refinement level is enforced globally as an
+  /// AMR policy.
   struct MaximumRefinementLevel {
     using type = size_t;
     static constexpr Options::String help = {
         "The maximum allowed refinement level."};
     static size_t upper_bound() { return ElementId<3>::max_refinement_level; }
   };
 
-  using options = tmpl::list<ChangeRefinementFraction, MaximumRefinementLevel>;
+  using options = tmpl::list<ProbabilityWeights, MaximumRefinementLevel>;
 
   static constexpr Options::String help = {
-      "Randomly h-refine (or coarsen) the grid"};
+      "Randomly refine (or coarsen) the grid"};
 
   Random() = default;
 
-  Random(const double do_something_fraction,
-         const size_t maximum_refinement_level);
+  explicit Random(
+      std::unordered_map<amr::Flag, size_t> probability_weights,
+      size_t maximum_refinement_level = std::numeric_limits<size_t>::max());
 
   /// \cond
   explicit Random(CkMigrateMessage* msg);
@@ -70,27 +77,29 @@ class Random : public Criterion {
 
   using argument_tags = tmpl::list<>;
 
-  template <typename Metavariables>
+  template <size_t Dim, typename Metavariables>
   auto operator()(Parallel::GlobalCache<Metavariables>& /*cache*/,
-                  const ElementId<Metavariables::volume_dim>& element_id) const;
+                  const ElementId<Dim>& element_id) const;
 
   void pup(PUP::er& p) override;
 
  private:
-  amr::Flag random_flag(size_t current_refinement_level) const;
-
-  double do_something_fraction_{0.0};
-  size_t maximum_refinement_level_{0};
+  std::unordered_map<amr::Flag, size_t> probability_weights_{};
+  size_t maximum_refinement_level_ = std::numeric_limits<size_t>::max();
 };
 
-template <typename Metavariables>
-auto Random::operator()(
-    Parallel::GlobalCache<Metavariables>& /*cache*/,
-    const ElementId<Metavariables::volume_dim>& element_id) const {
-  constexpr size_t volume_dim = Metavariables::volume_dim;
-  auto result = make_array<volume_dim>(amr::Flag::Undefined);
-  for (size_t d = 0; d < volume_dim; ++d) {
-    result[d] = random_flag(element_id.segment_ids()[d].refinement_level());
+template <size_t Dim, typename Metavariables>
+auto Random::operator()(Parallel::GlobalCache<Metavariables>& /*cache*/,
+                        const ElementId<Dim>& element_id) const {
+  auto result = make_array<Dim>(amr::Flag::Undefined);
+  for (size_t d = 0; d < Dim; ++d) {
+    result[d] = detail::random_flag(probability_weights_);
+    // Enforce max refinement level. Can be deleted once it's enforced globally.
+    if (result[d] == amr::Flag::Split and
+        element_id.segment_ids()[d].refinement_level() >=
+            maximum_refinement_level_) {
+      result[d] = amr::Flag::DoNothing;
+    }
   }
   return result;
 }

tests/InputFiles/ExampleExecutables/RandomAmr1D.yaml

@@ -18,7 +18,10 @@ Parallelization:
 Amr:
   Criteria:
     - Random:
-        ChangeRefinementFraction: 0.8
+        ProbabilityWeights:
+          Split: 2
+          Join: 1
+          DoNothing: 1
         MaximumRefinementLevel: 8
   Policies:
     Isotropy: Anisotropic

tests/InputFiles/ExportCoordinates/Input1D.yaml

@@ -17,7 +17,10 @@ Parallelization:
 Amr:
   Criteria:
     - Random:
-        ChangeRefinementFraction: 0.8
+        ProbabilityWeights:
+          Split: 2
+          Join: 1
+          DoNothing: 1
         MaximumRefinementLevel: 4
   Policies:
     Isotropy: Anisotropic

tests/Unit/ParallelAlgorithms/Amr/Actions/Test_EvaluateRefinementCriteria.cpp

@@ -5,7 +5,9 @@
 
 #include <array>
 #include <cstddef>
+#include <limits>
 #include <memory>
+#include <unordered_map>
 #include <unordered_set>
 
 #include "DataStructures/DataBox/DataBox.hpp"
@@ -37,15 +39,16 @@ namespace {
 
 // when called on the specified refinement level, this criteria
 // always will choose to join
-auto create_always_join(const size_t refinement_level) {
-  return std::make_unique<amr::Criteria::Random>(1.0, refinement_level);
+auto create_always_join() {
+  return std::make_unique<amr::Criteria::Random>(
+      std::unordered_map<amr::Flag, size_t>{{amr::Flag::Join, 1}});
 }
 
 // when called on any refinement level, this criteria always will choose to do
 // nothing
 auto create_always_do_nothing() {
   return std::make_unique<amr::Criteria::Random>(
-      0.0, ElementId<3>::max_refinement_level);
+      std::unordered_map<amr::Flag, size_t>{{amr::Flag::DoNothing, 1}});
 }
 
 template <typename Metavariables>
@@ -291,17 +294,17 @@ SPECTRE_TEST_CASE("Unit.Amr.Actions.EvaluateRefinementCriteria",
   // Run the test 3 times, twice with a single criterion that give known
   // decisions, and then once with two criteria, one of which always produces
   // flags of a higher priority than the other
-  criteria.emplace_back(create_always_join(1));
+  criteria.emplace_back(create_always_join());
   evaluate_criteria(std::move(criteria), std::array{amr::Flag::Join});
   criteria.clear();
   criteria.emplace_back(create_always_do_nothing());
   evaluate_criteria(std::move(criteria), std::array{amr::Flag::DoNothing});
   criteria.clear();
   criteria.emplace_back(create_always_do_nothing());
-  criteria.emplace_back(create_always_join(1));
+  criteria.emplace_back(create_always_join());
   evaluate_criteria(std::move(criteria), std::array{amr::Flag::DoNothing});
   criteria.clear();
-  criteria.emplace_back(create_always_join(1));
+  criteria.emplace_back(create_always_join());
   criteria.emplace_back(create_always_do_nothing());
   evaluate_criteria(std::move(criteria), std::array{amr::Flag::DoNothing});
   check_split_while_join_is_avoided();