Ptl add adaptive timestep (#25)

* adding the capability for an adaptive time step. If not provided in the config file, it defaults to false.

* updating readme

* updating readme

* updating giuh fxn with adaptive time step

* giuh now uses a fixed time step that is hardcoded to be the the same as the smallest adaptive time step for the rest of the model. If the adaptive time step is off, the giuh time step will be the same as the timestep input in the config file.

* updating readme in configs

* cleaned up code that was used for debugging

* updating Phillipsburg config file

* cleaned up code that was used for debugging

* After chatting with Ahmad, we determined that it is best if the adaptive time step strategy does not alter the function in giuh.c because this function is used in multiple models. So, the adaptive time step method has been altered such that the giuh is computed at the hourly time step rather than at the sub timestep. This allows for the giuh fxn to be unaltered, and has the added bonus of not reqiring resampling for the giuh based on time step.

* updating in response to comments, mostly focused on removing hardcoding from min and max time step values

* more changes in response to comments, including formatting of giuh.c, removing hardcoding for internal time steps and replacing with user specified minimum time step, and establishing a maximum time step equal to the forcing resolution

* removing code that was commented out

---------

Co-authored-by: Peter La Follette <[email protected]>
Co-authored-by: Peter La Follette <[email protected]>

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)

src/bmi_lgar.cxx

@@ -33,18 +33,19 @@ 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;
+  }
 
 }
 
@@ -76,7 +77,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 +110,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 +132,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++) {
@@ -358,16 +379,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,6 +429,13 @@ 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
 

src/lgar.cxx

@@ -213,6 +213,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
   // setting these options to false (defualt) 
   state->lgar_bmi_params.sft_coupled       = false;
   state->lgar_bmi_params.use_closed_form_G = false;
+  state->lgar_bmi_params.adaptive_timestep = false;
 
   bool is_layer_thickness_set       = false;
   bool is_initial_psi_set           = false;
@@ -398,6 +399,20 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
       continue;
     }
+    else if (param_key == "adaptive_timestep") { 
+      if ((param_value == "false") || (param_value == "0")) {
+        state->lgar_bmi_params.adaptive_timestep = false;
+      }
+      else if ( (param_value == "true") || (param_value == "1")) {
+        state->lgar_bmi_params.adaptive_timestep = true;
+      }
+      else {
+	std::cerr<<"Invalid option: adaptive_timestep must be true or false. \n";
+        abort();
+      }
+
+      continue;
+    }
     else if (param_key == "timestep") {
       state->lgar_bmi_params.timestep_h = stod(param_value);
 
@@ -411,6 +426,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
       assert (state->lgar_bmi_params.timestep_h > 0);
       is_timestep_set = true;
 
+      state->lgar_bmi_params.minimum_timestep_h = state->lgar_bmi_params.timestep_h;
+
       if (verbosity.compare("high") == 0) {
 	std::cerr<<"Model timestep [hours,seconds]: "<<state->lgar_bmi_params.timestep_h<<" , "
 		 <<state->lgar_bmi_params.timestep_h*3600<<"\n";
@@ -605,29 +622,28 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     throw runtime_error(errMsg.str());
   }
 
-
   if (is_giuh_ordinates_set) {
-    int factor = int(1.0/state->lgar_bmi_params.timestep_h);
+    int factor = int(1.0/state->lgar_bmi_params.forcing_resolution_h);
 
     state->lgar_bmi_params.num_giuh_ordinates = factor * (giuh_ordinates_temp.size() - 1);
     state->lgar_bmi_params.giuh_ordinates = new double[state->lgar_bmi_params.num_giuh_ordinates+1];
 
     for (int i=0; i<giuh_ordinates_temp.size()-1; i++) {
       for (int j=0; j<factor; j++) {
-	int index = j + i * factor + 1;
-	state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1]/double(factor);
+        int index = j + i * factor + 1;
+        state->lgar_bmi_params.giuh_ordinates[index] = giuh_ordinates_temp[i+1]/double(factor);
       }
     }
 
     if (verbosity.compare("high") == 0) {
       for (int i=1; i<=state->lgar_bmi_params.num_giuh_ordinates; i++)
-	std::cerr<<"GIUH ordinates (scaled) : "<<state->lgar_bmi_params.giuh_ordinates[i]<<"\n";
+	      std::cerr<<"GIUH ordinates (scaled) : "<<state->lgar_bmi_params.giuh_ordinates[i]<<"\n";
 
       std::cerr<<"          *****         \n";
     }
-
     giuh_ordinates_temp.clear();
   }
+
   else if (!is_giuh_ordinates_set) {
     stringstream errMsg;
     errMsg << "The configuration file \'" << config_file <<"\' does not set giuh_ordinates. \n";
@@ -1272,7 +1288,6 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
 	if (wf_free_drainage_demand == wf)
 	  prior_mass += precip_mass_to_add - (free_drainage_demand + actual_ET_demand);
   // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front
-
 	double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, AET_demand_cm,
 						   delta_thetas, delta_thickness, soil_type, soil_properties);
   actual_ET_demand = *AET_demand_cm;
@@ -2573,7 +2588,8 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
 
     // stop the loop if the error between the current and previous psi is less than 10^-15
     // 1. enough accuracy, 2. the algorithm can't improve the error further,
-    // 3. avoid infinite loop, 4. handles a corner case when prior mass is tiny (e.g., <1.E-5)
+    // 3. avoid infinite loop, 4. handles case where theta is very close to theta_r and convergence might be possible but would be extremely slow
+    // 5. handles a corner case when prior mass is tiny (e.g., <1.E-5)
     // printf("A1 = %.20f, %.18f %.18f %.18f %.18f \n ",fabs(psi_cm_loc - psi_cm_loc_prev) , psi_cm_loc, psi_cm_loc_prev, factor, delta_mass);
 
     if (fabs(psi_cm_loc - psi_cm_loc_prev) < 1E-15 && factor < 1E-13) break;
@@ -2588,6 +2604,11 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
     if (count_no_mass_change == break_no_mass_change)
       break;
 
+    if (psi_cm_loc > 1e7){//there are rare cases where theta is very close to theta_r, and delta_mass - delta_mass_prev will change extremely slowly. Convergence might be possible but the model will take hours to converge rather than seconds. 
+    //an alternative solution was to change the threshold in if (fabs(delta_mass - delta_mass_prev) < 1e-15) to 1e-11, but that solution is somewhat slow. 
+      break;
+    }
+
     // -ve pressure will return NAN, so terminate the loop if previous psi is way small and current psi is zero
     // the wetting front is almost saturated
     if (psi_cm_loc <= 0 && psi_cm_loc_prev < 0) break;

src/soil_funcs.cxx

@@ -154,7 +154,7 @@ double calc_theta_from_h(double h,double alpha, double m, double n, double theta
 /***********************************/
 double calc_Se_from_h(double h,double alpha, double m, double n)
 {
-  if(is_epsilon_less_than(h,1.0e-01)) return 1.0;  // this function doesn't work well ffor tiny h
+  if(is_epsilon_less_than(h,1.0e-10)) return 1.0;  // this function doesn't work well ffor tiny h
   else return(1.0/(pow(1.0+pow(alpha*h,n),m)));
 }