Allow random hp AMR for testing

Support triggering any of the AMR flags with a specified probability, rather than
only h-refinement. Simplifies the reasoning of the AMR criterion
and gives more control over the types of refinement.
Still enforces a maximum refinement level, which can be
deleted once it's enforced globally with an AMR policy.

src/ParallelAlgorithms/Amr/Criteria/Random.cpp

@@ -3,42 +3,53 @@
 
 #include "ParallelAlgorithms/Amr/Criteria/Random.hpp"
 
+#include <cstddef>
+#include <map>
 #include <random>
 
 #include "Domain/Amr/Flag.hpp"
 
 namespace amr::Criteria {
-Random::Random(const double do_something_fraction,
-               const size_t maximum_refinement_level)
-    : do_something_fraction_(do_something_fraction),
-      maximum_refinement_level_(maximum_refinement_level) {}
+Random::Random(std::map<amr::Flag, double> probabilities,
+               const size_t maximum_refinement_level,
+               const Options::Context& context)
+    : probabilities_(std::move(probabilities)),
+      maximum_refinement_level_(maximum_refinement_level) {
+  if (std::any_of(
+          probabilities_.begin(), probabilities_.end(),
+          [](const auto& p) { return p.second < 0. or p.second > 1.; })) {
+    PARSE_ERROR(context, "Probabilities must be in the range [0, 1].");
+  }
+}
 
 Random::Random(CkMigrateMessage* msg) : Criterion(msg) {}
 
-// NOLINTNEXTLINE(google-runtime-references)
 void Random::pup(PUP::er& p) {
   Criterion::pup(p);
-  p | do_something_fraction_;
+  p | probabilities_;
   p | maximum_refinement_level_;
 }
 
-amr::Flag Random::random_flag(const size_t current_refinement_level) const {
+namespace detail {
+amr::Flag random_flag(const std::map<amr::Flag, double>& probabilities) {
   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;
+  for (auto it = probabilities.rbegin(); it != probabilities.rend(); ++it) {
+    const auto& [flag, probability] = *it;
+    if (probability == 1.) {
+      return flag;
+    } else if (probability == 0.) {
+      continue;
+    } else if (distribution(generator) < probability) {
+      return flag;
+    }
   }
-  return amr::Flag::Split;
+  return amr::Flag::DoNothing;
 }
+}  // namespace detail
 
 PUP::able::PUP_ID Random::my_PUP_ID = 0;  // NOLINT
 }  // namespace amr::Criteria

src/ParallelAlgorithms/Amr/Criteria/Random.hpp

@@ -4,6 +4,7 @@
 #pragma once
 
 #include <cstddef>
+#include <map>
 #include <pup.h>
 
 #include "Domain/Amr/Flag.hpp"
@@ -16,49 +17,50 @@
 #include "Utilities/TMPL.hpp"
 
 namespace amr::Criteria {
+namespace detail {
+amr::Flag random_flag(const std::map<amr::Flag, double>& probabilities);
+}  // namespace detail
+
 /*!
- * \brief Randomly h-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.
+ * \brief Randomly refine (or coarsen) an Element in each dimension.
  *
- * \note This criterion is primarily useful for testing the mechanics of
- * h-refinement
+ * \details You can specify a probability for each possible `amr::Flag`. The
+ * probabilities are evaluated in the priority order defined by `amr::Flag`. The
+ * first one that triggers is used in that dimension. If none of the
+ * probabilities trigger, then the Element is not refined.
  */
 class Random : public Criterion {
  public:
-  /// The fraction of the time random refinement does changes the grid
-  struct ChangeRefinementFraction {
-    using type = double;
+  struct Probabilities {
+    using type = std::map<amr::Flag, double>;
     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. "
+        "The probabilities are evaluated in order of refinement type "
+        "priority: first Split, then IncreaseResolution, then "
+        "DecreaseResolution, then Join. If none of the probabilities "
+        "trigger, then the Element is not refined."};
   };
 
-  /// 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<Probabilities, MaximumRefinementLevel>;
 
   static constexpr Options::String help = {
       "Randomly h-refine (or coarsen) the grid"};
 
   Random() = default;
 
-  Random(const double do_something_fraction,
-         const size_t maximum_refinement_level);
+  Random(std::map<amr::Flag, double> probabilities,
+         const size_t maximum_refinement_level,
+         const Options::Context& context = {});
 
   /// \cond
   explicit Random(CkMigrateMessage* msg);
@@ -70,27 +72,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};
+  std::map<amr::Flag, double> probabilities_{};
   size_t maximum_refinement_level_{0};
 };
 
-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(probabilities_);
+    // 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,9 @@ Parallelization:
 Amr:
   Criteria:
     - Random:
-        ChangeRefinementFraction: 0.8
+        Probabilities:
+          Split: 0.6
+          Join: 0.2
         MaximumRefinementLevel: 8
   Verbosity: Verbose
 

tests/InputFiles/ExportCoordinates/Input1D.yaml

@@ -17,7 +17,9 @@ Parallelization:
 Amr:
   Criteria:
     - Random:
-        ChangeRefinementFraction: 0.8
+        Probabilities:
+          Split: 0.6
+          Join: 0.2
         MaximumRefinementLevel: 4
   Verbosity: Verbose
 

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

@@ -36,15 +36,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::map<amr::Flag, double>{{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::map<amr::Flag, double>{{amr::Flag::DoNothing, 1.}});
 }
 
 template <typename Metavariables>
@@ -289,17 +290,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();

tests/Unit/ParallelAlgorithms/Amr/Criteria/Test_Random.cpp

@@ -48,26 +48,17 @@ void test_criterion(const amr::Criterion& criterion) {
   ElementId<VolumeDim> root_id{0};
   auto flags = criterion.evaluate(empty_box, empty_cache, root_id);
   for (size_t d = 0; d < VolumeDim; ++d) {
-    CHECK((flags[d] == amr::Flag::Split or flags[d] == amr::Flag::DoNothing));
-  }
-  for (size_t level = 1; level < 5; ++level) {
-    ElementId<VolumeDim> id{0, make_array<VolumeDim>(SegmentId(level, 1))};
-    flags = criterion.evaluate(empty_box, empty_cache, id);
-    for (size_t d = 0; d < VolumeDim; ++d) {
-      CHECK((flags[d] == amr::Flag::Split or flags[d] == amr::Flag::DoNothing or
-             flags[d] == amr::Flag::Join));
-    }
-  }
-  ElementId<VolumeDim> id{0, make_array<VolumeDim>(SegmentId(5, 1))};
-  flags = criterion.evaluate(empty_box, empty_cache, id);
-  for (size_t d = 0; d < VolumeDim; ++d) {
-    CHECK((flags[d] == amr::Flag::Join or flags[d] == amr::Flag::DoNothing));
+    CHECK((flags[d] == amr::Flag::Split or flags[d] == amr::Flag::Join or
+           flags[d] == amr::Flag::DoNothing));
   }
 }
 
 template <size_t VolumeDim>
-void test_always_change_refinement() {
-  const amr::Criteria::Random criterion(1.0, 5);
+void test_always_split() {
+  const amr::Criteria::Random criterion{
+      {{amr::Flag::Split, 1.},
+       // Split has higher priority, so this should never be triggered
+       {amr::Flag::IncreaseResolution, 1.}}};
   Parallel::GlobalCache<Metavariables<VolumeDim>> empty_cache{};
   auto databox = db::create<db::AddSimpleTags<>>();
   auto empty_box =
@@ -78,23 +69,11 @@ void test_always_change_refinement() {
   for (size_t d = 0; d < VolumeDim; ++d) {
     CHECK(flags[d] == amr::Flag::Split);
   }
-  for (size_t level = 1; level < 5; ++level) {
-    ElementId<VolumeDim> id{0, make_array<VolumeDim>(SegmentId(level, 1))};
-    flags = criterion.evaluate(empty_box, empty_cache, id);
-    for (size_t d = 0; d < VolumeDim; ++d) {
-      CHECK((flags[d] == amr::Flag::Split or flags[d] == amr::Flag::Join));
-    }
-  }
-  ElementId<VolumeDim> id{0, make_array<VolumeDim>(SegmentId(5, 1))};
-  flags = criterion.evaluate(empty_box, empty_cache, id);
-  for (size_t d = 0; d < VolumeDim; ++d) {
-    CHECK(flags[d] == amr::Flag::Join);
-  }
 }
 
 template <size_t VolumeDim>
 void test_always_do_nothing() {
-  const amr::Criteria::Random criterion(0.0, 5);
+  const amr::Criteria::Random criterion{{}};
 
   Parallel::GlobalCache<Metavariables<VolumeDim>> empty_cache{};
   auto databox = db::create<db::AddSimpleTags<>>();
@@ -110,21 +89,43 @@ void test_always_do_nothing() {
   }
 }
 
+template <size_t VolumeDim>
+void test_h_or_p() {
+  const amr::Criteria::Random criterion{{{amr::Flag::Split, 0.5},
+                                         {amr::Flag::IncreaseResolution, 0.5},
+                                         {amr::Flag::DecreaseResolution, 0.5},
+                                         {amr::Flag::Join, 1.}}};
+  Parallel::GlobalCache<Metavariables<VolumeDim>> empty_cache{};
+  auto databox = db::create<db::AddSimpleTags<>>();
+  auto empty_box =
+      make_observation_box<db::AddComputeTags<>>(make_not_null(&databox));
+
+  ElementId<VolumeDim> root_id{0};
+  auto flags = criterion.evaluate(empty_box, empty_cache, root_id);
+  for (size_t d = 0; d < VolumeDim; ++d) {
+    CHECK((flags[d] != amr::Flag::DoNothing and
+           flags[d] != amr::Flag::Undefined));
+  }
+}
+
 template <size_t VolumeDim>
 void test() {
   register_factory_classes_with_charm<Metavariables<VolumeDim>>();
-  test_always_change_refinement<VolumeDim>();
+  test_always_split<VolumeDim>();
   test_always_do_nothing<VolumeDim>();
-  const amr::Criteria::Random random_criterion(0.8, 5);
+  test_h_or_p<VolumeDim>();
+  const amr::Criteria::Random random_criterion{
+      {{amr::Flag::Split, 0.8}, {amr::Flag::Join, 0.1}}};
   test_criterion<VolumeDim>(random_criterion);
   test_criterion<VolumeDim>(serialize_and_deserialize(random_criterion));
 
   const auto criterion =
       TestHelpers::test_creation<std::unique_ptr<amr::Criterion>,
                                  Metavariables<VolumeDim>>(
           "Random:\n"
-          "  ChangeRefinementFraction: 0.8\n"
-          "  MaximumRefinementLevel: 5\n");
+          "  Probabilities:\n"
+          "    Split: 0.5\n"
+          "    Join: 0.3\n");
   test_criterion<VolumeDim>(*criterion);
   test_criterion<VolumeDim>(*serialize_and_deserialize(criterion));
 }