v491 - Support for workflow vars with op regimes; :AllowSoilOverfill; recursive :RedirectToFile [BENCHMARKED]

Major changes control variables->workflow variables, with operating regime-based definitions (DemandExpressionHandling.cpp; DEmandOptimization.cpp/.h)
-new routine CDemandOptimizer::ParseCondition
-report workflow variables in ManagementOptimization.csv
-:DefineWorkflowVariable,:WorkflowVarDefinition,:OverrideStageDischargeCurve (ParseManagementFile.cpp)

nested :RedirectToFile (ParseHRUFile.cpp;ParseManagementFile.cpp;ParseTimeSeriesFile.cpp)

support for :AllowSoilOverfill (CapillaryRise.cpp;Infiltration.cpp;Percolation.cpp;Recharge.cpp;ParseInput.cpp)

Update CMakeLists.txt to support LPSOLVE option (with help from M. Krossovski)

Moved handling of :OverrideStageDischargeCurve OUTSIDE of :ManagementGoal statements

prelim support for reduced memory usage in advection.cpp

Report whether observations are available in SubbasinProperties.csv file (StandardOutput.cpp)

bug fixes:
-support for proper parameters of RAINSNOW_HSPF (found by R. Chlumsky)
-correct OW_PET with snow cover as well as PET (for frozen lakes) (OrographicCorrections.cpp)
-breaks in :GaugedSubBasinGroup command if invalid syntax (found by R. Chlumsky) (ParseHRUFile.cpp)

Author: James Craig

CMakeLists.txt

@@ -1,11 +1,13 @@
 cmake_minimum_required(VERSION 3.20)
 # made with help from Zhu Liang-Jun of IGSNRR in Beijing
 # modified for GitHub by Trevor James Smith of Ouranos in Montreal
+# modified for lpsolve support by Maxim Krassovski of DFO
 
 # optional cmake command line arguments (e.g, "cmake -D COMPILE_LIB=ON" .)
 option(COMPILE_LIB "If ON, will create a dynamic lib file (default: OFF)" OFF)
 option(COMPILE_EXE "If ON, will create a executable file (default: ON)" ON)
 option(PYTHON, "If ON, will create a share library for python (default: OFF)" OFF)
+option(LPSOLVE,"If ON, will link to lp_solve optimization library (default: OFF)" OFF)
 
 # Setup Project
 PROJECT(Raven CXX)
@@ -14,14 +16,15 @@ PROJECT(Raven CXX)
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake")
 
 # Find NetCDF
-find_package(NetCDF)
+find_package(NetCDF) #may also be 'netCDF'
+find_package(netCDF)
 # find header & source
 file(GLOB HEADER "src/*.h")
 file(GLOB SOURCE "src/*.cpp")
 
 # Create library with Python bindings
 # To work, install pybind11, making sure it comes with cmake files (conda install -c conda-forge pybind11)
-if (PYTHON)
+if(PYTHON)
   SET(PYBIND11_NEWPYTHON ON)
   find_package(Python COMPONENTS REQUIRED Interpreter Development)
   find_package(pybind11 CONFIG REQUIRED)
@@ -44,15 +47,25 @@ if(COMPILE_EXE)
   )
   target_compile_definitions(Raven PUBLIC STANDALONE)
   set_target_properties(Raven PROPERTIES LINKER_LANGUAGE CXX)
+  if(LPSOLVE)
+    target_link_directories(Raven PRIVATE lib/lp_solve)  # where liblpsolve55.so is located
+    target_link_libraries(Raven lpsolve55)
+    add_definitions(-D_LPSOLVE_)
+  endif()
 endif()
+
 source_group("Header Files" FILES ${HEADER})
 source_group("Source Files" FILES ${SOURCE})
 
-IF(NETCDF_FOUND)
+if(NETCDF_FOUND)
   add_definitions(-Dnetcdf)
   include_directories(${NetCDF_INCLUDE_DIRS})
   target_link_libraries(Raven NetCDF::NetCDF)
-ENDIF()
+elseif(netCDF_FOUND)
+  add_definitions(-Dnetcdf)
+  include_directories(${NetCDF_INCLUDE_DIRS})
+  target_link_libraries(Raven netcdf)
+endif()
 
 # unset cmake variables to avoid polluting the cache
 unset(COMPILE_LIB CACHE)

src/Advection.cpp

@@ -25,7 +25,9 @@ CmvAdvection::CmvAdvection(string constit_name,
   _constit_ind=pTransModel->GetConstituentIndex(constit_name);
 
   int nAdvConnections    =pTransModel->GetNumAdvConnections();
+  int nWaterStores       =pTransModel->GetNumWaterCompartments();
   CHydroProcessABC::DynamicSpecifyConnections(3*nAdvConnections+1);//+2 once GW added
+  //CHydroProcessABC::DynamicSpecifyConnections(nAdvConnections+nWaterStores+1);//+2 once GW added
 
   for (int q=0;q<nAdvConnections;q++) //for regular advection
   {
@@ -42,18 +44,29 @@ CmvAdvection::CmvAdvection(string constit_name,
     iFrom[2*nAdvConnections+q]=pModel->GetStateVarIndex(CONSTITUENT_SRC,_constit_ind);
     iTo  [2*nAdvConnections+q]=iTo[q];
   }
+  //for Dirichlet source corrections
+  //for (int ii = 0; ii < nWaterStores; ii++) {
+  //  iFrom[nAdvConnections+ii]=pModel->GetStateVarIndex(CONSTITUENT_SRC,_constit_ind);
+  //  iTo  [nAdvConnections+ii]=pTransModel->GetStorWaterIndex(ii);
+  // }
+
   //advection into surface water
   int iSW=pModel->GetStateVarIndex(SURFACE_WATER);
   int   m=pTransModel->GetLayerIndex(_constit_ind,iSW);
   iFrom[3*nAdvConnections]=pModel->GetStateVarIndex(CONSTITUENT,m);
   iTo  [3*nAdvConnections]=pModel->GetStateVarIndex(CONSTITUENT_SW,0);
 
+  //iFrom[nAdvConnections+nWaterStores]=pModel->GetStateVarIndex(CONSTITUENT,m);
+  //iTo  [nAdvConnections+nWaterStores]=pModel->GetStateVarIndex(CONSTITUENT_SW,0);
+
   //advection into groundwater
   //int iGW=pModel->GetStateVarIndex(GROUNDWATER);
   //int   m=pTransModel->GetLayerIndex(constit_ind,iGW);
   //iFrom[3*nAdvConnections+1]=pModel->GetStateVarIndex(CONSTITUENT,m);
   //iTo  [3*nAdvConnections+1]=pModel->GetStateVarIndex(CONSTITUENT_GW,0);
 
+  //iFrom[nAdvConnections+nWaterStores+1]=pModel->GetStateVarIndex(CONSTITUENT,m);
+  //iTo  [nAdvConnections+nWaterStores+1]=pModel->GetStateVarIndex(CONSTITUENT_SW,0);
 }
 
 //////////////////////////////////////////////////////////////////
@@ -223,6 +236,8 @@ void   CmvAdvection::GetRatesOfChange(const double      *state_vars,
       mass          =sv[iFrom[q]];
       dirichlet_mass=Cs*sv[iFromWater]; //Cs*V
       rates[nAdvConnections+q]+=(dirichlet_mass-mass)/Options.timestep; //From CONSTITUENT_SRC
+      //iiFromWater=pTransModel->GetWaterStorIndexFromSVIndex(iFromWater);
+      //rates[nAdvConnections+iiFromWater]+=(dirichlet_mass-mass)/Options.timestep; //From CONSTITUENT_SRC //JRC: Reduced memory edit (or _aStorIndex[iFromWater])
       sv[iFrom[q]]=dirichlet_mass;      //override previous mass/enthalpy
     }
 
@@ -237,6 +252,8 @@ void   CmvAdvection::GetRatesOfChange(const double      *state_vars,
       mass          =sv[iTo[q]];
       dirichlet_mass=Cs*sv[iToWater]; //C*V
       rates[2*nAdvConnections+q]+=(dirichlet_mass-mass)/Options.timestep; //From CONSTITUENT_SRC
+      //iiToWater=pTransModel->GetWaterStorIndexFromSVIndex(iToWater);
+      //rates[nAdvConnections+iiToWater]+=(dirichlet_mass-mass)/Options.timestep; //From CONSTITUENT_SRC //JRC: Reduced memory edit (or _aStorIndex[iFromWater] of size _nWaterStores)
       sv[iTo[q]]=dirichlet_mass;      //override previous mass/enthalpy
     }
 

src/CapillaryRise.cpp

@@ -167,10 +167,11 @@ void   CmvCapillaryRise::ApplyConstraints( const double     *storage,
   double min_stor=g_min_storage;
 
   //cant remove more than is there
-  rates[0]=threshMin(rates[0],max(storage[iFrom[0]],min_stor)/Options.timestep,0.0);
+  rates[0]=min(rates[0],max(storage[iFrom[0]],min_stor)/Options.timestep);
 
   //exceedance of max "to" compartment
   //water flow simply slows (or stops) so that receptor will not overfill during tstep
-  rates[0]=threshMin(rates[0],
-                     (pHRU->GetStateVarMax(iTo[0],storage,Options)-storage[iTo[0]])/Options.timestep,0.0);
+  if (!Options.allow_soil_overfill){
+    rates[0]=min(rates[0],(pHRU->GetStateVarMax(iTo[0],storage,Options)-storage[iTo[0]])/Options.timestep);
+  }
 }