diff --git a/CHANGELOG.md b/CHANGELOG.md index 53cf87782faa..ba5c56f56046 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[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 @@ -21,6 +22,7 @@ ### Changed (changing behavior/API/variables/...) - [[PR 974]](https://github.com/parthenon-hpc-lab/parthenon/pull/974) Change GetParentPointer to always return T* +- [[PR 976]](https://github.com/parthenon-hpc-lab/parthenon/pull/976) Move UserWorkBeforeLoop to be after first o - [[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 diff --git a/doc/sphinx/src/boundary_communication.rst b/doc/sphinx/src/boundary_communication.rst index 3df126de6f02..cc85fa8052c2 100644 --- a/doc/sphinx/src/boundary_communication.rst +++ b/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(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). diff --git a/doc/sphinx/src/figs/Curtis_et_al-ApJL-2023-1dhist.png b/doc/sphinx/src/figs/Curtis_et_al-ApJL-2023-1dhist.png new file mode 100644 index 000000000000..b7078512e10b Binary files /dev/null and b/doc/sphinx/src/figs/Curtis_et_al-ApJL-2023-1dhist.png differ diff --git a/doc/sphinx/src/TaskDiagram.png b/doc/sphinx/src/figs/TaskDiagram.png similarity index 100% rename from doc/sphinx/src/TaskDiagram.png rename to doc/sphinx/src/figs/TaskDiagram.png diff --git a/doc/sphinx/src/figs/yt_doc-2dhist.png b/doc/sphinx/src/figs/yt_doc-2dhist.png new file mode 100644 index 000000000000..ae9c6d1b2f73 Binary files /dev/null and b/doc/sphinx/src/figs/yt_doc-2dhist.png differ diff --git a/doc/sphinx/src/interface/state.rst b/doc/sphinx/src/interface/state.rst index 0841e8c5a1c7..e0760d4b8ff1 100644 --- a/doc/sphinx/src/interface/state.rst +++ b/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* rc)`` delgates to the ``std::function`` member ``FillDerivedBlock`` if set (defaults to ``nullptr`` and therefore a no-op) that allows an application to provide diff --git a/doc/sphinx/src/mesh/mesh.rst b/doc/sphinx/src/mesh/mesh.rst index 1c046bc14c87..aa891633cdde 100644 --- a/doc/sphinx/src/mesh/mesh.rst +++ b/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 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`). diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index 574d52f4a212..c01c957fdcf5 100644 --- a/doc/sphinx/src/outputs.rst +++ b/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 ```` block containing one simple and one complex +example might look like:: + + + 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 `_): + +.. 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 `_): + +.. 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) ----------------- diff --git a/doc/sphinx/src/tasks.rst b/doc/sphinx/src/tasks.rst index 77aa79325a97..d4c0b361b7f9 100644 --- a/doc/sphinx/src/tasks.rst +++ b/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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 512201459106..a3a0ad91c52c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -174,6 +174,7 @@ add_library(parthenon mesh/meshblock.cpp outputs/ascent.cpp + outputs/histogram.cpp outputs/history.cpp outputs/io_wrapper.cpp outputs/io_wrapper.hpp diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp index 03acf820e536..c0e21414e576 100644 --- a/src/driver/driver.cpp +++ b/src/driver/driver.cpp @@ -65,15 +65,6 @@ DriverStatus EvolutionDriver::Execute() { PreExecute(); InitializeBlockTimeStepsAndBoundaries(); SetGlobalTimeStep(); - OutputSignal signal = OutputSignal::none; - pouts->MakeOutputs(pmesh, pinput, &tm, signal); - pmesh->mbcnt = 0; - int perf_cycle_offset = - pinput->GetOrAddInteger("parthenon/time", "perf_cycle_offset", 0); - - // Output a text file of all parameters at this point - // Defaults must be set across all ranks - DumpInputParameters(); // Before loop do work // App input version @@ -81,13 +72,22 @@ DriverStatus EvolutionDriver::Execute() { if (app_input->UserWorkBeforeLoop != nullptr) { app_input->UserWorkBeforeLoop(pmesh, pinput, tm); } - // packages version for (auto &[name, pkg] : pmesh->packages.AllPackages()) { pkg->UserWorkBeforeLoop(pmesh, pinput, tm); } Kokkos::Profiling::popRegion(); // Driver_UserWorkBeforeLoop + OutputSignal signal = OutputSignal::none; + pouts->MakeOutputs(pmesh, pinput, &tm, signal); + pmesh->mbcnt = 0; + int perf_cycle_offset = + pinput->GetOrAddInteger("parthenon/time", "perf_cycle_offset", 0); + + // Output a text file of all parameters at this point + // Defaults must be set across all ranks + DumpInputParameters(); + Kokkos::Profiling::pushRegion("Driver_Main"); while (tm.KeepGoing()) { if (Globals::my_rank == 0) OutputCycleDiagnostics(); diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp new file mode 100644 index 000000000000..88c6ced2ecd3 --- /dev/null +++ b/src/outputs/histogram.cpp @@ -0,0 +1,570 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2023 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 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 +// 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. +//======================================================================================== +//! \file histogram.cpp +// \brief 1D and 2D histograms + +// options for building +#include "basic_types.hpp" +#include "config.hpp" +#include "globals.hpp" +#include "kokkos_abstraction.hpp" +#include "parameter_input.hpp" +#include "parthenon_array_generic.hpp" + +// Only proceed if HDF5 output enabled +#ifdef ENABLE_HDF5 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Parthenon headers +#include "coordinates/coordinates.hpp" +#include "defs.hpp" +#include "interface/variable_state.hpp" +#include "mesh/mesh.hpp" +#include "outputs/output_utils.hpp" +#include "outputs/outputs.hpp" +#include "outputs/parthenon_hdf5.hpp" +#include "utils/error_checking.hpp" +#include "utils/sort.hpp" // for upper_bound + +// ScatterView is not part of Kokkos core interface +#include "Kokkos_ScatterView.hpp" + +namespace parthenon { + +using namespace OutputUtils; + +namespace HistUtil { + +// Parse input for x and y vars from input +// std::tuple +auto ProcessVarInput(ParameterInput *pin, const std::string &block_name, + const std::string &prefix) { + auto var_name = pin->GetString(block_name, prefix + "variable"); + int var_component = -1; + VarType var_type; + if (var_name == "HIST_COORD_X1") { + var_type = VarType::X1; + } else if (var_name == "HIST_COORD_X2") { + var_type = VarType::X2; + } else if (var_name == "HIST_COORD_X3") { + var_type = VarType::X3; + } else if (var_name == "HIST_COORD_R") { + PARTHENON_REQUIRE_THROWS( + typeid(Coordinates_t) == typeid(UniformCartesian), + "Radial coordinate currently only works for uniform Cartesian coordinates."); + var_type = VarType::R; + } else { + var_type = VarType::Var; + var_component = pin->GetInteger(block_name, prefix + "variable_component"); + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(var_component >= 0, + "Negative component indices are not supported"); + } + + return std::make_tuple(var_name, var_component, var_type); +} + +// Parse edges from input parameters. Returns the edges themselves (to be used as list for +// arbitrary bins) as well as min and step sizes (potentially in log space) for direct +// indexing. +auto GetEdges(ParameterInput *pin, const std::string &block_name, + const std::string &prefix) { + std::vector edges_in; + auto edge_type = EdgeType::Undefined; + auto edge_min = std::numeric_limits::quiet_NaN(); + auto edge_dbin = std::numeric_limits::quiet_NaN(); + + const auto edge_type_str = pin->GetString(block_name, prefix + "type"); + if (edge_type_str == "lin" || edge_type_str == "log") { + edge_min = pin->GetReal(block_name, prefix + "min"); + const auto edge_max = pin->GetReal(block_name, prefix + "max"); + PARTHENON_REQUIRE_THROWS(edge_max > edge_min, + "Histogram max needs to be larger than min.") + + const auto edge_num_bins = pin->GetInteger(block_name, prefix + "num_bins"); + PARTHENON_REQUIRE_THROWS(edge_num_bins >= 1, "Need at least one bin for histogram."); + edges_in.reserve(edge_num_bins); + + if (edge_type_str == "lin") { + edge_type = EdgeType::Lin; + edge_dbin = (edge_max - edge_min) / (edge_num_bins); + for (int i = 0; i < edge_num_bins; i++) { + edges_in.emplace_back(edge_min + i * edge_dbin); + } + edges_in.emplace_back(edge_max); + } else if (edge_type_str == "log") { + edge_type = EdgeType::Log; + PARTHENON_REQUIRE_THROWS( + edge_min > 0.0 && edge_max > 0.0, + "Log binning for negative values not implemented. However, you can specify " + "arbitrary bin edges through the 'list' edge type.") + + // override start with log value for direct indexing in histogram kernel + edge_min = std::log10(edge_min); + edge_dbin = (std::log10(edge_max) - edge_min) / (edge_num_bins); + for (int i = 0; i < edge_num_bins; i++) { + edges_in.emplace_back(std::pow(10., edge_min + i * edge_dbin)); + } + edges_in.emplace_back(edge_max); + } else { + PARTHENON_FAIL("Not sure how I got here...") + } + + } else if (edge_type_str == "list") { + edge_type = EdgeType::List; + edges_in = pin->GetVector(block_name, prefix + "list"); + // required by binning index function + PARTHENON_REQUIRE_THROWS(std::is_sorted(edges_in.begin(), edges_in.end()), + "Bin edges must be in order."); + PARTHENON_REQUIRE_THROWS(edges_in.size() >= 2 && edges_in[1] > edges_in[0], + "Need at least one bin, i.e., two distinct edges."); + + } else { + PARTHENON_THROW( + "Unknown edge type for histogram. Supported types are lin, log, and list.") + } + auto edges = ParArray1D(prefix, edges_in.size()); + auto edges_h = edges.GetHostMirror(); + for (int i = 0; i < edges_in.size(); i++) { + edges_h(i) = edges_in[i]; + } + Kokkos::deep_copy(edges, edges_h); + + PARTHENON_REQUIRE_THROWS( + edge_type != EdgeType::Undefined, + "Edge type not set and it's unclear how this code was triggered..."); + return std::make_tuple(edges, edge_type, edge_min, edge_dbin); +} + +Histogram::Histogram(ParameterInput *pin, const std::string &block_name, + const std::string &name) { + name_ = name; + const auto prefix = name + "_"; + + ndim_ = pin->GetInteger(block_name, prefix + "ndim"); + PARTHENON_REQUIRE_THROWS(ndim_ == 1 || ndim_ == 2, "Histogram dim must be '1' or '2'"); + + std::tie(x_var_name_, x_var_component_, x_var_type_) = + ProcessVarInput(pin, block_name, prefix + "x_"); + + std::tie(x_edges_, x_edges_type_, x_edge_min_, x_edge_dbin_) = + GetEdges(pin, block_name, prefix + "x_edges_"); + + // For 1D profile default initalize y variables + y_var_name_ = ""; + y_var_component_ = -1; + y_var_type_ = VarType::Unused; + // and for 2D profile check if they're explicitly set (not default value) + if (ndim_ == 2) { + std::tie(y_var_name_, y_var_component_, y_var_type_) = + ProcessVarInput(pin, block_name, prefix + "y_"); + + std::tie(y_edges_, y_edges_type_, y_edge_min_, y_edge_dbin_) = + GetEdges(pin, block_name, prefix + "y_edges_"); + + } else { + y_edges_ = ParArray1D(prefix + "y_edges_unused", 0); + } + + binned_var_name_ = + pin->GetOrAddString(block_name, prefix + "binned_variable", "HIST_ONES"); + binned_var_component_ = -1; // implies that we're not binning a variable but count + if (binned_var_name_ != "HIST_ONES") { + binned_var_component_ = + pin->GetInteger(block_name, prefix + "binned_variable_component"); + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(binned_var_component_ >= 0, + "Negative component indices are not supported"); + } + + const auto nxbins = x_edges_.extent_int(0) - 1; + const auto nybins = ndim_ == 2 ? y_edges_.extent_int(0) - 1 : 1; + + result_ = ParArray2D(prefix + "result", nybins, nxbins); + scatter_result = + Kokkos::Experimental::ScatterView(result_.KokkosView()); + + accumulate_ = pin->GetOrAddBoolean(block_name, prefix + "accumulate", false); + weight_by_vol_ = pin->GetOrAddBoolean(block_name, prefix + "weight_by_volume", false); + + weight_var_name_ = + pin->GetOrAddString(block_name, prefix + "weight_variable", "HIST_ONES"); + weight_var_component_ = -1; // implies that weighting is not applied + if (weight_var_name_ != "HIST_ONES") { + weight_var_component_ = + pin->GetInteger(block_name, prefix + "weight_variable_component"); + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(weight_var_component_ >= 0, + "Negative component indices are not supported"); + } +} + +// Computes a 1D or 2D histogram with inclusive lower edges and inclusive rightmost edges. +// Function could in principle be templated on dimension, but it's currently not expected +// to be a performance concern (because it won't be called that often). +void Histogram::CalcHist(Mesh *pm) { + Kokkos::Profiling::pushRegion("Calculate single histogram"); + const auto x_var_component = x_var_component_; + const auto y_var_component = y_var_component_; + const auto binned_var_component = binned_var_component_; + const auto weight_var_component = weight_var_component_; + const auto x_var_type = x_var_type_; + const auto y_var_type = y_var_type_; + const auto x_edges = x_edges_; + const auto y_edges = y_edges_; + const auto x_edges_type = x_edges_type_; + const auto y_edges_type = y_edges_type_; + const auto x_edge_min = x_edge_min_; + const auto x_edge_dbin = x_edge_dbin_; + const auto y_edge_min = y_edge_min_; + const auto y_edge_dbin = y_edge_dbin_; + const auto hist_ndim = ndim_; + const auto weight_by_vol = weight_by_vol_; + const auto accumulate = accumulate_; + auto result = result_; + auto scatter = scatter_result; + + // Reset ScatterView from previous output + scatter.reset(); + // Also reset the histogram from previous call. + // Currently still required for consistent results between host and device backends, see + // https://github.com/kokkos/kokkos/issues/6363 + Kokkos::deep_copy(result, 0); + + const int num_partitions = pm->DefaultNumPartitions(); + + for (int p = 0; p < num_partitions; p++) { + auto &md = pm->mesh_data.GetOrAdd("base", p); + + const auto x_var_pack_string = x_var_type == VarType::Var + ? std::vector{x_var_name_} + : std::vector{}; + const auto x_var = md->PackVariables(x_var_pack_string); + + const auto y_var_pack_string = y_var_type == VarType::Var + ? std::vector{y_var_name_} + : std::vector{}; + const auto y_var = md->PackVariables(y_var_pack_string); + + const auto binned_var_pack_string = binned_var_component == -1 + ? std::vector{} + : std::vector{binned_var_name_}; + const auto binned_var = md->PackVariables(binned_var_pack_string); + + const auto weight_var_pack_string = weight_var_component == -1 + ? std::vector{} + : std::vector{weight_var_name_}; + const auto weight_var = md->PackVariables(weight_var_pack_string); + + const auto ib = md->GetBoundsI(IndexDomain::interior); + const auto jb = md->GetBoundsJ(IndexDomain::interior); + const auto kb = md->GetBoundsK(IndexDomain::interior); + + parthenon::par_for( + DEFAULT_LOOP_PATTERN, "CalcHist", DevExecSpace(), 0, md->NumBlocks() - 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) { + auto &coords = x_var.GetCoords(b); + auto x_val = std::numeric_limits::quiet_NaN(); + if (x_var_type == VarType::X1) { + x_val = coords.Xc<1>(k, j, i); + } else if (x_var_type == VarType::X2) { + x_val = coords.Xc<2>(k, j, i); + } else if (x_var_type == VarType::X3) { + x_val = coords.Xc<3>(k, j, i); + } else if (x_var_type == VarType::R) { + x_val = Kokkos::sqrt(SQR(coords.Xc<1>(k, j, i)) + SQR(coords.Xc<2>(k, j, i)) + + SQR(coords.Xc<3>(k, j, i))); + } else { + x_val = x_var(b, x_var_component, k, j, i); + } + + int x_bin = -1; + // First handle edge cases explicitly + if (x_val < x_edges(0)) { + if (accumulate) { + x_bin = 0; + } else { + return; + } + } else if (x_val > x_edges(x_edges.extent_int(0) - 1)) { + if (accumulate) { + x_bin = x_edges.extent_int(0) - 2; + } else { + return; + } + // if we're on the rightmost edge, directly set last bin + } else if (x_val == x_edges(x_edges.extent_int(0) - 1)) { + x_bin = x_edges.extent_int(0) - 2; + } else { + // for lin and log directly pick index + if (x_edges_type == EdgeType::Lin) { + x_bin = static_cast((x_val - x_edge_min) / x_edge_dbin); + } else if (x_edges_type == EdgeType::Log) { + x_bin = static_cast((Kokkos::log10(x_val) - x_edge_min) / x_edge_dbin); + // otherwise search + } else { + x_bin = upper_bound(x_edges, x_val) - 1; + } + } + PARTHENON_DEBUG_REQUIRE(x_bin >= 0, "Bin not found"); + + // needs to be zero as for the 1D histogram we need 0 as first index of the 2D + // result array + int y_bin = 0; + if (hist_ndim == 2) { + auto y_val = std::numeric_limits::quiet_NaN(); + if (y_var_type == VarType::X1) { + y_val = coords.Xc<1>(k, j, i); + } else if (y_var_type == VarType::X2) { + y_val = coords.Xc<2>(k, j, i); + } else if (y_var_type == VarType::X3) { + y_val = coords.Xc<3>(k, j, i); + } else if (y_var_type == VarType::R) { + y_val = + Kokkos::sqrt(SQR(coords.Xc<1>(k, j, i)) + SQR(coords.Xc<2>(k, j, i)) + + SQR(coords.Xc<3>(k, j, i))); + } else { + y_val = y_var(b, y_var_component, k, j, i); + } + + y_bin = -1; // reset to impossible value + // First handle edge cases explicitly + if (y_val < y_edges(0)) { + if (accumulate) { + y_bin = 0; + } else { + return; + } + } else if (y_val > y_edges(y_edges.extent_int(0) - 1)) { + if (accumulate) { + y_bin = y_edges.extent_int(0) - 2; + } else { + return; + } + // if we're on the rightmost edge, directly set last bin + } else if (y_val == y_edges(y_edges.extent_int(0) - 1)) { + y_bin = y_edges.extent_int(0) - 2; + } else { + // for lin and log directly pick index + if (y_edges_type == EdgeType::Lin) { + y_bin = static_cast((y_val - y_edge_min) / y_edge_dbin); + } else if (y_edges_type == EdgeType::Log) { + y_bin = + static_cast((Kokkos::log10(y_val) - y_edge_min) / y_edge_dbin); + // otherwise search + } else { + y_bin = upper_bound(y_edges, y_val) - 1; + } + } + PARTHENON_DEBUG_REQUIRE(y_bin >= 0, "Bin not found"); + } + auto res = scatter.access(); + const auto val_to_add = binned_var_component == -1 + ? 1 + : binned_var(b, binned_var_component, k, j, i); + auto weight = weight_by_vol ? coords.CellVolume(k, j, i) : 1.0; + weight *= weight_var_component == -1 + ? 1.0 + : weight_var(b, weight_var_component, k, j, i); + res(y_bin, x_bin) += val_to_add * weight; + }); + // "reduce" results from scatter view to original view. May be a no-op depending on + // backend. + Kokkos::Experimental::contribute(result.KokkosView(), scatter); + } + // Ensure all (implicit) reductions from contribute are done + Kokkos::fence(); // May not be required + + // Now reduce over ranks +#ifdef MPI_PARALLEL + if (Globals::my_rank == 0) { + PARTHENON_MPI_CHECK(MPI_Reduce(MPI_IN_PLACE, result.data(), result.size(), + MPI_PARTHENON_REAL, MPI_SUM, 0, MPI_COMM_WORLD)); + } else { + PARTHENON_MPI_CHECK(MPI_Reduce(result.data(), result.data(), result.size(), + MPI_PARTHENON_REAL, MPI_SUM, 0, MPI_COMM_WORLD)); + } +#endif + Kokkos::Profiling::popRegion(); // Calculate single histogram +} + +} // namespace HistUtil + +//---------------------------------------------------------------------------------------- +//! \fn void HistogramOutput:::SetupHistograms(ParameterInput *pin) +// \brief Process parameter input to setup persistent histograms +HistogramOutput::HistogramOutput(const OutputParameters &op, ParameterInput *pin) + : OutputType(op) { + + hist_names_ = pin->GetVector(op.block_name, "hist_names"); + + for (auto &hist_name : hist_names_) { + histograms_.emplace_back(pin, op.block_name, hist_name); + } +} + +std::string HistogramOutput::GenerateFilename_(ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal) { + using namespace HDF5; + + auto filename = std::string(output_params.file_basename); + filename.append("."); + filename.append(output_params.file_id); + filename.append(".histograms."); + if (signal == SignalHandler::OutputSignal::now) { + filename.append("now"); + } else if (signal == SignalHandler::OutputSignal::final && + output_params.file_label_final) { + filename.append("final"); + // default time based data dump + } else { + std::stringstream file_number; + file_number << std::setw(output_params.file_number_width) << std::setfill('0') + << output_params.file_number; + filename.append(file_number.str()); + } + filename.append(".hdf"); + + return filename; +} + +//---------------------------------------------------------------------------------------- +//! \fn void HistogramOutput:::WriteOutputFile(Mesh *pm) +// \brief Calculate histograms +void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal) { + Kokkos::Profiling::pushRegion("Calculate all histograms"); + for (auto &hist : histograms_) { + hist.CalcHist(pm); + } + Kokkos::Profiling::popRegion(); // Calculate all histograms + + Kokkos::Profiling::pushRegion("Dump histograms"); + // Given the expect size of histograms, we'll use serial HDF + if (Globals::my_rank == 0) { + using namespace HDF5; + H5P const pl_xfer = H5P::FromHIDCheck(H5Pcreate(H5P_DATASET_XFER)); + + // As we're reusing the interface from the existing hdf5 output, we have to define + // everything as 7D arrays. + // Counts will be set for each histogram individually below. + const std::array local_offset({0, 0, 0, 0, 0, 0, 0}); + std::array local_count({0, 0, 0, 0, 0, 0, 0}); + std::array global_count({0, 0, 0, 0, 0, 0, 0}); + + // create/open HDF5 file + const std::string filename = GenerateFilename_(pin, tm, signal); + + H5F file; + try { + file = H5F::FromHIDCheck( + H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)); + } catch (std::exception &ex) { + std::stringstream err; + err << "### ERROR: Failed to create HDF5 output file '" << filename + << "' with the following error:" << std::endl + << ex.what() << std::endl; + PARTHENON_THROW(err) + } + + const H5G info_group = MakeGroup(file, "/Info"); + if (tm != nullptr) { + HDF5WriteAttribute("NCycle", tm->ncycle, info_group); + HDF5WriteAttribute("Time", tm->time, info_group); + HDF5WriteAttribute("dt", tm->dt, info_group); + } + HDF5WriteAttribute("hist_names", hist_names_, info_group); + + for (auto &hist : histograms_) { + const H5G hist_group = MakeGroup(file, "/" + hist.name_); + HDF5WriteAttribute("ndim", hist.ndim_, hist_group); + HDF5WriteAttribute("x_var_name", hist.x_var_name_.c_str(), hist_group); + HDF5WriteAttribute("x_var_component", hist.x_var_component_, hist_group); + HDF5WriteAttribute("binned_var_name", hist.binned_var_name_.c_str(), hist_group); + HDF5WriteAttribute("binned_var_component", hist.binned_var_component_, hist_group); + + const auto x_edges_h = hist.x_edges_.GetHostMirrorAndCopy(); + local_count[0] = global_count[0] = x_edges_h.extent_int(0); + HDF5Write1D(hist_group, "x_edges", x_edges_h.data(), local_offset.data(), + local_count.data(), global_count.data(), pl_xfer); + + if (hist.ndim_ == 2) { + HDF5WriteAttribute("y_var_name", hist.y_var_name_.c_str(), hist_group); + HDF5WriteAttribute("y_var_component", hist.y_var_component_, hist_group); + + const auto y_edges_h = hist.y_edges_.GetHostMirrorAndCopy(); + local_count[0] = global_count[0] = y_edges_h.extent_int(0); + HDF5Write1D(hist_group, "y_edges", y_edges_h.data(), local_offset.data(), + local_count.data(), global_count.data(), pl_xfer); + } + + const auto hist_h = hist.result_.GetHostMirrorAndCopy(); + // Ensure correct output format (as the data in Parthenon may, in theory, vary by + // changing the default view layout) so that it matches the numpy output (row + // major, x first) + std::vector tmp_data(hist_h.size()); + int idx = 0; + for (int i = 0; i < hist_h.extent_int(1); ++i) { + for (int j = 0; j < hist_h.extent_int(0); ++j) { + tmp_data[idx++] = hist_h(j, i); + } + } + + local_count[0] = global_count[0] = hist_h.extent_int(1); + if (hist.ndim_ == 2) { + local_count[1] = global_count[1] = hist_h.extent_int(0); + + HDF5Write2D(hist_group, "data", tmp_data.data(), local_offset.data(), + local_count.data(), global_count.data(), pl_xfer); + } else { + // No y-dim for 1D histogram -- though unnecessary as it's not read anyway + local_count[1] = global_count[1] = 0; + HDF5Write1D(hist_group, "data", tmp_data.data(), local_offset.data(), + local_count.data(), global_count.data(), pl_xfer); + } + } + } + Kokkos::Profiling::popRegion(); // Dump histograms + + // advance file ids + if (signal == SignalHandler::OutputSignal::none) { + // After file has been opened with the current number, already advance output + // parameters so that for restarts the file is not immediatly overwritten again. + // Only applies to default time-based data dumps, so that writing "now" and "final" + // outputs does not change the desired output numbering. + output_params.file_number++; + output_params.next_time += output_params.dt; + pin->SetInteger(output_params.block_name, "file_number", output_params.file_number); + pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + } +} + +} // namespace parthenon +#endif // ifndef PARTHENON_ENABLE_ASCENT diff --git a/src/outputs/outputs.cpp b/src/outputs/outputs.cpp index de54d7827c59..33b35a24f284 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -198,7 +198,7 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { // set output variable and optional data format string used in formatted writes if ((op.file_type != "hst") && (op.file_type != "rst") && - (op.file_type != "ascent")) { + (op.file_type != "ascent") && (op.file_type != "histogram")) { op.variables = pin->GetOrAddVector(pib->block_name, "variables", std::vector()); // JMM: If the requested var isn't present for a given swarm, @@ -246,6 +246,17 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { pnew_type = new VTKOutput(op); } else if (op.file_type == "ascent") { pnew_type = new AscentOutput(op); + } else if (op.file_type == "histogram") { +#ifdef ENABLE_HDF5 + pnew_type = new HistogramOutput(op, pin); +#else + msg << "### FATAL ERROR in Outputs constructor" << std::endl + << "Executable not configured for HDF5 outputs, but HDF5 file format " + << "is requested in output/restart block '" << op.block_name << "'. " + << "You can disable this block without deleting it by setting a dt < 0." + << std::endl; + PARTHENON_FAIL(msg); +#endif // ifdef ENABLE_HDF5 } else if (is_hdf5_output) { const bool restart = (op.file_type == "rst"); if (restart) { diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 08c6676da8b4..4fd64236a19c 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -25,10 +25,13 @@ #include #include +#include "Kokkos_ScatterView.hpp" + #include "basic_types.hpp" #include "coordinates/coordinates.hpp" #include "interface/mesh_data.hpp" #include "io_wrapper.hpp" +#include "kokkos_abstraction.hpp" #include "parthenon_arrays.hpp" #include "utils/error_checking.hpp" @@ -214,6 +217,60 @@ class PHDF5Output : public OutputType { const IndexRange &kb, std::vector &x, std::vector &y, std::vector &z); }; + +//---------------------------------------------------------------------------------------- +//! \class HistogramOutput +// \brief derived OutputType class for histograms + +namespace HistUtil { + +enum class VarType { X1, X2, X3, R, Var, Unused }; +enum class EdgeType { Lin, Log, List, Undefined }; + +struct Histogram { + std::string name_; // name (id) of histogram + int ndim_; // 1D or 2D histogram + std::string x_var_name_, y_var_name_; // variable(s) for bins + VarType x_var_type_, y_var_type_; // type, e.g., coord related or actual field + int x_var_component_, y_var_component_; // components of bin variables (vector) + ParArray1D x_edges_, y_edges_; + EdgeType x_edges_type_, y_edges_type_; + // Lowest edge and difference between edges. + // Internally used to speed up lookup for log (and lin) bins as otherwise + // two more log10 calls would be required per index. + Real x_edge_min_, x_edge_dbin_, y_edge_min_, y_edge_dbin_; + bool accumulate_; // accumulate data outside binning range in outermost bins + std::string binned_var_name_; // variable name of variable to be binned + // component of variable to be binned. If -1 means no variable is binned but the + // histgram is a sample count. + int binned_var_component_; + bool weight_by_vol_; // use volume weighting + std::string weight_var_name_; // variable name of variable used as weight + // component of variable to be used as weight. If -1 means no weighting + int weight_var_component_; + ParArray2D result_; // resulting histogram + + // temp view for histogram reduction for better performance (switches + // between atomics and data duplication depending on the platform) + Kokkos::Experimental::ScatterView scatter_result; + Histogram(ParameterInput *pin, const std::string &block_name, const std::string &name); + void CalcHist(Mesh *pm); +}; + +} // namespace HistUtil + +class HistogramOutput : public OutputType { + public: + HistogramOutput(const OutputParameters &oparams, ParameterInput *pin); + void WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal) override; + + private: + std::string GenerateFilename_(ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal); + std::vector hist_names_; // names (used as id) for different histograms + std::vector histograms_; +}; #endif // ifdef ENABLE_HDF5 //---------------------------------------------------------------------------------------- diff --git a/src/utils/signal_handler.cpp b/src/utils/signal_handler.cpp index 940721e15433..7b54540f13fc 100644 --- a/src/utils/signal_handler.cpp +++ b/src/utils/signal_handler.cpp @@ -25,6 +25,9 @@ #include #include +#include FS_HEADER +namespace fs = FS_NAMESPACE; + #include "parthenon_mpi.hpp" #include "globals.hpp" @@ -61,10 +64,8 @@ void SignalHandlerInit() { OutputSignal CheckSignalFlags() { if (Globals::my_rank == 0) { - // TODO(the person bumping std to C++17): use std::filesystem::exists - struct stat buffer; // if file "output_now" exists - if (stat("output_now", &buffer) == 0) { + if (fs::exists("output_now")) { signalflag[nsignal] = 1; } } diff --git a/src/utils/sort.hpp b/src/utils/sort.hpp index 802cd6f56f52..97e9c77a88e4 100644 --- a/src/utils/sort.hpp +++ b/src/utils/sort.hpp @@ -31,6 +31,31 @@ namespace parthenon { +// Returns the upper bound (or the array size if value has not been found) +// Could/Should be replaced with a Kokkos std version once available (currently schedule +// for 4.2 release). +// Note, the API follows the std::upper_bound with the difference of taking an +// array/view as input rather than first and last Iterators, and returning an index +// rather than an Iterator. +template +KOKKOS_INLINE_FUNCTION int upper_bound(const T &arr, Real val) { + int l = 0; + int r = arr.extent_int(0); + int m; + while (l < r) { + m = l + (r - l) / 2; + if (val >= arr(m)) { + l = m + 1; + } else { + r = m; + } + } + if (l < arr.extent_int(0) && val >= arr(l)) { + l++; + } + return l; +} + template void sort(ParArray1D data, KeyComparator comparator, size_t min_idx, size_t max_idx) { diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index a11d41d6a1ff..efe7eca5b6ba 100644 --- a/tst/regression/test_suites/output_hdf5/output_hdf5.py +++ b/tst/regression/test_suites/output_hdf5/output_hdf5.py @@ -1,6 +1,6 @@ # ======================================================================================== # Parthenon performance portable AMR framework -# Copyright(C) 2020 The Parthenon collaboration +# Copyright(C) 2020-2023 The Parthenon collaboration # Licensed under the 3-clause BSD License, see LICENSE file for details # ======================================================================================== # (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved. @@ -20,6 +20,7 @@ import sys import os import utils.test_case +import h5py # To prevent littering up imported folders with .pyc files or __pycache_ folder sys.dont_write_bytecode = True @@ -92,8 +93,9 @@ def Analyse(self, parameters): try: import phdf_diff + import phdf except ModuleNotFoundError: - print("Couldn't find module to compare Parthenon hdf5 files.") + print("Couldn't find modules to read/compare Parthenon hdf5 files.") return False # TODO(pgrete) make sure this also works/doesn't fail for the user @@ -166,4 +168,75 @@ def Analyse(self, parameters): ) analyze_status = False + # Checking Parthenon histograms versus numpy ones + for dim in [2, 3]: + # 1D histogram with binning of a variable with bins defined by a var + data = phdf.phdf(f"advection_{dim}d.out0.final.phdf") + advected = data.Get("advected") + hist_np1d = np.histogram( + advected, [1e-9, 1e-4, 1e-1, 2e-1, 5e-1, 1e0], weights=advected + ) + with h5py.File( + f"advection_{dim}d.out2.histograms.final.hdf", "r" + ) as infile: + hist_parth = infile["hist0/data"][:] + all_close = np.allclose(hist_parth, hist_np1d[0]) + if not all_close: + print(f"1D variable-based hist for {dim}D setup don't match") + analyze_status = False + + # 2D histogram with binning of a variable with bins defined by one var and one coord + omadvected = data.Get("one_minus_advected_sq") + z, y, x = data.GetVolumeLocations() + hist_np2d = np.histogram2d( + x.flatten(), + omadvected.flatten(), + [[-0.5, -0.25, 0, 0.25, 0.5], [0, 0.5, 1]], + weights=advected.flatten(), + ) + with h5py.File( + f"advection_{dim}d.out2.histograms.final.hdf", "r" + ) as infile: + hist_parth = infile["name/data"][:] + # testing slices separately to ensure matching numpy convention + all_close = np.allclose(hist_parth[:, 0], hist_np2d[0][:, 0]) + all_close &= np.allclose(hist_parth[:, 1], hist_np2d[0][:, 1]) + if not all_close: + print(f"2D hist for {dim}D setup don't match") + analyze_status = False + + # 1D histogram (simple sampling) with bins defined by a var + hist_np1d = np.histogram(advected, np.logspace(-9, 0, 11, endpoint=True)) + with h5py.File( + f"advection_{dim}d.out2.histograms.final.hdf", "r" + ) as infile: + hist_parth = infile["other_name/data"][:] + all_close = np.allclose(hist_parth, hist_np1d[0]) + if not all_close: + print(f"1D sampling-based hist for {dim}D setup don't match") + analyze_status = False + + # 2D histogram with volume weighted binning of a variable with bins defined by coords + vols = np.einsum( + "ai,aj,ak->aijk", np.diff(data.zf), np.diff(data.yf), np.diff(data.xf) + ) + hist_np2d = np.histogram2d( + x.flatten(), + y.flatten(), + [[-0.5, -0.25, 0, 0.25, 0.5], [-0.5, -0.1, 0, 0.1, 0.5]], + weights=advected.flatten() * vols.flatten() * omadvected.flatten(), + ) + with h5py.File( + f"advection_{dim}d.out3.histograms.final.hdf", "r" + ) as infile: + hist_parth = infile["hist0/data"][:] + # testing slices separately to ensure matching numpy convention + all_close = np.allclose(hist_parth[:, 0], hist_np2d[0][:, 0]) + all_close &= np.allclose(hist_parth[:, 1], hist_np2d[0][:, 1]) + all_close &= np.allclose(hist_parth[:, 2], hist_np2d[0][:, 2]) + all_close &= np.allclose(hist_parth[:, 3], hist_np2d[0][:, 3]) + if not all_close: + print(f"2D vol-weighted hist for {dim}D setup don't match") + analyze_status = False + return analyze_status diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index e4ba9618c069..3061d2bff961 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -69,3 +69,68 @@ variables = advected, one_minus_advected, & # comments are ok file_type = hst dt = 0.25 + + +file_type = histogram +dt = 0.25 + +hist_names = hist0, name, other_name + +# 1D histogram of a variable, binned by a variable +hist0_ndim = 1 +hist0_x_variable = advected +hist0_x_variable_component = 0 +hist0_x_edges_type = list +hist0_x_edges_list = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 +hist0_binned_variable = advected +hist0_binned_variable_component = 0 + +# 2D histogram of a variable, binned by a coordinate and a different variable +name_ndim = 2 +name_x_variable = HIST_COORD_X1 +name_x_edges_type = lin +name_x_edges_num_bins = 4 +name_x_edges_min = -0.5 +name_x_edges_max = 0.5 +name_y_variable = one_minus_advected_sq +name_y_variable_component = 0 +name_y_edges_type = list +name_y_edges_list = 0, 0.5, 1.0 +name_binned_variable = advected +name_binned_variable_component = 0 + +# 1D histogram ("standard", i.e., counting occurance in bin) +other_name_ndim = 1 +other_name_x_variable = advected +other_name_x_variable_component = 0 +other_name_x_edges_type = log +other_name_x_edges_num_bins = 10 +other_name_x_edges_min = 1e-9 +other_name_x_edges_max = 1e0 +other_name_binned_variable = HIST_ONES + +# A second output block with different dt for histograms +# to double check that writing to different files works + +file_type = histogram +dt = 0.5 + +hist_names = hist0 + +# 2D histogram of volume weighted variable according to two coordinates +hist0_ndim = 2 +hist0_x_variable = HIST_COORD_X1 +hist0_x_edges_type = list +# Note that the coordinate edges are smaller than the domain extents on purpose +# to test the accumulation feature (as the reference histogram in the test calculated +# with numpy goes all the way out to the domain edges). +hist0_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.35 +hist0_y_variable = HIST_COORD_X2 +hist0_y_edges_type = list +hist0_y_edges_list = -0.25, -0.1, 0.0, 0.1, 0.5 +hist0_binned_variable = advected +hist0_binned_variable_component = 0 +hist0_weight_by_volume = true +hist0_weight_variable = one_minus_advected_sq +hist0_weight_variable_component = 0 +hist0_accumulate = true diff --git a/tst/unit/CMakeLists.txt b/tst/unit/CMakeLists.txt index 9efbdff6215f..95f7ca3ecebd 100644 --- a/tst/unit/CMakeLists.txt +++ b/tst/unit/CMakeLists.txt @@ -39,6 +39,7 @@ list(APPEND unit_tests_SOURCES test_partitioning.cpp test_state_descriptor.cpp test_unit_integrators.cpp + test_upper_bound.cpp ) add_executable(unit_tests "${unit_tests_SOURCES}") diff --git a/tst/unit/test_upper_bound.cpp b/tst/unit/test_upper_bound.cpp new file mode 100644 index 000000000000..20244ee95089 --- /dev/null +++ b/tst/unit/test_upper_bound.cpp @@ -0,0 +1,99 @@ +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2023 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== + +#include +#include + +#include + +#include "Kokkos_Core.hpp" +#include "utils/sort.hpp" + +TEST_CASE("upper_bound", "[between][out of bounds][on edges]") { + GIVEN("A sorted list") { + const std::vector data{-1, 0, 1e-2, 5, 10}; + + Kokkos::View arr("arr", data.size()); + auto arr_h = Kokkos::create_mirror_view(arr); + + for (int i = 0; i < data.size(); i++) { + arr_h(i) = data[i]; + } + + Kokkos::deep_copy(arr, arr_h); + + WHEN("a value between entries is given") { + int result; + double val = 0.001; + Kokkos::parallel_reduce( + "unit::upper_bound::between", 1, + KOKKOS_LAMBDA(int /*i*/, int &lres) { + lres = parthenon::upper_bound(arr, val); + }, + result); + THEN("then the next index is returned") { REQUIRE(result == 2); } + THEN("it matches the stl result") { + REQUIRE(result == std::upper_bound(data.begin(), data.end(), val) - data.begin()); + } + } + WHEN("a value below the lower bound is given") { + int result; + double val = -1.1; + Kokkos::parallel_reduce( + "unit::upper_bound::below", 1, + KOKKOS_LAMBDA(int /*i*/, int &lres) { + lres = parthenon::upper_bound(arr, val); + }, + result); + THEN("then the first index is returned") { REQUIRE(result == 0); } + THEN("it matches the stl result") { + REQUIRE(result == std::upper_bound(data.begin(), data.end(), val) - data.begin()); + } + } + WHEN("a value above the upper bound is given") { + int result; + double val = 10.01; + Kokkos::parallel_reduce( + "unit::upper_bound::above", 1, + KOKKOS_LAMBDA(int /*i*/, int &lres) { + lres = parthenon::upper_bound(arr, val); + }, + result); + THEN("then the length of the array is returned") { REQUIRE(result == data.size()); } + THEN("it matches the stl result") { + REQUIRE(result == std::upper_bound(data.begin(), data.end(), val) - data.begin()); + } + } + WHEN("a value on the left edge is given") { + int result; + double val = -1; + Kokkos::parallel_reduce( + "unit::upper_bound::left", 1, + KOKKOS_LAMBDA(int /*i*/, int &lres) { + lres = parthenon::upper_bound(arr, val); + }, + result); + THEN("then the second index is returned") { REQUIRE(result == 1); } + THEN("it matches the stl result") { + REQUIRE(result == std::upper_bound(data.begin(), data.end(), val) - data.begin()); + } + } + WHEN("a value on the right edge is given") { + int result; + double val = 10; + Kokkos::parallel_reduce( + "unit::upper_bound::right", 1, + KOKKOS_LAMBDA(int /*i*/, int &lres) { + lres = parthenon::upper_bound(arr, val); + }, + result); + THEN("then the length of the array is returned") { REQUIRE(result == data.size()); } + THEN("it matches the stl result") { + REQUIRE(result == std::upper_bound(data.begin(), data.end(), val) - data.begin()); + } + } + } +}