1
|
#include "EquivalenceModel.hxx"
|
2
|
#include "StringLine.hxx"
|
3
|
#include "CLASSMethod.hxx"
|
4
|
|
5
|
//________________________________________________________________________
|
6
|
EquivalenceModel::EquivalenceModel():CLASSObject()
|
7
|
{
|
8
|
fRelativMassPrecision = 5/10000.; // Mass precision
|
9
|
fMaxInterration = 100; // Max iterration in build fueld algorythum
|
10
|
fFirstGuessFissilContent = 0.02;
|
11
|
freaded = false;
|
12
|
EquivalenceModel::LoadKeyword();
|
13
|
}
|
14
|
//________________________________________________________________________
|
15
|
EquivalenceModel::EquivalenceModel(CLASSLogger* log):CLASSObject(log)
|
16
|
{
|
17
|
fRelativMassPrecision = 5/10000.; // Mass precision
|
18
|
fMaxInterration = 100; // Max iterration in build fueld algorythm
|
19
|
fFirstGuessFissilContent = 0.02;
|
20
|
freaded = false;
|
21
|
EquivalenceModel::LoadKeyword();
|
22
|
|
23
|
}
|
24
|
//________________________________________________________________________
|
25
|
EquivalenceModel::~EquivalenceModel()
|
26
|
{
|
27
|
|
28
|
}
|
29
|
//________________________________________________________________________
|
30
|
void EquivalenceModel::ReadNFO()
|
31
|
{
|
32
|
DBGL
|
33
|
ifstream NFO(fInformationFile.c_str());
|
34
|
|
35
|
if(!NFO)
|
36
|
{
|
37
|
ERROR << "Can't find/open file " << fInformationFile << endl;
|
38
|
exit(0);
|
39
|
}
|
40
|
|
41
|
do
|
42
|
{
|
43
|
string line;
|
44
|
getline(NFO,line);
|
45
|
|
46
|
EquivalenceModel::ReadLine(line);
|
47
|
|
48
|
} while(!NFO.eof());
|
49
|
|
50
|
DBGL
|
51
|
}
|
52
|
//________________________________________________________________________
|
53
|
void EquivalenceModel::ReadLine(string line)
|
54
|
{
|
55
|
DBGL
|
56
|
|
57
|
if (!freaded)
|
58
|
{
|
59
|
int pos = 0;
|
60
|
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
|
61
|
|
62
|
map<string, EQM_MthPtr>::iterator it = fKeyword.find(keyword);
|
63
|
|
64
|
if(it != fKeyword.end())
|
65
|
(this->*(it->second))( line );
|
66
|
|
67
|
freaded = true;
|
68
|
ReadLine(line);
|
69
|
|
70
|
}
|
71
|
|
72
|
freaded = false;
|
73
|
|
74
|
DBGL
|
75
|
}
|
76
|
//________________________________________________________________________
|
77
|
void EquivalenceModel::LoadKeyword()
|
78
|
{
|
79
|
DBGL
|
80
|
fKeyword.insert( pair<string, EQM_MthPtr>( "k_zail", & EquivalenceModel::ReadZAIlimits));
|
81
|
fKeyword.insert( pair<string, EQM_MthPtr>( "k_reactor", & EquivalenceModel::ReadType) );
|
82
|
fKeyword.insert( pair<string, EQM_MthPtr>( "k_fuel", & EquivalenceModel::ReadType) );
|
83
|
fKeyword.insert( pair<string, EQM_MthPtr>( "k_fissil", & EquivalenceModel::ReadFissil) );
|
84
|
fKeyword.insert( pair<string, EQM_MthPtr>( "k_fertil", & EquivalenceModel::ReadFertil) );
|
85
|
fKeyword.insert( pair<string, EQM_MthPtr>( "k_specpower", & EquivalenceModel::ReadSpecificPower));
|
86
|
DBGL
|
87
|
}
|
88
|
//________________________________________________________________________
|
89
|
void EquivalenceModel::PrintInfo()
|
90
|
{
|
91
|
INFO << "Reactor Type : "<< fDBRType << endl;
|
92
|
INFO << "Fuel Type : "<< fDBFType << endl;
|
93
|
INFO << "Specific Power [W/g]: "<< fSpecificPower << endl;
|
94
|
|
95
|
INFO << "Fissile Liste (Z A I) :" << endl;
|
96
|
map<ZAI ,double >::iterator it1;
|
97
|
map<ZAI ,double > fMap1 = fFissileList.GetIsotopicQuantity();
|
98
|
for(it1 = fMap1.begin() ; it1 != fMap1.end() ; it1++)
|
99
|
INFO << (*it1).first.Z() <<" "<< (*it1).first.A() <<" "<< (*it1).first.I() << endl;
|
100
|
|
101
|
INFO<<"Fertile Liste (Z A I Quantity) :"<<endl;
|
102
|
map<ZAI ,double >::iterator it2;
|
103
|
map<ZAI ,double > fMap2 = fFertileList.GetIsotopicQuantity();
|
104
|
for(it2 = fMap2.begin() ; it2 != fMap2.end() ; it2++)
|
105
|
INFO << (*it2).first.Z() <<" "<< (*it2).first.A() << " "<< (*it2).first.I() <<" "<< (*it2).second << endl;
|
106
|
|
107
|
|
108
|
INFO<<"ZAI limits (validity domain)[prop in fresh fuel] (Z A I min max) :"<<endl;
|
109
|
for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++)
|
110
|
{
|
111
|
double ThatZAIMin = Domain_it->second.first;
|
112
|
double ThatZAIMax = Domain_it->second.second;
|
113
|
int Z = Domain_it->first.Z();
|
114
|
int A = Domain_it->first.A();
|
115
|
int I = Domain_it->first.I();
|
116
|
|
117
|
INFO <<ThatZAIMin<<" < ZAI ("<< Z<< " " << A << " " << I<<")"<<" < "<<ThatZAIMax<< endl;
|
118
|
}
|
119
|
|
120
|
}
|
121
|
//________________________________________________________________________
|
122
|
void EquivalenceModel::ReadType(const string &line)
|
123
|
{
|
124
|
DBGL
|
125
|
int pos = 0;
|
126
|
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
|
127
|
if( keyword != "k_fuel" && keyword != "k_reactor" ) // Check the keyword
|
128
|
{
|
129
|
ERROR << " Bad keyword : " << keyword << " Not found !" << endl;
|
130
|
exit(1);
|
131
|
}
|
132
|
if( keyword == "k_fuel" )
|
133
|
fDBFType = StringLine::NextWord(line, pos, ' ');
|
134
|
else if( keyword == "k_reactor" )
|
135
|
fDBRType = StringLine::NextWord(line, pos, ' ');
|
136
|
|
137
|
DBGL
|
138
|
}
|
139
|
//________________________________________________________________________
|
140
|
void EquivalenceModel::ReadZAIlimits(const string &line)
|
141
|
{
|
142
|
DBGL
|
143
|
int pos = 0;
|
144
|
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
|
145
|
if( keyword != "k_zail" ) // Check the keyword
|
146
|
{
|
147
|
ERROR << " Bad keyword : \"k_zail\" not found !" << endl;
|
148
|
exit(1);
|
149
|
}
|
150
|
|
151
|
int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
152
|
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
153
|
int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
154
|
|
155
|
double downLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
|
156
|
double upLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
|
157
|
|
158
|
if (upLimit < downLimit)
|
159
|
{
|
160
|
double tmp = upLimit;
|
161
|
upLimit = downLimit;
|
162
|
downLimit = tmp;
|
163
|
}
|
164
|
fZAILimits.insert(pair<ZAI, pair<double, double> >(ZAI(Z,A,I), pair<double,double>(downLimit, upLimit)));
|
165
|
DBGL
|
166
|
}
|
167
|
//________________________________________________________________________
|
168
|
void EquivalenceModel::ReadFissil(const string &line)
|
169
|
{
|
170
|
DBGL
|
171
|
int pos = 0;
|
172
|
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
|
173
|
if( keyword != "k_fissil" ) // Check the keyword
|
174
|
{
|
175
|
ERROR << " Bad keyword : \"k_fissil\" not found !" << endl;
|
176
|
exit(1);
|
177
|
}
|
178
|
|
179
|
int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
180
|
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
181
|
int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
182
|
|
183
|
fFissileList.Add(Z, A, I, 1.0);
|
184
|
|
185
|
DBGL
|
186
|
}
|
187
|
//________________________________________________________________________
|
188
|
void EquivalenceModel::ReadFertil(const string &line)
|
189
|
{
|
190
|
DBGL
|
191
|
int pos = 0;
|
192
|
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
|
193
|
if( keyword != "k_fertil" ) // Check the keyword
|
194
|
{
|
195
|
ERROR << " Bad keyword : \"k_fertil\" not found !" << endl;
|
196
|
exit(1);
|
197
|
}
|
198
|
|
199
|
int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
200
|
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
201
|
int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
|
202
|
double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
|
203
|
|
204
|
fFertileList.Add(Z, A, I, Q);
|
205
|
|
206
|
|
207
|
DBGL
|
208
|
}
|
209
|
//________________________________________________________________________
|
210
|
void EquivalenceModel::ReadSpecificPower(const string &line)
|
211
|
{
|
212
|
DBGL
|
213
|
int pos = 0;
|
214
|
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
|
215
|
if( keyword != "k_specpower") // Check the keyword
|
216
|
{
|
217
|
ERROR << " Bad keyword : \"k_specpower\" Not found !" << endl;
|
218
|
exit(1);
|
219
|
}
|
220
|
|
221
|
fSpecificPower = atof(StringLine::NextWord(line, pos, ' ').c_str());
|
222
|
|
223
|
DBGL
|
224
|
}
|
225
|
//________________________________________________________________________
|
226
|
double EquivalenceModel::LAMBDA_TOT_FOR(double MassNeeded, vector<IsotopicVector> Stocks, string FisOrFer)
|
227
|
{
|
228
|
double Lambda_tot = 0;
|
229
|
|
230
|
// Calculating total mass of stock once and for all
|
231
|
if( fTotalFissileMassInStocks == 0 || fTotalFertileMassInStocks == 0 )
|
232
|
{
|
233
|
double TotalMassInStocks = 0;
|
234
|
for( int i = 0 ; i<(int)Stocks.size() ; i++ )
|
235
|
TotalMassInStocks += Stocks[i].GetTotalMass() ;
|
236
|
|
237
|
if(FisOrFer == "Fis")
|
238
|
fTotalFissileMassInStocks = TotalMassInStocks * 1e6; // in grams
|
239
|
else
|
240
|
fTotalFertileMassInStocks = TotalMassInStocks * 1e6; // in grams
|
241
|
}
|
242
|
|
243
|
double TotalMassInStocks = 0;
|
244
|
|
245
|
if(FisOrFer == "Fis")
|
246
|
TotalMassInStocks = fTotalFissileMassInStocks;
|
247
|
else
|
248
|
TotalMassInStocks = fTotalFertileMassInStocks;
|
249
|
|
250
|
// If there is not enought matter in stocks construction fails
|
251
|
if( MassNeeded > TotalMassInStocks )
|
252
|
{
|
253
|
WARNING << "Not enought " << FisOrFer << " material to build fuel" << endl;
|
254
|
WARNING << TotalMassInStocks << endl;
|
255
|
return -1;
|
256
|
}
|
257
|
|
258
|
for( int i = 0 ; i<(int)Stocks.size() ; i++ )
|
259
|
{
|
260
|
|
261
|
if( MassNeeded >= (Stocks[i].GetTotalMass()*1e6) )
|
262
|
{
|
263
|
Lambda_tot+= 1;
|
264
|
MassNeeded -= (Stocks[i].GetTotalMass()*1e6);
|
265
|
}
|
266
|
else
|
267
|
{
|
268
|
Lambda_tot+= MassNeeded/(Stocks[i].GetTotalMass()*1e6);
|
269
|
break;
|
270
|
}
|
271
|
}
|
272
|
|
273
|
return Lambda_tot;
|
274
|
}
|
275
|
//________________________________________________________________________
|
276
|
bool EquivalenceModel::Build_Fuel_According_Lambda(vector<double> &lambda,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray, double HMMass, IsotopicVector &Fissile, IsotopicVector &Fertile)
|
277
|
{
|
278
|
//Build the Plutonium vector from stocks
|
279
|
Fissile.Clear();
|
280
|
for( int i = 0; i < (int)FissilArray.size(); i++ )
|
281
|
Fissile += lambda[i] * FissilArray[i];
|
282
|
|
283
|
|
284
|
double AvailablePuMass = Fissile.GetTotalMass() * 1e6; // in grams
|
285
|
|
286
|
// Building complementary Fertile from stocks
|
287
|
double FertilMassNeeded = HMMass - AvailablePuMass;
|
288
|
//adrien
|
289
|
if( FertilMassNeeded < 0)
|
290
|
{
|
291
|
WARNING<<"Be Carefull, FertilMassNeeded negative and set to 0"<<endl;
|
292
|
FertilMassNeeded = 0;
|
293
|
}
|
294
|
double LAMBDA_FERTILE = LAMBDA_TOT_FOR( FertilMassNeeded , FertilArray , "Fer" );
|
295
|
//adrien
|
296
|
|
297
|
SetLambda(lambda, (int)FissilArray.size(), (int)lambda.size()-1, LAMBDA_FERTILE);
|
298
|
|
299
|
int j = -1;
|
300
|
Fertile.Clear();
|
301
|
for(int i = (int)FissilArray.size() ; i < (int)FissilArray.size()+(int)FertilArray.size() ; i++)
|
302
|
{
|
303
|
j++;
|
304
|
Fertile += lambda[i] * FertilArray[j];
|
305
|
}
|
306
|
|
307
|
if( fabs(Fertile.GetTotalMass()*1e6 - FertilMassNeeded) > FertilMassNeeded * 1e-6) // Not enought fertile in stocks
|
308
|
{
|
309
|
WARNING << "Not enought fertile material to build fuel" << endl;
|
310
|
return false;
|
311
|
}
|
312
|
|
313
|
return true;
|
314
|
}
|
315
|
//________________________________________________________________________
|
316
|
vector<double> EquivalenceModel::BuildFuel(double BurnUp, double HMMass, vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
|
317
|
{
|
318
|
|
319
|
DBGL
|
320
|
vector<double> lambda ; // vector of portion of stocks taken (fissile & fertil)
|
321
|
for(int i = 0 ; i < (int)FissilArray.size() + (int)FertilArray.size() ; i++ )
|
322
|
lambda.push_back(0);
|
323
|
|
324
|
/*** Test if there is a stock **/
|
325
|
if( (int)FissilArray.size() == 0 )
|
326
|
{
|
327
|
WARNING << " No fissile stocks available ! Fuel not build" << endl;
|
328
|
SetLambdaToErrorCode(lambda);
|
329
|
return lambda;
|
330
|
}
|
331
|
|
332
|
HMMass *= 1e6; // Unit onversion : tons to gram
|
333
|
|
334
|
/**** Some initializations **/
|
335
|
fTotalFissileMassInStocks = 0;
|
336
|
fTotalFertileMassInStocks = 0;
|
337
|
|
338
|
fActualFissileContent = GetBuildFuelFirstGuess();
|
339
|
|
340
|
IsotopicVector Fertile;
|
341
|
IsotopicVector Fissile;
|
342
|
|
343
|
double AvailablePuMass = 0;
|
344
|
double PuMassNeeded = HMMass * fActualFissileContent;
|
345
|
double WeightPuContent = 0;
|
346
|
int loopCount = 0;
|
347
|
|
348
|
do
|
349
|
{
|
350
|
double LAMBDA_NEEDED = LAMBDA_TOT_FOR(PuMassNeeded,FissilArray,"Fis");
|
351
|
if( LAMBDA_NEEDED == -1 ) // Check if previous lambda was well calculated
|
352
|
{
|
353
|
SetLambdaToErrorCode(lambda);
|
354
|
WARNING << "Not enought fissile material to build fuel" << endl;
|
355
|
return lambda;
|
356
|
}
|
357
|
|
358
|
SetLambda(lambda, 0, FissilArray.size()-1, LAMBDA_NEEDED );
|
359
|
|
360
|
bool succeed = Build_Fuel_According_Lambda(lambda, FissilArray, FertilArray, HMMass, Fissile, Fertile);
|
361
|
|
362
|
if(!succeed)
|
363
|
{
|
364
|
SetLambdaToErrorCode(lambda);
|
365
|
return lambda;
|
366
|
}
|
367
|
|
368
|
AvailablePuMass = Fissile.GetTotalMass() * 1e6; //in grams
|
369
|
|
370
|
if (loopCount > fMaxInterration)
|
371
|
{
|
372
|
ERROR << "Too much iterration in BuildFuel Method !";
|
373
|
ERROR << "Need improvement in fuel fabrication ! Ask for it or D.I.Y. !!" << endl;
|
374
|
exit(1);
|
375
|
}
|
376
|
|
377
|
/* Calcul the quantity of this composition needed to reach the burnup */
|
378
|
double MolarPuContent = GetFissileMolarFraction(Fissile, Fertile, BurnUp);
|
379
|
|
380
|
if( MolarPuContent < 0 || MolarPuContent > 1 )
|
381
|
{ SetLambdaToErrorCode(lambda);
|
382
|
WARNING << "GetFissileMolarFraction return negative or greater than one value";
|
383
|
return lambda;
|
384
|
}
|
385
|
|
386
|
double MeanMolarPu = Fissile.GetMeanMolarMass();
|
387
|
double MeanMolarDepletedU = Fertile.GetMeanMolarMass();
|
388
|
|
389
|
double MeanMolar = MeanMolarPu * MolarPuContent + (1-MolarPuContent) * MeanMolarDepletedU;
|
390
|
|
391
|
|
392
|
WeightPuContent = MolarPuContent * MeanMolarPu / MeanMolar;
|
393
|
fActualFissileContent = MolarPuContent; //fActualFissileContent can be accessed by a derivated EquivalenModel to accelerate GetFissileMolarFraction function (exemple in EQM_MLP_Kinf)
|
394
|
PuMassNeeded = WeightPuContent * HMMass ;
|
395
|
|
396
|
DBGV( "MolarPuContent " << MolarPuContent << " DeltaM " << PuMassNeeded - AvailablePuMass << " g" );
|
397
|
|
398
|
loopCount++;
|
399
|
|
400
|
}while( fabs( PuMassNeeded - AvailablePuMass )/HMMass > fRelativMassPrecision );
|
401
|
|
402
|
(*this).isIVInDomain(Fissile);
|
403
|
|
404
|
DBGV( "Weight percent fissile : " << PuMassNeeded/HMMass );
|
405
|
DBGV( "Lambda vector : " );
|
406
|
|
407
|
for(int i = 0; i < (int)FissilArray.size() + (int)FertilArray.size(); i++ )
|
408
|
DBGV(lambda[i]);
|
409
|
|
410
|
return lambda;
|
411
|
}
|
412
|
//________________________________________________________________________
|
413
|
void EquivalenceModel::SetLambda(vector<double>& lambda ,int FirstStockID, int LastStockID, double LAMBDA_TOT)
|
414
|
{
|
415
|
if( LAMBDA_TOT > LastStockID - FirstStockID + 1 )
|
416
|
{
|
417
|
ERROR << " FATAL ERROR " << endl;
|
418
|
exit(0);
|
419
|
}
|
420
|
|
421
|
for( int i = FirstStockID; i <= LastStockID; i++ ) //set to 0 all non touched value (to be sure)
|
422
|
lambda[i] = 0 ;
|
423
|
|
424
|
int IntegerPart = floor( LAMBDA_TOT );
|
425
|
double DecimalPart = LAMBDA_TOT - IntegerPart;
|
426
|
|
427
|
for( int i = FirstStockID; i < FirstStockID +IntegerPart; i++ )
|
428
|
lambda[i] = 1;
|
429
|
|
430
|
lambda[FirstStockID + IntegerPart] = DecimalPart;
|
431
|
}
|
432
|
|
433
|
//________________________________________________________________________
|
434
|
void EquivalenceModel::SetLambdaToErrorCode(vector<double>& lambda)
|
435
|
{
|
436
|
for(int i = 0 ; i < (int)lambda.size() ;i++ )
|
437
|
lambda[i] = -1;
|
438
|
}
|
439
|
//________________________________________________________________________
|
440
|
bool EquivalenceModel::isIVInDomain(IsotopicVector IV)
|
441
|
{
|
442
|
DBGL
|
443
|
bool IsInDomain = true;
|
444
|
|
445
|
if(fZAILimits.empty())
|
446
|
{
|
447
|
WARNING << "Fresh Fuel variation domain is not set" << endl;
|
448
|
WARNING << "CLASS has no clue if the computed evolution for this fresh fuel is correct" << endl;
|
449
|
WARNING << "Proceed finger crossed !!" << endl;
|
450
|
return true;
|
451
|
}
|
452
|
|
453
|
else
|
454
|
{
|
455
|
IsotopicVector IVNorm = IV /IV.GetSumOfAll();
|
456
|
for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++)
|
457
|
{
|
458
|
double ThatZAIProp = IVNorm.GetIsotopicQuantity()[Domain_it->first] ;
|
459
|
double ThatZAIMin = Domain_it->second.first;
|
460
|
double ThatZAIMax = Domain_it->second.second;
|
461
|
if( (ThatZAIProp > ThatZAIMax) || (ThatZAIProp < ThatZAIMin) )
|
462
|
{
|
463
|
IsInDomain = false;
|
464
|
|
465
|
WARNING << "Fresh fuel out of model range" << endl;
|
466
|
WARNING << "\t AT LEAST this ZAI is accused to be outrange :" << endl;
|
467
|
WARNING << "\t\t" << Domain_it->first.Z() << " " << Domain_it->first.A() << " " << Domain_it->first.I() << endl;
|
468
|
WARNING << "\t\t min = " << ThatZAIMin << " value = " << ThatZAIProp << " max = " << ThatZAIMax << endl;
|
469
|
WARNING << "\t IV accused :" << endl << endl;
|
470
|
WARNING << IVNorm.sPrint() << endl;
|
471
|
break;
|
472
|
}
|
473
|
}
|
474
|
}
|
475
|
DBGL
|
476
|
return IsInDomain;
|
477
|
|
478
|
}
|
479
|
|