Add standalone restart support

CMakeLists.txt

@@ -20,21 +20,21 @@ message(CMAKE_C_COMPILER " ${CMAKE_C_COMPILER}")
 message("CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 
 # standalone
-add_executable(lasam_standalone ./src/bmi_main_lgar.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx
+add_executable(lasam_standalone ./src/bmi_main_lgar.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx ./src/init_from_state.cxx
              ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx
              ./giuh/giuh.h ./giuh/giuh.c)
 target_link_libraries(lasam_standalone PRIVATE m)
 
 # unittest
-add_executable(lasam_unitest ./tests/main_unit_test_bmi.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx
+add_executable(lasam_unitest ./tests/main_unit_test_bmi.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx ./src/conceptual_reservoir.cxx ./src/init_from_state.cxx
              ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx ./giuh/giuh.h
              ./giuh/giuh.c)
 target_link_libraries(lasam_unitest PRIVATE m)
 
 # Make sure these are compiled with this directive
 add_compile_definitions(BMI_ACTIVE)
 
-add_library(lasambmi SHARED src/bmi_lgar.cxx src/lgar.cxx ./src/soil_funcs.cxx ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/conceptual_reservoir.cxx
+add_library(lasambmi SHARED src/bmi_lgar.cxx src/lgar.cxx ./src/soil_funcs.cxx ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/conceptual_reservoir.cxx ./src/init_from_state.cxx
         ./src/util_funcs.cxx ./src/aet.cxx ./giuh/giuh.c include/all.hxx ./giuh/giuh.h)
 
 target_compile_definitions(lasambmi PRIVATE NGEN)

configs/README.md

@@ -35,6 +35,51 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | spf_factor | double (scalar) | 0.1 <= spf_factor <= 1 | - | parameter for fluxes to nonlinear reservoir | storage that contributes directly to streamflow | Simple bypass of surface water to the nonlinear reservoir will occur when the most superficial wetting front achieves the theta_e value of its layer times spf_factor. When this occurs, the amount of water sent to the nonlinear reservoir is equal to the precipitation plus any ponded water times frac_to_CR. This is a rather simple representation of preferential flow that intends to simulate the episodic nature of streamflow events in arid or semi arid environments. Defaults to 0.98. |
 | allow_flux_caching | Boolean | true, false | - | trades a small amount of accuracy for a lot of speed | flux caching | During dry periods, it is often the case that wetting fronts will move very slowly and AET will be significantly less than PET. In these cases, in the context of streamflow simulation, it is not efficient to recompute fluxes and soil moisture dynamics for each time step. If this is set to true, then fluxes and wetting front movement will only be recomputed once every 24 hours, or when the conditions resulting in dry and slow wetting fronts and low AET cease. During the times for which fluxes are not recomputed, instead they are stored in a cache and fluxes for subsequent time steps are set using this cache. Sligtly different strategies are used for fluxes through the lower boundary and AET. Also note that because NextGen models should ideally provide output for each hour, simply setting an adaptive time step to be larger than one hour is not a preferred runtime reduction method here. Note that this can cause small mass balance errors when the lower boundary condition is set to free drainage. Defaults to false. |
 | log_mode | Boolean | true, false | - | helps calibration search space exploration | log transform of parameters | When this is set to true, then all inputs for the van Genuchten parameter alpha, saturated hydraulic conductivity, and the nonlinear reservoir parameter a must be input as their log values rather than the normal values. For example, if an saturated hydraulic conductivity of 0.1 cm/h is desired, then the input value must be -1 because 10^-1 = 0.1. The reasoning for this is that these parameters are not distributed normally in nature but rather are distributed log normally, such that simply sampling the parameter space normally during calibration will vastly undersample a big region of the parameter space in which we expect useful parameter sets to be. Defaults to false. |
+| init_state_path | string | - | - | filename | optional initialization | If set, CASAM will initialize by loading the wetting fronts from the last data row in this file. The expected format is the standalone-mode `data_layers.csv` restart file. |
+| init_non_vadose_state_path | string | - | - | filename | optional initialization | If set, CASAM will initialize non-vadose restart variables from the last data row in this file, including conceptual reservoir storage, `volon_timestep_cm`, `runoff_in_prev_step`, `precip_previous_timestep_cm`, and flux-caching state. The expected format is the standalone-mode `data_non_vadose_state.csv` restart file. If `allow_flux_caching=false`, saved flux-cache fields are ignored. |
+| init_giuh_state_path | string | - | - | filename | optional initialization | If set, CASAM will initialize the saved GIUH runoff queue from the last data row in this file. The expected format is the standalone-mode `data_giuh_state.csv` restart file. |
 | a_slow | double (scalar) | 1E-8 < a_slow < 1E-1 | cm^(1-b) h^-1 | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter a, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0.|
 | b_slow | double (scalar) | 0.01 < b_slow < 5 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter b, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0.|
-| frac_slow | double (scalar) | 0.0 < frac_slow <= 1 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This describes the partitioning of water to the two reservoris, where the the input to the slow reservoir is equal to the total input for the nonlinear reservoirs times frac_slow. Note that either all or none of a_slow, b_slow, and frac_slow must be specified. If none are specified then the model will not simulate a second nonlinear reservoir. Defaults to 0.|
+| frac_slow | double (scalar) | 0.0 < frac_slow <= 1 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This describes the partitioning of water to the two reservoris, where the the input to the slow reservoir is equal to the total input for the nonlinear reservoirs times frac_slow. Note that either all or none of a_slow, b_slow, and frac_slow must be specified. If none are specified then the model will not simulate a second nonlinear reservoir. Defaults to 0.|
+
+### Restart Files
+
+Restart files are currently written only by standalone mode. Restart initialization reads the last data row from each requested restart file. If an invalid soil type is specified, restart files are ignored because CASAM returns input precipitation as output discharge for invalid soil types.
+
+The standalone `data_layers.csv` restart file saves and loads the following wetting front variables:
+
+| Variable | Units | Description |
+| -------- | ----- | ----------- |
+| depth_cm | cm | wetting front depth |
+| theta | - | volumetric water content |
+| layer_num | - | 1-based soil layer number |
+| front_num | - | 1-based wetting front number |
+| to_bottom | - | whether the wetting front is in contact with a layer bottom |
+| psi_cm | cm | capillary pressure head |
+| dzdt_cm_per_h | cm/h | wetting front velocity |
+
+The standalone `data_non_vadose_state.csv` restart file saves and loads the following non-vadose and conceptual reservoir variables:
+
+| Variable | Units | Description |
+| -------- | ----- | ----------- |
+| CR_fast_storage_cm | cm | fast conceptual reservoir storage |
+| CR_slow_storage_cm | cm | slow conceptual reservoir storage |
+| volon_timestep_cm | cm | surface water available at the previous timestep |
+| runoff_in_prev_step | - | whether runoff occurred in the previous timestep |
+| precip_previous_timestep_cm | cm | precipitation from the previous timestep |
+| cache_fluxes | - | whether cached fluxes were active |
+| cache_count | - | number of timesteps using cached fluxes |
+| previous_AET | cm | cached actual evapotranspiration |
+| previous_PET | cm | cached potential evapotranspiration |
+| previous_recharge | cm | cached recharge |
+| accumulated_PET | cm | accumulated PET while using cached fluxes |
+| accumulated_free_drainage | cm | accumulated free drainage while using cached fluxes |
+
+If `allow_flux_caching=false`, saved values for `cache_fluxes`, `cache_count`, `previous_AET`, `previous_PET`, `previous_recharge`, `accumulated_PET`, and `accumulated_free_drainage` are ignored and initialized as inactive cache state.
+
+The standalone `data_giuh_state.csv` restart file saves and loads the following GIUH variables:
+
+| Variable | Units | Description |
+| -------- | ----- | ----------- |
+| num_giuh_ordinates | - | number of GIUH ordinates expected by the saved queue |
+| queue | cm | GIUH runoff queue values from index 0 through `num_giuh_ordinates` |

configs/config_lasam_Bushland.txt

@@ -1,13 +1,13 @@
-verbosity=none
+verbosity=high
 forcing_file=./forcing/forcing_data_resampled_uniform_Bushland.csv
 soil_params_file=./data/vG_default_params.dat
 layer_thickness=18.0,76.0,135.0[cm]
 initial_psi=2000.0[cm]
-timestep=3600[sec]
+timestep=300[sec]
 endtime=8333.0[hr]
 forcing_resolution=3600[sec]
 ponded_depth_max=0[cm]
-use_closed_form_G=true
+use_closed_form_G=false
 layer_soil_type=16,17,18
 max_valid_soil_types=25
 wilting_point_psi=15495.0[cm]
@@ -18,4 +18,7 @@ a=0.001
 b=3.0
 frac_to_CR=0.1
 spf_factor=0.9
-allow_flux_caching=true
+allow_flux_caching=true
+#init_state_path=./data_layers.csv
+#init_non_vadose_state_path=./data_non_vadose_state.csv
+#init_giuh_state_path=./data_giuh_state.csv

configs/config_lasam_Phillipsburg_with_reservoir.txt

@@ -19,3 +19,6 @@ a=0.001
 b=2.5
 spf_factor=0.9
 frac_to_CR=0.2
+#init_state_path=./data_layers.csv
+#init_non_vadose_state_path=./data_non_vadose_state.csv
+#init_giuh_state_path=./data_giuh_state.csv

include/all.hxx

@@ -132,6 +132,9 @@ 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 */
+  string init_state_path;                // optional path for saved wetting front state
+  string init_non_vadose_state_path;      // optional path for conceptual reservoir / restart state
+  string init_giuh_state_path;            // optional path for saved GIUH queue state
   bool   PET_affects_precip = false;     // set to true in config file if you want PET to be taken from precip 
   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.
   bool   free_drainage_enabled = false;  // free_drainage_enabled will specify whether the lower boundary condition is no flow (false), or free drainage (true). Defaults to false.
@@ -403,6 +406,21 @@ extern void InitFromConfigFile(string config_file, struct model_state *state);
 extern vector<double> ReadVectorData(string key);
 extern void InitializeWettingFronts(bool is_invalid_soil_type, int num_layers, double initial_psi_cm, int *layer_soil_type, double *cum_layer_thickness_cm,
 				    double *frozen_factor, struct wetting_front** head, struct soil_properties_ *soil_properties);
+extern void InitializeWettingFrontsFromCSV(
+    int num_layers,
+    const char *data_layers_csv_path,
+    int *layer_soil_type,
+    double *cum_layer_thickness_cm,
+    double *frozen_factor,
+    struct wetting_front **head,
+    struct soil_properties_ *soil_properties);
+extern void InitializenonvadoseStateFromCSV(
+    const char *non_vadose_state_csv_path,
+    struct model_state *state);
+extern void InitializeGIUHRunoffQueueFromCSV(
+    const char *giuh_state_csv_path,
+    double *giuh_runoff_queue,
+    int num_giuh_ordinates);
 
 /********************************************************************/
 /*Other function prototypes for doing hydrology calculations, etc.  */
@@ -423,6 +441,11 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run
 
 // writes full state of wetting fronts (depth, theta, no. of wetting front, no. of layer, dz/dt, psi) to a file at each time step
 extern void write_state(FILE *out, struct wetting_front* head);
+extern void write_non_vadose_state(FILE *out, struct model_state *state);
+extern void write_giuh_runoff_queue_state(
+    FILE *out,
+    const double *giuh_runoff_queue,
+    int num_giuh_ordinates);
 
 
 /********************************************************************/

include/bmi_lgar.hxx

@@ -147,6 +147,8 @@ public:
   void global_mass_balance();
   double update_calibratable_parameters();
   struct model_state* get_model();
+  double* get_giuh_runoff_queue();
+  int get_num_giuh_ordinates();
 
 private:
   void realloc_soil();

src/bmi_lgar.cxx

@@ -9,6 +9,7 @@
 #include <vector>
 #include <algorithm>
 #include <iostream>
+#include <iomanip>
 #include "../bmi/bmi.hxx"
 #include "../include/bmi_lgar.hxx"
 #include "../include/all.hxx"
@@ -106,8 +107,35 @@ Initialize (std::string config_file)
     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++){
-    giuh_runoff_queue[i] = 0.0;
+  if (!state->lgar_bmi_params.is_invalid_soil_type &&
+      !state->lgar_bmi_params.init_giuh_state_path.empty()) {
+    InitializeGIUHRunoffQueueFromCSV(
+        state->lgar_bmi_params.init_giuh_state_path.c_str(),
+        giuh_runoff_queue,
+        num_giuh_ordinates);
+
+    if (verbosity.compare("high") == 0) {
+      streamsize old_precision = std::cerr.precision();
+      ios::fmtflags old_flags = std::cerr.flags();
+
+      std::cerr << "GIUH runoff queue initialized from saved state file: "
+                << state->lgar_bmi_params.init_giuh_state_path << "\n";
+      std::cerr << "GIUH runoff queue values: [";
+      std::cerr << std::fixed << std::setprecision(12);
+      for (int i = 0; i <= num_giuh_ordinates; i++) {
+        std::cerr << giuh_runoff_queue[i];
+        if (i < num_giuh_ordinates) std::cerr << ",";
+      }
+      std::cerr << "]\n";
+      std::cerr.flags(old_flags);
+      std::cerr.precision(old_precision);
+      std::cerr << "          *****         \n";
+    }
+  }
+  else {
+    for (int i=0; i<=num_giuh_ordinates;i++){
+      giuh_runoff_queue[i] = 0.0;
+    }
   }
 
 }
@@ -269,6 +297,10 @@ Update()
     state->lgar_bmi_params.timestep_h = subtimestep_h;
   }
 
+  if (!state->lgar_bmi_params.allow_flux_caching) {
+    state->lgar_mass_balance.cache_fluxes = FALSE;
+  }
+
   bool caching_at_start = state->lgar_mass_balance.cache_fluxes;
   bool switch_caching = FALSE;
 
@@ -1663,4 +1695,14 @@ GetGridNodesPerFace(const int grid, int *nodes_per_face)
   throw bmi_lgar::NotImplemented();
 }
 
+double* BmiLGAR::get_giuh_runoff_queue()
+{
+  return giuh_runoff_queue;
+}
+
+int BmiLGAR::get_num_giuh_ordinates()
+{
+  return num_giuh_ordinates;
+}
+
 #endif