#737: SCHISM Coastal Formulation

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")

include/bmi/Bmi_Fortran_Adapter.hpp

@@ -11,6 +11,10 @@
 #include "State_Exception.hpp"
 #include "utilities/ExternalIntegrationException.hpp"
 
+#if NGEN_WITH_MPI
+#include <mpi.h>
+#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>(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;
 

include/forcing/DataProvider.hpp

@@ -96,6 +96,8 @@ namespace data_access
          */
         virtual std::vector<data_type> get_values(const selection_type& selector, ReSampleMethod m=SUM) = 0;
 
+        virtual void get_values(const selection_type& selector, boost::span<double> 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:

include/forcing/GenericDataProvider.hpp

@@ -3,10 +3,12 @@
 
 #include "DataProvider.hpp"
 #include "DataProviderSelectors.hpp"
+#include "MeshPointsSelectors.hpp"
 
 namespace data_access
 {
     using GenericDataProvider = DataProvider<double, CatchmentAggrDataSelector>;
+    using MeshPointsDataProvider = DataProvider<double, MeshPointsSelector>;
 }
 
 #endif

include/forcing/MeshPointsSelectors.hpp

@@ -0,0 +1,17 @@
+#pragma once
+
+#include <chrono>
+#include <string>
+#include <boost/variant.hpp>
+
+struct AllPoints {};
+static AllPoints all_points;
+
+struct MeshPointsSelector
+{
+    std::string variable_name;
+    std::chrono::time_point<std::chrono::system_clock> init_time;
+    std::chrono::seconds duration;
+    std::string output_units;
+    boost::variant<AllPoints, std::vector<int>> points;
+};

include/forcing/NetCDFMeshPointsDataProvider.hpp

@@ -0,0 +1,104 @@
+#pragma once
+
+#include <NGenConfig.h>
+
+#if NGEN_WITH_NETCDF
+
+#include "GenericDataProvider.hpp"
+#include "DataProviderSelectors.hpp"
+
+#include <string>
+#include <algorithm>
+#include <map>
+#include <memory>
+#include <string>
+#include <exception>
+
+namespace netCDF {
+    class NcVar;
+    class NcFile;
+}
+
+namespace data_access
+{
+    class NetCDFMeshPointsDataProvider : public MeshPointsDataProvider
+    {
+        public:
+
+        using time_point_type = std::chrono::time_point<std::chrono::system_clock>;
+
+        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 accessible by this data provider */
+        boost::span<const std::string> get_available_variable_names() const override;
+
+        /** 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 */
+        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.
+         */
+        data_type get_value(const selection_type& selector, ReSampleMethod m) override;
+
+        void get_values(const selection_type& selector, boost::span<data_type> data) override;
+
+        // And an implementation of the usual version using it
+        std::vector<data_type> get_values(const selection_type& selector, data_access::ReSampleMethod) override
+        {
+            throw std::runtime_error("Unimplemented");
+        }
+
+        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::vector<std::string> variable_names;
+        std::vector<time_point_type> time_vals;
+        std::chrono::seconds time_stride;                             // the amount of time between stored time values
+
+        std::shared_ptr<netCDF::NcFile> nc_file;
+
+        struct metadata_cache_entry;
+        std::map<std::string, metadata_cache_entry> ncvar_cache;
+    };
+}
+
+
+#endif // NGEN_WITH_NETCDF

include/realizations/coastal/CoastalFormulation.hpp

@@ -0,0 +1,50 @@
+#pragma once
+
+#include <forcing/GenericDataProvider.hpp>
+
+#include <string>
+#include <boost/core/span.hpp>
+
+class CoastalFormulation : public data_access::MeshPointsDataProvider
+{
+public:
+    CoastalFormulation(std::string id)
+        : id_(std::move(id))
+    { }
+
+    virtual ~CoastalFormulation() = default;
+
+    std::string get_id() const {
+        return this->id_;
+    }
+
+    virtual void initialize() = 0;
+    void finalize() override = 0;
+    virtual void update() = 0;
+
+    void get_values(const selection_type& selector, boost::span<double> data) override = 0;
+
+    // And an implementation of the usual version using it
+    std::vector<data_type> get_values(const selection_type& selector, data_access::ReSampleMethod) override
+    {
+        std::vector<data_type> 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 = boost::get<std::vector<int>>(&selector.points);
+        size_t size = points ? points->size() : this->mesh_size(selector.variable_name);
+        return size;
+    }
+
+private:
+    const std::string id_;
+};

include/realizations/coastal/SchismFormulation.hpp

@@ -0,0 +1,84 @@
+#pragma once
+
+#include <NGenConfig.h>
+
+#if NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI
+
+#include <realizations/coastal/CoastalFormulation.hpp>
+#include <bmi/Bmi_Fortran_Adapter.hpp>
+#include <memory>
+#include <map>
+#include <vector>
+
+class SchismFormulation final : public CoastalFormulation
+{
+public:
+    using ProviderType = data_access::MeshPointsDataProvider;
+
+    SchismFormulation(
+                      std::string const& id
+                      , std::string const& library_path
+                      , std::string const& init_config_path
+                      , MPI_Comm mpi_comm
+                      , std::shared_ptr<ProviderType> met_forcings
+                      , std::shared_ptr<ProviderType> offshore_boundary
+                      , std::shared_ptr<ProviderType> channel_flow_boundary
+                      );
+
+    ~SchismFormulation();
+
+    // Implementation of DataProvider
+    boost::span<const std::string> get_available_variable_names() const 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;
+
+    // Implementation of CoastalFormulation
+    void initialize() override;
+    void finalize() override;
+    void update() override;
+
+    void get_values(const selection_type& selector, boost::span<data_type> data) override;
+
+    // Visible only for testing use
+    enum ForcingSelector { METEO, OFFSHORE, CHANNEL_FLOW };
+    struct InputMapping { ForcingSelector selector; std::string name; };
+    static std::map<std::string, InputMapping> expected_input_variables_;
+
+protected:
+    size_t mesh_size(std::string const& variable_name) override;
+
+private:
+    std::unique_ptr<models::bmi::Bmi_Fortran_Adapter> bmi_;
+
+    std::map<std::string, std::string> input_variable_units_;
+    std::map<std::string, std::string> input_variable_type_;
+    std::map<std::string, size_t> input_variable_count_;
+
+    static std::vector<std::string> exported_output_variable_names_;
+    std::map<std::string, std::string> output_variable_units_;
+    std::map<std::string, std::string> output_variable_type_;
+    std::map<std::string, size_t> output_variable_count_;
+
+    std::chrono::time_point<std::chrono::system_clock> current_time_;
+    std::chrono::seconds time_step_length_;
+
+    // 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 make
+    // do with point values at the element centroids
+
+    std::shared_ptr<ProviderType> meteorological_forcings_provider_;
+    std::shared_ptr<ProviderType> offshore_boundary_provider_;
+    std::shared_ptr<ProviderType> channel_flow_boundary_provider_;
+
+    void set_inputs(int timestep_offset);
+};
+
+#endif // NGEN_WITH_BMI_FORTRAN && NGEN_WITH_MPI