From c0294b25b136ff7965c6a1587119d401e07c5b3f Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Thu, 2 May 2024 17:29:28 -0700 Subject: [PATCH 01/51] #737: Initial draft of Coastal Formulation interface, SCHISM implementation, and supporting classes --- include/forcing/MeshPointsSelectors.hpp | 16 ++++++ .../coastal/CoastalFormulation.hpp | 52 +++++++++++++++++++ .../coastal/SchismFormulation.hpp | 50 ++++++++++++++++++ 3 files changed, 118 insertions(+) create mode 100644 include/forcing/MeshPointsSelectors.hpp create mode 100644 include/realizations/coastal/CoastalFormulation.hpp create mode 100644 include/realizations/coastal/SchismFormulation.hpp diff --git a/include/forcing/MeshPointsSelectors.hpp b/include/forcing/MeshPointsSelectors.hpp new file mode 100644 index 0000000000..08a39fb36c --- /dev/null +++ b/include/forcing/MeshPointsSelectors.hpp @@ -0,0 +1,16 @@ +#pragma once + +#include +#include +#include + +struct all_points {}; + +struct MeshPointsSelector +{ + std::string variable_name; + std::chrono::time_point init_time; + std::chrono::duration duration; + std::string output_units; + boost::variant> points; +}; diff --git a/include/realizations/coastal/CoastalFormulation.hpp b/include/realizations/coastal/CoastalFormulation.hpp new file mode 100644 index 0000000000..9a4fc791f0 --- /dev/null +++ b/include/realizations/coastal/CoastalFormulation.hpp @@ -0,0 +1,52 @@ +#pragma once + +#include +#include + +#include +#include + +class CoastalFormulation : public data_access::DataProvider +{ +public: + CoastalFormulation(std::string id) + : id_(std::move(id)) + { } + + virtual ~CoastalFormulation() = default; + + std::string get_id() const { + return this->id; + } + + virtual void initialize() = 0; + virtual void finalize() = 0; + virtual void update() = 0; + + // The interface that DataProvider really should have + virtual void get_values(const selection_type& selector, boost::span data) = 0; + + // And an implementation of the usual version using it + std::vector get_values(const selection_type& selector, ReSampleMethod) override + { + std::vector output(selected_points_count(selector)); + get_values(selector, output); + return output; + } + +protected: + using selection_type = MeshPointsSelector; + using data_type = double; + + virtual size_t mesh_size(std::string const& variable_name) = 0; + + size_t selected_points_count(const selection_type& selector) + { + auto* points = get>(selector.points); + size_t size = points ? points->size() : this->mesh_size(selector.variable_name); + return size; + } + +private: + const std::string id_; +}; diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp new file mode 100644 index 0000000000..1809ce91ae --- /dev/null +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include +#include +#include + +class SchismFormulation : public CoastalFormulation +{ +public: + using MeshPointsDataProvider = data_access::DataProvider; + SchismFormulation( + std::string id + , std::shared_ptr met_forcings + , std::shared_ptr offshore_boundary + , std::shared_ptr inflow_boundary + ); + + ~SchismFormulation(); + + // Implementation of DataProvider + boost::span get_available_variable_names() override; + + long get_data_start_time() override; + long get_data_stop_time() override; + long record_duration() override; + size_t get_ts_index_for_time(const time_t &epoch_time) override; + + data_type get_value(const selection_type& selector, ReSampleMethod m) override; + std::vector get_values(const selection_type& selector, ReSampleMethod m) override; + + // Implementation of CoastalFormulation + void initialize() override; + void finalize() override; + void update() override; + + void get_values(const selection_type& selector, boost::span data) override; + +private: + std::unique_ptr bmi_; + + // TODO: Some of these maybe should be members of + // CoastalFormulation + + // TODO: We really want something that we can ask for + // area-averaged RAINRATE over elements, but we're going to makde + // do with point values at the element centroids + std::shared_ptr> meteorological_forcings_provider_; + std::shared_ptr> offshore_boundary_provider_; + std::shared_ptr> inflows_boundary_provider_; +}; From a57659f8b17e66410db4d62e05367eb359cdca49 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 24 Jun 2024 08:03:01 -0700 Subject: [PATCH 02/51] Clean up compiler issues --- include/forcing/MeshPointsSelectors.hpp | 4 ++-- .../realizations/coastal/CoastalFormulation.hpp | 8 ++++---- .../realizations/coastal/SchismFormulation.hpp | 8 ++++---- test/CMakeLists.txt | 15 +++++++++++++++ test/coastal/SchismFormulation_Test.cpp | 9 +++++++++ 5 files changed, 34 insertions(+), 10 deletions(-) create mode 100644 test/coastal/SchismFormulation_Test.cpp diff --git a/include/forcing/MeshPointsSelectors.hpp b/include/forcing/MeshPointsSelectors.hpp index 08a39fb36c..b27fb1ecd7 100644 --- a/include/forcing/MeshPointsSelectors.hpp +++ b/include/forcing/MeshPointsSelectors.hpp @@ -9,8 +9,8 @@ struct all_points {}; struct MeshPointsSelector { std::string variable_name; - std::chrono::time_point init_time; - std::chrono::duration duration; + std::chrono::time_point init_time; + std::chrono::duration duration; std::string output_units; boost::variant> points; }; diff --git a/include/realizations/coastal/CoastalFormulation.hpp b/include/realizations/coastal/CoastalFormulation.hpp index 9a4fc791f0..260ab15140 100644 --- a/include/realizations/coastal/CoastalFormulation.hpp +++ b/include/realizations/coastal/CoastalFormulation.hpp @@ -16,18 +16,18 @@ class CoastalFormulation : public data_access::DataProviderid; + return this->id_; } virtual void initialize() = 0; virtual void finalize() = 0; virtual void update() = 0; - + // The interface that DataProvider really should have virtual void get_values(const selection_type& selector, boost::span data) = 0; // And an implementation of the usual version using it - std::vector get_values(const selection_type& selector, ReSampleMethod) override + std::vector get_values(const selection_type& selector, data_access::ReSampleMethod) override { std::vector output(selected_points_count(selector)); get_values(selector, output); @@ -42,7 +42,7 @@ class CoastalFormulation : public data_access::DataProvider>(selector.points); + auto* points = boost::get>(selector.points); size_t size = points ? points->size() : this->mesh_size(selector.variable_name); return size; } diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 1809ce91ae..fb452410b6 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -1,7 +1,7 @@ #pragma once #include -#include +#include #include class SchismFormulation : public CoastalFormulation @@ -25,8 +25,8 @@ class SchismFormulation : public CoastalFormulation long record_duration() override; size_t get_ts_index_for_time(const time_t &epoch_time) override; - data_type get_value(const selection_type& selector, ReSampleMethod m) override; - std::vector get_values(const selection_type& selector, ReSampleMethod m) override; + data_type get_value(const selection_type& selector, data_access::ReSampleMethod m) override; + std::vector get_values(const selection_type& selector, data_access::ReSampleMethod m) override; // Implementation of CoastalFormulation void initialize() override; @@ -36,7 +36,7 @@ class SchismFormulation : public CoastalFormulation void get_values(const selection_type& selector, boost::span data) override; private: - std::unique_ptr bmi_; + std::unique_ptr bmi_; // TODO: Some of these maybe should be members of // CoastalFormulation diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 79fd3897cc..5faa54f02e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -134,6 +134,21 @@ endfunction() ########################## Individual tests ########################## +########################## Coastal Formulation Tests +ngen_add_test( + test_coastal + OBJECTS + coastal/SchismFormulation_Test.cpp + LIBRARIES + NGen::core + NGen::realizations_catchment + NGen::core_mediator + NGen::forcing + NGen::ngen_bmi + REQUIRES + NGEN_WITH_BMI_FORTRAN +) + ########################## GeoJSON Unit Tests ngen_add_test( test_geojson diff --git a/test/coastal/SchismFormulation_Test.cpp b/test/coastal/SchismFormulation_Test.cpp new file mode 100644 index 0000000000..c8e43f9dec --- /dev/null +++ b/test/coastal/SchismFormulation_Test.cpp @@ -0,0 +1,9 @@ +#include "gtest/gtest.h" +#include "realizations/coastal/SchismFormulation.hpp" + +struct Schism_Formulation_IT : public ::testing::Test +{ + +}; + + From 4d43b38e0c3fdb2c493a014513f0432f0db8f3f1 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 24 Jun 2024 08:06:26 -0700 Subject: [PATCH 03/51] Another compilation fix --- include/realizations/coastal/CoastalFormulation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/realizations/coastal/CoastalFormulation.hpp b/include/realizations/coastal/CoastalFormulation.hpp index 260ab15140..37b3fd04f4 100644 --- a/include/realizations/coastal/CoastalFormulation.hpp +++ b/include/realizations/coastal/CoastalFormulation.hpp @@ -42,7 +42,7 @@ class CoastalFormulation : public data_access::DataProvider>(selector.points); + auto* points = boost::get>(&selector.points); size_t size = points ? points->size() : this->mesh_size(selector.variable_name); return size; } From af3166f97d3041c94acdf5d8f205b886e2b48f22 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 25 Jun 2024 17:57:12 -0700 Subject: [PATCH 04/51] More work toward Schism Formulation implementation --- CMakeLists.txt | 1 + .../coastal/SchismFormulation.hpp | 5 +++- src/realizations/coastal/CMakeLists.txt | 22 +++++++++++++++ .../coastal/SchismFormulation.cpp | 27 +++++++++++++++++++ 4 files changed, 54 insertions(+), 1 deletion(-) create mode 100644 src/realizations/coastal/CMakeLists.txt create mode 100644 src/realizations/coastal/SchismFormulation.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 22a3af6b41..9d3c473497 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -316,6 +316,7 @@ add_subdirectory("src/core") add_subdirectory("src/geojson") add_subdirectory("src/bmi") add_subdirectory("src/realizations/catchment") +add_subdirectory("src/realizations/coastal") add_subdirectory("src/forcing") add_subdirectory("src/utilities/mdarray") add_subdirectory("src/utilities/mdframe") diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index fb452410b6..29f62c3bb7 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -1,5 +1,7 @@ #pragma once +#include + #include #include #include @@ -9,7 +11,8 @@ class SchismFormulation : public CoastalFormulation public: using MeshPointsDataProvider = data_access::DataProvider; SchismFormulation( - std::string id + std::string const& id + , std::string const& init_config_path , std::shared_ptr met_forcings , std::shared_ptr offshore_boundary , std::shared_ptr inflow_boundary diff --git a/src/realizations/coastal/CMakeLists.txt b/src/realizations/coastal/CMakeLists.txt new file mode 100644 index 0000000000..e60ba8baee --- /dev/null +++ b/src/realizations/coastal/CMakeLists.txt @@ -0,0 +1,22 @@ +include(${PROJECT_SOURCE_DIR}/cmake/dynamic_sourced_library.cmake) +dynamic_sourced_cxx_library(realizations_coastal "${CMAKE_CURRENT_SOURCE_DIR}") + +add_library(NGen::realizations_coastal ALIAS realizations_coastal) + +target_include_directories(realizations_coastal PUBLIC + ${PROJECT_SOURCE_DIR}/include/core + ${PROJECT_SOURCE_DIR}/include/realizations/coastal + ${PROJECT_SOURCE_DIR}/include/forcing + ${PROJECT_SOURCE_DIR}/include/simulation_time + ${PROJECT_SOURCE_DIR}/include/utilities + ${PROJECT_SOURCE_DIR}/include/bmi + ) + +target_link_libraries(realizations_coastal PUBLIC + ${CMAKE_DL_LIBS} + NGen::config_header + NGen::geojson + NGen::logging + NGen::ngen_bmi + ) + diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp new file mode 100644 index 0000000000..2f68715770 --- /dev/null +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -0,0 +1,27 @@ +#include + +const static auto library_path = "/path/to/built/libpschism.so"; +const static auto schism_registration_function = "schism_registration_function"; + +SchismFormulation::SchismFormulation( + std::string const& id + , std::string const& init_config_path + , std::shared_ptr met_forcings + , std::shared_ptr offshore_boundary + , std::shared_ptr inflow_boundary + ) + : CoastalFormulation(id) + , meteorological_forcings_provider_(met_forcings) + , offshore_boundary_provider_(offshore_boundary) + , inflows_boundary_provider_(inflow_boundary) +{ + bmi_ = std::make_unique + ( + "schism_formulation" + , library_path + , init_config_path + , true // allow_model_exceed_end_time + , true // model_time_step_fixed + , schism_registration_function + ); +} From d07ebf86f4593b420d6ae6c79878fcaa1ae39f67 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 26 Jun 2024 10:44:00 -0700 Subject: [PATCH 05/51] Reword --- src/realizations/coastal/SchismFormulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 2f68715770..07f7ba1432 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -17,7 +17,7 @@ SchismFormulation::SchismFormulation( { bmi_ = std::make_unique ( - "schism_formulation" + "schism_coastal_formulation" , library_path , init_config_path , true // allow_model_exceed_end_time From 17a8bfd685d5b1290d1cdcf4497d31d8493f8d16 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 26 Jun 2024 12:07:48 -0700 Subject: [PATCH 06/51] Adjust for refactor --- src/realizations/coastal/SchismFormulation.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 07f7ba1432..f936833b4c 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -20,7 +20,6 @@ SchismFormulation::SchismFormulation( "schism_coastal_formulation" , library_path , init_config_path - , true // allow_model_exceed_end_time , true // model_time_step_fixed , schism_registration_function ); From 28809e0bdcbc87e51dd7a923141689f2cb1b717b Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 1 Jul 2024 15:32:33 -0700 Subject: [PATCH 07/51] Typo fix --- include/realizations/coastal/SchismFormulation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 29f62c3bb7..39168adfc8 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -45,7 +45,7 @@ class SchismFormulation : public CoastalFormulation // CoastalFormulation // TODO: We really want something that we can ask for - // area-averaged RAINRATE over elements, but we're going to makde + // area-averaged RAINRATE over elements, but we're going to make // do with point values at the element centroids std::shared_ptr> meteorological_forcings_provider_; std::shared_ptr> offshore_boundary_provider_; From f345e6d5356b4dc5883992f15e834661cc7ab39b Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 8 Jul 2024 15:11:58 -0700 Subject: [PATCH 08/51] Parameterize SCHISM library path --- src/realizations/coastal/SchismFormulation.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index f936833b4c..3cf01acbd9 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -1,10 +1,10 @@ #include -const static auto library_path = "/path/to/built/libpschism.so"; -const static auto schism_registration_function = "schism_registration_function"; +const static auto s_schism_registration_function = "schism_registration_function"; SchismFormulation::SchismFormulation( std::string const& id + , std::string const& library_path , std::string const& init_config_path , std::shared_ptr met_forcings , std::shared_ptr offshore_boundary @@ -21,6 +21,6 @@ SchismFormulation::SchismFormulation( , library_path , init_config_path , true // model_time_step_fixed - , schism_registration_function + , s_schism_registration_function ); } From 25209b0beb1c1782d0b29c101d3bec8a96eeb9d1 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 7 Aug 2024 17:12:36 -0700 Subject: [PATCH 09/51] Update to match const-qualification changes --- include/realizations/coastal/SchismFormulation.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 39168adfc8..e832dfa4be 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -21,12 +21,12 @@ class SchismFormulation : public CoastalFormulation ~SchismFormulation(); // Implementation of DataProvider - boost::span get_available_variable_names() override; + boost::span get_available_variable_names() const override; - long get_data_start_time() override; - long get_data_stop_time() override; - long record_duration() override; - size_t get_ts_index_for_time(const time_t &epoch_time) override; + long get_data_start_time() const override; + long get_data_stop_time() const override; + long record_duration() const override; + size_t get_ts_index_for_time(const time_t &epoch_time) const override; data_type get_value(const selection_type& selector, data_access::ReSampleMethod m) override; std::vector get_values(const selection_type& selector, data_access::ReSampleMethod m) override; From e003e7602b1dec6478725dea6f3979d992c9c9c1 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 7 Aug 2024 17:12:52 -0700 Subject: [PATCH 10/51] Match added argument --- include/realizations/coastal/SchismFormulation.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index e832dfa4be..4f4c2f99a2 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -12,6 +12,7 @@ class SchismFormulation : public CoastalFormulation using MeshPointsDataProvider = data_access::DataProvider; SchismFormulation( std::string const& id + , std::string const& library_path , std::string const& init_config_path , std::shared_ptr met_forcings , std::shared_ptr offshore_boundary From 7478cd2ad8510dc81f51467df637292f8d79ae71 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 7 Aug 2024 17:13:04 -0700 Subject: [PATCH 11/51] Start outlining main implementation --- src/realizations/coastal/SchismFormulation.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 3cf01acbd9..204f030c8f 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -24,3 +24,18 @@ SchismFormulation::SchismFormulation( , s_schism_registration_function ); } + +void SchismFormulation::initialize() +{ + bmi_->Initialize(); +} + +void SchismFormulation::finalize() +{ + bmi_->Finalize(); +} + +void SchismFormulation::update() +{ + bmi_->Update(); +} From 68dc991d89d25c834805f885dfc2f38163028910 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 7 Aug 2024 17:27:05 -0700 Subject: [PATCH 12/51] Finalize forcings/boundary providers --- src/realizations/coastal/SchismFormulation.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 204f030c8f..3962d8ee28 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -32,6 +32,10 @@ void SchismFormulation::initialize() void SchismFormulation::finalize() { + meteorological_forcings_provider_->finalize(); + offshore_boundary_provider_->finalize(); + inflows_boundary_provider_->finalize(); + bmi_->Finalize(); } From 6fa3c1830e73042c60594479bf95777a13f2734d Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 7 Aug 2024 17:27:23 -0700 Subject: [PATCH 13/51] List variables to set up before Update --- src/realizations/coastal/SchismFormulation.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 3962d8ee28..730762f484 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -41,5 +41,16 @@ void SchismFormulation::finalize() void SchismFormulation::update() { + // Meteorological forcings + // RAINRATE - precipitation + // SFCPRS - surface atmospheric pressure + // SPFH2m - specific humidity at 2m + // TMP2m - temperature at 2m + // UU10m, VV10m - wind velocity components at 10m + + // ETA2_bnd - water surface elevation at the boundaries + // Q_bnd - flows at boundaries + + bmi_->Update(); } From 04c9cf51e4630a6c539038837cd708faeed5afc1 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 7 Aug 2024 17:29:34 -0700 Subject: [PATCH 14/51] Stylize argument-naming comment --- src/realizations/coastal/SchismFormulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 730762f484..5707d49131 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -20,7 +20,7 @@ SchismFormulation::SchismFormulation( "schism_coastal_formulation" , library_path , init_config_path - , true // model_time_step_fixed + , /* model_time_step_fixed = */ true , s_schism_registration_function ); } From 4d55b6ce4fb4e9fb016a1ec737e3e30d66107ffa Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Thu, 8 Aug 2024 14:41:36 -0700 Subject: [PATCH 15/51] Tweak definition of all_points to clarify call site usage --- include/forcing/MeshPointsSelectors.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/forcing/MeshPointsSelectors.hpp b/include/forcing/MeshPointsSelectors.hpp index b27fb1ecd7..bddae4c267 100644 --- a/include/forcing/MeshPointsSelectors.hpp +++ b/include/forcing/MeshPointsSelectors.hpp @@ -4,7 +4,7 @@ #include #include -struct all_points {}; +struct AllPoints {} all_points; struct MeshPointsSelector { @@ -12,5 +12,5 @@ struct MeshPointsSelector std::chrono::time_point init_time; std::chrono::duration duration; std::string output_units; - boost::variant> points; + boost::variant> points; }; From 67e17c8858a17e11656be1f5cb1dd2ef533b0fa3 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 16 Aug 2024 13:28:38 -0700 Subject: [PATCH 16/51] Conditionalize on Fortran support --- include/realizations/coastal/SchismFormulation.hpp | 4 ++++ src/realizations/coastal/SchismFormulation.cpp | 6 ++++++ 2 files changed, 10 insertions(+) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 4f4c2f99a2..295240d51c 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -2,6 +2,8 @@ #include +#if NGEN_WITH_BMI_FORTRAN + #include #include #include @@ -52,3 +54,5 @@ class SchismFormulation : public CoastalFormulation std::shared_ptr> offshore_boundary_provider_; std::shared_ptr> inflows_boundary_provider_; }; + +#endif // NGEN_WITH_BMI_FORTRAN diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 5707d49131..ff749aee79 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -1,3 +1,7 @@ +#include + +#if NGEN_WITH_BMI_FORTRAN + #include const static auto s_schism_registration_function = "schism_registration_function"; @@ -54,3 +58,5 @@ void SchismFormulation::update() bmi_->Update(); } + +#endif // NGEN_WITH_BMI_FORTRAN From c3660148cd52a1e86e04dc7b840d8b220f9c361f Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 16 Aug 2024 13:49:21 -0700 Subject: [PATCH 17/51] Start handling input variable setup --- .../coastal/SchismFormulation.hpp | 6 ++++ .../coastal/SchismFormulation.cpp | 31 +++++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 295240d51c..785dea5f6f 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -7,6 +7,7 @@ #include #include #include +#include class SchismFormulation : public CoastalFormulation { @@ -44,6 +45,11 @@ class SchismFormulation : public CoastalFormulation private: std::unique_ptr bmi_; + static std::set expected_input_variable_names_; + std::map input_variable_units_; + std::map input_variable_type_; + std::map input_variable_count_; + // TODO: Some of these maybe should be members of // CoastalFormulation diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index ff749aee79..ef5b88cb23 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -6,6 +6,18 @@ const static auto s_schism_registration_function = "schism_registration_function"; +std::set SchismFormulation::expected_input_variable_names_ = + { + "RAINRATE", + "SFCPRS", + "SPHF2m", + "TMP2m", + "UU10m", + "VV10m", + "ETA2_bnd", + "Q_bnd" + }; + SchismFormulation::SchismFormulation( std::string const& id , std::string const& library_path @@ -27,6 +39,25 @@ SchismFormulation::SchismFormulation( , /* model_time_step_fixed = */ true , s_schism_registration_function ); + + auto const& input_vars = bmi_->GetInputVarNames(); + + for (auto const& name : input_vars) { + if (expected_input_variable_names_.find(name) == expected_input_variable_names_.end()) { + throw std::runtime_error("SCHISM instance requests unexpected input variable '" + name + "'"); + } + + input_variable_units_[name] = bmi_->GetVarUnits(name); + input_variable_type_[name] = bmi_->GetVarType(name); + + auto nbytes = bmi_->GetVarNbytes(name); + auto itemsize = bmi_->GetVarItemsize(name); + if (nbytes % itemsize != 0) { + throw std::runtime_error("For SCHISM input variable '" + name + "', itemsize " + std::to_string(itemsize) + + " does not evenly divide nbytes " + std::to_string(nbytes)); + } + input_variable_count_[name] = nbytes / itemsize; + } } void SchismFormulation::initialize() From 706046b1a40582db521db67abb01b574f3317053 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 16 Aug 2024 13:49:44 -0700 Subject: [PATCH 18/51] Sketch time keeping members --- include/realizations/coastal/SchismFormulation.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 785dea5f6f..65fc5e3439 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -50,6 +50,9 @@ class SchismFormulation : public CoastalFormulation std::map input_variable_type_; std::map input_variable_count_; + time_t current_time_; + std::chrono::seconds time_step_length_; + // TODO: Some of these maybe should be members of // CoastalFormulation From 175112508c7e0884a843666ff611dfa1a95a2efa Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 16 Aug 2024 14:01:04 -0700 Subject: [PATCH 19/51] Call BMI Init function --- src/realizations/coastal/SchismFormulation.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index ef5b88cb23..12bb6e3eb4 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -40,6 +40,8 @@ SchismFormulation::SchismFormulation( , s_schism_registration_function ); + bmi_->Initialize(init_config_path); + auto const& input_vars = bmi_->GetInputVarNames(); for (auto const& name : input_vars) { From 2dba5650abe417791af8efa2a93cdc0bbef816dc Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 16 Aug 2024 14:01:58 -0700 Subject: [PATCH 20/51] Move input variable setup from ctor to Initialize --- src/realizations/coastal/SchismFormulation.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 12bb6e3eb4..f9e162c130 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -39,8 +39,11 @@ SchismFormulation::SchismFormulation( , /* model_time_step_fixed = */ true , s_schism_registration_function ); +} - bmi_->Initialize(init_config_path); +void SchismFormulation::initialize() +{ + bmi_->Initialize(); auto const& input_vars = bmi_->GetInputVarNames(); @@ -62,11 +65,6 @@ SchismFormulation::SchismFormulation( } } -void SchismFormulation::initialize() -{ - bmi_->Initialize(); -} - void SchismFormulation::finalize() { meteorological_forcings_provider_->finalize(); From d687db03314237038f62ddee20da3da5426d8549 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 19 Aug 2024 12:14:55 -0700 Subject: [PATCH 21/51] Match SCHISM BMI's registration/instantiaion function name --- src/realizations/coastal/SchismFormulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index f9e162c130..a61c77f01b 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -4,7 +4,7 @@ #include -const static auto s_schism_registration_function = "schism_registration_function"; +const static auto s_schism_registration_function = "register_bmi"; std::set SchismFormulation::expected_input_variable_names_ = { From b7d1ebada630aa1ed4e2929a1d19f1c52ecec0da Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 23 Aug 2024 09:16:09 -0700 Subject: [PATCH 22/51] Fix declaration of all_points to avoid linking errors on duplicate definition --- include/forcing/MeshPointsSelectors.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/forcing/MeshPointsSelectors.hpp b/include/forcing/MeshPointsSelectors.hpp index bddae4c267..0a230e34c3 100644 --- a/include/forcing/MeshPointsSelectors.hpp +++ b/include/forcing/MeshPointsSelectors.hpp @@ -4,7 +4,8 @@ #include #include -struct AllPoints {} all_points; +struct AllPoints {}; +static AllPoints all_points; struct MeshPointsSelector { From 9b5716e54bd752d642402a022ca5b0e2c58e7243 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 23 Aug 2024 09:17:15 -0700 Subject: [PATCH 23/51] Conditionalize SCHISM code on Fortran and MPI support --- include/realizations/coastal/SchismFormulation.hpp | 4 ++-- src/realizations/coastal/SchismFormulation.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 65fc5e3439..83a8cf395c 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -2,7 +2,7 @@ #include -#if NGEN_WITH_BMI_FORTRAN +#if NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI #include #include @@ -64,4 +64,4 @@ class SchismFormulation : public CoastalFormulation std::shared_ptr> inflows_boundary_provider_; }; -#endif // NGEN_WITH_BMI_FORTRAN +#endif // NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index a61c77f01b..a10e7210f2 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -1,6 +1,6 @@ #include -#if NGEN_WITH_BMI_FORTRAN +#if NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI #include @@ -90,4 +90,4 @@ void SchismFormulation::update() bmi_->Update(); } -#endif // NGEN_WITH_BMI_FORTRAN +#endif // NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI From fa6c223ab379ff507755b94af14ea0383e5a2730 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 23 Aug 2024 09:17:55 -0700 Subject: [PATCH 24/51] Square away virtual method implementation --- .../coastal/SchismFormulation.hpp | 6 ++- .../coastal/SchismFormulation.cpp | 52 +++++++++++++++++++ 2 files changed, 56 insertions(+), 2 deletions(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 83a8cf395c..b914a23298 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -9,7 +9,7 @@ #include #include -class SchismFormulation : public CoastalFormulation +class SchismFormulation final : public CoastalFormulation { public: using MeshPointsDataProvider = data_access::DataProvider; @@ -33,7 +33,6 @@ class SchismFormulation : public CoastalFormulation size_t get_ts_index_for_time(const time_t &epoch_time) const override; data_type get_value(const selection_type& selector, data_access::ReSampleMethod m) override; - std::vector get_values(const selection_type& selector, data_access::ReSampleMethod m) override; // Implementation of CoastalFormulation void initialize() override; @@ -42,6 +41,9 @@ class SchismFormulation : public CoastalFormulation void get_values(const selection_type& selector, boost::span data) override; +protected: + size_t mesh_size(std::string const& variable_name) override; + private: std::unique_ptr bmi_; diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index a10e7210f2..3e01324ee8 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -41,6 +41,8 @@ SchismFormulation::SchismFormulation( ); } +SchismFormulation::~SchismFormulation() = default; + void SchismFormulation::initialize() { bmi_->Initialize(); @@ -83,6 +85,8 @@ void SchismFormulation::update() // TMP2m - temperature at 2m // UU10m, VV10m - wind velocity components at 10m + //auto rain_points = MeshPointsSelector{"RAINRATE", current_time_, time_step_length_, input_variable_units_["RAINRATE"], all_points}; + // ETA2_bnd - water surface elevation at the boundaries // Q_bnd - flows at boundaries @@ -90,4 +94,52 @@ void SchismFormulation::update() bmi_->Update(); } +boost::span SchismFormulation::get_available_variable_names() const +{ + throw std::runtime_error(__func__); + return {}; +} + +long SchismFormulation::get_data_start_time() const +{ + throw std::runtime_error(__func__); + return 0; +} + +long SchismFormulation::get_data_stop_time() const +{ + throw std::runtime_error(__func__); + return 0; +} + +long SchismFormulation::record_duration() const +{ + throw std::runtime_error(__func__); + return 0; +} + +size_t SchismFormulation::get_ts_index_for_time(const time_t &epoch_time) const +{ + throw std::runtime_error(__func__); + return 0; +} + +SchismFormulation::data_type SchismFormulation::get_value(const selection_type& selector, data_access::ReSampleMethod m) +{ + throw std::runtime_error(__func__); + return 0.0; +} + +void SchismFormulation::get_values(const selection_type& selector, boost::span data) +{ + throw std::runtime_error(__func__); +} + +size_t SchismFormulation::mesh_size(std::string const& variable_name) +{ + throw std::runtime_error(__func__); + return 0; +} + + #endif // NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI From 56783a34ee0c4aca5e3e4c9d8059393ddc412955 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 23 Aug 2024 14:43:59 -0700 Subject: [PATCH 25/51] Typo fix in variable name --- src/realizations/coastal/SchismFormulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 3e01324ee8..991fed057e 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -10,7 +10,7 @@ std::set SchismFormulation::expected_input_variable_names_ = { "RAINRATE", "SFCPRS", - "SPHF2m", + "SPFH2m", "TMP2m", "UU10m", "VV10m", From 23725d169f12e550610a4de8bd093257948c06f2 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 23 Aug 2024 14:45:30 -0700 Subject: [PATCH 26/51] Bmi_Fortran_Adapter: Add MPI-specific constructor --- include/bmi/Bmi_Fortran_Adapter.hpp | 44 +++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/include/bmi/Bmi_Fortran_Adapter.hpp b/include/bmi/Bmi_Fortran_Adapter.hpp index 9f5bfe71bb..8abedef1b4 100644 --- a/include/bmi/Bmi_Fortran_Adapter.hpp +++ b/include/bmi/Bmi_Fortran_Adapter.hpp @@ -11,6 +11,10 @@ #include "State_Exception.hpp" #include "utilities/ExternalIntegrationException.hpp" +#if NGEN_WITH_MPI +#include +#endif // NGEN_WITH_MPI + // Forward declaration to provide access to protected items in testing class Bmi_Fortran_Adapter_Test; @@ -76,6 +80,46 @@ namespace models { } } +#if NGEN_WITH_MPI + // Special case constructor for formulation modules that + // implement the NGen BMI MPI protocol, in which the MPI + // Communicator is set through SetValue() between + // construction and Initialize() + Bmi_Fortran_Adapter(const std::string &type_name, + std::string library_file_path, + std::string bmi_init_config, + bool has_fixed_time_step, + std::string registration_func, + MPI_Comm comm) + : AbstractCLibBmiAdapter(type_name, + library_file_path, + bmi_init_config, + has_fixed_time_step, + registration_func + ) + { + try { + bmi_model = std::make_unique(Bmi_Fortran_Handle_Wrapper()); + dynamic_library_load(); + execModuleRegistration(); + + MPI_Fint comm_fortran = MPI_Comm_c2f(comm); + inner_set_value_int("bmi_mpi_comm_handle", &comm_fortran); + + int init_result = initialize(&bmi_model->handle, bmi_init_config.c_str()); + if (init_result != BMI_SUCCESS) { + init_exception_msg = "Failure when attempting to initialize " + model_name; + throw models::external::State_Exception(init_exception_msg); + } + model_initialized = true; + } + catch (...) { + model_initialized = true; + throw; + } + } +#endif // NGEN_WITH_MPI + Bmi_Fortran_Adapter(Bmi_Fortran_Adapter const&) = delete; Bmi_Fortran_Adapter(Bmi_Fortran_Adapter&&) = delete; From d618d800bf6edea9b196cfc0f8aedf194000f42e Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Fri, 23 Aug 2024 14:46:25 -0700 Subject: [PATCH 27/51] SchismFormulation: Match usage to Bmi_Fortran_Adapter constructor taking MPI_Comm (which calls BMI initialize) --- src/realizations/coastal/SchismFormulation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 991fed057e..02e591789f 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -3,6 +3,7 @@ #if NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI #include +#include const static auto s_schism_registration_function = "register_bmi"; @@ -38,6 +39,7 @@ SchismFormulation::SchismFormulation( , init_config_path , /* model_time_step_fixed = */ true , s_schism_registration_function + , MPI_COMM_WORLD ); } @@ -45,8 +47,6 @@ SchismFormulation::~SchismFormulation() = default; void SchismFormulation::initialize() { - bmi_->Initialize(); - auto const& input_vars = bmi_->GetInputVarNames(); for (auto const& name : input_vars) { From c7eaf3ec8b6e084760867580739b8e2f8912af07 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 11:42:50 -0700 Subject: [PATCH 28/51] Require MPI for the test --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5faa54f02e..26e83a963d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -147,6 +147,7 @@ ngen_add_test( NGen::ngen_bmi REQUIRES NGEN_WITH_BMI_FORTRAN + NGEN_WITH_MPI ) ########################## GeoJSON Unit Tests From 8b4d612743e68f5bb6dc99e3c19b1c1e9a4c0d19 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 11:51:29 -0700 Subject: [PATCH 29/51] Enumerate output variables --- include/realizations/coastal/SchismFormulation.hpp | 3 +++ src/realizations/coastal/SchismFormulation.cpp | 11 +++++++++-- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index b914a23298..588535e563 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -8,6 +8,7 @@ #include #include #include +#include class SchismFormulation final : public CoastalFormulation { @@ -52,6 +53,8 @@ class SchismFormulation final : public CoastalFormulation std::map input_variable_type_; std::map input_variable_count_; + static std::vector exported_output_variable_names_; + time_t current_time_; std::chrono::seconds time_step_length_; diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 02e591789f..495c58ec35 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -19,6 +19,14 @@ std::set SchismFormulation::expected_input_variable_names_ = "Q_bnd" }; +std::vector SchismFormulation::exported_output_variable_names_ = + { + "ETA2", + "VX", + "VY", + "BEDLEVEL" + }; + SchismFormulation::SchismFormulation( std::string const& id , std::string const& library_path @@ -96,8 +104,7 @@ void SchismFormulation::update() boost::span SchismFormulation::get_available_variable_names() const { - throw std::runtime_error(__func__); - return {}; + return exported_output_variable_names_; } long SchismFormulation::get_data_start_time() const From 90b33d98355e4baaa96d89bf5c1f9bba83c934e8 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 12:11:15 -0700 Subject: [PATCH 30/51] Implement mesh_size method --- src/realizations/coastal/SchismFormulation.cpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 495c58ec35..fbe906762d 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -64,14 +64,7 @@ void SchismFormulation::initialize() input_variable_units_[name] = bmi_->GetVarUnits(name); input_variable_type_[name] = bmi_->GetVarType(name); - - auto nbytes = bmi_->GetVarNbytes(name); - auto itemsize = bmi_->GetVarItemsize(name); - if (nbytes % itemsize != 0) { - throw std::runtime_error("For SCHISM input variable '" + name + "', itemsize " + std::to_string(itemsize) + - " does not evenly divide nbytes " + std::to_string(nbytes)); - } - input_variable_count_[name] = nbytes / itemsize; + input_variable_count_[name] = mesh_size(name); } } @@ -144,8 +137,13 @@ void SchismFormulation::get_values(const selection_type& selector, boost::spanGetVarNbytes(variable_name); + auto itemsize = bmi_->GetVarItemsize(variable_name); + if (nbytes % itemsize != 0) { + throw std::runtime_error("For SCHISM input variable '" + variable_name + "', itemsize " + std::to_string(itemsize) + + " does not evenly divide nbytes " + std::to_string(nbytes)); + } + return nbytes / itemsize; } From 9d74eb2cf458b4f6c616a60312e64d2a82bc6f35 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 12:16:23 -0700 Subject: [PATCH 31/51] Implement get_values --- src/realizations/coastal/SchismFormulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index fbe906762d..969d936108 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -132,7 +132,7 @@ SchismFormulation::data_type SchismFormulation::get_value(const selection_type& void SchismFormulation::get_values(const selection_type& selector, boost::span data) { - throw std::runtime_error(__func__); + bmi_->GetValue(selector.variable_name, data.data()); } size_t SchismFormulation::mesh_size(std::string const& variable_name) From b9c12374fde4f8f80b88753cb148fe348a10b51f Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 12:17:01 -0700 Subject: [PATCH 32/51] Generalize output string --- src/realizations/coastal/SchismFormulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 969d936108..9c6c2af04b 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -140,7 +140,7 @@ size_t SchismFormulation::mesh_size(std::string const& variable_name) auto nbytes = bmi_->GetVarNbytes(variable_name); auto itemsize = bmi_->GetVarItemsize(variable_name); if (nbytes % itemsize != 0) { - throw std::runtime_error("For SCHISM input variable '" + variable_name + "', itemsize " + std::to_string(itemsize) + + throw std::runtime_error("For SCHISM variable '" + variable_name + "', itemsize " + std::to_string(itemsize) + " does not evenly divide nbytes " + std::to_string(nbytes)); } return nbytes / itemsize; From ffb7057a7e5297aa85bcb69fe347b18bc19789ac Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 12:21:51 -0700 Subject: [PATCH 33/51] Position comments better --- src/realizations/coastal/SchismFormulation.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 9c6c2af04b..6210f35e3b 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -9,13 +9,23 @@ const static auto s_schism_registration_function = "register_bmi"; std::set SchismFormulation::expected_input_variable_names_ = { + /* Meteorological Forcings */ + // RAINRATE - precipitation "RAINRATE", + // SFCPRS - surface atmospheric pressure "SFCPRS", + // SPFH2m - specific humidity at 2m "SPFH2m", + // TMP2m - temperature at 2m "TMP2m", + // UU10m, VV10m - wind velocity components at 10m "UU10m", "VV10m", + + /* Input Boundary Conditions */ + // ETA2_bnd - water surface elevation at the boundaries "ETA2_bnd", + // Q_bnd - flows at boundaries "Q_bnd" }; @@ -79,7 +89,6 @@ void SchismFormulation::finalize() void SchismFormulation::update() { - // Meteorological forcings // RAINRATE - precipitation // SFCPRS - surface atmospheric pressure // SPFH2m - specific humidity at 2m From 5835aacd3ad07d666fcc2a7708997cfd9cf7dcb6 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 12:59:04 -0700 Subject: [PATCH 34/51] Square up time types --- include/forcing/MeshPointsSelectors.hpp | 2 +- include/realizations/coastal/SchismFormulation.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/include/forcing/MeshPointsSelectors.hpp b/include/forcing/MeshPointsSelectors.hpp index 0a230e34c3..680f882388 100644 --- a/include/forcing/MeshPointsSelectors.hpp +++ b/include/forcing/MeshPointsSelectors.hpp @@ -11,7 +11,7 @@ struct MeshPointsSelector { std::string variable_name; std::chrono::time_point init_time; - std::chrono::duration duration; + std::chrono::seconds duration; std::string output_units; boost::variant> points; }; diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 588535e563..c969f9365d 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -55,7 +55,7 @@ class SchismFormulation final : public CoastalFormulation static std::vector exported_output_variable_names_; - time_t current_time_; + std::chrono::time_point current_time_; std::chrono::seconds time_step_length_; // TODO: Some of these maybe should be members of From 686a706d6af2a9069eb14721ed3bc0c174abf64a Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 13:01:34 -0700 Subject: [PATCH 35/51] Implement much of forcings --- .../coastal/SchismFormulation.hpp | 15 +++-- .../coastal/SchismFormulation.cpp | 60 ++++++++++++------- 2 files changed, 47 insertions(+), 28 deletions(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index c969f9365d..9ed96e4735 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -7,7 +7,7 @@ #include #include #include -#include +#include #include class SchismFormulation final : public CoastalFormulation @@ -48,7 +48,8 @@ class SchismFormulation final : public CoastalFormulation private: std::unique_ptr bmi_; - static std::set expected_input_variable_names_; + enum ForcingSelector { METEO, OFFSHORE, INFLOW }; + static std::map expected_input_variables_; std::map input_variable_units_; std::map input_variable_type_; std::map input_variable_count_; @@ -64,9 +65,13 @@ class SchismFormulation final : public CoastalFormulation // TODO: We really want something that we can ask for // area-averaged RAINRATE over elements, but we're going to make // do with point values at the element centroids - std::shared_ptr> meteorological_forcings_provider_; - std::shared_ptr> offshore_boundary_provider_; - std::shared_ptr> inflows_boundary_provider_; + + using ProviderType = data_access::DataProvider; + std::shared_ptr meteorological_forcings_provider_; + std::shared_ptr offshore_boundary_provider_; + std::shared_ptr inflows_boundary_provider_; + + void set_inputs(); }; #endif // NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 6210f35e3b..abf4397c53 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -7,26 +7,28 @@ const static auto s_schism_registration_function = "register_bmi"; -std::set SchismFormulation::expected_input_variable_names_ = + + +std::map SchismFormulation::expected_input_variables_ = { /* Meteorological Forcings */ // RAINRATE - precipitation - "RAINRATE", + {"RAINRATE", SchismFormulation::METEO}, // SFCPRS - surface atmospheric pressure - "SFCPRS", + {"SFCPRS", SchismFormulation::METEO}, // SPFH2m - specific humidity at 2m - "SPFH2m", + {"SPFH2m", SchismFormulation::METEO}, // TMP2m - temperature at 2m - "TMP2m", + {"TMP2m", SchismFormulation::METEO}, // UU10m, VV10m - wind velocity components at 10m - "UU10m", - "VV10m", + {"UU10m", SchismFormulation::METEO}, + {"VV10m", SchismFormulation::METEO}, /* Input Boundary Conditions */ // ETA2_bnd - water surface elevation at the boundaries - "ETA2_bnd", + {"ETA2_bnd", SchismFormulation::OFFSHORE}, // Q_bnd - flows at boundaries - "Q_bnd" + {"Q_bnd", SchismFormulation::INFLOW}, }; std::vector SchismFormulation::exported_output_variable_names_ = @@ -68,7 +70,7 @@ void SchismFormulation::initialize() auto const& input_vars = bmi_->GetInputVarNames(); for (auto const& name : input_vars) { - if (expected_input_variable_names_.find(name) == expected_input_variable_names_.end()) { + if (expected_input_variables_.find(name) == expected_input_variables_.end()) { throw std::runtime_error("SCHISM instance requests unexpected input variable '" + name + "'"); } @@ -76,31 +78,43 @@ void SchismFormulation::initialize() input_variable_type_[name] = bmi_->GetVarType(name); input_variable_count_[name] = mesh_size(name); } + + //set_inputs(); } void SchismFormulation::finalize() { +#if 0 meteorological_forcings_provider_->finalize(); offshore_boundary_provider_->finalize(); inflows_boundary_provider_->finalize(); - +#endif bmi_->Finalize(); } -void SchismFormulation::update() +void SchismFormulation::set_inputs() { - // RAINRATE - precipitation - // SFCPRS - surface atmospheric pressure - // SPFH2m - specific humidity at 2m - // TMP2m - temperature at 2m - // UU10m, VV10m - wind velocity components at 10m - - //auto rain_points = MeshPointsSelector{"RAINRATE", current_time_, time_step_length_, input_variable_units_["RAINRATE"], all_points}; - - // ETA2_bnd - water surface elevation at the boundaries - // Q_bnd - flows at boundaries - + for (auto var : expected_input_variables_) { + auto& name = var.first; + auto selector = var.second; + auto points = MeshPointsSelector{name, current_time_, time_step_length_, input_variable_units_[name], all_points}; + + ProviderType* provider = [this, selector](){ + switch(selector) { + case METEO: return meteorological_forcings_provider_.get(); + case OFFSHORE: return offshore_boundary_provider_.get(); + case INFLOW: return inflows_boundary_provider_.get(); + default: throw std::runtime_error("Unknown SCHISM provider selector type"); + } + }(); + std::vector values = provider->get_values(points); + bmi_->SetValue(name, values.data()); + } +} +void SchismFormulation::update() +{ + set_inputs(); bmi_->Update(); } From bbe5c5faa770f1d4b16c4c36f5ee883c02b4b531 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 17:35:36 -0700 Subject: [PATCH 36/51] Split parallel_utils.h to header/source, and make bits usable from multiple calling objects --- include/utilities/parallel_utils.h | 361 ++--------------------------- src/NGen.cpp | 8 +- src/core/parallel_utils.cpp | 352 ++++++++++++++++++++++++++++ 3 files changed, 369 insertions(+), 352 deletions(-) create mode 100644 src/core/parallel_utils.cpp diff --git a/include/utilities/parallel_utils.h b/include/utilities/parallel_utils.h index 3e3a70a85d..509649e0dd 100644 --- a/include/utilities/parallel_utils.h +++ b/include/utilities/parallel_utils.h @@ -2,6 +2,10 @@ #define NGEN_PARALLEL_UTILS_H #include + +// Define in the non-MPI case so that we don't need to conditionally compile `if (mpi_rank == 0)` +extern int mpi_rank; + #if NGEN_WITH_MPI #ifndef MPI_HF_SUB_CODE_GOOD @@ -29,6 +33,8 @@ #include "python/HydrofabricSubsetter.hpp" #endif // NGEN_WITH_PYTHON +extern int mpi_num_procs; + namespace parallel { /** @@ -50,37 +56,7 @@ namespace parallel { * is not in the success/ready state. * @return Whether all ranks coordinating status had a success/ready status value. */ - bool mpiSyncStatusAnd(bool status, int mpi_rank, int mpi_num_procs, const std::string &taskDesc) { - - // Expect 0 is good and 1 is no good for goodCode - // TODO: assert this in constructor or somewhere, or maybe just in a unit test - unsigned short codeBuffer; - bool printMessage = !taskDesc.empty(); - // For the other ranks, start by properly setting the status code value in the buffer and send to rank 0 - if (mpi_rank != 0) { - codeBuffer = status ? MPI_HF_SUB_CODE_GOOD : MPI_HF_SUB_CODE_BAD; - MPI_Send(&codeBuffer, 1, MPI_UNSIGNED_SHORT, 0, 0, MPI_COMM_WORLD); - } - // In rank 0, the first step is to receive and process codes from the other ranks into unified global status - else { - for (int i = 1; i < mpi_num_procs; ++i) { - MPI_Recv(&codeBuffer, 1, MPI_UNSIGNED_SHORT, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - // If any is ever "not good", overwrite status to be "false" - if (codeBuffer != MPI_HF_SUB_CODE_GOOD) { - if (printMessage) { - std::cout << "Rank " << i << " not successful/ready after " << taskDesc << std::endl; - } - status = false; - } - } - // Rank 0 must also now prepare the codeBuffer value for broadcasting the global status - codeBuffer = status ? MPI_HF_SUB_CODE_GOOD : MPI_HF_SUB_CODE_BAD; - } - - // Execute broadcast of global status rooted at rank 0 - MPI_Bcast(&codeBuffer, 1, MPI_UNSIGNED_SHORT, 0, MPI_COMM_WORLD); - return codeBuffer == MPI_HF_SUB_CODE_GOOD; - } + bool mpiSyncStatusAnd(bool status, int mpi_rank, int mpi_num_procs, const std::string &taskDesc); /** * Convenience method for overloaded function when no message is needed, and thus no description param provided. @@ -91,7 +67,7 @@ namespace parallel { * @return Whether all ranks coordinating status had a success/ready status value. * @see mpiSyncStatusAnd(bool, const std::string&) */ - bool mpiSyncStatusAnd(bool status, int mpi_rank, int mpi_num_procs) { + inline bool mpiSyncStatusAnd(bool status, int mpi_rank, int mpi_num_procs) { return mpiSyncStatusAnd(status, mpi_rank, mpi_num_procs, ""); } @@ -113,34 +89,7 @@ namespace parallel { * @return Whether proprocessing has already been performed to divide the main hydrofabric into existing, individual * sub-hydrofabric files for each partition/process. */ - bool is_hydrofabric_subdivided(int mpi_rank, int mpi_num_procs, bool printMsg) { - std::string name = catchmentDataFile + "." + std::to_string(mpi_rank); - // Initialize isGood based on local state. Here, local file is "good" when it already exists. - // TODO: this isn't actually checking whether the files are right (just that they are present) so do we need to? - bool isGood = utils::FileChecker::file_is_readable(name); - - if (mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs)) { - if (printMsg) { std::cout << "Process " << mpi_rank << ": Hydrofabric already subdivided in " << mpi_num_procs << " files." << std::endl; } - return true; - } - else { - if (printMsg) { std::cout << "Process " << mpi_rank << ": Hydrofabric has not yet been subdivided." << std::endl; } - return false; - } - } - - /** - * Convenience overloaded method for when no supplemental output message is required. - * - * @param mpi_rank The rank of the current process. - * @param mpi_num_procs The total number of MPI processes. - * @return Whether proprocessing has already been performed to divide the main hydrofabric into existing, individual - * sub-hydrofabric files for each partition/process. - * @see is_hydrofabric_subdivided(bool) - */ - bool is_hydrofabric_subdivided(int mpi_rank, int mpi_num_procs) { - return is_hydrofabric_subdivided(mpi_rank, mpi_num_procs, false); - } + bool is_hydrofabric_subdivided(int mpi_rank, int mpi_num_procs, std::string const& catchmentDataFile, bool printMsg); /** * Set each rank's host "id" value in a provided host array. @@ -157,78 +106,7 @@ namespace parallel { * @param mpi_num_procs The number of ranks. * @param host_array A pointer to an allocated array of size ``mpi_num_procs``. */ - void get_hosts_array(int mpi_rank, int mpi_num_procs, int *host_array) { - const int ROOT_RANK = 0; - // These are the lengths of the (trimmed) C-string representations of the hostname for each rank - std::vector actualHostnameCStrLength(mpi_num_procs); - // Initialize to -1 to represent unknown - for (int i = 0; i < mpi_num_procs; ++i) { - actualHostnameCStrLength[i] = -1; - } - - // Get this rank's hostname (things should never be longer than 256) - char myhostname[256]; - gethostname(myhostname, 256); - - // Set the one for this rank - actualHostnameCStrLength[mpi_rank] = std::strlen(myhostname); - - // First, gather the hostname string lengths - MPI_Gather(&actualHostnameCStrLength[mpi_rank], 1, MPI_INT, actualHostnameCStrLength.data(), 1, MPI_INT, ROOT_RANK, - MPI_COMM_WORLD); - // Per-rank start offsets/displacements in our hostname strings gather buffer. - std::vector recvDisplacements(mpi_num_procs); - int totalLength = 0; - for (int i = 0; i < mpi_num_procs; ++i) { - // Displace each rank's string by the sum of the length of the previous rank strings - recvDisplacements[i] = totalLength; - // Adding the extra to make space for the \0 - actualHostnameCStrLength[i] += 1; - // Then update the total length - totalLength += actualHostnameCStrLength[i]; - - } - // Now we can create our buffer array and gather the hostname strings into it - char hostnames[totalLength]; - MPI_Gatherv(myhostname, actualHostnameCStrLength[mpi_rank], MPI_CHAR, hostnames, actualHostnameCStrLength.data(), - recvDisplacements.data(), MPI_CHAR, ROOT_RANK, MPI_COMM_WORLD); - - if (mpi_rank == ROOT_RANK) { - host_array[0] = 0; - int next_host_id = 1; - - int rank_with_matching_hostname; - char *checked_rank_hostname, *known_rank_hostname; - - for (int rank_being_check = 0; rank_being_check < mpi_num_procs; ++rank_being_check) { - // Set this as negative initially for each rank check to indicate no match found (at least yet) - rank_with_matching_hostname = -1; - // Get a C-string pointer for this rank's hostname, offset by the appropriate displacement - checked_rank_hostname = &hostnames[recvDisplacements[rank_being_check]]; - - // Assume that hostnames for any ranks less than the current rank being check are already known - for (int known_rank = 0; known_rank < rank_being_check; ++known_rank) { - // Get the right C-string pointer for the current known rank's hostname also - known_rank_hostname = &hostnames[recvDisplacements[known_rank]]; - // Compare the hostnames, setting and breaking if a match is found - if (std::strcmp(known_rank_hostname, checked_rank_hostname) == 0) { - rank_with_matching_hostname = known_rank; - break; - } - } - // This indicates this rank had no earlier rank with a matching hostname. - if (rank_with_matching_hostname < 0) { - // Assign new host id, then increment what the next id will be - host_array[rank_being_check] = next_host_id++; - } - else { - host_array[rank_being_check] = host_array[rank_with_matching_hostname]; - } - } - } - // Now, broadcast the results out - MPI_Bcast(host_array, mpi_num_procs, MPI_INT, 0, MPI_COMM_WORLD); - } + void get_hosts_array(int mpi_rank, int mpi_num_procs, int *host_array); /** * Send the contents of a text file to another MPI rank. @@ -240,58 +118,7 @@ namespace parallel { * @param destRank The MPI rank to which the file data should be sent. * @return Whether sending was successful. */ - bool mpi_send_text_file(const char *fileName, const int mpi_rank, const int destRank) { - int bufSize = 4096; - std::vector buf(bufSize); - int code; - // How much has been transferred so far - int totalNumTransferred = 0; - - FILE *file = fopen(fileName, "r"); - - // Transmit error code instead of expected size and return false if file can't be opened - if (file == NULL) { - // TODO: output error message - code = -1; - MPI_Send(&code, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - return false; - } - - // Send expected size to start - MPI_Send(&bufSize, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - - // Then get back expected size to infer other side is good to go - MPI_Recv(&code, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - if (code != bufSize) { - // TODO: output error message - fclose(file); - return false; - } - int continueCode = 1; - // Then while there is more of the file to read and send, read the next batch and ... - while (fgets(buf.data(), bufSize, file) != NULL) { - // Indicate we are ready to continue sending data - MPI_Send(&continueCode, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - - // Send this batch - MPI_Send(buf.data(), bufSize, MPI_CHAR, destRank, NGEN_MPI_DATA_TAG, MPI_COMM_WORLD); - - // Then get back a code, which will be -1 if bad and need to exit and otherwise good - MPI_Recv(&code, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - if (code < 0) { - // TODO: output error message - fclose(file); - return false; - } - } - // Once there is no more file to read and send, we should stop continuing - continueCode = 0; - MPI_Send(&continueCode, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - // Expect to get back a code of 0 - MPI_Recv(&code, 1, MPI_INT, destRank, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - fclose(file); - return code == 0; - } + bool mpi_send_text_file(const char *fileName, const int mpi_rank, const int destRank); /** * Receive text data from another MPI rank and write the contents to a text file. @@ -305,58 +132,7 @@ namespace parallel { * @param destRank The MPI rank to which the file data should be sent. * @return Whether sending was successful. */ - bool mpi_recv_text_file(const char *fileName, const int mpi_rank, const int srcRank) { - int bufSize, writeCode; - // Receive expected buffer size to start - MPI_Recv(&bufSize, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - - // If the sending side couldn't open the file, then immediately return false - if (bufSize == -1) { - // TODO: output error - return false; - } - - // Try to open recv file ... - FILE *file = fopen(fileName, "w"); - // ... and let sending size know whether this was successful by sending error code if not ... - if (file == NULL) { - // TODO: output error message - bufSize = -1; - MPI_Send(&bufSize, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - return false; - } - - // Send back the received buffer it if file opened, confirming things are good to go for transfer - MPI_Send(&bufSize, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - - // How much has been transferred so far - int totalNumTransferred = 0; - std::vector buf(bufSize); - - int continueCode; - - while (true) { - // Make sure the other side wants to continue sending data - MPI_Recv(&continueCode, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - if (continueCode <= 0) - break; - - MPI_Recv(buf.data(), bufSize, MPI_CHAR, srcRank, NGEN_MPI_DATA_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); - - writeCode = fputs(buf.data(), file); - MPI_Send(&writeCode, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - - if (writeCode < 0) { - fclose(file); - return false; - } - } - - fclose(file); - MPI_Send(&continueCode, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); - return true; - } - + bool mpi_recv_text_file(const char *fileName, const int mpi_rank, const int srcRank); /** * Distribute subdivided hydrofabric files to ranks on other hosts as needed. @@ -409,51 +185,7 @@ namespace parallel { */ bool distribute_subdivided_hydrofabric_files(const std::string &baseCatchmentFile, const std::string &baseNexusFile, const int sendingRank, const int mpi_rank, const int mpi_num_procs, - const int *hostIdForRank, bool syncReturnStatus, bool blockAll) - { - // Start with status as good - bool isGood = true; - // FIXME: For now, just have rank 0 send everything, but optimize with multiple procs or threads later - // Only need to process this if sending rank or a receiving ranks (i.e., not on same host as sending rank) - if (mpi_rank == sendingRank || hostIdForRank[mpi_rank] != hostIdForRank[sendingRank]) { - // Have the sending rank send out all files - if (mpi_rank == sendingRank) { - // In rank 0, for all the other ranks ... - for (int otherRank = 0; otherRank < mpi_num_procs; ++otherRank) { - // If another rank is on a different host (note that this covers otherRank == sendingRank case) ... - if (hostIdForRank[otherRank] != hostIdForRank[mpi_rank]) { - // ... then send that rank its rank-specific catchment and nexus files - std::string catFileToSend = baseCatchmentFile + "." + std::to_string(otherRank); - std::string nexFileToSend = baseNexusFile + "." + std::to_string(otherRank); - // Note that checking previous isGood is necessary here because of loop - isGood = isGood && mpi_send_text_file(catFileToSend.c_str(), mpi_rank, otherRank); - isGood = isGood && mpi_send_text_file(nexFileToSend.c_str(), mpi_rank, otherRank); - } - } - } - else { - // For a rank not on the same host as the sending rank, receive the transmitted file - std::string catFileToReceive = baseCatchmentFile + "." + std::to_string(mpi_rank); - std::string nexFileToReceive = baseNexusFile + "." + std::to_string(mpi_rank); - // Note that, unlike a bit earlier, don't need to check prior isGood in 1st receive, because not in loop - isGood = mpi_recv_text_file(catFileToReceive.c_str(), mpi_rank, sendingRank); - isGood = isGood && mpi_recv_text_file(nexFileToReceive.c_str(), mpi_rank, sendingRank); - } - } - - // Wait when appropriate - if (blockAll) { MPI_Barrier(MPI_COMM_WORLD); } - - // Sync status among the ranks also, if appropriate - if (syncReturnStatus) { - return mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "distributing subdivided hydrofabric files"); - } - // Otherwise, just return the local status value - else { - return isGood; - } - } - + const int *hostIdForRank, bool syncReturnStatus, bool blockAll); /** * Attempt to subdivide the passed hydrofabric files into a series of per-partition files. @@ -470,71 +202,8 @@ namespace parallel { * @return Whether subdividing was successful. */ bool subdivide_hydrofabric(int mpi_rank, int mpi_num_procs, const std::string &catchmentDataFile, - const std::string &nexusDataFile, const std::string &partitionConfigFile) - { - // Track whether things are good, meaning ok to continue and, at the end, whether successful - // Start with a value of true - bool isGood = true; - - #if !NGEN_WITH_PYTHON - // We can't be good to proceed with this, because Python is not active - isGood = false; - std::cerr << "Driver is unable to perform required hydrofabric subdividing when Python integration is not active." << std::endl; - - - // Sync with the rest of the ranks and bail if any aren't ready to proceed for any reason - if (!mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "initializing hydrofabric subdivider")) { - return false; - } - #else // i.e., #if NGEN_WITH_PYTHON - // Have rank 0 handle the generation task for all files/partitions - std::unique_ptr subdivider; - // Have rank 0 handle the generation task for all files/partitions - if (mpi_rank == 0) { - try { - subdivider = std::make_unique(catchmentDataFile, nexusDataFile, - partitionConfigFile); - } - catch (const std::exception &e) { - std::cerr << e.what() << std::endl; - // Set not good if the subdivider object couldn't be instantiated - isGood = false; - } - } - // Sync ranks and bail if any aren't ready to proceed for any reason - if (!mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "initializing hydrofabric subdivider")) { - return false; - } - - if (mpi_rank == 0) { - // Try to perform the subdividing - try { - isGood = subdivider->execSubdivision(); - } - catch (const std::exception &e) { - std::cerr << e.what() << std::endl; - // Set not good if the subdivider object couldn't be instantiated - isGood = false; - } - } - // Sync ranks again here on whether subdividing was successful, having them all exit at this point if not - if (!mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "executing hydrofabric subdivision")) { - return false; - } - // But if the subdividing went fine ... - else { - // ... figure out what ranks are on hosts with each other by getting an id for host of each rank - std::vector hostIdForRank(mpi_num_procs); - get_hosts_array(mpi_rank, mpi_num_procs, hostIdForRank.data()); - - // ... then (when necessary) transferring files around - return distribute_subdivided_hydrofabric_files(catchmentDataFile, nexusDataFile, 0, mpi_rank, - mpi_num_procs, hostIdForRank.data(), true, true); - - } - #endif // NGEN_WITH_PYTHON - } -} + const std::string &nexusDataFile, const std::string &partitionConfigFile); +} // namespace parallel #endif // NGEN_WITH_MPI diff --git a/src/NGen.cpp b/src/NGen.cpp index 1c0cf2bf28..d97145c7fc 100644 --- a/src/NGen.cpp +++ b/src/NGen.cpp @@ -38,8 +38,7 @@ std::string nexusDataFile = ""; std::string REALIZATION_CONFIG_PATH = ""; bool is_subdivided_hydrofabric_wanted = false; -// Define in the non-MPI case so that we don't need to conditionally compile `if (mpi_rank == 0)` -int mpi_rank = 0; +#include "parallel_utils.h" #if NGEN_WITH_MPI @@ -47,15 +46,12 @@ int mpi_rank = 0; #define MPI_HF_SUB_CLI_FLAG "--subdivided-hydrofabric" #endif -#include -#include "parallel_utils.h" #include "core/Partition_Parser.hpp" #include #include "core/Partition_One.hpp" std::string PARTITION_PATH = ""; -int mpi_num_procs; #endif // NGEN_WITH_MPI #include @@ -305,7 +301,7 @@ int main(int argc, char *argv[]) { // Do some extra steps if we expect to load a subdivided hydrofabric if (is_subdivided_hydrofabric_wanted) { // Ensure the hydrofabric is subdivided (either already or by doing it now), and then adjust these paths - if (parallel::is_hydrofabric_subdivided(mpi_rank, mpi_num_procs, true) || + if (parallel::is_hydrofabric_subdivided(mpi_rank, mpi_num_procs, catchmentDataFile, true) || parallel::subdivide_hydrofabric(mpi_rank, mpi_num_procs, catchmentDataFile, nexusDataFile, PARTITION_PATH)) { diff --git a/src/core/parallel_utils.cpp b/src/core/parallel_utils.cpp new file mode 100644 index 0000000000..72ce754384 --- /dev/null +++ b/src/core/parallel_utils.cpp @@ -0,0 +1,352 @@ +#include + +#include + +int mpi_rank = 0; + +#if NGEN_WITH_MPI + +int mpi_num_procs = 0; + +namespace parallel { + + bool mpiSyncStatusAnd(bool status, int mpi_rank, int mpi_num_procs, const std::string &taskDesc) { + // Expect 0 is good and 1 is no good for goodCode + // TODO: assert this in constructor or somewhere, or maybe just in a unit test + unsigned short codeBuffer; + bool printMessage = !taskDesc.empty(); + // For the other ranks, start by properly setting the status code value in the buffer and send to rank 0 + if (mpi_rank != 0) { + codeBuffer = status ? MPI_HF_SUB_CODE_GOOD : MPI_HF_SUB_CODE_BAD; + MPI_Send(&codeBuffer, 1, MPI_UNSIGNED_SHORT, 0, 0, MPI_COMM_WORLD); + } + // In rank 0, the first step is to receive and process codes from the other ranks into unified global status + else { + for (int i = 1; i < mpi_num_procs; ++i) { + MPI_Recv(&codeBuffer, 1, MPI_UNSIGNED_SHORT, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + // If any is ever "not good", overwrite status to be "false" + if (codeBuffer != MPI_HF_SUB_CODE_GOOD) { + if (printMessage) { + std::cout << "Rank " << i << " not successful/ready after " << taskDesc << std::endl; + } + status = false; + } + } + // Rank 0 must also now prepare the codeBuffer value for broadcasting the global status + codeBuffer = status ? MPI_HF_SUB_CODE_GOOD : MPI_HF_SUB_CODE_BAD; + } + + // Execute broadcast of global status rooted at rank 0 + MPI_Bcast(&codeBuffer, 1, MPI_UNSIGNED_SHORT, 0, MPI_COMM_WORLD); + return codeBuffer == MPI_HF_SUB_CODE_GOOD; + } + + bool is_hydrofabric_subdivided(int mpi_rank, int mpi_num_procs, std::string const& catchmentDataFile, bool printMsg) { + std::string name = catchmentDataFile + "." + std::to_string(mpi_rank); + // Initialize isGood based on local state. Here, local file is "good" when it already exists. + // TODO: this isn't actually checking whether the files are right (just that they are present) so do we need to? + bool isGood = utils::FileChecker::file_is_readable(name); + + if (mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs)) { + if (printMsg) { std::cout << "Process " << mpi_rank << ": Hydrofabric already subdivided in " << mpi_num_procs << " files." << std::endl; } + return true; + } + else { + if (printMsg) { std::cout << "Process " << mpi_rank << ": Hydrofabric has not yet been subdivided." << std::endl; } + return false; + } + } + + void get_hosts_array(int mpi_rank, int mpi_num_procs, int *host_array) { + const int ROOT_RANK = 0; + // These are the lengths of the (trimmed) C-string representations of the hostname for each rank + std::vector actualHostnameCStrLength(mpi_num_procs); + // Initialize to -1 to represent unknown + for (int i = 0; i < mpi_num_procs; ++i) { + actualHostnameCStrLength[i] = -1; + } + + // Get this rank's hostname (things should never be longer than 256) + char myhostname[256]; + gethostname(myhostname, 256); + + // Set the one for this rank + actualHostnameCStrLength[mpi_rank] = std::strlen(myhostname); + + // First, gather the hostname string lengths + MPI_Gather(&actualHostnameCStrLength[mpi_rank], 1, MPI_INT, actualHostnameCStrLength.data(), 1, MPI_INT, ROOT_RANK, + MPI_COMM_WORLD); + // Per-rank start offsets/displacements in our hostname strings gather buffer. + std::vector recvDisplacements(mpi_num_procs); + int totalLength = 0; + for (int i = 0; i < mpi_num_procs; ++i) { + // Displace each rank's string by the sum of the length of the previous rank strings + recvDisplacements[i] = totalLength; + // Adding the extra to make space for the \0 + actualHostnameCStrLength[i] += 1; + // Then update the total length + totalLength += actualHostnameCStrLength[i]; + + } + // Now we can create our buffer array and gather the hostname strings into it + char hostnames[totalLength]; + MPI_Gatherv(myhostname, actualHostnameCStrLength[mpi_rank], MPI_CHAR, hostnames, actualHostnameCStrLength.data(), + recvDisplacements.data(), MPI_CHAR, ROOT_RANK, MPI_COMM_WORLD); + + if (mpi_rank == ROOT_RANK) { + host_array[0] = 0; + int next_host_id = 1; + + int rank_with_matching_hostname; + char *checked_rank_hostname, *known_rank_hostname; + + for (int rank_being_check = 0; rank_being_check < mpi_num_procs; ++rank_being_check) { + // Set this as negative initially for each rank check to indicate no match found (at least yet) + rank_with_matching_hostname = -1; + // Get a C-string pointer for this rank's hostname, offset by the appropriate displacement + checked_rank_hostname = &hostnames[recvDisplacements[rank_being_check]]; + + // Assume that hostnames for any ranks less than the current rank being check are already known + for (int known_rank = 0; known_rank < rank_being_check; ++known_rank) { + // Get the right C-string pointer for the current known rank's hostname also + known_rank_hostname = &hostnames[recvDisplacements[known_rank]]; + // Compare the hostnames, setting and breaking if a match is found + if (std::strcmp(known_rank_hostname, checked_rank_hostname) == 0) { + rank_with_matching_hostname = known_rank; + break; + } + } + // This indicates this rank had no earlier rank with a matching hostname. + if (rank_with_matching_hostname < 0) { + // Assign new host id, then increment what the next id will be + host_array[rank_being_check] = next_host_id++; + } + else { + host_array[rank_being_check] = host_array[rank_with_matching_hostname]; + } + } + } + // Now, broadcast the results out + MPI_Bcast(host_array, mpi_num_procs, MPI_INT, 0, MPI_COMM_WORLD); + } + + bool mpi_send_text_file(const char *fileName, const int mpi_rank, const int destRank) { + int bufSize = 4096; + std::vector buf(bufSize); + int code; + // How much has been transferred so far + int totalNumTransferred = 0; + + FILE *file = fopen(fileName, "r"); + + // Transmit error code instead of expected size and return false if file can't be opened + if (file == NULL) { + // TODO: output error message + code = -1; + MPI_Send(&code, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + return false; + } + + // Send expected size to start + MPI_Send(&bufSize, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + + // Then get back expected size to infer other side is good to go + MPI_Recv(&code, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + if (code != bufSize) { + // TODO: output error message + fclose(file); + return false; + } + int continueCode = 1; + // Then while there is more of the file to read and send, read the next batch and ... + while (fgets(buf.data(), bufSize, file) != NULL) { + // Indicate we are ready to continue sending data + MPI_Send(&continueCode, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + + // Send this batch + MPI_Send(buf.data(), bufSize, MPI_CHAR, destRank, NGEN_MPI_DATA_TAG, MPI_COMM_WORLD); + + // Then get back a code, which will be -1 if bad and need to exit and otherwise good + MPI_Recv(&code, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + if (code < 0) { + // TODO: output error message + fclose(file); + return false; + } + } + // Once there is no more file to read and send, we should stop continuing + continueCode = 0; + MPI_Send(&continueCode, 1, MPI_INT, destRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + // Expect to get back a code of 0 + MPI_Recv(&code, 1, MPI_INT, destRank, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + fclose(file); + return code == 0; + } + + bool mpi_recv_text_file(const char *fileName, const int mpi_rank, const int srcRank) { + int bufSize, writeCode; + // Receive expected buffer size to start + MPI_Recv(&bufSize, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + + // If the sending side couldn't open the file, then immediately return false + if (bufSize == -1) { + // TODO: output error + return false; + } + + // Try to open recv file ... + FILE *file = fopen(fileName, "w"); + // ... and let sending size know whether this was successful by sending error code if not ... + if (file == NULL) { + // TODO: output error message + bufSize = -1; + MPI_Send(&bufSize, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + return false; + } + + // Send back the received buffer it if file opened, confirming things are good to go for transfer + MPI_Send(&bufSize, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + + // How much has been transferred so far + int totalNumTransferred = 0; + std::vector buf(bufSize); + + int continueCode; + + while (true) { + // Make sure the other side wants to continue sending data + MPI_Recv(&continueCode, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + if (continueCode <= 0) + break; + + MPI_Recv(buf.data(), bufSize, MPI_CHAR, srcRank, NGEN_MPI_DATA_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + + writeCode = fputs(buf.data(), file); + MPI_Send(&writeCode, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + + if (writeCode < 0) { + fclose(file); + return false; + } + } + + fclose(file); + MPI_Send(&continueCode, 1, MPI_INT, srcRank, NGEN_MPI_PROTOCOL_TAG, MPI_COMM_WORLD); + return true; + } + + bool distribute_subdivided_hydrofabric_files(const std::string &baseCatchmentFile, const std::string &baseNexusFile, + const int sendingRank, const int mpi_rank, const int mpi_num_procs, + const int *hostIdForRank, bool syncReturnStatus, bool blockAll) + { + // Start with status as good + bool isGood = true; + // FIXME: For now, just have rank 0 send everything, but optimize with multiple procs or threads later + // Only need to process this if sending rank or a receiving ranks (i.e., not on same host as sending rank) + if (mpi_rank == sendingRank || hostIdForRank[mpi_rank] != hostIdForRank[sendingRank]) { + // Have the sending rank send out all files + if (mpi_rank == sendingRank) { + // In rank 0, for all the other ranks ... + for (int otherRank = 0; otherRank < mpi_num_procs; ++otherRank) { + // If another rank is on a different host (note that this covers otherRank == sendingRank case) ... + if (hostIdForRank[otherRank] != hostIdForRank[mpi_rank]) { + // ... then send that rank its rank-specific catchment and nexus files + std::string catFileToSend = baseCatchmentFile + "." + std::to_string(otherRank); + std::string nexFileToSend = baseNexusFile + "." + std::to_string(otherRank); + // Note that checking previous isGood is necessary here because of loop + isGood = isGood && mpi_send_text_file(catFileToSend.c_str(), mpi_rank, otherRank); + isGood = isGood && mpi_send_text_file(nexFileToSend.c_str(), mpi_rank, otherRank); + } + } + } + else { + // For a rank not on the same host as the sending rank, receive the transmitted file + std::string catFileToReceive = baseCatchmentFile + "." + std::to_string(mpi_rank); + std::string nexFileToReceive = baseNexusFile + "." + std::to_string(mpi_rank); + // Note that, unlike a bit earlier, don't need to check prior isGood in 1st receive, because not in loop + isGood = mpi_recv_text_file(catFileToReceive.c_str(), mpi_rank, sendingRank); + isGood = isGood && mpi_recv_text_file(nexFileToReceive.c_str(), mpi_rank, sendingRank); + } + } + + // Wait when appropriate + if (blockAll) { MPI_Barrier(MPI_COMM_WORLD); } + + // Sync status among the ranks also, if appropriate + if (syncReturnStatus) { + return mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "distributing subdivided hydrofabric files"); + } + // Otherwise, just return the local status value + else { + return isGood; + } + } + + bool subdivide_hydrofabric(int mpi_rank, int mpi_num_procs, const std::string &catchmentDataFile, + const std::string &nexusDataFile, const std::string &partitionConfigFile) + { + // Track whether things are good, meaning ok to continue and, at the end, whether successful + // Start with a value of true + bool isGood = true; + + #if !NGEN_WITH_PYTHON + // We can't be good to proceed with this, because Python is not active + isGood = false; + std::cerr << "Driver is unable to perform required hydrofabric subdividing when Python integration is not active." << std::endl; + + + // Sync with the rest of the ranks and bail if any aren't ready to proceed for any reason + if (!mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "initializing hydrofabric subdivider")) { + return false; + } + #else // i.e., #if NGEN_WITH_PYTHON + // Have rank 0 handle the generation task for all files/partitions + std::unique_ptr subdivider; + // Have rank 0 handle the generation task for all files/partitions + if (mpi_rank == 0) { + try { + subdivider = std::make_unique(catchmentDataFile, nexusDataFile, + partitionConfigFile); + } + catch (const std::exception &e) { + std::cerr << e.what() << std::endl; + // Set not good if the subdivider object couldn't be instantiated + isGood = false; + } + } + // Sync ranks and bail if any aren't ready to proceed for any reason + if (!mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "initializing hydrofabric subdivider")) { + return false; + } + + if (mpi_rank == 0) { + // Try to perform the subdividing + try { + isGood = subdivider->execSubdivision(); + } + catch (const std::exception &e) { + std::cerr << e.what() << std::endl; + // Set not good if the subdivider object couldn't be instantiated + isGood = false; + } + } + // Sync ranks again here on whether subdividing was successful, having them all exit at this point if not + if (!mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "executing hydrofabric subdivision")) { + return false; + } + // But if the subdividing went fine ... + else { + // ... figure out what ranks are on hosts with each other by getting an id for host of each rank + std::vector hostIdForRank(mpi_num_procs); + get_hosts_array(mpi_rank, mpi_num_procs, hostIdForRank.data()); + + // ... then (when necessary) transferring files around + return distribute_subdivided_hydrofabric_files(catchmentDataFile, nexusDataFile, 0, mpi_rank, + mpi_num_procs, hostIdForRank.data(), true, true); + + } + #endif // NGEN_WITH_PYTHON + } +} // namespace parallel + +#endif From b850742dbd83d9dca0673a8af7e05a95a717d07b Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 17:52:13 -0700 Subject: [PATCH 37/51] Initialize input variables --- src/realizations/coastal/SchismFormulation.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index abf4397c53..f42bfd4780 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -3,12 +3,10 @@ #if NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI #include -#include +#include const static auto s_schism_registration_function = "register_bmi"; - - std::map SchismFormulation::expected_input_variables_ = { /* Meteorological Forcings */ @@ -79,7 +77,7 @@ void SchismFormulation::initialize() input_variable_count_[name] = mesh_size(name); } - //set_inputs(); + set_inputs(); } void SchismFormulation::finalize() From 0b0a237ec05dbae4dd27f640d766391c924d7ca6 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 17:56:52 -0700 Subject: [PATCH 38/51] Finalize input providers since mock is now set --- src/realizations/coastal/SchismFormulation.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index f42bfd4780..2e745032b8 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -82,11 +82,10 @@ void SchismFormulation::initialize() void SchismFormulation::finalize() { -#if 0 meteorological_forcings_provider_->finalize(); offshore_boundary_provider_->finalize(); inflows_boundary_provider_->finalize(); -#endif + bmi_->Finalize(); } From d9bcb52e1e191d99231981f6ddcf32d1cd0d39df Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 26 Aug 2024 17:57:43 -0700 Subject: [PATCH 39/51] Single-step test with mocked inputs, works for mpiexec -np 1 --- test/CMakeLists.txt | 31 +++++----- test/coastal/SchismFormulation_Test.cpp | 81 ++++++++++++++++++++++++- 2 files changed, 97 insertions(+), 15 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 26e83a963d..e82cb096b5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -135,20 +135,23 @@ endfunction() ########################## Individual tests ########################## ########################## Coastal Formulation Tests -ngen_add_test( - test_coastal - OBJECTS - coastal/SchismFormulation_Test.cpp - LIBRARIES - NGen::core - NGen::realizations_catchment - NGen::core_mediator - NGen::forcing - NGen::ngen_bmi - REQUIRES - NGEN_WITH_BMI_FORTRAN - NGEN_WITH_MPI -) +# ngen_add_test( +# test_coastal +# OBJECTS +# coastal/SchismFormulation_Test.cpp +# LIBRARIES +# NGen::core +# NGen::realizations_catchment +# NGen::core_mediator +# NGen::forcing +# NGen::ngen_bmi +# REQUIRES +# NGEN_WITH_BMI_FORTRAN +# NGEN_WITH_MPI +# ) + +add_executable(test_coastal coastal/SchismFormulation_Test.cpp) +target_link_libraries(test_coastal NGen::core NGen::realizations_catchment NGen::realizations_coastal NGen::core_mediator NGen::forcing NGen::ngen_bmi) ########################## GeoJSON Unit Tests ngen_add_test( diff --git a/test/coastal/SchismFormulation_Test.cpp b/test/coastal/SchismFormulation_Test.cpp index c8e43f9dec..38caccd23e 100644 --- a/test/coastal/SchismFormulation_Test.cpp +++ b/test/coastal/SchismFormulation_Test.cpp @@ -1,9 +1,88 @@ -#include "gtest/gtest.h" +//#include "gtest/gtest.h" +#include + #include "realizations/coastal/SchismFormulation.hpp" +const static std::string library_path = "/Users/phil/Code/noaa/BUILD/bmischism_2024-08-19-ngen/libtestbmifortranmodel.dylib"; +const static std::string init_config_path = "/Users/phil/Code/noaa/SCHISM_Lake_Champlain_BMI_test/namelist.input"; + +#if 0 struct Schism_Formulation_IT : public ::testing::Test { }; +#endif + +std::map input_variables_defaults = + { + /* Meteorological Forcings */ + // RAINRATE - precipitation + {"RAINRATE", 0.01}, + // SFCPRS - surface atmospheric pressure + {"SFCPRS", 101325.0}, + // SPFH2m - specific humidity at 2m + {"SPFH2m", 0.01}, + // TMP2m - temperature at 2m + {"TMP2m", 293}, + // UU10m, VV10m - wind velocity components at 10m + {"UU10m", 1.0}, + {"VV10m", 1.0}, + + /* Input Boundary Conditions */ + // ETA2_bnd - water surface elevation at the boundaries + {"ETA2_bnd", 30}, + // Q_bnd - flows at boundaries + {"Q_bnd", 0.1}, + }; + +struct MockProvider : data_access::DataProvider +{ + std::vector data; + + MockProvider() + : data(552697, 0.0) + {} + ~MockProvider() = default; + + // Implementation of DataProvider + std::vector variables; + boost::span get_available_variable_names() const override { return variables; } + + long get_data_start_time() const override { return 0; } + long get_data_stop_time() const override { return 0; } + long record_duration() const override { return 3600; } + size_t get_ts_index_for_time(const time_t &epoch_time) const override { return 1; } + + data_type get_value(const selection_type& selector, data_access::ReSampleMethod m) override { return data[0]; } + std::vector get_values(const selection_type& selector, data_access::ReSampleMethod) override + { + auto default_value = input_variables_defaults[selector.variable_name]; + for (auto& val : data) { + val = default_value; + } + return data; + } +}; + +int main(int argc, char **argv) +{ + MPI_Init(&argc, &argv); + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + + auto provider = std::make_shared(); + auto schism = std::make_unique(/*id=*/ "test_schism_formulation", + library_path, + init_config_path, + provider, + provider, + provider + ); + + schism->initialize(); + schism->update(); + schism->finalize(); + MPI_Finalize(); + return 0; +} From b1907972c98aa7771b466ab3020ff29d05cb1956 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 12:39:19 -0700 Subject: [PATCH 40/51] Annotate override --- include/realizations/coastal/CoastalFormulation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/realizations/coastal/CoastalFormulation.hpp b/include/realizations/coastal/CoastalFormulation.hpp index 37b3fd04f4..be4a32d0a6 100644 --- a/include/realizations/coastal/CoastalFormulation.hpp +++ b/include/realizations/coastal/CoastalFormulation.hpp @@ -20,7 +20,7 @@ class CoastalFormulation : public data_access::DataProvider Date: Tue, 17 Sep 2024 12:40:01 -0700 Subject: [PATCH 41/51] Firm up MeshPointsSelector --- include/forcing/GenericDataProvider.hpp | 2 ++ include/forcing/MeshPointsSelectors.hpp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/include/forcing/GenericDataProvider.hpp b/include/forcing/GenericDataProvider.hpp index d4c6bd7ae9..ed8dda9191 100644 --- a/include/forcing/GenericDataProvider.hpp +++ b/include/forcing/GenericDataProvider.hpp @@ -3,10 +3,12 @@ #include "DataProvider.hpp" #include "DataProviderSelectors.hpp" +#include "MeshPointsSelectors.hpp" namespace data_access { using GenericDataProvider = DataProvider; + using MeshPointsDataProvider = DataProvider; } #endif diff --git a/include/forcing/MeshPointsSelectors.hpp b/include/forcing/MeshPointsSelectors.hpp index 680f882388..55e9c639d2 100644 --- a/include/forcing/MeshPointsSelectors.hpp +++ b/include/forcing/MeshPointsSelectors.hpp @@ -10,7 +10,7 @@ static AllPoints all_points; struct MeshPointsSelector { std::string variable_name; - std::chrono::time_point init_time; + std::chrono::time_point init_time; std::chrono::seconds duration; std::string output_units; boost::variant> points; From eaac98751b428dcb5ab4195b656cba24ee58e5c4 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 12:41:21 -0700 Subject: [PATCH 42/51] Mostly implemented NetCDF unstructured mesh reader for coastal formulation use --- .../forcing/NetCDFMeshPointsDataProvider.hpp | 156 +++++++++++++ src/forcing/CMakeLists.txt | 5 +- src/forcing/NetCDFMeshPointsDataProvider.cpp | 221 ++++++++++++++++++ 3 files changed, 381 insertions(+), 1 deletion(-) create mode 100644 include/forcing/NetCDFMeshPointsDataProvider.hpp create mode 100644 src/forcing/NetCDFMeshPointsDataProvider.cpp diff --git a/include/forcing/NetCDFMeshPointsDataProvider.hpp b/include/forcing/NetCDFMeshPointsDataProvider.hpp new file mode 100644 index 0000000000..e8895c60c0 --- /dev/null +++ b/include/forcing/NetCDFMeshPointsDataProvider.hpp @@ -0,0 +1,156 @@ +#pragma once + +#include + +#if NGEN_WITH_NETCDF + +#include "GenericDataProvider.hpp" +#include "DataProviderSelectors.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include "assert.h" +#include +#include + +#include +#include + +#include "AorcForcing.hpp" + +namespace netCDF { + class NcVar; + class NcFile; +} + +namespace data_access +{ + class NetCDFMeshPointsDataProvider : public MeshPointsDataProvider + { + public: + + enum TimeUnit + { + TIME_HOURS, + TIME_MINUTES, + TIME_SECONDS, + TIME_MILLISECONDS, + TIME_MICROSECONDS, + TIME_NANOSECONDS + }; + + using time_point_type = std::chrono::time_point; + + /** + * @brief Factory method that creates or returns an existing provider for the provided path. + * @param input_path The path to a NetCDF file with lumped catchment forcing values. + * @param log_s An output log stream for messages from the underlying library. If a provider object for + * the given path already exists, this argument will be ignored. + */ + static std::shared_ptr get_shared_provider(std::string input_path, time_point_type sim_start, time_point_type sim_end); + + /** + * @brief Cleanup the shared providers cache, ensuring that the files get closed. + */ + static void cleanup_shared_providers(); + + NetCDFMeshPointsDataProvider(std::string input_path, + time_point_type sim_start, + time_point_type sim_end); + + // Default implementation defined in the .cpp file so that + // client code doesn't need to have the full definition of + // NcFile visible for the compiler to implicitly generate + // ~NetCDFMeshPointsDataProvider() = default; + // for every file that uses this class + ~NetCDFMeshPointsDataProvider(); + + void finalize() override; + + /** Return the variables that are accessable by this data provider */ + boost::span get_available_variable_names() const override; + + /** Return the first valid time for which data from the request variable can be requested */ + long get_data_start_time() const override; + + /** Return the last valid time for which data from the requested variable can be requested */ + long get_data_stop_time() const override; + + long record_duration() const override; + + /** + * Get the index of the data time step that contains the given point in time. + * + * An @ref std::out_of_range exception should be thrown if the time is not in any time step. + * + * @param epoch_time The point in time, as a seconds-based epoch time. + * @return The index of the forcing time step that contains the given point in time. + * @throws std::out_of_range If the given point is not in any time step. + */ + size_t get_ts_index_for_time(const time_t &epoch_time) const override; + + /** + * Get the value of a forcing property for an arbitrary time period, converting units if needed. + * + * An @ref std::out_of_range exception should be thrown if the data for the time period is not available. + * + * @param selector Data required to establish what subset of the stored data should be accessed + * @param m How data is to be resampled if there is a mismatch in data alignment or repeat rate + * @return The value of the forcing property for the described time period, with units converted if needed. + * @throws std::out_of_range If data for the time period is not available. + */ + double get_value(const selection_type& selector, ReSampleMethod m) override; + + // The interface that DataProvider really should have + virtual void get_values(const selection_type& selector, boost::span data); + + // And an implementation of the usual version using it + std::vector get_values(const selection_type& selector, data_access::ReSampleMethod) override + { + std::vector output(selected_points_count(selector)); + get_values(selector, output); + return output; + } + + private: + + size_t mesh_size(std::string const& variable_name); + + size_t selected_points_count(const selection_type& selector) + { + auto* points = boost::get>(&selector.points); + size_t size = points ? points->size() : this->mesh_size(selector.variable_name); + return size; + } + + time_point_type sim_start_date_time_epoch; + time_point_type sim_end_date_time_epoch; + std::chrono::seconds sim_to_data_time_offset; // Deliberately signed--sim should never start before data, yes? + + static std::mutex shared_providers_mutex; + static std::map> shared_providers; + + std::vector variable_names; + std::vector time_vals; + TimeUnit time_unit; // the unit that time was stored as in the file + std::chrono::seconds time_stride; // the amount of time between stored time values + + std::shared_ptr nc_file; + + struct metadata_cache_entry; + std::map ncvar_cache; + + boost::compute::detail::lru_cache>> value_cache; + size_t cache_slice_t_size = 1; + size_t cache_slice_c_size = 1; + }; +} + + +#endif // NGEN_WITH_NETCDF diff --git a/src/forcing/CMakeLists.txt b/src/forcing/CMakeLists.txt index d374032b2e..8cf40e198f 100644 --- a/src/forcing/CMakeLists.txt +++ b/src/forcing/CMakeLists.txt @@ -23,7 +23,10 @@ target_link_libraries(forcing PUBLIC target_sources(forcing PRIVATE "${CMAKE_CURRENT_LIST_DIR}/NullForcingProvider.cpp") if(NGEN_WITH_NETCDF) - target_sources(forcing PRIVATE "${CMAKE_CURRENT_LIST_DIR}/NetCDFPerFeatureDataProvider.cpp") + target_sources(forcing PRIVATE + "${CMAKE_CURRENT_LIST_DIR}/NetCDFPerFeatureDataProvider.cpp" + "${CMAKE_CURRENT_LIST_DIR}/NetCDFMeshPointsDataProvider.cpp" + ) target_link_libraries(forcing PUBLIC NetCDF) endif() diff --git a/src/forcing/NetCDFMeshPointsDataProvider.cpp b/src/forcing/NetCDFMeshPointsDataProvider.cpp new file mode 100644 index 0000000000..4356ab4020 --- /dev/null +++ b/src/forcing/NetCDFMeshPointsDataProvider.cpp @@ -0,0 +1,221 @@ +#include + +#if NGEN_WITH_NETCDF +#include "NetCDFMeshPointsDataProvider.hpp" +#include +#include + +std::mutex data_access::NetCDFMeshPointsDataProvider::shared_providers_mutex; +std::map> data_access::NetCDFMeshPointsDataProvider::shared_providers; + +namespace data_access { + +// Out-of-line class definition after forward-declaration so that the +// header doesn't need #include for NcVar to be a complete +// type there +struct NetCDFMeshPointsDataProvider::metadata_cache_entry { + netCDF::NcVar variable; + std::string units; + double scale_factor; + double offset; +}; + +std::shared_ptr NetCDFMeshPointsDataProvider::get_shared_provider(std::string input_path, time_point_type sim_start, time_point_type sim_end) +{ + const std::lock_guard lock(shared_providers_mutex); + std::shared_ptr p; + if(shared_providers.count(input_path) > 0){ + p = shared_providers[input_path]; + } else { + p = std::make_shared(input_path, sim_start, sim_end); + shared_providers[input_path] = p; + } + return p; +} + +void NetCDFMeshPointsDataProvider::cleanup_shared_providers() +{ + const std::lock_guard lock(shared_providers_mutex); + // First lets try just emptying the map... if all goes well, everything will destruct properly on its own... + shared_providers.clear(); +} + +NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_path, time_point_type sim_start, time_point_type sim_end) + : value_cache(20), + sim_start_date_time_epoch(sim_start), + sim_end_date_time_epoch(sim_end) +{ + nc_file = std::make_shared(input_path, netCDF::NcFile::read); + + auto var_set = nc_file->getVars(); + + // Populate the variable attributes cache + for (const auto& element : var_set) + { + std::string var_name = element.first; + auto ncvar = nc_file->getVar(var_name); + variable_names.push_back(var_name); + + std::string native_units; + auto units_att = ncvar.getAtt("units"); + units_att.getValues(native_units); + + double scale_factor = 1.0; + auto scale_var = ncvar.getAtt("scale_factor"); + if (!scale_var.isNull()) { + scale_var.getValues(&scale_factor); + } + + double offset = 0.0; + auto offset_var = ncvar.getAtt("add_offset"); + if (!offset_var.isNull()) { + offset_var.getValues(&offset); + } + + ncvar_cache[var_name] = {ncvar, native_units, scale_factor, offset}; + } + + auto num_times = nc_file->getDim("time").getSize(); + + std::vector>> raw_time(num_times); + + auto time_var = nc_file->getVar("Time"); + + if (time_var.getDimCount() != 1) { + throw "'Time' variable has dimension other than 1"; + } + + time_var.getVar(raw_time.data()); + + auto time_unit_att = time_var.getAtt("units"); + std::string time_unit_str; + + time_unit_att.getValues(time_unit_str); + if (time_unit_str != "minutes since 1970-01-01 00:00:00 UTC") { + throw "Time units not exactly as expected"; + } + + time_vals.reserve(num_times); + + for (int i = 0; i < num_times; ++i) { + // Assume that the system clock's epoch matches Unix time. + // This is guaranteed from C++20 onwards + time_vals.push_back(time_point_type(std::chrono::duration_cast(raw_time[i]))); + } + + time_stride = std::chrono::duration_cast(time_vals[1] - time_vals[0]); + +#if 0 + // verify the time stride + for( size_t i = 1; i < time_vals.size() -1; ++i) + { + auto tinterval = time_vals[i+1] - time_vals[i]; + + if ( tinterval - time_stride > 1us) + { + throw std::runtime_error("Time intervals in forcing file are not constant"); + } + } + + // determine start_time and stop_time; + start_time = epoch_start_time + time_vals[0]; + stop_time = epoch_start_time + time_vals.back() + time_stride; + + sim_to_data_time_offset = std::chrono::duration_cast(sim_start_date_time_epoch - start_time); +#endif +} + +NetCDFMeshPointsDataProvider::~NetCDFMeshPointsDataProvider() = default; + +void NetCDFMeshPointsDataProvider::finalize() +{ + ncvar_cache.clear(); + + if (nc_file != nullptr) { + nc_file->close(); + } + nc_file = nullptr; +} + +boost::span NetCDFMeshPointsDataProvider::get_available_variable_names() const +{ + return variable_names; +} + +/** Return the first valid time for which data from the request variable can be requested */ +long NetCDFMeshPointsDataProvider::get_data_start_time() const +{ + return std::chrono::system_clock::to_time_t(time_vals[0]); +#if 0 + //return start_time; + //FIXME: Matching behavior from CsvMeshPointsForcingProvider, but both are probably wrong! + return sim_start_date_time_epoch.time_since_epoch().count(); // return start_time + sim_to_data_time_offset; +#endif +} + +/** Return the last valid time for which data from the requested variable can be requested */ +long NetCDFMeshPointsDataProvider::get_data_stop_time() const +{ + return std::chrono::system_clock::to_time_t(time_vals.back()) + std::chrono::duration_cast(time_stride).count(); +#if 0 + //return stop_time; + //FIXME: Matching behavior from CsvMeshPointsForcingProvider, but both are probably wrong! + return sim_end_date_time_epoch.time_since_epoch().count(); // return end_time + sim_to_data_time_offset; +#endif +} + +long NetCDFMeshPointsDataProvider::record_duration() const +{ + return std::chrono::duration_cast(time_stride).count(); +} + +size_t NetCDFMeshPointsDataProvider::get_ts_index_for_time(const time_t &epoch_time_in) const +{ + auto start_time = time_vals.front(); + auto stop_time = time_vals.back() + time_stride; + + auto epoch_time = std::chrono::system_clock::from_time_t(epoch_time_in); + + if (start_time <= epoch_time && epoch_time < stop_time) + { + auto offset = epoch_time - start_time; + auto index = offset / time_stride; + return index; + } + else + { + //std::stringstream ss; + //ss << "The value " << (int)epoch_time << " was not in the range [" << (int)start_time << "," << (int)stop_time << ")\n" << SOURCE_LOC; + //throw std::out_of_range(ss.str().c_str()); + throw std::out_of_range(""); + } +} + +void NetCDFMeshPointsDataProvider::get_values(const selection_type& selector, boost::span data) +{ + if (!boost::get(&selector.points)) throw "Not implemented - only all_points"; + + auto const& metadata = ncvar_cache[selector.variable_name]; + std::string const& source_units = metadata.units; + + size_t time_index = get_ts_index_for_time(std::chrono::system_clock::to_time_t(selector.init_time)); + + metadata.variable.getVar({time_index, 0}, {1, data.size()}, data.data()); + + for (double& value : data) { + value = value * metadata.scale_factor + metadata.offset; + } + + UnitsHelper::convert_values(source_units, data.data(), selector.output_units, data.data(), data.size()); +} + +double NetCDFMeshPointsDataProvider::get_value(const selection_type& selector, ReSampleMethod m) +{ + throw std::runtime_error("Not implemented - access chunks of the mesh"); + + return 0.0; +} + +} + +#endif From aab3cdc73257ccbbc9abb699ab7c015b41a5c554 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 12:50:23 -0700 Subject: [PATCH 43/51] Drop unused member variable --- include/forcing/NetCDFMeshPointsDataProvider.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/forcing/NetCDFMeshPointsDataProvider.hpp b/include/forcing/NetCDFMeshPointsDataProvider.hpp index e8895c60c0..8078a7b6d8 100644 --- a/include/forcing/NetCDFMeshPointsDataProvider.hpp +++ b/include/forcing/NetCDFMeshPointsDataProvider.hpp @@ -138,7 +138,6 @@ namespace data_access std::vector variable_names; std::vector time_vals; - TimeUnit time_unit; // the unit that time was stored as in the file std::chrono::seconds time_stride; // the amount of time between stored time values std::shared_ptr nc_file; From 9fd9510a1359065633f945d448d59e0d6c1da504 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 14:40:02 -0700 Subject: [PATCH 44/51] Remove unused instance and value caching from NetCDF mesh provider --- .../forcing/NetCDFMeshPointsDataProvider.hpp | 20 --------------- src/forcing/NetCDFMeshPointsDataProvider.cpp | 25 +------------------ 2 files changed, 1 insertion(+), 44 deletions(-) diff --git a/include/forcing/NetCDFMeshPointsDataProvider.hpp b/include/forcing/NetCDFMeshPointsDataProvider.hpp index 8078a7b6d8..6c2aa100af 100644 --- a/include/forcing/NetCDFMeshPointsDataProvider.hpp +++ b/include/forcing/NetCDFMeshPointsDataProvider.hpp @@ -47,19 +47,6 @@ namespace data_access using time_point_type = std::chrono::time_point; - /** - * @brief Factory method that creates or returns an existing provider for the provided path. - * @param input_path The path to a NetCDF file with lumped catchment forcing values. - * @param log_s An output log stream for messages from the underlying library. If a provider object for - * the given path already exists, this argument will be ignored. - */ - static std::shared_ptr get_shared_provider(std::string input_path, time_point_type sim_start, time_point_type sim_end); - - /** - * @brief Cleanup the shared providers cache, ensuring that the files get closed. - */ - static void cleanup_shared_providers(); - NetCDFMeshPointsDataProvider(std::string input_path, time_point_type sim_start, time_point_type sim_end); @@ -133,9 +120,6 @@ namespace data_access time_point_type sim_end_date_time_epoch; std::chrono::seconds sim_to_data_time_offset; // Deliberately signed--sim should never start before data, yes? - static std::mutex shared_providers_mutex; - static std::map> shared_providers; - std::vector variable_names; std::vector time_vals; std::chrono::seconds time_stride; // the amount of time between stored time values @@ -144,10 +128,6 @@ namespace data_access struct metadata_cache_entry; std::map ncvar_cache; - - boost::compute::detail::lru_cache>> value_cache; - size_t cache_slice_t_size = 1; - size_t cache_slice_c_size = 1; }; } diff --git a/src/forcing/NetCDFMeshPointsDataProvider.cpp b/src/forcing/NetCDFMeshPointsDataProvider.cpp index 4356ab4020..ad477c9875 100644 --- a/src/forcing/NetCDFMeshPointsDataProvider.cpp +++ b/src/forcing/NetCDFMeshPointsDataProvider.cpp @@ -5,9 +5,6 @@ #include #include -std::mutex data_access::NetCDFMeshPointsDataProvider::shared_providers_mutex; -std::map> data_access::NetCDFMeshPointsDataProvider::shared_providers; - namespace data_access { // Out-of-line class definition after forward-declaration so that the @@ -20,28 +17,8 @@ struct NetCDFMeshPointsDataProvider::metadata_cache_entry { double offset; }; -std::shared_ptr NetCDFMeshPointsDataProvider::get_shared_provider(std::string input_path, time_point_type sim_start, time_point_type sim_end) -{ - const std::lock_guard lock(shared_providers_mutex); - std::shared_ptr p; - if(shared_providers.count(input_path) > 0){ - p = shared_providers[input_path]; - } else { - p = std::make_shared(input_path, sim_start, sim_end); - shared_providers[input_path] = p; - } - return p; -} - -void NetCDFMeshPointsDataProvider::cleanup_shared_providers() -{ - const std::lock_guard lock(shared_providers_mutex); - // First lets try just emptying the map... if all goes well, everything will destruct properly on its own... - shared_providers.clear(); -} - NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_path, time_point_type sim_start, time_point_type sim_end) - : value_cache(20), + : sim_start_date_time_epoch(sim_start), sim_end_date_time_epoch(sim_end) { From 5ed7fd5282e383776c8757ce60f1938ba99776dd Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 14:45:47 -0700 Subject: [PATCH 45/51] Clear out some residual cruft --- .../forcing/NetCDFMeshPointsDataProvider.hpp | 27 +++---------------- src/forcing/NetCDFMeshPointsDataProvider.cpp | 2 -- 2 files changed, 3 insertions(+), 26 deletions(-) diff --git a/include/forcing/NetCDFMeshPointsDataProvider.hpp b/include/forcing/NetCDFMeshPointsDataProvider.hpp index 6c2aa100af..0f63d7e404 100644 --- a/include/forcing/NetCDFMeshPointsDataProvider.hpp +++ b/include/forcing/NetCDFMeshPointsDataProvider.hpp @@ -35,16 +35,6 @@ namespace data_access { public: - enum TimeUnit - { - TIME_HOURS, - TIME_MINUTES, - TIME_SECONDS, - TIME_MILLISECONDS, - TIME_MICROSECONDS, - TIME_NANOSECONDS - }; - using time_point_type = std::chrono::time_point; NetCDFMeshPointsDataProvider(std::string input_path, @@ -60,10 +50,10 @@ namespace data_access void finalize() override; - /** Return the variables that are accessable by this data provider */ + /** Return the variables that are accessible by this data provider */ boost::span get_available_variable_names() const override; - /** Return the first valid time for which data from the request variable can be requested */ + /** Return the first valid time for which data from the requested variable can be requested */ long get_data_start_time() const override; /** Return the last valid time for which data from the requested variable can be requested */ @@ -100,22 +90,11 @@ namespace data_access // And an implementation of the usual version using it std::vector get_values(const selection_type& selector, data_access::ReSampleMethod) override { - std::vector output(selected_points_count(selector)); - get_values(selector, output); - return output; + throw std::runtime_error("Unimplemented"); } private: - size_t mesh_size(std::string const& variable_name); - - size_t selected_points_count(const selection_type& selector) - { - auto* points = boost::get>(&selector.points); - size_t size = points ? points->size() : this->mesh_size(selector.variable_name); - return size; - } - time_point_type sim_start_date_time_epoch; time_point_type sim_end_date_time_epoch; std::chrono::seconds sim_to_data_time_offset; // Deliberately signed--sim should never start before data, yes? diff --git a/src/forcing/NetCDFMeshPointsDataProvider.cpp b/src/forcing/NetCDFMeshPointsDataProvider.cpp index ad477c9875..385c4502ff 100644 --- a/src/forcing/NetCDFMeshPointsDataProvider.cpp +++ b/src/forcing/NetCDFMeshPointsDataProvider.cpp @@ -119,7 +119,6 @@ boost::span NetCDFMeshPointsDataProvider::get_available_varia return variable_names; } -/** Return the first valid time for which data from the request variable can be requested */ long NetCDFMeshPointsDataProvider::get_data_start_time() const { return std::chrono::system_clock::to_time_t(time_vals[0]); @@ -130,7 +129,6 @@ long NetCDFMeshPointsDataProvider::get_data_start_time() const #endif } -/** Return the last valid time for which data from the requested variable can be requested */ long NetCDFMeshPointsDataProvider::get_data_stop_time() const { return std::chrono::system_clock::to_time_t(time_vals.back()) + std::chrono::duration_cast(time_stride).count(); From 852c31cfa4a55fa105178dca1f79bd34f5b1f3fa Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 16:30:03 -0700 Subject: [PATCH 46/51] DataProvider: Add nicer interface to base class --- include/forcing/DataProvider.hpp | 2 ++ include/forcing/NetCDFMeshPointsDataProvider.hpp | 3 +-- include/realizations/coastal/CoastalFormulation.hpp | 3 +-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/forcing/DataProvider.hpp b/include/forcing/DataProvider.hpp index 32cbb36029..8d3da0ff5a 100644 --- a/include/forcing/DataProvider.hpp +++ b/include/forcing/DataProvider.hpp @@ -96,6 +96,8 @@ namespace data_access */ virtual std::vector get_values(const selection_type& selector, ReSampleMethod m=SUM) = 0; + virtual void get_values(const selection_type& selector, boost::span data) { throw std::runtime_error("DP::get_values(span) Unimplemented"); } + virtual bool is_property_sum_over_time_step(const std::string& name) const {return false; } private: diff --git a/include/forcing/NetCDFMeshPointsDataProvider.hpp b/include/forcing/NetCDFMeshPointsDataProvider.hpp index 0f63d7e404..e608a46b76 100644 --- a/include/forcing/NetCDFMeshPointsDataProvider.hpp +++ b/include/forcing/NetCDFMeshPointsDataProvider.hpp @@ -84,8 +84,7 @@ namespace data_access */ double get_value(const selection_type& selector, ReSampleMethod m) override; - // The interface that DataProvider really should have - virtual void get_values(const selection_type& selector, boost::span data); + void get_values(const selection_type& selector, boost::span data) override; // And an implementation of the usual version using it std::vector get_values(const selection_type& selector, data_access::ReSampleMethod) override diff --git a/include/realizations/coastal/CoastalFormulation.hpp b/include/realizations/coastal/CoastalFormulation.hpp index be4a32d0a6..8f2e4369dd 100644 --- a/include/realizations/coastal/CoastalFormulation.hpp +++ b/include/realizations/coastal/CoastalFormulation.hpp @@ -23,8 +23,7 @@ class CoastalFormulation : public data_access::DataProvider data) = 0; + void get_values(const selection_type& selector, boost::span data) override = 0; // And an implementation of the usual version using it std::vector get_values(const selection_type& selector, data_access::ReSampleMethod) override From 5d2ed334ffaad3efc6bba692b8e98bcababe3c64 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 16:58:46 -0700 Subject: [PATCH 47/51] NetCDF mesh data provider: fix up time keeping --- .../forcing/NetCDFMeshPointsDataProvider.hpp | 13 +- src/forcing/NetCDFMeshPointsDataProvider.cpp | 122 ++++++++++-------- 2 files changed, 71 insertions(+), 64 deletions(-) diff --git a/include/forcing/NetCDFMeshPointsDataProvider.hpp b/include/forcing/NetCDFMeshPointsDataProvider.hpp index e608a46b76..fe4530c07c 100644 --- a/include/forcing/NetCDFMeshPointsDataProvider.hpp +++ b/include/forcing/NetCDFMeshPointsDataProvider.hpp @@ -12,17 +12,7 @@ #include #include #include -#include #include -#include -#include "assert.h" -#include -#include - -#include -#include - -#include "AorcForcing.hpp" namespace netCDF { class NcVar; @@ -94,9 +84,10 @@ namespace data_access private: + void cache_variable(std::string const& var_name); + time_point_type sim_start_date_time_epoch; time_point_type sim_end_date_time_epoch; - std::chrono::seconds sim_to_data_time_offset; // Deliberately signed--sim should never start before data, yes? std::vector variable_names; std::vector time_vals; diff --git a/src/forcing/NetCDFMeshPointsDataProvider.cpp b/src/forcing/NetCDFMeshPointsDataProvider.cpp index 385c4502ff..c7f4c88a21 100644 --- a/src/forcing/NetCDFMeshPointsDataProvider.cpp +++ b/src/forcing/NetCDFMeshPointsDataProvider.cpp @@ -2,56 +2,31 @@ #if NGEN_WITH_NETCDF #include "NetCDFMeshPointsDataProvider.hpp" -#include +#include + #include +#include +#include + namespace data_access { // Out-of-line class definition after forward-declaration so that the // header doesn't need #include for NcVar to be a complete // type there struct NetCDFMeshPointsDataProvider::metadata_cache_entry { - netCDF::NcVar variable; + netCDF::NcVar ncVar; std::string units; double scale_factor; double offset; }; NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_path, time_point_type sim_start, time_point_type sim_end) - : - sim_start_date_time_epoch(sim_start), - sim_end_date_time_epoch(sim_end) + : sim_start_date_time_epoch(sim_start) + , sim_end_date_time_epoch(sim_end) { nc_file = std::make_shared(input_path, netCDF::NcFile::read); - auto var_set = nc_file->getVars(); - - // Populate the variable attributes cache - for (const auto& element : var_set) - { - std::string var_name = element.first; - auto ncvar = nc_file->getVar(var_name); - variable_names.push_back(var_name); - - std::string native_units; - auto units_att = ncvar.getAtt("units"); - units_att.getValues(native_units); - - double scale_factor = 1.0; - auto scale_var = ncvar.getAtt("scale_factor"); - if (!scale_var.isNull()) { - scale_var.getValues(&scale_factor); - } - - double offset = 0.0; - auto offset_var = ncvar.getAtt("add_offset"); - if (!offset_var.isNull()) { - offset_var.getValues(&offset); - } - - ncvar_cache[var_name] = {ncvar, native_units, scale_factor, offset}; - } - auto num_times = nc_file->getDim("time").getSize(); std::vector>> raw_time(num_times); @@ -59,7 +34,7 @@ NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_pat auto time_var = nc_file->getVar("Time"); if (time_var.getDimCount() != 1) { - throw "'Time' variable has dimension other than 1"; + throw std::runtime_error("'Time' variable has dimension other than 1"); } time_var.getVar(raw_time.data()); @@ -69,7 +44,7 @@ NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_pat time_unit_att.getValues(time_unit_str); if (time_unit_str != "minutes since 1970-01-01 00:00:00 UTC") { - throw "Time units not exactly as expected"; + throw std::runtime_error("Time units not exactly as expected"); } time_vals.reserve(num_times); @@ -82,24 +57,16 @@ NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_pat time_stride = std::chrono::duration_cast(time_vals[1] - time_vals[0]); -#if 0 // verify the time stride for( size_t i = 1; i < time_vals.size() -1; ++i) { auto tinterval = time_vals[i+1] - time_vals[i]; - if ( tinterval - time_stride > 1us) + if ( tinterval - time_stride > std::chrono::microseconds(1) ) { throw std::runtime_error("Time intervals in forcing file are not constant"); } } - - // determine start_time and stop_time; - start_time = epoch_start_time + time_vals[0]; - stop_time = epoch_start_time + time_vals.back() + time_stride; - - sim_to_data_time_offset = std::chrono::duration_cast(sim_start_date_time_epoch - start_time); -#endif } NetCDFMeshPointsDataProvider::~NetCDFMeshPointsDataProvider() = default; @@ -146,11 +113,12 @@ long NetCDFMeshPointsDataProvider::record_duration() const size_t NetCDFMeshPointsDataProvider::get_ts_index_for_time(const time_t &epoch_time_in) const { + // time_t in simulation engine's time frame - i.e. seconds, starting at 0 + auto epoch_time = sim_start_date_time_epoch + std::chrono::seconds(epoch_time_in); + auto start_time = time_vals.front(); auto stop_time = time_vals.back() + time_stride; - auto epoch_time = std::chrono::system_clock::from_time_t(epoch_time_in); - if (start_time <= epoch_time && epoch_time < stop_time) { auto offset = epoch_time - start_time; @@ -159,29 +127,43 @@ size_t NetCDFMeshPointsDataProvider::get_ts_index_for_time(const time_t &epoch_t } else { - //std::stringstream ss; - //ss << "The value " << (int)epoch_time << " was not in the range [" << (int)start_time << "," << (int)stop_time << ")\n" << SOURCE_LOC; - //throw std::out_of_range(ss.str().c_str()); - throw std::out_of_range(""); + std::stringstream ss; + ss << "The value " << (int)epoch_time_in << " was not in the range [" + << std::chrono::system_clock::to_time_t(start_time) << ", " + << std::chrono::system_clock::to_time_t(stop_time) << ")\n" + << SOURCE_LOC; + throw std::out_of_range(ss.str().c_str()); } } void NetCDFMeshPointsDataProvider::get_values(const selection_type& selector, boost::span data) { - if (!boost::get(&selector.points)) throw "Not implemented - only all_points"; + if (!boost::get(&selector.points)) throw std::runtime_error("Not implemented - only all_points"); + + cache_variable(selector.variable_name); auto const& metadata = ncvar_cache[selector.variable_name]; std::string const& source_units = metadata.units; size_t time_index = get_ts_index_for_time(std::chrono::system_clock::to_time_t(selector.init_time)); - metadata.variable.getVar({time_index, 0}, {1, data.size()}, data.data()); + metadata.ncVar.getVar({time_index, 0}, {1, data.size()}, data.data()); for (double& value : data) { value = value * metadata.scale_factor + metadata.offset; } - UnitsHelper::convert_values(source_units, data.data(), selector.output_units, data.data(), data.size()); + // These mass and and volume flux density units are very close to + // numerically identical for liquid water at atmospheric + // conditions, and so we currently treat them as interchangeable + bool RAINRATE_equivalence = + selector.variable_name == "RAINRATE" && + source_units == "mm s^-1" && + selector.output_units == "kg m-2 s-1"; + + if (!RAINRATE_equivalence) { + UnitsHelper::convert_values(source_units, data.data(), selector.output_units, data.data(), data.size()); + } } double NetCDFMeshPointsDataProvider::get_value(const selection_type& selector, ReSampleMethod m) @@ -191,6 +173,40 @@ double NetCDFMeshPointsDataProvider::get_value(const selection_type& selector, R return 0.0; } +void NetCDFMeshPointsDataProvider::cache_variable(std::string const& var_name) +{ + if (ncvar_cache.find(var_name) != ncvar_cache.end()) return; + + auto ncvar = nc_file->getVar(var_name); + variable_names.push_back(var_name); + + std::string native_units; + auto units_att = ncvar.getAtt("units"); + units_att.getValues(native_units); + + double scale_factor = 1.0; + try { + auto scale_var = ncvar.getAtt("scale_factor"); + if (!scale_var.isNull()) { + scale_var.getValues(&scale_factor); + } + } catch (...) { + // Assume it's just not present, and so keeps the default value + } + + double offset = 0.0; + try { + auto offset_var = ncvar.getAtt("add_offset"); + if (!offset_var.isNull()) { + offset_var.getValues(&offset); + } + } catch (...) { + // Assume it's just not present, and so keeps the default value + } + + ncvar_cache[var_name] = {ncvar, native_units, scale_factor, offset}; +} + } #endif From 1b29181878d3c9b5b039ea2f7bfb4c08b140dd3b Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Tue, 17 Sep 2024 17:02:05 -0700 Subject: [PATCH 48/51] Adjust Schism formulation and test to exploit new NetCDF mesh data provider --- .../coastal/SchismFormulation.hpp | 8 ++- .../coastal/SchismFormulation.cpp | 49 ++++++++----- test/coastal/SchismFormulation_Test.cpp | 72 +++++++++++++++++-- 3 files changed, 107 insertions(+), 22 deletions(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 9ed96e4735..22d63d6571 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -49,14 +49,18 @@ class SchismFormulation final : public CoastalFormulation std::unique_ptr bmi_; enum ForcingSelector { METEO, OFFSHORE, INFLOW }; - static std::map expected_input_variables_; + struct InputMapping { ForcingSelector selector; std::string name; }; + static std::map expected_input_variables_; std::map input_variable_units_; std::map input_variable_type_; std::map input_variable_count_; static std::vector exported_output_variable_names_; + std::map output_variable_units_; + std::map output_variable_type_; + std::map output_variable_count_; - std::chrono::time_point current_time_; + std::chrono::time_point current_time_; std::chrono::seconds time_step_length_; // TODO: Some of these maybe should be members of diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 2e745032b8..590aae08b1 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -7,26 +7,26 @@ const static auto s_schism_registration_function = "register_bmi"; -std::map SchismFormulation::expected_input_variables_ = +std::map SchismFormulation::expected_input_variables_ = { /* Meteorological Forcings */ // RAINRATE - precipitation - {"RAINRATE", SchismFormulation::METEO}, + {"RAINRATE", { SchismFormulation::METEO, "RAINRATE"}}, // SFCPRS - surface atmospheric pressure - {"SFCPRS", SchismFormulation::METEO}, + {"SFCPRS", { SchismFormulation::METEO, "PSFC"}}, // SPFH2m - specific humidity at 2m - {"SPFH2m", SchismFormulation::METEO}, + {"SPFH2m", { SchismFormulation::METEO, "Q2D"}}, // TMP2m - temperature at 2m - {"TMP2m", SchismFormulation::METEO}, + {"TMP2m", { SchismFormulation::METEO, "T2D"}}, // UU10m, VV10m - wind velocity components at 10m - {"UU10m", SchismFormulation::METEO}, - {"VV10m", SchismFormulation::METEO}, + {"UU10m", { SchismFormulation::METEO, "U2D"}}, + {"VV10m", { SchismFormulation::METEO, "V2D"}}, /* Input Boundary Conditions */ // ETA2_bnd - water surface elevation at the boundaries - {"ETA2_bnd", SchismFormulation::OFFSHORE}, + {"ETA2_bnd", { SchismFormulation::OFFSHORE, "ETA2_bnd"}}, // Q_bnd - flows at boundaries - {"Q_bnd", SchismFormulation::INFLOW}, + {"Q_bnd", { SchismFormulation::INFLOW, "Q_bnd"}}, }; std::vector SchismFormulation::exported_output_variable_names_ = @@ -57,7 +57,7 @@ SchismFormulation::SchismFormulation( , init_config_path , /* model_time_step_fixed = */ true , s_schism_registration_function - , MPI_COMM_WORLD + , MPI_COMM_SELF ); } @@ -77,6 +77,20 @@ void SchismFormulation::initialize() input_variable_count_[name] = mesh_size(name); } + auto const& output_vars = bmi_->GetOutputVarNames(); + + for (auto const& name : output_vars) { + //if (expected_output_variables_.find(name) == expected_output_variables_.end()) { + // throw std::runtime_error("SCHISM instance requests unexpected output variable '" + name + "'"); + //} + + output_variable_units_[name] = bmi_->GetVarUnits(name); + output_variable_type_[name] = bmi_->GetVarType(name); + output_variable_count_[name] = mesh_size(name); + } + + time_step_length_ = std::chrono::seconds((long long)bmi_->GetTimeStep()); + set_inputs(); } @@ -92,9 +106,11 @@ void SchismFormulation::finalize() void SchismFormulation::set_inputs() { for (auto var : expected_input_variables_) { - auto& name = var.first; - auto selector = var.second; - auto points = MeshPointsSelector{name, current_time_, time_step_length_, input_variable_units_[name], all_points}; + auto const& name = var.first; + auto const& mapping = var.second; + auto selector = mapping.selector; + auto const& source_name = mapping.name; + auto points = MeshPointsSelector{source_name, current_time_, time_step_length_, input_variable_units_[name], all_points}; ProviderType* provider = [this, selector](){ switch(selector) { @@ -104,13 +120,15 @@ void SchismFormulation::set_inputs() default: throw std::runtime_error("Unknown SCHISM provider selector type"); } }(); - std::vector values = provider->get_values(points); - bmi_->SetValue(name, values.data()); + std::vector buffer(mesh_size(name)); + provider->get_values(points, buffer); + bmi_->SetValue(name, buffer.data()); } } void SchismFormulation::update() { + current_time_ += time_step_length_; set_inputs(); bmi_->Update(); } @@ -166,5 +184,4 @@ size_t SchismFormulation::mesh_size(std::string const& variable_name) return nbytes / itemsize; } - #endif // NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI diff --git a/test/coastal/SchismFormulation_Test.cpp b/test/coastal/SchismFormulation_Test.cpp index 38caccd23e..e342b91047 100644 --- a/test/coastal/SchismFormulation_Test.cpp +++ b/test/coastal/SchismFormulation_Test.cpp @@ -2,9 +2,13 @@ #include #include "realizations/coastal/SchismFormulation.hpp" +#include +#include +#include const static std::string library_path = "/Users/phil/Code/noaa/BUILD/bmischism_2024-08-19-ngen/libtestbmifortranmodel.dylib"; const static std::string init_config_path = "/Users/phil/Code/noaa/SCHISM_Lake_Champlain_BMI_test/namelist.input"; +const static std::string met_forcing_netcdf_path = "/Users/phil/Code/noaa/NextGen_Forcings_Engine_SCHISM_Lake_Champlain_2015120100.nc"; #if 0 struct Schism_Formulation_IT : public ::testing::Test @@ -54,13 +58,13 @@ struct MockProvider : data_access::DataProvider size_t get_ts_index_for_time(const time_t &epoch_time) const override { return 1; } data_type get_value(const selection_type& selector, data_access::ReSampleMethod m) override { return data[0]; } - std::vector get_values(const selection_type& selector, data_access::ReSampleMethod) override + std::vector get_values(const selection_type& selector, data_access::ReSampleMethod m) override { throw ""; return data; } + void get_values(const selection_type& selector, boost::span out_data) override { auto default_value = input_variables_defaults[selector.variable_name]; - for (auto& val : data) { + for (auto& val : out_data) { val = default_value; } - return data; } }; @@ -70,17 +74,77 @@ int main(int argc, char **argv) MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); auto provider = std::make_shared(); + + std::tm start_time_tm{}; + start_time_tm.tm_year = 2015 - 1900; + start_time_tm.tm_mon = 12 - 1; + start_time_tm.tm_mday = 1; + auto start_time_t = std::mktime(&start_time_tm); + + std::tm stop_time_tm{}; + stop_time_tm.tm_year = 2015 - 1900; + stop_time_tm.tm_mon = 12 - 1; + stop_time_tm.tm_mday = 2; + auto stop_time_t = std::mktime(&stop_time_tm); + + auto netcdf_met_provider = std::make_shared(met_forcing_netcdf_path, + std::chrono::system_clock::from_time_t(start_time_t), + std::chrono::system_clock::from_time_t(stop_time_t)); auto schism = std::make_unique(/*id=*/ "test_schism_formulation", library_path, init_config_path, - provider, + netcdf_met_provider, provider, provider ); schism->initialize(); + for (int i = 0; i < 3; ++i) + schism->update(); + + using namespace std::chrono_literals; + + auto report = [](std::vector const& data, std::string name) { + double min = 10e6, max = -10e6; + for (int i = 0; i < data.size(); ++i) { + double val = data[i]; + if (std::isnan(val)) { + std::cout << "Nan found at " << i << std::endl; + break; + } + + min = std::min(val, min); + max = std::max(val, max); + } + std::cout << name << " ranges from " << min << " to " << max << std::endl; + }; + + std::vector bedlevel(278784, std::numeric_limits::quiet_NaN()); + MeshPointsSelector bedlevel_selector{"BEDLEVEL", std::chrono::system_clock::now(), 3600s, "m", all_points}; + schism->get_values(bedlevel_selector, bedlevel); + report(bedlevel, "BEDLEVEL"); + + std::vector eta2(278784, std::numeric_limits::quiet_NaN()); + MeshPointsSelector eta2_selector{"ETA2", std::chrono::system_clock::now(), 3600s, "m", all_points}; + schism->get_values(eta2_selector, eta2); + report(eta2, "ETA2"); + + std::vector vx(278784, std::numeric_limits::quiet_NaN()); + MeshPointsSelector vx_selector{"VX", std::chrono::system_clock::now(), 3600s, "m s-1", all_points}; + schism->get_values(vx_selector, vx); + report(vx, "VX"); + + std::vector vy(278784, std::numeric_limits::quiet_NaN()); + MeshPointsSelector vy_selector{"VY", std::chrono::system_clock::now(), 3600s, "m s-1", all_points}; + schism->get_values(vy_selector, vy); + report(vy, "VY"); + schism->update(); + schism->get_values(vx_selector, vx); + schism->get_values(vy_selector, vy); + report(vx, "VX"); + report(vy, "VY"); schism->finalize(); MPI_Finalize(); From 6a2a3dd35f23d4e35440b375128ec553d8b26415 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 18 Sep 2024 14:39:03 -0700 Subject: [PATCH 49/51] Time logic cleanup --- src/forcing/NetCDFMeshPointsDataProvider.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/forcing/NetCDFMeshPointsDataProvider.cpp b/src/forcing/NetCDFMeshPointsDataProvider.cpp index c7f4c88a21..61c8ad653f 100644 --- a/src/forcing/NetCDFMeshPointsDataProvider.cpp +++ b/src/forcing/NetCDFMeshPointsDataProvider.cpp @@ -28,17 +28,12 @@ NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_pat nc_file = std::make_shared(input_path, netCDF::NcFile::read); auto num_times = nc_file->getDim("time").getSize(); - - std::vector>> raw_time(num_times); - auto time_var = nc_file->getVar("Time"); if (time_var.getDimCount() != 1) { throw std::runtime_error("'Time' variable has dimension other than 1"); } - time_var.getVar(raw_time.data()); - auto time_unit_att = time_var.getAtt("units"); std::string time_unit_str; @@ -47,8 +42,10 @@ NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_pat throw std::runtime_error("Time units not exactly as expected"); } - time_vals.reserve(num_times); + std::vector>> raw_time(num_times); + time_var.getVar(raw_time.data()); + time_vals.reserve(num_times); for (int i = 0; i < num_times; ++i) { // Assume that the system clock's epoch matches Unix time. // This is guaranteed from C++20 onwards @@ -62,7 +59,8 @@ NetCDFMeshPointsDataProvider::NetCDFMeshPointsDataProvider(std::string input_pat { auto tinterval = time_vals[i+1] - time_vals[i]; - if ( tinterval - time_stride > std::chrono::microseconds(1) ) + if ( tinterval - time_stride > std::chrono::microseconds(1) || + time_stride - tinterval > std::chrono::microseconds(1) ) { throw std::runtime_error("Time intervals in forcing file are not constant"); } @@ -98,7 +96,7 @@ long NetCDFMeshPointsDataProvider::get_data_start_time() const long NetCDFMeshPointsDataProvider::get_data_stop_time() const { - return std::chrono::system_clock::to_time_t(time_vals.back()) + std::chrono::duration_cast(time_stride).count(); + return std::chrono::system_clock::to_time_t(time_vals.back() + time_stride); #if 0 //return stop_time; //FIXME: Matching behavior from CsvMeshPointsForcingProvider, but both are probably wrong! From b22f210d8fea0cfbadb763a2118f477de949bb66 Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Wed, 18 Sep 2024 16:14:55 -0700 Subject: [PATCH 50/51] SCHISM Formulation: Pass MPI communicator as constructor argument --- include/realizations/coastal/SchismFormulation.hpp | 1 + src/realizations/coastal/SchismFormulation.cpp | 3 ++- test/coastal/SchismFormulation_Test.cpp | 1 + 3 files changed, 4 insertions(+), 1 deletion(-) diff --git a/include/realizations/coastal/SchismFormulation.hpp b/include/realizations/coastal/SchismFormulation.hpp index 22d63d6571..04f1df3c15 100644 --- a/include/realizations/coastal/SchismFormulation.hpp +++ b/include/realizations/coastal/SchismFormulation.hpp @@ -18,6 +18,7 @@ class SchismFormulation final : public CoastalFormulation std::string const& id , std::string const& library_path , std::string const& init_config_path + , MPI_Comm mpi_comm , std::shared_ptr met_forcings , std::shared_ptr offshore_boundary , std::shared_ptr inflow_boundary diff --git a/src/realizations/coastal/SchismFormulation.cpp b/src/realizations/coastal/SchismFormulation.cpp index 590aae08b1..064ce728d2 100644 --- a/src/realizations/coastal/SchismFormulation.cpp +++ b/src/realizations/coastal/SchismFormulation.cpp @@ -41,6 +41,7 @@ SchismFormulation::SchismFormulation( std::string const& id , std::string const& library_path , std::string const& init_config_path + , MPI_Comm mpi_comm , std::shared_ptr met_forcings , std::shared_ptr offshore_boundary , std::shared_ptr inflow_boundary @@ -57,7 +58,7 @@ SchismFormulation::SchismFormulation( , init_config_path , /* model_time_step_fixed = */ true , s_schism_registration_function - , MPI_COMM_SELF + , mpi_comm ); } diff --git a/test/coastal/SchismFormulation_Test.cpp b/test/coastal/SchismFormulation_Test.cpp index e342b91047..921a285883 100644 --- a/test/coastal/SchismFormulation_Test.cpp +++ b/test/coastal/SchismFormulation_Test.cpp @@ -93,6 +93,7 @@ int main(int argc, char **argv) auto schism = std::make_unique(/*id=*/ "test_schism_formulation", library_path, init_config_path, + MPI_COMM_SELF, netcdf_met_provider, provider, provider From d680b07e960353bbe8a9e637805656ec11638dbb Mon Sep 17 00:00:00 2001 From: Phil Miller Date: Mon, 30 Sep 2024 14:41:56 -0700 Subject: [PATCH 51/51] Fix missing headers and conditional return in non-Python path --- src/core/parallel_utils.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/parallel_utils.cpp b/src/core/parallel_utils.cpp index 72ce754384..d0d0f87a5a 100644 --- a/src/core/parallel_utils.cpp +++ b/src/core/parallel_utils.cpp @@ -1,6 +1,8 @@ #include #include +#include +#include int mpi_rank = 0; @@ -294,11 +296,9 @@ namespace parallel { isGood = false; std::cerr << "Driver is unable to perform required hydrofabric subdividing when Python integration is not active." << std::endl; - // Sync with the rest of the ranks and bail if any aren't ready to proceed for any reason - if (!mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "initializing hydrofabric subdivider")) { - return false; - } + return mpiSyncStatusAnd(isGood, mpi_rank, mpi_num_procs, "initializing hydrofabric subdivider"); + #else // i.e., #if NGEN_WITH_PYTHON // Have rank 0 handle the generation task for all files/partitions std::unique_ptr subdivider;