Add features to history output

CHANGELOG.md

@@ -6,9 +6,11 @@
 
 
 ### 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/...)
+- [[PR 1024]](https://github.com/parthenon-hpc-lab/parthenon/pull/1024) Add features to history output
 - [[PR 1031]](https://github.com/parthenon-hpc-lab/parthenon/pull/1031) Fix bug in non-cell centered AMR
 
 ### Infrastructure (changes irrelevant to downstream codes)

doc/sphinx/src/interface/state.rst

@@ -159,7 +159,7 @@ functions that are all called at the interval according to the input
 parameters, see :ref:`output documention <output hist files>`.
 
 To enroll functions create a list of callback function with the
-appropriate reduction operation:
+appropriate reduction operations for either scalar data:
 
 .. code:: cpp
 
@@ -172,6 +172,19 @@ appropriate reduction operation:
    // add callbacks for HST output identified by the `hist_param_key`
    pkg->AddParam<>(parthenon::hist_param_key, hst_vars);
 
+or vector data:
+
+.. code:: cpp
+
+   // List (vector) of HistoryOutputVar that will all be enrolled as output variables
+   parthenon::HstVec_list hst_vecs = {};
+
+   // Add a callback function
+   hst_vecs.emplace_back(parthenon::HistoryOutputVec(UserHistoryOperation::sum, MyHstVecFunction, "my vector label"));
+
+   // add callbacks for HST output identified by the `hist_vec_param_key`
+   pkg->AddParam<>(parthenon::hist_vec_param_key, hst_vecs);
+
 Here, ``HistoryOutputVar`` is a ``struct`` containing the global (over
 all blocks of all ranks) reduction operation, ``MyHstFunction`` is a
 callback function (see below), and ``"my label"`` is the string to be
@@ -194,6 +207,7 @@ Callback functions need to have the following signature
 .. code:: cpp
 
    Real MyHstFunction(MeshData<Real> *md);
+   std::vector<Real> MyHstVecFunction(MeshData<Real> *md);
 
 i.e., they will always work on ``MeshData``. *Note*, currently history
 output will always be calculated for the “base” container. More

doc/sphinx/src/outputs.rst

@@ -149,18 +149,22 @@ History Files
 
 In the input file, include a ``<parthenon/output*>`` block and specify
 ``file_type = hst``. A ``dt`` parameter controls the frequency of
-outputs for simulations involving evolution. A ``<parthenon/output*>``
+outputs for simulations involving evolution. The default behavior is to provide
+all enrolled history outputs, but output can be limited to a specific set of
+packages with an optional comma-separated list argument
+``packages =  package_a, package_b``. A ``<parthenon/output*>``
 block might look like
 
 ::
 
    <parthenon/output8>
    file_type = hst
    dt = 1.0
+   packages = advection_app
 
 This will produce a text file (``.hst``) output file every 1 units of
 simulation time. The content of the file is determined by the functions
-enrolled by a specific package, see :ref:`state history output`.
+enrolled by specific packages, see :ref:`state history output`.
 
 Histograms
 ----------
@@ -171,7 +175,7 @@ Currently supported are
 
 - 1D and 2D histograms (see examples below)
 - binning by variable or coordinate (x1, x2, x3 and radial distance)
-- counting samples and or summing a variable  
+- counting samples and or summing a variable
 - weighting by volume and/or variable
 
 The output format follows ``numpy`` convention, so that plotting data
@@ -322,7 +326,7 @@ The following is a minimal example to plot a 1D and 2D histogram from the output
       y = infile["other_name/y_edges"][:]
       z = infile["other_name/data"][:].T   # note the transpose here (so that the data matches the axis for the pcolormesh)
       plt.pcolormesh(x,y,z,)
-      plt.show()   
+      plt.show()
 
 Ascent (optional)
 -----------------

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

example/sparse_advection/sparse_advection_driver.cpp

@@ -122,12 +122,12 @@ TaskCollection SparseAdvectionDriver::MakeTaskCollection(BlockList_t &blocks,
                              mdudt.get(), beta * dt, mc1.get());
 
     // do boundary exchange
-    auto restrict =
+    auto prolong_restrict =
         parthenon::AddBoundaryExchangeTasks(update, tl, mc1, pmesh->multilevel);
 
     // if this is the last stage, check if we can deallocate any sparse variables
     if (stage == integrator->nstages) {
-      tl.AddTask(restrict, SparseDealloc, mc1.get());
+      tl.AddTask(prolong_restrict, SparseDealloc, mc1.get());
     }
   }