From 3f3226eee642b3820d30b80593494ad1e80976e4 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 12 Oct 2023 21:39:26 +0200 Subject: [PATCH 01/22] Add output skeleton for histogram outputs --- src/CMakeLists.txt | 1 + src/outputs/outputs.cpp | 13 ++++++++++++- src/outputs/outputs.hpp | 11 +++++++++++ .../test_suites/output_hdf5/parthinput.advection | 5 +++++ 4 files changed, 29 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 54a226c553831..567b3ec800f7c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -172,6 +172,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/outputs/outputs.cpp b/src/outputs/outputs.cpp index de54d7827c59e..f1084d4a360d1 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); +#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 08c6676da8b4e..66cffc56bdb3e 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -214,6 +214,17 @@ class PHDF5Output : public OutputType { const IndexRange &kb, std::vector &x, std::vector &y, std::vector &z); }; + +//---------------------------------------------------------------------------------------- +//! \class HistogramOutput +// \brief derived OutputType class for histograms + +class HistogramOutput : public OutputType { + public: + explicit HistogramOutput(const OutputParameters &oparams) : OutputType(oparams) {} + void WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal) override; +}; #endif // ifdef ENABLE_HDF5 //---------------------------------------------------------------------------------------- diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index e4ba9618c069b..42472be2d012d 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -69,3 +69,8 @@ variables = advected, one_minus_advected, & # comments are ok file_type = hst dt = 0.25 + + +file_type = histogram +dt = 0.25 + From 67434b3ab6482b1945833de401dc485fc282e540 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 12 Oct 2023 23:23:07 +0200 Subject: [PATCH 02/22] histogram input processing --- src/outputs/histogram.cpp | 156 ++++++++++++++++++ src/outputs/outputs.cpp | 2 +- src/outputs/outputs.hpp | 5 +- .../output_hdf5/parthinput.advection | 4 + 4 files changed, 165 insertions(+), 2 deletions(-) create mode 100644 src/outputs/histogram.cpp diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp new file mode 100644 index 0000000000000..0791f4133a2aa --- /dev/null +++ b/src/outputs/histogram.cpp @@ -0,0 +1,156 @@ +//======================================================================================== +// 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 "config.hpp" +#include "globals.hpp" +#include "kokkos_abstraction.hpp" +#include "parameter_input.hpp" +#include "utils/error_checking.hpp" +#include +#include +#include + +// Only proceed if HDF5 output enabled +#ifdef ENABLE_HDF5 + +#include +#include +#include +#include +#include +#include +#include + +// Parthenon headers +#include "coordinates/coordinates.hpp" +#include "defs.hpp" +#include "globals.hpp" +#include "interface/variable_state.hpp" +#include "mesh/mesh.hpp" +#include "outputs/output_utils.hpp" +#include "outputs/outputs.hpp" +#include "utils/error_checking.hpp" + +// Ascent headers +#ifdef PARTHENON_ENABLE_ASCENT +#include "ascent.hpp" +#include "conduit_blueprint.hpp" +#include "conduit_relay_io.hpp" +#include "conduit_relay_io_blueprint.hpp" +#endif // ifdef PARTHENON_ENABLE_ASCENT + +namespace parthenon { + +using namespace OutputUtils; + +namespace HistUtil { + +struct Histogram { + int ndim; // 1D or 2D histogram + std::array bin_var_names; // variable(s) for bins + std::array bin_var_components; // components of bin variables (vector) + ParArray2D bin_edges; + std::string val_var_name; // variable name of variable to be binned + int val_var_component; // component of variable to be binned + ParArray2D hist; // resulting histogram + + Histogram(ParameterInput *pin, const std::string & block_name, const std::string & prefix) { + ndim = pin->GetInteger(block_name, prefix + "ndim"); + PARTHENON_REQUIRE_THROWS(ndim == 1 || ndim == 2, "Histogram dim must be '1' or '2'"); + + const auto x_var_name = pin->GetString(block_name, prefix + "x_variable"); + const auto x_var_component = + pin->GetInteger(block_name, prefix + "x_variable_component"); + const auto x_edges = pin->GetVector(block_name, prefix + "x_edges"); + + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(x_var_component >= 0, + "Negative component indices are not supported"); + // required by binning index function + PARTHENON_REQUIRE_THROWS(std::is_sorted(x_edges.begin(), x_edges.end()), + "Bin edges must be in order."); + + // For 1D profile default initalize y variables + std::string y_var_name = ""; + int y_var_component = -1; + auto y_edges = std::vector(); + // and for 2D profile check if they're explicitly set (not default value) + if (ndim == 2) { + y_var_name = pin->GetString(block_name, prefix + "y_variable"); + y_var_component = pin->GetInteger(block_name, prefix + "y_variable_component"); + y_edges = pin->GetVector(block_name, prefix + "y_edges"); + + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(y_var_component >= 0, + "Negative component indices are not supported"); + // required by binning index function + PARTHENON_REQUIRE_THROWS(std::is_sorted(y_edges.begin(), y_edges.end()), + "Bin edges must be in order."); + } + + bin_var_names = {x_var_name, y_var_name}; + bin_var_components = {x_var_component, y_var_component}; + + bin_edges = ParArray2D(prefix + "bin_edges", 2); // TODO split these... + + + val_var_name = pin->GetString(block_name, prefix + "val_variable"); + val_var_component = + pin->GetInteger(block_name, prefix + "val_variable_component"); + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(val_var_component >= 0, + "Negative component indices are not supported"); + + } +}; + +} // 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) { + + num_histograms_ = pin->GetOrAddInteger(op.block_name, "num_histograms", 0); + + std::vector histograms_; // TODO make private class member + + for (int i = 0; i < num_histograms_; i++) { + const auto prefix = "hist" + std::to_string(i) + "_"; + histograms_.emplace_back(pin, op.block_name, prefix); + } +} + +//---------------------------------------------------------------------------------------- +//! \fn void HistogramOutput:::WriteOutputFile(Mesh *pm) +// \brief Calculate histograms +void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal) { + + // advance output parameters + 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 f1084d4a360d1..33b35a24f284f 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -248,7 +248,7 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { pnew_type = new AscentOutput(op); } else if (op.file_type == "histogram") { #ifdef ENABLE_HDF5 - pnew_type = new HistogramOutput(op); + 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 " diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 66cffc56bdb3e..fa949d2cd51a5 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -221,9 +221,12 @@ class PHDF5Output : public OutputType { class HistogramOutput : public OutputType { public: - explicit HistogramOutput(const OutputParameters &oparams) : OutputType(oparams) {} + HistogramOutput(const OutputParameters &oparams, ParameterInput *pin); void WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, const SignalHandler::OutputSignal signal) override; + + private: + int num_histograms_; // number of different histograms to compute }; #endif // ifdef ENABLE_HDF5 diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index 42472be2d012d..0f6161606a7fa 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -74,3 +74,7 @@ dt = 0.25 file_type = histogram dt = 0.25 +num_histograms = 1 + + + From 6cfd0bbf81cd708a8fc1896a0395996eaeb4939a Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Fri, 13 Oct 2023 16:17:26 +0200 Subject: [PATCH 03/22] Add CalcHist function --- src/outputs/histogram.cpp | 165 +++++++++++++++++++++++++++++++------- 1 file changed, 137 insertions(+), 28 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 0791f4133a2aa..2183c1e68b6d8 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -18,10 +18,12 @@ // \brief 1D and 2D histograms // options for building +#include "Kokkos_ScatterView.hpp" #include "config.hpp" #include "globals.hpp" #include "kokkos_abstraction.hpp" #include "parameter_input.hpp" +#include "parthenon_array_generic.hpp" #include "utils/error_checking.hpp" #include #include @@ -63,63 +65,170 @@ using namespace OutputUtils; namespace HistUtil { struct Histogram { - int ndim; // 1D or 2D histogram - std::array bin_var_names; // variable(s) for bins - std::array bin_var_components; // components of bin variables (vector) - ParArray2D bin_edges; - std::string val_var_name; // variable name of variable to be binned - int val_var_component; // component of variable to be binned - ParArray2D hist; // resulting histogram - - Histogram(ParameterInput *pin, const std::string & block_name, const std::string & prefix) { + int ndim; // 1D or 2D histogram + std::string x_var_name, y_var_name; // variable(s) for bins + int x_var_component, y_var_component; // components of bin variables (vector) + ParArray1D x_edges, y_edges; + std::string binned_var_name; // variable name of variable to be binned + int binned_var_component; // component of variable to be binned + 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 &prefix) { ndim = pin->GetInteger(block_name, prefix + "ndim"); PARTHENON_REQUIRE_THROWS(ndim == 1 || ndim == 2, "Histogram dim must be '1' or '2'"); - const auto x_var_name = pin->GetString(block_name, prefix + "x_variable"); - const auto x_var_component = - pin->GetInteger(block_name, prefix + "x_variable_component"); - const auto x_edges = pin->GetVector(block_name, prefix + "x_edges"); + x_var_name = pin->GetString(block_name, prefix + "x_variable"); + x_var_component = pin->GetInteger(block_name, prefix + "x_variable_component"); // would add additional logic to pick it from a pack... PARTHENON_REQUIRE_THROWS(x_var_component >= 0, "Negative component indices are not supported"); + + const auto x_edges_in = pin->GetVector(block_name, prefix + "x_edges"); // required by binning index function - PARTHENON_REQUIRE_THROWS(std::is_sorted(x_edges.begin(), x_edges.end()), + PARTHENON_REQUIRE_THROWS(std::is_sorted(x_edges_in.begin(), x_edges_in.end()), "Bin edges must be in order."); + PARTHENON_REQUIRE_THROWS(x_edges_in.size() >= 2, + "Need at least one bin, i.e., two edges."); + x_edges = ParArray1D(prefix + "x_edges", x_edges_in.size()); + auto x_edges_h = x_edges.GetHostMirror(); + for (int i = 0; i < x_edges_in.size(); i++) { + x_edges_h(i) = x_edges_in[i]; + } + Kokkos::deep_copy(x_edges, x_edges_h); // For 1D profile default initalize y variables std::string y_var_name = ""; int y_var_component = -1; - auto y_edges = std::vector(); // and for 2D profile check if they're explicitly set (not default value) if (ndim == 2) { y_var_name = pin->GetString(block_name, prefix + "y_variable"); - y_var_component = pin->GetInteger(block_name, prefix + "y_variable_component"); - y_edges = pin->GetVector(block_name, prefix + "y_edges"); + y_var_component = pin->GetInteger(block_name, prefix + "y_variable_component"); // would add additional logic to pick it from a pack... PARTHENON_REQUIRE_THROWS(y_var_component >= 0, "Negative component indices are not supported"); + + const auto y_edges_in = pin->GetVector(block_name, prefix + "y_edges"); // required by binning index function - PARTHENON_REQUIRE_THROWS(std::is_sorted(y_edges.begin(), y_edges.end()), + PARTHENON_REQUIRE_THROWS(std::is_sorted(y_edges_in.begin(), y_edges_in.end()), "Bin edges must be in order."); + PARTHENON_REQUIRE_THROWS(y_edges_in.size() >= 2, + "Need at least one bin, i.e., two edges."); + y_edges = ParArray1D(prefix + "y_edges", y_edges_in.size()); + auto y_edges_h = y_edges.GetHostMirror(); + for (int i = 0; i < y_edges_in.size(); i++) { + y_edges_h(i) = y_edges_in[i]; + } + Kokkos::deep_copy(y_edges, y_edges_h); + } else { + y_edges = ParArray1D(prefix + "y_edges_unused", 0); } - bin_var_names = {x_var_name, y_var_name}; - bin_var_components = {x_var_component, y_var_component}; + binned_var_name = pin->GetString(block_name, prefix + "binned_variable"); + 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"); - bin_edges = ParArray2D(prefix + "bin_edges", 2); // TODO split these... + 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()); + } +}; - val_var_name = pin->GetString(block_name, prefix + "val_variable"); - val_var_component = - pin->GetInteger(block_name, prefix + "val_variable_component"); - // would add additional logic to pick it from a pack... - PARTHENON_REQUIRE_THROWS(val_var_component >= 0, - "Negative component indices are not supported"); +// Returns the lower 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). +// TODO add unit test +KOKKOS_INLINE_FUNCTION int lower_bound(const ParArray1D &arr, Real val) { + int l = 0; + int r = arr.GetDim(0); + int m; + while (l < r) { + m = l + (r - l) / 2; + if (val <= arr(m)) { + r = m; + } else { + l = m + 1; + } + } + return l; +} + +// Computes a 1D or 2D histogram with inclusive lower edges (and exclusive right ones). +// 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 CalcHist(Mesh *pm, const Histogram &hist) { + const auto x_var_component = hist.x_var_component; + const auto y_var_component = hist.y_var_component; + const auto binned_var_component = hist.binned_var_component; + const auto x_edges = hist.x_edges; + const auto y_edges = hist.y_edges; + const auto hist_ndim = hist.ndim; + auto result = hist.result; + auto scatter = hist.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 + result.Reset(); + 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 = md->PackVariables(std::vector{hist.x_var_name}); + const auto y_var = md->PackVariables(std::vector{hist.y_var_name}); + const auto binned_var = + md->PackVariables(std::vector{hist.binned_var_name}); + 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, x_var.GetDim(5) - 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) { + const auto &x_val = x_var(b, x_var_component, k, j, i); + if (x_val < x_edges(0) || x_val >= x_edges(x_edges.extent_int(0))) { + return; + } + // No further check for x_bin required as the preceeding if-statement guarantees + // x_val to fall in one bin. + const auto x_bin = lower_bound(x_edges, x_val); + + int y_bin = 0; + if (hist_ndim == 2) { + const auto &y_val = y_var(b, y_var_component, k, j, i); + if (y_val < y_edges(0) || y_val >= y_edges(y_edges.extent_int(0))) { + return; + } + // No further check for y_bin required as the preceeding if-statement + // guarantees y_val to fall in one bin. + y_bin = lower_bound(y_edges, y_val); + } + auto res = scatter.access(); + res(y_bin, x_bin) += binned_var(b, binned_var_component, k, j, i); + }); + // "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 +} } // namespace HistUtil From a6d718f84954317aad03669127c132b775fa8913 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sat, 14 Oct 2023 18:55:34 +0200 Subject: [PATCH 04/22] Actually calc histograms --- src/outputs/histogram.cpp | 176 +++++++++--------- src/outputs/outputs.hpp | 22 +++ .../output_hdf5/parthinput.advection | 7 +- 3 files changed, 121 insertions(+), 84 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 2183c1e68b6d8..ce71154f7bbeb 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -19,6 +19,7 @@ // options for building #include "Kokkos_ScatterView.hpp" +#include "basic_types.hpp" #include "config.hpp" #include "globals.hpp" #include "kokkos_abstraction.hpp" @@ -64,107 +65,96 @@ using namespace OutputUtils; namespace HistUtil { -struct Histogram { - int ndim; // 1D or 2D histogram - std::string x_var_name, y_var_name; // variable(s) for bins - int x_var_component, y_var_component; // components of bin variables (vector) - ParArray1D x_edges, y_edges; - std::string binned_var_name; // variable name of variable to be binned - int binned_var_component; // component of variable to be binned - 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 &prefix) { - ndim = pin->GetInteger(block_name, prefix + "ndim"); - PARTHENON_REQUIRE_THROWS(ndim == 1 || ndim == 2, "Histogram dim must be '1' or '2'"); +Histogram::Histogram(ParameterInput *pin, const std::string &block_name, + const std::string &prefix) { + ndim = pin->GetInteger(block_name, prefix + "ndim"); + PARTHENON_REQUIRE_THROWS(ndim == 1 || ndim == 2, "Histogram dim must be '1' or '2'"); + + x_var_name = pin->GetString(block_name, prefix + "x_variable"); + + x_var_component = pin->GetInteger(block_name, prefix + "x_variable_component"); + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(x_var_component >= 0, + "Negative component indices are not supported"); + + const auto x_edges_in = pin->GetVector(block_name, prefix + "x_edges"); + // required by binning index function + PARTHENON_REQUIRE_THROWS(std::is_sorted(x_edges_in.begin(), x_edges_in.end()), + "Bin edges must be in order."); + PARTHENON_REQUIRE_THROWS(x_edges_in.size() >= 2, + "Need at least one bin, i.e., two edges."); + x_edges = ParArray1D(prefix + "x_edges", x_edges_in.size()); + auto x_edges_h = x_edges.GetHostMirror(); + for (int i = 0; i < x_edges_in.size(); i++) { + x_edges_h(i) = x_edges_in[i]; + } + Kokkos::deep_copy(x_edges, x_edges_h); - x_var_name = pin->GetString(block_name, prefix + "x_variable"); + // For 1D profile default initalize y variables + y_var_name = ""; + y_var_component = -1; + // and for 2D profile check if they're explicitly set (not default value) + if (ndim == 2) { + y_var_name = pin->GetString(block_name, prefix + "y_variable"); - x_var_component = pin->GetInteger(block_name, prefix + "x_variable_component"); + y_var_component = pin->GetInteger(block_name, prefix + "y_variable_component"); // would add additional logic to pick it from a pack... - PARTHENON_REQUIRE_THROWS(x_var_component >= 0, + PARTHENON_REQUIRE_THROWS(y_var_component >= 0, "Negative component indices are not supported"); - const auto x_edges_in = pin->GetVector(block_name, prefix + "x_edges"); + const auto y_edges_in = pin->GetVector(block_name, prefix + "y_edges"); // required by binning index function - PARTHENON_REQUIRE_THROWS(std::is_sorted(x_edges_in.begin(), x_edges_in.end()), + PARTHENON_REQUIRE_THROWS(std::is_sorted(y_edges_in.begin(), y_edges_in.end()), "Bin edges must be in order."); - PARTHENON_REQUIRE_THROWS(x_edges_in.size() >= 2, + PARTHENON_REQUIRE_THROWS(y_edges_in.size() >= 2, "Need at least one bin, i.e., two edges."); - x_edges = ParArray1D(prefix + "x_edges", x_edges_in.size()); - auto x_edges_h = x_edges.GetHostMirror(); - for (int i = 0; i < x_edges_in.size(); i++) { - x_edges_h(i) = x_edges_in[i]; - } - Kokkos::deep_copy(x_edges, x_edges_h); - - // For 1D profile default initalize y variables - std::string y_var_name = ""; - int y_var_component = -1; - // and for 2D profile check if they're explicitly set (not default value) - if (ndim == 2) { - y_var_name = pin->GetString(block_name, prefix + "y_variable"); - - y_var_component = pin->GetInteger(block_name, prefix + "y_variable_component"); - // would add additional logic to pick it from a pack... - PARTHENON_REQUIRE_THROWS(y_var_component >= 0, - "Negative component indices are not supported"); - - const auto y_edges_in = pin->GetVector(block_name, prefix + "y_edges"); - // required by binning index function - PARTHENON_REQUIRE_THROWS(std::is_sorted(y_edges_in.begin(), y_edges_in.end()), - "Bin edges must be in order."); - PARTHENON_REQUIRE_THROWS(y_edges_in.size() >= 2, - "Need at least one bin, i.e., two edges."); - y_edges = ParArray1D(prefix + "y_edges", y_edges_in.size()); - auto y_edges_h = y_edges.GetHostMirror(); - for (int i = 0; i < y_edges_in.size(); i++) { - y_edges_h(i) = y_edges_in[i]; - } - Kokkos::deep_copy(y_edges, y_edges_h); - } else { - y_edges = ParArray1D(prefix + "y_edges_unused", 0); + y_edges = ParArray1D(prefix + "y_edges", y_edges_in.size()); + auto y_edges_h = y_edges.GetHostMirror(); + for (int i = 0; i < y_edges_in.size(); i++) { + y_edges_h(i) = y_edges_in[i]; } + Kokkos::deep_copy(y_edges, y_edges_h); + } else { + y_edges = ParArray1D(prefix + "y_edges_unused", 0); + } - binned_var_name = pin->GetString(block_name, prefix + "binned_variable"); - 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"); + binned_var_name = pin->GetString(block_name, prefix + "binned_variable"); + 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; + 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()); - } -}; + result = ParArray2D(prefix + "result", nybins, nxbins); + scatter_result = Kokkos::Experimental::ScatterView(result.KokkosView()); +} -// Returns the lower bound (or the array size if value has not been found) +// 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). // TODO add unit test -KOKKOS_INLINE_FUNCTION int lower_bound(const ParArray1D &arr, Real val) { +KOKKOS_INLINE_FUNCTION int upper_bound(const ParArray1D &arr, Real val) { int l = 0; - int r = arr.GetDim(0); + int r = arr.extent_int(0); int m; while (l < r) { m = l + (r - l) / 2; - if (val <= arr(m)) { - r = m; - } else { + if (val >= arr(m)) { l = m + 1; + } else { + r = m; } } + if (l < arr.extent_int(0) && val >= arr(l)) { + l++; + } return l; } -// Computes a 1D or 2D histogram with inclusive lower edges (and exclusive right ones). +// 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 CalcHist(Mesh *pm, const Histogram &hist) { @@ -182,7 +172,7 @@ void CalcHist(Mesh *pm, const Histogram &hist) { // 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 - result.Reset(); + Kokkos::deep_copy(result, 0); const int num_partitions = pm->DefaultNumPartitions(); @@ -202,22 +192,22 @@ void CalcHist(Mesh *pm, const Histogram &hist) { kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { const auto &x_val = x_var(b, x_var_component, k, j, i); - if (x_val < x_edges(0) || x_val >= x_edges(x_edges.extent_int(0))) { + if (x_val < x_edges(0) || x_val > x_edges(x_edges.extent_int(0) - 1)) { return; } // No further check for x_bin required as the preceeding if-statement guarantees // x_val to fall in one bin. - const auto x_bin = lower_bound(x_edges, x_val); + const auto x_bin = upper_bound(x_edges, x_val) - 1; int y_bin = 0; if (hist_ndim == 2) { const auto &y_val = y_var(b, y_var_component, k, j, i); - if (y_val < y_edges(0) || y_val >= y_edges(y_edges.extent_int(0))) { + if (y_val < y_edges(0) || y_val > y_edges(y_edges.extent_int(0) - 1)) { return; } // No further check for y_bin required as the preceeding if-statement // guarantees y_val to fall in one bin. - y_bin = lower_bound(y_edges, y_val); + y_bin = upper_bound(y_edges, y_val) - 1; } auto res = scatter.access(); res(y_bin, x_bin) += binned_var(b, binned_var_component, k, j, i); @@ -228,6 +218,17 @@ void CalcHist(Mesh *pm, const Histogram &hist) { } // 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 } } // namespace HistUtil @@ -240,8 +241,6 @@ HistogramOutput::HistogramOutput(const OutputParameters &op, ParameterInput *pin num_histograms_ = pin->GetOrAddInteger(op.block_name, "num_histograms", 0); - std::vector histograms_; // TODO make private class member - for (int i = 0; i < num_histograms_; i++) { const auto prefix = "hist" + std::to_string(i) + "_"; histograms_.emplace_back(pin, op.block_name, prefix); @@ -253,7 +252,18 @@ HistogramOutput::HistogramOutput(const OutputParameters &op, ParameterInput *pin // \brief Calculate histograms void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, const SignalHandler::OutputSignal signal) { - + for (auto &hist : histograms_) { + CalcHist(pm, hist); + + if (Globals::my_rank == 0) { + const auto hist_h = hist.result.GetHostMirrorAndCopy(); + std::cout << "Hist result: "; + for (int i = 0; i < hist_h.extent_int(1); i++) { + std::cout << hist_h(0, i) << " "; + } + std::cout << "\n"; + } + } // advance output parameters output_params.file_number++; output_params.next_time += output_params.dt; diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index fa949d2cd51a5..cadac050da595 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -25,6 +25,8 @@ #include #include +#include "Kokkos_ScatterView.hpp" + #include "basic_types.hpp" #include "coordinates/coordinates.hpp" #include "interface/mesh_data.hpp" @@ -219,6 +221,25 @@ class PHDF5Output : public OutputType { //! \class HistogramOutput // \brief derived OutputType class for histograms +namespace HistUtil { +struct Histogram { + int ndim; // 1D or 2D histogram + std::string x_var_name, y_var_name; // variable(s) for bins + int x_var_component, y_var_component; // components of bin variables (vector) + ParArray1D x_edges, y_edges; + std::string binned_var_name; // variable name of variable to be binned + int binned_var_component; // component of variable to be binned + 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 &prefix); +}; + +} // namespace HistUtil + class HistogramOutput : public OutputType { public: HistogramOutput(const OutputParameters &oparams, ParameterInput *pin); @@ -227,6 +248,7 @@ class HistogramOutput : public OutputType { private: int num_histograms_; // number of different histograms to compute + std::vector histograms_; }; #endif // ifdef ENABLE_HDF5 diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index 0f6161606a7fa..bdaf6025f8f98 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -76,5 +76,10 @@ dt = 0.25 num_histograms = 1 - +hist0_ndim = 1 +hist0_x_variable = advected +hist0_x_variable_component = 0 +hist0_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.00001 +hist0_binned_variable = advected +hist0_binned_variable_component = 0 From 37eaac2154a5e093d8f00de6d0fffba984868eb0 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sat, 14 Oct 2023 21:46:09 +0200 Subject: [PATCH 05/22] First attempt at dumping to hdf5. existing groups are a problem... --- src/outputs/histogram.cpp | 96 ++++++++++++++++++++++++++++++++------- 1 file changed, 80 insertions(+), 16 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index ce71154f7bbeb..7a4cb9e8ebcf4 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -18,7 +18,6 @@ // \brief 1D and 2D histograms // options for building -#include "Kokkos_ScatterView.hpp" #include "basic_types.hpp" #include "config.hpp" #include "globals.hpp" @@ -26,13 +25,12 @@ #include "parameter_input.hpp" #include "parthenon_array_generic.hpp" #include "utils/error_checking.hpp" -#include -#include -#include // Only proceed if HDF5 output enabled #ifdef ENABLE_HDF5 +#include +#include #include #include #include @@ -40,6 +38,7 @@ #include #include #include +#include // Parthenon headers #include "coordinates/coordinates.hpp" @@ -49,15 +48,14 @@ #include "mesh/mesh.hpp" #include "outputs/output_utils.hpp" #include "outputs/outputs.hpp" +#include "outputs/parthenon_hdf5.hpp" #include "utils/error_checking.hpp" -// Ascent headers -#ifdef PARTHENON_ENABLE_ASCENT -#include "ascent.hpp" -#include "conduit_blueprint.hpp" -#include "conduit_relay_io.hpp" -#include "conduit_relay_io_blueprint.hpp" -#endif // ifdef PARTHENON_ENABLE_ASCENT +// ScatterView is not part of Kokkos core interface +#include "Kokkos_ScatterView.hpp" + +#include FS_HEADER +namespace fs = FS_NAMESPACE; namespace parthenon { @@ -158,6 +156,7 @@ KOKKOS_INLINE_FUNCTION int upper_bound(const ParArray1D &arr, Real val) { // 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 CalcHist(Mesh *pm, const Histogram &hist) { + Kokkos::Profiling::pushRegion("Calculate single histogram"); const auto x_var_component = hist.x_var_component; const auto y_var_component = hist.y_var_component; const auto binned_var_component = hist.binned_var_component; @@ -229,6 +228,7 @@ void CalcHist(Mesh *pm, const Histogram &hist) { MPI_PARTHENON_REAL, MPI_SUM, 0, MPI_COMM_WORLD)); } #endif + Kokkos::Profiling::popRegion(); // Calculate single histogram } } // namespace HistUtil @@ -252,10 +252,66 @@ HistogramOutput::HistogramOutput(const OutputParameters &op, ParameterInput *pin // \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_) { CalcHist(pm, hist); + } + Kokkos::Profiling::popRegion(); // Calculate all histograms + + Kokkos::Profiling::pushRegion("Dump histograms"); + if (Globals::my_rank == 0) { + using namespace HDF5; + // create/open HDF5 file + const std::string filename = "histogram.hdf"; + H5F file; + try { + if (fs::exists(filename)) { + file = H5F::FromHIDCheck(H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT)); + } else { + 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 open/create HDF5 output file '" << filename + << "' with the following error:" << std::endl + << ex.what() << std::endl; + PARTHENON_THROW(err) + } + + std::string out_label; + if (signal == SignalHandler::OutputSignal::now) { + out_label = "now"; + } else if (signal == SignalHandler::OutputSignal::final && + output_params.file_label_final) { + out_label = "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; + out_label = file_number.str(); + } + + const H5G all_hist_group = MakeGroup(file, "/" + out_label); + if (tm != nullptr) { + HDF5WriteAttribute("NCycle", tm->ncycle, all_hist_group); + HDF5WriteAttribute("Time", tm->time, all_hist_group); + HDF5WriteAttribute("dt", tm->dt, all_hist_group); + } + HDF5WriteAttribute("num_histograms", num_histograms_, all_hist_group); + + for (int i = 0; i < num_histograms_; i++) { + auto &hist = histograms_[i]; + const H5G hist_group = MakeGroup(all_hist_group, "/" + std::to_string(i)); + HDF5WriteAttribute("x_var_name", hist.x_var_name, hist_group); + HDF5WriteAttribute("x_var_component", hist.x_var_component, hist_group); + HDF5WriteAttribute("y_var_name", hist.y_var_name, hist_group); + HDF5WriteAttribute("y_var_component", hist.y_var_component, hist_group); + HDF5WriteAttribute("binned_var_name", hist.binned_var_name, hist_group); + HDF5WriteAttribute("binned_var_component", hist.binned_var_component, hist_group); - if (Globals::my_rank == 0) { const auto hist_h = hist.result.GetHostMirrorAndCopy(); std::cout << "Hist result: "; for (int i = 0; i < hist_h.extent_int(1); i++) { @@ -264,11 +320,19 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm std::cout << "\n"; } } + Kokkos::Profiling::popRegion(); // Dump histograms + // advance output parameters - 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); + 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 From b030779fd7c724cfc05a31340561654e6400eaa0 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sun, 15 Oct 2023 13:46:35 +0200 Subject: [PATCH 06/22] Use separate hdf5 files for histogram outputs --- src/outputs/histogram.cpp | 70 +++++++++++++++++++----------------- src/outputs/outputs.hpp | 2 ++ src/utils/signal_handler.cpp | 7 ++-- 3 files changed, 44 insertions(+), 35 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 7a4cb9e8ebcf4..e5ad8683a278c 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -54,9 +54,6 @@ // ScatterView is not part of Kokkos core interface #include "Kokkos_ScatterView.hpp" -#include FS_HEADER -namespace fs = FS_NAMESPACE; - namespace parthenon { using namespace OutputUtils; @@ -247,6 +244,31 @@ HistogramOutput::HistogramOutput(const OutputParameters &op, ParameterInput *pin } } +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 @@ -260,51 +282,35 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm 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; // create/open HDF5 file - const std::string filename = "histogram.hdf"; + const std::string filename = GenerateFilename_(pin, tm, signal); + H5F file; try { - if (fs::exists(filename)) { - file = H5F::FromHIDCheck(H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT)); - } else { - file = H5F::FromHIDCheck( - H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)); - } + 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 open/create HDF5 output file '" << filename + err << "### ERROR: Failed to create HDF5 output file '" << filename << "' with the following error:" << std::endl << ex.what() << std::endl; PARTHENON_THROW(err) } - std::string out_label; - if (signal == SignalHandler::OutputSignal::now) { - out_label = "now"; - } else if (signal == SignalHandler::OutputSignal::final && - output_params.file_label_final) { - out_label = "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; - out_label = file_number.str(); - } - - const H5G all_hist_group = MakeGroup(file, "/" + out_label); + const H5G info_group = MakeGroup(file, "/Info"); if (tm != nullptr) { - HDF5WriteAttribute("NCycle", tm->ncycle, all_hist_group); - HDF5WriteAttribute("Time", tm->time, all_hist_group); - HDF5WriteAttribute("dt", tm->dt, all_hist_group); + HDF5WriteAttribute("NCycle", tm->ncycle, info_group); + HDF5WriteAttribute("Time", tm->time, info_group); + HDF5WriteAttribute("dt", tm->dt, info_group); } - HDF5WriteAttribute("num_histograms", num_histograms_, all_hist_group); + HDF5WriteAttribute("num_histograms", num_histograms_, info_group); for (int i = 0; i < num_histograms_; i++) { auto &hist = histograms_[i]; - const H5G hist_group = MakeGroup(all_hist_group, "/" + std::to_string(i)); + const H5G hist_group = MakeGroup(file, "/" + std::to_string(i)); HDF5WriteAttribute("x_var_name", hist.x_var_name, hist_group); HDF5WriteAttribute("x_var_component", hist.x_var_component, hist_group); HDF5WriteAttribute("y_var_name", hist.y_var_name, hist_group); @@ -322,7 +328,7 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm } Kokkos::Profiling::popRegion(); // Dump histograms - // advance output parameters + // 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. diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index cadac050da595..089eff04c9760 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -247,6 +247,8 @@ class HistogramOutput : public OutputType { const SignalHandler::OutputSignal signal) override; private: + std::string GenerateFilename_(ParameterInput *pin, SimTime *tm, + const SignalHandler::OutputSignal signal); int num_histograms_; // number of different histograms to compute std::vector histograms_; }; diff --git a/src/utils/signal_handler.cpp b/src/utils/signal_handler.cpp index 940721e154332..7b54540f13fcb 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; } } From eacd5cb678e084adbdb67df2d66d7503bbb7673e Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sun, 15 Oct 2023 21:02:11 +0200 Subject: [PATCH 07/22] Fix hdf5 output --- src/outputs/histogram.cpp | 57 +++++++++++++++++++++++++++++++++++---- 1 file changed, 52 insertions(+), 5 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index e5ad8683a278c..8728a758ba0c8 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -285,6 +285,15 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm // 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); @@ -308,17 +317,55 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm } HDF5WriteAttribute("num_histograms", num_histograms_, info_group); - for (int i = 0; i < num_histograms_; i++) { - auto &hist = histograms_[i]; - const H5G hist_group = MakeGroup(file, "/" + std::to_string(i)); + for (int h = 0; h < num_histograms_; h++) { + auto &hist = histograms_[h]; + const H5G hist_group = MakeGroup(file, "/" + std::to_string(h)); + HDF5WriteAttribute("ndim", hist.ndim, hist_group); HDF5WriteAttribute("x_var_name", hist.x_var_name, hist_group); HDF5WriteAttribute("x_var_component", hist.x_var_component, hist_group); - HDF5WriteAttribute("y_var_name", hist.y_var_name, hist_group); - HDF5WriteAttribute("y_var_component", hist.y_var_component, hist_group); HDF5WriteAttribute("binned_var_name", hist.binned_var_name, 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, 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 + local_count[1] = global_count[1] = 0; + HDF5Write2D(hist_group, "data", tmp_data.data(), local_offset.data(), + local_count.data(), global_count.data(), pl_xfer); + } + std::cout << "Hist result: "; for (int i = 0; i < hist_h.extent_int(1); i++) { std::cout << hist_h(0, i) << " "; From 47065bca5d458f8a206c831c24055c6b9bdd046e Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sun, 15 Oct 2023 23:06:46 +0200 Subject: [PATCH 08/22] Move upper_bound to utils and add unit test --- src/outputs/histogram.cpp | 23 +------------------ src/utils/sort.hpp | 22 ++++++++++++++++++ .../output_hdf5/parthinput.advection | 12 +++++++++- tst/unit/CMakeLists.txt | 1 + 4 files changed, 35 insertions(+), 23 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 8728a758ba0c8..4d2d653c58ade 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -50,6 +50,7 @@ #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" @@ -127,28 +128,6 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, scatter_result = Kokkos::Experimental::ScatterView(result.KokkosView()); } -// 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). -// TODO add unit test -KOKKOS_INLINE_FUNCTION int upper_bound(const ParArray1D &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; -} - // 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). diff --git a/src/utils/sort.hpp b/src/utils/sort.hpp index 802cd6f56f52f..aed0deeff473c 100644 --- a/src/utils/sort.hpp +++ b/src/utils/sort.hpp @@ -31,6 +31,28 @@ 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). +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/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index bdaf6025f8f98..52ad9c180752e 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -74,7 +74,7 @@ dt = 0.25 file_type = histogram dt = 0.25 -num_histograms = 1 +num_histograms = 2 hist0_ndim = 1 hist0_x_variable = advected @@ -83,3 +83,13 @@ hist0_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.00001 hist0_binned_variable = advected hist0_binned_variable_component = 0 +hist1_ndim = 2 +hist1_x_variable = advected +hist1_x_variable_component = 0 +hist1_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0000001 +hist1_y_variable = one_minus_advected_sq +hist1_y_variable_component = 0 +hist1_y_edges = 0, 0.5, 1.00001 +hist1_binned_variable = advected +hist1_binned_variable_component = 0 + diff --git a/tst/unit/CMakeLists.txt b/tst/unit/CMakeLists.txt index 9efbdff6215ff..95f7ca3ecebd3 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}") From 32a25d57a008306226657f45ee61580f51c955fd Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Sun, 15 Oct 2023 23:21:16 +0200 Subject: [PATCH 09/22] Fix edge bin calc --- src/outputs/histogram.cpp | 25 ++++++++----------- .../output_hdf5/parthinput.advection | 4 +-- 2 files changed, 13 insertions(+), 16 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 4d2d653c58ade..bcc41c1affb54 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -170,9 +170,11 @@ void CalcHist(Mesh *pm, const Histogram &hist) { if (x_val < x_edges(0) || x_val > x_edges(x_edges.extent_int(0) - 1)) { return; } - // No further check for x_bin required as the preceeding if-statement guarantees - // x_val to fall in one bin. - const auto x_bin = upper_bound(x_edges, x_val) - 1; + + // if we're on the rightmost edge, directly set last bin, otherwise search + const auto x_bin = x_val == x_edges(x_edges.extent_int(0) - 1) + ? x_edges.extent_int(0) - 2 + : upper_bound(x_edges, x_val) - 1; int y_bin = 0; if (hist_ndim == 2) { @@ -180,9 +182,10 @@ void CalcHist(Mesh *pm, const Histogram &hist) { if (y_val < y_edges(0) || y_val > y_edges(y_edges.extent_int(0) - 1)) { return; } - // No further check for y_bin required as the preceeding if-statement - // guarantees y_val to fall in one bin. - y_bin = upper_bound(y_edges, y_val) - 1; + // if we're on the rightmost edge, directly set last bin, otherwise search + const auto y_bin = y_val == y_edges(y_edges.extent_int(0) - 1) + ? y_edges.extent_int(0) - 2 + : upper_bound(y_edges, y_val) - 1; } auto res = scatter.access(); res(y_bin, x_bin) += binned_var(b, binned_var_component, k, j, i); @@ -339,17 +342,11 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm 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 + // No y-dim for 1D histogram -- though unnecessary as it's not read anyway local_count[1] = global_count[1] = 0; - HDF5Write2D(hist_group, "data", tmp_data.data(), local_offset.data(), + HDF5Write1D(hist_group, "data", tmp_data.data(), local_offset.data(), local_count.data(), global_count.data(), pl_xfer); } - - std::cout << "Hist result: "; - for (int i = 0; i < hist_h.extent_int(1); i++) { - std::cout << hist_h(0, i) << " "; - } - std::cout << "\n"; } } Kokkos::Profiling::popRegion(); // Dump histograms diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index 52ad9c180752e..a9e21c6581b6f 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -86,10 +86,10 @@ hist0_binned_variable_component = 0 hist1_ndim = 2 hist1_x_variable = advected hist1_x_variable_component = 0 -hist1_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0000001 +hist1_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.00 hist1_y_variable = one_minus_advected_sq hist1_y_variable_component = 0 -hist1_y_edges = 0, 0.5, 1.00001 +hist1_y_edges = 0, 0.5, 1.0 hist1_binned_variable = advected hist1_binned_variable_component = 0 From c80c146844260c87fe499be7cb9375751c816a9f Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 16 Oct 2023 00:09:27 +0200 Subject: [PATCH 10/22] Add histogram regression test --- tst/regression/CMakeLists.txt | 1 + .../test_suites/output_hdf5/output_hdf5.py | 40 ++++++++++++++++++- .../output_hdf5/parthinput.advection | 2 +- 3 files changed, 40 insertions(+), 3 deletions(-) diff --git a/tst/regression/CMakeLists.txt b/tst/regression/CMakeLists.txt index bf902381bbe79..96ac7c583c678 100644 --- a/tst/regression/CMakeLists.txt +++ b/tst/regression/CMakeLists.txt @@ -127,6 +127,7 @@ endif() # Any external modules that are required by python can be added to REQUIRED_PYTHON_MODULES # list variable, before including TestSetup.cmake. list(APPEND REQUIRED_PYTHON_MODULES numpy) +list(APPEND REQUIRED_PYTHON_MODULES h5py) list(APPEND DESIRED_PYTHON_MODULES matplotlib) # Include test setup functions, and check for python interpreter and modules diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index a11d41d6a1ff1..74f90e61fe3b3 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,38 @@ def Analyse(self, parameters): ) analyze_status = False + # Checking Parthenon histograms versus numpy ones + for dim in [2, 3]: + 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["0/data"][:] + all_close = np.allclose(hist_parth, hist_np1d[0]) + if not all_close: + print(f"1D hist for {dim}D setup don't match") + analyze_status = False + + omadvected = data.Get("one_minus_advected_sq") + hist_np2d = np.histogram2d( + advected.flatten(), + omadvected.flatten(), + [[1e-9, 1e-4, 1e-1, 2e-1, 5e-1, 1e0], [0, 0.5, 1]], + weights=advected.flatten(), + ) + with h5py.File( + f"advection_{dim}d.out2.histograms.final.hdf", "r" + ) as infile: + hist_parth = infile["1/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 + return analyze_status diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index a9e21c6581b6f..812c4ef756050 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -79,7 +79,7 @@ num_histograms = 2 hist0_ndim = 1 hist0_x_variable = advected hist0_x_variable_component = 0 -hist0_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.00001 +hist0_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 hist0_binned_variable = advected hist0_binned_variable_component = 0 From e801cf48e91bc2bdea1bff82073b8894ea42b72f Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 16 Oct 2023 00:12:05 +0200 Subject: [PATCH 11/22] Add missing unit test file --- tst/unit/test_upper_bound.cpp | 100 ++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) create mode 100644 tst/unit/test_upper_bound.cpp diff --git a/tst/unit/test_upper_bound.cpp b/tst/unit/test_upper_bound.cpp new file mode 100644 index 0000000000000..dff514de8e342 --- /dev/null +++ b/tst/unit/test_upper_bound.cpp @@ -0,0 +1,100 @@ +//======================================================================================== +// 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()); + } + } + } +} From 26719323777ec6d20f873aed1410346e3754e953 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 16 Oct 2023 14:10:07 +0200 Subject: [PATCH 12/22] Allow coordinate based binning --- src/outputs/histogram.cpp | 101 ++++++++++++++---- src/outputs/outputs.hpp | 3 + src/utils/sort.hpp | 2 +- .../test_suites/output_hdf5/output_hdf5.py | 7 +- .../output_hdf5/parthinput.advection | 5 +- 5 files changed, 93 insertions(+), 25 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index bcc41c1affb54..a7217079911aa 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -31,6 +31,7 @@ #include #include +#include #include #include #include @@ -67,11 +68,25 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, PARTHENON_REQUIRE_THROWS(ndim == 1 || ndim == 2, "Histogram dim must be '1' or '2'"); x_var_name = pin->GetString(block_name, prefix + "x_variable"); - - x_var_component = pin->GetInteger(block_name, prefix + "x_variable_component"); - // would add additional logic to pick it from a pack... - PARTHENON_REQUIRE_THROWS(x_var_component >= 0, - "Negative component indices are not supported"); + x_var_component = -1; + if (x_var_name == "COORD_X1") { + x_var_type = VarType::X1; + } else if (x_var_name == "COORD_X2") { + x_var_type = VarType::X2; + } else if (x_var_name == "COORD_X3") { + x_var_type = VarType::X3; + } else if (x_var_name == "COORD_R") { + PARTHENON_REQUIRE_THROWS( + typeid(Coordinates_t) == typeid(UniformCartesian), + "Radial coordinate currently only works for uniform Cartesian coordinates."); + x_var_type = VarType::R; + } else { + x_var_type = VarType::Var; + x_var_component = pin->GetInteger(block_name, prefix + "x_variable_component"); + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(x_var_component >= 0, + "Negative component indices are not supported"); + } const auto x_edges_in = pin->GetVector(block_name, prefix + "x_edges"); // required by binning index function @@ -89,14 +104,28 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, // 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) { y_var_name = pin->GetString(block_name, prefix + "y_variable"); - - y_var_component = pin->GetInteger(block_name, prefix + "y_variable_component"); - // would add additional logic to pick it from a pack... - PARTHENON_REQUIRE_THROWS(y_var_component >= 0, - "Negative component indices are not supported"); + if (y_var_name == "COORD_X1") { + y_var_type = VarType::X1; + } else if (y_var_name == "COORD_X2") { + y_var_type = VarType::X2; + } else if (y_var_name == "COORD_X3") { + y_var_type = VarType::X3; + } else if (y_var_name == "COORD_R") { + PARTHENON_REQUIRE_THROWS( + typeid(Coordinates_t) == typeid(UniformCartesian), + "Radial coordinate currently only works for uniform Cartesian coordinates."); + y_var_type = VarType::R; + } else { + y_var_type = VarType::Var; + y_var_component = pin->GetInteger(block_name, prefix + "y_variable_component"); + // would add additional logic to pick it from a pack... + PARTHENON_REQUIRE_THROWS(y_var_component >= 0, + "Negative component indices are not supported"); + } const auto y_edges_in = pin->GetVector(block_name, prefix + "y_edges"); // required by binning index function @@ -136,6 +165,8 @@ void CalcHist(Mesh *pm, const Histogram &hist) { const auto x_var_component = hist.x_var_component; const auto y_var_component = hist.y_var_component; const auto binned_var_component = hist.binned_var_component; + const auto x_var_type = hist.x_var_type; + const auto y_var_type = hist.y_var_type; const auto x_edges = hist.x_edges; const auto y_edges = hist.y_edges; const auto hist_ndim = hist.ndim; @@ -154,8 +185,14 @@ void CalcHist(Mesh *pm, const Histogram &hist) { for (int p = 0; p < num_partitions; p++) { auto &md = pm->mesh_data.GetOrAdd("base", p); - const auto x_var = md->PackVariables(std::vector{hist.x_var_name}); - const auto y_var = md->PackVariables(std::vector{hist.y_var_name}); + const auto x_var_pack_string = x_var_type == VarType::Var + ? std::vector{hist.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{hist.y_var_name} + : std::vector{}; + const auto y_var = md->PackVariables(y_var_pack_string); const auto binned_var = md->PackVariables(std::vector{hist.binned_var_name}); const auto ib = md->GetBoundsI(IndexDomain::interior); @@ -163,10 +200,23 @@ void CalcHist(Mesh *pm, const Histogram &hist) { const auto kb = md->GetBoundsK(IndexDomain::interior); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "CalcHist", DevExecSpace(), 0, x_var.GetDim(5) - 1, kb.s, + 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) { - const auto &x_val = x_var(b, x_var_component, k, j, 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); + } if (x_val < x_edges(0) || x_val > x_edges(x_edges.extent_int(0) - 1)) { return; } @@ -178,14 +228,28 @@ void CalcHist(Mesh *pm, const Histogram &hist) { int y_bin = 0; if (hist_ndim == 2) { - const auto &y_val = y_var(b, y_var_component, k, j, i); + 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); + } + if (y_val < y_edges(0) || y_val > y_edges(y_edges.extent_int(0) - 1)) { return; } // if we're on the rightmost edge, directly set last bin, otherwise search - const auto y_bin = y_val == y_edges(y_edges.extent_int(0) - 1) - ? y_edges.extent_int(0) - 2 - : upper_bound(y_edges, y_val) - 1; + y_bin = y_val == y_edges(y_edges.extent_int(0) - 1) + ? y_edges.extent_int(0) - 2 + : upper_bound(y_edges, y_val) - 1; } auto res = scatter.access(); res(y_bin, x_bin) += binned_var(b, binned_var_component, k, j, i); @@ -256,7 +320,6 @@ std::string HistogramOutput::GenerateFilename_(ParameterInput *pin, SimTime *tm, // \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_) { CalcHist(pm, hist); diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 089eff04c9760..2cdac27732fd4 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -222,9 +222,12 @@ class PHDF5Output : public OutputType { // \brief derived OutputType class for histograms namespace HistUtil { + +enum class VarType { X1, X2, X3, R, Var, Unused }; struct 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; std::string binned_var_name; // variable name of variable to be binned diff --git a/src/utils/sort.hpp b/src/utils/sort.hpp index aed0deeff473c..9662cdbc58356 100644 --- a/src/utils/sort.hpp +++ b/src/utils/sort.hpp @@ -34,7 +34,7 @@ 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). -template +template KOKKOS_INLINE_FUNCTION int upper_bound(const T &arr, Real val) { int l = 0; int r = arr.extent_int(0); diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index 74f90e61fe3b3..da14f6c3e1831 100644 --- a/tst/regression/test_suites/output_hdf5/output_hdf5.py +++ b/tst/regression/test_suites/output_hdf5/output_hdf5.py @@ -170,6 +170,7 @@ def Analyse(self, parameters): # 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( @@ -184,11 +185,13 @@ def Analyse(self, parameters): print(f"1D 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( - advected.flatten(), + x.flatten(), omadvected.flatten(), - [[1e-9, 1e-4, 1e-1, 2e-1, 5e-1, 1e0], [0, 0.5, 1]], + [[-0.5, -0.25, 0, 0.25, 0.5], [0, 0.5, 1]], weights=advected.flatten(), ) with h5py.File( diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index 812c4ef756050..d495a94126c3b 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -84,9 +84,8 @@ hist0_binned_variable = advected hist0_binned_variable_component = 0 hist1_ndim = 2 -hist1_x_variable = advected -hist1_x_variable_component = 0 -hist1_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.00 +hist1_x_variable = COORD_X1 +hist1_x_edges = -0.5, -0.25, 0.0, 0.25, 0.5 hist1_y_variable = one_minus_advected_sq hist1_y_variable_component = 0 hist1_y_edges = 0, 0.5, 1.0 From b3d144df0512313796c5e4d2e377f88dff4f4bdd Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 16 Oct 2023 14:46:12 +0200 Subject: [PATCH 13/22] Add sampling based hist --- src/outputs/histogram.cpp | 46 ++++++++++++------- src/outputs/outputs.hpp | 5 +- .../test_suites/output_hdf5/output_hdf5.py | 13 +++++- .../output_hdf5/parthinput.advection | 10 +++- 4 files changed, 52 insertions(+), 22 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index a7217079911aa..42c12ac281129 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -69,13 +69,13 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, x_var_name = pin->GetString(block_name, prefix + "x_variable"); x_var_component = -1; - if (x_var_name == "COORD_X1") { + if (x_var_name == "HIST_COORD_X1") { x_var_type = VarType::X1; - } else if (x_var_name == "COORD_X2") { + } else if (x_var_name == "HIST_COORD_X2") { x_var_type = VarType::X2; - } else if (x_var_name == "COORD_X3") { + } else if (x_var_name == "HIST_COORD_X3") { x_var_type = VarType::X3; - } else if (x_var_name == "COORD_R") { + } else if (x_var_name == "HIST_COORD_R") { PARTHENON_REQUIRE_THROWS( typeid(Coordinates_t) == typeid(UniformCartesian), "Radial coordinate currently only works for uniform Cartesian coordinates."); @@ -108,13 +108,13 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, // and for 2D profile check if they're explicitly set (not default value) if (ndim == 2) { y_var_name = pin->GetString(block_name, prefix + "y_variable"); - if (y_var_name == "COORD_X1") { + if (y_var_name == "HIST_COORD_X1") { y_var_type = VarType::X1; - } else if (y_var_name == "COORD_X2") { + } else if (y_var_name == "HIST_COORD_X2") { y_var_type = VarType::X2; - } else if (y_var_name == "COORD_X3") { + } else if (y_var_name == "HIST_COORD_X3") { y_var_type = VarType::X3; - } else if (y_var_name == "COORD_R") { + } else if (y_var_name == "HIST_COORD_R") { PARTHENON_REQUIRE_THROWS( typeid(Coordinates_t) == typeid(UniformCartesian), "Radial coordinate currently only works for uniform Cartesian coordinates."); @@ -143,12 +143,16 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, y_edges = ParArray1D(prefix + "y_edges_unused", 0); } - binned_var_name = pin->GetString(block_name, prefix + "binned_variable"); - 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"); + 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; @@ -189,12 +193,17 @@ void CalcHist(Mesh *pm, const Histogram &hist) { ? std::vector{hist.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{hist.y_var_name} : std::vector{}; const auto y_var = md->PackVariables(y_var_pack_string); - const auto binned_var = - md->PackVariables(std::vector{hist.binned_var_name}); + + const auto binned_var_pack_string = + binned_var_component == -1 ? std::vector{} + : std::vector{hist.binned_var_name}; + const auto binned_var = md->PackVariables(binned_var_pack_string); + const auto ib = md->GetBoundsI(IndexDomain::interior); const auto jb = md->GetBoundsJ(IndexDomain::interior); const auto kb = md->GetBoundsK(IndexDomain::interior); @@ -252,7 +261,10 @@ void CalcHist(Mesh *pm, const Histogram &hist) { : upper_bound(y_edges, y_val) - 1; } auto res = scatter.access(); - res(y_bin, x_bin) += binned_var(b, binned_var_component, k, j, i); + const auto to_add = binned_var_component == -1 + ? 1 + : binned_var(b, binned_var_component, k, j, i); + res(y_bin, x_bin) += to_add; }); // "reduce" results from scatter view to original view. May be a no-op depending on // backend. diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 2cdac27732fd4..64f0a96ff0452 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -231,8 +231,9 @@ struct Histogram { int x_var_component, y_var_component; // components of bin variables (vector) ParArray1D x_edges, y_edges; std::string binned_var_name; // variable name of variable to be binned - int binned_var_component; // component of variable to be binned - ParArray2D result; // resulting histogram + int binned_var_component; // component of variable to be binned. If -1 means no variable + // is binned but the histgram is a sample count. + ParArray2D result; // resulting histogram // temp view for histogram reduction for better performance (switches // between atomics and data duplication depending on the platform) diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index da14f6c3e1831..96a5673778337 100644 --- a/tst/regression/test_suites/output_hdf5/output_hdf5.py +++ b/tst/regression/test_suites/output_hdf5/output_hdf5.py @@ -182,7 +182,7 @@ def Analyse(self, parameters): hist_parth = infile["0/data"][:] all_close = np.allclose(hist_parth, hist_np1d[0]) if not all_close: - print(f"1D hist for {dim}D setup don't match") + 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 @@ -205,4 +205,15 @@ def Analyse(self, parameters): 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, [1e-9, 1e-4, 1e-1, 2e-1, 5e-1, 1e0]) + with h5py.File( + f"advection_{dim}d.out2.histograms.final.hdf", "r" + ) as infile: + hist_parth = infile["2/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 + return analyze_status diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index d495a94126c3b..e31568c762a5b 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -74,7 +74,7 @@ dt = 0.25 file_type = histogram dt = 0.25 -num_histograms = 2 +num_histograms = 3 hist0_ndim = 1 hist0_x_variable = advected @@ -84,7 +84,7 @@ hist0_binned_variable = advected hist0_binned_variable_component = 0 hist1_ndim = 2 -hist1_x_variable = COORD_X1 +hist1_x_variable = HIST_COORD_X1 hist1_x_edges = -0.5, -0.25, 0.0, 0.25, 0.5 hist1_y_variable = one_minus_advected_sq hist1_y_variable_component = 0 @@ -92,3 +92,9 @@ hist1_y_edges = 0, 0.5, 1.0 hist1_binned_variable = advected hist1_binned_variable_component = 0 +hist2_ndim = 1 +hist2_x_variable = advected +hist2_x_variable_component = 0 +hist2_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 +hist2_binned_variable = HIST_ONES + From 4a4d6ae8e66bdc69b8f20214119d66f727793557 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Mon, 16 Oct 2023 15:08:48 +0200 Subject: [PATCH 14/22] Add volume weighting --- src/outputs/histogram.cpp | 12 ++++++---- src/outputs/outputs.hpp | 8 ++++--- .../test_suites/output_hdf5/output_hdf5.py | 23 +++++++++++++++++++ .../output_hdf5/parthinput.advection | 14 ++++++++++- 4 files changed, 49 insertions(+), 8 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 42c12ac281129..9b164a28cbf26 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -159,6 +159,8 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, result = ParArray2D(prefix + "result", nybins, nxbins); scatter_result = Kokkos::Experimental::ScatterView(result.KokkosView()); + + weight_by_vol = pin->GetOrAddBoolean(block_name, prefix + "weight_by_volume", false); } // Computes a 1D or 2D histogram with inclusive lower edges and inclusive rightmost edges. @@ -174,6 +176,7 @@ void CalcHist(Mesh *pm, const Histogram &hist) { const auto x_edges = hist.x_edges; const auto y_edges = hist.y_edges; const auto hist_ndim = hist.ndim; + const auto weight_by_vol = hist.weight_by_vol; auto result = hist.result; auto scatter = hist.scatter_result; @@ -261,10 +264,11 @@ void CalcHist(Mesh *pm, const Histogram &hist) { : upper_bound(y_edges, y_val) - 1; } auto res = scatter.access(); - const auto to_add = binned_var_component == -1 - ? 1 - : binned_var(b, binned_var_component, k, j, i); - res(y_bin, x_bin) += to_add; + const auto val_to_add = binned_var_component == -1 + ? 1 + : binned_var(b, binned_var_component, k, j, i); + const auto weight = weight_by_vol ? coords.CellVolume(k, j, i) : 1.0; + 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. diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 64f0a96ff0452..acba841dd3ff1 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -231,9 +231,11 @@ struct Histogram { int x_var_component, y_var_component; // components of bin variables (vector) ParArray1D x_edges, y_edges; std::string binned_var_name; // variable name of variable to be binned - int binned_var_component; // component of variable to be binned. If -1 means no variable - // is binned but the histgram is a sample count. - ParArray2D result; // resulting histogram + // 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 + ParArray2D result; // resulting histogram // temp view for histogram reduction for better performance (switches // between atomics and data duplication depending on the platform) diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index 96a5673778337..857f528c578fc 100644 --- a/tst/regression/test_suites/output_hdf5/output_hdf5.py +++ b/tst/regression/test_suites/output_hdf5/output_hdf5.py @@ -216,4 +216,27 @@ def Analyse(self, parameters): 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(), + ) + with h5py.File( + f"advection_{dim}d.out2.histograms.final.hdf", "r" + ) as infile: + hist_parth = infile["3/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 e31568c762a5b..fd12704472809 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -74,8 +74,9 @@ dt = 0.25 file_type = histogram dt = 0.25 -num_histograms = 3 +num_histograms = 4 +# 1D histogram of a variable, binned by a variable hist0_ndim = 1 hist0_x_variable = advected hist0_x_variable_component = 0 @@ -83,6 +84,7 @@ hist0_x_edges = 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 hist1_ndim = 2 hist1_x_variable = HIST_COORD_X1 hist1_x_edges = -0.5, -0.25, 0.0, 0.25, 0.5 @@ -92,9 +94,19 @@ hist1_y_edges = 0, 0.5, 1.0 hist1_binned_variable = advected hist1_binned_variable_component = 0 +# 1D histogram ("standard", i.e., counting occurance in bin) hist2_ndim = 1 hist2_x_variable = advected hist2_x_variable_component = 0 hist2_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 hist2_binned_variable = HIST_ONES +# 2D histogram of volume weighted variable according to two coordinates +hist3_ndim = 2 +hist3_x_variable = HIST_COORD_X1 +hist3_x_edges = -0.5, -0.25, 0.0, 0.25, 0.5 +hist3_y_variable = HIST_COORD_X2 +hist3_y_edges = -0.5, -0.1, 0.0, 0.1, 0.5 +hist3_binned_variable = advected +hist3_binned_variable_component = 0 +hist3_weight_by_volume = true From 3eb6df90067dcb5445d306a78ce9afc6a8e83df1 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 17 Oct 2023 14:23:22 +0200 Subject: [PATCH 15/22] Fix Scatterview layout --- src/outputs/histogram.cpp | 3 ++- src/outputs/outputs.hpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 9b164a28cbf26..1d0855a74960b 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -158,7 +158,8 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, 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()); + scatter_result = + Kokkos::Experimental::ScatterView(result.KokkosView()); weight_by_vol = pin->GetOrAddBoolean(block_name, prefix + "weight_by_volume", false); } diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index acba841dd3ff1..6fe08bb3af0fc 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -31,6 +31,7 @@ #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" @@ -239,7 +240,7 @@ struct Histogram { // temp view for histogram reduction for better performance (switches // between atomics and data duplication depending on the platform) - Kokkos::Experimental::ScatterView scatter_result; + Kokkos::Experimental::ScatterView scatter_result; Histogram(ParameterInput *pin, const std::string &block_name, const std::string &prefix); }; From 5961d213127ef3dbd075b089a00953bfd493fe9a Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 17 Oct 2023 15:00:38 +0200 Subject: [PATCH 16/22] Separate edge parsing --- src/outputs/histogram.cpp | 58 +++++++++++-------- .../output_hdf5/parthinput.advection | 18 ++++-- 2 files changed, 46 insertions(+), 30 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 1d0855a74960b..94158b55635e0 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -62,6 +62,37 @@ using namespace OutputUtils; namespace HistUtil { +ParArray1D GetEdges(ParameterInput *pin, const std::string &block_name, + const std::string &prefix) { + + std::vector edges_in; + + const auto edge_type_str = pin->GetString(block_name, prefix + "type"); + if (edge_type_str == "lin") { + + } else if (edge_type_str == "log") { + + } else if (edge_type_str == "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, + "Need at least one bin, i.e., two 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); + return edges; +} + Histogram::Histogram(ParameterInput *pin, const std::string &block_name, const std::string &prefix) { ndim = pin->GetInteger(block_name, prefix + "ndim"); @@ -88,18 +119,7 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, "Negative component indices are not supported"); } - const auto x_edges_in = pin->GetVector(block_name, prefix + "x_edges"); - // required by binning index function - PARTHENON_REQUIRE_THROWS(std::is_sorted(x_edges_in.begin(), x_edges_in.end()), - "Bin edges must be in order."); - PARTHENON_REQUIRE_THROWS(x_edges_in.size() >= 2, - "Need at least one bin, i.e., two edges."); - x_edges = ParArray1D(prefix + "x_edges", x_edges_in.size()); - auto x_edges_h = x_edges.GetHostMirror(); - for (int i = 0; i < x_edges_in.size(); i++) { - x_edges_h(i) = x_edges_in[i]; - } - Kokkos::deep_copy(x_edges, x_edges_h); + x_edges = GetEdges(pin, block_name, prefix + "x_edges_"); // For 1D profile default initalize y variables y_var_name = ""; @@ -127,18 +147,8 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, "Negative component indices are not supported"); } - const auto y_edges_in = pin->GetVector(block_name, prefix + "y_edges"); - // required by binning index function - PARTHENON_REQUIRE_THROWS(std::is_sorted(y_edges_in.begin(), y_edges_in.end()), - "Bin edges must be in order."); - PARTHENON_REQUIRE_THROWS(y_edges_in.size() >= 2, - "Need at least one bin, i.e., two edges."); - y_edges = ParArray1D(prefix + "y_edges", y_edges_in.size()); - auto y_edges_h = y_edges.GetHostMirror(); - for (int i = 0; i < y_edges_in.size(); i++) { - y_edges_h(i) = y_edges_in[i]; - } - Kokkos::deep_copy(y_edges, y_edges_h); + y_edges = GetEdges(pin, block_name, prefix + "y_edges_"); + } else { y_edges = ParArray1D(prefix + "y_edges_unused", 0); } diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index fd12704472809..677a2ebf2e18c 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -80,17 +80,20 @@ num_histograms = 4 hist0_ndim = 1 hist0_x_variable = advected hist0_x_variable_component = 0 -hist0_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.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 hist1_ndim = 2 hist1_x_variable = HIST_COORD_X1 -hist1_x_edges = -0.5, -0.25, 0.0, 0.25, 0.5 +hist1_x_edges_type = list +hist1_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5 hist1_y_variable = one_minus_advected_sq hist1_y_variable_component = 0 -hist1_y_edges = 0, 0.5, 1.0 +hist1_y_edges_type = list +hist1_y_edges_list = 0, 0.5, 1.0 hist1_binned_variable = advected hist1_binned_variable_component = 0 @@ -98,15 +101,18 @@ hist1_binned_variable_component = 0 hist2_ndim = 1 hist2_x_variable = advected hist2_x_variable_component = 0 -hist2_x_edges = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 +hist2_x_edges_type = list +hist2_x_edges_list = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 hist2_binned_variable = HIST_ONES # 2D histogram of volume weighted variable according to two coordinates hist3_ndim = 2 hist3_x_variable = HIST_COORD_X1 -hist3_x_edges = -0.5, -0.25, 0.0, 0.25, 0.5 +hist3_x_edges_type = list +hist3_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5 hist3_y_variable = HIST_COORD_X2 -hist3_y_edges = -0.5, -0.1, 0.0, 0.1, 0.5 +hist3_y_edges_type = list +hist3_y_edges_list = -0.5, -0.1, 0.0, 0.1, 0.5 hist3_binned_variable = advected hist3_binned_variable_component = 0 hist3_weight_by_volume = true From e515f7c3dc848a88c6053aa2d97e0c2bf9e87b21 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Tue, 17 Oct 2023 15:24:26 +0200 Subject: [PATCH 17/22] Add support for calculating lin and log bin edges --- src/outputs/histogram.cpp | 31 +++++++++++++++++-- .../output_hdf5/parthinput.advection | 6 ++-- 2 files changed, 32 insertions(+), 5 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 94158b55635e0..046dcaec3fbc4 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -64,13 +64,38 @@ namespace HistUtil { ParArray1D GetEdges(ParameterInput *pin, const std::string &block_name, const std::string &prefix) { - std::vector edges_in; const auto edge_type_str = pin->GetString(block_name, prefix + "type"); - if (edge_type_str == "lin") { + if (edge_type_str == "lin" || edge_type_str == "log") { + const auto 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->GetReal(block_name, prefix + "num_bins"); + PARTHENON_REQUIRE_THROWS(edge_num_bins >= 1, "Need at least one bin for histogram."); + + if (edge_type_str == "lin") { + auto dbin = (edge_max - edge_min) / (edge_num_bins); + for (int i = 0; i < edge_num_bins; i++) { + edges_in.emplace_back(edge_min + i * dbin); + } + edges_in.emplace_back(edge_max); + } else if (edge_type_str == "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.") - } else if (edge_type_str == "log") { + auto dbin = (std::log10(edge_max) - std::log10(edge_min)) / (edge_num_bins); + for (int i = 0; i < edge_num_bins; i++) { + edges_in.emplace_back(std::pow(10., std::log10(edge_min) + i * dbin)); + } + edges_in.emplace_back(edge_max); + } else { + PARTHENON_FAIL("Not sure how I got here...") + } } else if (edge_type_str == "list") { edges_in = pin->GetVector(block_name, prefix + "list"); diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index 677a2ebf2e18c..347ff6d75644d 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -88,8 +88,10 @@ hist0_binned_variable_component = 0 # 2D histogram of a variable, binned by a coordinate and a different variable hist1_ndim = 2 hist1_x_variable = HIST_COORD_X1 -hist1_x_edges_type = list -hist1_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5 +hist1_x_edges_type = lin +hist1_x_edges_num_bins = 4 +hist1_x_edges_min = -0.5 +hist1_x_edges_max = 0.5 hist1_y_variable = one_minus_advected_sq hist1_y_variable_component = 0 hist1_y_edges_type = list From 102baa058648666255ae62cf2c0f0d31c1b900ed Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 19 Oct 2023 09:43:16 +0200 Subject: [PATCH 18/22] Write string as strings to HDF5 --- src/outputs/histogram.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 046dcaec3fbc4..6215338ffe475 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -418,9 +418,9 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm auto &hist = histograms_[h]; const H5G hist_group = MakeGroup(file, "/" + std::to_string(h)); HDF5WriteAttribute("ndim", hist.ndim, hist_group); - HDF5WriteAttribute("x_var_name", hist.x_var_name, 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, 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(); @@ -429,7 +429,7 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm local_count.data(), global_count.data(), pl_xfer); if (hist.ndim == 2) { - HDF5WriteAttribute("y_var_name", hist.y_var_name, hist_group); + 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(); From 6b12bc4ceeddc909b4c5ff6bcc1b17d47bade743 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 19 Oct 2023 09:58:41 +0200 Subject: [PATCH 19/22] Add support for variable weights to histogram --- src/outputs/histogram.cpp | 22 ++++++++++++- src/outputs/outputs.hpp | 5 ++- .../test_suites/output_hdf5/output_hdf5.py | 6 ++-- .../output_hdf5/parthinput.advection | 32 ++++++++++++------- 4 files changed, 49 insertions(+), 16 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 6215338ffe475..8b2a04bd20527 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -197,6 +197,17 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, Kokkos::Experimental::ScatterView(result.KokkosView()); 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. @@ -207,6 +218,7 @@ void CalcHist(Mesh *pm, const Histogram &hist) { const auto x_var_component = hist.x_var_component; const auto y_var_component = hist.y_var_component; const auto binned_var_component = hist.binned_var_component; + const auto weight_var_component = hist.weight_var_component; const auto x_var_type = hist.x_var_type; const auto y_var_type = hist.y_var_type; const auto x_edges = hist.x_edges; @@ -243,6 +255,11 @@ void CalcHist(Mesh *pm, const Histogram &hist) { : std::vector{hist.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{hist.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); @@ -303,7 +320,10 @@ void CalcHist(Mesh *pm, const Histogram &hist) { const auto val_to_add = binned_var_component == -1 ? 1 : binned_var(b, binned_var_component, k, j, i); - const auto weight = weight_by_vol ? coords.CellVolume(k, j, i) : 1.0; + 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 diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 6fe08bb3af0fc..00793e0f2dfe1 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -235,7 +235,10 @@ struct Histogram { // 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 + 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 diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index 857f528c578fc..64e04b2c65847 100644 --- a/tst/regression/test_suites/output_hdf5/output_hdf5.py +++ b/tst/regression/test_suites/output_hdf5/output_hdf5.py @@ -224,12 +224,12 @@ def Analyse(self, parameters): 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(), + weights=advected.flatten() * vols.flatten() * omadvected.flatten(), ) with h5py.File( - f"advection_{dim}d.out2.histograms.final.hdf", "r" + f"advection_{dim}d.out3.histograms.final.hdf", "r" ) as infile: - hist_parth = infile["3/data"][:] + hist_parth = infile["0/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]) diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index 347ff6d75644d..e5981dfd77dbf 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -74,7 +74,7 @@ dt = 0.25 file_type = histogram dt = 0.25 -num_histograms = 4 +num_histograms = 3 # 1D histogram of a variable, binned by a variable hist0_ndim = 1 @@ -107,14 +107,24 @@ hist2_x_edges_type = list hist2_x_edges_list = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 hist2_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 + +num_histograms = 1 + # 2D histogram of volume weighted variable according to two coordinates -hist3_ndim = 2 -hist3_x_variable = HIST_COORD_X1 -hist3_x_edges_type = list -hist3_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5 -hist3_y_variable = HIST_COORD_X2 -hist3_y_edges_type = list -hist3_y_edges_list = -0.5, -0.1, 0.0, 0.1, 0.5 -hist3_binned_variable = advected -hist3_binned_variable_component = 0 -hist3_weight_by_volume = true +hist0_ndim = 2 +hist0_x_variable = HIST_COORD_X1 +hist0_x_edges_type = list +hist0_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5 +hist0_y_variable = HIST_COORD_X2 +hist0_y_edges_type = list +hist0_y_edges_list = -0.5, -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 From ead3aecc63315c721ad8d3445d366b4533dd8dc2 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 19 Oct 2023 10:53:24 +0200 Subject: [PATCH 20/22] Use direct indexing for lin and log bins --- src/outputs/histogram.cpp | 85 ++++++++++++++++++++++++++++++--------- src/outputs/outputs.hpp | 7 ++++ 2 files changed, 73 insertions(+), 19 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index 8b2a04bd20527..e976884f25491 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -62,13 +62,19 @@ using namespace OutputUtils; namespace HistUtil { -ParArray1D GetEdges(ParameterInput *pin, const std::string &block_name, - const std::string &prefix) { +// 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. +std::tuple, EdgeType, Real, Real> +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") { - const auto edge_min = pin->GetReal(block_name, prefix + "min"); + 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.") @@ -77,20 +83,24 @@ ParArray1D GetEdges(ParameterInput *pin, const std::string &block_name, PARTHENON_REQUIRE_THROWS(edge_num_bins >= 1, "Need at least one bin for histogram."); if (edge_type_str == "lin") { - auto dbin = (edge_max - edge_min) / (edge_num_bins); + 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 * dbin); + 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.") - auto dbin = (std::log10(edge_max) - std::log10(edge_min)) / (edge_num_bins); + // 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., std::log10(edge_min) + i * dbin)); + edges_in.emplace_back(std::pow(10., edge_min + i * edge_dbin)); } edges_in.emplace_back(edge_max); } else { @@ -98,6 +108,7 @@ ParArray1D GetEdges(ParameterInput *pin, const std::string &block_name, } } 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()), @@ -115,7 +126,11 @@ ParArray1D GetEdges(ParameterInput *pin, const std::string &block_name, edges_h(i) = edges_in[i]; } Kokkos::deep_copy(edges, edges_h); - return edges; + + PARTHENON_REQUIRE_THROWS( + edge_type != EdgeType::Undefined, + "Edge type not set and it's unclear how this code was triggered..."); + return {edges, edge_type, edge_min, edge_dbin}; } Histogram::Histogram(ParameterInput *pin, const std::string &block_name, @@ -144,7 +159,8 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, "Negative component indices are not supported"); } - x_edges = GetEdges(pin, block_name, prefix + "x_edges_"); + 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 = ""; @@ -172,7 +188,8 @@ Histogram::Histogram(ParameterInput *pin, const std::string &block_name, "Negative component indices are not supported"); } - y_edges = GetEdges(pin, block_name, prefix + "y_edges_"); + 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); @@ -223,6 +240,12 @@ void CalcHist(Mesh *pm, const Histogram &hist) { const auto y_var_type = hist.y_var_type; const auto x_edges = hist.x_edges; const auto y_edges = hist.y_edges; + const auto x_edges_type = hist.x_edges_type; + const auto y_edges_type = hist.y_edges_type; + const auto x_edge_min = hist.x_edge_min; + const auto x_edge_dbin = hist.x_edge_dbin; + const auto y_edge_min = hist.y_edge_min; + const auto y_edge_dbin = hist.y_edge_dbin; const auto hist_ndim = hist.ndim; const auto weight_by_vol = hist.weight_by_vol; auto result = hist.result; @@ -286,12 +309,24 @@ void CalcHist(Mesh *pm, const Histogram &hist) { return; } - // if we're on the rightmost edge, directly set last bin, otherwise search - const auto x_bin = x_val == x_edges(x_edges.extent_int(0) - 1) - ? x_edges.extent_int(0) - 2 - : upper_bound(x_edges, x_val) - 1; + int x_bin = -1; + // if we're on the rightmost edge, directly set last bin + 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"); - int y_bin = 0; + int y_bin = -1; if (hist_ndim == 2) { auto y_val = std::numeric_limits::quiet_NaN(); if (y_var_type == VarType::X1) { @@ -311,10 +346,22 @@ void CalcHist(Mesh *pm, const Histogram &hist) { if (y_val < y_edges(0) || y_val > y_edges(y_edges.extent_int(0) - 1)) { return; } - // if we're on the rightmost edge, directly set last bin, otherwise search - y_bin = y_val == y_edges(y_edges.extent_int(0) - 1) - ? y_edges.extent_int(0) - 2 - : upper_bound(y_edges, y_val) - 1; + // if we're on the rightmost edge, directly set last bin + 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 diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 00793e0f2dfe1..eb975ef27e513 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -225,12 +225,19 @@ class PHDF5Output : public OutputType { namespace HistUtil { enum class VarType { X1, X2, X3, R, Var, Unused }; +enum class EdgeType { Lin, Log, List, Undefined }; + struct 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; 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. From 2ac1c9aa98c379e087b12f5b8a5132cdc2ed15ea Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 19 Oct 2023 11:05:49 +0200 Subject: [PATCH 21/22] Fix test case for direct log indexed bins --- src/outputs/histogram.cpp | 6 +++++- tst/regression/test_suites/output_hdf5/output_hdf5.py | 2 +- tst/regression/test_suites/output_hdf5/parthinput.advection | 6 ++++-- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index e976884f25491..34779dbf65151 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -326,7 +326,9 @@ void CalcHist(Mesh *pm, const Histogram &hist) { } PARTHENON_DEBUG_REQUIRE(x_bin >= 0, "Bin not found"); - int y_bin = -1; + // 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) { @@ -346,6 +348,8 @@ void CalcHist(Mesh *pm, const Histogram &hist) { if (y_val < y_edges(0) || y_val > y_edges(y_edges.extent_int(0) - 1)) { return; } + + y_bin = -1; // reset to impossible value // if we're on the rightmost edge, directly set last bin if (y_val == y_edges(y_edges.extent_int(0) - 1)) { y_bin = y_edges.extent_int(0) - 2; diff --git a/tst/regression/test_suites/output_hdf5/output_hdf5.py b/tst/regression/test_suites/output_hdf5/output_hdf5.py index 64e04b2c65847..1de94b678448f 100644 --- a/tst/regression/test_suites/output_hdf5/output_hdf5.py +++ b/tst/regression/test_suites/output_hdf5/output_hdf5.py @@ -206,7 +206,7 @@ def Analyse(self, parameters): analyze_status = False # 1D histogram (simple sampling) with bins defined by a var - hist_np1d = np.histogram(advected, [1e-9, 1e-4, 1e-1, 2e-1, 5e-1, 1e0]) + 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: diff --git a/tst/regression/test_suites/output_hdf5/parthinput.advection b/tst/regression/test_suites/output_hdf5/parthinput.advection index e5981dfd77dbf..8442d082e8ec1 100644 --- a/tst/regression/test_suites/output_hdf5/parthinput.advection +++ b/tst/regression/test_suites/output_hdf5/parthinput.advection @@ -103,8 +103,10 @@ hist1_binned_variable_component = 0 hist2_ndim = 1 hist2_x_variable = advected hist2_x_variable_component = 0 -hist2_x_edges_type = list -hist2_x_edges_list = 1e-9, 1e-4,1e-1, 2e-1, 5e-1 ,1.0 +hist2_x_edges_type = log +hist2_x_edges_num_bins = 10 +hist2_x_edges_min = 1e-9 +hist2_x_edges_max = 1e0 hist2_binned_variable = HIST_ONES # A second output block with different dt for histograms From bc2684f33ac09732e050a5b64f883353cda54f7c Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Thu, 19 Oct 2023 17:00:14 +0200 Subject: [PATCH 22/22] Add doc --- CHANGELOG.md | 1 + doc/sphinx/src/interface/state.rst | 2 +- doc/sphinx/src/outputs.rst | 135 ++++++++++++++++++++++++++++- 3 files changed, 135 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 49d929c46ab5c..0a94116ccd52e 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 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 diff --git a/doc/sphinx/src/interface/state.rst b/doc/sphinx/src/interface/state.rst index b2a91fa8e4cc2..ae183c1cb2949 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/outputs.rst b/doc/sphinx/src/outputs.rst index 574d52f4a212e..2e1d24fdb3d6b 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,134 @@ 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 +- 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 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 + num_histograms = 2 # required, specifies how many histograms are defined in this block + + # 1D histogram ("standard", i.e., counting occurance in bin) + hist0_ndim = 1 + hist0_x_variable = advected + hist0_x_variable_component = 0 + hist0_x_edges_type = log + hist0_x_edges_num_bins = 10 + hist0_x_edges_min = 1e-9 + hist0_x_edges_max = 1e0 + hist0_binned_variable = HIST_ONES + + # 2D histogram of volume weighted variable according to two coordinates + hist1_ndim = 2 + hist1_x_variable = HIST_COORD_X1 + hist1_x_edges_type = list + hist1_x_edges_list = -0.5, -0.25, 0.0, 0.25, 0.5 + hist1_y_variable = HIST_COORD_X2 + hist1_y_edges_type = list + hist1_y_edges_list = -0.5, -0.1, 0.0, 0.1, 0.5 + hist1_binned_variable = advected + hist1_binned_variable_component = 0 + hist1_weight_by_volume = true + hist1_weight_variable = one_minus_advected_sq + hist1_weight_variable_component = 0 + +with the following parameters + +- ``num_histograms=INT`` + The number of histograms defined in this block. + All histogram definitions need to be prefix with ``hist#_`` where ``#`` is the + histogram number starting to count from ``0``. + All histograms will be written to the same output file with the "group" in the + output corresponding to the histogram number. +- ``hist#_ndim=INT`` (either ``1`` or ``2``) + Dimensionality of the histogram. +- ``hist#_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. +- ``hist#_x_variable_component=INT`` + Component index of the binning variable. + Used/required only if a non-coordinate variable is used for binning. +- ``hist#_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. +- ``hist#_x_edges_min=FLOAT`` + Minimum value (inclusive) of the bins in the first dim. + Used/required only for ``lin`` and ``log`` edge type. +- ``hist#_x_edges_max=FLOAT`` + Maximum value (inclusive) of the bins in the first dim. + Used/required only for ``lin`` and ``log`` edge type. +- ``hist#_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. +- ``hist#_x_edges_list=FLOAT,FLOAT,FLOAT,...`` (comma separated list of increasing values) + Arbitrary definition of edge values. + Used/required only for ``list`` edge type. +- ``hist#_y_edges...`` + Same as the ``hist#_x_edges...`` parameters except for being used in the second + dimension for ``ndim=2`` histograms. +- ``hist#_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. +- ``hist#_binned_variable_component=INT`` + Component index of the variable to be binned. + Used/required only if a variable is binned and not ``HIST_ONES``. +- ``hist#_weight_by_volume=BOOL`` (``true`` or ``false``) + Apply volume weighting to the binned variable. Can be used simultaneously with binning + by a different variable. +- ``hist#_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. +- ``hist#_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. + +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["0/x_edges"][:] + y = infile["0/data"][:] + plt.plot(x, y) + plt.show() + + # 2D histogram + x = infile["1/x_edges"][:] + y = infile["1/y_edges"][:] + z = infile["1/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) -----------------