Skip to content

Commit 851e42f

Browse files
James CraigJames Craig
James Craig
authored and
James Craig
committed
v491 - Invasive overhaul - SBID and demandID as long long integer; DemandExpression [BENCHMARKED]
SBID and demandID as long long integer - many files; fixed some likely casting errors on the way -now supports SBIDs with more than 9 digits Prep work for ECCC Data assimilation - _aDAQadjust (Assimilate.cpp; Model.h/Model.cpp) support for :DemandExpression, :DemandFract (WaterDemands.cpp) Better warnings in Model initialize for time of conc, etc. (ModelInitialize.cpp)
1 parent f1ced19 commit 851e42f

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

42 files changed

+429
-401
lines changed

‎src/Assimilate.cpp‎

Lines changed: 35 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2021 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------*/
55
#include "RavenInclude.h"
66
#include "Model.h"
77

88
/*****************************************************************
99
Model Streamflow Assimilation Routines
1010
*****************************************************************/
11-
bool IsContinuousFlowObs2(const CTimeSeriesABC *pObs,long SBID)
11+
bool IsContinuousFlowObs2(const CTimeSeriesABC *pObs,long long SBID)
1212
{
1313
// clears up terribly ugly repeated if statements
1414
if(pObs==NULL) { return false; }
@@ -25,19 +25,21 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
2525
// Initialize Streamflow assimilation
2626
if(Options.assimilate_flow)
2727
{
28-
_aDAscale =new double[_nSubBasins];
29-
_aDAlength =new double[_nSubBasins];
30-
_aDAtimesince=new double[_nSubBasins];
31-
_aDAoverride =new bool [_nSubBasins];
32-
_aDAobsQ =new double[_nSubBasins];
33-
_aDAlast =new double[_nSubBasins];
28+
_aDAscale =new double[_nSubBasins];
29+
_aDAscale_last=new double[_nSubBasins];
30+
_aDAQadjust =new double[_nSubBasins];
31+
_aDAlength =new double[_nSubBasins];
32+
_aDAtimesince =new double[_nSubBasins];
33+
_aDAoverride =new bool [_nSubBasins];
34+
_aDAobsQ =new double[_nSubBasins];
3435
for(int p=0; p<_nSubBasins; p++) {
35-
_aDAscale [p]=1.0;
36-
_aDAlength [p]=0.0;
37-
_aDAtimesince[p]=0.0;
38-
_aDAoverride [p]=false;
39-
_aDAobsQ [p]=0.0;
40-
_aDAlast [p]=1.0;
36+
_aDAscale [p]=1.0;
37+
_aDAscale_last[p]=1.0;
38+
_aDAQadjust [p]=0.0;
39+
_aDAlength [p]=0.0;
40+
_aDAtimesince [p]=0.0;
41+
_aDAoverride [p]=false;
42+
_aDAobsQ [p]=0.0;
4143
}
4244
int count=0;
4345
for(int p=0; p<_nSubBasins; p++) {
@@ -93,17 +95,20 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
9395
//Option B: instantaneous flow
9496
Qmod = _pSubBasins[p]->GetOutflowRate();
9597
if(Qmod>PRETTY_SMALL) {
96-
_aDAscale[p]=1.0+alpha*((Qobs-Qmod)/Qmod);
98+
_aDAscale [p]=1.0+alpha*((Qobs-Qmod)/Qmod); //if alpha = 1, Q=Qobs in observation basin
99+
_aDAQadjust[p]=alpha*(Qobs-Qmod);
97100
}
98101
else {
99-
_aDAscale[p]=1.0;
102+
_aDAscale [p]=1.0;
103+
_aDAQadjust[p]=0.0;
100104
}
101-
_aDAlast[p]=1.0;//no need to scale previous
105+
_aDAscale_last[p]=1.0;//no need to scale previous
102106
}
103107

104108
// actually scale flows
105109
//---------------------------------------------------------------------
106-
double mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAlast[p],_aDAoverride[p],Options.timestep,tt.model_time);
110+
double mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time);
111+
//double mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAadjust[p],_aDAoverride[p],Options.timestep,tt.model_time);
107112
if(mass_added>0.0){_CumulInput +=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
108113
else {_CumulOutput-=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
109114

@@ -119,17 +124,17 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
119124
//
120125
void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
121126
{
122-
if (!Options.assimilate_flow) {return;}
123-
if(tt.model_time<(Options.assimilation_start+Options.timestep/2.0)){return;}//assimilates all data after assimilation date
127+
if (!Options.assimilate_flow) {return;}
128+
if (tt.model_time<(Options.assimilation_start-Options.timestep/2.0)){return;}//assimilates all data after or on assimilation date
124129

125130
int p,pdown;
126131
double Qobs;
127132
double t_observationsOFF=ALMOST_INF;//Only used for debugging - keep as ALMOST_INF otherwise
128133

129-
int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep);//current timestep index
134+
int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep)+1;//end of timestep index
130135

131136
for(p=0; p<_nSubBasins; p++) {
132-
_aDAlast[p]=_aDAscale[p];
137+
_aDAscale_last[p]=_aDAscale[p];
133138
}
134139

135140
for(int pp=_nSubBasins-1; pp>=0; pp--)//downstream to upstream
@@ -145,22 +150,24 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
145150
{
146151
if(IsContinuousFlowObs2(_pObservedTS[i],_pSubBasins[p]->GetID()))//flow observation is available and linked to this subbasin
147152
{
148-
Qobs = _pObservedTS[i]->GetSampledValue(nn+1);//+1: correction for period ending storage of hydrograph
153+
Qobs = _pObservedTS[i]->GetSampledValue(nn); //end of timestep flow
149154

150155
//bool fakeblank=((tt.model_time>30) && (tt.model_time<40)) || ((tt.model_time>45) && (tt.model_time<47));//TMP DEBUG
151156
//if (fakeblank){Qobs=RAV_BLANK_DATA;}
152157

153158
if((Qobs!=RAV_BLANK_DATA) && (tt.model_time<t_observationsOFF))
154159
{
155160
//_aDAscale[p] calculated live in AssimilationOverride when up-to-date modelled flow available
161+
//same with _aDAQadjust[p]
156162
_aDAlength [p]=0.0;
157163
_aDAtimesince[p]=0.0;
158164
_aDAoverride [p]=true;
159165
_aDAobsQ [p]=Qobs;
160166
}
161167
else
162168
{ //found a blank or zero flow value
163-
_aDAscale [p]=_aDAscale[p];
169+
_aDAscale [p]=_aDAscale[p];//same adjustment as before - scaling persists
170+
_aDAQadjust [p]=_aDAQadjust[p];//same adjustment as before - flow magnitude persists
164171
_aDAtimesince[p]+=Options.timestep;
165172
_aDAlength [p]=0.0;
166173
_aDAoverride [p]=false;
@@ -178,12 +185,14 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
178185
//if ((pdown!=DOESNT_EXIST) && (!_aDAoverride[pdown])){ //alternate - allow information to pass through reservoirs
179186
if((pdown!=DOESNT_EXIST) && (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[pdown])) {
180187
_aDAscale [p]= _aDAscale [pdown];
188+
_aDAQadjust [p]= _aDAQadjust [pdown] * (_pSubBasins[p]->GetDrainageArea() / _pSubBasins[pdown]->GetDrainageArea());
181189
_aDAlength [p]+=_pSubBasins [pdown]->GetReachLength();
182190
_aDAtimesince [p]= _aDAtimesince[pdown];
183191
_aDAoverride [p]=false;
184192
}
185193
else{ //Nothing downstream or reservoir present in this basin, no assimilation
186194
_aDAscale [p]=1.0;
195+
_aDAQadjust [p]=0.0;
187196
_aDAlength [p]=0.0;
188197
_aDAtimesince[p]=0.0;
189198
_aDAoverride [p]=false;
@@ -197,6 +206,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
197206
double distfact = _pGlobalParams->GetParams()->assim_upstream_decay/M_PER_KM; //[1/km]->[1/m]
198207
for(p=0; p<_nSubBasins; p++)
199208
{
200-
_aDAscale[p] =1.0+(_aDAscale[p]-1.0)*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
209+
_aDAscale [p] =1.0+(_aDAscale[p]-1.0)*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
210+
_aDAQadjust[p] =_aDAQadjust[p]*exp(-distfact*_aDAlength[p])*exp(-time_fact*_aDAtimesince[p]);
201211
}
202212
}

‎src/ChannelXSect.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2022 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------*/
55
#include "ChannelXSect.h"
66
void TestManningsInfluence(const CChannelXSect *pChan,const double &Qref);
@@ -403,7 +403,7 @@ double CChannelXSect::GetDiffusivity(const double &Q, const double &SB_slope, co
403403
/// \param SB_n [in] subbasin mannings (or AUTO_COMPUTE, if channel mannings to be used)
404404
/// \param SBID [in] subbasin identifier
405405
//
406-
void CChannelXSect::CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long SBID) const
406+
void CChannelXSect::CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long long SBID) const
407407
{
408408
string warn;
409409
double junk,Q_mult;

‎src/ChannelXSect.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2022 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------
55
ChannelXSect.h
66
------------------------------------------------------------------
@@ -104,6 +104,6 @@ class CChannelXSect
104104
double GetCelerity (const double &Qref, const double &SB_slope,const double &SB_n) const;
105105
double GetDiffusivity(const double &Q, const double &SB_slope, const double &SB_n) const;
106106

107-
void CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long SBID) const;
107+
void CheckReferenceFlow(const double& Qref,const double& SB_slope,const double& SB_n,const long long SBID) const;
108108
};
109109
#endif

‎src/ConstituentModel.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*----------------------------------------------------------------
22
Raven Library Source Code
3-
Copyright (c) 2008-2022 the Raven Development Team
3+
Copyright (c) 2008-2025 the Raven Development Team
44
----------------------------------------------------------------
55
Constituent Transport/Tracer Model class
66
coordinates information about constituent storage
@@ -12,7 +12,7 @@ coordinates information about constituent storage
1212
#include "Transport.h"
1313
#include "EnergyTransport.h"
1414

15-
bool IsContinuousConcObs(const CTimeSeriesABC *pObs,const long SBID,const int c); //Defined in StandardOutput.cpp
15+
bool IsContinuousConcObs(const CTimeSeriesABC *pObs,const long long SBID,const int c); //Defined in StandardOutput.cpp
1616
void WriteNetCDFGlobalAttributes(const int out_ncid,const optStruct& Options,const string descript);
1717
int NetCDFAddMetadata (const int fileid,const int time_dimid,string shortname,string longname,string units);
1818
int NetCDFAddMetadata2D (const int fileid,const int time_dimid,int nbasins_dimid,string shortname,string longname,string units);

‎src/ControlStructures.cpp‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111
/// \param downID [in] downstream subbasin
1212
/// \notes assumes SBID basin has valid reservoir
1313
//
14-
CControlStructure::CControlStructure(string name, const long SBID,const long downID) {
14+
CControlStructure::CControlStructure(string name, const long long SBID,const long long downID) {
1515
_name=name;
1616
_SBID=SBID;
1717
_target_SBID=downID;
@@ -37,13 +37,13 @@ string CControlStructure::GetName () const {return _name;}
3737
/// \brief returns target subbasin ID to which outflow is directed
3838
/// \returns target subbasin ID
3939
//
40-
long CControlStructure::GetTargetBasinID() const {return _target_SBID;}
40+
long long CControlStructure::GetTargetBasinID() const {return _target_SBID;}
4141

4242
//////////////////////////////////////////////////////////////////
4343
/// \brief sets target subbasin ID to which outflow is directed
4444
/// \param SBID [in] - valid subbasin ID (or DOESNT_EXIST, if flow is directed outside model)
4545
//
46-
void CControlStructure::SetTargetBasin(const long SBID)
46+
void CControlStructure::SetTargetBasin(const long long SBID)
4747
{
4848
_target_SBID=SBID;
4949
}
@@ -228,7 +228,7 @@ bool COutflowRegime::AreConditionsMet(const time_struct& tt) const
228228
for (int i = 0; i < _nConditions; i++)
229229
{
230230
comparison comp=_pConditions[i]->compare;
231-
long SBID=_pConditions[i]->basinID;
231+
long long SBID=_pConditions[i]->basinID;
232232
double v1=_pConditions[i]->compare_val1;
233233
double v2=_pConditions[i]->compare_val2;
234234
string var=_pConditions[i]->variable;
@@ -286,7 +286,7 @@ bool COutflowRegime::AreConditionsMet(const time_struct& tt) const
286286
/// \param Q_start [in] control structure outflow at start of timestep [m3/s]
287287
/// \returns outflow, in [m3/s]
288288
//
289-
double COutflowRegime::GetOutflow(const double &h, const double &h_start, const double &Q_start, const long &target_SBID, const double &drefelev) const
289+
double COutflowRegime::GetOutflow(const double &h, const double &h_start, const double &Q_start, const long long &target_SBID, const double &drefelev) const
290290
{
291291
double tstep = _pModel->GetOptStruct()->timestep;
292292
double rivdepth = _pModel->GetSubBasinByID(target_SBID)->GetRiverDepth();

‎src/ControlStructures.h‎

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,8 @@
1717
*****************************************************************/
1818
struct RegimeCondition
1919
{
20-
string variable; ///< one of {DATE, DOY, STAGE, FLOW, ...}
21-
long basinID; ///< defaults to -1 (this basin)
20+
string variable; ///< one of {DATE, DOY, STAGE, FLOW, ...}
21+
long long basinID; ///< defaults to -1 (this basin)
2222
comparison compare; ///< e.g., greater than, between, etc.
2323
double compare_val1; ///< primary comparison value
2424
double compare_val2; ///< secondary comparison value (used for between)
@@ -78,7 +78,7 @@ class COutflowRegime
7878
void AddRegimeConstraint(RegimeConstraint *pCond);
7979

8080
bool AreConditionsMet(const time_struct &tt) const;
81-
double GetOutflow(const double &h, const double &hstart, const double &Qstart, const long &target_SBID, const double &drefelev) const;
81+
double GetOutflow(const double &h, const double &hstart, const double &Qstart, const long long &target_SBID, const double &drefelev) const;
8282
};
8383

8484

@@ -91,28 +91,28 @@ class CControlStructure
9191
{
9292
private:/*-------------------------------------------------------*/
9393
string _name; ///< structure name
94-
long _SBID; ///< subbasin ID
95-
long _target_SBID; ///< outflow directed to this subbasin ID
94+
long long _SBID; ///< subbasin ID
95+
long long _target_SBID; ///< outflow directed to this subbasin ID
9696
double _dRefElev; ///< downstream reference elevation [m] (used in some structures for backwater, limiting flow, etc.)
9797

9898
int _nRegimes; ///< number of flow regimes
9999
COutflowRegime **_aRegimes; ///< array of pointers to flow regimes with unique Q=f(h,Q,...)
100100

101101
public:/*-------------------------------------------------------*/
102102
//Constructors:
103-
CControlStructure(const string name,const long SBID, const long downID);
103+
CControlStructure(const string name,const long long SBID, const long long downID);
104104
~CControlStructure();
105105

106-
void AddRegime (COutflowRegime *pRegime);
107-
void SetTargetBasin (const long SBID);
106+
void AddRegime (COutflowRegime *pRegime);
107+
void SetTargetBasin (const long long SBID);
108108

109-
void SetDownstreamRefElevation (const double dRefElev);
109+
void SetDownstreamRefElevation (const double dRefElev);
110110

111-
string GetName () const;
112-
long GetTargetBasinID() const;
111+
string GetName () const;
112+
long long GetTargetBasinID() const;
113113

114-
double GetOutflow(const double &stage,const double &stage_start,const double &Q_start,const time_struct &tt) const;
114+
double GetOutflow(const double &stage,const double &stage_start,const double &Q_start,const time_struct &tt) const;
115115

116-
string GetCurrentRegimeName(const time_struct &tt) const;
116+
string GetCurrentRegimeName(const time_struct &tt) const;
117117
};
118118
#endif

‎src/DDS.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ void CDDSEnsemble::SetPerturbationValue(const double &perturb) {
6767
//////////////////////////////////////////////////////////////////
6868
/// \brief set DDS calibration target
6969
//
70-
void CDDSEnsemble::SetCalibrationTarget(const long SBID,const diag_type object_diag,const string period)
70+
void CDDSEnsemble::SetCalibrationTarget(const long long SBID,const diag_type object_diag,const string period)
7171
{
7272
_calib_SBID =SBID;
7373
_calib_Obj =object_diag;

‎src/DemandExpressionHandling.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -231,7 +231,7 @@ int CDemandOptimizer::GetIndexFromDVString(string s) const //String in format !Q
231231
return DOESNT_EXIST;
232232
}
233233
else{
234-
long SBID=s_to_l(s.substr(2).c_str());
234+
long long SBID=s_to_ll(s.substr(2).c_str());
235235
return _pModel->GetSubBasinIndex(SBID);
236236
}
237237
}
@@ -644,7 +644,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t
644644
else if (s.substr(0, 8) == "@SB_var(") // SubBasin state var (e.g., @SB_var(SNOW,[id])
645645
{
646646
string sv_name;
647-
long SBID;
647+
long long SBID;
648648
size_t is = s.find("@SB_var(");
649649
size_t ie = s.find(",",is);
650650
size_t ip = s.find_last_of(")");
@@ -660,7 +660,7 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t
660660
{
661661
bool found=false;
662662
sv_name = s.substr(is+8,ie-(is+8));
663-
SBID = s_to_l(s.substr(ie+1,ip-(ie+1)).c_str());
663+
SBID = s_to_ll(s.substr(ie+1,ip-(ie+1)).c_str());
664664

665665
int lay=DOESNT_EXIST;
666666
sv_type sv=_pModel->GetStateVarInfo()->StringToSVType(sv_name,lay,false);
@@ -1002,7 +1002,7 @@ exp_condition* CDemandOptimizer::ParseCondition(const char** s, const int Len, c
10021002
char code=pCond->dv_name[1];
10031003
if ((code=='Q') || (code=='h') || (code=='I')) //subbasin state decision variable
10041004
{
1005-
long SBID=s_to_l(tmp2.c_str());
1005+
long long SBID=s_to_ll(tmp2.c_str());
10061006
if (_pModel->GetSubBasinByID(SBID) == NULL) {
10071007
ExitGracefully("ParseManagementFile: Subbasin ID in :Condition statement is invalid.",BAD_DATA_WARN);
10081008
}

‎src/DemandGroups.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -50,10 +50,10 @@ int CDemandGroup::GetGlobalIndex () const {return _global_ii;}
5050
///
5151
/// \return true if demand with ID demandID is in group
5252
//
53-
bool CDemandGroup::IsInGroup (const int demandID) const
53+
bool CDemandGroup::IsInGroup (const long long demandID) const
5454
{
5555
for (int p=0;p<_nDemands; p++){
56-
if (_pDemands[p]->GetID() == demandID) { return true; }
56+
if (_pDemands[p]->GetDemandID() == demandID) { return true; }
5757
}
5858
return false;
5959
}

0 commit comments

Comments
 (0)