Project

General

Profile

round off problem in EquivalenceModel.cxx ยป EquivalenceModel.cxx

bidaud adrien, 04/28/2016 12:00 PM

 
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

    
    (1-1/1)