Reset Time Signal (NGWPC-9621)

include/bmi_serialization.h

@@ -8,7 +8,7 @@ extern "C" {
 #include "bmi.h"
 
 const int serialize_topmodel(Bmi* bmi);
-const int deserialize_topmodel(Bmi* bmi, const char* buffer);
+const int deserialize_topmodel(Bmi* bmi, char* buffer);
 
 #ifdef __cplusplus
 }

include/topmodel.h

@@ -26,28 +26,28 @@
 
 /*** Function/subroutine prototypes ***/
 extern void convert_dist_to_histords(const double * const dist_from_outlet, const int num_channels,
-					const double * const chv, const double * const rv, const double dt, double* const tch);
+					const double chv, const double rv, const double dt, double* const tch);
 
 extern void calc_time_delay_histogram(int num_channels, double area,
-		double* tch, double *cum_dist_area_with_dist,
+		const double *tch, const double *cum_dist_area_with_dist,
 		int *num_time_delay_histo_ords, int *num_delay,	double **time_delay_histogram);
 
 extern void init_water_balance(
-				int num_topodex_values, double dt, double *sr0, 
-				double *szm, double *Q0, double *t0, double tl,
+				int num_topodex_values, double dt, double sr0, 
+				double szm, double Q0, double t0, double tl,
 				double **stor_unsat_zone, double *szq,
 				double **deficit_local, double **deficit_root_zone, 
 				double *sbar, double *bal);
 
-extern void init_discharge_array(int stand_alone, int *num_delay, double *Q0, double area, 
-			int *num_time_delay_histo_ords, double **time_delay_histogram,
+extern void init_discharge_array(int stand_alone, int num_delay, double Q0, double area, 
+			int num_time_delay_histo_ords, double **time_delay_histogram,
                         double **Q);
 
 extern int init(FILE *in_param_fptr, FILE *output_fptr, char *subcat, int stand_alone,
 	      int num_channels, int num_topodex_values, int yes_print_output,
 	      double area, double **time_delay_histogram,
 	      double *cum_dist_area_with_dist, double dt, 
-        double tl, double *dist_from_outlet, 
+        double tl, const double *dist_from_outlet, 
 	      int *num_time_delay_histo_ords,int *num_delay,
 	      double *szm, double *t0, double *chv, double *rv, double *td, double *srmax, 
         double *Q0,double *sr0, int *infex, double *xk0, double *hf, double *dth,
@@ -61,11 +61,11 @@ extern int inputs(FILE *input_fptr, int *nstep, double *dt, double **rain,
 extern void topmod(FILE *output_fptr, int nstep, int num_topodex_values,
                 int yes_print_output,int infex, double dt, double szm,
 	        double *stor_unsat_zone, double *deficit_root_zone,
-                double *deficit_local, double *pe, double *rain,double xk0,double hf, 
-                double *dist_area_lnaotb, double tl, double *lnaotb, double td,
+                double *deficit_local, const double *pe, const double *rain,double xk0,double hf, 
+                const double *dist_area_lnaotb, double tl, const double *lnaotb, double td,
                 double srmax, double *contrib_area, double szq, double *Qout, 
                 int num_time_delay_histo_ords,double *Q,
-                double *time_delay_histogram,char *subcat,double *bal,
+                const double *time_delay_histogram,
                 double *sbar,int num_delay, int current_time_step, int stand_alone,
                 double *sump, double *sumae, double *sumq, double *sumrz, double *sumuz,
                 double *quz, double *qb, double *qof, double *p, double *ep);

include/vecbuf.hpp

@@ -28,14 +28,17 @@ class vecbuf : public std::basic_streambuf<CharT> {
     // Forwarder for std::vector::clear()
     constexpr void clear() { vector_.clear(); }
 
+    // Forwarder for std::vector::resize(size)
+    constexpr void resize(size_type size) { vector_.resize(size); }
+
     // Forwarder for std::vector::reserve
     constexpr void reserve(size_type capacity) { vector_.reserve(capacity); setp_from_vector(); }
 
     // Increase the capacity of the buffer by reserving the current_size + additional_capacity
     constexpr void reserve_additional(size_type additional_capacity) { reserve(size() + additional_capacity); }
 
     // Forwarder for std::vector::data
-    constexpr const value_type* data() const { return vector_.data(); }
+    constexpr value_type* data() { return vector_.data(); }
 
     // Forwarder for std::vector::size
     constexpr size_type size() const { return vector_.size(); }
@@ -112,4 +115,11 @@ class vecbuf : public std::basic_streambuf<CharT> {
 
 };
 
+class membuf : public std::streambuf {
+public:
+    membuf(char *begin, size_t size) {
+        this->setg(begin, begin, begin + size);
+    }
+};
+
 #endif

src/bmi_serialization.cpp

@@ -32,6 +32,12 @@ class TopmodelSerializer {
 template<class Archive>
 void TopmodelSerializer::serialize(Archive& ar, const unsigned int version) {
     topmodel_model* model = this->model;
+    if (model->stand_alone == TRUE) {
+        // the number of timesteps makes hindcasting nigh imposible when stand alone
+        auto error = "Topmodel serialization is not currently implemented when running stand alone.";
+        Log(SEVERE, error);
+        throw std::runtime_error(error);
+    }
     ar & model->current_time_step;
 
     // data summed between runs
@@ -62,12 +68,37 @@ void TopmodelSerializer::serialize(Archive& ar, const unsigned int version) {
     ar & boost::serialization::make_array(
         model->deficit_local, num_topodex_values
     );
+
+    // nsteps will always be 1 for non-stand-alone models
     ar & boost::serialization::make_array(
-        model->contrib_area, num_topodex_values
+        model->contrib_area, model->nstep + 1
     );
+
+    // copy the current sizes to detect changes, then archive the model value
+    int num_time_delay_histo_ords = model->num_time_delay_histo_ords;
+    ar & model->num_time_delay_histo_ords;
+    int num_delay = model->num_delay;
+    ar & model->num_delay;
+    size_t num_Q = model->num_delay + model->num_time_delay_histo_ords + 1;
+    if (Archive::is_loading::value) {
+        // if loading and array size has changed, reallocate
+        if (num_time_delay_histo_ords != model->num_time_delay_histo_ords) {
+            if (model->time_delay_histogram != NULL)
+                free(model->time_delay_histogram);
+            model->time_delay_histogram = (double *)malloc(
+                (model->num_time_delay_histo_ords + 1) * sizeof(double)
+            );
+        }
+        if (num_delay != model->num_delay || num_time_delay_histo_ords != model->num_time_delay_histo_ords) {
+            if (model->Q != NULL)
+                free(model->Q);
+            model->Q = (double *)malloc(num_Q * sizeof(double));
+        }
+    }
     ar & boost::serialization::make_array(
-        model->Q, model->num_time_delay_histo_ords + 1
+        model->time_delay_histogram, model->num_time_delay_histo_ords + 1
     );
+    ar & boost::serialization::make_array(model->Q, num_Q);
 }
 
 
@@ -99,15 +130,17 @@ const int serialize_topmodel(Bmi* bmi) {
         free(model->serialized);
     }
     // set size and allocate memory
-    model->serialized_length = stream.size();
-    model->serialized = (char*)malloc(sizeof(char) * model->serialized_length);
+    uint64_t serialized_size = stream.size();
+    model->serialized_length = serialized_size + sizeof(uint64_t);
+    model->serialized = (char*)malloc(model->serialized_length);
     // make sure memory could be allocated
     if (model->serialized == NULL) {
         model->serialized_length = 0;
         return BMI_FAILURE;
     }
     // copy stream data to new allocation
-    memcpy(model->serialized, stream.data(), model->serialized_length);
+    memcpy(model->serialized, &serialized_size, sizeof(uint64_t));
+    memcpy(model->serialized + sizeof(uint64_t), stream.data(), serialized_size);
     return BMI_SUCCESS;
 }
 
@@ -118,9 +151,13 @@ const int serialize_topmodel(Bmi* bmi) {
   * @param buffer Start of data that wil be read as previously serialized state
   * @return int signifiying whether the serialization process completed successfully.
   */
-const int deserialize_topmodel(Bmi* bmi, const char* buffer) {
+const int deserialize_topmodel(Bmi* bmi, char* buffer) {
     TopmodelSerializer serializer(bmi);
-    std::istringstream stream(buffer);
+    // copy size of data out of header
+    uint64_t size;
+    memcpy(&size, buffer, sizeof(uint64_t));
+    // create stream from data after header
+    membuf stream(buffer + sizeof(uint64_t), size);
     boost::archive::binary_iarchive archive(stream);
     try {
         archive >> serializer;

src/bmi_topmodel.c

@@ -481,8 +481,6 @@ static int Update(Bmi *self) {
         topmodel->num_time_delay_histo_ords,
         topmodel->Q,
         topmodel->time_delay_histogram,
-        topmodel->subcat,
-        &topmodel->bal,
         &topmodel->sbar,
         topmodel->num_delay,
         topmodel->current_time_step,
@@ -652,6 +650,9 @@ static int Get_var_type(Bmi *self, const char *name, char *type) {
     } else if (strcmp(name, "serialization_free") == 0) {
         strncpy(type, "int", BMI_MAX_TYPE_NAME);
         return BMI_SUCCESS;
+    } else if (strcmp(name, "reset_time") == 0) {
+        strncpy(type, "double", BMI_MAX_TYPE_NAME);
+        return BMI_SUCCESS;
     }
     // If we get here, it means the variable name wasn't recognized
     type[0] = '\0';
@@ -776,7 +777,10 @@ static int Get_var_nbytes(Bmi *self, const char *name, int *nbytes) {
     }
     // special cases for save state
     if (item_count < 1) {
-        if (strcmp(name, "serialization_create") == 0 || strcmp(name, "serialization_size") == 0 || strcmp(name, "serialization_free") == 0) {
+        if (strcmp(name, "serialization_create") == 0
+            || strcmp(name, "serialization_size") == 0
+            || strcmp(name, "serialization_free") == 0
+            || strcmp(name, "reset_time") == 0) {
             item_count = 1;
         } else if (strcmp(name, "serialization_state") == 0) {
             topmodel_model* model = (topmodel_model*)self->data;
@@ -1063,6 +1067,11 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         } else {
             return BMI_FAILURE;
         }
+    } else if (strcmp(name, "reset_time") == 0) {
+        topmodel_model* model = (topmodel_model *)self->data;
+        // current_time_step is mainly used for indexing into config data, so should be safe to reset and nothing else
+        model->current_time_step = 0;
+        return BMI_SUCCESS;
     }
 
     if (self->get_value_ptr(self, name, &dest) == BMI_FAILURE)
@@ -1168,8 +1177,8 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         convert_dist_to_histords(
             topmodel->dist_from_outlet,
             topmodel->num_channels,
-            &topmodel->chv,
-            &topmodel->rv,
+            topmodel->chv,
+            topmodel->rv,
             topmodel->dt,
             tch
         );
@@ -1188,10 +1197,10 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         // Reinitialise discharge array
         init_discharge_array(
             topmodel->stand_alone,
-            &topmodel->num_delay,
-            &topmodel->Q0,
+            topmodel->num_delay,
+            topmodel->Q0,
             topmodel->area,
-            &topmodel->num_time_delay_histo_ords,
+            topmodel->num_time_delay_histo_ords,
             &topmodel->time_delay_histogram,
             &topmodel->Q
         );
@@ -1211,10 +1220,10 @@ static int Set_value(Bmi *self, const char *name, void *array) {
         init_water_balance(
             topmodel->num_topodex_values,
             topmodel->dt,
-            &topmodel->sr0,
-            &topmodel->szm,
-            &topmodel->Q0,
-            &topmodel->t0,
+            topmodel->sr0,
+            topmodel->szm,
+            topmodel->Q0,
+            topmodel->t0,
             topmodel->tl,
             &topmodel->stor_unsat_zone,
             &topmodel->szq,