BACKUP COMMIT DO NOT MERGE

include/bmi/Bmi_Adapter.hpp

@@ -32,9 +32,10 @@ namespace models {
             {
                 // This replicates a lot of Initialize, but it's necessary to be able to do it separately to support
                 // "initializing" on construction, given using Initialize requires use of virtual functions
+                errno = 0;
                 if (!utils::FileChecker::file_is_readable(this->bmi_init_config)) {
                     init_exception_msg = "Cannot create and initialize " + this->model_name + " using unreadable file '"
-                            + this->bmi_init_config + "'";
+                            + this->bmi_init_config + "'. Error: "+std::strerror(errno);
                     throw std::runtime_error(init_exception_msg);
                 }
             }
@@ -125,6 +126,7 @@ namespace models {
             void Initialize() {
                 // If there was previous init attempt but w/ failure exception, throw runtime error and include previous
                 // message
+                errno = 0;
                 if (model_initialized && !init_exception_msg.empty()) {
                     throw std::runtime_error(
                             "Previous " + model_name + " init attempt had exception: \n\t" + init_exception_msg);
@@ -135,7 +137,7 @@ namespace models {
                 }
                 else if (!utils::FileChecker::file_is_readable(bmi_init_config)) {
                     init_exception_msg = "Cannot initialize " + model_name + " using unreadable file '"
-                            + bmi_init_config + "'";
+                            + bmi_init_config + "'. Error: "+std::strerror(errno);;
                     throw std::runtime_error(init_exception_msg);
                 }
                 else {

include/core/DomainLayer.hpp

@@ -0,0 +1,50 @@
+#ifndef __NGEN_DOMAIN_LAYER__
+#define __NGEN_DOMAIN_LAYER__
+
+#include "Catchment_Formulation.hpp"
+#include "Layer.hpp"
+
+namespace ngen
+{
+    class DomainLayer : public Layer
+    {
+        public:
+        // DomainLayer() = delete;
+        DomainLayer(
+                const LayerDescription& desc,
+                const Simulation_Time& s_t,
+                feature_type& features,
+                long idx,
+                std::shared_ptr<realization::Catchment_Formulation> formulation):
+                    //description(desc), features(features), output_time_index(0), simulation_time(s_t), formulation(formulation)
+                    Layer(desc, s_t, features, idx), formulation(formulation)
+        {
+            formulation->write_output("Time Step,""Time,"+formulation->get_output_header_line(",")+"\n");
+        }
+
+        /***
+         * @brief Run one simulation timestep for each model in this layer
+        */
+        void update_models() override{
+            formulation->get_response(output_time_index, simulation_time.get_output_interval_seconds());
+            std::string current_timestamp = simulation_time.get_timestamp(output_time_index);
+            std::string output = std::to_string(output_time_index)+","+current_timestamp+","+
+            formulation->get_output_line_for_timestep(output_time_index)+"\n";
+            formulation->write_output(output);
+            ++output_time_index;
+            if ( output_time_index < simulation_time.get_total_output_times() )
+            {
+               simulation_time.advance_timestep();
+            }
+
+        }
+
+        private:
+        // LayerDescription description;
+        // Simulation_Time simulation_time;
+        // long output_time_index;
+        std::shared_ptr<realization::Catchment_Formulation> formulation;
+    };
+}
+
+#endif

include/core/HY_Features.hpp

@@ -79,6 +79,7 @@ namespace hy_features {
          * @param id 
          * @return std::shared_ptr<HY_CatchmentRealization> 
          */
+        // std::shared_ptr<realization::Formulation> catchment_at(std::string id)
         std::shared_ptr<HY_CatchmentRealization> catchment_at(std::string id)
         {
           if( _catchments.find(id) != _catchments.end() )

include/core/Layer.hpp

@@ -40,6 +40,19 @@ namespace ngen
 
         }
 
+        Layer(
+                const LayerDescription& desc, 
+                const Simulation_Time& s_t, 
+                feature_type& f,
+                long idx) :
+            description(desc),
+            simulation_time(s_t),
+            features(f),
+            output_time_index(idx)
+        {
+
+        }
+
         virtual ~Layer() {}
 
         /***

include/core/LayerData.hpp

@@ -51,6 +51,17 @@ namespace ngen
 
         void put_layer(const LayerDescription& desc, int id) 
         {
+            // check to see if this layer was allready defined
+            if ( exists(id) )
+            {
+                std::string message = "A layer with id = ";
+                message += std::to_string(id);
+                message += " was defined more than once";
+
+                std::runtime_error r_error(message);
+
+                throw r_error;
+            }
              stored_layers[id] = desc;
              keys.push_back(id);
         }

include/core/catchment/HY_Catchment.hpp

@@ -7,6 +7,7 @@
 
 #include "HY_CatchmentRealization.hpp"
 #include "HY_HydroFeature.hpp"
+#include "Formulation.hpp"
 
 class HY_HydroNexus;
 
@@ -21,6 +22,7 @@ class HY_Catchment : public HY_HydroFeature
     public:
 
     HY_Catchment();
+    // HY_Catchment(std::string id, Nexuses inflows, Nexuses outflows, std::shared_ptr<realization::Formulation> realization, long lyr = 0):
     HY_Catchment(std::string id, Nexuses inflows, Nexuses outflows, std::shared_ptr<HY_CatchmentRealization> realization, long lyr = 0):
     id(std::move(id)),
     inflows(std::move(inflows)),
@@ -34,6 +36,7 @@ class HY_Catchment : public HY_HydroFeature
     layer(lyr){}
     virtual ~HY_Catchment();
     const Nexuses& get_outflow_nexuses(){ return outflows; }
+    // std::shared_ptr<realization::Formulation> realization;
     std::shared_ptr<HY_CatchmentRealization> realization;
 
     /*! \brief get the hydrofabric layer of this catchment

include/core/catchment/HY_CatchmentRealization.hpp

@@ -5,6 +5,7 @@
 #include <string>
 #include <AorcForcing.hpp>
 #include "GenericDataProvider.hpp"
+#include "Formulation.hpp"
 
 using std::shared_ptr;
 
@@ -15,6 +16,7 @@ class HY_Catchment;
 typedef long time_step_t;
 //TODO template<forcing>
 //TODO template<et_datatype>
+// class HY_CatchmentRealization: public realization::Formulation
 class HY_CatchmentRealization
 {
     public:

include/core/network.hpp

@@ -136,7 +136,48 @@ namespace network {
    * 
    */
   using IndexPair = std::pair< NetworkIndexT::const_iterator, NetworkIndexT::const_iterator>;
+struct detect_loops : public boost::dfs_visitor<>
+{
+    using colormap = std::map<Graph::vertex_descriptor, boost::default_color_type>;
+    colormap vertex_coloring;
+
+    using edgeColorMap = std::map<Graph::edge_descriptor, boost::default_color_type>;
+    edgeColorMap  edge_coloring;
+   using vertex_t = Graph::vertex_descriptor;
+
+   template <class Edge, class Graph>
+  void tree_edge(Edge e, const Graph& g) {
+      //std::cout << "tree_edge: " << boost::source(e, g) << " --> " << boost::target(e, g) << std::endl;
+
+    //edgeVisited.push(e);
+    if (vertexVisited.empty()) {
+        vertexVisited.push(boost::source(e, g));
+    }
+    vertexVisited.push(boost::target(e, g));
+  }
 
+   template <class Edge, class Graph>
+   void back_edge(Edge e, const Graph& g) {
+      std::cout << source(e, g)
+        << " -- "
+        << target(e, g) << "\n";
+      std::cout<< get(boost::vertex_name, g)[ source(e, g) ] <<
+        " -- " << get(boost::vertex_name, g)[ target(e,g)  ]<<"\n";
+       vertex_t v2;
+       std::cout << "Cycle end= " << boost::target(e, g) << std::endl;
+       while ( vertexVisited.top() != boost::target(e, g) )
+        {
+            //std::cout << " Cycle middle=" << vertexVisited.top() << std::endl;
+            v2 = vertexVisited.top();
+            vertexVisited.pop();
+        }
+        std::cout << "Cycle starting= " << vertexVisited.top() << std::endl;
+        vertexVisited.push(v2);
+      cycles.push_back(std::make_pair(source(e,g), target(e,g)));
+   }
+  std::stack<Graph::vertex_descriptor> vertexVisited;
+  std::vector<std::pair<Graph::vertex_descriptor, Graph::vertex_descriptor>> cycles;
+};
     /**
      * @brief A lightweight, graph based index of hydrologic features.
      * 

include/forcing/DataProviderSelectors.hpp

@@ -9,20 +9,15 @@
  * 
  */
 
-class CatchmentAggrDataSelector
-{
-    public:
+class DataSelector{
 
-    /**
-     * @brief Construct a new Catchment Aggregate Data Selector object with default values
-     * 
-     */
-    CatchmentAggrDataSelector() : 
+  public:
+
+    DataSelector():
         variable_name(), 
         init_time(0), 
         duration_s(1), 
-        output_units(),
-        id_str()
+        output_units()
     {}
 
     /**
@@ -34,8 +29,7 @@ class CatchmentAggrDataSelector
      * @param units The units to output the result in
      * @param id the id of the associated catchment
      */
-    CatchmentAggrDataSelector(std::string id, std::string var, time_t start, long dur, std::string units) : 
-        id_str(id),
+    DataSelector(std::string var, time_t start, long dur, std::string units) : 
         variable_name(var), 
         init_time(start), 
         duration_s(dur), 
@@ -98,6 +92,41 @@ class CatchmentAggrDataSelector
      */
     void set_output_units(std::string units) { output_units = units; }
 
+  private:
+
+    std::string variable_name; //!< The variable name that should be queried
+    time_t init_time; //!< The inital time to query the requested variable
+    long duration_s;  //!< The duration of the query
+    std::string output_units; //!< required units for the result to be return in
+}; 
+
+class CatchmentAggrDataSelector: public DataSelector
+{
+    public:
+
+    /**
+     * @brief Construct a new Catchment Aggregate Data Selector object with default values
+     * 
+     */
+    CatchmentAggrDataSelector() :
+        DataSelector(),
+        id_str()
+    {}
+
+    /**
+     * @brief Construct a new Data Selector object with inital values
+     * 
+     * @param var The variable to be accessed by the selector   
+     * @param start THe start time for this selector
+     * @param dur The duration for this selector   
+     * @param units The units to output the result in
+     * @param id the id of the associated catchment
+     */
+    CatchmentAggrDataSelector(std::string id, std::string var, time_t start, long dur, std::string units) : 
+        DataSelector(var, start, dur, units),
+        id_str(id)
+    {}
+
     /**
      * @brief Get the id string for this NetCDF Data Selector
      * 
@@ -114,10 +143,6 @@ class CatchmentAggrDataSelector
 
     private:
 
-    std::string variable_name; //!< The variable name that should be queried
-    time_t init_time; //!< The inital time to query the requested variable
-    long duration_s;  //!< The duration of the query
-    std::string output_units; //!< required units for the result to be return in
     std::string id_str; //< the catchment to access data for
 };