Merge branch 'lroberts36/multigrid-example-update' into lroberts36/mg…

…-riot-updates

CHANGELOG.md

@@ -4,8 +4,11 @@
 
 ### Added (new features/APIs/variables/...)
 - [[PR 968]](https://github.com/parthenon-hpc-lab/parthenon/pull/968) Add per package registration of boundary conditions
-- [[PR 948]](https://github.com/parthenon-hpc-lab/parthenon/pull/948) Add solver interface and update Poisson geometric multi-grid example 
+- [[PR 948]](https://github.com/parthenon-hpc-lab/parthenon/pull/948) Add solver interface and update Poisson geometric multi-grid example
+- [[PR 962]](https://github.com/parthenon-hpc-lab/parthenon/pull/962) Add support for in-situ histograms/profiles
 - [[PR 911]](https://github.com/parthenon-hpc-lab/parthenon/pull/911) Add infrastructure for geometric multi-grid
+- [[PR 971]](https://github.com/parthenon-hpc-lab/parthenon/pull/971) Add UserWorkBeforeLoop
+- [[PR 907]](https://github.com/parthenon-hpc-lab/parthenon/pull/907) PEP1: Allow subclassing StateDescriptor
 - [[PR 932]](https://github.com/parthenon-hpc-lab/parthenon/pull/932) Add GetOrAddFlag to metadata
 - [[PR 931]](https://github.com/parthenon-hpc-lab/parthenon/pull/931) Allow SparsePacks with subsets of blocks
 - [[PR 921]](https://github.com/parthenon-hpc-lab/parthenon/pull/921) Add more flexible ways of adding and using MeshData/MeshBlockData objects to DataCollections
@@ -20,6 +23,9 @@
 - [[PR 868]](https://github.com/parthenon-hpc-lab/parthenon/pull/868) Add block-local face, edge, and nodal fields and allow for packing
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 975]](https://github.com/parthenon-hpc-lab/parthenon/pull/975) Construct staged integrators via arbitrary name
+- [[PR 976]](https://github.com/parthenon-hpc-lab/parthenon/pull/976) Move UserWorkBeforeLoop to be after first output
+- [[PR 965]](https://github.com/parthenon-hpc-lab/parthenon/pull/965) Allow leading whitespace in input parameters
 - [[PR 926]](https://github.com/parthenon-hpc-lab/parthenon/pull/926) Internal refinement op registration
 - [[PR 897]](https://github.com/parthenon-hpc-lab/parthenon/pull/897) Deflate compression filter is not called any more if compression is soft disabled
 - [[PR 896]](https://github.com/parthenon-hpc-lab/parthenon/pull/896) Update Kokkos integration to support installed version. Use `serial` (flat MPI) host parallelization by default (instead of OpenMP)
@@ -38,6 +44,7 @@
 - [[PR 890]](https://github.com/parthenon-hpc-lab/parthenon/pull/890) Fix bugs in sparse communication and prolongation
 
 ### Infrastructure (changes irrelevant to downstream codes)
+- [[PR 967]](https://github.com/parthenon-hpc-lab/parthenon/pull/967) Change INLINE to FORCEINLINE on par_for_inner overloads
 - [[PR 938]](https://github.com/parthenon-hpc-lab/parthenon/pull/938) Restructure buffer packing/unpacking kernel hierarchical parallelism
 - [[PR 944]](https://github.com/parthenon-hpc-lab/parthenon/pull/944) Move sparse pack identifier creation to descriptor
 - [[PR 904]](https://github.com/parthenon-hpc-lab/parthenon/pull/904) Move to prolongation/restriction in one for AMR and communicate non-cell centered fields

CMakeLists.txt

@@ -62,9 +62,9 @@ include(cmake/Format.cmake)
 include(cmake/Lint.cmake)
 
 # regression test reference data
-set(REGRESSION_GOLD_STANDARD_VER 19 CACHE STRING "Version of gold standard to download and use")
+set(REGRESSION_GOLD_STANDARD_VER 20 CACHE STRING "Version of gold standard to download and use")
 set(REGRESSION_GOLD_STANDARD_HASH
-  "SHA512=e1d1a06b9cf9b761d42d0b6b241056ac75658db90138b6b867b1ca7ead4308af4f980285af092b40aee1dbbfb68b4e8cb15efcc9b83d7930c18bf992ae95c729"
+  "SHA512=e5e421f3c0be01e4708965542bb8b1b79b5c96de97091e46972e375c7616588d026a9a8e29226d9c7ef75346bc859fd9af72acdc7e95e0d783b5ef29aa4630b1"
   CACHE STRING "Hash of default gold standard file to download")
 option(REGRESSION_GOLD_STANDARD_SYNC "Automatically sync gold standard files." ON)
 
@@ -353,6 +353,10 @@ if (PARTHENON_ENABLE_ASCENT)
   find_package(Ascent REQUIRED NO_DEFAULT_PATH)
 endif()
 
+# Installation configuration
+include(GNUInstallDirs)
+set(CMAKE_INSTALL_INCLUDEDIR "${CMAKE_INSTALL_INCLUDEDIR}/parthenon")
+
 add_subdirectory(src)
 add_subdirectory(example)
 add_subdirectory(benchmarks)

cmake/parthenonConfig.cmake.in

@@ -11,6 +11,14 @@
 # the public, perform publicly and display publicly, and to permit others to do so.
 #=========================================================================================
 
+if(NOT PARTHENON_CMAKE)
+  cmake_path(SET PARTHENON_CMAKE_BASE_DIR NORMALIZE "${CMAKE_CURRENT_LIST_DIR}/..")
+  message(STATUS "Appending parthenon cmake module directory: " ${PARTHENON_CMAKE_BASE_DIR})
+  list(APPEND CMAKE_MODULE_PATH ${PARTHENON_CMAKE_BASE_DIR})
+  set(PARTHENON_CMAKE TRUE)
+endif()
+
+
 # Favor using the kokkos package that was built with parthenon, if one has not been specified
 if (@PARTHENON_IMPORT_KOKKOS@)
   find_package(Kokkos 3 REQUIRED PATHS @Kokkos_DIR@ NO_DEFAULT_PATH)
@@ -37,4 +45,6 @@ if(@OpenMP_FOUND@)
   find_package(OpenMP REQUIRED COMPONENTS CXX)
 endif()
 
+find_package(Filesystem REQUIRED COMPONENTS Experimental Final)
+
 include("${CMAKE_CURRENT_LIST_DIR}/parthenonTargets.cmake")

doc/sphinx/src/boundary_communication.rst

@@ -43,13 +43,14 @@ subset, and the columns point to the indices in the ``bnd_info`` array
 containing the subset of sub-halos you wish to operate on.
 
 To communicate across a particular boundary type, the templated 
-boundary communication routines (see :boundary_comm_tasks:`boundary_comm_tasks`.) 
+boundary communication routines (see :ref:`boundary_comm_tasks`) 
 should be instantiated with the desired ``BoundaryType``, i.e. 
 
 .. code:: cpp 
+
   SendBoundBufs<BoundaryType::gmg_restrict_send>(md); 
 
-The different ``BoundaryType``s are: 
+The different ``BoundaryType``\ s are: 
 
 - ``any``: Communications are performed between all leaf blocks (i.e. the
   standard Parthenon grid that does not include multi-grid related blocks). 

doc/sphinx/src/interface/state.rst

@@ -81,7 +81,7 @@ several useful features and functions.
   with specified fields to the ``DataCollection`` objects in ``Mesh`` and
   ``MeshBlock``.  For convenience, the ``Mesh`` class also provides this
   function, which provides a list of variables gathered from all the package
-  ``StateDescriptor``s.
+  ``StateDescriptor``\s.
 - ``void FillDerivedBlock(MeshBlockData<Real>* rc)`` delgates to the
   ``std::function`` member ``FillDerivedBlock`` if set (defaults to
   ``nullptr`` and therefore a no-op) that allows an application to provide
@@ -112,6 +112,13 @@ several useful features and functions.
   deletgates to the ``std::function`` member ``PostStepDiagnosticsMesh``
   if set (defaults to ``nullptr`` an therefore a no-op) to print
   diagnostics after the time-integration advance
+- ``void UserWorkBeforeLoopMesh(Mesh *, ParameterInput *pin, SimTime
+  &tm)`` performs a per-package, mesh-wide calculation after the mesh
+  has been generated, and problem generators called, but before any
+  time evolution. This work is done both on first initialization and
+  on restart. If you would like to avoid doing the work upon restart,
+  you can check for the const ``is_restart`` member field of the ``Mesh``
+  object.
 
 The reasoning for providing ``FillDerived*`` and ``EstimateTimestep*``
 function pointers appropriate for usage with both ``MeshData`` and

doc/sphinx/src/mesh/mesh.rst

@@ -85,6 +85,7 @@ To work with these GMG levels, ``MeshData`` objects containing these blocks can
 be recovered from a ``Mesh`` pointer using 
 
 .. code:: c++
+
   auto &md = pmesh->gmg_mesh_data[level].GetOrAdd(level, "base", partition_idx);
 
 This ``MeshData`` will include blocks at the current level and possibly some 
@@ -94,6 +95,7 @@ communication). To make packs containing only a subset of blocks from a
 GMG ``MeshData`` pointer ``md``, one would use 
 
 .. code:: c++
+
   int nblocks = md->NumBlocks();
   std::vector<bool> include_block(nblocks, true);
   for (int b = 0; b < nblocks; ++b)
@@ -104,10 +106,10 @@ GMG ``MeshData`` pointer ``md``, one would use
   auto pack = desc.GetPack(md.get(), include_block);
 
 In addition to creating the ``LogicalLocation`` and block lists for the GMG levels, 
-``Mesh`` fills neigbor arrays in ``MeshBlock`` for intra- and inter-GMG block list 
+``Mesh`` fills neighbor arrays in ``MeshBlock`` for intra- and inter-GMG block list 
 communication (i.e. boundary communication and internal prolongation/restriction, 
 respectively). Communication within and between GMG levels can be done by calling 
 boundary communication routines with the boundary tags ``gmg_same``, 
 ``gmg_restrict_send``, ``gmg_restrict_recv``, ``gmg_prolongate_send``, 
-``gmg_prolongate_recv`` (see :boundary_communication:`boundary_communication`). 
+``gmg_prolongate_recv`` (see :ref:`boundary_comm_tasks`). 
 

doc/sphinx/src/outputs.rst

@@ -10,7 +10,7 @@ To disable an output block without removing it from the intput file set
 the block's ``dt < 0.0``.
 
 In addition to time base outputs, two additional options to trigger
-outputs (applies to HDF5 and restart outputs) exist.
+outputs (applies to HDF5, restart and histogram outputs) exist.
 
 -  Signaling: If ``Parthenon`` catches a signal, e.g., ``SIGALRM`` which
    is often sent by schedulers such as Slurm to signal a job of
@@ -28,7 +28,10 @@ outputs (applies to HDF5 and restart outputs) exist.
 
 Note, in both cases the original numbering of the output will be
 unaffected and the ``final`` and ``now`` files will be overwritten each
-time without warning. ## HDF5
+time without warning.
+
+HDF5
+----
 
 Parthenon allows users to select which fields are captured in the HDF5
 (``.phdf``) dumps at runtime. In the input file, include a
@@ -158,6 +161,168 @@ 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`.
 
+Histograms
+----------
+
+Parthenon supports calculating flexible 1D and 2D histograms in-situ that
+are written to disk in HDF5 format.
+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  
+- weighting by volume and/or variable
+
+The output format follows ``numpy`` convention, so that plotting data
+with Python based machinery should be straightfoward (see example below).
+In other words, 2D histograms use C-ordering corresponding to ``[x,y]``
+indexing with ``y`` being the fast index.
+In general, histograms are calculated using inclusive left bin edges and
+data equal to the rightmost edge is also included in the last bin.
+
+A ``<parthenon/output*>`` block containing one simple and one complex
+example might look like::
+
+   <parthenon/output8>
+   file_type = histogram            # required, sets the output type
+   dt = 1.0                         # required, sets the output interval
+   hist_names = myname, other_name  # required, specifies the names of the histograms
+                                    # in this block (used a prefix below and in the output)
+
+   # 1D histogram ("standard", i.e., counting occurance in bin)
+   myname_ndim = 1
+   myname_x_variable = advected
+   myname_x_variable_component = 0
+   myname_x_edges_type = log
+   myname_x_edges_num_bins = 10
+   myname_x_edges_min = 1e-9
+   myname_x_edges_max = 1e0
+   myname_binned_variable = HIST_ONES
+
+   # 2D histogram of volume weighted variable according to two coordinates
+   other_name_ndim = 2
+   other_name_x_variable = HIST_COORD_X1
+   other_name_x_edges_type = list
+   other_name_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5
+   other_name_y_variable = HIST_COORD_X2
+   other_name_y_edges_type = list
+   other_name_y_edges_list = -0.5, -0.1, 0.0, 0.1, 0.5
+   other_name_binned_variable = advected
+   other_name_binned_variable_component = 0
+   other_name_weight_by_volume = true
+   other_name_weight_variable = one_minus_advected_sq
+   other_name_weight_variable_component = 0
+
+with the following parameters
+
+- ``hist_names=STRING, STRING, STRING, ...`` (comma separated names)
+   The names of the histograms in this block.
+   Will be used as preifx in the block as well as in the output file.
+   All histograms will be written to the same output file with the "group" in the
+   output corresponding to the histogram name.
+- ``NAME_ndim=INT`` (either ``1`` or ``2``)
+   Dimensionality of the histogram.
+- ``NAME_x_variable=STRING`` (variable name or special coordinate string ``HIST_COORD_X1``, ``HIST_COORD_X2``, ``HIST_COORD_X3`` or ``HIST_COORD_R``)
+   Variable to be used as bin. If a variable name is given a component has to be specified, too,
+   see next parameter.
+   For a scalar variable, the component needs to be specified as ``0`` anyway.
+   If binning should be done by coordinate the special strings allow to bin by either one
+   of the three dimensions or by radial distance from the origin.
+- ``NAME_x_variable_component=INT``
+   Component index of the binning variable.
+   Used/required only if a non-coordinate variable is used for binning.
+- ``NAME_x_edges_type=STRING`` (``lin``, ``log``, or ``list``)
+   How the bin edges are defined in the first dimension.
+   For ``lin`` and ``log`` direct indexing is used to determine the bin, which is significantly
+   faster than specifying the edges via a ``list`` as the latter requires a binary search.
+- ``NAME_x_edges_min=FLOAT``
+   Minimum value (inclusive) of the bins in the first dim.
+   Used/required only for ``lin`` and ``log`` edge type.
+- ``NAME_x_edges_max=FLOAT``
+   Maximum value (inclusive) of the bins in the first dim.
+   Used/required only for ``lin`` and ``log`` edge type.
+- ``NAME_x_edges_num_bins=INT`` (must be ``>=1``)
+   Number of equally spaced bins between min and max value in the first dim.
+   Used/required only for ``lin`` and ``log`` edge type.
+- ``NAME_x_edges_list=FLOAT,FLOAT,FLOAT,...`` (comma separated list of increasing values)
+   Arbitrary definition of edge values with inclusive innermost and outermost edges.
+   Used/required only for ``list`` edge type.
+- ``NAME_y_edges...``
+   Same as the ``NAME_x_edges...`` parameters except for being used in the second
+   dimension for ``ndim=2`` histograms.
+- ``NAME_accumulate=BOOL`` (``true`` or ``false`` default)
+   Accumulate data that is outside the binning range in the outermost bins.
+- ``NAME_binned_variable=STRING`` (variable name or ``HIST_ONES``)
+   Variable to be binned. If a variable name is given a component has to be specified, too,
+   see next parameter.
+   For a scalar variable, the component needs to be specified as ``0`` anyway.
+   If sampling (i.e., counting the number of value inside a bin) is to be used the special
+   string ``HIST_ONES`` can be set.
+- ``NAME_binned_variable_component=INT``
+   Component index of the variable to be binned.
+   Used/required only if a variable is binned and not ``HIST_ONES``.
+- ``NAME_weight_by_volume=BOOL`` (``true`` or ``false``)
+   Apply volume weighting to the binned variable. Can be used simultaneously with binning
+   by a different variable. Note that this does *not* include any normalization
+   (e.g., by total volume or the sum of the weight variable in question) and is left to
+   the user during post processing.
+- ``NAME_weight_variable=STRING``
+   Variable to be used as weight.
+   Can be used together with volume weighting.
+   For a scalar variable, the component needs to be specified as ``0`` anyway.
+- ``NAME_weight_variable_component=INT``
+   Component index of the variable to be used as weight.
+
+Note, weighting by volume and variable simultaneously might seem counterintuitive, but
+easily allows for, e.g., mass-weighted profiles, by enabling weighting by volume and
+using a mass density field as additional weight variable.
+
+In practice, a 1D histogram in the astrophysical context may look like (top panel from
+Fig 4 in `Curtis et al 2023 ApJL 945 L13 <https://dx.doi.org/10.3847/2041-8213/acba16>`_):
+
+.. figure:: figs/Curtis_et_al-ApJL-2023-1dhist.png
+   :alt: 1D histogram example from Fig 2 in Curtis et al 2023 ApJL 945 L13
+
+Translating this to the notation used for Parthenon histogram outputs means specifying
+for each histogram
+
+- the field containing the Electron fraction as ``x_variable``\ ,
+- the field containing the traced mass density as ``binned_variable``\ , and
+- enable ``weight_by_volume`` (to get the total traced mass).
+
+Similarly, a 2D histogram (also referred to as phase plot) example may look like
+(from the `yt Project documentation <https://yt-project.org/doc/visualizing/plots.html#d-phase-plots>`_):
+
+.. figure:: figs/yt_doc-2dhist.png
+   :alt: 2D histogram example from the yt documentation
+
+Translating this to the notation used for Parthenon histogram outputs means using
+
+- the field containing the density as ``x_variable``\ ,
+- the field containing the temperature as ``y_variable``\ ,
+- the field containing the mass density as ``binned_variable``\ , and
+- enable ``weight_by_volume`` (to get the total mass).
+
+
+
+The following is a minimal example to plot a 1D and 2D histogram from the output file:
+
+.. code:: python
+
+   with h5py.File("parthenon.out8.histograms.00040.hdf", "r") as infile:
+     # 1D histogram
+      x = infile["myname/x_edges"][:]
+      y = infile["myname/data"][:]
+      plt.plot(x, y)
+      plt.show()
+
+      # 2D histogram
+      x = infile["other_name/x_edges"][:]
+      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()   
+
 Ascent (optional)
 -----------------
 

doc/sphinx/src/tasks.rst

@@ -117,7 +117,7 @@ finally another round of asynchronous work.
 A diagram illustrating the relationship between these different classes
 is shown below.
 
-.. figure:: TaskDiagram.png
+.. figure:: figs/TaskDiagram.png
    :alt: Task Diagram
 
 ``TaskCollection`` provides two member functions, ``AddRegion`` and

example/advection/advection_package.cpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2023. 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
@@ -13,12 +13,14 @@
 
 #include <algorithm>
 #include <cmath>
+#include <iostream>
 #include <limits>
 #include <memory>
 #include <string>
 #include <vector>
 
 #include <coordinates/coordinates.hpp>
+#include <globals.hpp>
 #include <parthenon/package.hpp>
 
 #include "advection_package.hpp"
@@ -215,10 +217,19 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   }
   pkg->CheckRefinementBlock = CheckRefinement;
   pkg->EstimateTimestepBlock = EstimateTimestepBlock;
+  pkg->UserWorkBeforeLoopMesh = AdvectionGreetings;
 
   return pkg;
 }
 
+void AdvectionGreetings(Mesh *pmesh, ParameterInput *pin, parthenon::SimTime &tm) {
+  if (parthenon::Globals::my_rank == 0) {
+    std::cout << "Hello from the advection package in the advection example!\n"
+              << "This run is a restart: " << pmesh->is_restart << "\n"
+              << std::endl;
+  }
+}
+
 AmrTag CheckRefinement(MeshBlockData<Real> *rc) {
   // refine on advected, for example.  could also be a derived quantity
   auto pmb = rc->GetBlockPointer();