Merge branch 'main' into multiplatform-builds

Author: Zeitsperre

.pre-commit-config.yaml

@@ -26,7 +26,7 @@ repos:
 #      - id: include-what-you-use
 
   - repo: https://github.com/python-jsonschema/check-jsonschema
-    rev: 0.29.3
+    rev: 0.30.0
     hooks:
       - id: check-github-workflows
   - repo: meta

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,15 +16,16 @@ 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)
@@ -45,13 +48,22 @@ 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()
 
-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()
 
 source_group("Header Files" FILES ${HEADER})
 source_group("Source Files" FILES ${SOURCE})

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);
+  }
 }

src/CommonFunctions.cpp

@@ -50,7 +50,7 @@ string GetProcessName(process_type p)
   case(SNOWSQUEEZE):        {name="Liquid snow release";      break;}
   case(REFREEZE):           {name="Snow Refreeze";            break;}
   case(SUBLIMATION):        {name="Sublimation";              break;}
-  case(SNOW_BALANCE):       {name="Snow Melt & Refreeze";     break;}
+  case(SNOW_BALANCE):       {name="Snow Balance";             break;}
   case(GLACIER_MELT):       {name="Glacier Melt";             break;}
   case(GLACIER_RELEASE):    {name="Glacier Release";          break;}
   case(GLACIER_INFIL):      {name="Glacier Infiltration";     break;}

src/ConstituentModel.cpp

@@ -730,7 +730,7 @@ void CConstituentModel::WriteOutputFileHeaders(const optStruct &Options)
     _OUTPUT<<", sensible (r) "<<kg<<", conductive (r) " << kg << ", latent (r) "<<kg<<", GW mixing (r) "<<kg<<", radiant (r) "<<kg<<", friction (r)"<<kg<< endl;//TMP DEBUG
   }
 
-  //Pollutograph / stream temperatures file
+  // Pollutograph / stream temperatures file
   //--------------------------------------------------------------------
   if(!_is_passive) {
     if(_type!=ENTHALPY) { filename=_name+"Pollutographs.csv"; }
@@ -1262,14 +1262,14 @@ void CConstituentModel::WriteMinorOutput(const optStruct &Options,const time_str
 
   double latent_flux=0;
   double Q_sens(0.0),Q_cond(0.0),Q_lat(0.0),Q_GW(0.0),Q_rad(0.0),Q_fric(0.0);
-  double Qs,Qc,Ql,Qg,Qr,Qlw,Qlw_in,Qf,Tave;
+  double Qs,Qc,Ql,Qg,Qr,Qlw,Qlw_in,Qltrl,Qf,Tave;
   if(_type==ENTHALPY)
   {
     pEnthalpyModel=(CEnthalpyModel*)(this);
     latent_flux =pEnthalpyModel->GetAvgLatentHeatFlux()*(area*M2_PER_KM2)*Options.timestep; // [MJ] (loss term)t)
 
     for(int p=0;p<_pModel->GetNumSubBasins();p++) {
-      pEnthalpyModel->GetEnergyLossesFromReach(p,Qs,Qc,Ql,Qg,Qr,Qlw_in,Qlw,Qf,Tave);
+      pEnthalpyModel->GetEnergyLossesFromReach(p,Qs,Qc,Ql,Qg,Qr,Qlw_in,Qlw,Qltrl,Qf,Tave);
       Q_sens+=Qs;Q_cond+=Qc;Q_lat+=Ql;Q_GW+=Qg;Q_rad+=Qr+Qlw+Qlw_in; Q_fric+=Qf;
     }
   }
@@ -1529,7 +1529,7 @@ void CConstituentModel::WriteNetCDFMinorOutput(const optStruct& Options,const ti
 
   CEnthalpyModel *pEnthalpyModel=NULL; //Necessary evil to reduce code duplication rather than having an entirely separate WriteMinorOutput for Enthalpy
   if (_type==ENTHALPY){
-    CEnthalpyModel *pEnthalpyModel=(CEnthalpyModel*)(this);
+    pEnthalpyModel=(CEnthalpyModel*)(this);
   }
 
   if((Options.suppressICs) && (tt.model_time==0.0)) { return; }

src/ControlStructures.cpp

@@ -283,7 +283,7 @@ bool      COutflowRegime::AreConditionsMet(const time_struct& tt) const
 //////////////////////////////////////////////////////////////////
 /// \brief returns outflow from control structure for given stage, with prioritized constraints applied
 /// \param h [in] stage, in [m]
-/// \param Q_start [in] outflow at start of timestep [m3/s]
+/// \param Q_start [in] control structure outflow at start of timestep [m3/s]
 /// \returns outflow, in [m3/s]
 //
 double          COutflowRegime::GetOutflow(const double &h, const double &h_start, const double &Q_start, const long &target_SBID, const double &drefelev) const
@@ -292,7 +292,8 @@ double          COutflowRegime::GetOutflow(const double &h, const double &h_star
   double rivdepth = _pModel->GetSubBasinByID(target_SBID)->GetRiverDepth();
 
   double Q = _pCurve->GetDischarge(h, h_start, Q_start, rivdepth, drefelev);
-  double dQdt = (Q-Q_start)/tstep;
+
+  double dQdt = (Q-Q_start)/tstep; //m3/s/d
 
   //apply constraints here
   for (int j = _nConstraints-1; j >=0; j--) //priority from first to last

src/ControlStructures.h

@@ -78,7 +78,7 @@ class COutflowRegime
   void   AddRegimeConstraint(RegimeConstraint *pCond);
 
   bool      AreConditionsMet(const time_struct &tt) const;
-  double          GetOutflow(const double &h, const double &Qstart, const double &hstart, const long &target_SBID, const double &drefelev) const;
+  double          GetOutflow(const double &h, const double &hstart, const double &Qstart, const long &target_SBID, const double &drefelev) const;
 };