partitioning subset hydrofabrics (NGWPC-6553)

CMakeLists.txt

@@ -403,6 +403,7 @@ message(STATUS "----------------------------------------------------------------
 ngen_multiline_message(
 "NGen version: ${ngen_VERSION}"
 "Build configuration summary:"
+"  Source directory: ${CMAKE_SOURCE_DIR}"
 "  Generator: ${CMAKE_GENERATOR}"
 "  Build type: ${CMAKE_BUILD_TYPE}"
 "  System: ${CMAKE_SYSTEM_NAME}"

include/geojson/FeatureVisitor.hpp

@@ -9,6 +9,7 @@ namespace geojson {
     class MultiLineStringFeature;
     class MultiPolygonFeature;
     class CollectionFeature;
+    class SentinelFeature;
 
     /**
      * An abstract class used to operate on the children of FeatureBase
@@ -27,6 +28,7 @@ namespace geojson {
             virtual void visit(MultiLineStringFeature *feature) = 0;
             virtual void visit(MultiPolygonFeature *feature) = 0;
             virtual void visit(CollectionFeature* feature) = 0;
+            virtual void visit(SentinelFeature *feature) = 0;
 
             virtual ~FeatureVisitor() = default;
     };

include/geojson/features/FeatureBase.hpp

@@ -41,16 +41,21 @@ namespace geojson {
 
     /**
      *  Describes a type of features
+     *
+     *  These are numbered in accordance with the GeoPackage
+     *  'well-known binary' (WKB) feature types, other than the
+     *  additional 'Sentinel' value defined here
      */
     enum class FeatureType {
-        None,                   /*!< Represents an empty feature with no sort of geometry */
-        Point,                  /*!< Represents a feature that contains a single Point geometry */
-        LineString,             /*!< Represents a feature that is represented by a series of interconnected points */
-        Polygon,                /*!< Represents a feature that is represented by a defined area */
-        MultiPoint,             /*!< Represents a feature that is represented by many points */
-        MultiLineString,        /*!< Represents a feature that is represented by multiple series of interconnected points */
-        MultiPolygon,           /*!< Represents a feature that is represented by multiple areas */
-        GeometryCollection      /*!< Represents a feature that contains a collection of different types of geometry */
+        None                 = 0,  /*!< Represents an empty feature with no sort of geometry */
+        Point                = 1,  /*!< Represents a feature that contains a single Point geometry */
+        LineString           = 2,  /*!< Represents a feature that is represented by a series of interconnected points */
+        Polygon              = 3,  /*!< Represents a feature that is represented by a defined area */
+        MultiPoint           = 4,  /*!< Represents a feature that is represented by many points */
+        MultiLineString      = 5,  /*!< Represents a feature that is represented by multiple series of interconnected points */
+        MultiPolygon         = 6,  /*!< Represents a feature that is represented by multiple areas */
+        GeometryCollection   = 7,  /*!< Represents a feature that contains a collection of different types of geometry */
+        Sentinel             = 100 /*!< Represents a 'dummy' feature included for computational consistency */
     };
 
     /**

include/geojson/features/Features.hpp

@@ -9,5 +9,6 @@
 #include "MultiLineStringFeature.hpp"
 #include "MultiPolygonFeature.hpp"
 #include "CollectionFeature.hpp"
+#include "SentinelFeature.hpp"
 
-#endif // GEOJSON_FEATURE_H
+#endif // GEOJSON_FEATURE_H

include/geojson/features/SentinelFeature.hpp

@@ -0,0 +1,29 @@
+#ifndef GEOJSON_SENTINEL_FEATURE_H
+#define GEOJSON_SENTINEL_FEATURE_H
+
+#include "FeatureBase.hpp"
+#include <FeatureVisitor.hpp>
+#include <JSONGeometry.hpp>
+
+#include <string>
+#include <vector>
+#include <map>
+#include <exception>
+
+namespace geojson {
+
+    class SentinelFeature : public FeatureBase {
+        public:
+            SentinelFeature(
+                std::string new_id
+            ) : FeatureBase(std::move(new_id)) {
+                this->type = geojson::FeatureType::Sentinel;
+            }
+
+            void visit(FeatureVisitor& visitor) override {
+                visitor.visit(this);
+            }
+    };
+}
+
+#endif // GEOJSON_SENTINEL_FEATURE_H

src/NGen.cpp

@@ -444,6 +444,13 @@ int main(int argc, char* argv[]) {
     // TODO: Instead of iterating through a collection of FeatureBase objects mapping to catchments,
     // we instead want to iterate through HY_Catchment objects
     geojson::GeoJSON catchment_collection;
+    // As part of the fix for NOAA-OWP/ngen#284 / NGWPC-6553,
+    // partitioning may insert sentinel flowpaths downstream of
+    // terminal nexuses. Those sentinels will not exist in the
+    // catchmentDataFile. Their listing in catchment_subset_ids works
+    // because the respective geoFOO::read() functions return the
+    // intersection of features in the file and the specified subset,
+    // rather than erroring on missing features.
     if (boost::algorithm::ends_with(catchmentDataFile, "gpkg")) {
 #if NGEN_WITH_SQLITE3
         try {

src/partitionGenerator.cpp

@@ -107,15 +107,16 @@ void write_remote_connections(const PartitionVSet& catchment_part, const Partiti
  * 
  * @param network 
  * @param num_partitions 
- * @param num_catchments
- * @param catchment_part 
+ * @param catchment_part  - which catchments will be simulation on each partition
+ * @param nexus_part - which nexuses have contributing catchments on each partition
  */
-void generate_partitions(network::Network& network, const int& num_partitions, const int& num_catchments, PartitionVSet& catchment_part,
-     PartitionVSet& nexus_part)
+void generate_partitions(network::Network& network, const int& num_partitions, PartitionVSet& catchment_part, PartitionVSet& nexus_part)
 {
+    auto catchments = network.filter("cat", network::SortOrder::TransposedDepthFirstPreorder);
+
     int partition = 0;
     int counter = 0;
-    int total = num_catchments;
+    int total = size(catchments);
     int partition_size = total/num_partitions;
     int partition_size_norm = partition_size;
     int remainder;
@@ -134,9 +135,7 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
     nexus_set.reserve(partition_size);
     std::string part_id, partition_str;
 
-    std::string up_nexus;
-    std::string down_nexus;
-    for(const auto& catchment : network.filter("cat", network::SortOrder::TransposedDepthFirstPreorder)){
+    for(const auto& catchment : catchments){
             if (partition < remainder)
                 partition_size = partition_size_plus1;
             else
@@ -150,9 +149,7 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
             }
             if(nexus_set.size() == 0){
                 partgen_ss <<"Error: Catchment "<<catchment<<" has no destination nexus.\n";
-                LOG(partgen_ss.str(), LogLevel::FATAL); partgen_ss.str("");
-
-                exit(1);
+                LOG(partgen_ss.str(), LogLevel::WARNING); partgen_ss.str("");
             }
             for( auto upstream : network.get_origination_ids(catchment) ){
                 nexus_set.emplace(upstream);
@@ -161,18 +158,12 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
 
             //keep track of all the features in this partition
             catchment_set.emplace(catchment);
+            LOG(catchment + " placed in partition " + std::to_string(partition), LogLevel::DEBUG);
+
             counter++;
             if(counter == partition_size)
             {
-                //partgen_ss<<"nexus "<<nexus<<" is remote DOWN on partition "<<partition<<std::endl;
                 //FIXME partitioning shouldn't have to assume dendritic network
-                std::vector<std::string> destinations = network.get_destination_ids(catchment);
-                if(destinations.size() == 0){
-                    partgen_ss <<"Error: Catchment "<<catchment<<" has no destination nexus.\n";
-                    LOG(partgen_ss.str(), LogLevel::FATAL); partgen_ss.str("");
-                    exit(1);
-                }
-                down_nexus = destinations[0];
 
                 part_id = std::to_string(partition);  // Is id used?
                 partition_str = std::to_string(partition);
@@ -187,12 +178,6 @@ void generate_partitions(network::Network& network, const int& num_partitions, c
 
                 partition++;
                 counter = 0;
-                //partgen_ss<<"\nnexus "<<nexus<<" is remote UP on partition "<<partition<<std::endl;
-
-                //this nexus overlaps partitions
-                //Handled above by ensure all up/down stream nexuses are recorded
-                up_nexus = down_nexus;
-                //partgen_ss<<"\nin partition "<<partition<<":"<<std::endl;
             }
     }
 
@@ -474,11 +459,8 @@ int main(int argc, char* argv[])
 
     std::string link_key = "toid";
 
-    Network catchment_network(catchment_collection, &link_key);
     //Assumes dendritic, can add check in network if needed.
     PartitionVSet catchment_part, nexus_part;
-
-    //catchment_network.print_network();
 
     //build the remote connections from network
     // read the nexus hydrofabric, reuse the catchments
@@ -517,17 +499,64 @@ int main(int argc, char* argv[])
     //Do this before linking features so that the alt ids can lookup the correct feature
     global_nexus_collection->update_ids("id");
     global_nexus_collection->link_features_from_property(nullptr, &link_key);
+
+    // ngen misbehaves in gathering and storing output when a terminal
+    // nexus is fed by multiple catchments partitioned to different
+    // processes. This was recorded as NOAA-OWP/ngen#284 and NGWPC-6553.
+    //
+    // Address that here by inserting sentinel flowpaths downstream of
+    // those nexuses. Those sentinels will be assigned to specific
+    // processes, which will then properly receive all of the flow for
+    // the nexus in question.
+    //
+    // These features will not exist in the hydrofabric
+    // GeoJSON/GeoPackage. As implemented, that means that they will
+    // simply be ignored when ngen figures out what features to
+    // simulate on each process.
+    //
+    // Store the sentinels separately to avoid iterator invalidation
+    // from inserting them eagerly
+    std::vector<std::shared_ptr<geojson::FeatureBase>> sentinels;
+    for (auto& feature : *global_nexus_collection)
+    {
+        auto id = feature->get_id();
+        auto type = id.substr(0,3);
+        if (hy_features::identifiers::isNexus(type)) {
+            if (feature->get_number_of_destination_features() == 0) {
+                std::string sentinel_id = "wb-TERMINAL_SENTINEL-" + feature->get_id();
+                geojson::Feature sentinel_feature = std::make_shared<geojson::SentinelFeature>(sentinel_id);
+                sentinels.push_back(sentinel_feature);
+                feature->add_destination_feature(sentinel_feature.get());
+                LOG("Nexus " + feature->get_id() + " has no destination features; adding " + sentinel_id + " below it", LogLevel::INFO);
+            }
+        }
+    }
+    for (auto& sentinel : sentinels)
+    {
+        global_nexus_collection->add_feature(sentinel);
+    }
+
     // make a global network
     Network global_network(global_nexus_collection);
 
     //Generate the partitioning
-    generate_partitions(global_network, num_partitions, num_catchments, catchment_part, nexus_part);
+    generate_partitions(global_network, num_partitions, catchment_part, nexus_part);
 
     //global_network.print_network();
 
     //The container holding all remote_connections
     std::vector<RemoteConnectionVec> remote_connections_vec;
 
+    for (int i = 0; i < num_partitions; ++i) {
+        partgen_ss << "Partition " << i << " catchments: " << catchment_part[i].size() << "\n";
+        for (auto& c : catchment_part[i])
+            partgen_ss << c << std::endl;
+        partgen_ss << "nexuses " << nexus_part[i].size() << "\n";
+        for (auto& n : nexus_part[i])
+            partgen_ss << n << std::endl;
+        LOG(partgen_ss.str(), LogLevel::DEBUG); partgen_ss.str("");
+    }
+
     int total_remotes = 0;
     // loop over all partitions by partition id
     for (int ipart=0; ipart < catchment_part.size(); ++ipart)

test/geojson/FeatureCollection_Test.cpp

@@ -55,6 +55,10 @@ class Visitor : public geojson::FeatureVisitor {
             this->types.push_back("CollectionFeature");
         }
 
+        void visit(geojson::SentinelFeature *feature) override {
+            this->types.push_back("SentinelFeature");
+        }
+
         std::string get(int index) {
             if( index >= 0 && index < types.size() )
               return this->types[index];