Implementing a BMI mass balance checking protocol

CMakeLists.txt

@@ -321,6 +321,7 @@ add_subdirectory("src/forcing")
 add_subdirectory("src/utilities/mdarray")
 add_subdirectory("src/utilities/mdframe")
 add_subdirectory("src/utilities/logging")
+add_subdirectory("src/utilities/bmi")
 
 target_link_libraries(ngen
     PUBLIC
@@ -333,6 +334,7 @@ target_link_libraries(ngen
         NGen::forcing
         NGen::core_mediator
         NGen::logging
+        NGen::bmi_protocols
 )
 
 if(NGEN_WITH_SQLITE)

extern/test_bmi_c/include/bmi_test_bmi_c.h

@@ -1,6 +1,11 @@
 #ifndef BMI_TEST_BMI_C_H
 #define BMI_TEST_BMI_C_H
 
+#define NGEN_MASS_IN "ngen::mass_in"
+#define NGEN_MASS_OUT "ngen::mass_out"
+#define NGEN_MASS_STORED "ngen::mass_stored"
+#define NGEN_MASS_LEAKED "ngen::mass_leaked"
+
 #if defined(__cplusplus)
 extern "C" {
 #endif

extern/test_bmi_c/include/test_bmi_c.h

@@ -31,6 +31,9 @@ struct test_bmi_c_model {
     int param_var_1;
     double param_var_2;
     double* param_var_3;
+
+    double mass_stored; // Mass balance variable, for testing purposes
+    double mass_leaked; //Mass balance variable, for testing purposes
 };
 typedef struct test_bmi_c_model test_bmi_c_model;
 

extern/test_bmi_c/src/bmi_test_bmi_c.c

@@ -9,6 +9,7 @@
 #define INPUT_VAR_NAME_COUNT 2
 #define OUTPUT_VAR_NAME_COUNT 2
 #define PARAM_VAR_NAME_COUNT 3
+#define MASS_BALANCE_VAR_NAME_COUNT 4
 
 // Don't forget to update Get_value/Get_value_at_indices (and setter) implementation if these are adjusted
 static const char *output_var_names[OUTPUT_VAR_NAME_COUNT] = { "OUTPUT_VAR_1", "OUTPUT_VAR_2" };
@@ -34,6 +35,13 @@ static const int param_var_item_count[PARAM_VAR_NAME_COUNT] = { 1, 1, 2 };
 static const char *param_var_grids[PARAM_VAR_NAME_COUNT] = { 0, 0, 0 };
 static const char *param_var_locations[PARAM_VAR_NAME_COUNT] = { "node", "node", "node" };
 
+static const char *mass_balance_var_names[MASS_BALANCE_VAR_NAME_COUNT] = { NGEN_MASS_IN, NGEN_MASS_OUT, NGEN_MASS_STORED, NGEN_MASS_LEAKED};
+static const char *mass_balance_var_types[MASS_BALANCE_VAR_NAME_COUNT] = { "double", "double", "double", "double"};
+static const char *mass_balance_var_units[MASS_BALANCE_VAR_NAME_COUNT] = { "m", "m", "m", "m" };
+static const int mass_balance_var_item_count[MASS_BALANCE_VAR_NAME_COUNT] = { 1, 1, 1, 1};
+static const char *mass_balance_var_grids[MASS_BALANCE_VAR_NAME_COUNT] = { 0, 0, 0, 0 };
+static const char *mass_balance_var_locations[MASS_BALANCE_VAR_NAME_COUNT] = { "node", "node", "node", "node" };
+
 static int Finalize (Bmi *self)
 {
     // Function assumes everything that is needed is retrieved from the model before Finalize is called.
@@ -387,6 +395,23 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest)
         *dest = ((test_bmi_c_model *)(self->data))->param_var_3;
         return BMI_SUCCESS;
     }
+
+    if (strcmp (name, NGEN_MASS_IN) == 0) {
+        *dest = ((test_bmi_c_model *)(self->data))->input_var_1;
+        return BMI_SUCCESS;
+    }
+    if (strcmp (name, NGEN_MASS_OUT) == 0) {
+        *dest = ((test_bmi_c_model *)(self->data))->output_var_1;
+        return BMI_SUCCESS;
+    }
+    if (strcmp (name, NGEN_MASS_STORED) == 0) {
+        *dest = &((test_bmi_c_model *)(self->data))->mass_stored;
+        return BMI_SUCCESS;
+    }
+    if (strcmp (name, NGEN_MASS_LEAKED) == 0) {
+        *dest = &((test_bmi_c_model *)(self->data))->mass_leaked;
+        return BMI_SUCCESS;
+    }
     return BMI_FAILURE;
 }
 
@@ -483,6 +508,14 @@ static int Get_var_nbytes (Bmi *self, const char *name, int * nbytes)
             }
         }
     }
+    if (item_count < 1) {
+        for (i = 0; i < MASS_BALANCE_VAR_NAME_COUNT; i++) {
+            if (strcmp(name, mass_balance_var_names[i]) == 0) {
+                item_count = mass_balance_var_item_count[i];
+                break;
+            }
+        }
+    }
     if (item_count < 1)
         item_count = ((test_bmi_c_model *) self->data)->num_time_steps;
 
@@ -515,6 +548,13 @@ static int Get_var_type (Bmi *self, const char *name, char * type)
             return BMI_SUCCESS;
         }
     }
+    // Finally check to see if in mass balance array
+    for (i = 0; i < MASS_BALANCE_VAR_NAME_COUNT; i++) {
+        if (strcmp(name, mass_balance_var_names[i]) == 0) {
+            snprintf(type, BMI_MAX_TYPE_NAME, "%s", mass_balance_var_types[i]);
+            return BMI_SUCCESS;
+        }
+    }
     // If we get here, it means the variable name wasn't recognized
     type[0] = '\0';
     return BMI_FAILURE;
@@ -538,6 +578,13 @@ static int Get_var_units (Bmi *self, const char *name, char * units)
             return BMI_SUCCESS;
         }
     }
+    //Check for mass balance
+    for (i = 0; i < MASS_BALANCE_VAR_NAME_COUNT; i++) {
+        if (strcmp(name, mass_balance_var_names[i]) == 0) {
+            snprintf(units, BMI_MAX_UNITS_NAME, "%s", mass_balance_var_units[i]);
+            return BMI_SUCCESS;
+        }
+    }
     // If we get here, it means the variable name wasn't recognized
     units[0] = '\0';
     return BMI_FAILURE;
@@ -560,6 +607,9 @@ static int Initialize (Bmi *self, const char *file)
     else
         model = (test_bmi_c_model *) self->data;
 
+    model->mass_stored = 0.0;
+    model->mass_leaked = 0.0;
+
     if (read_init_config(file, model) == BMI_FAILURE)
         return BMI_FAILURE;
 

extern/test_bmi_c/src/test_bmi_c.c

@@ -19,5 +19,7 @@ extern int run(test_bmi_c_model* model, long dt)
     }
     model->current_model_time += (double)dt;
 
+    model->mass_stored = *model->output_var_1 - *model->input_var_1;
+    model->mass_leaked = 0;
     return 0;
 }

include/core/Layer.hpp

@@ -119,6 +119,8 @@ namespace ngen
                 double response(0.0);
                 try{
                     response = r_c->get_response(output_time_index, simulation_time.get_output_interval_seconds());
+                    // Check mass balance if able
+                    r_c->check_mass_balance(output_time_index, simulation_time.get_total_output_times(), current_timestamp);
                 }
                 catch(models::external::State_Exception& e){
                     std::string msg = e.what();
@@ -127,6 +129,13 @@ namespace ngen
                              +" at feature id "+id;
                     throw models::external::State_Exception(msg);
                 }
+                catch(std::exception& e){
+                    std::string msg = e.what();
+                    msg = msg+" at timestep "+std::to_string(output_time_index)
+                             +" ("+current_timestamp+")"
+                             +" at feature id "+id;
+                    throw std::runtime_error(msg);
+                }
                 std::string output = std::to_string(output_time_index)+","+current_timestamp+","+
                                     r_c->get_output_line_for_timestep(output_time_index)+"\n";
                 r_c->write_output(output);

include/realizations/catchment/Bmi_Module_Formulation.hpp

@@ -7,6 +7,7 @@
 #include "Bmi_Adapter.hpp"
 #include <DataProvider.hpp>
 #include "bmi_utilities.hpp"
+#include "bmi/protocols.hpp"
 
 using data_access::MEAN;
 using data_access::SUM;
@@ -250,6 +251,12 @@ namespace realization {
         const std::vector<std::string> get_bmi_input_variables() const override;
         const std::vector<std::string> get_bmi_output_variables() const override;
 
+        virtual void check_mass_balance(const int& iteration, const int& total_steps, const std::string& timestamp) const override {
+            //Create the protocol context, each member is const, and cannot change during the check
+            models::bmi::protocols::Context ctx{iteration, total_steps, timestamp, id};
+            bmi_protocols.mass_balance.run(ctx);
+        }
+
     protected:
 
         /**
@@ -419,6 +426,7 @@ namespace realization {
         int next_time_step_index = 0;
 
     private:
+        models::bmi::protocols::NgenBmiProtocols bmi_protocols;
         /**
          * Whether model ``Update`` calls are allowed and handled in some way by the backing model for time steps after
          * the model's ``end_time``.

include/realizations/catchment/Bmi_Multi_Formulation.hpp

@@ -42,6 +42,15 @@ namespace realization {
 
         virtual ~Bmi_Multi_Formulation() {};
 
+        virtual void check_mass_balance(const int& iteration, const int& total_steps, const std::string& timestamp) const final {
+            for( const auto &module : modules ) {
+                // TODO may need to check on outputs form each module indepdently???
+                // Right now, the assumption is that if each component is mass balanced
+                // then the entire formulation is mass balanced
+                module->check_mass_balance(iteration, total_steps, timestamp);
+            }
+        };
+
         /**
          * Convert a time value from the model to an epoch time in seconds.
          *

include/realizations/catchment/Formulation.hpp

@@ -44,6 +44,7 @@ namespace realization {
             virtual void create_formulation(boost::property_tree::ptree &config, geojson::PropertyMap *global = nullptr) = 0;
             virtual void create_formulation(geojson::PropertyMap properties) = 0;
 
+            virtual void check_mass_balance(const int& iteration, const int& total_steps, const std::string& timestamp) const = 0;
         protected:
 
             virtual const std::vector<std::string>& get_required_parameters() const = 0;