Memory Leak Fixes

configs/README.md

@@ -12,7 +12,7 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | initial_psi | double (scalar)| >=0 | cm | capillary head | - | used to initialize layers with a constant head |
 | ponded_depth_max | double (scalar)| >=0 | cm | maximum surface ponding | - | the maximum amount of water unavailable for surface drainage, default is set to zero |
 | timestep | double (scalar)| >0 | sec/min/hr | temporal resolution | - | timestep of the model |
-| forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data |
+| forcing_resolution | double (scalar)| - | sec/min/hr | temporal resolution | - | timestep of the forcing data. Recommended value of 3600 seconds. |
 | endtime | double (scalar)| >0 | sec, min, hr, d | simulation duration | - | time at which model simulation ends |
 | layer_soil_type | int (1D array) | - | - | state variable | - | layer soil type (read from the database file soil_params_file) |
 | max_soil_types | int | >1 | - | - | - | maximum number of soil types read from the file soil_params_file (default is set to 15) |
@@ -24,3 +24,4 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | sft_coupled | Boolean | true, false | - | model coupling | impacts hydraulic conductivity | couples LASAM to SFT. Coupling to SFT reduces hydraulic conducitivity, and hence infiltration, when soil is frozen|
 | soil_z | double (1D array) | - | cm | spatial resolution | - | vertical resolution of the soil column (computational domain of the SFT model) |
 | calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content |
+| adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep, and the above timestep is used as a minimum timestep (recommended value of 300 seconds). The adaptive timestep will never be larger than the forcing resolution. If set to false, LGAR will use the above specified timestep as a fixed timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. |

configs/config_lasam_Bushland.txt

@@ -13,3 +13,4 @@ max_soil_types=25
 wilting_point_psi=15495.0[cm]
 field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
+adaptive_timestep=true

configs/config_lasam_Phillipsburg.txt

@@ -13,3 +13,4 @@ max_soil_types=25
 wilting_point_psi=15495.0[cm]
 field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
+adaptive_timestep=true

configs/config_lasam_sft_ngen.txt

@@ -15,3 +15,4 @@ field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
 sft_coupled=true
 soil_z=10,20,30,40,50,60,70,80,90,100.0,110.,120,130.,140.,150.,160.,170.,180.,190.,200.0[cm]
+adaptive_timestep=true

include/all.hxx

@@ -114,6 +114,7 @@ struct lgar_bmi_parameters
   double initial_psi_cm;           // model initial (psi) condition
   double timestep_h;               // model timestep in hours
   double forcing_resolution_h;     // forcing resolution in hours
+  double minimum_timestep_h;       // minimum time step in hours, only used if adaptive_timestep is true
   int    forcing_interval;         // = forcing_resolution_h/timestep_h
   int    num_soil_types;           // number of soil types; must be less than or equal to MAX_NUM_SOIL_TYPES
   double AET_cm;                   // actual evapotranspiration in cm
@@ -131,6 +132,7 @@ struct lgar_bmi_parameters
   double  field_capacity_psi_cm;          // field capacity represented as a capillary head. Note that both wilting point and field capacity are specified for the whole model domain with single values
   bool   use_closed_form_G = false;      /* true if closed form of capillary drive calculation is desired, false if numeric integral
 					    for capillary drive calculation is desired */
+  bool   adaptive_timestep = false;      // if set to true, model uses adaptive timestep. In this case, the minimum timestep is the timestep specified in the config file. The maximum time step will be equal to the forcing resolution.
   double ponded_depth_cm;                // amount of water on the surface unavailable for surface runoff
   double ponded_depth_max_cm;            // maximum amount of water on the surface unavailable for surface runoff
   double precip_previous_timestep_cm;    // amount of rainfall (previous time step)
@@ -239,7 +241,7 @@ extern void                     listReverseOrder(struct wetting_front** head_ref
 extern bool                     listFindLayer(struct wetting_front* link, int num_layers, double *cum_layer_thickness_cm,
 					      int *lives_in_layer, bool *extends_to_bottom_flag);
 extern struct wetting_front*    listCopy(struct wetting_front* current, struct wetting_front* state_previous=NULL);
-
+extern void listDelete(struct wetting_front* head);
 
 
 

include/bmi_lgar.hxx

@@ -29,7 +29,8 @@ class NotImplemented : public std::logic_error {
 
 class BmiLGAR : public bmi::Bmi {
 public:
-  BmiLGAR() {
+  ~BmiLGAR();
+  BmiLGAR():giuh_ordinates(nullptr), giuh_runoff_queue(nullptr) {
     this->input_var_names[0] = "precipitation_rate";
     this->input_var_names[1] = "potential_evapotranspiration_rate";
     this->input_var_names[2] = "soil_temperature_profile";
@@ -131,6 +132,7 @@ public:
   struct model_state* get_model();
 
 private:
+  void realloc_soil();
   struct model_state* state;
   static const int input_var_name_count  = 3;
   static const int output_var_name_count = 15;

src/bmi_lgar.cxx

@@ -17,6 +17,14 @@
 // default verbosity is set to 'none' other option 'high' or 'low' needs to be specified in the config file
 string verbosity="none";
 
+/**
+ * @brief Delete dynamic arrays allocated in Initialize() and held by this object
+ * 
+ */
+BmiLGAR::~BmiLGAR(){
+  delete [] giuh_ordinates;
+  delete [] giuh_runoff_queue;
+}
 
 /* The `head` pointer stores the address in memory of the first member of the linked list containing
    all the wetting fronts. The contents of struct wetting_front are defined in "all.h" */
@@ -33,19 +41,33 @@ Initialize (std::string config_file)
 
   num_giuh_ordinates = state->lgar_bmi_params.num_giuh_ordinates;
 
-
   /* giuh ordinates are static and read in the lgar.cxx, and we need to have a copy of it to pass to
      giuh.cxx, so allocating/copying here*/
-  
+
   giuh_ordinates = new double[num_giuh_ordinates];
   giuh_runoff_queue = new double[num_giuh_ordinates+1];
 
-  for (int i=0; i<num_giuh_ordinates;i++)
+  for (int i=0; i<num_giuh_ordinates;i++){
     giuh_ordinates[i] = state->lgar_bmi_params.giuh_ordinates[i+1]; // note lgar uses 1-indexing
+  }
 
-  for (int i=0; i<=num_giuh_ordinates;i++)
+  for (int i=0; i<=num_giuh_ordinates;i++){
     giuh_runoff_queue[i] = 0.0;
+  }
+
+}
 
+/**
+ * @brief Allocate (or reallocate) storage for soil parameters
+ * 
+ */
+void BmiLGAR::realloc_soil(){
+
+  delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
+  delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
+
+  state->lgar_bmi_params.soil_depth_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts];
+  state->lgar_bmi_params.soil_moisture_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts];
 }
 
 /*
@@ -76,7 +98,7 @@ Update()
   double mm_to_cm = 0.1; // unit conversion
 
   // local variables for readibility
-  int subcycles = state->lgar_bmi_params.forcing_interval;
+  int subcycles;
   int num_layers = state->lgar_bmi_params.num_layers;
 
   // local variables for a full timestep (i.e., timestep of the forcing data)
@@ -109,14 +131,14 @@ Update()
   double volrech_subtimestep_cm;
   double surface_runoff_subtimestep_cm; // direct surface runoff
   double precip_previous_subtimestep_cm;
-  double volrunoff_giuh_subtimestep_cm;
   double volQ_gw_subtimestep_cm = 0.0; // fix it for non-zero values after adding groundwater reservoir
 
   double subtimestep_h = state->lgar_bmi_params.timestep_h;
   int nint = state->lgar_bmi_params.nint;
   double wilting_point_psi_cm = state->lgar_bmi_params.wilting_point_psi_cm;
   double field_capacity_psi_cm = state->lgar_bmi_params.field_capacity_psi_cm;
   bool use_closed_form_G = state->lgar_bmi_params.use_closed_form_G; 
+  bool adaptive_timestep = state->lgar_bmi_params.adaptive_timestep;
 
   // constant value used in the AET function
   double AET_thresh_Theta = 0.85;    // scaled soil moisture (0-1) above which AET=PET (fix later!)
@@ -131,6 +153,26 @@ Update()
 
   assert (state->lgar_bmi_input_params->precipitation_mm_per_h >= 0.0);
   assert(state->lgar_bmi_input_params->PET_mm_per_h >=0.0);
+
+  // adaptive time step is set 
+  if (adaptive_timestep) {
+    subtimestep_h = state->lgar_bmi_params.forcing_resolution_h;
+    if (state->lgar_bmi_input_params->precipitation_mm_per_h > 10.0 || volon_timestep_cm > 0.0 ) {
+      subtimestep_h = state->lgar_bmi_params.minimum_timestep_h;  //case where precip > 1 cm/h, or there is ponded head from the last time step
+    }
+    else if (state->lgar_bmi_input_params->precipitation_mm_per_h > 0.0) {
+      subtimestep_h = state->lgar_bmi_params.minimum_timestep_h * 2.0;  //case where precip is less than 1 cm/h but greater than 0, and there is no ponded head 
+    }
+    subtimestep_h = fmin(subtimestep_h, state->lgar_bmi_params.forcing_resolution_h);  //just in case the user has specified a minimum time step that would make the subtimestep_h greater than the forcing resolution 
+    state->lgar_bmi_params.timestep_h = subtimestep_h;
+  }
+
+  state->lgar_bmi_params.forcing_interval = int(state->lgar_bmi_params.forcing_resolution_h/state->lgar_bmi_params.timestep_h+1.0e-08); // add 1.0e-08 to prevent truncation error
+  subcycles = state->lgar_bmi_params.forcing_interval;
+
+  if (verbosity.compare("high") == 0) {
+    printf("time step size in hours: %lf \n", state->lgar_bmi_params.timestep_h);
+  }
 
   // subcycling loop (loop over model's timestep)
   for (int cycle=1; cycle <= subcycles; cycle++) {
@@ -143,7 +185,10 @@ Update()
       std::cerr<<"BMI Update |Timesteps = "<< state->lgar_bmi_params.timesteps<<", Time [h] = "<<this->state->lgar_bmi_params.time_s / 3600.<<", Subcycle = "<< cycle <<" of "<<subcycles<<std::endl;
     }
 
-    state->state_previous = NULL;
+    if( state->state_previous != NULL ){
+      listDelete(state->state_previous);
+      state->state_previous = NULL;
+    }
     state->state_previous = listCopy(state->head);
 
     // ensure precip and PET are non-negative
@@ -270,7 +315,10 @@ Update()
         listPrint(state->head);
       }
 
-      state->state_previous = NULL;
+      if(state->state_previous != NULL ){
+        listDelete(state->state_previous);
+        state->state_previous = NULL;
+      }
       state->state_previous = listCopy(state->head);
 
       volin_timestep_cm += volin_subtimestep_cm;
@@ -358,16 +406,9 @@ Update()
 
 
     /*----------------------------------------------------------------------*/
-    // compute giuh runoff for the subtimestep
+    // increment runoff for the subtimestep
     surface_runoff_subtimestep_cm = volrunoff_subtimestep_cm;
-    volrunoff_giuh_subtimestep_cm = giuh_convolution_integral(volrunoff_subtimestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
-
     surface_runoff_timestep_cm += surface_runoff_subtimestep_cm ;
-    volrunoff_giuh_timestep_cm += volrunoff_giuh_subtimestep_cm;
-
-    // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well.
-
-    volQ_timestep_cm += volrunoff_giuh_subtimestep_cm;
 
     // adding groundwater flux to stream channel (note: this will be updated/corrected after adding the groundwater reservoir)
     volQ_gw_timestep_cm += volQ_gw_subtimestep_cm;
@@ -415,15 +456,21 @@ Update()
 
   } // end of subcycling
 
+  //update giuh at the time step level (was previously updated at the sub time step level)
+  volrunoff_giuh_timestep_cm = giuh_convolution_integral(volrunoff_timestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue);
+
+  // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well.
+  // when groundwater component is added, it should probably happen inside of the subcycling loop.
+  volQ_timestep_cm = volrunoff_giuh_timestep_cm;
+
   /*----------------------------------------------------------------------*/
   // Everything related to lgar state is done at this point, now time to update some dynamic variables
 
   // update number of wetting fronts
   state->lgar_bmi_params.num_wetting_fronts = listLength(state->head);
 
   // allocate new memory based on updated wetting fronts; we could make it conditional i.e. create only if no. of wf are changed
-  state->lgar_bmi_params.soil_depth_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts];
-  state->lgar_bmi_params.soil_moisture_wetting_fronts = new double[state->lgar_bmi_params.num_wetting_fronts];
+  realloc_soil();
 
   // update thickness/depth and soil moisture of wetting fronts (used for state coupling)
   struct wetting_front *current = state->head;
@@ -583,6 +630,30 @@ void BmiLGAR::
 Finalize()
 {
   global_mass_balance();
+  listDelete(state->head);
+  listDelete(state->state_previous);
+
+  delete [] state->soil_properties;
+
+  delete [] state->lgar_bmi_params.soil_depth_wetting_fronts;
+  delete [] state->lgar_bmi_params.soil_moisture_wetting_fronts;
+
+  delete [] state->lgar_bmi_params.soil_temperature;
+  delete [] state->lgar_bmi_params.soil_temperature_z;
+  delete [] state->lgar_bmi_params.layer_soil_type;
+
+  delete [] state->lgar_calib_params.theta_e;
+  delete [] state->lgar_calib_params.theta_r;
+  delete [] state->lgar_calib_params.vg_n;
+  delete [] state->lgar_calib_params.vg_alpha;
+  delete [] state->lgar_calib_params.Ksat;
+
+  delete [] state->lgar_bmi_params.layer_thickness_cm;
+  delete [] state->lgar_bmi_params.cum_layer_thickness_cm;
+  delete [] state->lgar_bmi_params.giuh_ordinates;
+  delete [] state->lgar_bmi_params.frozen_factor;
+  delete state->lgar_bmi_input_params;
+  delete state;
 }
 
 

src/bmi_main_lgar.cxx

@@ -155,13 +155,14 @@ int main(int argc, char *argv[])
       // write layers data to file
       fprintf(outlayer_fptr,"# Timestep = %d, %s \n", i, time[i].c_str());
       write_state(outlayer_fptr, model_state.get_model()->head);
+      delete [] soil_moisture_wetting_front;
+      delete [] soil_thickness_wetting_front;
     }
 
   }
 
-  // do final mass balance
-  model_state.global_mass_balance();
-
+  // do final mass balance ( inside Finalize() ) and finish the simulation
+  model_state.Finalize();
   if (outdata_fptr) {
     fclose(outdata_fptr);
     fclose(outlayer_fptr);