Merge pull request #1996 from nilsdeppe/domain_cleanups

Domain cleanups, operator== for functions of time, serialization of unique_ptr Base when nullptr

src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp

@@ -105,6 +105,12 @@ void PiecewisePolynomial<MaxDeriv>::DerivInfo::pup(PUP::er& p) noexcept {
   p | derivs_coefs;
 }
 
+template <size_t MaxDeriv>
+bool PiecewisePolynomial<MaxDeriv>::DerivInfo::operator==(
+    const PiecewisePolynomial<MaxDeriv>::DerivInfo& rhs) const noexcept {
+  return time == rhs.time and derivs_coefs == rhs.derivs_coefs;
+}
+
 template <size_t MaxDeriv>
 const typename PiecewisePolynomial<MaxDeriv>::DerivInfo&
 PiecewisePolynomial<MaxDeriv>::deriv_info_from_upper_bound(const double t) const
@@ -140,16 +146,37 @@ void PiecewisePolynomial<MaxDeriv>::pup(PUP::er& p) {
   p | deriv_info_at_update_times_;
 }
 
+template <size_t MaxDeriv>
+bool operator==(const PiecewisePolynomial<MaxDeriv>& lhs,
+                const PiecewisePolynomial<MaxDeriv>& rhs) noexcept {
+  return lhs.deriv_info_at_update_times_ == rhs.deriv_info_at_update_times_;
+}
+
+template <size_t MaxDeriv>
+bool operator!=(const PiecewisePolynomial<MaxDeriv>& lhs,
+                const PiecewisePolynomial<MaxDeriv>& rhs) noexcept {
+  return not(lhs == rhs);
+}
+
 // do explicit instantiation of MaxDeriv = {2,3,4}
 // along with all combinations of MaxDerivReturned = {0,...,MaxDeriv}
 /// \cond
-template class PiecewisePolynomial<2_st>;
-template class PiecewisePolynomial<3_st>;
-template class PiecewisePolynomial<4_st>;
-
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 #define DIMRETURNED(data) BOOST_PP_TUPLE_ELEM(1, data)
 
+#define INSTANTIATE(_, data)                                       \
+  template bool operator==                                         \
+      <DIM(data)>(const PiecewisePolynomial<DIM(data)>&,           \
+                  const PiecewisePolynomial<DIM(data)>&) noexcept; \
+  template class PiecewisePolynomial<DIM(data)>;                   \
+  template bool operator!=                                         \
+      <DIM(data)>(const PiecewisePolynomial<DIM(data)>&,           \
+                  const PiecewisePolynomial<DIM(data)>&) noexcept;
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (2, 3, 4))
+
+#undef INSTANTIATE
+
 #define INSTANTIATE(_, data)                                          \
   template std::array<DataVector, DIMRETURNED(data) + 1>              \
   PiecewisePolynomial<DIM(data)>::func_and_derivs<DIMRETURNED(data)>( \

src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp

@@ -67,6 +67,11 @@ class PiecewisePolynomial : public FunctionOfTime {
   void pup(PUP::er& p) override;
 
  private:
+  template <size_t LocalMaxDeriv>
+  friend bool operator==(  // NOLINT(readability-redundant-declaration)
+      const PiecewisePolynomial<LocalMaxDeriv>& lhs,
+      const PiecewisePolynomial<LocalMaxDeriv>& rhs) noexcept;
+
   /// Returns the function and `MaxDerivReturned` derivatives at
   /// an arbitrary time `t`.
   /// The function has multiple components.
@@ -93,6 +98,8 @@ class PiecewisePolynomial : public FunctionOfTime {
 
     // NOLINTNEXTLINE(google-runtime-references)
     void pup(PUP::er& p) noexcept;
+
+    bool operator==(const DerivInfo& rhs) const noexcept;
   };
 
   /// Returns a DerivInfo corresponding to the closest element in the range of
@@ -105,6 +112,10 @@ class PiecewisePolynomial : public FunctionOfTime {
   std::vector<DerivInfo> deriv_info_at_update_times_;
 };
 
+template <size_t MaxDeriv>
+bool operator!=(const PiecewisePolynomial<MaxDeriv>& lhs,
+                const PiecewisePolynomial<MaxDeriv>& rhs) noexcept;
+
 /// \cond
 template <size_t MaxDeriv>
 PUP::able::PUP_ID PiecewisePolynomial<MaxDeriv>::my_PUP_ID = 0;  // NOLINT

src/Domain/FunctionsOfTime/SettleToConstant.cpp

@@ -66,6 +66,18 @@ void SettleToConstant::pup(PUP::er& p) {
   p | inv_decay_time_;
 }
 
+bool operator==(const SettleToConstant& lhs,
+                const SettleToConstant& rhs) noexcept {
+  return lhs.coef_a_ == rhs.coef_a_ and lhs.coef_b_ == rhs.coef_b_ and
+         lhs.coef_c_ == rhs.coef_c_ and lhs.match_time_ == rhs.match_time_ and
+         lhs.inv_decay_time_ == rhs.inv_decay_time_;
+}
+
+bool operator!=(const SettleToConstant& lhs,
+                const SettleToConstant& rhs) noexcept {
+  return not(lhs == rhs);
+}
+
 /// \cond
 PUP::able::PUP_ID SettleToConstant::my_PUP_ID = 0;  // NOLINT
 

src/Domain/FunctionsOfTime/SettleToConstant.hpp

@@ -71,6 +71,9 @@ class SettleToConstant : public FunctionOfTime {
   void pup(PUP::er& p) override;
 
  private:
+  friend bool operator==(const SettleToConstant& lhs,
+                         const SettleToConstant& rhs) noexcept;
+
   template <size_t MaxDerivReturned = 2>
   std::array<DataVector, MaxDerivReturned + 1> func_and_derivs(double t) const
       noexcept;
@@ -79,5 +82,8 @@ class SettleToConstant : public FunctionOfTime {
   double match_time_{std::numeric_limits<double>::signaling_NaN()};
   double inv_decay_time_{std::numeric_limits<double>::signaling_NaN()};
 };
+
+bool operator!=(const SettleToConstant& lhs,
+                const SettleToConstant& rhs) noexcept;
 }  // namespace FunctionsOfTime
 }  // namespace domain

src/Parallel/PupStlCpp11.hpp

@@ -185,13 +185,17 @@ inline void pup(PUP::er& p, std::unique_ptr<T>& t) {  // NOLINT
 
 template <typename T, Requires<std::is_base_of<PUP::able, T>::value> = nullptr>
 inline void pup(PUP::er& p, std::unique_ptr<T>& t) {  // NOLINT
-  T* t1 = nullptr;
-  if (p.isUnpacking()) {
-    p | t1;
-    t = std::unique_ptr<T>(t1);
-  } else {
-    t1 = t.get();
-    p | t1;
+  bool is_nullptr = (nullptr == t);
+  p | is_nullptr;
+  if (not is_nullptr) {
+    T* t1 = nullptr;
+    if (p.isUnpacking()) {
+      p | t1;
+      t = std::unique_ptr<T>(t1);
+    } else {
+      t1 = t.get();
+      p | t1;
+    }
   }
 }
 // @}

tests/Unit/Domain/CoordinateMaps/TestMapHelpers.hpp

@@ -37,7 +37,7 @@ template <typename Map>
 bool are_maps_equal(
     const Map& map,
     const domain::CoordinateMapBase<Frame::Logical, Frame::Inertial, Map::dim>&
-        map_base) {
+        map_base) noexcept {
   const auto* map_derived = dynamic_cast<const Map*>(&map_base);
   return map_derived == nullptr ? false : (*map_derived == map);
 }
@@ -50,7 +50,7 @@ void check_if_maps_are_equal(
     const domain::CoordinateMapBase<SourceFrame, TargetFrame, VolumeDim>&
         map_one,
     const domain::CoordinateMapBase<SourceFrame, TargetFrame, VolumeDim>&
-        map_two) {
+        map_two) noexcept {
   MAKE_GENERATOR(gen);
   std::uniform_real_distribution<> real_dis(-1, 1);
 
@@ -72,7 +72,7 @@ void check_if_maps_are_equal(
 /// \brief Given a coordinate map, check that this map is equal to the identity
 /// by evaluating the map at a random set of points.
 template <typename Map>
-void check_if_map_is_identity(const Map& map) {
+void check_if_map_is_identity(const Map& map) noexcept {
   using IdentityMap = domain::CoordinateMaps::Identity<Map::dim>;
   check_if_maps_are_equal(
       domain::make_coordinate_map<Frame::Inertial, Frame::Grid>(IdentityMap{}),
@@ -87,7 +87,7 @@ void check_if_map_is_identity(const Map& map) {
  */
 template <typename Map>
 void test_jacobian(const Map& map,
-                   const std::array<double, Map::dim>& test_point) {
+                   const std::array<double, Map::dim>& test_point) noexcept {
   // Our default approx value is too stringent for this test
   Approx local_approx = Approx::custom().epsilon(1e-10).scale(1.0);
   const double dx = 1e-4;
@@ -107,8 +107,8 @@ void test_jacobian(const Map& map,
  * multiply together to produce the identity matrix
  */
 template <typename Map>
-void test_inv_jacobian(const Map& map,
-                       const std::array<double, Map::dim>& test_point) {
+void test_inv_jacobian(
+    const Map& map, const std::array<double, Map::dim>& test_point) noexcept {
   const auto jacobian = map.jacobian(test_point);
   const auto inv_jacobian = map.inv_jacobian(test_point);
 
@@ -140,7 +140,7 @@ void test_inv_jacobian(const Map& map,
  * the template parameter to the `CoordinateMap` type.
  */
 template <typename Map, typename... Args>
-void test_coordinate_map_implementation(const Map& map) {
+void test_coordinate_map_implementation(const Map& map) noexcept {
   const auto coord_map =
       domain::make_coordinate_map<Frame::Logical, Frame::Grid>(map);
   MAKE_GENERATOR(gen);
@@ -290,7 +290,7 @@ void test_coordinate_map_argument_types(
  */
 template <typename Map, typename T>
 void test_inverse_map(const Map& map,
-                      const std::array<T, Map::dim>& test_point) {
+                      const std::array<T, Map::dim>& test_point) noexcept {
   CHECK_ITERABLE_APPROX(test_point, map.inverse(map(test_point)).get());
 }
 
@@ -302,7 +302,7 @@ void test_inverse_map(const Map& map,
  * the origin.  The map is expected to be valid on the boundaries of the cube.
  */
 template <typename Map>
-void test_suite_for_map_on_unit_cube(const Map& map) {
+void test_suite_for_map_on_unit_cube(const Map& map) noexcept {
   // Set up random number generator
   MAKE_GENERATOR(gen);
   std::uniform_real_distribution<> real_dis(-1.0, 1.0);
@@ -353,9 +353,9 @@ void test_suite_for_map_on_unit_cube(const Map& map) {
  * This test works only in 3 dimensions.
  */
 template <typename Map>
-void test_suite_for_map_on_sphere(const Map& map,
-                                  const bool include_origin = true,
-                                  const double radius_of_sphere = 1.0) {
+void test_suite_for_map_on_sphere(
+    const Map& map, const bool include_origin = true,
+    const double radius_of_sphere = 1.0) noexcept {
   static_assert(Map::dim == 3, "Works only for a 3d map");
 
   // Set up random number generator
@@ -467,7 +467,7 @@ class OrientationMapIterator {
  * \brief Wedge OrientationMap in each of the six directions used in the
  * Shell and Sphere domain creators.
  */
-inline std::array<OrientationMap<3>, 6> all_wedge_directions() {
+inline std::array<OrientationMap<3>, 6> all_wedge_directions() noexcept {
   const OrientationMap<3> upper_zeta_rotation{};
   const OrientationMap<3> lower_zeta_rotation(std::array<Direction<3>, 3>{
       {Direction<3>::upper_xi(), Direction<3>::lower_eta(),

tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp

@@ -21,6 +21,9 @@ void test(const gsl::not_null<FunctionsOfTime::FunctionOfTime*> f_of_t,
           const gsl::not_null<FunctionsOfTime::PiecewisePolynomial<DerivOrder>*>
               f_of_t_derived,
           double t, const double dt, const double final_time) noexcept {
+  const FunctionsOfTime::PiecewisePolynomial<DerivOrder> f_of_t_derived_copy =
+      *f_of_t_derived;
+  CHECK(*f_of_t_derived == f_of_t_derived_copy);
   while (t < final_time) {
     const auto lambdas0 = f_of_t->func_and_2_derivs(t);
     CHECK(approx(lambdas0[0][0]) == cube(t));
@@ -42,6 +45,7 @@ void test(const gsl::not_null<FunctionsOfTime::FunctionOfTime*> f_of_t,
 
     t += dt;
     f_of_t_derived->update(t, {6.0, 0.0});
+    CHECK(*f_of_t_derived != f_of_t_derived_copy);
   }
   // test time_bounds function
   const auto t_bounds = f_of_t->time_bounds();
@@ -55,9 +59,13 @@ void test_non_const_deriv(
     const gsl::not_null<FunctionsOfTime::PiecewisePolynomial<DerivOrder>*>
         f_of_t_derived,
     double t, const double dt, const double final_time) noexcept {
+  const FunctionsOfTime::PiecewisePolynomial<DerivOrder> f_of_t_derived_copy =
+      *f_of_t_derived;
+  CHECK(*f_of_t_derived == f_of_t_derived_copy);
   while (t < final_time) {
     t += dt;
     f_of_t_derived->update(t, {3.0 + t});
+    CHECK(*f_of_t_derived != f_of_t_derived_copy);
   }
   const auto lambdas0 = f_of_t->func_and_2_derivs(t);
   CHECK(approx(lambdas0[0][0]) == 33.948);

tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstant.cpp

@@ -80,4 +80,52 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.SettleToConstant",
   const std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime> f_of_t3 =
       f_of_t->get_clone();
   test(f_of_t3, match_time, f_t0, dtf_t0, d2tf_t0, A);
+
+  {
+    INFO("Test operator==");
+    CHECK(SettleToConstant{init_func, match_time, decay_time} ==
+          SettleToConstant{init_func, match_time, decay_time});
+    CHECK_FALSE(SettleToConstant{init_func, match_time, decay_time} !=
+                SettleToConstant{init_func, match_time, decay_time});
+
+    CHECK(SettleToConstant{init_func, match_time, decay_time} !=
+          SettleToConstant{init_func, match_time, 2.0 * decay_time});
+    CHECK_FALSE(SettleToConstant{init_func, match_time, decay_time} ==
+                SettleToConstant{init_func, match_time, 2.0 * decay_time});
+
+    CHECK(SettleToConstant{init_func, match_time, decay_time} !=
+          SettleToConstant{init_func, 2.0 * match_time, decay_time});
+    CHECK_FALSE(SettleToConstant{init_func, match_time, decay_time} ==
+                SettleToConstant{init_func, 2.0 * match_time, decay_time});
+
+    CHECK(SettleToConstant{init_func, match_time, decay_time} !=
+          SettleToConstant{{{init_func[0] + 1.0, init_func[1], init_func[2]}},
+                           match_time,
+                           decay_time});
+    CHECK_FALSE(
+        SettleToConstant{init_func, match_time, decay_time} ==
+        SettleToConstant{{{init_func[0] + 1.0, init_func[1], init_func[2]}},
+                         match_time,
+                         decay_time});
+
+    CHECK(SettleToConstant{init_func, match_time, decay_time} !=
+          SettleToConstant{{{init_func[0], init_func[1] + 1.0, init_func[2]}},
+                           match_time,
+                           decay_time});
+    CHECK_FALSE(
+        SettleToConstant{init_func, match_time, decay_time} ==
+        SettleToConstant{{{init_func[0], init_func[1] + 1.0, init_func[2]}},
+                         match_time,
+                         decay_time});
+
+    CHECK(SettleToConstant{init_func, match_time, decay_time} !=
+          SettleToConstant{{{init_func[0], init_func[1], init_func[2] + 1.0}},
+                           match_time,
+                           decay_time});
+    CHECK_FALSE(
+        SettleToConstant{init_func, match_time, decay_time} ==
+        SettleToConstant{{{init_func[0], init_func[1], init_func[2] + 1.0}},
+                         match_time,
+                         decay_time});
+  }
 }

tests/Unit/Domain/Test_Domain.cpp

@@ -133,11 +133,9 @@ void test_1d_domains() {
                              expected_maps);
   }
 }
-}  // namespace
 
-SPECTRE_TEST_CASE("Unit.Domain.Domain", "[Domain][Unit]") { test_1d_domains(); }
-SPECTRE_TEST_CASE("Unit.Domain.Domain.Rectilinear1D1", "[Domain][Unit]") {
-  SECTION("Aligned domain.") {
+void test_1d_rectilinear_domains() {
+  INFO("Aligned domain.") {
     const auto domain = rectilinear_domain<1>(
         Index<1>{3}, std::array<std::vector<double>, 1>{{{0.0, 1.0, 2.0, 3.0}}},
         {}, {}, {{false}}, {}, true);
@@ -160,7 +158,7 @@ SPECTRE_TEST_CASE("Unit.Domain.Domain.Rectilinear1D1", "[Domain][Unit]") {
       CHECK(domain.blocks()[i].neighbors() == expected_block_neighbors[i]);
     }
   }
-  SECTION("Antialigned domain.") {
+  INFO("Antialigned domain.") {
     const OrientationMap<1> aligned{};
     const OrientationMap<1> antialigned{
         std::array<Direction<1>, 1>{{Direction<1>::lower_xi()}}};
@@ -188,7 +186,7 @@ SPECTRE_TEST_CASE("Unit.Domain.Domain.Rectilinear1D1", "[Domain][Unit]") {
   }
 }
 
-SPECTRE_TEST_CASE("Unit.Domain.Domain.Rectilinear2D", "[Domain][Unit]") {
+void test_2d_rectilinear_domains() {
   const OrientationMap<2> half_turn{std::array<Direction<2>, 2>{
       {Direction<2>::lower_xi(), Direction<2>::lower_eta()}}};
   const OrientationMap<2> quarter_turn_cw{std::array<Direction<2>, 2>{
@@ -231,7 +229,7 @@ SPECTRE_TEST_CASE("Unit.Domain.Domain.Rectilinear2D", "[Domain][Unit]") {
   }
 }
 
-SPECTRE_TEST_CASE("Unit.Domain.Domain.Rectilinear3D", "[Domain][Unit]") {
+void test_3d_rectilinear_domains() {
   const OrientationMap<3> aligned{};
   const OrientationMap<3> quarter_turn_cw_xi{std::array<Direction<3>, 3>{
       {Direction<3>::upper_xi(), Direction<3>::upper_zeta(),
@@ -265,7 +263,14 @@ SPECTRE_TEST_CASE("Unit.Domain.Domain.Rectilinear3D", "[Domain][Unit]") {
     CHECK(domain.blocks()[i].neighbors() == expected_block_neighbors[i]);
   }
 }
+}  // namespace
 
+SPECTRE_TEST_CASE("Unit.Domain.Domain", "[Domain][Unit]") {
+  test_1d_domains();
+  test_1d_rectilinear_domains();
+  test_2d_rectilinear_domains();
+  test_3d_rectilinear_domains();
+}
 // [[OutputRegex, Must pass same number of maps as block corner sets]]
 [[noreturn]] SPECTRE_TEST_CASE("Unit.Domain.Domain.BadArgs", "[Domain][Unit]") {
   ASSERTION_TEST();