Skip to content

Commit f9c77a2

Browse files
James CraigJames Craig
James Craig
authored and
James Craig
committed
v529 - fixed direct insertion assimilation issues
adjusted assimilation to aim for 0.5*(Qobs^n+Qobs^n+1) rather than exactly match Qobs- oscillation issue -adjuested in DemandOptimization.cpp, SubBasin.cpp, Assimilate.cpp -new DA variables (Model.h/.cpp) -handling of flow from Reservoirs (Reservoir.cpp) Proper handling of workflow variables in non-linear loop added WETLAND to HRUTYPE conditionals in :PopulateHRUGroup (ParseHRUFile.cpp) improved QA/QC for non-linear management variables (DemandExpressionHandling.cpp) bug: fixed erroneous check for special soil horizions (ParsePropertyFile.cpp)
1 parent b788f64 commit f9c77a2

File tree

11 files changed

+115
-66
lines changed

11 files changed

+115
-66
lines changed

‎src/Assimilate.cpp‎

Lines changed: 62 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,8 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
3434
_aDAtimesince =new double[_nSubBasins];
3535
_aDAoverride =new bool [_nSubBasins];
3636
_aDAobsQ =new double[_nSubBasins];
37+
_aDAobsQ2 =new double[_nSubBasins];
38+
_aDASinceLastBlank=new double [_nSubBasins];
3739
for(int p=0; p<_nSubBasins; p++) {
3840
_aDAscale [p]=1.0;
3941
_aDAscale_last[p]=1.0;
@@ -44,6 +46,8 @@ void CModel::InitializeDataAssimilation(const optStruct &Options)
4446
_aDAtimesince [p]=0.0;
4547
_aDAoverride [p]=false;
4648
_aDAobsQ [p]=0.0;
49+
_aDAobsQ2 [p]=0.0;
50+
_aDASinceLastBlank[p]=0.0; //initialized to zero to ramp from IC (default)
4751
}
4852
int count=0;
4953
for(int p=0; p<_nSubBasins; p++) {
@@ -71,17 +75,28 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
7175

7276
//Get current, up-to-date flow and calculate scaling factor
7377
//---------------------------------------------------------------------
78+
double Qobs;
7479
if(_aDAoverride[p])
7580
{
76-
double Qobs,Qmod,Qmodlast;
81+
double Qobs2,Qmod,Qmodlast;
7782
double alpha = _pGlobalParams->GetParams()->assimilation_fact;
7883

84+
//alpha*=(1.0-exp(-0.06*_aDASinceLastBlank[p])); //shock/oscillation prevention
85+
7986
Qobs = _aDAobsQ[p];
87+
Qobs2 = _aDAobsQ2[p];
8088
Qmod = _pSubBasins[p]->GetOutflowRate();
8189
Qmodlast= _pSubBasins[p]->GetLastOutflowRate();
8290
if(Qmod>PRETTY_SMALL) {
8391
_aDAscale [p]=1.0+alpha*((Qobs-Qmod)/Qmod); //if alpha = 1, Q=Qobs in observation basin
84-
_aDAQadjust[p]=0.5*alpha*(2.0*Qobs-Qmodlast-Qmod);//Option B: mean flow
92+
//_aDAQadjust[p]=alpha*(Qobs-Qmod); //Option A: end of time step flow
93+
//_aDAQadjust[p]=0.5*alpha*(2.0*Qobs-Qmodlast-Qmod);//Option B: mean flow - rapidly oscillatory Q - Qmean is perfect
94+
if (Qobs2 == RAV_BLANK_DATA) { //Option C: aim for midpoint of two observation flows (this is the way)
95+
_aDAQadjust[p]=alpha*(Qobs-Qmod);
96+
}
97+
else {
98+
_aDAQadjust[p]=alpha*(0.5*(Qobs2+Qobs)-Qmod);
99+
}
85100
}
86101
else {
87102
_aDAscale [p]=1.0;
@@ -92,17 +107,16 @@ void CModel::AssimilationOverride(const int p,const optStruct& Options,const tim
92107

93108
// actually scale flows
94109
//---------------------------------------------------------------------
95-
double mass_added;
110+
double mass_added=0.0;
96111
if (Options.assim_method==DA_RAVEN_DEFAULT){
97112
mass_added=_pSubBasins[p]->ScaleAllFlows(_aDAscale[p]/_aDAscale_last[p],_aDAoverride[p],Options.timestep,tt.model_time);
98113
}
99114
else if (Options.assim_method==DA_ECCC) {
100115
if(_aDAoverride[p]){
101-
mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],true,Options.timestep,tt.model_time);
116+
mass_added=_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],true,true,Options.timestep,tt.model_time);
102117
}
103118
}
104119

105-
//
106120
if(mass_added>0.0){_CumulInput +=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
107121
else {_CumulOutput-=mass_added/(_WatershedArea*M2_PER_KM2)*MM_PER_METER;}
108122

@@ -122,14 +136,17 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
122136
if (tt.model_time<(Options.assimilation_start-Options.timestep/2.0)){return;}//assimilates all data after or on assimilation date
123137

124138
int p,pdown;
125-
double Qobs;
139+
double Qobs,Qobs2;
126140
double t_observationsOFF=ALMOST_INF;//Only used for debugging - keep as ALMOST_INF otherwise
141+
bool ObsExists; //observation available in THIS basin
127142

128143
for(p=0; p<_nSubBasins; p++) {
129144
_aDAscale_last[p]=_aDAscale[p];
130145
_aDADrainSum[p]=0.0;
131146
}
132147

148+
// check for observations in each basin
149+
//----------------------------------------------------------------
133150
int nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep)+1;//end of timestep index
134151

135152
for(int pp=_nSubBasins-1; pp>=0; pp--)//downstream to upstream
@@ -138,7 +155,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
138155

139156
pdown=GetDownstreamBasin(p);
140157

141-
bool ObsExists=false; //observation available in THIS basin
158+
ObsExists=false;
142159

143160
if(_pSubBasins[p]->UseInFlowAssimilation())
144161
{
@@ -147,6 +164,7 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
147164
if(IsContinuousFlowObs2(_pObservedTS[i],_pSubBasins[p]->GetID()))//flow observation is available and linked to this subbasin
148165
{
149166
Qobs = _pObservedTS[i]->GetSampledValue(nn); //mean timestep flow
167+
Qobs2 = _pObservedTS[i]->GetSampledValue(nn+1); //mean timestep flow
150168
ObsExists=true;
151169
break; //avoids duplicate observations
152170
}
@@ -165,6 +183,8 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
165183
_aDAtimesince[p]=0.0;
166184
_aDAoverride [p]=true;
167185
_aDAobsQ [p]=Qobs;
186+
_aDAobsQ2 [p]=Qobs2;
187+
_aDASinceLastBlank[p]+=1.0;
168188
}
169189
else
170190
{ //found a blank or zero flow value
@@ -173,14 +193,19 @@ void CModel::PrepareAssimilation(const optStruct &Options,const time_struct &tt)
173193
_aDAtimesince[p]+=Options.timestep;
174194
_aDAoverride [p]=false;
175195
_aDAobsQ [p]=0.0;
196+
_aDAobsQ2 [p]=0.0;
197+
_aDASinceLastBlank[p]+=0.0;
176198
}
177199
}
178200
// no observations in this basin, get scaling from downstream
179201
//----------------------------------------------------------------
180202
else if(!ObsExists) //observations may be downstream, propagate scaling upstream
181203
{
182204
//if ((pdown!=DOESNT_EXIST) && (!_aDAoverride[pdown])){ //alternate - allow information to pass through reservoirs
183-
if( (pdown!=DOESNT_EXIST) && (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[p]) && (_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled())) {
205+
if( (pdown!=DOESNT_EXIST) &&
206+
(_pSubBasins[p]->GetReservoir()==NULL) &&
207+
(!_aDAoverride[p]) &&
208+
(_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled())) {
184209
_aDAscale [p]= _aDAscale [pdown];
185210
_aDAlength [p]+=_pSubBasins [pdown]->GetReachLength();
186211
_aDAtimesince [p]= _aDAtimesince[pdown];
@@ -258,54 +283,55 @@ void CModel::AssimilationBackPropagate(const optStruct &Options,const time_struc
258283

259284
if (!Options.assimilate_flow) {return;}
260285

286+
// Determine flow adjustment magnitude
287+
//----------------------------------------------------------------
261288
for(int pp=_nSubBasins-1; pp>=0; pp--)//downstream to upstream
262289
{
263-
p=GetOrderedSubBasinIndex(pp);
264-
265-
pdown=DOESNT_EXIST;
266-
if (_pSubBasins[p]->GetDownstreamID()!=DOESNT_EXIST){
267-
pdown = GetSubBasinByID(_pSubBasins[p]->GetDownstreamID())->GetGlobalIndex();
268-
}
269-
270-
// no observations in this basin, get scaling from downstream
271-
//----------------------------------------------------------------
290+
p =GetOrderedSubBasinIndex(pp);
272291
pdown=GetDownstreamBasin(p);
273-
long long SBID=_pSubBasins[p]->GetID();
274-
if(!_aDAoverride[p]) { //observations may be downstream, propagate scaling upstream
275-
if( (pdown!=DOESNT_EXIST) && (_pSubBasins[p]->GetReservoir()==NULL) && (!_aDAoverride[p]) && (_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled())) {
276-
_aDAQadjust [p]= _aDAQadjust [pdown] * (_pSubBasins[p]->GetDrainageArea() / _pSubBasins[pdown]->GetDrainageArea());
277-
//cout<<"PROPAGATING "<<SBID<<": " <<setprecision(3)<< _aDAQadjust[p] << " from " << _aDAQadjust[pdown] << endl;
292+
293+
if(!_aDAoverride[p]) { //observations may be downstream, get scaling from downstream
294+
if( (pdown!=DOESNT_EXIST ) &&
295+
(!_aDAoverride[p] ) &&
296+
(_pSubBasins[p]->GetReservoir()==NULL) &&
297+
(_pSubBasins[p]->IsEnabled()) && (_pSubBasins[pdown]->IsEnabled()))
298+
{
299+
_aDAQadjust[p]= _aDAQadjust [pdown] * (_pSubBasins[p]->GetDrainageArea() / _pSubBasins[pdown]->GetDrainageArea());
278300
}
279-
else{ //Nothing downstream or reservoir present in this basin, no assimilation
280-
_aDAQadjust [p]=0.0;
301+
else{ //Nothing downstream or reservoir present in this basin, no adjustment
302+
_aDAQadjust[p]= 0.0;
281303
}
282304
}
283-
else if (_pSubBasins[p]->IsEnabled()) {
284-
//cout<<"ASSIMILATING AT "<<SBID<<": " << _aDAQadjust[p] << endl;
285-
}
286-
} // end downstream to upstream
305+
}
287306

288307
// Apply time and space correction factors
289308
//----------------------------------------------------------------
290309
double distfact = _pGlobalParams->GetParams()->assim_upstream_decay/M_PER_KM; //[1/km]->[1/m]
291310
double ECCCwt;
292311
for(p=0; p<_nSubBasins; p++)
293312
{
294-
_aDAQadjust[p] *=exp(-distfact*_aDAlength[p]);
295-
ECCCwt=1.0;
296-
if (_aDADrainSum[p]!=0.0){
297-
ECCCwt = (_pSubBasins[p]->GetDrainageArea() - _aDADrainSum[p])/(_aDADownSum[p]-_aDADrainSum[p]);
298-
}
313+
if(!_aDAoverride[p])
314+
{
315+
pdown=GetDownstreamBasin(p);
316+
if (pdown!=DOESNT_EXIST){
317+
_aDAQadjust[p] *=exp(-distfact*_aDAlength[pdown]);
318+
}
299319

300-
if (!_aDAoverride[p]) { //no scaling for stations being overridden, only upstream
301-
_aDAQadjust[p] = _aDAQadjust[p]*ECCCwt;
320+
ECCCwt=1.0;
321+
if (_aDADrainSum[p]!=0.0){
322+
ECCCwt = (_pSubBasins[p]->GetDrainageArea() - _aDADrainSum[p])/(_aDADownSum[p]-_aDADrainSum[p]);
323+
}
324+
325+
if (!_aDAoverride[p]) { //no scaling for stations being overridden, only upstream
326+
_aDAQadjust[p] = _aDAQadjust[p]*ECCCwt;
327+
}
302328
}
303329
}
304330

305331
// Actually update flows
306332
//----------------------------------------------------------------
307333
for(p=0; p<_nSubBasins; p++)
308334
{
309-
_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],false,Options.timestep,tt.model_time);
335+
_pSubBasins[p]->AdjustAllFlows(_aDAQadjust[p],false,_aDAoverride[p],Options.timestep,tt.model_time);
310336
}
311337
}

‎src/DemandExpressionHandling.cpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -521,6 +521,11 @@ bool CDemandOptimizer::ConvertToExpressionTerm(const string s, expressionTerm* t
521521
else if (s[0] == '?') {
522522
term->type=TERM_ITER;
523523
term->DV_ind=GetNLIndexFromGuessString(s);
524+
if (term->DV_ind == DOESNT_EXIST) {
525+
warn="ConvertToExpressionTerm: non-linear expression in expression " +warnstring +" was not defined: goal/constraint will be ignored";
526+
WriteWarning(warn.c_str(),true);
527+
return false;
528+
}
524529
}
525530
//----------------------------------------------------------------------
526531
else if (s.substr(0, 4) == "@ts(")//time series (e.g., @ts(my_time_series,n)

‎src/DemandOptimization.cpp‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1709,7 +1709,7 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
17091709
int p2,id;
17101710
dv_type dvtype;
17111711
bool assimilating;
1712-
double Qobs;
1712+
double Qobs,Qobs2;
17131713

17141714
for (int pp = 0; pp<pModel->GetNumSubBasins(); pp++)
17151715
{
@@ -1729,7 +1729,8 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
17291729
SBID=pSB->GetID();
17301730

17311731
i=0;
1732-
Qobs=_pModel->GetObservedFlow(p,nn);
1732+
Qobs =_pModel->GetObservedFlow(p,nn );
1733+
Qobs2=_pModel->GetObservedFlow(p,nn+1);
17331734
assimilating = ((Options.assimilate_flow) && (pSB->UseInFlowAssimilation()) && (Qobs!=RAV_BLANK_DATA) && (pSB->GetReservoir()==NULL));
17341735

17351736
// outflow term =============================
@@ -1799,7 +1800,8 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
17991800
}
18001801
else //assimilating - skip mass balance
18011802
{
1802-
RHS=Qobs;
1803+
if (Qobs2 != RAV_BLANK_DATA) {RHS=0.5*(Qobs+Qobs2);}
1804+
else {RHS=Qobs; }
18031805
}
18041806

18051807
retval = lp_lib::add_constraintex(pLinProg,i,row_val,col_ind,ROWTYPE_EQ,RHS);
@@ -2392,6 +2394,10 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
23922394
AddConstraintToLP( i, _pGoals[i]->active_regime, pLinProg, tt, col_ind, row_val,true,lpgoalrow[i]);
23932395
}
23942396
}
2397+
// adjust workflow variables to respond to non-linear DVs
2398+
// ----------------------------------------------------------------------------
2399+
UpdateWorkflowVariables(tt,Options);
2400+
23952401
iter++;
23962402
} while ((iter<_maxIterations));/*end iteration loop*/
23972403

‎src/Model.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,8 @@ CModel::CModel(const int nsoillayers,
6767
_aDAtimesince =NULL;
6868
_aDAoverride =NULL;
6969
_aDAobsQ =NULL;
70+
_aDAobsQ2 =NULL;
71+
_aDASinceLastBlank=NULL;
7072

7173
_pOptStruct = &Options;
7274
_pGlobalParams = new CGlobalParams();
@@ -223,6 +225,8 @@ CModel::~CModel()
223225
delete [] _aDAtimesince; _aDAtimesince=NULL;
224226
delete [] _aDAoverride; _aDAoverride=NULL;
225227
delete [] _aDAobsQ; _aDAobsQ=NULL;
228+
delete [] _aDAobsQ2; _aDAobsQ2=NULL;
229+
delete [] _aDASinceLastBlank; _aDASinceLastBlank=NULL;
226230

227231
this->DestroyAllLanduseClasses();
228232
this->DestroyAllSoilClasses();

‎src/Model.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -151,6 +151,8 @@ class CModel: public CModelABC
151151
double *_aDAtimesince; ///< array of downstream time since most recent downstream DA observation [size: _nSubBasins] (NULL w/o DA)
152152
bool *_aDAoverride; ///< array of booleans indicating if observation data is available for assimilation at basin p's outlet [size: _nSubBasins] (NULL w/o DA)
153153
double *_aDAobsQ; ///< array of observed flow values in basins [size: _nSubBasins] (NULL w/o DA)
154+
double *_aDAobsQ2; ///< array of observed flow values in next time step [size: _nSubBasins] (NULL w/o DA)
155+
double *_aDASinceLastBlank; ///< array of time steps since last blank value occurred [size: _nSubBasins] (NULL w/o DA)
154156
double *_aDAscale_last; ///< array of scale factors from previous time step [size: _nSubBasins] (NULL w/o DA)
155157

156158
force_perturb**_pPerturbations; ///< array of pointers to perturbation data; defines which forcing functions to perturb and how [size: _nPerturbations]

‎src/ParseHRUFile.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -628,6 +628,7 @@ bool ParseHRUPropsFile(CModel *&pModel, const optStruct &Options, bool terrain_r
628628
if((pModel->GetHydroUnit(k)->GetHRUType()==HRU_GLACIER) && (!strcmp(s[4],"GLACIER"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
629629
if((pModel->GetHydroUnit(k)->GetHRUType()==HRU_LAKE) && (!strcmp(s[4],"LAKE"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
630630
if((pModel->GetHydroUnit(k)->GetHRUType()==HRU_WATER) && (!strcmp(s[4],"WATER"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
631+
if((pModel->GetHydroUnit(k)->GetHRUType()==HRU_WETLAND) && (!strcmp(s[4],"WETLAND"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
631632
}
632633
}
633634
if(!strcmp(s[4],"DOESNTEQUAL")){
@@ -638,6 +639,7 @@ bool ParseHRUPropsFile(CModel *&pModel, const optStruct &Options, bool terrain_r
638639
if((pModel->GetHydroUnit(k)->GetHRUType()!=HRU_GLACIER) && (!strcmp(s[4],"GLACIER"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
639640
if((pModel->GetHydroUnit(k)->GetHRUType()!=HRU_LAKE) && (!strcmp(s[4],"LAKE"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
640641
if((pModel->GetHydroUnit(k)->GetHRUType()==HRU_WATER) && (!strcmp(s[4],"WATER"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
642+
if((pModel->GetHydroUnit(k)->GetHRUType()==HRU_WETLAND) && (!strcmp(s[4],"WETLAND"))){ pHRUGrp->AddHRU(pModel->GetHydroUnit(k)); }
641643
}
642644
}
643645
}

‎src/ParsePropertyFile.cpp‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -460,12 +460,13 @@ bool ParseClassPropertiesFile(CModel *&pModel,
460460
//zero is a valid entry for lakes & glaciers (zero horizons)
461461
ExitGracefullyIf(Len!=(nhoriz*2+2),
462462
"ParseClassPropertiesFile: :SoilProfiles invalid command length",BAD_DATA);
463-
bool is_special=((!string(s[0]).compare("LAKE")) ||
464-
(!string(s[0]).compare("GLACIER")) ||
465-
(!string(s[0]).compare("PAVEMENT")) ||
466-
(!string(s[0]).compare("WATER")) ||
467-
(!string(s[0]).compare("WETLAND")) ||
468-
(!string(s[0]).compare("ROCK")));
463+
bool is_special= (!string(s[0]).substr(0,4).compare("LAKE" )) ||
464+
(!string(s[0]).substr(0,5).compare("WATER" )) ||
465+
(!string(s[0]).substr(0,7).compare("GLACIER" )) ||
466+
(!string(s[0]).substr(0,4).compare("ROCK" )) ||
467+
(!string(s[0]).substr(0,8).compare("PAVEMENT")) ||
468+
(!string(s[0]).substr(0,7).compare("WETLAND" ));
469+
469470
ExitGracefullyIf((nhoriz==0) && (!is_special),
470471
"ParseClassPropertiesFile: only special soil profiles (LAKE,WATER,GLACIER,PAVEMENT,WETLAND, or ROCK) can have zero horizons",BAD_DATA);
471472

‎src/Reservoir.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1235,10 +1235,9 @@ double CReservoir::AdjustFlow(const double& Qadjust, const bool overriding, cons
12351235

12361236
_DAscale=1.0;
12371237
_DAscale_last=1.0;
1238-
return 0.0; //scaling always stops at reservoir or only adjusts inflow
12391238

12401239
double scale=(_Qout+Qadjust)/_Qout;
1241-
if (_Qout=0){scale=1.0;}
1240+
if (_Qout==0.0){scale=1.0;}
12421241

12431242
double va=0.0; //volume added
12441243
double sf=(scale-1.0)/scale;

‎src/StandardOutput.cpp‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -620,12 +620,6 @@ void CModel::WriteOutputFileHeaders(const optStruct &Options)
620620
}
621621
}
622622
ASSIM<<endl;
623-
for (int p = 0; p < _nSubBasins; p++) {
624-
if((_pSubBasins[p]->IsGauged()) && (_pSubBasins[p]->IsEnabled())){
625-
ASSIM<<", ";
626-
}
627-
}
628-
ASSIM<<endl;
629623
ASSIM.close();
630624
}
631625

0 commit comments

Comments
 (0)