Add support for "analysis" restart

doc/sphinx/src/boundary_conditions.rst

@@ -101,8 +101,8 @@ you like. Reference implementations of the standard boundary conditions
 are available `here <https://github.com/parthenon-hpc-lab/parthenon/blob/develop/src/bvals/boundary_conditions.cpp>`__.
 
 
-Per package user-defined boundary conditions.
----------------------------------
+Per package user-defined boundary conditions
+--------------------------------------------
 
 In addition to user defined *global* boundary conditions, Parthenon also supports 
 registration of boundary conditions at the *per package* level. These per package 

doc/sphinx/src/outputs.rst

@@ -143,6 +143,28 @@ selected as all the variables that have either the ``Independent`` or
 ``Restart`` ``Metadata`` flags specified. No other intervention is
 required by the developer.
 
+Postprocessing/native analysis
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+A rudimentary postprocessing option is available to analyze data in
+restart files within the downstream/Parthenon framework, i.e., making use of
+native compute capabilities.
+
+To trigger this kind of analysis, restart with ``-a restart.rhdf`` (instead of
+``-r``).
+This will launch the standard driver as if a simulation would be restarted,
+but never call the main time loop.
+Only the callbacks ``UserWorkBeforeLoop`` are executed.
+Afterwards, all output blocks that have the parameter ``analysis_output=true`` are
+going to be processed including ``UserWorkBeforeOutput`` callbacks.
+
+Note, the standard modifications to the original input parameters via the command line
+or via an input file apply.
+A typical use case is, for example, to calculate histograms a posteriori, i.e., a new
+output block is specified in a file called ``sample_hist.in`` with all necessary details
+(particularly the ``analysis_output=true`` parameter within said block) and then the
+data can be processed via ``-a restart.rhdf -i sample_hist.in``.
+
 .. _output hist files:
 
 History Files

doc/sphinx/src/tasks.rst

@@ -5,9 +5,9 @@ Tasks
 
 Parthenon's tasking infrastructure is how downstream applications describe 
 and execute their work.  Tasks are organized into a hierarchy of objects.
-``TaskCollection``s have one or more ``TaskRegion``s, ``TaskRegion``s have
-one or more ``TaskList``s, and ``TaskList``s can have one or more sublists
-(that are themselves ``TaskList``s).
+``TaskCollection``\s have one or more ``TaskRegion``\s, ``TaskRegion``\s have
+one or more ``TaskList``\s, and ``TaskList``\s can have one or more sublists
+(that are themselves ``TaskList``\s).
 
 Task
 ----
@@ -18,7 +18,7 @@ that stores the necessary data to invoke a downstream code's functions with
 the desired arguments.  Importantly, however, it also stores information that
 relates itself to other tasks, namely the tasks that must be complete before
 it should execute and the tasks that may be available to run after it completes.
-In other words, ``Task``s are nodes in a directed (possibly cyclic) graph, and
+In other words, ``Task``\s are nodes in a directed (possibly cyclic) graph, and
 include the edges that connect to it and emerge from it.
 
 TaskList
@@ -52,10 +52,13 @@ and maximum number of times the sublist should execute.  Passing something like
 leading to a sublist that never cycles.  ``AddSublist``
 returns a ``std::pair<TaskList&, TaskID>`` which is conveniently accessed via
 a structured binding, e.g.
+
 .. code:: cpp
+
   TaskID none;
   auto [child_list, child_list_id] = parent_list.AddSublist(dependencies, {1,3});
   auto task_id = child_list.AddTask(none, SomeFunction, arg1, arg2);
+
 In the above example, passing ``none`` as the dependency for the task added to
 ``child_list`` does not imply that this task can execute at any time since
 ``child_list`` itself has dependencies that must be satisfied before any of its
@@ -65,22 +68,22 @@ TaskRegion
 ----------
 
 Under the hood, a ``TaskRegion`` is a directed, possibly cyclic graph.  The graph
-is built up incrementally as tasks are added to the ``TaskList``s within the 
-``TaskRegion``, and it's construction is completed upon the first time it's
-executed.  ``TaskRegion``s can have one or more ``TaskList``s.  The primary reason
+is built up incrementally as tasks are added to the ``TaskList``\s within the
+``TaskRegion``\, and it's construction is completed upon the first time it's
+executed.  ``TaskRegion``\s can have one or more ``TaskList``\s.  The primary reason
 for this is to allow flexibility in how work is broken up into tasks (and
 eventually kernels).  A region with many lists will produce many small
 tasks/kernels, but may expose more asynchrony (e.g. MPI communication).  A region
 with fewer lists will produce more work per kernel (which may be good for GPUs,
 for example), but may limit asynchrony.  Typically, each list is tied to a unique
 partition of the mesh blocks owned by a rank.  ``TaskRegion`` only provides a few
 public facing functions:
-- ``TaskListStatus Execute(ThreadPool &pool)``: ``TaskRegion``s can be executed, requiring a
+- ``TaskListStatus Execute(ThreadPool &pool)``\: ``TaskRegion``\s can be executed, requiring a
 ``ThreadPool`` be provided by the caller.  In practice, ``Execute`` is usually
 called from the ``Execute`` member function of ``TaskCollection``.
-- ``TaskList& operator[](const int i)``: return a reference to the ``i``th
+- ``TaskList& operator[](const int i)``\: return a reference to the ``i``\th
 ``TaskList`` in the region.
-- ``size_t size()``: return the number of ``TaskList``s in the region.
+- ``size_t size()``\: return the number of ``TaskList``\s in the region.
 
 TaskCollection
 --------------

src/argument_parser.hpp

@@ -53,6 +53,12 @@ class ArgParse {
           res_flag = 1;
           restart_filename = argv[++i];
           break;
+        case 'a': // -a <restart_file>
+          invalid = invalid_arg();
+          res_flag = 1;
+          analysis_flag = true;
+          restart_filename = argv[++i];
+          break;
         case 'd': // -d <run_directory>
           invalid = invalid_arg();
           prundir = argv[++i];
@@ -80,6 +86,7 @@ class ArgParse {
             std::cout << "Options:" << std::endl;
             std::cout << "  -i <file>       specify input file [athinput]\n";
             std::cout << "  -r <file>       restart with this file\n";
+            std::cout << "  -a <file>       analyze/postprocess this file\n";
             std::cout << "  -d <directory>  specify run dir [current dir]\n";
             std::cout << "  -n              parse input file and quit\n";
             std::cout << "  -c              show configuration and quit\n";
@@ -118,6 +125,7 @@ class ArgParse {
   char *input_filename = nullptr;
   char *restart_filename = nullptr;
   char *prundir = nullptr;
+  bool analysis_flag = false;
   int res_flag = 0;
   int narg_flag = 0;
   int mesh_flag = 0;

src/driver/driver.cpp

@@ -79,7 +79,9 @@ DriverStatus EvolutionDriver::Execute() {
     }
   } // UserWorkBeforeLoop
 
-  OutputSignal signal = OutputSignal::none;
+  OutputSignal signal = pinput->GetBoolean("parthenon/job", "run_only_analysis")
+                            ? OutputSignal::analysis
+                            : OutputSignal::none;
   pouts->MakeOutputs(pmesh, pinput, &tm, signal);
   pmesh->mbcnt = 0;
   int perf_cycle_offset =
@@ -91,7 +93,7 @@ DriverStatus EvolutionDriver::Execute() {
 
   { // Main t < tmax loop region
     PARTHENON_INSTRUMENT
-    while (tm.KeepGoing()) {
+    while (tm.KeepGoing() && signal != OutputSignal::analysis) {
       if (Globals::my_rank == 0) OutputCycleDiagnostics();
 
       pmesh->PreStepUserWorkInLoop(pmesh, pinput, tm);
@@ -141,8 +143,11 @@ DriverStatus EvolutionDriver::Execute() {
   pmesh->UserWorkAfterLoop(pmesh, pinput, tm);
 
   DriverStatus status = DriverStatus::complete;
-
-  pouts->MakeOutputs(pmesh, pinput, &tm, OutputSignal::final);
+  // Do *not* write the "final" output, if this is analysis run.
+  // The analysis output itself has already been written above before the main loop.
+  if (signal != OutputSignal::analysis) {
+    pouts->MakeOutputs(pmesh, pinput, &tm, OutputSignal::final);
+  }
   PostExecute(status);
   return status;
 }

src/outputs/outputs.cpp

@@ -76,6 +76,7 @@
 #include "parameter_input.hpp"
 #include "parthenon_arrays.hpp"
 #include "utils/error_checking.hpp"
+#include "utils/utils.hpp"
 
 namespace parthenon {
 
@@ -146,6 +147,8 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) {
       // read cartesian mapping option
       op.cartesian_vector = false;
 
+      op.analysis_flag = pin->GetOrAddBoolean(op.block_name, "analysis_output", false);
+
       // read single precision output option
       const bool is_hdf5_output = (op.file_type == "rst") || (op.file_type == "hdf5");
 
@@ -425,7 +428,10 @@ void Outputs::MakeOutputs(Mesh *pm, ParameterInput *pin, SimTime *tm,
     if ((tm == nullptr) ||
         ((ptype->output_params.dt >= 0.0) &&
          ((tm->ncycle == 0) || (tm->time >= ptype->output_params.next_time) ||
-          (tm->time >= tm->tlim) || (signal != SignalHandler::OutputSignal::none)))) {
+          (tm->time >= tm->tlim) || (signal == SignalHandler::OutputSignal::now) ||
+          (signal == SignalHandler::OutputSignal::final) ||
+          (signal == SignalHandler::OutputSignal::analysis &&
+           ptype->output_params.analysis_flag)))) {
       if (first && ptype->output_params.file_type != "hst") {
         pm->ApplyUserWorkBeforeOutput(pm, pin, *tm);
         first = false;

src/outputs/outputs.hpp

@@ -51,6 +51,7 @@ struct OutputParameters {
   std::string file_basename;
   int file_number_width;
   bool file_label_final;
+  bool analysis_flag; // write this output for analysis/postprocessing restarts
   std::string file_id;
   std::vector<std::string> variables;
   std::vector<std::string> component_labels;

src/parthenon_manager.cpp

@@ -113,7 +113,7 @@ ParthenonStatus ParthenonManager::ParthenonInitEnv(int argc, char *argv[]) {
     std::istringstream is(inputString);
     pinput->LoadFromStream(is);
   }
-  // If an input was provided
+  // If an input file was provided
   if (arg.input_filename != nullptr) {
     // Modify info read from restart file
     if (arg.res_flag != 0) {
@@ -130,6 +130,16 @@ ParthenonStatus ParthenonManager::ParthenonInitEnv(int argc, char *argv[]) {
 
   // Modify based on command line inputs
   pinput->ModifyFromCmdline(argc, argv);
+
+  PARTHENON_REQUIRE_THROWS(
+      !pinput->DoesParameterExist("parthenon/job", "run_only_analysis") ||
+          pinput->GetBoolean("parthenon/job", "run_only_analysis") == false,
+      "'parthenon/job/run_only_analysis=true' input parameter was found indicating "
+      "manual modification or restarting from an output written during analysis, which "
+      "is undefined behavior. If you don't know how this was triggered, please contact "
+      "the Parthenon developers.");
+  pinput->SetBoolean("parthenon/job", "run_only_analysis", arg.analysis_flag);
+
   // Set the global number of ghost zones
   Globals::nghost = pinput->GetOrAddInteger("parthenon/mesh", "nghost", 2);
 

src/utils/utils.hpp

@@ -80,7 +80,7 @@ constexpr auto get_array_from_tuple(tuple_t &&tuple) {
 //  \brief static data and functions that implement a simple signal handling system
 namespace SignalHandler {
 
-enum class OutputSignal { none, now, final };
+enum class OutputSignal { none, now, final, analysis };
 constexpr int nsignal = 3;
 // using the +1 for signaling based on "output_now" trigger
 static volatile int signalflag[nsignal + 1];