add history of burgers mass and block counts

benchmarks/burgers/README.md

@@ -74,6 +74,10 @@ To build for execution on a single GPU, it should be sufficient to add the follo
 ```
 where `Kokkos_ARCH` should be set appropriately for the machine (see [here](https://kokkos.github.io/kokkos-core-wiki/keywords.html)).
 
+### Diagnostics
+
+Parthenon-VIBE prints to a history file (default name `burgers.hst`) a time series of the sum of squares of evolved variables integrated over volume for each octant of the domain, as well as the total number of meshblocks in the simulation at that time. To compare these quantities between runs, we provide the `burgers_diff.py` program in the benchmark folder. This will diff two history files and report when the relative difference is greater than some tolerance.
+
 ### Memory Usage
 
 The dominant memory usage in Parthenon-VIBE is for storage of the solution, for which two copies are required to support second order time stepping, for storing the update for a integrator stage (essentially the flux divergence), the intercell fluxes of each variable, for intermediate values of each solution variable on each side of every face, and for a derived quantity that we compute from the evolved solution.  From this we can construct a simple model for the memory usage $M$ as 
@@ -110,4 +114,3 @@ For the GPU, we measure throughput on a single-level mesh ("parthenon/mesh/numle
 
 <p style="text-align:center;"><img src="data/pvibe_gpu_throughput.png" alt="Plot showing throughput on an A100 at different mesh and block sizes" style=width:50%><br />Figure 3: Throughput for different mesh and block sizes on a single 40 GB A100 GPU.</p>
 
-

benchmarks/burgers/burgers.pin

@@ -67,6 +67,11 @@ file_type = hdf5
 dt = -0.4
 variables = U, derived
 
+<parthenon/output1>
+file_type = hst
+data_format = %.14e
+dt = 0.01
+
 <burgers>
 cfl = 0.8
 recon = weno5

benchmarks/burgers/burgers_diff.py

@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+import numpy as np
+from argparse import ArgumentParser
+
+parser = ArgumentParser(
+    prog="burgers_diff.py",
+    description="Compute difference between two history solvers parthenon VIBE",
+)
+parser.add_argument("file1", type=str, help="First file in diff")
+parser.add_argument("file2", type=str, help="Second fiel in diff")
+parser.add_argument(
+    "-t", "--tolerance", type=float, default=1e-8, help="Relative tolerance for diff"
+)
+
+
+def matprint(mat, fmt="14g"):
+    "Print a numpy arraay prettily"
+    col_maxes = [max([len(("{:" + fmt + "}").format(x)) for x in col]) for col in mat.T]
+    for x in mat:
+        for i, y in enumerate(x):
+            print(("{:" + str(col_maxes[i]) + fmt + "}").format(y), end="  ")
+        print("")
+
+
+def get_rel_diff(d1, d2):
+    "Get relative difference between two numpy arrays"
+    return 2 * np.abs(d1 - d2) / (d1 + d2 + 1e-20)
+
+
+if __name__ == "__main__":
+    args = parser.parse_args()
+    d1 = np.loadtxt(args.file1)
+    d2 = np.loadtxt(args.file2)
+    diffs = get_rel_diff(d1, d2)
+    mask = diffs > args.tolerance
+    if np.any(mask):
+        print("Diffs found!")
+        indices = np.transpose(np.nonzero(mask))
+        print("Diff locations (row, column) =", indices)
+        print("Diffs =", diffs[mask])
+    else:
+        print("No diffs found!")

benchmarks/burgers/burgers_package.cpp

@@ -93,6 +93,46 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy});
   pkg->AddField("derived", m);
 
+  // Compute the octants
+  std::vector<Region> octants;
+  std::vector<Real> mesh_mins, mesh_maxs, mesh_mids;
+  for (int d = 1; d <= 3; ++d) {
+    mesh_mins.push_back(pin->GetReal("parthenon/mesh", "x" + std::to_string(d) + "min"));
+    mesh_maxs.push_back(pin->GetReal("parthenon/mesh", "x" + std::to_string(d) + "max"));
+    mesh_mids.push_back(0.5 * (mesh_mins.back() + mesh_maxs.back()));
+  }
+  for (int side1 = 0; side1 < 2; ++side1) {
+    Region r;
+    r.xmin[0] = side1 ? mesh_mids[0] : mesh_mins[0];
+    r.xmax[0] = side1 ? mesh_maxs[0] : mesh_mids[0];
+    for (int side2 = 0; side2 < 2; ++side2) {
+      r.xmin[1] = side2 ? mesh_mids[1] : mesh_mins[1];
+      r.xmax[1] = side2 ? mesh_maxs[1] : mesh_mids[1];
+      for (int side3 = 0; side3 < 2; ++side3) {
+        r.xmin[2] = side3 ? mesh_mids[2] : mesh_mins[2];
+        r.xmax[2] = side3 ? mesh_maxs[2] : mesh_mids[2];
+        octants.push_back(r);
+      }
+    }
+  }
+
+  // Histories
+  auto HstSum = parthenon::UserHistoryOperation::sum;
+  using parthenon::HistoryOutputVar;
+  parthenon::HstVar_list hst_vars = {};
+  int i_octant = 0;
+  for (auto &octant : octants) {
+    auto ReduceMass = [=](MeshData<Real> *md) {
+      return MassHistory(md, octant.xmin[0], octant.xmax[0], octant.xmin[1],
+                         octant.xmax[1], octant.xmin[2], octant.xmax[2]);
+    };
+    hst_vars.emplace_back(HstSum, ReduceMass,
+                          "MS Mass " + std::to_string(i_octant));
+    i_octant++;
+  }
+  hst_vars.emplace_back(HstSum, MeshCountHistory, "Meshblock count");
+  pkg->AddParam(parthenon::hist_param_key, hst_vars);
+
   pkg->EstimateTimestepMesh = EstimateTimestepMesh;
   pkg->FillDerivedMesh = CalculateDerived;
 
@@ -364,4 +404,34 @@ TaskStatus CalculateFluxes(MeshData<Real> *md) {
   return TaskStatus::complete;
 }
 
+Real MassHistory(MeshData<Real> *md, const Real x1min, const Real x1max, const Real x2min,
+                 const Real x2max, const Real x3min, const Real x3max) {
+  const auto ib = md->GetBoundsI(IndexDomain::interior);
+  const auto jb = md->GetBoundsJ(IndexDomain::interior);
+  const auto kb = md->GetBoundsK(IndexDomain::interior);
+
+  std::vector<std::string> vars = {"U"};
+  const auto pack = md->PackVariables(vars);
+
+  Real result = 0.0;
+  parthenon::par_reduce(
+      parthenon::LoopPatternMDRange(), "MassHistory", DevExecSpace(), 0,
+      pack.GetDim(5) - 1, 0, pack.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int v, const int k, const int j, const int i,
+                    Real &lresult) {
+        const auto &coords = pack.GetCoords(b);
+        const Real vol = coords.CellVolume(k, j, i);
+        const Real x1 = coords.Xc<X1DIR>(k, j, i);
+        const Real x2 = coords.Xc<X2DIR>(k, j, i);
+        const Real x3 = coords.Xc<X3DIR>(k, j, i);
+        const Real mask = (x1min <= x1) && (x1 <= x1max) && (x2min <= x2) &&
+                          (x2 <= x2max) && (x3min <= x3) && (x3 <= x3max);
+        lresult += mask * pack(b, v, k, j, i) * pack(b, v, k, j, i) * vol;
+      },
+      Kokkos::Sum<Real>(result));
+  return result;
+}
+
+Real MeshCountHistory(MeshData<Real> *md) { return md->NumBlocks(); }
+
 } // namespace burgers_package

benchmarks/burgers/burgers_package.hpp

@@ -25,6 +25,9 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
 void CalculateDerived(MeshData<Real> *md);
 Real EstimateTimestepMesh(MeshData<Real> *md);
 TaskStatus CalculateFluxes(MeshData<Real> *md);
+Real MassHistory(MeshData<Real> *md, const Real x1min, const Real x1max, const Real x2min,
+                 const Real x2max, const Real x3min, const Real x3max);
+Real MeshCountHistory(MeshData<Real> *md);
 
 // compute the hll flux for Burgers' equation
 KOKKOS_INLINE_FUNCTION
@@ -40,6 +43,10 @@ void lr_to_flux(const Real uxl, const Real uxr, const Real uyl, const Real uyr,
   fuz = 0.5 * (sr * uzl * upl - sl * uzr * upr + sl * sr * (uzr - uzl)) * islsr;
 }
 
+struct Region {
+  std::array<Real, 3> xmin, xmax;
+};
+
 } // namespace burgers_package
 
 #endif // BENCHMARKS_BURGERS_BURGERS_PACKAGE_HPP_