diff --git a/source/trunk/include/CLASSHeaders.hxx b/source/trunk/include/CLASSHeaders.hxx
index dbd173fd7b449f83e29dd7e7011f93000f351cfa..2a9e3fe3f17e81ec5320e451092a9a3efacc70b4 100755
--- a/source/trunk/include/CLASSHeaders.hxx
+++ b/source/trunk/include/CLASSHeaders.hxx
@@ -20,6 +20,4 @@
#include "PhysicsModels.hxx"
#include "DecayDataBank.hxx"
-#include "ZAIMass.hxx"
-
#endif
diff --git a/source/trunk/include/EquivalenceModel.hxx b/source/trunk/include/EquivalenceModel.hxx
index 60ecb6b2e51e55a71dbbb27fa2f96b1ec1767136..b005d87133296a071b04ac12b71b866d1edf296c 100644
--- a/source/trunk/include/EquivalenceModel.hxx
+++ b/source/trunk/include/EquivalenceModel.hxx
@@ -53,6 +53,11 @@ class EquivalenceModel : public CLASSObject
virtual ~EquivalenceModel(); //!< Destructor
//@}
+ /*!
+ \name Fuel Construction Method
+ */
+ //@{
+
//{
/// BuildFuel function.
/*!
@@ -70,17 +75,35 @@ class EquivalenceModel : public CLASSObject
int FirstStockID, int LastStockID, double DeltaM,
vector<IsotopicVector> Stocks, double HMMass); //!< Guess a portion of IsotopicVectors to take (dichotomy)
- virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp) = 0; //!< Return the molar fraction of fissile element in the fuel according to the burnup, and a given fuel composition (this is the heart of the equivalence model)
-
- void SetBuildFuelFirstGuess(double FirstGuess){fFirstGuessFissilContent = FirstGuess;} //!<set the initialization value for BuildFuel algorithm
- double GetBuildFuelFirstGuess(){return fFirstGuessFissilContent;} //!<Get the initialization value for BuildFuel algorithm
+ //@}
+ /*!
+ \name Get/Set Method
+ */
+ //@{
IsotopicVector GetFertileList() {return fFertileList;} //!<return the fertile list
IsotopicVector GetFissileList() {return fFissileList;} //!<return the fissile list
+ double GetBuildFuelFirstGuess(){return fFirstGuessFissilContent;} //!<Get the initialization value for BuildFuel algorithm
+ virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp) {return 0;} //!< Return the molar fraction of fissile element in the fuel according to the burnup, and a given fuel composition (this is the heart of the equivalence model)
+
+ double GetRelativMassPrecision() const { return fRelativMassPrecision; } //!< Mass precision
+ int GetMaxInterration() const { return fMaxInterration; } //!< Max iterration in build fueld algorythm
+
+
+
void SetFertileList(IsotopicVector IV) {fFertileList = IV;}//!<set the fertile list
void SetFissileList(IsotopicVector IV) {fFissileList = IV;}//!<set the fissile list
+ void SetBuildFuelFirstGuess(double FirstGuess){fFirstGuessFissilContent = FirstGuess;} //!<set the initialization value for BuildFuel algorithm
+ void SetRelativMassPrecision( double val) { fRelativMassPrecision = val; } //!< Mass precision
+ void SetMaxInterration(int val) { fMaxInterration = val; } //!< Max iterration in build fueld algorythm
+
+
+ //@}
+
+
+
protected :
@@ -93,14 +116,32 @@ class EquivalenceModel : public CLASSObject
double fFirstGuessFissilContent;//!< fissile content for BuildFuel initialization (in weight proportion)
+
+
+
private :
- void SetLambda(vector<double>& lambda ,int FirstStockID, int LastStockID, double LAMBDA_TOT); //!< Set individual lambda according to the LAMBDA_TOT (lambda of all stocks)
- double FindLambdaMax( vector<IsotopicVector> Stocks, double HMMass); //!< Find the maximum LAMBDA_TOT of Stocks (ie lambda to reach HMass)
- double fOld_Lambda_Tot;//!<The old (old iteration) guessed lambda_tot (guessed from GuessLambda)
- double fLambda_max;//!<Value calculated by FindLambdaMax
+ double fOld_Lambda_Tot; //!< The old (old iteration) guessed lambda_tot (guessed from GuessLambda)
+ double fLambda_max; //!< Value calculated by FindLambdaMax
+
+ double fRelativMassPrecision; //!< Mass precision
+ int fMaxInterration; //!< Max iterration in build fueld algorythm
+
+
+ double FindLambdaMax( vector<IsotopicVector> Stocks, double HMMass); //!< Find the maximum LAMBDA_TOT of Stocks (ie lambda to reach HMass)
+ void SetLambda(vector<double>& lambda ,int FirstStockID, int LastStockID, double LAMBDA_TOT); //!< Set individual lambda according to the LAMBDA_TOT (lambda of all stocks)
+
+
};
#endif
+
+
+
+
+
+
+
+
diff --git a/source/trunk/include/IsotopicVector.hxx b/source/trunk/include/IsotopicVector.hxx
index 685e865f78c8b45b0c151857ed4894fa11f41ec8..2e4d9af1a3306ed98eaa7d22f040b453153c1adf 100755
--- a/source/trunk/include/IsotopicVector.hxx
+++ b/source/trunk/include/IsotopicVector.hxx
@@ -150,9 +150,9 @@ public :
protected :
- map<ZAI ,double> fIsotopicQuantity; //!< Isotopic vector composition in atomes Number
+ map<ZAI ,double> fIsotopicQuantity; ///< Isotopic vector composition in atomes Number
- map<ZAI ,double> fIsotopicQuantityNeeded; //!< map where negative value are saved
+ map<ZAI ,double> fIsotopicQuantityNeeded; ///< map where negative value are saved
ClassDef(IsotopicVector,1);
};
diff --git a/source/trunk/include/Scenario.hxx b/source/trunk/include/Scenario.hxx
index 7023ad167b75aada6a45bc6055a28fa4f66ca50d..276a1c6b9240f70504a8d6c8e80c2e874c5179de 100755
--- a/source/trunk/include/Scenario.hxx
+++ b/source/trunk/include/Scenario.hxx
@@ -285,7 +285,7 @@ class Scenario : public CLASSObject
IsotopicVector fWaste; ///< Waste IV
IsotopicVector fTotalStorage; ///< Sum of all IV in Storage IV
- IsotopicVector fOutIncome; ///< OutIncomeIncome IV
+ IsotopicVector fOutIncome; ///< OutIncomeIncome IV
IsotopicVector fTotalCooling; ///< Sum of all IV in Cooling IV
IsotopicVector fFuelFabrication; ///< Sum of all IV in Fabrication IV
IsotopicVector fTotalInReactor; ///< Sum of all IV in Reactor IV
diff --git a/source/trunk/src/EquivalenceModel.cxx b/source/trunk/src/EquivalenceModel.cxx
index 9b0337426fb58722f9d0c5e6f8a490100c36fbbb..b73d2536de6753cd06b6846508b8a1a77d230704 100644
--- a/source/trunk/src/EquivalenceModel.cxx
+++ b/source/trunk/src/EquivalenceModel.cxx
@@ -2,16 +2,19 @@
-
-
-
EquivalenceModel::EquivalenceModel():CLASSObject()
{
+ fRelativMassPrecision = 1/1000.; //Mass precision
+ fMaxInterration = 100; // Max iterration in build fueld algorythum
+ fFirstGuessFissilContent = 0.04;
}
EquivalenceModel::EquivalenceModel(CLASSLogger* log):CLASSObject(log)
{
+ fRelativMassPrecision = 1/1000.; //Mass precision
+ fMaxInterration = 100; // Max iterration in build fueld algorythm
+ fFirstGuessFissilContent = 0.04;
}
@@ -21,7 +24,7 @@ EquivalenceModel::~EquivalenceModel()
}
//________________________________________________________________________
vector<double> EquivalenceModel::BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
-{
+{
DBGL
HMMass*=1e6;//tons to gram
@@ -44,22 +47,36 @@ DBGL
IsotopicVector Fertile;
IsotopicVector Fissile;
- double AvailablePuMass=0;
- double PuMassNeeded=1000; //At present time I have no clue what is the requiered Pu mass, I assume at least 1kg is needed
- double WeightPuContent=0;
- double MassPrecision=100; //Mass precision is 100 grams
+ double AvailablePuMass = 0;
+ double PuMassNeeded = HMMass * 0.01; //At present time I have no clue what is the requiered Pu mass, I assume at least 1kg is needed
+ double WeightPuContent = 0;
- while( fabs(PuMassNeeded - AvailablePuMass) > MassPrecision )
- { //Increase the portion of the stock stokID taken, according to the followings variables
+ int loopCount = 0;
+
+ while( fabs(PuMassNeeded - AvailablePuMass)/HMMass > fRelativMassPrecision )
+ {
+ //Increase the portion of the stock stokID taken, according to the followings variables
double DeltaM=PuMassNeeded-AvailablePuMass;
- GuessLambda( lambda,0,FissilArray.size()-1, DeltaM, FissilArray ,HMMass);
- if(lambda[0]==-1)
+ GuessLambda( lambda, 0, (int)FissilArray.size()-1, DeltaM, FissilArray ,HMMass);
+
+ if(lambda[0]==-1) // Check if previous lambda was well calculated
{
for(int i=0 ; i < (int)FissilArray.size() + (int)FertilArray.size();i++ )
- lambda[i ]= -1;
+ lambda[i]= -1;
+
WARNING<<"Not enought fissile material to build fuel"<<endl;
return lambda;
}
+ else if (loopCount > fMaxInterration)
+ {
+ for(int i=0 ; i < (int)FissilArray.size() + (int)FertilArray.size();i++ )
+ lambda[i]= -1;
+
+ WARNING << "Too much iterration in BuildFuel Method ! Fuel fabrication fail !!" << endl;
+ return lambda;
+ }
+
+
//Build the Plutonium vector from stocks
Fissile.Clear();
@@ -70,7 +87,9 @@ DBGL
//Building complementary Fertile from stocks
double FertilMassNeeded = HMMass - AvailablePuMass;
double LAMBDA_FERTILE = FindLambdaMax(FertilArray, FertilMassNeeded);
- SetLambda(lambda,FissilArray.size(), lambda.size()-1,LAMBDA_FERTILE);
+
+ SetLambda(lambda, (int)FissilArray.size(), (int)lambda.size()-1,LAMBDA_FERTILE);
+
int j=-1;
Fertile.Clear();
for(int i = (int)FissilArray.size() ; i < (int)FissilArray.size()+(int)FertilArray.size() ; i++)
@@ -85,14 +104,19 @@ DBGL
WARNING<<"Not enought fertile material to build fuel"<<endl;
return lambda;
}
+
/*Calcul the quantity of this composition needed to reach the burnup*/
double MolarPuContent = GetFissileMolarFraction(Fissile, Fertile, BurnUp);
double MeanMolarPu = Fissile.MeanMolar();
double MeanMolarDepletedU = Fertile.MeanMolar();
double MeanMolar = MeanMolarPu * MolarPuContent + (1-MolarPuContent)*MeanMolarDepletedU;
+
DBGV("MolarPuContent "<<MolarPuContent);
+ DBGV("DeltaM "<<DeltaM << " g");
+
WeightPuContent = MolarPuContent * MeanMolarPu / MeanMolar ;
PuMassNeeded = WeightPuContent * HMMass ;
+ loopCount++;
}
@@ -155,8 +179,7 @@ void EquivalenceModel::GuessLambda(vector<double>& lambda,int FirstStockID, int
break;
}
}
-
- SetLambda(lambda,FirstStockID,LastStockID,LAMBDA_TOT );
+ SetLambda(lambda,FirstStockID,LastStockID,LAMBDA_TOT );
IsotopicVector test;
int j=0;
@@ -176,15 +199,8 @@ void EquivalenceModel::GuessLambda(vector<double>& lambda,int FirstStockID, int
fOld_Lambda_Tot = LAMBDA_TOT;
LAMBDA_TOT += (fLambda_max - LAMBDA_TOT)/2.;
- if(LAMBDA_TOT > 0.99*fLambda_max) //if we get close to the total of the stocks
- {
- double MasseTot=0;
- for(int i=0;i<(int)Stocks.size();i++)
- MasseTot+=Stocks[i].GetTotalMass()*1e6;
-
- if( (fLambda_max -LAMBDA_TOT )*MasseTot < 5 ) //If it remains less than 5gram in stocks => the stocks are not enought
+ if(fLambda_max -LAMBDA_TOT )/fLambda_max < 0.01) //if we get close to the total of the stocks
lambda[0]=-1;//error code;
- }
//cout<<"DM>0 LAMBDA_TOT "<<LAMBDA_TOT<<endl;
else
SetLambda(lambda,FirstStockID,LastStockID,LAMBDA_TOT );
diff --git a/source/trunk/src/FabricationPlant.cxx b/source/trunk/src/FabricationPlant.cxx
index cca660549e4b25c04a9b23a9e1d6a10e0805fac5..441ea890200726fafdf27643b5f6d83ff42f6070 100644
--- a/source/trunk/src/FabricationPlant.cxx
+++ b/source/trunk/src/FabricationPlant.cxx
@@ -191,7 +191,7 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
pair<CLASSFuel, double > FuelBU = GetParc()->GetReactor()[ReactorId]->GetFuelPlan()->GetFuelAt(t+GetCycleTime()) ;
PhysicsModels FuelType = *FuelBU.first.GetPhysicsModels();
double R_BU = FuelBU.second;
-
+
fFissileList = FuelType.GetEquivalenceModel()->GetFissileList();
BuildFissileArray();
@@ -206,6 +206,7 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
fFertileArray.push_back( fFertileList / fFertileList.GetTotalMass() * R_HM_Mass );
}
+
vector<double> LambdaArray = FuelType.GetEquivalenceModel()->BuildFuel(R_BU, R_HM_Mass, fFissileArray, fFertileArray);
double LambdaSum = 0;
@@ -290,9 +291,11 @@ DBGL
void FabricationPlant::BuildFissileArray()
{
-
+DBGL
+
for(int i = 0; i < (int)fFissileStorage.size(); i++)
{
+
vector<IsotopicVector> IVArray = fFissileStorage[i]->GetIVArray();
for(int j = 0; j < (int)IVArray.size(); j++)
@@ -314,12 +317,13 @@ void FabricationPlant::BuildFissileArray()
}
SortArray(0);
+DBGL
}
void FabricationPlant::BuildFertileArray()
{
-
+DBGL
for(int i = 0; i < (int)fFertileStorage.size(); i++)
{
@@ -341,7 +345,7 @@ void FabricationPlant::BuildFertileArray()
}
SortArray(1);
-
+DBGL
}
void FabricationPlant::SortArray(int i)
diff --git a/source/trunk/src/Scenario.cxx b/source/trunk/src/Scenario.cxx
index 28b7b6e8a81b09b218fa019d2f9847ac4cd41abf..7c6b1c327ffc6f51fd2ff49348c485ab67ec065d 100755
--- a/source/trunk/src/Scenario.cxx
+++ b/source/trunk/src/Scenario.cxx
@@ -734,8 +734,6 @@ void Scenario::UpdateParc()
fIVTotal = fWaste + fTotalStorage + fTotalCooling + fFuelFabrication + fTotalInReactor;
fIVInCycleTotal = fTotalStorage + fTotalCooling + fFuelFabrication + fTotalInReactor;
-
-
}