stopping point

example/count_cells/count_cells.cpp

@@ -20,6 +20,12 @@
 
 #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();
@@ -42,13 +48,13 @@ 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) && !SufficientlyRefined(pkg,coords)) {
+  if (BlockInRegion(pkg.get(), pmb.get()) && !SufficientlyRefined(pkg.get(), coords)) {
     return AmrTag::refine;
   }
   return AmrTag::same;
 }
 
-bool BlockInRegion(const StateDescriptor *pkg, const MeshBlock *pmb) {
+bool BlockInRegion(const StateDescriptor *pkg, MeshBlock *pmb) {
   const auto radius = pkg->Param<Real>("radius");
   const auto coords = pmb->coords;
 
@@ -82,33 +88,35 @@ bool SufficientlyRefined(const StateDescriptor *pkg, const Coordinates_t &coords
   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;
+    if (coords.Dx(d) > dx_target[d - 1]) return false;
   }
   return true;
 }
 
 void CountCells(Mesh *pmesh) {
   // 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 = ncells_total - ncells_interior;
+  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*ncells_total + ncells_ghost;
+  const size_t mb_ncells_with_extra_buffs = 5 * mb_ncells_total + mb_ncells_ghost;
 
   const size_t num_blocks = pmesh->block_list.size();
 
-  size_t ncells_interior = num_blocks * mb_ncells_interior;
-  size_t ncells_total = num_blocks * mb_ncells_total;
-  size_t ncells_ghost = num_blocks * mb_ncells_ghost;
-  size_t ncells_with_extra_buffs = mb_ncells_with_extra_buffs;
-
-  std::cout << "num blocks         = " << num_blocks << "\n"
-            << "num cells interior = " << ncells_interior << "\n"
-            << "num cells total    = " << ncells_total << "\n"
-            << "num ghosts         = " << ncells_ghost << "\n"
-            << "num with comms etc = " << ncells_with_extra_buffs
-            << std::endl;
+  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;
+
+  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;
 }
+} // namespace count_cells

example/count_cells/count_cells.hpp

@@ -10,7 +10,7 @@
 // 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.
 //========================================================================================
-#infdef EXAMPLE_COUNT_CELLS_HPP_
+#ifndef EXAMPLE_COUNT_CELLS_HPP_
 #define EXAMPLE_COUNT_CELLS_HPP_
 
 // Parthenon Includes
@@ -26,7 +26,7 @@ using parthenon::MeshBlock;
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
 AmrTag CheckRefinement(MeshBlockData<Real> *rc);
-bool BlockInRegion(const StateDescriptor *pkg, const MeshBlock *pmb);
+bool BlockInRegion(const StateDescriptor *pkg, MeshBlock *pmb);
 bool SufficientlyRefined(const StateDescriptor *pkg, const Coordinates_t &coords);
 void CountCells(Mesh *pmesh);
 

example/count_cells/count_cells_main.cpp

@@ -13,21 +13,26 @@
 
 // C++ includes
 #include <iostream>
+#include <memory>
 
 // 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;') {
+int main(int argc, char *argv[]) {
   ParthenonManager pman;
   pman.app_input->ProcessPackages = ProcessPackages;
 
@@ -42,10 +47,16 @@ int main(int argc, char *argv;') {
     pman.ParthenonFinalize();
     return 1;
   }
-
-  count_cells::CountCells(pman.pmesh.get());
-
+  { // scoped so unique pointers cleaned up
+    // count cells
+    count_cells::CountCells(pman.pmesh.get());
+
+    // 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,51 @@
+# ========================================================================================
+#  (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 = 5
+
+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>
+sphere_radius = 0.5
+dx1_target = 0.001
+dx2_target = 0.001
+dx3_target = 0.001
+
+<parthenon/output0>
+file_type = rst

scripts/python/packages/parthenon_tools/parthenon_tools/plot_mesh.py

@@ -0,0 +1,54 @@
+# =========================================================================================
+# (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.
+# =========================================================================================
+
+from argparse import ArgumentParser
+import numpy as np
+import h5py
+import matplotlib as mpl
+from matplotlib import pyplot as plt
+import matplotlib.patches as mpatches
+
+parser = ArgumentParser("Plot x3 slice of Parthenon mesh")
+parser.add_argument('file', type=str,
+                    help='rhdf or phdf file to use for mesh')
+parser.add_argument('--slice', type=float, default=None,
+                    help='Slice through X3 to plot. Default is middle.')
+parser.add_argument('--save', type=str, default='mesh.png',
+                    help='File to save figure as. Defaults to mesh.png')
+args = parser.parse_args()
+
+with h5py.File(args.file,'r') as f:
+    domain = f['Info'].attrs['RootGridDomain'].reshape(3,3)
+    NB = f['Info'].attrs['NumMeshBlocks']
+    if args.slice is not None:
+        zslice = args.slice
+    else:
+        zslice = 0.5(domain[2,0] + domain[2,1])
+    fig,ax = plt.subplots()
+    for b in range(NB):
+        if f['Locations/z'][b,0] <= zslice <= f['Locations/z'][b,-1]:
+            xb = (f['Locations/x'][b,0], f['Locations/x'][b,-1])
+            yb = (f['Locations/y'][b,0], f['Locations/y'][b,-1])
+            rect = mpatches.Rectangle((xb[0],yb[0]),
+                                      xb[1]-xb[0],
+                                      yb[1]-yb[0],
+                                      linewidth=0.225,
+                                      alpha=1,
+                                      edgecolor='k',
+                                      facecolor='none',
+                                      linestyle='solid')
+            ax.add_patch(rect)
+    plt.xlim(domain[0,0],domain[0,1])
+    plt.ylim(domain[1,0],domain[1,1])
+    plt.savefig(args.save,bbox_inches='tight')
+