From 8e55e3173f8fde17b69773bc8cfe70a111f483a9 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Fri, 11 Jul 2025 11:19:20 -0700 Subject: [PATCH 01/15] added histogram diagnostic that is passing initial tests --- .../eamxx/src/diagnostics/CMakeLists.txt | 1 + .../eamxx/src/diagnostics/histogram.cpp | 94 +++++++++ .../eamxx/src/diagnostics/histogram.hpp | 41 ++++ .../src/diagnostics/register_diagnostics.hpp | 2 + .../src/diagnostics/tests/CMakeLists.txt | 5 +- .../src/diagnostics/tests/histogram_test.cpp | 189 ++++++++++++++++++ 6 files changed, 331 insertions(+), 1 deletion(-) create mode 100644 components/eamxx/src/diagnostics/histogram.cpp create mode 100644 components/eamxx/src/diagnostics/histogram.hpp create mode 100644 components/eamxx/src/diagnostics/tests/histogram_test.cpp diff --git a/components/eamxx/src/diagnostics/CMakeLists.txt b/components/eamxx/src/diagnostics/CMakeLists.txt index 0cf2178ade05..5d7faaa944fd 100644 --- a/components/eamxx/src/diagnostics/CMakeLists.txt +++ b/components/eamxx/src/diagnostics/CMakeLists.txt @@ -9,6 +9,7 @@ set(DIAGNOSTIC_SRCS field_at_height.cpp field_at_level.cpp field_at_pressure_level.cpp + histogram.cpp horiz_avg.cpp longwave_cloud_forcing.cpp number_path.cpp diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp new file mode 100644 index 000000000000..90c163060cd4 --- /dev/null +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -0,0 +1,94 @@ +#include "diagnostics/histogram.hpp" + +namespace scream { + + +HistogramDiag::HistogramDiag(const ekat::Comm &comm, const ekat::ParameterList ¶ms) + : AtmosphereDiagnostic(comm, params) { + const auto &field_name = m_params.get("field_name"); + m_bin_configuration = params.get("bin_configuration"); + m_diag_name = field_name + "_histogram_" + m_bin_configuration; +} + +void HistogramDiag::set_grids(const std::shared_ptr grids_manager) { + const auto &field_name = m_params.get("field_name"); + const auto &grid_name = m_params.get("grid_name"); + add_field(field_name, grid_name); +} + +void HistogramDiag::initialize_impl(const RunType /*run_type*/) { + using FieldIdentifier = FieldHeader::identifier_type; + using FieldLayout = FieldIdentifier::layout_type; + using ShortFieldTagsNames::CMP; + const Field &field = get_fields_in().front(); + const FieldIdentifier &field_id = field.get_header().get_identifier(); + + // extract bin values from configuration + std::istringstream stream(m_bin_configuration); + std::string token; + std::vector bin_values; + while (std::getline(stream, token, '_')) { + bin_values.push_back(std::stod(token)); + } + int num_bins = bin_values.size()-1; + + // allocate histogram field + FieldLayout diagnostic_layout({CMP}, {num_bins}, {"bin"}); + FieldIdentifier diagnostic_id(m_diag_name, diagnostic_layout, field_id.get_units(), + field_id.get_grid_name(), DataType::IntType); + m_diagnostic_output = Field(diagnostic_id); + m_diagnostic_output.allocate_view(); + + // allocate field for bin values + FieldLayout bin_values_layout({CMP}, {num_bins+1}, {"bin"}); + FieldIdentifier bin_values_id(m_diag_name + "_bin_values", bin_values_layout, + field_id.get_units(), field_id.get_grid_name()); + m_bin_values = Field(bin_values_id); + m_bin_values.allocate_view(); + auto bin_values_view = m_bin_values.get_view(); + using RangePolicy = Kokkos::RangePolicy; + Kokkos::parallel_for("store_histogram_bin_values_" + field.name(), + RangePolicy(0, bin_values_layout.dim(0)), [&] (int i) { + bin_values_view(i) = bin_values[i]; + }); +} + +void HistogramDiag::compute_diagnostic_impl() { + const auto &field = get_fields_in().front(); + auto field_layout = field.get_header().get_identifier().get_layout(); + auto histogram_layout = m_diagnostic_output.get_header().get_identifier().get_layout(); + const int num_bins = histogram_layout.dim(0); + const int field_size = field_layout.size(); + + auto bin_values_view = m_bin_values.get_view(); + auto histogram_view = m_diagnostic_output.get_view(); + auto field_view = field.get_view(); + using KT = ekat::KokkosTypes; + using TeamPolicy = Kokkos::TeamPolicy; + using TeamMember = typename TeamPolicy::member_type; + using ESU = ekat::ExeSpaceUtils; + TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, field_size); + Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, + KOKKOS_LAMBDA(const TeamMember &tm) { + const int bin_i = tm.league_rank(); + const Real bin_lower = bin_values_view(bin_i); + const Real bin_upper = bin_values_view(bin_i+1); + Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, field_size), + [&](int i, Int &val) { + if ((bin_lower <= field_view(i)) && (field_view(i) < bin_upper)) + val++; + }, + histogram_view(bin_i)); + }); + + // TODO: use device-side MPI calls + // TODO: the dev ptr causes problems; revisit this later + // TODO: doing cuda-aware MPI allreduce would be ~10% faster + Kokkos::fence(); + m_diagnostic_output.sync_to_host(); + m_comm.all_reduce(m_diagnostic_output.template get_internal_view_data(), + histogram_layout.size(), MPI_SUM); + m_diagnostic_output.sync_to_dev(); +} + +} // namespace scream diff --git a/components/eamxx/src/diagnostics/histogram.hpp b/components/eamxx/src/diagnostics/histogram.hpp new file mode 100644 index 000000000000..2e05543c32bb --- /dev/null +++ b/components/eamxx/src/diagnostics/histogram.hpp @@ -0,0 +1,41 @@ +#ifndef EAMXX_HISTOTGRAM_HPP +#define EAMXX_HISTOGRAM_HPP + +#include "share/atm_process/atmosphere_diagnostic.hpp" + +namespace scream { +/* + * This diagnostic will calculate a histogram of a field across all dimensions + * producing a one dimensional field, with CMP tag dimension named "bin", that + * indicates how many times a field value in the specified range occurred. + */ + +class HistogramDiag : public AtmosphereDiagnostic { + +public: + // Constructors + HistogramDiag(const ekat::Comm &comm, const ekat::ParameterList ¶ms); + + // The name of the diagnostic + std::string name() const { return m_diag_name; } + + // Set the grid + void set_grids(const std::shared_ptr grids_manager); + +protected: +#ifdef KOKKOS_ENABLE_CUDA +public: +#endif + void initialize_impl(const RunType /*run_type*/); + void compute_diagnostic_impl(); + +protected: + std::string m_diag_name; + std::string m_bin_configuration; + + Field m_bin_values; +}; + +} // namespace scream + +#endif // EAMXX_HISTOGRAM_HPP diff --git a/components/eamxx/src/diagnostics/register_diagnostics.hpp b/components/eamxx/src/diagnostics/register_diagnostics.hpp index b7e4a3aa4d6b..3c78a0088672 100644 --- a/components/eamxx/src/diagnostics/register_diagnostics.hpp +++ b/components/eamxx/src/diagnostics/register_diagnostics.hpp @@ -30,6 +30,7 @@ #include "diagnostics/zonal_avg.hpp" #include "diagnostics/conditional_sampling.hpp" #include "diagnostics/binary_ops.hpp" +#include "diagnostics/histogram.hpp" namespace scream { @@ -63,6 +64,7 @@ inline void register_diagnostics () { diag_factory.register_product("ZonalAvgDiag",&create_atmosphere_diagnostic); diag_factory.register_product("ConditionalSampling",&create_atmosphere_diagnostic); diag_factory.register_product("BinaryOpsDiag", &create_atmosphere_diagnostic); + diag_factory.register_product("HistogramDiag",&create_atmosphere_diagnostic); } } // namespace scream diff --git a/components/eamxx/src/diagnostics/tests/CMakeLists.txt b/components/eamxx/src/diagnostics/tests/CMakeLists.txt index 0dd755b1e27b..0516aadba818 100644 --- a/components/eamxx/src/diagnostics/tests/CMakeLists.txt +++ b/components/eamxx/src/diagnostics/tests/CMakeLists.txt @@ -82,10 +82,13 @@ CreateDiagTest(vert_contract "vert_contract_test.cpp" MPI_RANKS 1 ${SCREAM_TEST_ CreateDiagTest(vert_derivative "vert_derivative_test.cpp") # Test zonal averaging -CreateDiagTest(zonal_avg zonal_avg_test.cpp MPI_RANKS 1 ${SCREAM_TEST_MAX_RANKS}) +CreateDiagTest(zonal_avg "zonal_avg_test.cpp" MPI_RANKS 1 ${SCREAM_TEST_MAX_RANKS}) # Test conditional sampling CreateDiagTest(conditional_sampling "conditional_sampling_test.cpp") # Test binary ops CreateDiagTest(binary_ops "binary_ops_test.cpp") + +# Test histogram +CreateDiagTest(histogram "histogram_test.cpp" MPI_RANKS 1 ${SCREAM_TEST_MAX_RANKS}) diff --git a/components/eamxx/src/diagnostics/tests/histogram_test.cpp b/components/eamxx/src/diagnostics/tests/histogram_test.cpp new file mode 100644 index 000000000000..5901fc7b0f5a --- /dev/null +++ b/components/eamxx/src/diagnostics/tests/histogram_test.cpp @@ -0,0 +1,189 @@ +#include "catch2/catch.hpp" +#include "diagnostics/register_diagnostics.hpp" +#include "share/field/field_utils.hpp" +#include "share/grid/mesh_free_grids_manager.hpp" +#include "share/util/eamxx_setup_random_test.hpp" +#include "share/util/eamxx_universal_constants.hpp" + +namespace scream { + +std::shared_ptr create_gm(const ekat::Comm &comm, const int ncols, const int nlevs) { + const int num_global_cols = ncols * comm.size(); + + using vos_t = std::vector; + ekat::ParameterList gm_params; + gm_params.set("grids_names", vos_t{"Point Grid"}); + auto &pl = gm_params.sublist("Point Grid"); + pl.set("type", "point_grid"); + pl.set("aliases", vos_t{"Physics"}); + pl.set("number_of_global_columns", num_global_cols); + pl.set("number_of_vertical_levels", nlevs); + + auto gm = create_mesh_free_grids_manager(comm, gm_params); + gm->build_grids(); + + return gm; +} + +TEST_CASE("histogram") { + using namespace ShortFieldTagsNames; + using namespace ekat::units; + + // A numerical tolerance + auto tol = std::numeric_limits::epsilon() * 100; + + // A world comm + ekat::Comm comm(MPI_COMM_WORLD); + + // A time stamp + util::TimeStamp t0({2024, 1, 1}, {0, 0, 0}); + + // Create a grids manager - single column for these tests + constexpr int nlevs = 3; + constexpr int dim3 = 4; + const int ngcols = 6 * comm.size(); + + auto gm = create_gm(comm, ngcols, nlevs); + auto grid = gm->get_grid("Physics"); + + // Specify histogram bins + const std::string bin_configuration = "0_50_100_150_200"; + const std::vector bin_values = {0.0, 50.0, 100.0, 150.0, 200.0}; + + // Input (randomized) qc + FieldLayout scalar1d_layout{{COL}, {ngcols}}; + FieldLayout scalar2d_layout{{COL, LEV}, {ngcols, nlevs}}; + FieldLayout scalar3d_layout{{COL, CMP, LEV}, {ngcols, dim3, nlevs}}; + + FieldIdentifier qc1_id("qc", scalar1d_layout, kg / kg, grid->name()); + FieldIdentifier qc2_fid("qc", scalar2d_layout, kg / kg, grid->name()); + FieldIdentifier qc3_fid("qc", scalar3d_layout, kg / kg, grid->name()); + + Field qc1(qc1_id); + Field qc2(qc2_fid); + Field qc3(qc3_fid); + + qc1.allocate_view(); + qc2.allocate_view(); + qc3.allocate_view(); + + // Construct random number generator stuff + using RPDF = std::uniform_real_distribution; + RPDF pdf(sp(0.0), sp(200.0)); + auto engine = scream::setup_random_test(); + + // Set time for qc and randomize its values + qc1.get_header().get_tracking().update_time_stamp(t0); + qc2.get_header().get_tracking().update_time_stamp(t0); + qc3.get_header().get_tracking().update_time_stamp(t0); + randomize(qc1, engine, pdf); + randomize(qc2, engine, pdf); + randomize(qc3, engine, pdf); + + // Construct the Diagnostics + std::map> diags; + auto &diag_factory = AtmosphereDiagnosticFactory::instance(); + register_diagnostics(); + + // Create and set up the diagnostic + ekat::ParameterList params; + REQUIRE_THROWS(diag_factory.create("HistogramDiag", comm, + params)); // Bad construction + + params.set("grid_name", grid->name()); + REQUIRE_THROWS(diag_factory.create("HistogramDiag", comm, + params)); // Still no field_name + + params.set("field_name", "qc"); + REQUIRE_THROWS(diag_factory.create("HistogramDiag", comm, + params)); // Still no bin configuration + + params.set("bin_configuration", bin_configuration); + // Now we should be good to go... + auto diag1 = diag_factory.create("HistogramDiag", comm, params); + auto diag2 = diag_factory.create("HistogramDiag", comm, params); + auto diag3 = diag_factory.create("HistogramDiag", comm, params); + diag1->set_grids(gm); + diag2->set_grids(gm); + diag3->set_grids(gm); + + // Test the zonal average of qc1 + diag1->set_required_field(qc1); + diag1->initialize(t0, RunType::Initial); + diag1->compute_diagnostic(); + auto diag1_field = diag1->get_diagnostic(); + + // Manual calculation + const int num_bins = bin_values.size()-1; + const std::string bin_dim_name = diag1_field.get_header().get_identifier().get_layout().name(0); + FieldLayout diag0_layout({CMP}, {num_bins}, {bin_dim_name}); + FieldIdentifier diag0_id("qc_histogram_manual", diag0_layout, kg / kg, grid->name(), DataType::IntType); + Field diag0_field(diag0_id); + diag0_field.allocate_view(); + + // calculate the histogram + auto qc1_view_h = qc1.get_view(); + auto diag0_view_h = diag0_field.get_view(); + for (int bin_i = 0; bin_i < num_bins; bin_i++) { + for (int col_i = 0; col_i < ngcols; col_i++) { + if (bin_values[bin_i] <= qc1_view_h(col_i) && qc1_view_h(col_i) < bin_values[bin_i+1]) + diag0_view_h(bin_i) += 1; + } + } + diag0_field.sync_to_dev(); + + // Compare + REQUIRE(views_are_equal(diag1_field, diag0_field)); +/* + // Try other known cases + // Set qc1_v to 1.0 to get zonal averages of 1.0/nlats + const Real zavg1 = sp(1.0); + qc1.deep_copy(zavg1); + diag1->compute_diagnostic(); + auto diag1_view_host = diag1_field.get_view(); + for (int nlat = 0; nlat < nlats; nlat++) { + REQUIRE_THAT(diag1_view_host(nlat), Catch::Matchers::WithinRel(zavg1, tol)); + } + + // other diags + // Set qc2_v to 5.0 to get weighted average of 5.0 + const Real zavg2 = sp(5.0); + qc2.deep_copy(zavg2); + diag2->set_required_field(qc2); + diag2->initialize(t0, RunType::Initial); + diag2->compute_diagnostic(); + auto diag2_field = diag2->get_diagnostic(); + + auto diag2_view_host = diag2_field.get_view(); + for (int i = 0; i < nlevs; ++i) { + for (int nlat = 0; nlat < nlats; nlat++) { + REQUIRE_THAT(diag2_view_host(nlat, i), Catch::Matchers::WithinRel(zavg2, tol)); + } + } + + // Try a random case with qc3 + FieldLayout diag3m_layout({CMP, CMP, LEV}, {nlats, dim3, nlevs}, + {bin_dim_name, e2str(CMP), e2str(LEV)}); + FieldIdentifier diag3m_id("qc_zonal_avg_manual", diag3m_layout, kg / kg, grid->name()); + Field diag3m_field(diag3m_id); + diag3m_field.allocate_view(); + auto qc3_view_h = qc3.get_view(); + auto diag3m_view_h = diag3m_field.get_view(); + for (int i = 0; i < ngcols; i++) { + const int nlat = i % nlats; + for (int j = 0; j < dim3; j++) { + for (int k = 0; k < nlevs; k++) { + diag3m_view_h(nlat, j, k) += area_view_h(i) / zonal_areas[nlat] * qc3_view_h(i, j, k); + } + } + } + diag3m_field.sync_to_dev(); + diag3->set_required_field(qc3); + diag3->initialize(t0, RunType::Initial); + diag3->compute_diagnostic(); + auto diag3_field = diag3->get_diagnostic(); + REQUIRE(views_are_equal(diag3_field, diag3m_field)); + */ +} + +} // namespace scream From ab4e2e57a803bb46356d57d4ad53e50476c11e2b Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Fri, 11 Jul 2025 12:36:16 -0700 Subject: [PATCH 02/15] histogram diagnostic now supports up to rank 3 and passes unit tests --- .../eamxx/src/diagnostics/histogram.cpp | 94 +++++++++++++++---- .../src/diagnostics/tests/histogram_test.cpp | 54 +++++------ 2 files changed, 104 insertions(+), 44 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index 90c163060cd4..62a32ac19076 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -22,6 +22,14 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { using ShortFieldTagsNames::CMP; const Field &field = get_fields_in().front(); const FieldIdentifier &field_id = field.get_header().get_identifier(); + const FieldLayout &field_layout = field_id.get_layout(); + EKAT_REQUIRE_MSG(field_layout.rank() >= 1 && field_layout.rank() <= 3, + "Error! Field rank not supported by HistogramDiag.\n" + " - field name: " + + field_id.name() + + "\n" + " - field layout: " + + field_layout.to_string() + "\n"); // extract bin values from configuration std::istringstream stream(m_bin_configuration); @@ -34,7 +42,8 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { // allocate histogram field FieldLayout diagnostic_layout({CMP}, {num_bins}, {"bin"}); - FieldIdentifier diagnostic_id(m_diag_name, diagnostic_layout, field_id.get_units(), + FieldIdentifier diagnostic_id(m_diag_name, diagnostic_layout, + FieldIdentifier::Units::nondimensional(), field_id.get_grid_name(), DataType::IntType); m_diagnostic_output = Field(diagnostic_id); m_diagnostic_output.allocate_view(); @@ -58,28 +67,81 @@ void HistogramDiag::compute_diagnostic_impl() { auto field_layout = field.get_header().get_identifier().get_layout(); auto histogram_layout = m_diagnostic_output.get_header().get_identifier().get_layout(); const int num_bins = histogram_layout.dim(0); - const int field_size = field_layout.size(); auto bin_values_view = m_bin_values.get_view(); auto histogram_view = m_diagnostic_output.get_view(); - auto field_view = field.get_view(); using KT = ekat::KokkosTypes; using TeamPolicy = Kokkos::TeamPolicy; using TeamMember = typename TeamPolicy::member_type; using ESU = ekat::ExeSpaceUtils; - TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, field_size); - Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, - KOKKOS_LAMBDA(const TeamMember &tm) { - const int bin_i = tm.league_rank(); - const Real bin_lower = bin_values_view(bin_i); - const Real bin_upper = bin_values_view(bin_i+1); - Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, field_size), - [&](int i, Int &val) { - if ((bin_lower <= field_view(i)) && (field_view(i) < bin_upper)) - val++; - }, - histogram_view(bin_i)); - }); + switch (field_layout.rank()) + { + case 1: { + const int d1 = field_layout.dim(0); + auto field_view = field.get_view(); + TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, d1); + Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, + KOKKOS_LAMBDA(const TeamMember &tm) { + const int bin_i = tm.league_rank(); + const Real bin_lower = bin_values_view(bin_i); + const Real bin_upper = bin_values_view(bin_i+1); + Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1), + [&](int i, Int &val) { + if ((bin_lower <= field_view(i)) && (field_view(i) < bin_upper)) + val++; + }, + histogram_view(bin_i)); + }); + } break; + case 2: { + const int d1 = field_layout.dim(0); + const int d2 = field_layout.dim(1); + auto field_view = field.get_view(); + TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, d1*d2); + Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, + KOKKOS_LAMBDA(const TeamMember &tm) { + const int bin_i = tm.league_rank(); + const Real bin_lower = bin_values_view(bin_i); + const Real bin_upper = bin_values_view(bin_i+1); + histogram_view(bin_i) = 0; + Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1*d2), + [&](int ind, Int &val) { + const int i1 = ind / d2; + const int i2 = ind % d2; + if ((bin_lower <= field_view(i1,i2)) && (field_view(i1,i2) < bin_upper)) + val += 1; + }, + histogram_view(bin_i)); + }); + } break; + case 3: { + const int d1 = field_layout.dim(0); + const int d2 = field_layout.dim(1); + const int d3 = field_layout.dim(2); + auto field_view = field.get_view(); + TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, d1*d2*d3); + Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, + KOKKOS_LAMBDA(const TeamMember &tm) { + const int bin_i = tm.league_rank(); + const Real bin_lower = bin_values_view(bin_i); + const Real bin_upper = bin_values_view(bin_i+1); + histogram_view(bin_i) = 0; + Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1*d2*d3), + [&](int ind, Int &val) { + const int i1 = ind / (d2*d3); + const int ind2 = ind % (d2*d3); + const int i2 = ind2 / d3; + const int i3 = ind2 % d3; + if ((bin_lower <= field_view(i1,i2,i3)) && (field_view(i1,i2,i3) < bin_upper)) + val += 1; + }, + histogram_view(bin_i)); + }); + } break; + + default: + EKAT_ERROR_MSG("Error! Unsupported field rank for histogram.\n"); + } // TODO: use device-side MPI calls // TODO: the dev ptr causes problems; revisit this later diff --git a/components/eamxx/src/diagnostics/tests/histogram_test.cpp b/components/eamxx/src/diagnostics/tests/histogram_test.cpp index 5901fc7b0f5a..a7656651bb2c 100644 --- a/components/eamxx/src/diagnostics/tests/histogram_test.cpp +++ b/components/eamxx/src/diagnostics/tests/histogram_test.cpp @@ -29,9 +29,6 @@ TEST_CASE("histogram") { using namespace ShortFieldTagsNames; using namespace ekat::units; - // A numerical tolerance - auto tol = std::numeric_limits::epsilon() * 100; - // A world comm ekat::Comm comm(MPI_COMM_WORLD); @@ -117,7 +114,8 @@ TEST_CASE("histogram") { const int num_bins = bin_values.size()-1; const std::string bin_dim_name = diag1_field.get_header().get_identifier().get_layout().name(0); FieldLayout diag0_layout({CMP}, {num_bins}, {bin_dim_name}); - FieldIdentifier diag0_id("qc_histogram_manual", diag0_layout, kg / kg, grid->name(), DataType::IntType); + FieldIdentifier diag0_id("qc_histogram_manual", diag0_layout, + FieldIdentifier::Units::nondimensional(), grid->name(), DataType::IntType); Field diag0_field(diag0_id); diag0_field.allocate_view(); @@ -134,46 +132,47 @@ TEST_CASE("histogram") { // Compare REQUIRE(views_are_equal(diag1_field, diag0_field)); -/* + // Try other known cases - // Set qc1_v to 1.0 to get zonal averages of 1.0/nlats - const Real zavg1 = sp(1.0); + // Set qc1_v to so histogram is all entries in first bin + const Real zavg1 = sp(0.5*(bin_values[0]+bin_values[1])); qc1.deep_copy(zavg1); diag1->compute_diagnostic(); - auto diag1_view_host = diag1_field.get_view(); - for (int nlat = 0; nlat < nlats; nlat++) { - REQUIRE_THAT(diag1_view_host(nlat), Catch::Matchers::WithinRel(zavg1, tol)); + auto diag1_view_host = diag1_field.get_view(); + REQUIRE(diag1_view_host(0) == qc1.get_header().get_identifier().get_layout().size()); + for (int bin_i = 1; bin_i < num_bins; bin_i++) { + REQUIRE(diag1_view_host(bin_i) == 0); } // other diags - // Set qc2_v to 5.0 to get weighted average of 5.0 - const Real zavg2 = sp(5.0); + // Set qc2_v so histogram is all entries in last bin + const Real zavg2 = sp(0.5*(bin_values[num_bins-1]+bin_values[num_bins])); qc2.deep_copy(zavg2); diag2->set_required_field(qc2); diag2->initialize(t0, RunType::Initial); diag2->compute_diagnostic(); auto diag2_field = diag2->get_diagnostic(); - - auto diag2_view_host = diag2_field.get_view(); - for (int i = 0; i < nlevs; ++i) { - for (int nlat = 0; nlat < nlats; nlat++) { - REQUIRE_THAT(diag2_view_host(nlat, i), Catch::Matchers::WithinRel(zavg2, tol)); - } + auto diag2_view_host = diag2_field.get_view(); + REQUIRE(diag2_view_host(num_bins-1) == qc2.get_header().get_identifier().get_layout().size()); + for (int bin_i = num_bins-2; bin_i >=0; bin_i--) { + REQUIRE(diag2_view_host(bin_i) == 0); } // Try a random case with qc3 - FieldLayout diag3m_layout({CMP, CMP, LEV}, {nlats, dim3, nlevs}, - {bin_dim_name, e2str(CMP), e2str(LEV)}); - FieldIdentifier diag3m_id("qc_zonal_avg_manual", diag3m_layout, kg / kg, grid->name()); + FieldLayout diag3m_layout({CMP}, {num_bins}, {bin_dim_name}); + FieldIdentifier diag3m_id("qc_zonal_avg_manual", diag3m_layout, + FieldIdentifier::Units::nondimensional(), grid->name(), DataType::IntType); Field diag3m_field(diag3m_id); diag3m_field.allocate_view(); auto qc3_view_h = qc3.get_view(); - auto diag3m_view_h = diag3m_field.get_view(); - for (int i = 0; i < ngcols; i++) { - const int nlat = i % nlats; - for (int j = 0; j < dim3; j++) { - for (int k = 0; k < nlevs; k++) { - diag3m_view_h(nlat, j, k) += area_view_h(i) / zonal_areas[nlat] * qc3_view_h(i, j, k); + auto diag3m_view_h = diag3m_field.get_view(); + for (int bin_i = 0; bin_i < num_bins; bin_i++) { + for (int i = 0; i < ngcols; i++) { + for (int j = 0; j < dim3; j++) { + for (int k = 0; k < nlevs; k++) { + if (bin_values[bin_i] <= qc3_view_h(i,j,k) && qc3_view_h(i,j,k) < bin_values[bin_i+1]) + diag3m_view_h(bin_i) += 1; + } } } } @@ -183,7 +182,6 @@ TEST_CASE("histogram") { diag3->compute_diagnostic(); auto diag3_field = diag3->get_diagnostic(); REQUIRE(views_are_equal(diag3_field, diag3m_field)); - */ } } // namespace scream From 283fdafda7d0a6a89cbca78ac552c3a5aa933979 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 21 Jul 2025 15:39:26 -0700 Subject: [PATCH 03/15] updated histogram diag test to support multiple MPI ranks --- .../eamxx/src/diagnostics/histogram.cpp | 2 +- .../src/diagnostics/tests/histogram_test.cpp | 28 ++++++++++--------- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index 62a32ac19076..afc0c487838e 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -38,7 +38,7 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { while (std::getline(stream, token, '_')) { bin_values.push_back(std::stod(token)); } - int num_bins = bin_values.size()-1; + const int num_bins = bin_values.size()-1; // allocate histogram field FieldLayout diagnostic_layout({CMP}, {num_bins}, {"bin"}); diff --git a/components/eamxx/src/diagnostics/tests/histogram_test.cpp b/components/eamxx/src/diagnostics/tests/histogram_test.cpp index a7656651bb2c..70915407869e 100644 --- a/components/eamxx/src/diagnostics/tests/histogram_test.cpp +++ b/components/eamxx/src/diagnostics/tests/histogram_test.cpp @@ -7,8 +7,7 @@ namespace scream { -std::shared_ptr create_gm(const ekat::Comm &comm, const int ncols, const int nlevs) { - const int num_global_cols = ncols * comm.size(); +std::shared_ptr create_gm(const ekat::Comm &comm, const int ngcols, const int nlevs) { using vos_t = std::vector; ekat::ParameterList gm_params; @@ -16,7 +15,7 @@ std::shared_ptr create_gm(const ekat::Comm &comm, const int ncols, auto &pl = gm_params.sublist("Point Grid"); pl.set("type", "point_grid"); pl.set("aliases", vos_t{"Physics"}); - pl.set("number_of_global_columns", num_global_cols); + pl.set("number_of_global_columns", ngcols); pl.set("number_of_vertical_levels", nlevs); auto gm = create_mesh_free_grids_manager(comm, gm_params); @@ -38,19 +37,20 @@ TEST_CASE("histogram") { // Create a grids manager - single column for these tests constexpr int nlevs = 3; constexpr int dim3 = 4; - const int ngcols = 6 * comm.size(); + const int ncols = 6; - auto gm = create_gm(comm, ngcols, nlevs); - auto grid = gm->get_grid("Physics"); + const int ngcols = ncols * comm.size(); + auto gm = create_gm(comm, ngcols, nlevs); + auto grid = gm->get_grid("Physics"); // Specify histogram bins const std::string bin_configuration = "0_50_100_150_200"; const std::vector bin_values = {0.0, 50.0, 100.0, 150.0, 200.0}; // Input (randomized) qc - FieldLayout scalar1d_layout{{COL}, {ngcols}}; - FieldLayout scalar2d_layout{{COL, LEV}, {ngcols, nlevs}}; - FieldLayout scalar3d_layout{{COL, CMP, LEV}, {ngcols, dim3, nlevs}}; + FieldLayout scalar1d_layout{{COL}, {ncols}}; + FieldLayout scalar2d_layout{{COL, LEV}, {ncols, nlevs}}; + FieldLayout scalar3d_layout{{COL, CMP, LEV}, {ncols, dim3, nlevs}}; FieldIdentifier qc1_id("qc", scalar1d_layout, kg / kg, grid->name()); FieldIdentifier qc2_fid("qc", scalar2d_layout, kg / kg, grid->name()); @@ -123,11 +123,12 @@ TEST_CASE("histogram") { auto qc1_view_h = qc1.get_view(); auto diag0_view_h = diag0_field.get_view(); for (int bin_i = 0; bin_i < num_bins; bin_i++) { - for (int col_i = 0; col_i < ngcols; col_i++) { + for (int col_i = 0; col_i < ncols; col_i++) { if (bin_values[bin_i] <= qc1_view_h(col_i) && qc1_view_h(col_i) < bin_values[bin_i+1]) diag0_view_h(bin_i) += 1; } } + comm.all_reduce(diag0_field.template get_internal_view_data(), diag0_layout.size(), MPI_SUM); diag0_field.sync_to_dev(); // Compare @@ -139,7 +140,7 @@ TEST_CASE("histogram") { qc1.deep_copy(zavg1); diag1->compute_diagnostic(); auto diag1_view_host = diag1_field.get_view(); - REQUIRE(diag1_view_host(0) == qc1.get_header().get_identifier().get_layout().size()); + REQUIRE(diag1_view_host(0) == ngcols); for (int bin_i = 1; bin_i < num_bins; bin_i++) { REQUIRE(diag1_view_host(bin_i) == 0); } @@ -153,7 +154,7 @@ TEST_CASE("histogram") { diag2->compute_diagnostic(); auto diag2_field = diag2->get_diagnostic(); auto diag2_view_host = diag2_field.get_view(); - REQUIRE(diag2_view_host(num_bins-1) == qc2.get_header().get_identifier().get_layout().size()); + REQUIRE(diag2_view_host(num_bins-1) == ngcols*nlevs); for (int bin_i = num_bins-2; bin_i >=0; bin_i--) { REQUIRE(diag2_view_host(bin_i) == 0); } @@ -167,7 +168,7 @@ TEST_CASE("histogram") { auto qc3_view_h = qc3.get_view(); auto diag3m_view_h = diag3m_field.get_view(); for (int bin_i = 0; bin_i < num_bins; bin_i++) { - for (int i = 0; i < ngcols; i++) { + for (int i = 0; i < ncols; i++) { for (int j = 0; j < dim3; j++) { for (int k = 0; k < nlevs; k++) { if (bin_values[bin_i] <= qc3_view_h(i,j,k) && qc3_view_h(i,j,k) < bin_values[bin_i+1]) @@ -176,6 +177,7 @@ TEST_CASE("histogram") { } } } + comm.all_reduce(diag3m_field.template get_internal_view_data(), diag3m_layout.size(), MPI_SUM); diag3m_field.sync_to_dev(); diag3->set_required_field(qc3); diag3->initialize(t0, RunType::Initial); From 16d27a2ad929484c6638ffff6084d1be5929f123 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 21 Jul 2025 15:40:09 -0700 Subject: [PATCH 04/15] added regex to eamxx_io_utils for histogram diag --- components/eamxx/src/diagnostics/tests/histogram_test.cpp | 7 ++++--- components/eamxx/src/share/io/eamxx_io_utils.cpp | 8 +++++++- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/components/eamxx/src/diagnostics/tests/histogram_test.cpp b/components/eamxx/src/diagnostics/tests/histogram_test.cpp index 70915407869e..8828ed195466 100644 --- a/components/eamxx/src/diagnostics/tests/histogram_test.cpp +++ b/components/eamxx/src/diagnostics/tests/histogram_test.cpp @@ -8,7 +8,6 @@ namespace scream { std::shared_ptr create_gm(const ekat::Comm &comm, const int ngcols, const int nlevs) { - using vos_t = std::vector; ekat::ParameterList gm_params; gm_params.set("grids_names", vos_t{"Point Grid"}); @@ -128,7 +127,8 @@ TEST_CASE("histogram") { diag0_view_h(bin_i) += 1; } } - comm.all_reduce(diag0_field.template get_internal_view_data(), diag0_layout.size(), MPI_SUM); + comm.all_reduce(diag0_field.template get_internal_view_data(), + diag0_layout.size(), MPI_SUM); diag0_field.sync_to_dev(); // Compare @@ -177,7 +177,8 @@ TEST_CASE("histogram") { } } } - comm.all_reduce(diag3m_field.template get_internal_view_data(), diag3m_layout.size(), MPI_SUM); + comm.all_reduce(diag3m_field.template get_internal_view_data(), + diag3m_layout.size(), MPI_SUM); diag3m_field.sync_to_dev(); diag3->set_required_field(qc3); diag3->initialize(t0, RunType::Initial); diff --git a/components/eamxx/src/share/io/eamxx_io_utils.cpp b/components/eamxx/src/share/io/eamxx_io_utils.cpp index 85074ef8f8ba..d93e34e80f62 100644 --- a/components/eamxx/src/share/io/eamxx_io_utils.cpp +++ b/components/eamxx/src/share/io/eamxx_io_utils.cpp @@ -131,7 +131,7 @@ create_diagnostic (const std::string& diag_field_name, // Note: the number for field_at_p/h can match positive integer/floating-point numbers // Start with a generic for a field name allowing for all letters, all numbers, dash, dot, plus, minus, product, and division // Escaping all the special ones just in case - std::string generic_field = "([A-Za-z0-9_.+\\-\\*\\÷]+)"; + std::string generic_field = "([A-Za-z0-9_.+\\-\\*\\÷]+)"; std::regex field_at_l (R"()" + generic_field + R"(_at_(lev_(\d+)|model_(top|bot))$)"); std::regex field_at_p (R"()" + generic_field + R"(_at_(\d+(\.\d+)?)(hPa|mb|Pa)$)"); std::regex field_at_h (R"()" + generic_field + R"(_at_(\d+(\.\d+)?)(m)_above_(sealevel|surface)$)"); @@ -149,6 +149,7 @@ create_diagnostic (const std::string& diag_field_name, std::regex conditional_sampling (R"()" + generic_field + R"(_where_)" + generic_field + R"(_(gt|ge|eq|ne|le|lt)_([+-]?\d+(?:\.\d+)?)$)"); std::regex binary_ops (generic_field + "_" "(plus|minus|times|over)" + "_" + generic_field + "$"); std::regex vert_derivative (generic_field + "_(p|z)vert_derivative$"); + std::regex histogram (R"()" + generic_field + R"(_histogram_(\d+(\.\d+)?(_\d+(\.\d+)?)+)$)"); std::string diag_name; std::smatch matches; @@ -246,6 +247,11 @@ create_diagnostic (const std::string& diag_field_name, params.set("condition_field", matches[2].str()); params.set("condition_operator", matches[3].str()); params.set("condition_value", matches[4].str()); + else if (std::regex_search(diag_field_name,matches,histogram)) { + diag_name = "HistogramDiag"; + params.set("grid_name", grid->name()); + params.set("field_name", matches[1].str()); + params.set("bin_configuration", matches[2].str()); } else if (std::regex_search(diag_field_name,matches,binary_ops)) { diag_name = "BinaryOpsDiag"; From ece25fd1981408b674b882d4c22f866e15c13353 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 21 Jul 2025 18:24:07 -0700 Subject: [PATCH 05/15] converted histogram diag to use Real instead of Int --- .../eamxx/src/diagnostics/histogram.cpp | 20 +++++------ .../src/diagnostics/tests/histogram_test.cpp | 34 +++++++++++-------- 2 files changed, 29 insertions(+), 25 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index afc0c487838e..209a8ebf12e0 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -2,7 +2,6 @@ namespace scream { - HistogramDiag::HistogramDiag(const ekat::Comm &comm, const ekat::ParameterList ¶ms) : AtmosphereDiagnostic(comm, params) { const auto &field_name = m_params.get("field_name"); @@ -44,7 +43,7 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { FieldLayout diagnostic_layout({CMP}, {num_bins}, {"bin"}); FieldIdentifier diagnostic_id(m_diag_name, diagnostic_layout, FieldIdentifier::Units::nondimensional(), - field_id.get_grid_name(), DataType::IntType); + field_id.get_grid_name()); m_diagnostic_output = Field(diagnostic_id); m_diagnostic_output.allocate_view(); @@ -68,8 +67,9 @@ void HistogramDiag::compute_diagnostic_impl() { auto histogram_layout = m_diagnostic_output.get_header().get_identifier().get_layout(); const int num_bins = histogram_layout.dim(0); + m_diagnostic_output.deep_copy(sp(0.0)); auto bin_values_view = m_bin_values.get_view(); - auto histogram_view = m_diagnostic_output.get_view(); + auto histogram_view = m_diagnostic_output.get_view(); using KT = ekat::KokkosTypes; using TeamPolicy = Kokkos::TeamPolicy; using TeamMember = typename TeamPolicy::member_type; @@ -86,9 +86,9 @@ void HistogramDiag::compute_diagnostic_impl() { const Real bin_lower = bin_values_view(bin_i); const Real bin_upper = bin_values_view(bin_i+1); Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1), - [&](int i, Int &val) { + [&](int i, Real &val) { if ((bin_lower <= field_view(i)) && (field_view(i) < bin_upper)) - val++; + val += sp(1.0); }, histogram_view(bin_i)); }); @@ -105,11 +105,11 @@ void HistogramDiag::compute_diagnostic_impl() { const Real bin_upper = bin_values_view(bin_i+1); histogram_view(bin_i) = 0; Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1*d2), - [&](int ind, Int &val) { + [&](int ind, Real &val) { const int i1 = ind / d2; const int i2 = ind % d2; if ((bin_lower <= field_view(i1,i2)) && (field_view(i1,i2) < bin_upper)) - val += 1; + val += sp(1.0); }, histogram_view(bin_i)); }); @@ -127,13 +127,13 @@ void HistogramDiag::compute_diagnostic_impl() { const Real bin_upper = bin_values_view(bin_i+1); histogram_view(bin_i) = 0; Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1*d2*d3), - [&](int ind, Int &val) { + [&](int ind, Real &val) { const int i1 = ind / (d2*d3); const int ind2 = ind % (d2*d3); const int i2 = ind2 / d3; const int i3 = ind2 % d3; if ((bin_lower <= field_view(i1,i2,i3)) && (field_view(i1,i2,i3) < bin_upper)) - val += 1; + val += sp(1.0); }, histogram_view(bin_i)); }); @@ -148,7 +148,7 @@ void HistogramDiag::compute_diagnostic_impl() { // TODO: doing cuda-aware MPI allreduce would be ~10% faster Kokkos::fence(); m_diagnostic_output.sync_to_host(); - m_comm.all_reduce(m_diagnostic_output.template get_internal_view_data(), + m_comm.all_reduce(m_diagnostic_output.template get_internal_view_data(), histogram_layout.size(), MPI_SUM); m_diagnostic_output.sync_to_dev(); } diff --git a/components/eamxx/src/diagnostics/tests/histogram_test.cpp b/components/eamxx/src/diagnostics/tests/histogram_test.cpp index 8828ed195466..509782757b22 100644 --- a/components/eamxx/src/diagnostics/tests/histogram_test.cpp +++ b/components/eamxx/src/diagnostics/tests/histogram_test.cpp @@ -27,6 +27,9 @@ TEST_CASE("histogram") { using namespace ShortFieldTagsNames; using namespace ekat::units; + // A numerical tolerance + const auto tol = std::numeric_limits::epsilon() * 100; + // A world comm ekat::Comm comm(MPI_COMM_WORLD); @@ -114,20 +117,21 @@ TEST_CASE("histogram") { const std::string bin_dim_name = diag1_field.get_header().get_identifier().get_layout().name(0); FieldLayout diag0_layout({CMP}, {num_bins}, {bin_dim_name}); FieldIdentifier diag0_id("qc_histogram_manual", diag0_layout, - FieldIdentifier::Units::nondimensional(), grid->name(), DataType::IntType); + FieldIdentifier::Units::nondimensional(), grid->name()); Field diag0_field(diag0_id); diag0_field.allocate_view(); // calculate the histogram + diag0_field.deep_copy(sp(0.0)); auto qc1_view_h = qc1.get_view(); - auto diag0_view_h = diag0_field.get_view(); + auto diag0_view_h = diag0_field.get_view(); for (int bin_i = 0; bin_i < num_bins; bin_i++) { for (int col_i = 0; col_i < ncols; col_i++) { if (bin_values[bin_i] <= qc1_view_h(col_i) && qc1_view_h(col_i) < bin_values[bin_i+1]) - diag0_view_h(bin_i) += 1; + diag0_view_h(bin_i) += sp(1.0); } } - comm.all_reduce(diag0_field.template get_internal_view_data(), + comm.all_reduce(diag0_field.template get_internal_view_data(), diag0_layout.size(), MPI_SUM); diag0_field.sync_to_dev(); @@ -139,10 +143,10 @@ TEST_CASE("histogram") { const Real zavg1 = sp(0.5*(bin_values[0]+bin_values[1])); qc1.deep_copy(zavg1); diag1->compute_diagnostic(); - auto diag1_view_host = diag1_field.get_view(); - REQUIRE(diag1_view_host(0) == ngcols); + auto diag1_view_host = diag1_field.get_view(); + REQUIRE_THAT(diag1_view_host(0), Catch::Matchers::WithinRel(ngcols, tol)); for (int bin_i = 1; bin_i < num_bins; bin_i++) { - REQUIRE(diag1_view_host(bin_i) == 0); + REQUIRE(diag1_view_host(bin_i) == sp(0.0)); } // other diags @@ -153,31 +157,31 @@ TEST_CASE("histogram") { diag2->initialize(t0, RunType::Initial); diag2->compute_diagnostic(); auto diag2_field = diag2->get_diagnostic(); - auto diag2_view_host = diag2_field.get_view(); - REQUIRE(diag2_view_host(num_bins-1) == ngcols*nlevs); + auto diag2_view_host = diag2_field.get_view(); + REQUIRE_THAT(diag2_view_host(num_bins-1), Catch::Matchers::WithinRel(ngcols*nlevs, tol)); for (int bin_i = num_bins-2; bin_i >=0; bin_i--) { - REQUIRE(diag2_view_host(bin_i) == 0); + REQUIRE(diag2_view_host(bin_i) == sp(0.0)); } // Try a random case with qc3 FieldLayout diag3m_layout({CMP}, {num_bins}, {bin_dim_name}); FieldIdentifier diag3m_id("qc_zonal_avg_manual", diag3m_layout, - FieldIdentifier::Units::nondimensional(), grid->name(), DataType::IntType); + FieldIdentifier::Units::nondimensional(), grid->name()); Field diag3m_field(diag3m_id); diag3m_field.allocate_view(); - auto qc3_view_h = qc3.get_view(); - auto diag3m_view_h = diag3m_field.get_view(); + auto qc3_view_h = qc3.get_view(); + auto diag3m_view_h = diag3m_field.get_view(); for (int bin_i = 0; bin_i < num_bins; bin_i++) { for (int i = 0; i < ncols; i++) { for (int j = 0; j < dim3; j++) { for (int k = 0; k < nlevs; k++) { if (bin_values[bin_i] <= qc3_view_h(i,j,k) && qc3_view_h(i,j,k) < bin_values[bin_i+1]) - diag3m_view_h(bin_i) += 1; + diag3m_view_h(bin_i) += sp(1.0); } } } } - comm.all_reduce(diag3m_field.template get_internal_view_data(), + comm.all_reduce(diag3m_field.template get_internal_view_data(), diag3m_layout.size(), MPI_SUM); diag3m_field.sync_to_dev(); diag3->set_required_field(qc3); From ab483e124c0080f4472fb304b79759d5ff064c90 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 18 Aug 2025 13:57:43 -0700 Subject: [PATCH 06/15] addressed some issues in histogram after rebase --- components/eamxx/src/diagnostics/histogram.cpp | 11 +++++++---- components/eamxx/src/share/io/eamxx_io_utils.cpp | 11 ++++++----- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index 209a8ebf12e0..243c016f2bcd 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -1,5 +1,8 @@ #include "diagnostics/histogram.hpp" +#include + + namespace scream { HistogramDiag::HistogramDiag(const ekat::Comm &comm, const ekat::ParameterList ¶ms) @@ -73,13 +76,13 @@ void HistogramDiag::compute_diagnostic_impl() { using KT = ekat::KokkosTypes; using TeamPolicy = Kokkos::TeamPolicy; using TeamMember = typename TeamPolicy::member_type; - using ESU = ekat::ExeSpaceUtils; + using TPF = ekat::TeamPolicyFactory; switch (field_layout.rank()) { case 1: { const int d1 = field_layout.dim(0); auto field_view = field.get_view(); - TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, d1); + TeamPolicy team_policy = TPF::get_default_team_policy(num_bins, d1); Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, KOKKOS_LAMBDA(const TeamMember &tm) { const int bin_i = tm.league_rank(); @@ -97,7 +100,7 @@ void HistogramDiag::compute_diagnostic_impl() { const int d1 = field_layout.dim(0); const int d2 = field_layout.dim(1); auto field_view = field.get_view(); - TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, d1*d2); + TeamPolicy team_policy = TPF::get_default_team_policy(num_bins, d1*d2); Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, KOKKOS_LAMBDA(const TeamMember &tm) { const int bin_i = tm.league_rank(); @@ -119,7 +122,7 @@ void HistogramDiag::compute_diagnostic_impl() { const int d2 = field_layout.dim(1); const int d3 = field_layout.dim(2); auto field_view = field.get_view(); - TeamPolicy team_policy = ESU::get_default_team_policy(num_bins, d1*d2*d3); + TeamPolicy team_policy = TPF::get_default_team_policy(num_bins, d1*d2*d3); Kokkos::parallel_for("compute_histogram_" + field.name(), team_policy, KOKKOS_LAMBDA(const TeamMember &tm) { const int bin_i = tm.league_rank(); diff --git a/components/eamxx/src/share/io/eamxx_io_utils.cpp b/components/eamxx/src/share/io/eamxx_io_utils.cpp index d93e34e80f62..72ee9a032713 100644 --- a/components/eamxx/src/share/io/eamxx_io_utils.cpp +++ b/components/eamxx/src/share/io/eamxx_io_utils.cpp @@ -247,11 +247,6 @@ create_diagnostic (const std::string& diag_field_name, params.set("condition_field", matches[2].str()); params.set("condition_operator", matches[3].str()); params.set("condition_value", matches[4].str()); - else if (std::regex_search(diag_field_name,matches,histogram)) { - diag_name = "HistogramDiag"; - params.set("grid_name", grid->name()); - params.set("field_name", matches[1].str()); - params.set("bin_configuration", matches[2].str()); } else if (std::regex_search(diag_field_name,matches,binary_ops)) { diag_name = "BinaryOpsDiag"; @@ -260,6 +255,12 @@ create_diagnostic (const std::string& diag_field_name, params.set("field_2", matches[3].str()); params.set("binary_op", matches[2].str()); } + else if (std::regex_search(diag_field_name,matches,histogram)) { + diag_name = "HistogramDiag"; + params.set("grid_name", grid->name()); + params.set("field_name", matches[1].str()); + params.set("bin_configuration", matches[2].str()); + } else { // No existing special regex matches, so we assume that the diag field name IS the diag name. From 3b748f651bf66d96b17752d29f24d81e7adf483c Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 18 Aug 2025 13:58:04 -0700 Subject: [PATCH 07/15] added documentation for histogram diagnostic --- .../eamxx/docs/user/diags/field_contraction.md | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/components/eamxx/docs/user/diags/field_contraction.md b/components/eamxx/docs/user/diags/field_contraction.md index 33ad504c8b41..ec7efbd1a18c 100644 --- a/components/eamxx/docs/user/diags/field_contraction.md +++ b/components/eamxx/docs/user/diags/field_contraction.md @@ -3,7 +3,7 @@ In EAMxx, we can automatically calculate field reductions across the horizontal columns and across the model vertical levels. We call these horizontal and vertical reductions. -We can also automatically calculate zonal averages. +We can also automatically calculate zonal averages and histograms. We have *experimental* support for composing diagnostics; below, the horizontal and vertical reductions can be composed sequentially, but if using dp- or dz-weighted reductions, @@ -85,6 +85,19 @@ And so on... | --------- | ------ | ----------- | | `X_zonal_avg_Y_bins` | Area fraction | Average across the zonal direction | +## Histograms + +We currently have a utility to calculate histograms online. +To select the histogram, you need to suffix a field name `X` with +`_histogram_` followed by the values specifying the edges of the bins, +separated by `_`. For example, the histogram specified by +`T_mid_histogram_250_270_290_310` has 3 bins defined by +[250, 270), [270, 290), and [290, 310]. + +| Reduction | Weight | Description | +| --------- | ------ | ----------- | +| `X_histogram_V0_V1_..._VN` | 1 or 0 | Count of field values within each range | + ## Example ```yaml @@ -99,7 +112,7 @@ fields: # in this example, we use T_mid in units of K - T_mid_horiz_avg # K - T_mid_vert_avg_dp_weighted # K - - T_mid_vert_sum_dp_weighted # K * Pa * s / (m * m) + - T_mid_vert_sum_dp_weighted # K * Pa * s / (m * m) - T_mid_vert_avg_dz_weighted # K - T_mid_vert_sum_dz_weighted # K * m - T_mid_vert_avg # K @@ -110,6 +123,7 @@ fields: - T_mid_vert_avg_dp_weighted_horiz_avg # K - T_mid_zonal_avg_180_bins # K - T_mid_zonal_avg_90_bins # K + - T_mid_histogram_250_270_290_310 # count output_control: frequency: 1 frequency_units: nmonths From 06f7173c1ede6cf659a8efc2d0c1f1cb26c2f20a Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Thu, 21 Aug 2025 16:24:19 -0700 Subject: [PATCH 08/15] removed unnecessary initializations in histogram diagnostic --- components/eamxx/src/diagnostics/histogram.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index 243c016f2bcd..18cee9b6b8bd 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -70,7 +70,6 @@ void HistogramDiag::compute_diagnostic_impl() { auto histogram_layout = m_diagnostic_output.get_header().get_identifier().get_layout(); const int num_bins = histogram_layout.dim(0); - m_diagnostic_output.deep_copy(sp(0.0)); auto bin_values_view = m_bin_values.get_view(); auto histogram_view = m_diagnostic_output.get_view(); using KT = ekat::KokkosTypes; @@ -106,7 +105,6 @@ void HistogramDiag::compute_diagnostic_impl() { const int bin_i = tm.league_rank(); const Real bin_lower = bin_values_view(bin_i); const Real bin_upper = bin_values_view(bin_i+1); - histogram_view(bin_i) = 0; Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1*d2), [&](int ind, Real &val) { const int i1 = ind / d2; @@ -128,7 +126,6 @@ void HistogramDiag::compute_diagnostic_impl() { const int bin_i = tm.league_rank(); const Real bin_lower = bin_values_view(bin_i); const Real bin_upper = bin_values_view(bin_i+1); - histogram_view(bin_i) = 0; Kokkos::parallel_reduce(Kokkos::TeamVectorRange(tm, d1*d2*d3), [&](int ind, Real &val) { const int i1 = ind / (d2*d3); From 69f3d952e240da357f426ce169e061d08bd4ec16 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 8 Sep 2025 13:01:36 -0700 Subject: [PATCH 09/15] addressed missing includes --- .../ml_correction/eamxx_ml_correction_process_interface.cpp | 1 + .../single-process/ml_correction/ml_correction_standalone.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/components/eamxx/src/physics/ml_correction/eamxx_ml_correction_process_interface.cpp b/components/eamxx/src/physics/ml_correction/eamxx_ml_correction_process_interface.cpp index 79dbbd07cfde..e3f182070b36 100644 --- a/components/eamxx/src/physics/ml_correction/eamxx_ml_correction_process_interface.cpp +++ b/components/eamxx/src/physics/ml_correction/eamxx_ml_correction_process_interface.cpp @@ -5,6 +5,7 @@ #include "share/field/field_utils.hpp" #include +#include #include #include #include diff --git a/components/eamxx/tests/single-process/ml_correction/ml_correction_standalone.cpp b/components/eamxx/tests/single-process/ml_correction/ml_correction_standalone.cpp index 071c30e2a521..2004d8cfea9d 100644 --- a/components/eamxx/tests/single-process/ml_correction/ml_correction_standalone.cpp +++ b/components/eamxx/tests/single-process/ml_correction/ml_correction_standalone.cpp @@ -4,6 +4,7 @@ #include "physics/register_physics.hpp" #include "share/grid/mesh_free_grids_manager.hpp" +#include #include #include From 837c068a356ba202348017e3730088dd0d26a870 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 8 Sep 2025 13:02:11 -0700 Subject: [PATCH 10/15] leveraged ekat::split utility in histogram diagnostic --- components/eamxx/src/diagnostics/histogram.cpp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index 18cee9b6b8bd..7eb80896f274 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -1,6 +1,7 @@ #include "diagnostics/histogram.hpp" #include +#include namespace scream { @@ -34,13 +35,8 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { field_layout.to_string() + "\n"); // extract bin values from configuration - std::istringstream stream(m_bin_configuration); - std::string token; - std::vector bin_values; - while (std::getline(stream, token, '_')) { - bin_values.push_back(std::stod(token)); - } - const int num_bins = bin_values.size()-1; + std::vector bin_strings = ekat::split(m_bin_configuration, "_"); + const int num_bins = bin_strings.size()-1; // allocate histogram field FieldLayout diagnostic_layout({CMP}, {num_bins}, {"bin"}); @@ -60,7 +56,7 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { using RangePolicy = Kokkos::RangePolicy; Kokkos::parallel_for("store_histogram_bin_values_" + field.name(), RangePolicy(0, bin_values_layout.dim(0)), [&] (int i) { - bin_values_view(i) = bin_values[i]; + bin_values_view(i) = std::stod(bin_strings[i]); }); } From 886d74759000ea7ccea25a1623c4f12270fa9070 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 8 Sep 2025 16:12:41 -0700 Subject: [PATCH 11/15] added infinite extents to histogram diagnostics, along with check for monotonic bin values (test updated too) --- .../eamxx/src/diagnostics/histogram.cpp | 27 +++++++++++++------ .../eamxx/src/diagnostics/histogram.hpp | 3 +-- .../src/diagnostics/tests/histogram_test.cpp | 9 +++++-- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index 7eb80896f274..b77af9e8ccd4 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -8,9 +8,23 @@ namespace scream { HistogramDiag::HistogramDiag(const ekat::Comm &comm, const ekat::ParameterList ¶ms) : AtmosphereDiagnostic(comm, params) { - const auto &field_name = m_params.get("field_name"); - m_bin_configuration = params.get("bin_configuration"); - m_diag_name = field_name + "_histogram_" + m_bin_configuration; + const auto &field_name = m_params.get("field_name"); + const std::string bin_config = m_params.get("bin_configuration"); + m_diag_name = field_name + "_histogram_" + bin_config; + + // extract bin values from configuration, append end values, and check + const std::vector bin_strings = ekat::split(bin_config, "_"); + m_bin_reals.resize(bin_strings.size()+2); + m_bin_reals[0] = std::numeric_limits::lowest(); + for (int i=1; i < m_bin_reals.size()-1; i++) + { + m_bin_reals[i] = std::stod(bin_strings[i-1]); + EKAT_REQUIRE_MSG(m_bin_reals[i] > m_bin_reals[i-1], + "Error! HistogramDiag bin values must be monotonically " + "increasing.\n" + " - bin configuration: " + bin_config + "\n"); + } + m_bin_reals.back() = std::numeric_limits::max(); } void HistogramDiag::set_grids(const std::shared_ptr grids_manager) { @@ -34,11 +48,8 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { " - field layout: " + field_layout.to_string() + "\n"); - // extract bin values from configuration - std::vector bin_strings = ekat::split(m_bin_configuration, "_"); - const int num_bins = bin_strings.size()-1; - // allocate histogram field + const int num_bins = m_bin_reals.size()-1; FieldLayout diagnostic_layout({CMP}, {num_bins}, {"bin"}); FieldIdentifier diagnostic_id(m_diag_name, diagnostic_layout, FieldIdentifier::Units::nondimensional(), @@ -56,7 +67,7 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { using RangePolicy = Kokkos::RangePolicy; Kokkos::parallel_for("store_histogram_bin_values_" + field.name(), RangePolicy(0, bin_values_layout.dim(0)), [&] (int i) { - bin_values_view(i) = std::stod(bin_strings[i]); + bin_values_view(i) = m_bin_reals[i]; }); } diff --git a/components/eamxx/src/diagnostics/histogram.hpp b/components/eamxx/src/diagnostics/histogram.hpp index 2e05543c32bb..dcfdbd8335e9 100644 --- a/components/eamxx/src/diagnostics/histogram.hpp +++ b/components/eamxx/src/diagnostics/histogram.hpp @@ -31,8 +31,7 @@ class HistogramDiag : public AtmosphereDiagnostic { protected: std::string m_diag_name; - std::string m_bin_configuration; - + std::vector m_bin_reals; Field m_bin_values; }; diff --git a/components/eamxx/src/diagnostics/tests/histogram_test.cpp b/components/eamxx/src/diagnostics/tests/histogram_test.cpp index 509782757b22..8c615e947c02 100644 --- a/components/eamxx/src/diagnostics/tests/histogram_test.cpp +++ b/components/eamxx/src/diagnostics/tests/histogram_test.cpp @@ -46,8 +46,8 @@ TEST_CASE("histogram") { auto grid = gm->get_grid("Physics"); // Specify histogram bins - const std::string bin_configuration = "0_50_100_150_200"; - const std::vector bin_values = {0.0, 50.0, 100.0, 150.0, 200.0}; + const std::string bin_configuration = "50_100_150"; + const std::vector bin_values = {-100.0, 50.0, 100.0, 150.0, 1000.0}; // Input (randomized) qc FieldLayout scalar1d_layout{{COL}, {ncols}}; @@ -97,7 +97,12 @@ TEST_CASE("histogram") { REQUIRE_THROWS(diag_factory.create("HistogramDiag", comm, params)); // Still no bin configuration + params.set("bin_configuration", "100_25"); + REQUIRE_THROWS(diag_factory.create("HistogramDiag", comm, + params)); // Non-monotonic bin configuation + params.set("bin_configuration", bin_configuration); + // Now we should be good to go... auto diag1 = diag_factory.create("HistogramDiag", comm, params); auto diag2 = diag_factory.create("HistogramDiag", comm, params); From 07b15e47835b26e263f6bd1f280cb232c80bbd74 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 8 Sep 2025 16:16:30 -0700 Subject: [PATCH 12/15] updated diagnostic documentation for histogram --- components/eamxx/docs/user/diags/field_contraction.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/eamxx/docs/user/diags/field_contraction.md b/components/eamxx/docs/user/diags/field_contraction.md index ec7efbd1a18c..722f9e8ff261 100644 --- a/components/eamxx/docs/user/diags/field_contraction.md +++ b/components/eamxx/docs/user/diags/field_contraction.md @@ -91,8 +91,8 @@ We currently have a utility to calculate histograms online. To select the histogram, you need to suffix a field name `X` with `_histogram_` followed by the values specifying the edges of the bins, separated by `_`. For example, the histogram specified by -`T_mid_histogram_250_270_290_310` has 3 bins defined by -[250, 270), [270, 290), and [290, 310]. +`T_mid_histogram_250_270_290_310` has 5 bins effectively defined by +[$-\infty$, 250), [250, 270), [270, 290), [290, 310), and [310, $\infty$). | Reduction | Weight | Description | | --------- | ------ | ----------- | From 6600ff37b978cfc197e0fea4ba23d76b65160284 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 8 Sep 2025 16:47:09 -0700 Subject: [PATCH 13/15] correcting rebase issue in histogram diagnostics --- components/eamxx/src/share/io/eamxx_io_utils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eamxx/src/share/io/eamxx_io_utils.cpp b/components/eamxx/src/share/io/eamxx_io_utils.cpp index 72ee9a032713..79f80c0406eb 100644 --- a/components/eamxx/src/share/io/eamxx_io_utils.cpp +++ b/components/eamxx/src/share/io/eamxx_io_utils.cpp @@ -148,8 +148,8 @@ create_diagnostic (const std::string& diag_field_name, std::regex zonal_avg (R"()" + generic_field + R"(_zonal_avg_(\d+)_bins$)"); std::regex conditional_sampling (R"()" + generic_field + R"(_where_)" + generic_field + R"(_(gt|ge|eq|ne|le|lt)_([+-]?\d+(?:\.\d+)?)$)"); std::regex binary_ops (generic_field + "_" "(plus|minus|times|over)" + "_" + generic_field + "$"); - std::regex vert_derivative (generic_field + "_(p|z)vert_derivative$"); std::regex histogram (R"()" + generic_field + R"(_histogram_(\d+(\.\d+)?(_\d+(\.\d+)?)+)$)"); + std::regex vert_derivative (generic_field + "_(p|z)vert_derivative$"); std::string diag_name; std::smatch matches; From c479bfe28467cc148d2e9a86ee0ab04b877c45d9 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Fri, 19 Sep 2025 18:19:20 -0700 Subject: [PATCH 14/15] fixed issue in histogram diagnostic where device view was used instead of host view --- components/eamxx/src/diagnostics/histogram.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/components/eamxx/src/diagnostics/histogram.cpp b/components/eamxx/src/diagnostics/histogram.cpp index b77af9e8ccd4..92dd141c98c8 100644 --- a/components/eamxx/src/diagnostics/histogram.cpp +++ b/components/eamxx/src/diagnostics/histogram.cpp @@ -16,7 +16,7 @@ HistogramDiag::HistogramDiag(const ekat::Comm &comm, const ekat::ParameterList & const std::vector bin_strings = ekat::split(bin_config, "_"); m_bin_reals.resize(bin_strings.size()+2); m_bin_reals[0] = std::numeric_limits::lowest(); - for (int i=1; i < m_bin_reals.size()-1; i++) + for (long unsigned int i=1; i < m_bin_reals.size()-1; i++) { m_bin_reals[i] = std::stod(bin_strings[i-1]); EKAT_REQUIRE_MSG(m_bin_reals[i] > m_bin_reals[i-1], @@ -63,12 +63,13 @@ void HistogramDiag::initialize_impl(const RunType /*run_type*/) { field_id.get_units(), field_id.get_grid_name()); m_bin_values = Field(bin_values_id); m_bin_values.allocate_view(); + + // copy bin values into field auto bin_values_view = m_bin_values.get_view(); - using RangePolicy = Kokkos::RangePolicy; - Kokkos::parallel_for("store_histogram_bin_values_" + field.name(), - RangePolicy(0, bin_values_layout.dim(0)), [&] (int i) { - bin_values_view(i) = m_bin_reals[i]; - }); + auto bin_values_view_host = Kokkos::create_mirror_view(bin_values_view); + for (long unsigned int i=0; i < bin_values_layout.dim(0); i++) + bin_values_view_host(i) = m_bin_reals[i]; + Kokkos::deep_copy(bin_values_view,bin_values_view_host); } void HistogramDiag::compute_diagnostic_impl() { From 21a991af8481bc8642db4931bfde3752e2b8ef62 Mon Sep 17 00:00:00 2001 From: Chris Vogl Date: Mon, 22 Sep 2025 11:03:24 -0700 Subject: [PATCH 15/15] fixed header in histogram diag after move of eamxx_setup_random_test --- components/eamxx/src/diagnostics/tests/histogram_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eamxx/src/diagnostics/tests/histogram_test.cpp b/components/eamxx/src/diagnostics/tests/histogram_test.cpp index 8c615e947c02..81074aad3c67 100644 --- a/components/eamxx/src/diagnostics/tests/histogram_test.cpp +++ b/components/eamxx/src/diagnostics/tests/histogram_test.cpp @@ -2,7 +2,7 @@ #include "diagnostics/register_diagnostics.hpp" #include "share/field/field_utils.hpp" #include "share/grid/mesh_free_grids_manager.hpp" -#include "share/util/eamxx_setup_random_test.hpp" +#include "share/core/eamxx_setup_random_test.hpp" #include "share/util/eamxx_universal_constants.hpp" namespace scream {