Add kerr evolution test mc

docs/DevGuide/OptionParsing.md

@@ -62,8 +62,8 @@ input file as usual, but a warning will be issued if the specified
 value does not match the suggestion.
 
 Examples:
-\snippet Test_Options.cpp options_example_scalar_struct
-\snippet Test_Options.cpp options_example_vector_struct
+\snippet Options/Test_Options.cpp options_example_scalar_struct
+\snippet Options/Test_Options.cpp options_example_vector_struct
 
 The option type can be any type understood natively by yaml-cpp
 (fundamentals, `std::string`, and `std::map`, `std::vector`,
@@ -78,7 +78,7 @@ struct. The alias should refer to a type that, like option tags, defines a help
 string and may override a static `name()` function.
 
 Example:
-\snippet Test_Options.cpp options_example_group
+\snippet Options/Test_Options.cpp options_example_group
 
 ## Constructible classes
 
@@ -105,15 +105,15 @@ options ensuring that the error message will have a full backtrace so it is easy
 for the user to diagnose.
 
 Example:
-\snippet Test_CustomTypeConstruction.cpp class_creation_example
+\snippet Options/Test_CustomTypeConstruction.cpp class_creation_example
 
 Classes may use the Metavariables struct, which is effectively the compile time
 input file, in their parsing by templating the `options` type alias or by taking
 the Metavariables as a final argument to the constructor (after the
 `Options::Context`).
 
 Example:
-\snippet Test_CustomTypeConstruction.cpp class_creation_example_with_metavariables
+\snippet Options/Test_CustomTypeConstruction.cpp class_creation_example_with_metavariables
 
 ## Factory
 
@@ -123,7 +123,7 @@ from `Base`.  The list of creatable derived classes is specified in
 the `factory_creation` struct in the metavariables, which must contain
 a `factory_classes` type alias that is a `tmpl::map` from base classes
 to lists of derived classes:
-\snippet Test_Factory.cpp factory_creation
+\snippet Options/Test_Factory.cpp factory_creation
 
 When a `std::unique_ptr<Base>` is requested, the factory will expect a
 single YAML argument specifying the name of the class (as given by a
@@ -133,8 +133,8 @@ string, otherwise it must be given as a single key-value pair, with
 the key the name of the class.  The value portion of this pair is then
 used to create the requested derived class in the same way as an
 explicitly constructible class.  Examples:
-\snippet Test_Factory.cpp factory_without_arguments
-\snippet Test_Factory.cpp factory_with_arguments
+\snippet Options/Test_Factory.cpp factory_without_arguments
+\snippet Options/Test_Factory.cpp factory_with_arguments
 
 \anchor custom-parsing
 ## Custom parsing
@@ -154,7 +154,7 @@ The create function can perform any operations required to construct
 the object.
 
 Example of using a specialization to parse an enum:
-\snippet Test_CustomTypeConstruction.cpp enum_creation_example
+\snippet Options/Test_CustomTypeConstruction.cpp enum_creation_example
 
 Note that in the case where the `create` function does *not* need to use the
 `Metavariables` it is recommended that a general implementation forward to an
@@ -167,8 +167,8 @@ metavariables-independent case efficiently. As a concrete example, the general
 definition and forward declaration of the `void` specialization in the header
 file would be:
 
-\snippet Test_CustomTypeConstruction.cpp enum_void_creation_header_example
+\snippet Options/Test_CustomTypeConstruction.cpp enum_void_creation_header_example
 
 while in the `cpp` file the definition of the `void` specialization is:
 
-\snippet Test_CustomTypeConstruction.cpp enum_void_creation_cpp_example
+\snippet Options/Test_CustomTypeConstruction.cpp enum_void_creation_cpp_example

src/Evolution/DiscontinuousGalerkin/AtomicInboxBoundaryData.cpp

@@ -3,15 +3,39 @@
 
 #include "Evolution/DiscontinuousGalerkin/AtomicInboxBoundaryData.hpp"
 
+#include <atomic>
 #include <cstddef>
 
 #include "Domain/Structure/Direction.hpp"
 #include "Domain/Structure/DirectionalId.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Utilities/ConstantExpressions.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/System/Abort.hpp"
 
 namespace evolution::dg {
+template <size_t Dim>
+AtomicInboxBoundaryData<Dim>::AtomicInboxBoundaryData(
+    AtomicInboxBoundaryData<Dim>&& rhs) noexcept {
+  if (rhs.message_count.load(std::memory_order_acquire) != 0) {
+    sys::abort(
+        "You cannot move an AtomicInboxBoundaryData with non-zero message "
+        "count.");
+  }
+  for (size_t i = 0; i < rhs.boundary_data_in_directions.size(); ++i) {
+    // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
+    if (not rhs.boundary_data_in_directions[i].empty()) {
+      sys::abort(
+          "You cannot move an AtomicInboxBoundaryData with data in "
+          "boundary_data_in_directions.");
+    }
+  }
+  message_count.store(0, std::memory_order_release);
+  number_of_neighbors.store(
+      rhs.number_of_neighbors.load(std::memory_order_acquire),
+      std::memory_order_release);
+}
+
 template <size_t Dim>
 size_t AtomicInboxBoundaryData<Dim>::index(
     const DirectionalId<Dim>& neighbor_directional_id) {
@@ -47,6 +71,35 @@ size_t AtomicInboxBoundaryData<Dim>::index(
   }
 }
 
+template <size_t Dim>
+void AtomicInboxBoundaryData<Dim>::pup(PUP::er& p) {
+  if (UNLIKELY(number_of_neighbors.load(std::memory_order_acquire) != 0)) {
+    ERROR(
+        "Can only serialize AtomicInboxBoundaryData if there are no messages. "
+        "We need to be very careful about serializing atomics since "
+        "serialization requires strong synchronization like a lock.");
+  }
+  for (size_t i = 0; i < boundary_data_in_directions.size(); ++i) {
+    if (UNLIKELY(not gsl::at(boundary_data_in_directions, i).empty())) {
+      ERROR(
+          "We can only serialize empty StaticSpscQueues but the queue in "
+          "element "
+          << i << " is not empty.");
+    }
+  }
+  if (p.isUnpacking()) {
+    std::atomic_uint::value_type number_of_neighbors_to_serialize = 0;
+    p | number_of_neighbors_to_serialize;
+    number_of_neighbors.store(number_of_neighbors_to_serialize,
+                              std::memory_order_release);
+    message_count.store(0, std::memory_order_release);
+  } else {
+    std::atomic_uint::value_type number_of_neighbors_to_serialize =
+        number_of_neighbors.load(std::memory_order_acquire);
+    p | number_of_neighbors_to_serialize;
+  }
+}
+
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
 
 #define INSTANTIATION(r, data) \

src/Evolution/DiscontinuousGalerkin/AtomicInboxBoundaryData.hpp

@@ -33,10 +33,22 @@ namespace evolution::dg {
  * messages. Note that some additional logic is needed also for supporting
  * local time stepping, since not every message entry "counts" since it
  * depends on the time level of neighboring elements.
+ *
+ * \warning Only `AtomicInboxBoundaryData` with zero messages can be move
+ * constructed. A non-zero number of neighbors is allowed. This is necessary
+ * in order to be able to serialize a
+ * `std::unordered_map<Key,AtomicInboxBoundaryData>`.
  */
 template <size_t Dim>
 struct AtomicInboxBoundaryData {
   using stored_type = evolution::dg::BoundaryData<Dim>;
+  AtomicInboxBoundaryData() = default;
+  AtomicInboxBoundaryData(const AtomicInboxBoundaryData&) = delete;
+  AtomicInboxBoundaryData& operator=(const AtomicInboxBoundaryData&) = delete;
+  AtomicInboxBoundaryData(AtomicInboxBoundaryData&& rhs) noexcept;
+  AtomicInboxBoundaryData& operator=(AtomicInboxBoundaryData&&) noexcept =
+      delete;
+  ~AtomicInboxBoundaryData() = default;
 
   /*!
    * Computes the 1d index into the `boundary_data_in_directions` array
@@ -62,7 +74,9 @@ struct AtomicInboxBoundaryData {
    */
   static size_t index(const DirectionalId<Dim>& directional_id);
 
-  // We use 20 entiries in the SPSC under the assumption that each neighbor
+  void pup(PUP::er& p);
+
+  // We use 20 entries in the SPSC under the assumption that each neighbor
   // will never insert more than 20 entries before the element uses
   // them. While in practice a smaller buffer could be used, this is to
   // safeguard against future features.

src/Evolution/Particles/MonteCarlo/EvolvePackets.cpp

@@ -5,6 +5,7 @@
 
 #include "Evolution/Particles/MonteCarlo/Packet.hpp"
 #include "Evolution/Particles/MonteCarlo/Scattering.hpp"
+#include "Parallel/Printf/Printf.hpp"
 
 namespace Particles::MonteCarlo {
 
@@ -61,6 +62,14 @@ void evolve_single_packet_on_geodesic(
   // Calculate time derivative of 3-momentum at beginning of time step
   detail::time_derivative_momentum_geodesic(&dpdt, *packet, lapse, d_lapse,
                                             d_shift, d_inv_spatial_metric);
+  // Parallel::printf("Initial coord: %.5f %.5f %.5f\n", packet->coordinates[0],
+  //                  packet->coordinates[1], packet->coordinates[2]);
+  // Parallel::printf("Initial momentum: %.5f %.5f %.5f\n", packet->momentum[0],
+  //                  packet->momentum[1], packet->momentum[2]);
+  // Parallel::printf("dpdt: %.5f %.5f %.5f\n", gsl::at(dpdt, 0), gsl::at(dpdt,
+  // 1),
+  //                  gsl::at(dpdt, 2));
+
   // Take half-step (momentum only, as time derivative is independent of
   // position)
   for (size_t i = 0; i < 3; i++) {
@@ -74,6 +83,12 @@ void evolve_single_packet_on_geodesic(
   // Calculate time derivative of 3-momentum and position at half-step
   detail::time_derivative_momentum_geodesic(&dpdt, *packet, lapse, d_lapse,
                                             d_shift, d_inv_spatial_metric);
+  // Parallel::printf("Half step momentum: %.5f %.5f %.5f\n",
+  // packet->momentum[0],
+  //                 packet->momentum[1], packet->momentum[2]);
+  // Parallel::printf("dpdt: %.5f %.5f %.5f\n", gsl::at(dpdt, 0), gsl::at(dpdt,
+  // 1),
+  //                 gsl::at(dpdt, 2));
   for (size_t i = 0; i < 3; i++) {
     gsl::at(dxdt_inertial, i) = (-1.0) * shift.get(i)[closest_point_index];
     if (mesh_velocity.has_value()) {
@@ -86,19 +101,29 @@ void evolve_single_packet_on_geodesic(
           packet->momentum[j] / packet->momentum_upper_t;
     }
   }
+  // Parallel::printf("dxdt_inert: %.5f %.5f %.5f\n", gsl::at(dxdt_inertial, 0),
+  //                 gsl::at(dxdt_inertial, 1), gsl::at(dxdt_inertial, 2));
+
   for (size_t i = 0; i < 3; i++) {
     for (size_t j = 0; j < 3; j++) {
       gsl::at(dxdt_logical, i) +=
           gsl::at(dxdt_inertial, j) *
           inverse_jacobian_logical_to_inertial.get(i, j)[closest_point_index];
     }
   }
+  // Parallel::printf("dxdt_inert: %.5f %.5f %.5f\n", gsl::at(dxdt_logical, 0),
+  //                 gsl::at(dxdt_logical, 1), gsl::at(dxdt_logical, 2));
+
   // Take full time step
   for (size_t i = 0; i < 3; i++) {
     packet->momentum[i] = gsl::at(p0, i) + gsl::at(dpdt, i) * time_step;
     packet->coordinates[i] += gsl::at(dxdt_logical, i) * time_step;
   }
   packet->time += time_step;
+  // Parallel::printf("Final coord: %.5f %.5f %.5f\n", packet->coordinates[0],
+  //                 packet->coordinates[1], packet->coordinates[2]);
+  // Parallel::printf("Final momentum: %.5f %.5f %.5f\n", packet->momentum[0],
+  //                  packet->momentum[1], packet->momentum[2]);
 }
 
 }  // namespace Particles::MonteCarlo

src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp

@@ -144,7 +144,7 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
   double dt_min = -1.0;
   double initial_time = -1.0;
   // Loop over packets
-  const size_t n_packets = packets->size();
+  size_t n_packets = packets->size();
   for (size_t p = 0; p < n_packets; p++) {
     Packet& packet = (*packets)[p];
 
@@ -224,6 +224,7 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
         std::swap((*packets)[p], (*packets)[n_packets - 1]);
         packets->pop_back();
         p--;
+        n_packets--;
         break;
       }
       // If the next event was a scatter, perform that scatter and
@@ -293,6 +294,7 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
           (*packets)[p] = (*packets)[n_packets - 1];
           packets->pop_back();
           p--;
+          n_packets--;
           break;
         }
       }

src/Options/Protocols/FactoryCreation.hpp

@@ -22,7 +22,7 @@ namespace Options::protocols {
  *
  * Here's an example for a class conforming to this protocol:
  *
- * \snippet Test_Factory.cpp factory_creation
+ * \snippet Options/Test_Factory.cpp factory_creation
  */
 struct FactoryCreation {
   template <typename ConformingType>

src/Time/TimeSteppers/Factory.hpp

@@ -21,21 +21,19 @@
 namespace TimeSteppers {
 /// Typelist of available TimeSteppers
 using time_steppers =
-    tmpl::list<TimeSteppers::AdamsBashforth,
-               TimeSteppers::AdamsMoultonPc<false>,
-               TimeSteppers::AdamsMoultonPc<true>,
-               TimeSteppers::ClassicalRungeKutta4, TimeSteppers::DormandPrince5,
-               TimeSteppers::Heun2, TimeSteppers::Rk3HesthavenSsp,
-               TimeSteppers::Rk3Kennedy, TimeSteppers::Rk3Owren,
-               TimeSteppers::Rk3Pareschi, TimeSteppers::Rk4Kennedy,
-               TimeSteppers::Rk4Owren, TimeSteppers::Rk5Owren,
-               TimeSteppers::Rk5Tsitouras>;
+    tmpl::list<AdamsBashforth, AdamsMoultonPc<false>, AdamsMoultonPc<true>,
+               ClassicalRungeKutta4, DormandPrince5, Heun2, Rk3HesthavenSsp,
+               Rk3Kennedy, Rk3Owren, Rk3Pareschi, Rk4Kennedy, Rk4Owren,
+               Rk5Owren, Rk5Tsitouras>;
 
 /// Typelist of available LtsTimeSteppers
-using lts_time_steppers = tmpl::list<TimeSteppers::AdamsBashforth>;
+using lts_time_steppers = tmpl::list<AdamsBashforth>;
 
 /// Typelist of available ImexTimeSteppers
 using imex_time_steppers =
-    tmpl::list<TimeSteppers::Heun2, TimeSteppers::Rk3Kennedy,
-               TimeSteppers::Rk3Pareschi, TimeSteppers::Rk4Kennedy>;
+    tmpl::list<Heun2, Rk3Kennedy, Rk3Pareschi, Rk4Kennedy>;
+
+/// Typelist of TimeSteppers whose substep times are strictly increasing
+using increasing_substep_time_steppers =
+    tmpl::list<AdamsBashforth, Rk3Owren, Rk4Owren>;
 }  // namespace TimeSteppers

support/Pipelines/Bbh/__init__.py

@@ -3,22 +3,48 @@
 
 import click
 
-from .FindHorizon import find_horizon_command
-from .InitialData import generate_id_command
-from .Inspiral import start_inspiral_command
-from .Ringdown import start_ringdown_command
 
+# Load subcommands lazily, i.e., only import the module when the subcommand is
+# invoked. This is important so the CLI responds quickly.
+class Bbh(click.MultiCommand):
+    def list_commands(self, ctx):
+        return [
+            "find-horizon",
+            "generate-id",
+            "start-inspiral",
+            "start-ringdown",
+        ]
 
-@click.group(name="bbh")
+    def get_command(self, ctx, name):
+        if name == "find-horizon":
+            from .FindHorizon import find_horizon_command
+
+            return find_horizon_command
+        if name == "generate-id":
+            from .InitialData import generate_id_command
+
+            return generate_id_command
+        elif name == "start-inspiral":
+            from .Inspiral import start_inspiral_command
+
+            return start_inspiral_command
+        elif name == "start-ringdown":
+            from .Ringdown import start_ringdown_command
+
+            return start_ringdown_command
+
+        available_commands = " " + "\n ".join(self.list_commands(ctx))
+        raise click.UsageError(
+            f"The command '{name}' is not implemented. "
+            f"Available commands are:\n{available_commands}"
+        )
+
+
+@click.group(name="bbh", cls=Bbh)
 def bbh_pipeline():
     """Pipeline for binary black hole simulations."""
     pass
 
 
-bbh_pipeline.add_command(generate_id_command)
-bbh_pipeline.add_command(find_horizon_command)
-bbh_pipeline.add_command(start_inspiral_command)
-bbh_pipeline.add_command(start_ringdown_command)
-
 if __name__ == "__main__":
     bbh_pipeline(help_option_names=["-h", "--help"])