Implementing a BMI mass balance checking protocol

.github/workflows/test_and_validate.yml

@@ -278,6 +278,35 @@ jobs:
       - name: Clean Up Unit Test Build
         uses: ./.github/actions/clean-build
 
+  # Run BMI protocol tests in linux/unix environment
+  test_bmi_protocols:
+    # The type of runner that the job will run on
+    strategy:
+      matrix:
+        os: [ubuntu-22.04, macos-12]
+      fail-fast: false
+    runs-on: ${{ matrix.os }}
+
+    # Steps represent a sequence of tasks that will be executed as part of the job
+    steps:
+      # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
+      - uses: actions/checkout@v4
+
+      - name: Build BMI Protocol Unit Tests
+        uses: ./.github/actions/ngen-build
+        with:
+          targets: "test_bmi_protocols"
+          build-cores: ${{ env.LINUX_NUM_PROC_CORES }}
+
+      - name: run_bmi_protocol_tests
+        run: |
+          cd ./cmake_build/test/
+          ./test_bmi_protocols
+          cd ../../
+        timeout-minutes: 15
+
+      - name: Clean Up BMI Protocol Unit Test Build
+        uses: ./.github/actions/clean-build
 
   # TODO: fails due to compilation error, at least in large part due to use of POSIX functions not supported on Windows.
   # TODO: Need to determine whether Windows support (in particular, development environment support) is necessary.

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;
 }

extern/test_bmi_cpp/include/test_bmi_cpp.hpp

@@ -25,6 +25,11 @@
 #define BMI_TYPE_NAME_SHORT "short"
 #define BMI_TYPE_NAME_LONG "long"
 
+#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"
+
 class TestBmiCpp : public bmi::Bmi {
     public:
         /**
@@ -179,6 +184,11 @@ class TestBmiCpp : public bmi::Bmi {
         std::vector<std::string> output_var_locations = { "node", "node" };
         std::vector<std::string> model_var_locations = {};
 
+        std::vector<std::string> mass_balance_var_names = { NGEN_MASS_IN, NGEN_MASS_OUT, NGEN_MASS_STORED, NGEN_MASS_LEAKED};
+        std::vector<std::string> mass_balance_var_types = { "double", "double", "double", "double"};
+        std::vector<std::string> mass_balance_var_units = { "m", "m", "m", "m" };
+        std::vector<std::string> mass_balance_var_locations = { "node", "node", "node", "node"};
+
         std::vector<int> input_var_item_count = { 1, 1 };
         std::vector<int> output_var_item_count = { 1, 1 };
         std::vector<int> model_var_item_count = {};
@@ -223,6 +233,9 @@ class TestBmiCpp : public bmi::Bmi {
         std::unique_ptr<double> model_var_1 = nullptr;
         std::unique_ptr<double> model_var_2 = nullptr;
 
+        double mass_stored = 0.0;
+        double mass_leaked = 0.0;
+
         /**
         * Read the BMI initialization config file and use its contents to set the state of the model.
         *

extern/test_bmi_cpp/src/test_bmi_cpp.cpp

@@ -165,6 +165,19 @@ void* TestBmiCpp::GetValuePtr(std::string name){
     }
   }
 
+  if (name == NGEN_MASS_STORED) {
+    return &this->mass_stored;
+  }
+  if (name == NGEN_MASS_LEAKED) {
+    return &this->mass_leaked;
+  }
+  if (name == NGEN_MASS_IN) {
+    return this->input_var_1.get();
+  }
+  if (name == NGEN_MASS_OUT) {
+    return this->output_var_1.get();
+  }
+
   throw std::runtime_error("GetValuePtr called for unknown variable: "+name);
 }
 
@@ -212,6 +225,10 @@ int TestBmiCpp::GetVarNbytes(std::string name){
   if(iter != this->model_var_names.end()){
     item_count = this->model_var_item_count[iter - this->model_var_names.begin()];
   }
+  iter = std::find(this->mass_balance_var_names.begin(), this->mass_balance_var_names.end(), name);
+  if(iter != this->mass_balance_var_names.end()){
+    item_count = 1;
+  }
   if(item_count == -1){
     // This is probably impossible to reach--the same conditions above failing will cause a throw
     // in GetVarItemSize --> GetVarType (called earlier) instead.
@@ -233,6 +250,10 @@ std::string TestBmiCpp::GetVarType(std::string name){
   if(iter != this->model_var_names.end()){
     return this->model_var_types[iter - this->model_var_names.begin()];
   }
+  iter = std::find(this->mass_balance_var_names.begin(), this->mass_balance_var_names.end(), name);
+  if(iter != this->mass_balance_var_names.end()){
+    return this->mass_balance_var_types[iter - this->mass_balance_var_names.begin()];
+  }
   throw std::runtime_error("GetVarType called for non-existent variable: "+name+"" SOURCE_LOC );
 }
 
@@ -249,6 +270,10 @@ std::string TestBmiCpp::GetVarUnits(std::string name){
   if(iter != this->model_var_names.end()){
     return this->model_var_types[iter - this->model_var_names.begin()];
   }
+  iter = std::find(this->mass_balance_var_names.begin(), this->mass_balance_var_names.end(), name);
+  if(iter != this->mass_balance_var_names.end()){
+    return this->mass_balance_var_units[iter - this->mass_balance_var_names.begin()];
+  }
   throw std::runtime_error("GetVarUnits called for non-existent variable: "+name+"" SOURCE_LOC);
 }
 
@@ -517,4 +542,6 @@ void TestBmiCpp::run(long dt)
         *this->output_var_5 = *this->model_var_2 * 1.0;
     }
     this->current_model_time += (double)dt;
+    this->mass_stored = *this->output_var_1 - *this->input_var_1;
+    this->mass_leaked = 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.run(models::bmi::protocols::Protocol::MASS_BALANCE, 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;