[WIP]: Cell Count Example

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 734]](https://github.com/lanl/parthenon/pull/734) Count cells example
 - [[PR 729]](https://github.com/lanl/parthenon/pull/729) Optional modifications to output format
 - [[PR 717]](https://github.com/lanl/parthenon/pull/717) Add ghost zone plotting capability to phdf.py and movie2d.py
 - [[PR 712]](https://github.com/lanl/parthenon/pull/712) Allow to add params from cmdline

docs/outputs.md

@@ -53,6 +53,24 @@ be used to set the compression level (between 1 and 9, default is 5). Compressio
 altogether with the CMake build option `PARTHENON_DISABLE_HDF5_COMPRESSION`. See the [build
 doc](building.md) for more details.
 
+## Gridh5 Files
+
+Gridh5 files are a special, lightweight subclass of HDF5 files, where
+only the grid structure is output. No variables are output, nor are
+coordinate values. Instead, the bounding box of each meshblock is
+output, along with all the metadata normally in a `phdf` file. A
+relevant block might look like:
+```
+<parthenon/output2>
+file_type = gridh5
+dt = 1.0
+```
+All options (except for the list of output variables) supported by
+`hdf5` are supported. The `plot_mesh.py` script in the
+`parthenon_tools` library found in the `scripts/python` folder can
+generate a constant `X3` slice of the mesh based on an `hdf5` or
+`gridh5` output file.
+
 ## Tuning HDF5 Performance
 
 Tuning IO parameters can be passed to Parthenon through the use of environment variables. Available environment variables are:

example/CMakeLists.txt

@@ -1,5 +1,5 @@
 #=========================================================================================
-# (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+# (C) (or copyright) 2020-2022. 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 @@
 add_subdirectory(stochastic_subgrid)
 add_subdirectory(advection)
 add_subdirectory(calculate_pi)
+add_subdirectory(count_cells)
 add_subdirectory(face_fields)
 add_subdirectory(kokkos_pi)
 add_subdirectory(particles)

example/count_cells/CMakeLists.txt

@@ -0,0 +1,24 @@
+#=========================================================================================
+# (C) (or copyright) 2022. 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
+# for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+# in the program are reserved by Triad National Security, LLC, and the U.S. Department
+# of Energy/National Nuclear Security Administration. The Government is granted for
+# itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+# license in this material to reproduce, prepare derivative works, distribute copies to
+# the public, perform publicly and display publicly, and to permit others to do so.
+#=========================================================================================
+
+get_property(DRIVER_LIST GLOBAL PROPERTY DRIVERS_USED_IN_TESTS)
+if( "count_cells" IN_LIST DRIVER_LIST OR NOT PARTHENON_DISABLE_EXAMPLES)
+  add_executable(
+      count-cells
+      count_cells_main.cpp
+      count_cells.hpp
+      count_cells.cpp
+  )
+  target_link_libraries(count-cells PRIVATE Parthenon::parthenon)
+  lint_target(count-cells)
+endif()

example/count_cells/count_cells.cpp

@@ -0,0 +1,143 @@
+//========================================================================================
+// (C) (or copyright) 2022. 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
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+// C++ includes
+#include <array>
+#include <fstream>
+#include <iostream>
+#include <memory>
+#include <string>
+
+// Parthenon Includes
+#include <globals.hpp>
+#include <interface/state_descriptor.hpp>
+#include <parthenon/package.hpp>
+
+#include "count_cells.hpp"
+
+using parthenon::ParameterInput;
+using parthenon::Params;
+using parthenon::Real;
+using parthenon::StateDescriptor;
+
+namespace count_cells {
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
+  auto package = std::make_shared<StateDescriptor>("count_cells");
+  Params &params = package->AllParams();
+
+  Real radius = pin->GetOrAddReal("count_cells", "sphere_radius", 1.0);
+  params.Add("radius", radius);
+
+  Real dx1_target = pin->GetOrAddReal("count_cells", "dx1_target", 0.1);
+  Real dx2_target = pin->GetOrAddReal("count_cells", "dx2_target", 0.1);
+  Real dx3_target = pin->GetOrAddReal("count_cells", "dx3_target", 0.1);
+  std::array<Real, 3> dx_target{dx1_target, dx2_target, dx3_target};
+  params.Add("dx_target", dx_target);
+
+  package->CheckRefinementBlock = CheckRefinement;
+
+  return package;
+}
+
+AmrTag CheckRefinement(MeshBlockData<Real> *rc) {
+  auto pmb = rc->GetBlockPointer();
+  auto pkg = pmb->packages.Get("count_cells");
+  const auto &coords = pmb->coords;
+  if (BlockInRegion(pkg.get(), pmb.get()) && !SufficientlyRefined(pkg.get(), coords)) {
+    return AmrTag::refine;
+  }
+  return AmrTag::same;
+}
+
+bool BlockInRegion(const StateDescriptor *pkg, MeshBlock *pmb) {
+  const auto radius = pkg->Param<Real>("radius");
+  const auto coords = pmb->coords;
+
+  auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  // Loop through all cell centers in the block to determine if block
+  // intersects with region. Technically this could still miss some.
+  // Could improve fidelity by also looping through cell face and nodes.
+  // Note this doesn't require actually allocating any cells, but you
+  // do have to loop through them.
+  int num_intersections = 0;
+  pmb->par_reduce(
+      "Check intersections of cells with regions", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(int k, int j, int i, int &n) {
+        Real x1 = coords.x1v(k, j, i);
+        Real x2 = coords.x2v(k, j, i);
+        Real x3 = coords.x3v(k, j, i);
+        if (std::sqrt(x1 * x1 + x2 * x2 + x3 * x3) <= radius) {
+          n += 1;
+        }
+      },
+      num_intersections);
+
+  return num_intersections > 0;
+}
+
+bool SufficientlyRefined(const StateDescriptor *pkg, const Coordinates_t &coords) {
+  const auto &dx_target = pkg->Param<std::array<Real, 3>>("dx_target");
+  for (int d = 1; d <= 3; ++d) {
+    // assumes uniform cartesian coordinates
+    // which have constant Dx accross whole meshblock
+    if (coords.Dx(d) > dx_target[d - 1]) return false;
+  }
+  return true;
+}
+
+void CountCells(const std::string &name, Mesh *pmesh, bool save_to_file) {
+  // a representative meshblock
+  auto pmb = pmesh->block_list[0];
+
+  const size_t mb_ncells_interior = pmb->cellbounds.GetTotal(IndexDomain::interior);
+  const size_t mb_ncells_total = pmb->cellbounds.GetTotal(IndexDomain::entire);
+  const size_t mb_ncells_ghost = mb_ncells_total - mb_ncells_interior;
+
+  // includes 3 flux buffers + coarse buffer + comm buffers
+  const size_t mb_ncells_with_extra_buffs = 5 * mb_ncells_total + mb_ncells_ghost;
+
+  const size_t num_blocks = pmesh->block_list.size();
+
+  Real ncells_interior = num_blocks * mb_ncells_interior;
+  Real ncells_total = num_blocks * mb_ncells_total;
+  Real ncells_ghost = num_blocks * mb_ncells_ghost;
+  Real ncells_with_extra_buffs = num_blocks * mb_ncells_with_extra_buffs;
+
+  if (parthenon::Globals::my_rank == 0) {
+    std::cout << std::scientific << "num blocks         = " << std::setw(14) << num_blocks
+              << "\n"
+              << "num cells interior = " << std::setw(14) << ncells_interior << "\n"
+              << "num cells total    = " << std::setw(14) << ncells_total << "\n"
+              << "num ghosts         = " << std::setw(14) << ncells_ghost << "\n"
+              << "num with comms etc = " << std::setw(14) << ncells_with_extra_buffs
+              << std::endl;
+
+    if (save_to_file) {
+      std::fstream fs;
+      fs.open(name, std::fstream::out);
+      if (!fs.is_open()) PARTHENON_THROW("Unable to write to summary file");
+      fs << std::scientific << "num blocks         = " << std::setw(14) << num_blocks
+         << "\n"
+         << "num cells interior = " << std::setw(14) << ncells_interior << "\n"
+         << "num cells total    = " << std::setw(14) << ncells_total << "\n"
+         << "num ghosts         = " << std::setw(14) << ncells_ghost << "\n"
+         << "num with comms etc = " << std::setw(14) << ncells_with_extra_buffs
+         << std::endl;
+      fs.close();
+    }
+  }
+}
+} // namespace count_cells

example/count_cells/count_cells.hpp

@@ -0,0 +1,39 @@
+//========================================================================================
+// (C) (or copyright) 2022. 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
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_COUNT_CELLS_COUNT_CELLS_HPP_
+#define EXAMPLE_COUNT_CELLS_COUNT_CELLS_HPP_
+
+// C++ Includes
+#include <memory>
+#include <string>
+
+// Parthenon Includes
+#include <coordinates/coordinates.hpp>
+#include <interface/state_descriptor.hpp>
+#include <parthenon/package.hpp>
+
+namespace count_cells {
+using namespace parthenon::package::prelude;
+using parthenon::Coordinates_t;
+using parthenon::Mesh;
+using parthenon::MeshBlock;
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
+AmrTag CheckRefinement(MeshBlockData<Real> *rc);
+bool BlockInRegion(const StateDescriptor *pkg, MeshBlock *pmb);
+bool SufficientlyRefined(const StateDescriptor *pkg, const Coordinates_t &coords);
+void CountCells(const std::string &name, Mesh *pmesh, bool save_to_file);
+
+} // namespace count_cells
+
+#endif // EXAMPLE_COUNT_CELLS_COUNT_CELLS_HPP_

example/count_cells/count_cells_main.cpp

@@ -0,0 +1,67 @@
+//========================================================================================
+// (C) (or copyright) 2022. 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
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+// C++ includes
+#include <iostream>
+#include <memory>
+#include <string>
+
+// Parthenon includes
+#include <basic_types.hpp>
+#include <outputs/outputs.hpp>
+#include <parthenon/driver.hpp>
+
+#include "count_cells.hpp"
+
+using namespace parthenon::driver::prelude;
+using parthenon::Outputs;
+using parthenon::SignalHandler::OutputSignal;
+
+Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
+  Packages_t packages;
+  packages.Add(count_cells::Initialize(pin.get()));
+  return packages;
+}
+
+int main(int argc, char *argv[]) {
+  ParthenonManager pman;
+  pman.app_input->ProcessPackages = ProcessPackages;
+
+  // Generates the mesh, which should refine appropriately
+  // despite no variables being set.
+  auto manager_status = pman.ParthenonInit(argc, argv);
+  if (manager_status == ParthenonStatus::complete) {
+    pman.ParthenonFinalize();
+    return 0;
+  }
+  if (manager_status == ParthenonStatus::error) {
+    pman.ParthenonFinalize();
+    return 1;
+  }
+  { // scoped so unique pointers cleaned up
+    // count cells
+    bool save_count_to_file =
+        pman.pinput->GetOrAddBoolean("count_cells", "save_to_file", false);
+    std::string countfile_name =
+        pman.pinput->GetOrAddString("count_cells", "filename", "cell_count.txt");
+    count_cells::CountCells(countfile_name, pman.pmesh.get(), save_count_to_file);
+
+    // Dump grid
+    std::unique_ptr<Outputs> pouts =
+        std::make_unique<Outputs>(pman.pmesh.get(), pman.pinput.get());
+    OutputSignal signal = OutputSignal::none;
+    pouts->MakeOutputs(pman.pmesh.get(), pman.pinput.get(), nullptr, signal);
+  }
+  pman.ParthenonFinalize();
+  return 0;
+}

example/count_cells/parthinput.count

@@ -0,0 +1,59 @@
+# ========================================================================================
+#  (C) (or copyright) 2022. 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
+#  for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+#  in the program are reserved by Triad National Security, LLC, and the U.S. Department
+#  of Energy/National Nuclear Security Administration. The Government is granted for
+#  itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+#  license in this material to reproduce, prepare derivative works, distribute copies to
+#  the public, perform publicly and display publicly, and to permit others to do so.
+# ========================================================================================
+
+<parthenon/job>
+problem_id = count_cells
+
+<parthenon/mesh>
+refinement = adaptive
+numlevel = 6
+
+nx1 = 64
+x1min = -1
+x1max = 1
+ix1_bc = periodic
+ox1_bc = periodic
+
+nx2 = 64
+x2min = -1
+x2max = 1
+ix2_bc = periodic
+ox2_bc = periodic
+
+nx3 = 64
+x3min = -1
+x3max = 1
+ix3_bc = periodic
+ox3_bc = periodic
+
+<parthenon/meshblock>
+nx1 = 32
+nx2 = 32
+nx3 = 32
+
+<count_cells>
+# AMR tries to resolve a sphere of a given radius
+# with a dx of dx_target in each direction.
+# It will do the best it can given the number of levels
+# requested.
+sphere_radius = 0.25
+dx1_target = 0.001
+dx2_target = 0.001
+dx3_target = 0.001
+
+# Change to true to save cell count to a file.
+save_to_file = false
+filename = cell_count.txt
+
+<parthenon/output0>
+file_type = gridh5