Add features to history output

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1024]](https://github.com/parthenon-hpc-lab/parthenon/pull/1024) Add features to history output
 - [[PR 852]](https://github.com/parthenon-hpc-lab/parthenon/pull/852) Add Mesh version of UserWorkBeforeOutput
 - [[PR 998]](https://github.com/parthenon-hpc-lab/parthenon/pull/998) tensor indices added to sparse pack
 - [[PR 999]](https://github.com/parthenon-hpc-lab/parthenon/pull/999) Add a post-initialization hook
@@ -15,6 +16,7 @@
 - [[PR 996]](https://github.com/parthenon-hpc-lab/parthenon/pull/996) Remove dynamic allocations from swarm particle creation
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 1024]](https://github.com/parthenon-hpc-lab/parthenon/pull/1024) Add .outN. to history output filenames
 
 ### Fixed (not changing behavior/API/variables/...)
 - [[PR1023]](https://github.com/parthenon-hpc-lab/parthenon/pull/1023) Fix broken param of a scalar bool

example/advection/advection_package.cpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -213,8 +213,15 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       UserHistoryOperation::min, AdvectionHst<Kokkos::Min<Real, HostExecSpace>>,
       "min_advected"));
 
+  // Enroll example vector history output
+  parthenon::HstVec_list hst_vecs = {};
+  hst_vecs.emplace_back(parthenon::HistoryOutputVec(
+      UserHistoryOperation::sum, AdvectionVecHst<Kokkos::Sum<Real, HostExecSpace>>,
+      "advected_powers"));
+
   // add callbacks for HST output identified by the `hist_param_key`
   pkg->AddParam<>(parthenon::hist_param_key, hst_vars);
+  pkg->AddParam<>(parthenon::hist_vec_param_key, hst_vecs);
 
   if (fill_derived) {
     pkg->FillDerivedBlock = SquareIt;
@@ -375,6 +382,45 @@ void PostFill(MeshBlockData<Real> *rc) {
   }
 }
 
+template <typename T>
+std::vector<Real> AdvectionVecHst(MeshData<Real> *md) {
+  auto pmb = md->GetBlockData(0)->GetBlockPointer();
+
+  const auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  const auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  const auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  // Packing variable over MeshBlock as the function is called for MeshData, i.e., a
+  // collection of blocks
+  const auto &advected_pack = md->PackVariables(std::vector<std::string>{"advected"});
+
+  const int nvec = 3;
+
+  std::vector<Real> result(nvec, 0);
+
+  // We choose to apply volume weighting when using the sum reduction.
+  // Downstream this choice will be done on a variable by variable basis and volume
+  // weighting needs to be applied in the reduction region.
+  const bool volume_weighting = std::is_same<T, Kokkos::Sum<Real, HostExecSpace>>::value;
+
+  for (int n = 0; n < nvec; n++) {
+    T reducer(result[n]);
+    parthenon::par_reduce(
+        parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), 0,
+        advected_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lresult) {
+          const auto &coords = advected_pack.GetCoords(b);
+          // `join` is a function of the Kokkos::ReducerConecpt that allows to use the
+          // same call for different reductions
+          const Real vol = volume_weighting ? coords.CellVolume(k, j, i) : 1.0;
+          reducer.join(lresult, pow(advected_pack(b, 0, k, j, i), n + 1) * vol);
+        },
+        reducer);
+  }
+
+  return result;
+}
+
 // Example of how to enroll a history function.
 // Templating is *NOT* required and just implemented here to reuse this function
 // for testing of the UserHistoryOperations curently available in Parthenon (Sum, Min,
@@ -400,9 +446,9 @@ Real AdvectionHst(MeshData<Real> *md) {
   // weighting needs to be applied in the reduction region.
   const bool volume_weighting = std::is_same<T, Kokkos::Sum<Real, HostExecSpace>>::value;
 
-  pmb->par_reduce(
-      PARTHENON_AUTO_LABEL, 0, advected_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s,
-      ib.e,
+  parthenon::par_reduce(
+      parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), 0,
+      advected_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lresult) {
         const auto &coords = advected_pack.GetCoords(b);
         // `join` is a function of the Kokkos::ReducerConecpt that allows to use the same

example/advection/advection_package.hpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -14,6 +14,7 @@
 #define EXAMPLE_ADVECTION_ADVECTION_PACKAGE_HPP_
 
 #include <memory>
+#include <vector>
 
 #include <parthenon/package.hpp>
 
@@ -30,6 +31,8 @@ Real EstimateTimestepBlock(MeshBlockData<Real> *rc);
 TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc);
 template <typename T>
 Real AdvectionHst(MeshData<Real> *md);
+template <typename T>
+std::vector<Real> AdvectionVecHst(MeshData<Real> *md);
 } // namespace advection_package
 
 #endif // EXAMPLE_ADVECTION_ADVECTION_PACKAGE_HPP_

example/advection/parthinput.advection

@@ -3,7 +3,7 @@
 #  Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
 #  Licensed under the 3-clause BSD License, see LICENSE file for details
 # ========================================================================================
-#  (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved.
+#  (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 #
 #  This program was produced under U.S. Government contract 89233218CNA000001 for Los
 #  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -83,6 +83,11 @@ file_type = hst
 dt = 0.05
 
 <parthenon/output4>
+file_type = hst
+dt = 0.1
+packages = advection_app
+
+<parthenon/output5>
 file_type = ascent
 dt = -0.05   # soft disabled by default, as Ascent is an optional dependency
 actions_file = custom_ascent_actions.yaml

src/outputs/history.cpp

@@ -3,7 +3,7 @@
 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
-// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -49,68 +49,113 @@ void HistoryOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm,
   if (signal == SignalHandler::OutputSignal::now) {
     return;
   }
-  std::vector<std::string> all_labels = {};
-  std::vector<Real> all_results = {};
+
+  std::map<UserHistoryOperation, std::vector<Real>> results;
+  std::map<UserHistoryOperation, std::vector<std::string>> labels;
+  std::vector<UserHistoryOperation> ops = {
+      UserHistoryOperation::sum, UserHistoryOperation::max, UserHistoryOperation::min};
+  for (auto &op : ops) {
+    results[op] = {};
+    labels[op] = {};
+  }
+
+  // If "packages" not provided, output history for all packages
+  auto &requested_packages = output_params.packages;
+  Dictionary<std::shared_ptr<StateDescriptor>> packages;
+  if (requested_packages.empty()) {
+    packages = pm->packages.AllPackages();
+  } else {
+    const auto &all_packages = pm->packages.AllPackages();
+    for (const auto &pkg_name : requested_packages) {
+      const auto &it = all_packages.find(pkg_name);
+      if (it != all_packages.end()) {
+        packages[(*it).first] = (*it).second;
+      }
+    }
+  }
 
   // Loop over all packages of the application
-  for (const auto &pkg : pm->packages.AllPackages()) {
-    // Check if the package has enrolled functions which are stored in the
+  for (const auto &pkg : packages) {
+    // Check if the package has enrolled scalar history functions which are stored in the
     // Params under the `hist_param_key` name.
     const auto &params = pkg.second->AllParams();
     if (!params.hasKey(hist_param_key)) {
       continue;
     }
     const auto &hist_vars = params.Get<HstVar_list>(hist_param_key);
 
+    // Get "base" MeshData, which always exists but may not be populated yet
+    auto &md_base = pm->mesh_data.Get();
+    // Populated with all blocks
+    if (md_base->NumBlocks() == 0) {
+      md_base->Set(pm->block_list, pm);
+    } else if (md_base->NumBlocks() != pm->block_list.size()) {
+      PARTHENON_WARN(
+          "Resetting \"base\" MeshData to contain all blocks. This indicates that "
+          "the \"base\" MeshData container has been modified elsewhere. Double check "
+          "that the modification was intentional and is compatible with this reset.")
+      md_base->Set(pm->block_list, pm);
+    }
+
     for (const auto &hist_var : hist_vars) {
-      // Get "base" MeshData, which always exists but may not be populated yet
-      auto &md_base = pm->mesh_data.Get();
-      // Populated with all blocks
-      if (md_base->NumBlocks() == 0) {
-        md_base->Set(pm->block_list, pm);
-      } else if (md_base->NumBlocks() != pm->block_list.size()) {
-        PARTHENON_WARN(
-            "Resetting \"base\" MeshData to contain all blocks. This indicates that "
-            "the \"base\" MeshData container has been modified elsewhere. Double check "
-            "that the modification was intentional and is compatible with this reset.")
-        md_base->Set(pm->block_list, pm);
-      }
       auto result = hist_var.hst_fun(md_base.get());
-#ifdef MPI_PARALLEL
-      // need fence so that the result is ready prior to the MPI call
-      Kokkos::fence();
-      // apply separate chosen operations to each user-defined history output
-      MPI_Op usr_op;
-      switch (hist_var.hst_op) {
-      case UserHistoryOperation::sum:
-        usr_op = MPI_SUM;
-        break;
-      case UserHistoryOperation::max:
-        usr_op = MPI_MAX;
-        break;
-      case UserHistoryOperation::min:
-        usr_op = MPI_MIN;
-        break;
-      }
-      if (Globals::my_rank == 0) {
-        PARTHENON_MPI_CHECK(MPI_Reduce(MPI_IN_PLACE, &result, 1, MPI_PARTHENON_REAL,
-                                       usr_op, 0, MPI_COMM_WORLD));
-      } else {
-        PARTHENON_MPI_CHECK(MPI_Reduce(&result, &result, 1, MPI_PARTHENON_REAL, usr_op, 0,
-                                       MPI_COMM_WORLD));
+      results[hist_var.hst_op].push_back(result);
+      labels[hist_var.hst_op].push_back(hist_var.label);
+    }
+
+    // Check if the package has enrolled vector history functions which are stored in the
+    // Params under the `hist_vec_param_key` name.
+    if (!params.hasKey(hist_vec_param_key)) {
+      continue;
+    }
+    const auto &hist_vecs = params.Get<HstVec_list>(hist_vec_param_key);
+    for (const auto &hist_vec : hist_vecs) {
+      auto result = hist_vec.hst_vec_fun(md_base.get());
+      for (int n = 0; n < result.size(); n++) {
+        std::string label = hist_vec.label + std::to_string(n);
+        results[hist_vec.hst_op].push_back(result[n]);
+        labels[hist_vec.hst_op].push_back(label);
       }
-#endif
+    }
+  }
 
-      all_results.emplace_back(result);
-      all_labels.emplace_back(hist_var.label);
+  // Need fence so result is ready prior to MPI call
+  Kokkos::fence();
+
+#ifdef MPI_PARALLEL
+  for (auto &op : ops) {
+    if (results[op].size() == 0) {
+      continue;
+    }
+
+    MPI_Op usr_op;
+    if (op == UserHistoryOperation::sum) {
+      usr_op = MPI_SUM;
+    } else if (op == UserHistoryOperation::max) {
+      usr_op = MPI_MAX;
+    } else if (op == UserHistoryOperation::min) {
+      usr_op = MPI_MIN;
+    }
+
+    if (Globals::my_rank == 0) {
+      PARTHENON_MPI_CHECK(MPI_Reduce(MPI_IN_PLACE, &(results[op])[0], results[op].size(),
+                                     MPI_PARTHENON_REAL, usr_op, 0, MPI_COMM_WORLD));
+    } else {
+      PARTHENON_MPI_CHECK(MPI_Reduce(&(results[op])[0], &(results[op])[0],
+                                     results[op].size(), MPI_PARTHENON_REAL, usr_op, 0,
+                                     MPI_COMM_WORLD));
     }
   }
+#endif // MPI_PARALLEL
 
-  // only the master rank writes the file
-  // create filename: "file_basename" + ".hst".  There is no file number.
+  // Only the master rank writes the file.
+  // Create filename: "file_basename.out" + [output number] + ".hst". There is no file
+  // number.
   if (Globals::my_rank == 0) {
     std::string fname;
     fname.assign(output_params.file_basename);
+    fname.append(".");
+    fname.append(output_params.file_id);
     fname.append(".hst");
 
     // open file for output
@@ -124,23 +169,26 @@ void HistoryOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm,
 
     // If this is the first output, write header
     if (output_params.file_number == 0) {
-      // NEW_OUTPUT_TYPES:
-
+      // TODO(BRR) optionally overwrite file if this is the first write and not a restart?
       int iout = 1;
       std::fprintf(pfile, "#  History data\n"); // descriptor is first line
       std::fprintf(pfile, "# [%d]=time     ", iout++);
       std::fprintf(pfile, "[%d]=dt       ", iout++);
-      for (auto &lbl : all_labels) {
-        std::fprintf(pfile, "[%d]=%-8s", iout++, lbl.c_str());
+      for (auto &op : ops) {
+        for (auto &label : labels[op]) {
+          std::fprintf(pfile, "[%d]=%-8s", iout++, label.c_str());
+        }
       }
       std::fprintf(pfile, "\n"); // terminate line
     }
 
     // write history variables
     std::fprintf(pfile, output_params.data_format.c_str(), tm->time);
     std::fprintf(pfile, output_params.data_format.c_str(), tm->dt);
-    for (auto &val : all_results) {
-      std::fprintf(pfile, output_params.data_format.c_str(), val);
+    for (auto &op : ops) {
+      for (auto &result : results[op]) {
+        std::fprintf(pfile, output_params.data_format.c_str(), result);
+      }
     }
     std::fprintf(pfile, "\n"); // terminate line
     std::fclose(pfile);

src/outputs/outputs.cpp

@@ -7,7 +7,7 @@
 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
-// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -196,6 +196,11 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) {
         }
       }
 
+      if (op.file_type == "hst") {
+        op.packages = pin->GetOrAddVector<std::string>(pib->block_name, "packages",
+                                                       std::vector<std::string>());
+      }
+
       // set output variable and optional data format string used in formatted writes
       if ((op.file_type != "hst") && (op.file_type != "rst") &&
           (op.file_type != "ascent") && (op.file_type != "histogram")) {
@@ -291,11 +296,10 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) {
     pib = pib->pnext; // move to next input block name
   }
 
-  // check there were no more than one history or restart files requested
-  if (num_hst_outputs > 1 || num_rst_outputs > 1) {
+  // check there were no more than one restart file requested
+  if (num_rst_outputs > 1) {
     msg << "### FATAL ERROR in Outputs constructor" << std::endl
-        << "More than one history or restart output block detected in input file"
-        << std::endl;
+        << "More than one restart output block detected in input file" << std::endl;
     PARTHENON_FAIL(msg);
   }