From 2492f56c5344d58a803cac7f5f5ca6997d95de78 Mon Sep 17 00:00:00 2001
From: Baptiste LENIAU <baptiste.leniau@subatech.in2p3.fr>
Date: Tue, 19 Apr 2016 14:19:04 +0000
Subject: [PATCH] commit all change from funny courtin
git-svn-id: svn+ssh://svn.in2p3.fr/class@813 0e7d625b-0364-4367-a6be-d5be4a48d228
---
.../Model/Equivalence/EQM_FBR_MLP_Keff.cxx | 43 +-
.../Model/Equivalence/EQM_FBR_MLP_Keff.hxx | 19 +-
.../SL3/Model/Equivalence/EQM_MLP_Kinf.cxx | 185 ++--
.../SL3/Model/Equivalence/EQM_MLP_Kinf.hxx | 23 +-
.../Model/Equivalence/EQM_MLP_PWR_MOxEUS.cxx | 659 ++++++++++++++
.../Model/Equivalence/EQM_MLP_PWR_MOxEUS.hxx | 67 ++
.../Equivalence/EQM_PWR_FixedContent.cxx | 554 ++++++++++++
.../Equivalence/EQM_PWR_FixedContent.hxx | 75 ++
.../SL3/Model/Equivalence/EQM_PWR_LIN_MOX.cxx | 6 +-
.../SL3/Model/Equivalence/EQM_PWR_LIN_MOX.hxx | 2 +-
.../SL3/Model/Equivalence/EQM_PWR_MLP_MOX.cxx | 116 +--
.../SL3/Model/Equivalence/EQM_PWR_MLP_MOX.hxx | 15 +-
.../Model/Equivalence/EQM_PWR_MLP_MOX_Am.cxx | 17 +-
.../Model/Equivalence/EQM_PWR_MLP_MOX_Am.hxx | 15 +-
.../SL3/Model/Equivalence/EQM_PWR_POL_UO2.cxx | 9 +-
.../SL3/Model/Equivalence/EQM_PWR_POL_UO2.hxx | 2 +-
.../Model/Equivalence/EQM_PWR_QUAD_MOX.cxx | 13 +-
.../Model/Equivalence/EQM_PWR_QUAD_MOX.hxx | 2 +-
.../branches/SL3/Model/Irradiation/IM_RK4.cxx | 41 +-
source/branches/SL3/src/EquivalenceModel.cxx | 492 ++++++-----
source/branches/SL3/src/EvolutionData.cxx | 81 --
source/branches/SL3/src/FabricationPlant.cxx | 809 +++++++-----------
source/branches/SL3/src/Reactor.cxx | 83 +-
source/branches/SL3/src/Scenario.cxx | 10 +-
source/branches/SL3/src/XSModel.cxx | 7 +-
25 files changed, 2330 insertions(+), 1015 deletions(-)
create mode 100755 source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.cxx
create mode 100755 source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.hxx
create mode 100755 source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.cxx
create mode 100755 source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.hxx
diff --git a/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.cxx b/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.cxx
index c6998f89e..8c5a5926f 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.cxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.cxx
@@ -68,14 +68,6 @@ EQM_FBR_MLP_Keff::EQM_FBR_MLP_Keff(string TMVAWeightPath, double keff_target, st
INFO << "\tThis model is based on the prediction of keff at a specific time" << endl;
INFO << "\tThe TMVA (weight | information) files are :" << endl;
INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fInformationFile << " )" << endl;
- INFO << "Maximal fissile content (molar proportion) : "<<fMaximalContent<<endl;
- EquivalenceModel::PrintInfo();
-
- if(fMapOfTMVAVariableNames.empty() || fFertileList.GetIsotopicQuantity().empty() || fFissileList.GetIsotopicQuantity().empty())
- {
- ERROR<<"Missing information file in : "<<fInformationFile<<endl;
- exit(1);
- }
DBGL
}
@@ -114,16 +106,10 @@ EQM_FBR_MLP_Keff::EQM_FBR_MLP_Keff(CLASSLogger* log, string TMVAWeightPath, doub
INFO << "\tThis model is based on the prediction of keff at a specific time" << endl;
INFO << "\tThe TMVA (weight | information) files are :" << endl;
INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fInformationFile << " )" << endl;
- EquivalenceModel::PrintInfo();
-
- if(fMapOfTMVAVariableNames.empty() || fFertileList.GetIsotopicQuantity().empty() || fFissileList.GetIsotopicQuantity().empty())
- {
- ERROR<<"Missing information file in : "<<fInformationFile<<endl;
- exit(1);
- }
DBGL
}
+
//________________________________________________________________________
TTree* EQM_FBR_MLP_Keff::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime)
{
@@ -232,6 +218,7 @@ double EQM_FBR_MLP_Keff::ExecuteTMVA(TTree* InputTree, bool IsTimeDependent)
DBGL
return (double)val; //return k_{eff}(t = Time)
}
+
//________________________________________________________________________
void EQM_FBR_MLP_Keff::LoadKeyword()
{
@@ -243,11 +230,11 @@ void EQM_FBR_MLP_Keff::LoadKeyword()
DBGL
}
+
//________________________________________________________________________
void EQM_FBR_MLP_Keff::ReadZAIName(const string &line)
{
DBGL
-
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
if( keyword != "k_zainame" ) // Check the keyword
@@ -260,12 +247,12 @@ void EQM_FBR_MLP_Keff::ReadZAIName(const string &line)
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- string name = StringLine::NextWord(line, pos, ' ');
+ fFissileList.Add(Z, A, I, 1.0);
+ fStreamList.push_back(fFissileList);
- fMapOfTMVAVariableNames.insert( pair<ZAI,string>( ZAI(Z, A, I), name ) );
-
- DBGL
+ DBGL
}
+
//________________________________________________________________________
void EQM_FBR_MLP_Keff::ReadMaxFisContent(const string &line)
{
@@ -282,6 +269,7 @@ void EQM_FBR_MLP_Keff::ReadMaxFisContent(const string &line)
DBGL
}
+
//________________________________________________________________________
void EQM_FBR_MLP_Keff::ReadLine(string line)
{
@@ -297,14 +285,17 @@ void EQM_FBR_MLP_Keff::ReadLine(string line)
DBGL
}
+
//________________________________________________________________________
-double EQM_FBR_MLP_Keff::GetFissileMolarFraction(IsotopicVector Fissile,IsotopicVector Fertile,double TargetBU)
+map < string , double> EQM_FBR_MLP_Keff::GetMolarFraction(vector <IsotopicVector> IVStream,double TargetBU)
{
DBGL
if(TargetBU != 0)
WARNING << "The third arguement : Burnup has no effect here.";
+ IsotopicVector Fissile = IVStream["Fissile"];
+ IsotopicVector Fertile = IVStream["Fertile"];
//initialization
double FissileContent = GetActualFissileContent();
@@ -347,7 +338,13 @@ double EQM_FBR_MLP_Keff::GetFissileMolarFraction(IsotopicVector Fissile,Isotopic
}while(fabs(fTargetKeff-PredictedKeff)>Precision);
DBGV( "Predicted keff " << PredictedKeff << " FissileContent " << FissileContent << endl);
- return FissileContent;
+
+ map < string , double> MolarFraction;
+ MolarFraction["Fissile"] = FissileContent;
+ MolarFraction["Fertile"] = 1.- FissileContent;
+
+ return MolarFraction; //return Molar content of each component in the fuel
+
}
-//________________________________________________________________________
\ No newline at end of file
+//________________________________________________________________________
diff --git a/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.hxx b/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.hxx
index 9a12a9d59..133f184da 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.hxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_FBR_MLP_Keff.hxx
@@ -80,7 +80,7 @@ class EQM_FBR_MLP_Keff : public EquivalenceModel
\param Fertil : The composition of the Fertil matter
\param BurnUp : Maximum achievable burn up envisaged
*/
- virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp = 0);
+ virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream, double BurnUp = 0);
//}
/*!
@@ -137,7 +137,7 @@ class EQM_FBR_MLP_Keff : public EquivalenceModel
//{
/// ReadMaxFisContent : read a guessed (very overestimated) maximum fissile content (purpose : algorithm initialization)
/*!
- \param line : line suppossed to contain the ZAI name starts with "k_maxfiscontent" keyword
+ \param line : line suppossed to contain the ZAI name starts with "k_zainame" keyword
*/
void ReadMaxFisContent(const string &line);
//}
@@ -166,6 +166,7 @@ class EQM_FBR_MLP_Keff : public EquivalenceModel
map<ZAI,string> fMapOfTMVAVariableNames;//!< List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step
vector<double> fMLP_Time; //!< Time (in seconds) when the MLP(t) = keff(t) has been trained.
+ double fSpecificPower; //!< The specific power in W/gHM (HM: heavy Metal)
double fMaximalContent; //!< The approx. maximum fissile content reachable by the MLP model
@@ -178,18 +179,14 @@ class EQM_FBR_MLP_Keff : public EquivalenceModel
double fTargetKeff; //!< Use for Varying Fissile content to reach fTargetKeff at time used in the MLP Training
- /*!
- \name keff prediction methods & keff averaging
- */
- //@{
+
+ double GetKeffAtFixedTime(IsotopicVector FreshFuel){return ExecuteTMVA( CreateTMVAInputTree(FreshFuel,-1), false );} //!<time independant since the MLP is trained for 1 time
- double GetKeffAtFixedTime(IsotopicVector FreshFuel){TTree* Input = CreateTMVAInputTree(FreshFuel,-1); double Keff = ExecuteTMVA( Input, false ); delete Input; return Keff;} //!<time independant since the MLP is trained for 1 time
+ TGraph* BuildKeffGraph(IsotopicVector FreshFuel);
+ TGraph* BuildAverageKeffGraph(TGraph* GRAPH_KEFF);
- TGraph* BuildKeffGraph(IsotopicVector FreshFuel);
- TGraph* BuildAverageKeffGraph(TGraph* GRAPH_KEFF);
- double GetKeffAt(TGraph* GRAPH_KEFF, int Step);
+ double GetKeffAt(TGraph* GRAPH_KEFF, int Step);
- //@}
diff --git a/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.cxx b/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.cxx
index 24a8eb53f..3b35edf62 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.cxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.cxx
@@ -31,17 +31,17 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(string WeightPathAlpha0, string WeightPathAlpha1, str
fTMVAWeightPath.push_back(WeightPathAlpha1);
fTMVAWeightPath.push_back(WeightPathAlpha2);
- fMaximalBU = 0;
- fMaximalContent = 0;
- fInformationFile = InformationFile;
+ fMaximalBU = 0;
+ fMaximalContent = 0;
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold ;
+
+ SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+
+ fInformationFile = InformationFile;
LoadKeyword();
ReadNFO();//Getting information from fInformationFile
- fNumberOfBatch = NumOfBatch;
- fKThreshold = CriticalityThreshold ;
- SetBurnUpPrecision(0.005);//1 % of the targeted burnup
- SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
-
INFO << "__An equivalence model has been define__" << endl;
INFO << "\tThis model is based on the prediction of kinf" << endl;
INFO << "\tThe TMVA weight files are :" << endl;
@@ -54,7 +54,23 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(string WeightPathAlpha0, string WeightPathAlpha1, str
INFO << "Maximal burnup (GWd/tHM) : " << fMaximalBU << endl;
EquivalenceModel::PrintInfo();
- if(fMapOfTMVAVariableNames.empty() || fFertileList.GetIsotopicQuantity().empty() || fFissileList.GetIsotopicQuantity().empty() || fMaximalBU == 0 || fMaximalContent == 0 )
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ bool EmptyList = false;
+ bool FirstContentNULL = false;
+
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ if(fStreamList[(* it_s_IV).first].GetIsotopicQuantity().empty()){EmptyList=true;}
+ }
+
+ for( it_s_D = fFirstGuessContent .begin(); it_s_D != fFirstGuessContent.end(); it_s_D++)
+ {
+ if(fFirstGuessContent[(* it_s_D).first] ==0){FirstContentNULL=true;}
+ }
+
+
+ if(fMapOfTMVAVariableNames.empty() || EmptyList==true || FirstContentNULL==true || fMaximalBU == 0 || fMaximalContent == 0 )
{
ERROR<<"Missing information file in : "<<fInformationFile<<endl;
exit(1);
@@ -69,17 +85,17 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(CLASSLogger* log, string WeightPathAlpha0, string Wei
fTMVAWeightPath.push_back(WeightPathAlpha1);
fTMVAWeightPath.push_back(WeightPathAlpha2);
- fMaximalBU = 0;
- fMaximalContent = 0;
- fInformationFile = InformationFile;
+ fMaximalBU = 0;
+ fMaximalContent = 0;
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold ;
+
+ SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+
+ fInformationFile = InformationFile;
LoadKeyword();
ReadNFO();//Getting information from fInformationFile
- fNumberOfBatch = NumOfBatch;
- fKThreshold = CriticalityThreshold ;
- SetBurnUpPrecision(0.005);//1 % of the targeted burnup
- SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
-
INFO << "__An equivalence model has been define__" << endl;
INFO << "\tThis model is based on the prediction of kinf" << endl;
INFO << "\tThe TMVA weight files are :" << endl;
@@ -92,7 +108,23 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(CLASSLogger* log, string WeightPathAlpha0, string Wei
INFO << "Maximal burnup (GWd/tHM) : " << fMaximalBU << endl;
EquivalenceModel::PrintInfo();
- if(fMapOfTMVAVariableNames.empty() || fFertileList.GetIsotopicQuantity().empty() || fFissileList.GetIsotopicQuantity().empty() || fMaximalBU == 0 || fMaximalContent == 0 )
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ bool EmptyList = false;
+ bool FirstContentNULL = false;
+
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ if(fStreamList[(* it_s_IV).first].GetIsotopicQuantity().empty()){EmptyList=true;}
+ }
+
+ for( it_s_D = fFirstGuessContent .begin(); it_s_D != fFirstGuessContent.end(); it_s_D++)
+ {
+ if(fFirstGuessContent[(* it_s_D).first] ==0){FirstContentNULL=true;}
+ }
+
+
+ if(fMapOfTMVAVariableNames.empty() || EmptyList==true || FirstContentNULL==true || fMaximalBU == 0 || fMaximalContent == 0 )
{
ERROR<<"Missing information file in : "<<fInformationFile<<endl;
exit(1);
@@ -110,18 +142,18 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(string TMVAWeightPath, int NumOfBatch, string Inform
if(InformationFile == "")
InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
- fMaximalBU = 0;
- fMaximalContent = 0;
- fInformationFile = InformationFile;
+ fMaximalBU = 0;
+ fMaximalContent = 0;
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold ;
+
+ SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+ SetPCMPrecision(10);
+
+ fInformationFile = InformationFile;
LoadKeyword();
ReadNFO();//Getting information from fInformationFile
- fNumberOfBatch = NumOfBatch;
- fKThreshold = CriticalityThreshold ;
- SetBurnUpPrecision(0.005);//1 % of the targeted burnup
- SetPCMprecision(10);
- SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
-
INFO << "__An equivalence model has been define__" << endl;
INFO << "\tThis model is based on the prediction of kinf" << endl;
INFO << "\tThe TMVA (weight | information) files are :" << endl;
@@ -130,7 +162,23 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(string TMVAWeightPath, int NumOfBatch, string Inform
INFO << "Maximal burnup (GWd/tHM) : " << fMaximalBU << endl;
EquivalenceModel::PrintInfo();
- if(fMapOfTMVAVariableNames.empty() || fFertileList.GetIsotopicQuantity().empty() || fFissileList.GetIsotopicQuantity().empty() || fMaximalBU == 0 || fMaximalContent == 0 )
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ bool EmptyList = false;
+ bool FirstContentNULL = false;
+
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ if(fStreamList[(* it_s_IV).first].GetIsotopicQuantity().empty()){EmptyList=true;}
+ }
+
+ for( it_s_D = fFirstGuessContent .begin(); it_s_D != fFirstGuessContent.end(); it_s_D++)
+ {
+ if(fFirstGuessContent[(* it_s_D).first] ==0){FirstContentNULL=true;}
+ }
+
+
+ if(fMapOfTMVAVariableNames.empty() || EmptyList==true || FirstContentNULL==true || fMaximalBU == 0 || fMaximalContent == 0 )
{
ERROR<<"Missing information file in : "<<fInformationFile<<endl;
exit(1);
@@ -147,17 +195,16 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(CLASSLogger* log, string TMVAWeightPath, int NumOfBa
if(InformationFile == "")
InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
- fMaximalBU = 0;
- fMaximalContent = 0;
- fInformationFile = InformationFile;
+ fMaximalBU = 0;
+ fMaximalContent = 0;
+ fInformationFile = InformationFile;
LoadKeyword();
ReadNFO();//Getting information from fMLPInformationFile
- fNumberOfBatch = NumOfBatch;
- fKThreshold = CriticalityThreshold ;
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold ;
SetBurnUpPrecision(0.005);//1 % of the targeted burnup
- SetPCMprecision(10);
- SetBuildFuelFirstGuess(0.04);//First fissile content guess for the EquivalenceModel::BuildFuel algorithm
+ SetPCMPrecision(10);
INFO << "__An equivalence model has been define__" << endl;
INFO << "\tThis model is based on the prediction of kinf" << endl;
@@ -166,10 +213,25 @@ EQM_MLP_Kinf::EQM_MLP_Kinf(CLASSLogger* log, string TMVAWeightPath, int NumOfBa
INFO << "Maximal fissile content (molar proportion) : " << fMaximalContent << endl;
INFO << "Maximal burnup (GWd/tHM) : " << fMaximalBU << endl;
-
EquivalenceModel::PrintInfo();
- if(fMapOfTMVAVariableNames.empty() || fFertileList.GetIsotopicQuantity().empty() || fFissileList.GetIsotopicQuantity().empty() || fMaximalBU == 0 || fMaximalContent == 0 )
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ bool EmptyList = false;
+ bool FirstContentNULL = false;
+
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ if(fStreamList[(* it_s_IV).first].GetIsotopicQuantity().empty()){EmptyList=true;}
+ }
+
+ for( it_s_D = fFirstGuessContent .begin(); it_s_D != fFirstGuessContent.end(); it_s_D++)
+ {
+ if(fFirstGuessContent[(* it_s_D).first] ==0){FirstContentNULL=true;}
+ }
+
+
+ if(fMapOfTMVAVariableNames.empty() || EmptyList==true || FirstContentNULL==true || fMaximalBU == 0 || fMaximalContent == 0 )
{
ERROR<<"Missing information file in : "<<fInformationFile<<endl;
exit(1);
@@ -202,7 +264,7 @@ void EQM_MLP_Kinf::ReadZAIName(const string &line)
int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
string name = StringLine::NextWord(line, pos, ' ');
@@ -384,7 +446,7 @@ double EQM_MLP_Kinf::GetMaximumBurnUp_MLP(IsotopicVector TheFuel, double TargetB
{
if(count > MaximumLoopCount )
{
- ERROR << "CRITICAL ! Can't manage to predict burnup\nHint : Try to increase the precision on k effective using :\n YourEQM_MLP_Kinf->SetPCMprecision(pcm); with pcm the precision in pcm (default 10) REDUCE IT\n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
+ ERROR << "CRITICAL ! Can't manage to predict burnup\nHint : Try to increase the precision on k effective using :\n YourEQM_MLP_Kinf->SetPCMPrecision(pcm); with pcm the precision in pcm (default 10) REDUCE IT\n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
exit(1);
}
@@ -414,12 +476,12 @@ double EQM_MLP_Kinf::GetMaximumBurnUp_MLP(IsotopicVector TheFuel, double TargetB
delete InputTree;
}
k_av/= fNumberOfBatch;
-
+ //cout<<SecondToBurnup(TheFinalTime)<<" ";
OldPredictedk_av = k_av;
count++;
- } while( fabs(OldPredictedk_av-fKThreshold) > GetPCMprecision() ) ;
-
+ } while( fabs(OldPredictedk_av-fKThreshold) > GetPCMPrecision() ) ;
+ //cout<<endl;
return SecondToBurnup(TheFinalTime);
}
//________________________________________________________________________
@@ -517,20 +579,25 @@ double EQM_MLP_Kinf::GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU)
return TheBurnUp;
}
//________________________________________________________________________
-double EQM_MLP_Kinf::GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double TargetBU)
+map < string , double> EQM_MLP_Kinf::GetMolarFraction(map <string , IsotopicVector > IVStream, double TargetBU)
{
-
//initialization
- double FissileContent = GetActualFissileContent();
- double OldFissileContentMinus = 0;
- double OldFissileContentPlus = fMaximalContent;
- double PredictedBU = 0 ;
- IsotopicVector FreshFuel = (1-FissileContent)*(Fertil/Fertil.GetSumOfAll()) + FissileContent*(Fissil/Fissil.GetSumOfAll());
- double OldPredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
-
- double Precision = GetBurnUpPrecision()*TargetBU; //1 % of the targeted burnup
- int count = 0;
- int MaximumLoopCount = 500;
+ IsotopicVector Fissil = IVStream["Fissile"];
+ IsotopicVector Fertil = IVStream["Fertile"];
+
+ map < string , double> MolarFraction ;
+
+ double FissileContent = fActualMolarContentInFuel["Fissile"];
+ double OldFissileContentMinus = 0;
+ double OldFissileContentPlus = fMaximalContent;
+
+ double PredictedBU = 0 ;
+ IsotopicVector FreshFuel = (1-FissileContent)*(Fertil/Fertil.GetSumOfAll()) + FissileContent*(Fissil/Fissil.GetSumOfAll());
+ double OldPredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+
+ double Precision = GetBurnUpPrecision()*TargetBU; //1 % of the targeted burnup
+ int count = 0;
+ int MaximumLoopCount = 500;
do
{
if(count > MaximumLoopCount )
@@ -557,15 +624,17 @@ double EQM_MLP_Kinf::GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVecto
IsotopicVector FreshFuel = (1-FissileContent)*(Fertil/Fertil.GetSumOfAll()) + FissileContent*(Fissil/Fissil.GetSumOfAll());
PredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
-
+
OldPredictedBU = PredictedBU;
count ++;
}while(fabs(TargetBU-PredictedBU)>Precision);
-
-
DBGV("Predicted BU " << PredictedBU << " FissileContent " << FissileContent);
- return FissileContent;
+
+ MolarFraction["Fissile"] = FissileContent;
+ MolarFraction["Fertile"] = 1 - FissileContent;
+
+ return MolarFraction;
}
//________________________________________________________________________
\ No newline at end of file
diff --git a/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.hxx b/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.hxx
index 48e0296bb..c793b4270 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.hxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_MLP_Kinf.hxx
@@ -37,7 +37,7 @@ in non simulated devices such as control rods and mixing grid.
class EQM_MLP_Kinf : public EquivalenceModel
{
- public:
+ public :
/*!
\name Constructor
*/
@@ -94,7 +94,7 @@ class EQM_MLP_Kinf : public EquivalenceModel
\param InformationFile : Total path to the file containing time steps, fissile and ferile list (ante and post fabrication time cooling). Default is the same total path as TMVAWeightPath but extension is replaced by .nfo
\param CriticalityThreshold : Threshold for the @f$k_{\infty}@f$ (see detailed description)
*/
- EQM_MLP_Kinf(CLASSLogger* log, string TMVAWeightPath,int NumOfBatch, string InformationFile = "", double CriticalityThreshold = 1.01);
+ EQM_MLP_Kinf(CLASSLogger* log, string TMVAWeightPath, int NumOfBatch, string InformationFile = "", double CriticalityThreshold = 1.01);
//}
//@}
@@ -105,7 +105,7 @@ class EQM_MLP_Kinf : public EquivalenceModel
\param Fertil : The composition of the Fertil matter
\param BurnUp : Maximum achievable burn up envisaged
*/
- virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);
+ map < string , double> GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp);
//}
/*!
@@ -114,9 +114,9 @@ class EQM_MLP_Kinf : public EquivalenceModel
//@{
void SetBurnUpPrecision(double prop){fBurnUpPrecision = prop;} //!< Set the precision on Burnup : proportion of the targeted burnup
- void SetPCMprecision(double pcm){fPCMprecision = pcm;} //!< Set the precision on @f$\langle k \rangle@f$ prediction [pcm]. Neural network predictor constructors
+ void SetPCMPrecision(double pcm){fPCMprecision = pcm;} //!< Set the precision on @f$\langle k \rangle@f$ prediction [pcm]. Neural network predictor constructors
double GetBurnUpPrecision(){return fBurnUpPrecision;}//!< Get the precision on Burnup : proportion of the targeted burnup
- double GetPCMprecision(){return fPCMprecision/1e5;}//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors
+ double GetPCMPrecision(){return fPCMprecision/1e5;}//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors
double GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU);//!<Get the maximum reachable burnup according the freshfuel composition
void GetModelInformation();//!<Read the fMLPInformationFile and fill containers and variables
@@ -185,6 +185,9 @@ class EQM_MLP_Kinf : public EquivalenceModel
void ReadLine(string line);
//}
+ double GetMaximumBurnUp_MLP(IsotopicVector TheFuel, double TargetBU);
+ double GetMaximumBurnUp_Pol2(IsotopicVector TheFuel, double TargetBU);
+
//@}
private :
@@ -201,17 +204,15 @@ class EQM_MLP_Kinf : public EquivalenceModel
*/
//@{
- double GetMaximumBurnUp_MLP(IsotopicVector TheFuel, double TargetBU);
- double GetMaximumBurnUp_Pol2(IsotopicVector TheFuel, double TargetBU);
//@}
int fNumberOfBatch; //!< The number of batches for the loading plan
- double fKThreshold; //!< The @f$k_{Threshold}@f$
+ double fKThreshold; //!< The @f$k_{Threshold}@f$
double fMaximalBU; //!< The approx. maximum burnup reachable by the MLP model
- double fMaximalContent; //!< The approx. maximum fissile content reachable by the MLP model
- double fBurnUpPrecision; //!< precision on Burnup
- double fPCMprecision; //!< precision on @f$\langle k \rangle@f$ prediction [pcm]
+ double fMaximalContent; //!< The approx. maximum fissile content reachable by the MLP model
+ double fBurnUpPrecision; //!< precision on Burnup
+ double fPCMprecision; //!< precision on @f$\langle k \rangle@f$ prediction [pcm]
};
diff --git a/source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.cxx b/source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.cxx
new file mode 100755
index 000000000..f6dae0776
--- /dev/null
+++ b/source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.cxx
@@ -0,0 +1,659 @@
+#include "EquivalenceModel.hxx"
+#include "EQM_MLP_PWR_MOxEUS.hxx"
+#include "CLASSLogger.hxx"
+#include "StringLine.hxx"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+
+#include "TSystem.h"
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+
+//________________________________________________________________________
+//
+// EQM_MLP_PWR_MOxEUS
+//
+// Equivalenve Model based on multi layer perceptron from TMVA (root cern)
+// For PWR MOxEUS use
+// -> MOxEUS building
+// - if PuMassContent calculated < Maximal Pu value => MOx fuel
+// - if PuMassContent calculated > Maximal Pu value => Uranium enrichment adjusted to reach BU
+//
+//________________________________________________________________________
+
+EQM_MLP_PWR_MOxEUS::EQM_MLP_PWR_MOxEUS(string TMVAWeightPath, int NumOfBatch, double CriticalityThreshold, double MaximalPuMassContent):EquivalenceModel(new CLASSLogger("EQM_MLP_PWR_MOxEUS.log"))
+{
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold ;
+ fMaximalPuMassContent = MaximalPuMassContent;
+
+ fTMVAWeightPath.push_back(TMVAWeightPath);
+
+ SetMinimalU5Enrichment(0.003);
+
+ ZAI U8(92,238,0);
+ ZAI U5(92,235,0);
+ ZAI Pu8(94,238,0);
+ ZAI Pu9(94,239,0);
+ ZAI Pu0(94,240,0);
+ ZAI Pu1(94,241,0);
+ ZAI Pu2(94,242,0);
+
+ fStreamList["PuList"] = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+ fStreamList["EnrichmentList"] = U5*1;
+ fStreamList["FertileList"] = U8*1;
+
+ fFirstGuessContent["PuList"] = 0.02;
+ fFirstGuessContent["EnrichmentList"] = (1-fFirstGuessContent["PuList"])*GetMinimalU5Enrichment();
+ fFirstGuessContent["FertileList"] = 1 - fFirstGuessContent["PuList"] - fFirstGuessContent["EnrichmentList"];
+
+ fSpecificPower = 34.24;
+ fMaximalPuMolarContent = 0.50;
+ fMaximalBU = 100;
+
+ SetRelativMassPrecision(0.0008);
+ SetMaximalU5Enrichment(0.10);
+ SetBurnUpPrecision(0.008);//1 % of the targeted burnup
+ SetPCMPrecision(10);
+
+ INFO<<"__An equivalence model of PWR MOxEUS has been define__"<<endl;
+ INFO<<"\tThis model is based on a multi layer perceptron"<<endl;
+ INFO<<"\t\tThe TMVA weight file is :"<<endl;
+ INFO<<"\t\t\t"<<fTMVAWeightPath[0]<<endl;
+
+}
+
+//________________________________________________________________________
+EQM_MLP_PWR_MOxEUS::EQM_MLP_PWR_MOxEUS(CLASSLogger* log, string TMVAWeightPath, int NumOfBatch, double CriticalityThreshold, double MaximalPuMassContent):EquivalenceModel(log)
+{
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold ;
+ fMaximalPuMassContent = MaximalPuMassContent;
+
+ fTMVAWeightPath.push_back(TMVAWeightPath);
+
+ SetMinimalU5Enrichment(0.003);
+
+ ZAI U8(92,238,0);
+ ZAI U5(92,235,0);
+ ZAI Pu8(94,238,0);
+ ZAI Pu9(94,239,0);
+ ZAI Pu0(94,240,0);
+ ZAI Pu1(94,241,0);
+ ZAI Pu2(94,242,0);
+
+ fStreamList["PuList"] = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+ fStreamList["EnrichmentList"] = U5*1;
+ fStreamList["FertileList"] = U8*1;
+
+ fFirstGuessContent["PuList"] = 0.02;
+ fFirstGuessContent["EnrichmentList"] = (1-fFirstGuessContent["PuList"])*GetMinimalU5Enrichment();
+ fFirstGuessContent["FertileList"] = 1 - fFirstGuessContent["PuList"] - fFirstGuessContent["EnrichmentList"];
+
+ fSpecificPower = 34.24;
+ fMaximalPuMolarContent = 1.0;
+ fMaximalBU = 100;
+
+ SetRelativMassPrecision(0.0008);
+ SetMaximalU5Enrichment(0.10);
+ SetBurnUpPrecision(0.008);//1 % of the targeted burnup
+ SetPCMPrecision(10);
+
+ INFO<<"__An equivalence model of PWR MOxEUS has been define__"<<endl;
+ INFO<<"\tThis model is based on a multi layer perceptron"<<endl;
+ INFO<<"\t\tThe TMVA weight file is :"<<endl;
+ INFO<<"\t\t\t"<<fTMVAWeightPath[0]<<endl;
+
+}
+
+//________________________________________________________________________
+map <string , vector<double> > EQM_MLP_PWR_MOxEUS::BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray)
+{
+
+ map <string , vector<double> > lambda ; // map containing name of the list and associated vector of proportions taken from stocks
+
+ //Iterators declaration
+ map < string , vector <IsotopicVector> >::iterator it_s_vIV;
+ map < string , vector <double> >::iterator it_s_vD;
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ map < string , bool >::iterator it_s_B;
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i<StreamArray[(*it_s_vIV).first].size(); i++)
+ {
+ lambda[(*it_s_vIV).first].push_back(0);
+ }
+ }
+
+ /*** Test if there is stocks **/
+ bool BreakReturnLambda = false;
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ if(StreamArray[(*it_s_vIV).first].size() == 0)
+ {
+ WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<". Fuel not build." << endl;
+ SetLambdaToErrorCode(lambda[(*it_s_vIV).first]);
+ BreakReturnLambda = true;
+ }
+ }
+ if(BreakReturnLambda) { return lambda;}
+ HMMass *= 1e6; //Unit onversion : tons to gram
+
+ /**** Some initializations **/
+
+ StocksTotalMassCalculation(StreamArray);
+
+ fActualMassContentInFuel = GetBuildFuelFirstGuess();
+ int loopCount = 0;
+ bool FuelBuiltCorrectly = false ;
+
+ map <string , IsotopicVector > IVStream;
+ map <string , double > MaterialMolarContent;
+ map <string , double > MaterialMassContent;
+ map <string , double > MaterialMassNeeded;
+ map <string , double > LambdaNeeded;
+ map <string , double > DeltaMass;
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ LambdaNeeded[(*it_s_vIV).first] = 0;
+ DeltaMass[(*it_s_vIV).first] = 0;
+ }
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ MaterialMassNeeded[(*it_s_vIV).first] = HMMass*fActualMassContentInFuel[(*it_s_vIV).first];
+ DeltaMass[(*it_s_vIV).first] = - MaterialMassNeeded[(*it_s_vIV).first];
+ }
+
+ do
+ {
+ map < string , double > MaterialMeanMolar;
+ map < string , double > MaterialMassAvailable;
+ map < string , bool > CheckOnMass;
+
+ map < string , double > LambdaPreviousStep;
+
+ double FuelMeanMolar = 0;
+
+ bool NotEnoughPuInStock = false;
+ BreakReturnLambda = false;
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ LambdaPreviousStep[(*it_s_D ).first] =LambdaNeeded[(*it_s_D ).first];
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ LambdaNeeded[(*it_s_vIV).first] = LambdaCalculation((*it_s_vIV).first, LambdaPreviousStep[(*it_s_vIV).first], MaterialMassNeeded[(*it_s_vIV).first], DeltaMass[(*it_s_vIV).first], StreamArray[(*it_s_vIV).first]);
+
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ {
+ if( LambdaNeeded[(*it_s_D).first] == -1 )
+ {
+ SetLambdaToErrorCode(lambda[(*it_s_D).first]);
+ WARNING << "Not enough : "<< (*it_s_D).first <<" to build fuel." << endl;
+ BreakReturnLambda = true;
+ if((*it_s_D).first=="PuList")
+ NotEnoughPuInStock = true;
+ }
+ }
+
+ if(BreakReturnLambda && !NotEnoughPuInStock){ return lambda;}
+
+ if(BreakReturnLambda && NotEnoughPuInStock)
+ {
+ MaterialMolarContent["PuList"] = -1;
+ break;
+ }
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ SetLambda(lambda[(*it_s_D).first], LambdaNeeded[(*it_s_D).first]);
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ IVStream[(*it_s_vIV).first].Clear();
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i<StreamArray[(*it_s_vIV).first].size(); i++)
+ IVStream[(*it_s_vIV).first] += lambda[(*it_s_vIV).first][i] * StreamArray[(*it_s_vIV).first][i];
+ }
+
+ if (loopCount > fMaxInterration)
+ {
+ ERROR << "Too much iterration in BuildFuel Method !";
+ ERROR << "Need improvement in fuel fabrication ! Ask for it or D.I.Y. !!" << endl;
+ exit(1);
+ }
+
+ /* Calcul the quantity of this composition needed to reach the burnup */
+ MaterialMolarContent = GetMolarFraction(IVStream, BurnUp);
+ //cout<<MaterialMolarContent["PuList"]<<endl;
+ for( it_s_D = MaterialMolarContent.begin(); it_s_D != MaterialMolarContent.end(); it_s_D++)
+ {
+ if( MaterialMolarContent[(*it_s_D).first] < 0 || MaterialMolarContent[(*it_s_D).first] > 1 )
+ {
+ SetLambdaToErrorCode(lambda[(*it_s_D).first]);
+ WARNING << "GetMolarFraction return negative or greater than one value, at least for one material : "<<(*it_s_D).first;
+ return lambda;
+ }
+ }
+ for( it_s_IV = IVStream.begin(); it_s_IV != IVStream.end(); it_s_IV++)
+ MaterialMassAvailable[(*it_s_IV).first] = IVStream[(*it_s_IV).first].GetTotalMass()*1e06;
+
+ for( it_s_D = MaterialMolarContent.begin(); it_s_D != MaterialMolarContent.end(); it_s_D++)
+ {
+ MaterialMeanMolar[(*it_s_D).first] = IVStream[(*it_s_D).first].GetMeanMolarMass();
+ FuelMeanMolar += MaterialMolarContent[(*it_s_D).first] * MaterialMeanMolar[(*it_s_D).first];
+ }
+
+ for( it_s_D = MaterialMolarContent.begin(); it_s_D != MaterialMolarContent.end(); it_s_D++)
+ MaterialMassContent [(*it_s_D).first] = MaterialMolarContent[(*it_s_D).first] * MaterialMeanMolar[(*it_s_D).first] / FuelMeanMolar;
+
+ fActualMolarContentInFuel = MaterialMolarContent; //fActualMolarContentInFuel is used in GetMolarFraction
+
+ for( it_s_D = MaterialMassContent.begin(); it_s_D != MaterialMassContent.end(); it_s_D++)
+ {
+ MaterialMassNeeded[(*it_s_D).first] = HMMass*MaterialMassContent[(*it_s_D).first];
+ DeltaMass[(*it_s_D).first] = MaterialMassAvailable[(*it_s_D).first] - MaterialMassNeeded[(*it_s_D).first];
+ }
+
+ for( it_s_D = MaterialMassNeeded.begin(); it_s_D != MaterialMassNeeded.end(); it_s_D++)
+ {
+ double DeltaM = fabs(MaterialMassNeeded[(*it_s_D).first] - MaterialMassAvailable[(*it_s_D).first]) / HMMass ;
+ if(DeltaM<fRelativMassPrecision) {CheckOnMass[(*it_s_D).first] = true;}
+ else{CheckOnMass[(*it_s_D).first] = false;}
+ }
+
+ FuelBuiltCorrectly = true;
+ for( it_s_B = CheckOnMass.begin(); it_s_B != CheckOnMass.end(); it_s_B++)
+ {
+ if(!CheckOnMass[(*it_s_B).first]) {FuelBuiltCorrectly = false;}
+ }
+
+ loopCount++;
+
+ }while(!FuelBuiltCorrectly);
+
+ if(MaterialMassContent["PuList"]>fMaximalPuMassContent || MaterialMolarContent["PuList"]==-1)
+ {
+ for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++)
+ {
+ for(int i=0; i<lambda[(*it_s_vD).first].size(); i++)
+ {
+ lambda[(*it_s_vD).first][i] = 0;
+ }
+ }
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ LambdaNeeded[(*it_s_D).first] = 0;
+
+ for( it_s_IV = IVStream.begin(); it_s_IV != IVStream.end(); it_s_IV++)
+ IVStream[(*it_s_IV).first] = IsotopicVector();
+
+ if(MaterialMassContent["PuList"]>fMaximalPuMassContent)
+ {
+ MaterialMolarContent["PuList"] = 0;
+ MaterialMassNeeded["PuList"] = fMaximalPuMassContent*HMMass;
+ LambdaNeeded["PuList"] = LambdaCalculation("PuList", 0., MaterialMassNeeded["PuList"], 1., StreamArray["PuList"]);
+ }
+
+ if(MaterialMolarContent["PuList"]==-1)
+ {
+ LambdaNeeded["PuList"] = LambdaCalculation("PuList", 0., fTotalMassInStocks["PuList"], 1., StreamArray["PuList"]);
+ }
+
+ SetLambda(lambda["PuList"], LambdaNeeded["PuList"]);
+
+ for( int i = 0 ; i < (int)lambda["PuList"].size() ; i++ )
+ IVStream["PuList"] +=lambda["PuList"][i] * StreamArray["PuList"][i];
+
+ fActualMassContentInFuel["PuList"] = (IVStream["PuList"].GetTotalMass()*1e6)/HMMass;
+
+ LambdaNeeded["EnrichmentList"] = LambdaCalculation("EnrichmentList", 0., HMMass, 1., StreamArray["EnrichmentList"]);
+ LambdaNeeded["FertileList"] = LambdaCalculation("FertileList", 0., HMMass, 1., StreamArray["FertileList"]);
+
+ SetLambda(lambda["EnrichmentList"], LambdaNeeded["EnrichmentList"]);
+ SetLambda(lambda["FertileList"], LambdaNeeded["FertileList"]);
+
+ for( int i = 0 ; i < (int)lambda["EnrichmentList"].size() ; i++ )
+ IVStream["EnrichmentList"] +=lambda["EnrichmentList"][i] * StreamArray["EnrichmentList"][i];
+
+ for( int i = 0 ; i < (int)lambda["FertileList"].size() ; i++ )
+ IVStream["FertileList"] +=lambda["FertileList"][i] * StreamArray["FertileList"][i];
+
+ double U5Enrichment = GetU5Enrichment(IVStream, BurnUp);
+ double MeanMolarUdepleted = U5Enrichment*IVStream["EnrichmentList"].GetMeanMolarMass() + (1 - U5Enrichment)*IVStream["FertileList"].GetMeanMolarMass();
+
+ double UdepletedMassNeeded = (HMMass - IVStream["PuList"].GetTotalMass() * 1e6);
+ double EnrichedSupportMassNeeded = U5Enrichment*UdepletedMassNeeded*IVStream["EnrichmentList"].GetMeanMolarMass()/IVStream["FertileList"].GetMeanMolarMass();
+ double FertilMassNeeded = UdepletedMassNeeded - EnrichedSupportMassNeeded;
+
+ LambdaNeeded["EnrichmentList"] = LambdaCalculation("EnrichmentList", 0., EnrichedSupportMassNeeded, 1., StreamArray["EnrichmentList"]);
+ LambdaNeeded["FertileList"] = LambdaCalculation("FertileList", 0., FertilMassNeeded, 1., StreamArray["FertileList"]);
+
+ BreakReturnLambda = false;
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ {
+ if( LambdaNeeded[(*it_s_D).first] == -1 )
+ {
+ SetLambdaToErrorCode(lambda[(*it_s_D).first]);
+ WARNING << "Not enough : "<< (*it_s_D).first <<" to build fuel." << endl;
+ BreakReturnLambda = true;
+ }
+ }
+ if(BreakReturnLambda) { return lambda;}
+
+ SetLambda(lambda["EnrichmentList"], LambdaNeeded["EnrichmentList"]);
+ SetLambda(lambda["FertileList"], LambdaNeeded["FertileList"]);
+
+ IVStream["EnrichmentList"] = IsotopicVector();
+ IVStream["FertileList"] = IsotopicVector();
+
+ for( int i = 0 ; i < (int)lambda["EnrichmentList"].size() ; i++ )
+ IVStream["EnrichmentList"] +=lambda["EnrichmentList"][i] * StreamArray["EnrichmentList"][i];
+
+ for( int i = 0 ; i < (int)lambda["FertileList"].size() ; i++ )
+ IVStream["FertileList"] +=lambda["FertileList"][i] * StreamArray["FertileList"][i];
+ }
+
+ for( it_s_D = MaterialMassNeeded.begin(); it_s_D != MaterialMassNeeded.end(); it_s_D++)
+ DBGV( "Weight percent : "<<(*it_s_D).first<<" "<< MaterialMassNeeded[(*it_s_D).first]/HMMass);
+
+ for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++)
+ {
+ DBGV( "Lambda vector : "<<(*it_s_vD).first );
+
+ for(int i=0; i<lambda[(*it_s_vD).first].size(); i++)
+ {
+ DBGV(lambda[(*it_s_vD).first][i]);
+ }
+ }
+
+ return lambda;
+}
+
+//________________________________________________________________________
+TTree* EQM_MLP_PWR_MOxEUS::CreateTMVAInputTree(IsotopicVector TheFuel, double ThisTime)
+{
+ TTree* InputTree = new TTree("EQTMP", "EQTMP");
+
+ float Pu8 = 0;
+ float Pu9 = 0;
+ float Pu10 = 0;
+ float Pu11 = 0;
+ float Pu12 = 0;
+ float Am1 = 0;
+ float U5 = 0;
+ float U8 = 0;
+ float Time = 0;
+
+ InputTree->Branch( "U5" ,&U5 ,"U5/F" );
+ InputTree->Branch( "U8" ,&U8 ,"U8/F" );
+ InputTree->Branch( "Pu8" ,&Pu8 ,"Pu8/F" );
+ InputTree->Branch( "Pu9" ,&Pu9 ,"Pu9/F" );
+ InputTree->Branch( "Pu10" ,&Pu10 ,"Pu10/F" );
+ InputTree->Branch( "Pu11" ,&Pu11 ,"Pu11/F" );
+ InputTree->Branch( "Pu12" ,&Pu12 ,"Pu12/F" );
+ InputTree->Branch( "Am1" ,&Am1 ,"Am1/F" );
+ InputTree->Branch( "Time" ,&Time ,"Time/F" );
+
+ double Ntot = TheFuel.GetSumOfAll();
+
+ U5 = TheFuel.GetZAIIsotopicQuantity(92,235,0)/Ntot;
+ U8 = TheFuel.GetZAIIsotopicQuantity(92,238,0)/Ntot;
+ Pu8 = TheFuel.GetZAIIsotopicQuantity(94,238,0)/Ntot;
+ Pu9 = TheFuel.GetZAIIsotopicQuantity(94,239,0)/Ntot;
+ Pu10 = TheFuel.GetZAIIsotopicQuantity(94,240,0)/Ntot;
+ Pu11 = TheFuel.GetZAIIsotopicQuantity(94,241,0)/Ntot;
+ Pu12 = TheFuel.GetZAIIsotopicQuantity(94,242,0)/Ntot;
+ Am1 = TheFuel.GetZAIIsotopicQuantity(95,241,0)/Ntot;
+
+ Time=ThisTime;
+
+ InputTree->Fill();
+ return InputTree;
+}
+//________________________________________________________________________
+double EQM_MLP_PWR_MOxEUS::ExecuteTMVA(TTree* theTree, string WeightPath)
+{
+ // --- Create the Reader object
+ TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+ // Create a set of variables and declare them to the reader
+ // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+ Float_t U5, U8, Pu8,Pu9,Pu10,Pu11,Pu12,Am1, Time;
+
+ reader->AddVariable( "U5" ,&U5 );
+ reader->AddVariable( "U8" ,&U8 );
+ reader->AddVariable( "Pu8" ,&Pu8 );
+ reader->AddVariable( "Pu9" ,&Pu9 );
+ reader->AddVariable( "Pu10" ,&Pu10 );
+ reader->AddVariable( "Pu11" ,&Pu11 );
+ reader->AddVariable( "Pu12" ,&Pu12 );
+ reader->AddVariable( "Am1" ,&Am1 );
+ reader->AddVariable( "Time" ,&Time );
+
+ // --- Book the MVA methods
+
+ // Book method MLP
+ TString methodName = "MLP method";
+ reader->BookMVA( methodName, fTMVAWeightPath[0] );
+ theTree->SetBranchAddress( "U5" ,&U5 );
+ theTree->SetBranchAddress( "U8" ,&U8 );
+ theTree->SetBranchAddress( "Pu8" ,&Pu8 );
+ theTree->SetBranchAddress( "Pu9" ,&Pu9 );
+ theTree->SetBranchAddress( "Pu10" ,&Pu10 );
+ theTree->SetBranchAddress( "Pu11" ,&Pu11 );
+ theTree->SetBranchAddress( "Pu12" ,&Pu12 );
+ theTree->SetBranchAddress( "Am1" ,&Am1 );
+ theTree->SetBranchAddress( "Time" ,&Time );
+
+ theTree->GetEntry(0);
+
+ Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+ delete reader;
+ return (double)val; //return k_eff
+}
+
+//________________________________________________________________________
+double EQM_MLP_PWR_MOxEUS::GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU)
+{
+ /**************************************************************************/
+ //With a dichotomy, the maximal irradiation time (TheFinalTime) is calculated
+ //When average Kinf is very close (according "Precision") to the threshold
+ //then the corresponding irradiation time is convert in burnup and returned
+ /**************************************************************************/
+ //Algorithm initialization
+
+ double TheFinalTime = BurnupToSecond(TargetBU);
+ double OldFinalTimeMinus = 0;
+ double MaximumBU = fMaximalBU;
+ double OldFinalTimePlus = BurnupToSecond(MaximumBU);
+ double k_av = 0; //average kinf
+ double OldPredictedk_av = 0;
+
+ for(int b=0;b<fNumberOfBatch;b++)
+ {
+ float TheTime = (b+1)*TheFinalTime/fNumberOfBatch;
+ TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
+ OldPredictedk_av += ExecuteTMVA(InputTree,fTMVAWeightPath[0]);
+ delete InputTree;
+ }
+ OldPredictedk_av/=fNumberOfBatch;
+
+ //Algorithm control
+ int count=0;
+ int MaximumLoopCount = 500;
+ do
+ {
+ if(count > MaximumLoopCount )
+ {
+ ERROR<<"CRITICAL ! Can't manage to predict burnup\nHint : Try to increase the precision on k effective using :\n YourEQM_MLP_Kinf->SetPCMprecision(pcm); with pcm the precision in pcm (default 10) REDUCE IT\n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr "<<endl;
+ exit(1);
+ }
+
+ if( (OldPredictedk_av-fKThreshold) > 0) //The burnup can be increased
+ {
+ OldFinalTimeMinus = TheFinalTime;
+ TheFinalTime = TheFinalTime + fabs(OldFinalTimePlus - TheFinalTime)/2.;
+
+ if(SecondToBurnup(TheFinalTime) >= (MaximumBU-MaximumBU*GetBurnUpPrecision() ) )
+ return MaximumBU;
+ }
+
+ else if( (OldPredictedk_av-fKThreshold) < 0)//The burnup is too high
+ {
+ OldFinalTimePlus = TheFinalTime;
+ TheFinalTime = TheFinalTime - fabs(OldFinalTimeMinus - TheFinalTime)/2.;
+ if( SecondToBurnup(TheFinalTime) < TargetBU*GetBurnUpPrecision() )
+ return 0;
+ }
+
+ k_av = 0;
+ for(int b=0;b<fNumberOfBatch;b++)
+ {
+ float TheTime = (b+1)*TheFinalTime/fNumberOfBatch;
+ TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
+ k_av += ExecuteTMVA(InputTree,fTMVAWeightPath[0]);
+ delete InputTree;
+ }
+ k_av/=fNumberOfBatch;
+
+ OldPredictedk_av = k_av;
+ count++;
+ } while( fabs(OldPredictedk_av-fKThreshold) > GetPCMPrecision() ) ;
+
+ return SecondToBurnup(TheFinalTime);
+}
+
+//________________________________________________________________________
+map < string , double> EQM_MLP_PWR_MOxEUS::GetMolarFraction(map <string , IsotopicVector > IVStream, double TargetBU)
+{
+
+ IsotopicVector Pu = IVStream["PuList"];
+ IsotopicVector EnrichedSupport = IVStream["EnrichmentList"];
+ IsotopicVector Fertil = IVStream["FertileList"];
+
+ map < string , double> MolarFraction ;
+
+ double OldPuContentMinus = 0;
+ double OldPuContentPlus = fMaximalPuMolarContent;
+
+ double PredictedBU = 0 ;
+ double PuContent = fActualMolarContentInFuel["PuList"];
+ double U5Enrichment = fMinimalU5Enrichment;
+
+ IsotopicVector FreshFuel = PuContent*(Pu/Pu.GetSumOfAll()) + U5Enrichment*(1-PuContent) *(EnrichedSupport/EnrichedSupport.GetSumOfAll()) + (1 - U5Enrichment)*(1-PuContent)*(Fertil/Fertil.GetSumOfAll());
+ double OldPredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+
+ double Precision = GetBurnUpPrecision()*TargetBU; //1 % of the targeted burnup
+ int count = 0;
+ int MaximumLoopCount = 500;
+ do
+ {
+ if(count > MaximumLoopCount )
+ {
+ ERROR<<"CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on burnup using :\nYourEQM_MLP_Kinf->SetBurnUpPrecision(prop); with prop the precision (default 0.5percent : 0.005) INCREASE IT\nIf this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr "<<endl;
+ exit(1);
+ }
+ if( (OldPredictedBU - TargetBU) < 0 ) //The Content can be increased
+ {
+ OldPuContentMinus = PuContent;
+ PuContent = PuContent + fabs(OldPuContentPlus-PuContent)/2.;
+ }
+ else if( (OldPredictedBU - TargetBU) > 0) //The Content is too high
+ {
+ OldPuContentPlus = PuContent;
+ PuContent = PuContent - fabs(OldPuContentMinus-PuContent)/2.;
+ }
+
+ IsotopicVector FreshFuel= PuContent*(Pu/Pu.GetSumOfAll()) + U5Enrichment*(1-PuContent) *(EnrichedSupport/EnrichedSupport.GetSumOfAll()) + (1 - U5Enrichment)*(1-PuContent)*(Fertil/Fertil.GetSumOfAll());
+ PredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+
+ OldPredictedBU = PredictedBU;
+ count ++;
+INFO<<"count "<<count<<" => Target BU "<<TargetBU<<" - Predicted BU "<<PredictedBU<<endl;
+
+ }while(fabs(TargetBU-PredictedBU)>Precision);
+
+ DBGV("Predicted BU "<<PredictedBU<<" PuContent "<<PuContent);
+
+ MolarFraction["PuList"] = PuContent;
+ MolarFraction["EnrichmentList"] = (1 - MolarFraction["PuList"])*fMinimalU5Enrichment;
+ MolarFraction["FertileList"] = 1 - MolarFraction["PuList"] - MolarFraction["EnrichmentList"];
+
+return MolarFraction;
+}
+//________________________________________________________________________
+
+//________________________________________________________________________
+double EQM_MLP_PWR_MOxEUS::GetU5Enrichment(map <string , IsotopicVector > IVStream, double TargetBU)
+{
+ IsotopicVector Pu = IVStream["PuList"];
+ IsotopicVector EnrichedSupport = IVStream["EnrichmentList"];
+ IsotopicVector Fertil = IVStream["FertileList"];
+
+ double OldU5EnrichmentMinus = fMinimalU5Enrichment;
+ double OldU5EnrichmentPlus = fMaximalU5Enrichment;
+ double U5Enrichment = fMinimalU5Enrichment;
+
+ double PredictedBU = 0 ;
+ double USupportMassNeeded = (Pu.GetTotalMass()*1e6)/fActualMassContentInFuel["PuList"] - (Pu.GetTotalMass()*1e6);
+
+ IsotopicVector USupport = U5Enrichment*EnrichedSupport + (1 - U5Enrichment)*Fertil;
+ double USupportMassContent = USupportMassNeeded/(USupport.GetTotalMass()*1e6);
+
+ IsotopicVector FreshFuel = Pu + USupportMassContent*USupport;
+ FreshFuel = FreshFuel/FreshFuel.GetSumOfAll();
+
+ double OldPredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+
+ double Precision = GetBurnUpPrecision()*TargetBU; //1 % of the targeted burnup
+ int count = 0;
+ int MaximumLoopCount = 500;
+
+ do
+ {
+
+ if(count > MaximumLoopCount )
+ {
+ ERROR<<"CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on burnup using :\nYourEQM_MLP_Kinf->SetBurnUpPrecision(prop); with prop the precision (default 0.5percent : 0.005) INCREASE IT\nIf this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr "<<endl;
+ exit(1);
+ }
+ if( (OldPredictedBU - TargetBU) < 0 ) //The Content can be increased
+ {
+ OldU5EnrichmentMinus = U5Enrichment;
+ U5Enrichment = U5Enrichment + fabs(OldU5EnrichmentPlus-U5Enrichment)/2.;
+ }
+ else if( (OldPredictedBU - TargetBU) > 0) //The Content is too high
+ {
+ OldU5EnrichmentPlus = U5Enrichment;
+ U5Enrichment = U5Enrichment - fabs(OldU5EnrichmentMinus-U5Enrichment)/2.;
+ }
+
+ IsotopicVector USupport = U5Enrichment*EnrichedSupport + (1 - U5Enrichment)*Fertil;
+
+ double USupportMassContent = USupportMassNeeded/(USupport.GetTotalMass()*1e6);
+ IsotopicVector FreshFuel = Pu + USupportMassContent*USupport;
+ FreshFuel = FreshFuel/FreshFuel.GetSumOfAll();
+
+ PredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+
+ OldPredictedBU = PredictedBU;
+ count ++;
+ }while(fabs(TargetBU-PredictedBU)>Precision);
+
+DBGV("Predicted BU "<<PredictedBU<<" U5Enrichment "<<U5Enrichment);
+return U5Enrichment;
+
+}
diff --git a/source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.hxx b/source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.hxx
new file mode 100755
index 000000000..e4a138cf3
--- /dev/null
+++ b/source/branches/SL3/Model/Equivalence/EQM_MLP_PWR_MOxEUS.hxx
@@ -0,0 +1,67 @@
+#ifndef _EQM_MLP_PWR_MOxEUS_HXX
+#define _EQM_MLP_PWR_MOxEUS_HXX
+
+#include "EquivalenceModel.hxx"
+#include "TTree.h"
+
+/*!
+ \file
+ \brief Header file for EQM_MLP_PWR_MOxEUS
+
+
+ @author FC
+ @version 1.0
+ */
+class EQM_MLP_PWR_MOxEUS : public EquivalenceModel
+{
+ public :
+
+ EQM_MLP_PWR_MOxEUS(string TMVAWeightPath, int NumOfBatch, double CriticalityThreshold, double MaximalPuMassContent);
+ EQM_MLP_PWR_MOxEUS(CLASSLogger* log, string TMVAWeightPath, int NumOfBatch, double CriticalityThreshold, double MaximalPuMassContent);
+
+ virtual map <string , vector<double> > BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray);
+
+ void SetSpecificPower(double SpecificPower) {fSpecificPower = SpecificPower;}
+ void SetMinimalU5Enrichment(double MinEU5) {fMinimalU5Enrichment= MinEU5;}
+ void SetMaximalU5Enrichment(double MaxEU5){fMaximalU5Enrichment= MaxEU5;}
+ void SetBurnUpPrecision(double precision) {fBurnUpPrecision= precision;}
+ void SetPCMPrecision(double prop) {fPCMPrecision= prop;}
+
+ double GetBurnUpPrecision(){return fBurnUpPrecision;}
+ double GetPCMPrecision(){return fPCMPrecision/1e5;}
+ double GetMinimalU5Enrichment(){return fMinimalU5Enrichment;};
+
+ TTree* CreateTMVAInputTree(IsotopicVector TheFuel, double ThisTime);
+ double ExecuteTMVA(TTree* theTree, string WeightPath);
+ double GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU);
+ double GetU5Enrichment(map <string , IsotopicVector > IVStream, double TargetBU);
+
+ map < string , double> GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp);
+
+ private:
+
+ double BurnupToSecond(double BurnUp){return BurnUp/fSpecificPower*(24*3.6e6);}
+ double SecondToBurnup(double Second){return Second*fSpecificPower/(24*3.6e6);}
+
+ int fNumberOfBatch ;
+ double fKThreshold ;
+ double fSpecificPower ;
+ double fMaximalBU;
+ double fBurnUpPrecision;
+ double fPCMPrecision;
+
+ double fActualPuMolarContent;
+ double fActualPuMassContent;
+ double fMaximalPuMassContent ;
+ double fMaximalPuMolarContent;
+ double fMinimalU5Enrichment;
+ double fMaximalU5Enrichment;
+
+ vector <string> fTMVAWeightPath; //!<The weight needed by TMVA to construct and execute the multilayer perceptron
+
+
+
+
+};
+
+#endif
\ No newline at end of file
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.cxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.cxx
new file mode 100755
index 000000000..eea784b6c
--- /dev/null
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.cxx
@@ -0,0 +1,554 @@
+#include "EQM_PWR_FixedContent.hxx"
+#include "CLASSLogger.hxx"
+#include "CLASSMethod.hxx"
+#include "StringLine.hxx"
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+#include <cassert>
+
+#include "TSystem.h"
+#include "TMVA/Reader.h"
+#include "TMVA/Tools.h"
+#include "TMVA/MethodCuts.h"
+
+
+
+//________________________________________________________________________
+EQM_PWR_FixedContent::EQM_PWR_FixedContent(string TMVAWeightPath, int NumOfBatch, string InformationFile, double CriticalityThreshold):EquivalenceModel(new CLASSLogger("EQM_PWR_FixedContent.log"))
+{
+ /**The information file and tmva weight*/
+ fTMVAWeightPath.push_back(TMVAWeightPath);
+
+ /* INFORMATION FILE HANDLING */
+
+ if(InformationFile == "")
+ InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+
+ fMaximalBU = 0;
+ fMaximalContent = 0;
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold ;
+
+ SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+ SetPCMPrecision(10);
+
+ fInformationFile = InformationFile;
+
+ LoadKeyword();
+ ReadNFO();//Getting information from fInformationFile
+
+ INFO << "__An equivalence model has been define__" << endl;
+ INFO << "\tThis model is based on the prediction of kinf" << endl;
+ INFO << "\tThe TMVA (weight | information) files are :" << endl;
+ INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fInformationFile << " )" << endl;
+ INFO << "Maximal fissile content (molar proportion) : " << fMaximalContent << endl;
+ INFO << "Maximal burnup (GWd/tHM) : " << fMaximalBU << endl;
+
+ EquivalenceModel::PrintInfo();
+
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ bool EmptyList = false;
+ bool FirstContentNULL = false;
+ bool FixedContentNULL = false;
+
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ if(fStreamList[(* it_s_IV).first].GetIsotopicQuantity().empty()){EmptyList=true;}
+ }
+
+ for( it_s_D = fFirstGuessContent .begin(); it_s_D != fFirstGuessContent.end(); it_s_D++)
+ {
+ if(fFirstGuessContent[(* it_s_D).first] ==0){FirstContentNULL=true;}
+ }
+
+ for( it_s_D = fFixedMassContent .begin(); it_s_D != fFixedMassContent.end(); it_s_D++)
+ {
+ if(fFixedMassContent[(* it_s_D).first] ==0){FixedContentNULL=true;}
+ }
+
+ if(fMapOfTMVAVariableNames.empty() || EmptyList==true || FirstContentNULL==true || FixedContentNULL==true || fMaximalBU == 0 || fMaximalContent == 0 )
+ {
+ ERROR<<"Missing information file in : "<<fInformationFile<<endl;
+ exit(1);
+ }
+
+}
+//________________________________________________________________________
+EQM_PWR_FixedContent::EQM_PWR_FixedContent(CLASSLogger* log, string TMVAWeightPath, int NumOfBatch, string InformationFile, double CriticalityThreshold):EquivalenceModel(log)
+{
+
+ /**The information file and tmva weight*/
+ fTMVAWeightPath.push_back(TMVAWeightPath);
+
+ if(InformationFile == "")
+ InformationFile = StringLine::ReplaceAll(TMVAWeightPath,".xml",".nfo");
+
+ fMaximalBU = 0;
+ fMaximalContent = 0;
+ fNumberOfBatch = NumOfBatch;
+ fKThreshold = CriticalityThreshold;
+
+ SetBurnUpPrecision(0.005);//1 % of the targeted burnup
+ SetPCMPrecision(10);
+
+ fInformationFile = InformationFile;
+
+ LoadKeyword();
+ ReadNFO();//Getting information from fMLPInformationFile
+
+ INFO << "__An equivalence model has been define__" << endl;
+ INFO << "\tThis model is based on the prediction of kinf" << endl;
+ INFO << "\tThe TMVA (weight | information) files are :" << endl;
+ INFO << "\t" << "( " << fTMVAWeightPath[0] << " | " << fInformationFile << " )" << endl;
+ INFO << "Maximal fissile content (molar proportion) : " << fMaximalContent << endl;
+ INFO << "Maximal burnup (GWd/tHM) : " << fMaximalBU << endl;
+
+ EquivalenceModel::PrintInfo();
+
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+
+ bool EmptyList = false;
+ bool FirstContentNULL = false;
+ bool FixedContentNULL = false;
+
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ if(fStreamList[(* it_s_IV).first].GetIsotopicQuantity().empty()){EmptyList=true;}
+ }
+
+ for( it_s_D = fFirstGuessContent .begin(); it_s_D != fFirstGuessContent.end(); it_s_D++)
+ {
+ if(fFirstGuessContent[(* it_s_D).first] ==0){FirstContentNULL=true;}
+ }
+
+ for( it_s_D = fFixedMassContent .begin(); it_s_D != fFixedMassContent.end(); it_s_D++)
+ {
+ if(fFixedMassContent[(* it_s_D).first] ==0){FixedContentNULL=true;}
+ }
+
+ if(fMapOfTMVAVariableNames.empty() || EmptyList==true || FirstContentNULL==true || FixedContentNULL==true || fMaximalBU == 0 || fMaximalContent == 0 )
+ {
+ ERROR<<"Missing information file in : "<<fInformationFile<<endl;
+ exit(1);
+ }
+}
+//________________________________________________________________________
+void EQM_PWR_FixedContent::LoadKeyword()
+{
+ DBGL
+
+ fDKeyword.insert( pair<string, PWR_Fixed_DMthPtr>( "k_zainame" , &EQM_PWR_FixedContent::ReadZAIName) );
+ fDKeyword.insert( pair<string, PWR_Fixed_DMthPtr>( "k_maxburnup" , &EQM_PWR_FixedContent::ReadMaxBurnUp) );
+ fDKeyword.insert( pair<string, PWR_Fixed_DMthPtr>( "k_maxfiscontent" , &EQM_PWR_FixedContent::ReadMaxFisContent) );
+ fDKeyword.insert( pair<string, PWR_Fixed_DMthPtr>( "k_fixedmasscontent" , &EQM_PWR_FixedContent::ReadFixedMassContent) );
+ DBGL
+}
+//________________________________________________________________________
+void EQM_PWR_FixedContent::ReadZAIName(const string &line)
+{
+ DBGL
+
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+ if( keyword != "k_zainame" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"k_zainame\" not found !" << endl;
+ exit(1);
+ }
+
+ int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+
+ string name = StringLine::NextWord(line, pos, ' ');
+
+ fMapOfTMVAVariableNames.insert( pair<ZAI,string>( ZAI(Z, A, I), name ) );
+
+ DBGL
+}
+//________________________________________________________________________
+void EQM_PWR_FixedContent::ReadMaxBurnUp(const string &line)
+{
+ DBGL
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+ if( keyword != "k_maxburnup" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"k_maxburnup\" not found !" << endl;
+ exit(1);
+ }
+
+ fMaximalBU = atof(StringLine::NextWord(line, pos, ' ').c_str());
+
+ DBGL
+}
+//________________________________________________________________________
+void EQM_PWR_FixedContent::ReadMaxFisContent(const string &line)
+{
+ DBGL
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+ if( keyword != "k_maxfiscontent" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"k_maxfiscontent\" not found !" << endl;
+ exit(1);
+ }
+
+ fMaximalContent = atof(StringLine::NextWord(line, pos, ' ').c_str());
+
+ DBGL
+}
+//________________________________________________________________________
+void EQM_PWR_FixedContent::ReadLine(string line)
+{
+ DBGL
+
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+
+ map<string, PWR_Fixed_DMthPtr>::iterator it = fDKeyword.find(keyword);
+
+ if(it != fDKeyword.end())
+ (this->*(it->second))( line );
+
+ DBGL
+}
+//________________________________________________________________________
+//________________________________________________________________________
+void EQM_PWR_FixedContent::ReadFixedMassContent(const string &line)
+{
+ DBGL
+ int pos = 0;
+ string keyword = tlc(StringLine::NextWord(line, pos, ' '));
+ if( keyword != "k_fixedmasscontent" ) // Check the keyword
+ {
+ ERROR << " Bad keyword : \"k_fixedmasscontent\" not found !" << endl;
+ exit(1);
+ }
+ string ListName = StringLine::NextWord(line, pos, ' ');
+ double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
+ fFixedMassContent[ListName] = Q;
+
+ DBGL
+}
+
+//________________________________________________________________________
+map <string , vector<double> > EQM_PWR_FixedContent::BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray)
+{
+ map <string , vector<double> > lambda ; // map containing name of the list and associated vector of proportions taken from stocks
+
+ //Iterators declaration
+ map < string , vector <IsotopicVector> >::iterator it_s_vIV;
+ map < string , vector <double> >::iterator it_s_vD;
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ map < string , bool >::iterator it_s_B;
+
+ map <string , IsotopicVector > IVStream;
+ map <string , double > MaterialMolarContent;
+ map <string , double > MaterialMassContent;
+ map <string , double > MaterialMassNeeded;
+ map <string , double > LambdaNeeded;
+ map <string , double > DeltaMass;
+
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i<StreamArray[(*it_s_vIV).first].size(); i++)
+ {
+ lambda[(*it_s_vIV).first].push_back(0);
+ }
+ }
+
+ /*** Test if there is stocks **/
+ bool BreakReturnLambda = false;
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ if(StreamArray[(*it_s_vIV).first].size() == 0)
+ {
+ WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<". Fuel not build." << endl;
+ SetLambdaToErrorCode(lambda[(*it_s_vIV).first]);
+ BreakReturnLambda = true;
+ }
+ }
+ if(BreakReturnLambda) { return lambda;}
+ HMMass *= 1e6; //Unit onversion : tons to gram
+
+ /**** Some initializations **/
+ StocksTotalMassCalculation(StreamArray);
+
+ MaterialMassContent = GetFixedMassContent();
+ bool FuelBuiltCorrectly = false ;
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ LambdaNeeded[(*it_s_vIV).first] = 0;
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ MaterialMassNeeded[(*it_s_vIV).first] = HMMass*MaterialMassContent[(*it_s_vIV).first];
+ DeltaMass[(*it_s_vIV).first] = fTotalMassInStocks[(* it_s_vIV).first] - MaterialMassNeeded[(*it_s_vIV).first];
+ }
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ LambdaNeeded[(*it_s_vIV).first] = LambdaCalculation((*it_s_vIV).first, 0., MaterialMassNeeded[(*it_s_vIV).first], DeltaMass[(*it_s_vIV).first], StreamArray[(*it_s_vIV).first]);
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ SetLambda(lambda[(*it_s_D).first], LambdaNeeded[(*it_s_D).first]);
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i<StreamArray[(*it_s_vIV).first].size(); i++)
+ IVStream[(*it_s_vIV).first] += lambda[(*it_s_vIV).first][i] * StreamArray[(*it_s_vIV).first][i];
+ }
+
+ IsotopicVector FreshFuel;
+
+ for( it_s_IV = IVStream.begin(); it_s_IV != IVStream.end(); it_s_IV++)
+ FreshFuel += MaterialMassContent[(*it_s_IV).first]*(IVStream[(*it_s_IV).first]/IVStream[(*it_s_IV).first].GetSumOfAll());
+
+ double PredictedBU = GetMaximumBurnUp(FreshFuel,BurnUp);
+
+ if(PredictedBU<BurnUp)
+ {
+ WARNING << " Fixed contents are not adequate. Predicted BU is lower than target BU. "<< endl;
+ WARNING << " Predicted BU : "<<PredictedBU<<" Target BU: "<<BurnUp<< endl;
+ }
+
+ if(PredictedBU>BurnUp)
+ {
+ WARNING << " Fixed contents are not adequate. Predicted BU is greater than target BU. "<< endl;
+ WARNING << " Predicted BU : "<<PredictedBU<<" Target BU: "<<BurnUp<< endl;
+ }
+
+ return lambda;
+}
+
+TTree* EQM_PWR_FixedContent::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime)
+{
+ /******Create Input data tree to be interpreted by TMVA::Reader***/
+ TTree* InputTree = new TTree("InTMPKinf", "InTMPKinf");
+
+ vector<float> InputTMVA;
+ for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+ InputTMVA.push_back(0);
+
+ float Time = 0;
+
+ IsotopicVector IVInputTMVA;
+ map<ZAI ,string >::iterator it;
+ int j = 0;
+
+ for( it = fMapOfTMVAVariableNames.begin() ; it != fMapOfTMVAVariableNames.end() ; it++)
+ {
+ InputTree->Branch( ((*it).second).c_str() ,&InputTMVA[j], ((*it).second + "/F").c_str());
+ IVInputTMVA+= ((*it).first)*1;
+ j++;
+ }
+
+ if(ThisTime != -1)
+ InputTree->Branch( "Time" ,&Time ,"Time/F" );
+
+ IsotopicVector IVAccordingToUserInfoFile = TheFreshfuel.GetThisComposition(IVInputTMVA);
+
+ double Ntot = IVAccordingToUserInfoFile.GetSumOfAll();
+
+ IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot;
+
+ j = 0;
+ map<ZAI ,string >::iterator it2;
+
+ for( it2 = fMapOfTMVAVariableNames.begin() ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+ {
+ InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it2).first ) ;
+ j++;
+ }
+
+ Time = ThisTime;
+
+ InputTree->Fill();
+
+ return InputTree;
+
+}
+//________________________________________________________________________
+double EQM_PWR_FixedContent::ExecuteTMVA(TTree* InputTree,string WeightPath)
+{
+
+ // --- Create the Reader object
+ TMVA::Reader *reader = new TMVA::Reader( "Silent" );
+
+ // Create a set of variables and declare them to the reader
+ // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
+ vector<float> InputTMVA;
+ for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++)
+ InputTMVA.push_back(0);
+ Float_t Time;
+
+ map<ZAI ,string >::iterator it;
+ int j = 0;
+ for( it = fMapOfTMVAVariableNames.begin() ; it != fMapOfTMVAVariableNames.end() ; it++)
+ { reader->AddVariable( ( (*it).second ).c_str(),&InputTMVA[j]);
+ j++;
+ }
+
+ reader->AddVariable( "Time" ,&Time);
+
+ // --- Book the MVA methods
+
+ // Book method MLP
+ TString methodName = "MLP method";
+ reader->BookMVA( methodName, WeightPath );
+
+ map<ZAI ,string >::iterator it2;
+ j = 0;
+ for( it2 = fMapOfTMVAVariableNames.begin() ; it2 != fMapOfTMVAVariableNames.end() ; it2++)
+ {
+ InputTree->SetBranchAddress(( (*it2).second ).c_str(),&InputTMVA[j]);
+ j++;
+ }
+
+ InputTree->SetBranchAddress( "Time" ,&Time );
+
+ InputTree->GetEntry(0);
+ Float_t val = (reader->EvaluateRegression( methodName ))[0];
+
+ delete reader;
+
+ return (double)val;//retourn k_{inf}(t = Time)
+}
+//________________________________________________________________________
+double EQM_PWR_FixedContent::GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU)
+{
+ /**************************************************************************/
+ //With a dichotomy, the maximal irradiation time (TheFinalTime) is calculated
+ //When average Kinf is very close (according "Precision") to the threshold
+ //then the corresponding irradiation time is convert in burnup and returned
+ /**************************************************************************/
+ //Algorithm initialization
+ double TheFinalTime = BurnupToSecond(TargetBU);
+ double OldFinalTimeMinus = 0;
+ double MaximumBU = fMaximalBU;
+ double OldFinalTimePlus = BurnupToSecond(MaximumBU);
+ double k_av = 0; //average kinf
+ double OldPredictedk_av = 0;
+ for(int b = 0;b<fNumberOfBatch;b++)
+ {
+ float TheTime = (b+1)*TheFinalTime/fNumberOfBatch;
+ TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
+ OldPredictedk_av += ExecuteTMVA(InputTree,fTMVAWeightPath[0]);
+ delete InputTree;
+ }
+ OldPredictedk_av/= fNumberOfBatch;
+
+ //Algorithm control
+ int count = 0;
+ int MaximumLoopCount = 500;
+ do
+ {
+ if(count > MaximumLoopCount )
+ {
+ ERROR << "CRITICAL ! Can't manage to predict burnup\nHint : Try to increase the precision on k effective using :\n YourEQM_PWR_FixedContent->SetPCMPrecision(pcm); with pcm the precision in pcm (default 10) REDUCE IT\n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
+ exit(1);
+ }
+
+ if( (OldPredictedk_av-fKThreshold) > 0) //The burnup can be increased
+ {
+ OldFinalTimeMinus = TheFinalTime;
+ TheFinalTime = TheFinalTime + fabs(OldFinalTimePlus - TheFinalTime)/2.;
+
+ if(SecondToBurnup(TheFinalTime) >= (MaximumBU-MaximumBU*GetBurnUpPrecision() ) )
+ return MaximumBU;
+ }
+
+ else if( (OldPredictedk_av-fKThreshold) < 0)//The burnup is too high
+ {
+ OldFinalTimePlus = TheFinalTime;
+ TheFinalTime = TheFinalTime - fabs(OldFinalTimeMinus-TheFinalTime)/2.;
+ if( SecondToBurnup(TheFinalTime) < TargetBU*GetBurnUpPrecision() )
+ return 0;
+ }
+
+ k_av = 0;
+ for(int b = 0;b<fNumberOfBatch;b++)
+ {
+ float TheTime = (b+1)*TheFinalTime/fNumberOfBatch;
+ TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime);
+ k_av += ExecuteTMVA(InputTree,fTMVAWeightPath[0]);
+ delete InputTree;
+ }
+ k_av/= fNumberOfBatch;
+
+ OldPredictedk_av = k_av;
+ count++;
+
+ } while( fabs(OldPredictedk_av-fKThreshold) > GetPCMPrecision() ) ;
+
+ return SecondToBurnup(TheFinalTime);
+}
+//________________________________________________________________________
+map < string , double> EQM_PWR_FixedContent::GetMolarFraction(map <string , IsotopicVector > IVStream, double TargetBU)
+{
+ //initialization
+ IsotopicVector Fissil = IVStream["Fissile"];
+ IsotopicVector Fertil = IVStream["Fertile"];
+
+ map < string , double> MolarFraction ;
+
+ double FissileContent = fActualMolarContentInFuel["Fissile"];
+ double OldFissileContentMinus = 0;
+ double OldFissileContentPlus = fMaximalContent;
+
+ double PredictedBU = 0 ;
+ IsotopicVector FreshFuel = (1-FissileContent)*(Fertil/Fertil.GetSumOfAll()) + FissileContent*(Fissil/Fissil.GetSumOfAll());
+ double OldPredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+
+ double Precision = GetBurnUpPrecision()*TargetBU; //1 % of the targeted burnup
+ int count = 0;
+ int MaximumLoopCount = 500;
+ do
+ {
+ if(count > MaximumLoopCount )
+ {
+ ERROR << "CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on burnup using :\nYourEQM_PWR_FixedContent->SetBurnUpPrecision(prop); with prop the precision (default 0.5percent : 0.005) INCREASE IT\nIf this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl;
+ ERROR << "Targeted Burnup :" <<TargetBU<<endl;
+ ERROR << "Last calculated Burnup :" <<OldPredictedBU<<endl;
+ ERROR << "Last Fresh fuel composition :" <<endl;
+ ERROR << FreshFuel.sPrint()<<endl;
+
+ exit(1);
+ }
+
+ if( (OldPredictedBU - TargetBU) < 0 ) //The Content can be increased
+ {
+ OldFissileContentMinus = FissileContent;
+ FissileContent = FissileContent + fabs(OldFissileContentPlus-FissileContent)/2.;
+ }
+ else if( (OldPredictedBU - TargetBU) > 0) //The Content is too high
+ {
+ OldFissileContentPlus = FissileContent;
+ FissileContent = FissileContent - fabs(OldFissileContentMinus-FissileContent)/2.;
+ }
+
+ IsotopicVector FreshFuel = (1-FissileContent)*(Fertil/Fertil.GetSumOfAll()) + FissileContent*(Fissil/Fissil.GetSumOfAll());
+ PredictedBU = GetMaximumBurnUp(FreshFuel,TargetBU);
+
+ OldPredictedBU = PredictedBU;
+ count ++;
+
+ }while(fabs(TargetBU-PredictedBU)>Precision);
+
+ DBGV("Predicted BU " << PredictedBU << " FissileContent " << FissileContent);
+
+ MolarFraction["Fissile"] = FissileContent;
+ MolarFraction["Fertile"] = 1 - FissileContent;
+
+ return MolarFraction;
+}
+
+
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.hxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.hxx
new file mode 100755
index 000000000..1185f234e
--- /dev/null
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_FixedContent.hxx
@@ -0,0 +1,75 @@
+#ifndef _EQM_PWR_FixedContent_HXX
+#define _EQM_PWR_FixedContent_HXX
+
+#include "EquivalenceModel.hxx"
+#include "TTree.h"
+#include <map>
+
+using namespace std;
+
+class EQM_PWR_FixedContent;
+#ifndef __CINT__
+typedef void (EQM_PWR_FixedContent::*PWR_Fixed_DMthPtr)( const string & ) ;
+#endif
+
+/*!
+ \file
+ \brief Header file for EQM_PWR_FixedContent
+
+
+ @author FC
+ @version 1.0
+ */
+class EQM_PWR_FixedContent : public EquivalenceModel
+{
+ public :
+
+ EQM_PWR_FixedContent(string TMVAWeightPath, int NumOfBatch, string InformationFile = "", double CriticalityThreshold=1.01);
+ EQM_PWR_FixedContent(CLASSLogger* log, string TMVAWeightPath, int NumOfBatch, string InformationFile = "", double CriticalityThreshold=1.01);
+
+ virtual map <string , vector<double> > BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray);
+
+ void SetSpecificPower(double SpecificPower) {fSpecificPower = SpecificPower;}
+ void SetBurnUpPrecision(double precision) {fBurnUpPrecision= precision;}
+ void SetPCMPrecision(double prop) {fPCMPrecision= prop;}
+
+ map < string, double> GetFixedMassContent(){return fFixedMassContent;}
+ double GetBurnUpPrecision(){return fBurnUpPrecision;}
+ double GetPCMPrecision(){return fPCMPrecision/1e5;}
+
+ TTree* CreateTMVAInputTree(IsotopicVector TheFuel, double ThisTime);
+ double ExecuteTMVA(TTree* theTree, string WeightPath);
+ double GetMaximumBurnUp(IsotopicVector TheFuel, double TargetBU);
+
+ void GetModelInformation();//!<Read the fMLPInformationFile and fill containers and variables
+ void LoadKeyword();
+ void ReadZAIName(const string &line);
+ void ReadMaxBurnUp(const string &line);
+ void ReadMaxFisContent(const string &line);
+ void ReadFixedMassContent(const string &line);
+ void ReadLine(string line);
+
+ map < string , double> GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp);
+
+ double BurnupToSecond(double BurnUp){return BurnUp/fSpecificPower*(24*3.6e6);}
+ double SecondToBurnup(double Second){return Second*fSpecificPower/(24*3.6e6);}
+
+ private :
+ vector <string> fTMVAWeightPath; //!<The weight needed by TMVA to construct and execute the multilayer perceptron
+
+#ifndef __CINT__
+ map<string, PWR_Fixed_DMthPtr> fDKeyword;
+#endif
+
+ int fNumberOfBatch ;
+ double fKThreshold ;
+ double fMaximalBU;
+ double fBurnUpPrecision;
+ double fPCMPrecision;
+ map < string, double> fFixedMassContent;
+
+ map<ZAI,string> fMapOfTMVAVariableNames;//!< List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step
+
+};
+
+#endif
\ No newline at end of file
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.cxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.cxx
index 6a1a76d03..a9a9033df 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.cxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.cxx
@@ -57,8 +57,6 @@ EQM_PWR_LIN_MOX::EQM_PWR_LIN_MOX(string WeightPath):EquivalenceModel(new CLASSLo
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
-
- EquivalenceModel::PrintInfo();
}
@@ -121,7 +119,7 @@ EQM_PWR_LIN_MOX::~EQM_PWR_LIN_MOX()
}
//________________________________________________________________________
-vector<double> EQM_PWR_LIN_MOX::BuildFuel(double BurnUp, double HMMass,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
+vector<double> EQM_PWR_LIN_MOX::BuildFuel(double BurnUp, double HMMass,vector <IsotopicVector> IVStream)
{
//-----------------------------------------------------------------------------//
@@ -135,6 +133,8 @@ vector<double> EQM_PWR_LIN_MOX::BuildFuel(double BurnUp, double HMMass,vector<Is
//-----------------------------------------------------------------------------//
//-----------------------------------------------------------------------------//
+ vector<IsotopicVector> FissilArray = StreamArray[0];
+ vector<IsotopicVector> FertilArray = StreamArray[1];
vector<double> lambda;
for(int i = 0; i < (int) (FissilArray.size() + FertilArray.size()); i++)
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.hxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.hxx
index 4946ad7d7..7f31adfe0 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.hxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_LIN_MOX.hxx
@@ -58,7 +58,7 @@ class EQM_PWR_LIN_MOX : public EquivalenceModel
~EQM_PWR_LIN_MOX();
//@}
- virtual vector<double> BuildFuel(double BurnUp, double HMMass, vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray );
+ virtual vector<double> BuildFuel(double BurnUp, double HMMass, vector < vector <IsotopicVector> > StreamArray);
private :
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.cxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.cxx
index af3b7c820..98ea3fa66 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.cxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.cxx
@@ -39,16 +39,16 @@ EQM_PWR_MLP_MOX::EQM_PWR_MLP_MOX(string TMVAWeightPath):EquivalenceModel(new CLA
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
- fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
- fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
-
- SetBuildFuelFirstGuess(0.04);
+ fStreamList["Fissile"] = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+ fStreamList["Fertile"] = U5*U5_enrich + U8*(1-U5_enrich);
+
+ fFirstGuessContent["Fissile"] = 0.02;
+ fFirstGuessContent["Fertile"] = 1 - fFirstGuessContent["Fissile"];
INFO << "__An equivalence model of PWR MOX has been define__" << endl;
INFO << "\tThis model is based on a multi layer perceptron" << endl;
INFO << "\t\tThe TMVA weight file is :" << endl;
INFO << "\t\t\t" << fTMVAWeightPath << endl;
- EquivalenceModel::PrintInfo();
}
@@ -67,21 +67,21 @@ EQM_PWR_MLP_MOX::EQM_PWR_MLP_MOX(CLASSLogger* log, string TMVAWeightPath):Equiva
ZAI Pu1(94,241,0);
ZAI Pu2(94,242,0);
- fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
- fFertileList = U5*U5_enrich + U8*(1-U5_enrich);
-
- SetBuildFuelFirstGuess(0.04);
+ fStreamList["Fissile"] = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
+ fStreamList["Fertile"] = U5*U5_enrich + U8*(1-U5_enrich);
+
+ fFirstGuessContent["Fissile"] = 0.02;
+ fFirstGuessContent["Fertile"] = 1 - fFirstGuessContent["Fissile"];
INFO << "__An equivalence model of PWR MOX has been define__" << endl;
INFO << "\tThis model is based on a multi layer perceptron" << endl;
INFO << "\t\tThe TMVA weight file is :" << endl;
INFO << "\t\t\t" << fTMVAWeightPath << endl;
- EquivalenceModel::PrintInfo();
}
//________________________________________________________________________
-TTree* EQM_PWR_MLP_MOX::CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
+TTree* EQM_PWR_MLP_MOX::CreateTMVAInputTree(map < string , IsotopicVector> IVStream, double BurnUp)
{
TTree* InputTree = new TTree("EQTMP", "EQTMP");
float Pu8 = 0;
@@ -90,40 +90,40 @@ TTree* EQM_PWR_MLP_MOX::CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector
float Pu11 = 0;
float Pu12 = 0;
float Am1 = 0;
- float U5_enrichment = 0;
- float BU = 0;
-
- InputTree->Branch( "Pu8" ,&Pu8 ,"Pu8/F" );
- InputTree->Branch( "Pu9" ,&Pu9 ,"Pu9/F" );
- InputTree->Branch( "Pu10" ,&Pu10 ,"Pu10/F" );
- InputTree->Branch( "Pu11" ,&Pu11 ,"Pu11/F" );
- InputTree->Branch( "Pu12" ,&Pu12 ,"Pu12/F" );
- InputTree->Branch( "Am1" ,&Am1 ,"Am1/F" );
+ float U5_enrichment = 0;
+ float BU = 0;
+
+ InputTree->Branch( "Pu8" ,&Pu8 ,"Pu8/F" );
+ InputTree->Branch( "Pu9" ,&Pu9 ,"Pu9/F" );
+ InputTree->Branch( "Pu10" ,&Pu10 ,"Pu10/F" );
+ InputTree->Branch( "Pu11" ,&Pu11 ,"Pu11/F" );
+ InputTree->Branch( "Pu12" ,&Pu12 ,"Pu12/F" );
+ InputTree->Branch( "Am1" ,&Am1 ,"Am1/F" );
InputTree->Branch( "U5_enrichment" ,&U5_enrichment ,"U5_enrichment/F" );
- InputTree->Branch( "BU" ,&BU ,"BU/F" );
+ InputTree->Branch( "BU" ,&BU ,"BU/F" );
- float U8 = Fertil.GetZAIIsotopicQuantity(92,238,0);
- float U5 = Fertil.GetZAIIsotopicQuantity(92,235,0);
- float U4 = Fertil.GetZAIIsotopicQuantity(92,234,0);
+ float U8 = IVStream["Fertile"].GetZAIIsotopicQuantity(92,238,0);
+ float U5 = IVStream["Fertile"].GetZAIIsotopicQuantity(92,235,0);
+ float U4 = IVStream["Fertile"].GetZAIIsotopicQuantity(92,234,0);
float UTOT = U8 + U5 + U4;
- Pu8 = Fissil.GetZAIIsotopicQuantity(94,238,0);
- Pu9 = Fissil.GetZAIIsotopicQuantity(94,239,0);
- Pu10 = Fissil.GetZAIIsotopicQuantity(94,240,0);
- Pu11 = Fissil.GetZAIIsotopicQuantity(94,241,0);
- Pu12 = Fissil.GetZAIIsotopicQuantity(94,242,0);
- Am1 = Fissil.GetZAIIsotopicQuantity(95,241,0);
+ Pu8 = IVStream["Fissile"].GetZAIIsotopicQuantity(94,238,0);
+ Pu9 = IVStream["Fissile"].GetZAIIsotopicQuantity(94,239,0);
+ Pu10 = IVStream["Fissile"].GetZAIIsotopicQuantity(94,240,0);
+ Pu11 = IVStream["Fissile"].GetZAIIsotopicQuantity(94,241,0);
+ Pu12 = IVStream["Fissile"].GetZAIIsotopicQuantity(94,242,0);
+ Am1 = IVStream["Fissile"].GetZAIIsotopicQuantity(95,241,0);
double TOTPU = (Pu8+Pu9+Pu10+Pu11+Pu12+Am1);
- Pu8 = Pu8 / TOTPU;
- Pu9 = Pu9 / TOTPU;
- Pu10 = Pu10 / TOTPU;
- Pu11 = Pu11 / TOTPU;
- Pu12 = Pu12 / TOTPU;
- Am1 = Am1 / TOTPU;
+ Pu8 = Pu8 / TOTPU ;
+ Pu9 = Pu9 / TOTPU ;
+ Pu10 = Pu10 / TOTPU ;
+ Pu11 = Pu11 / TOTPU ;
+ Pu12 = Pu12 / TOTPU ;
+ Am1 = Am1 / TOTPU ;
U5_enrichment = U5 / UTOT;
@@ -139,47 +139,51 @@ TTree* EQM_PWR_MLP_MOX::CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector
return InputTree;
}
//________________________________________________________________________
-double EQM_PWR_MLP_MOX::ExecuteTMVA(TTree* theTree)
+map < string , double> EQM_PWR_MLP_MOX::ExecuteTMVA(TTree* theTree)
{
+ map < string , double> MolarFraction;
// --- Create the Reader object
TMVA::Reader *reader = new TMVA::Reader( "Silent" );
// Create a set of variables and declare them to the reader
// - the variable names MUST corresponds in name and type to those given in the weight file(s) used
Float_t Pu8,Pu9,Pu10,Pu11,Pu12,Am1,BU,U5_enrichment;
- reader->AddVariable( "BU" ,&BU );
- reader->AddVariable( "U5_enrichment",&U5_enrichment );
- reader->AddVariable( "Pu8" ,&Pu8 );
- reader->AddVariable( "Pu9" ,&Pu9 );
- reader->AddVariable( "Pu10" ,&Pu10);
- reader->AddVariable( "Pu11" ,&Pu11);
- reader->AddVariable( "Pu12" ,&Pu12);
- reader->AddVariable( "Am1" ,&Am1 );
+ reader->AddVariable( "BU" ,&BU );
+ reader->AddVariable( "U5_enrichment" ,&U5_enrichment );
+ reader->AddVariable( "Pu8" ,&Pu8 );
+ reader->AddVariable( "Pu9" ,&Pu9 );
+ reader->AddVariable( "Pu10" ,&Pu10 );
+ reader->AddVariable( "Pu11" ,&Pu11 );
+ reader->AddVariable( "Pu12" ,&Pu12 );
+ reader->AddVariable( "Am1" ,&Am1 );
// --- Book the MVA methods
// Book method MLP
TString methodName = "MLP method";
reader->BookMVA( methodName, fTMVAWeightPath );
- theTree->SetBranchAddress( "BU" ,&BU );
- theTree->SetBranchAddress( "U5_enrichment" ,&U5_enrichment ) ;
- theTree->SetBranchAddress( "Pu8" ,&Pu8 );
- theTree->SetBranchAddress( "Pu9" ,&Pu9 );
- theTree->SetBranchAddress( "Pu10" ,&Pu10 );
- theTree->SetBranchAddress( "Pu11" ,&Pu11 );
- theTree->SetBranchAddress( "Pu12" ,&Pu12 );
- theTree->SetBranchAddress( "Am1" ,&Am1 );
+ theTree->SetBranchAddress( "BU" ,&BU );
+ theTree->SetBranchAddress( "U5_enrichment" ,&U5_enrichment );
+ theTree->SetBranchAddress( "Pu8" ,&Pu8 );
+ theTree->SetBranchAddress( "Pu9" ,&Pu9 );
+ theTree->SetBranchAddress( "Pu10" ,&Pu10 );
+ theTree->SetBranchAddress( "Pu11" ,&Pu11 );
+ theTree->SetBranchAddress( "Pu12" ,&Pu12 );
+ theTree->SetBranchAddress( "Am1" ,&Am1 );
theTree->GetEntry(0);
Float_t val = (reader->EvaluateRegression( methodName ))[0];
delete reader;
delete theTree;
+
+ MolarFraction["Fissile"] = (double)val;
+ MolarFraction["Fertile"] = 1.- (double)val;
- return (double)val; //retourne teneur
+ return MolarFraction; //return Molar content of each component in the fuel
}
//________________________________________________________________________
-double EQM_PWR_MLP_MOX::GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
+map < string , double> EQM_PWR_MLP_MOX::GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp)
{DBGL
- return ExecuteTMVA(CreateTMVAInputTree(Fissil,Fertil,BurnUp));
+ return ExecuteTMVA(CreateTMVAInputTree(IVStream,BurnUp));
}
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.hxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.hxx
index f0a9bae11..976cc239b 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.hxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX.hxx
@@ -53,20 +53,13 @@ class EQM_PWR_MLP_MOX : public EquivalenceModel
\param Fertil : The composition of the Fertil matter
\param BurnUp : Maximum achievable burn up envisaged
*/
- virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);
+ map < string , double> GetMolarFraction(map < string , IsotopicVector> IVStream, double BurnUp);
//}
- /*!
- \name TMVA related methods
- */
- //@{
-
- TTree* CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);//!<Create input tmva tree to be read by ExecuteTMVA
- double ExecuteTMVA(TTree* theTree);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
-
- //@}
-
private :
+
+ TTree* CreateTMVAInputTree(map < string , IsotopicVector> IVStream, double BurnUp);//!<Create input tmva tree to be read by ExecuteTMVA
+ map < string , double> ExecuteTMVA(TTree* theTree);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
string fTMVAWeightPath;;//!<The weight needed by TMVA to construct and execute the multilayer perceptron
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.cxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.cxx
index d6e914c4a..4bb7537a8 100755
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.cxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.cxx
@@ -47,11 +47,13 @@ EQM_PWR_MLP_MOX_AM::EQM_PWR_MLP_MOX_AM(string TMVAWeightPath):EquivalenceModel(n
SetBuildFuelFirstGuess(0.04);
+ fStreamList.push_back(FissileList);
+ fStreamList.push_back(FertileList);
+
INFO << "__An equivalence model of PWR MOX has been define__" << endl;
INFO << "\tThis model is based on a multi layer perceptron" << endl;
INFO << "\t\tThe TMVA weight file is :" << endl;
INFO << "\t\t\t"<<fTMVAWeightPath << endl;
- EquivalenceModel::PrintInfo();
}
@@ -78,17 +80,22 @@ EQM_PWR_MLP_MOX_AM::EQM_PWR_MLP_MOX_AM(CLASSLogger* log, string TMVAWeightPath):
SetBuildFuelFirstGuess(0.04);
+ fStreamList.push_back(FissileList);
+ fStreamList.push_back(FertileList);
+
INFO << "__An equivalence model of PWR MOX has been define__" << endl;
INFO << "\tThis model is based on a multi layer perceptron" << endl;
INFO << "\t\tThe TMVA weight file is :" << endl;
INFO << "\t\t\t"<<fTMVAWeightPath << endl;
- EquivalenceModel::PrintInfo();
}
//________________________________________________________________________
-TTree* EQM_PWR_MLP_MOX_AM::CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
+TTree* EQM_PWR_MLP_MOX_AM::CreateTMVAInputTree(vector <IsotopicVector> IVStream,double BurnUp)
{
+ IsotopicVector Fissil = IVStream[0];
+ IsotopicVector Fertil = IVStream[1];
+
TTree* InputTree = new TTree("EQTMP", "EQTMP");
float Pu8 = 0;
float Pu9 = 0;
@@ -197,7 +204,7 @@ double EQM_PWR_MLP_MOX_AM::ExecuteTMVA(TTree* theTree)
return (double)val; //retourne teneur
}
//________________________________________________________________________
-double EQM_PWR_MLP_MOX_AM::GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp)
+double EQM_PWR_MLP_MOX_AM::GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp)
{DBGL
- return ExecuteTMVA(CreateTMVAInputTree(Fissil,Fertil,BurnUp));
+ return ExecuteTMVA(CreateTMVAInputTree(IVStream,BurnUp));
}
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.hxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.hxx
index f742c9220..eca891033 100755
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.hxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_MLP_MOX_Am.hxx
@@ -53,21 +53,14 @@ class EQM_PWR_MLP_MOX_AM : public EquivalenceModel
\param Fertil : The composition of the Fertil matter
\param BurnUp : Maximum achievable burn up envisaged
*/
- virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);
+ virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp);
//}
- /*!
- \name TMVA related methods
- */
- //@{
-
- TTree* CreateTMVAInputTree(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);//!<Create input tmva tree to be read by ExecuteTMVA
- double ExecuteTMVA(TTree* theTree);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
-
- //@}
-
private :
+ TTree* CreateTMVAInputTree(vector <IsotopicVector> IVStream,double BurnUp);//!<Create input tmva tree to be read by ExecuteTMVA
+ double ExecuteTMVA(TTree* theTree);//!<Execute the MLP according to the input tree created by CreateTMVAInputTree
+
string fTMVAWeightPath;;//!<The weight needed by TMVA to construct and execute the multilayer perceptron
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.cxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.cxx
index 0c25731ef..63810e0f4 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.cxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.cxx
@@ -25,7 +25,6 @@ EQM_PWR_POL_UO2::EQM_PWR_POL_UO2(string PathToWeightFile):EquivalenceModel(new C
fFissileList = U5*1;
ReadWeightFile(PathToWeightFile);
- EquivalenceModel::PrintInfo();
}
// _______________________________________________________________________
@@ -41,7 +40,6 @@ EQM_PWR_POL_UO2::EQM_PWR_POL_UO2(CLASSLogger* log,string PathToWeightFile):Equiv
fFissileList = U5*1;
ReadWeightFile(PathToWeightFile);
- EquivalenceModel::PrintInfo();
}
// _______________________________________________________________________
@@ -68,11 +66,8 @@ void EQM_PWR_POL_UO2::ReadWeightFile(string PathToWeightFile)
INFO << "\t U enrichment = " << fParam_Bu_0 << " + " << fParam_Bu << "*Burnup + " << fParam_BuSquare << "*Burnup*Burnup" << endl;
}
// _______________________________________________________________________
-double EQM_PWR_POL_UO2::GetFissileMolarFraction ( IsotopicVector Fissil , IsotopicVector Fertil , double BurnUp )
+double EQM_PWR_POL_UO2::GetFissileMolarFraction ( vector <IsotopicVector> IVStream , double BurnUp )
{
- double Fraction = fParam_Bu_0 + fParam_Bu*BurnUp + fParam_BuSquare*BurnUp*BurnUp;
- DBGV("Fissile molar fraction : "<<Fraction)
-
- return Fraction ;
+ return fParam_Bu_0 + fParam_Bu*BurnUp + fParam_BuSquare*BurnUp*BurnUp ;
}
\ No newline at end of file
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.hxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.hxx
index d6b855dd5..b56d1aab2 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.hxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_POL_UO2.hxx
@@ -53,7 +53,7 @@ public:
//@}
/**This function IS the equivalence model**/
- double GetFissileMolarFraction(IsotopicVector Fissil, IsotopicVector Fertil,double BurnUp) ; // !<Return the molar fraction of fissile element
+ double GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp) ; // !<Return the molar fraction of fissile element
private:
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.cxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.cxx
index 5f984d08f..8ff918b8d 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.cxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.cxx
@@ -45,8 +45,9 @@ EQM_PWR_QUAD_MOX::EQM_PWR_QUAD_MOX(string WeightPath):EquivalenceModel(new CLASS
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
SetBuildFuelFirstGuess(0.04);
- EquivalenceModel::PrintInfo();
-
+
+ fStreamList.push_back(FissileList);
+ fStreamList.push_back(FertileList);
}
@@ -87,7 +88,9 @@ EQM_PWR_QUAD_MOX::EQM_PWR_QUAD_MOX(CLASSLogger* log, string WeightPath):Equivale
fFissileList = Pu8*1+Pu9*1+Pu0*1+Pu1*1+Pu2*1;
SetBuildFuelFirstGuess(0.04);
- EquivalenceModel::PrintInfo();
+
+ fStreamList.push_back(FissileList);
+ fStreamList.push_back(FertileList);
}
@@ -101,9 +104,11 @@ EQM_PWR_QUAD_MOX::~EQM_PWR_QUAD_MOX()
-double EQM_PWR_QUAD_MOX::GetFissileMolarFraction(IsotopicVector Fissile,IsotopicVector Fertile,double BurnUp)
+double EQM_PWR_QUAD_MOX::GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp)
{
+ IsotopicVector Fissile = IVStream[0];
+ IsotopicVector Fertile = IVStream[1];
ZAI ZAIList[6] = {ZAI(94,238,0), ZAI(94,239,0), ZAI(94,240,0), ZAI(94,241,0), ZAI(94,242,0), ZAI(95,241,0) };
diff --git a/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.hxx b/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.hxx
index 94bbe4c65..5117faf4f 100644
--- a/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.hxx
+++ b/source/branches/SL3/Model/Equivalence/EQM_PWR_QUAD_MOX.hxx
@@ -57,7 +57,7 @@ class EQM_PWR_QUAD_MOX : public EquivalenceModel
~EQM_PWR_QUAD_MOX();
//@}
- virtual double GetFissileMolarFraction(IsotopicVector Fissil,IsotopicVector Fertil,double BurnUp);
+ virtual double GetFissileMolarFraction(vector <IsotopicVector> IVStream,double BurnUp);
private :
diff --git a/source/branches/SL3/Model/Irradiation/IM_RK4.cxx b/source/branches/SL3/Model/Irradiation/IM_RK4.cxx
index 0aaa26893..f0fec721d 100644
--- a/source/branches/SL3/Model/Irradiation/IM_RK4.cxx
+++ b/source/branches/SL3/Model/Irradiation/IM_RK4.cxx
@@ -40,7 +40,7 @@ IM_RK4::IM_RK4():IrradiationModel(new CLASSLogger("IM_RK4.log")), DynamicalSyste
}
-//________________________________________________________________________
+
IM_RK4::IM_RK4(CLASSLogger* log):IrradiationModel(log), DynamicalSystem()
{
fTheNucleiVector = 0;
@@ -51,6 +51,13 @@ IM_RK4::IM_RK4(CLASSLogger* log):IrradiationModel(log), DynamicalSystem()
}
+
+
+
+
+
+
+
//________________________________________________________________________
/* Evolution Calculation */
//________________________________________________________________________
@@ -105,9 +112,9 @@ EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, Evolu
}
- //-------------------------//
- //--- Perform Evolution ---//
- //-------------------------//
+ //---------------------------------//
+ //--- Perform Evolution ----//
+ //--------------------------------//
ReactorType = XSSet.GetReactorType();
double M_ref = XSSet.GetHeavyMetalMass();
@@ -148,7 +155,7 @@ EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, Evolu
vector< TMatrixT<double> > FissionXSMatrix; // Store The Fisison XS Matrix
vector< TMatrixT<double> > CaptureXSMatrix; // Store The Capture XS Matrix
- vector< TMatrixT<double> > n2nXSMatrix; // Store The n2N XS Matrix
+ vector< TMatrixT<double> > n2nXSMatrix; // Store The n2N XS Matrix
DBGL
for(int i = 0; i < NStep-1; i++)
{
@@ -241,11 +248,10 @@ EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, Evolu
}
GeneratedDB.NucleiInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NMatrix.size(), timevector, ZAIQuantity)));
- /* GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, FissionXS)));
+/* GeneratedDB.FissionXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, FissionXS)));
GeneratedDB.CaptureXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, CaptureXS)));
GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, n2nXS)));
- */
- }
+*/ }
DBGL
GeneratedDB.SetPower(Power );
GeneratedDB.SetHeavyMetalMass(M);
@@ -269,7 +275,7 @@ EvolutionData IM_RK4::GenerateEvolutionData(IsotopicVector isotopicvector, Evolu
}
-//________________________________________________________________________
+
void IM_RK4::ResetTheMatrix()
{
@@ -282,7 +288,6 @@ void IM_RK4::ResetTheMatrix()
fTheMatrix = 0;
}
-//________________________________________________________________________
void IM_RK4::SetTheMatrixToZero()
{
ResetTheMatrix();
@@ -373,3 +378,19 @@ TMatrixT<double> IM_RK4::GetTheNucleiVector()
return NEvolutionMatrix;
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/source/branches/SL3/src/EquivalenceModel.cxx b/source/branches/SL3/src/EquivalenceModel.cxx
index 460e5f285..6047e5111 100644
--- a/source/branches/SL3/src/EquivalenceModel.cxx
+++ b/source/branches/SL3/src/EquivalenceModel.cxx
@@ -5,28 +5,27 @@
//________________________________________________________________________
EquivalenceModel::EquivalenceModel():CLASSObject()
{
- fRelativMassPrecision = 5/10000.; // Mass precision
- fMaxInterration = 100; // Max iterration in build fueld algorythum
- fFirstGuessFissilContent = 0.02;
- freaded = false;
+ fRelativMassPrecision = 5/10000.; //Mass precision
+ fMaxInterration = 10000; // Max iterration in build fueld algorythum
+ freaded = false;
+
EquivalenceModel::LoadKeyword();
}
-//________________________________________________________________________
+
EquivalenceModel::EquivalenceModel(CLASSLogger* log):CLASSObject(log)
{
- fRelativMassPrecision = 5/10000.; // Mass precision
- fMaxInterration = 100; // Max iterration in build fueld algorythm
- fFirstGuessFissilContent = 0.02;
- freaded = false;
+ fRelativMassPrecision = 5/10000.; //Mass precision
+ fMaxInterration = 10000; // Max iterration in build fueld algorythm
+ freaded = false;
+
EquivalenceModel::LoadKeyword();
-
}
-//________________________________________________________________________
+
EquivalenceModel::~EquivalenceModel()
{
-
+
}
-//________________________________________________________________________
+
void EquivalenceModel::ReadNFO()
{
DBGL
@@ -77,12 +76,12 @@ void EquivalenceModel::ReadLine(string line)
void EquivalenceModel::LoadKeyword()
{
DBGL
- fKeyword.insert( pair<string, EQM_MthPtr>( "k_zail", & EquivalenceModel::ReadZAIlimits));
- fKeyword.insert( pair<string, EQM_MthPtr>( "k_reactor", & EquivalenceModel::ReadType) );
- fKeyword.insert( pair<string, EQM_MthPtr>( "k_fuel", & EquivalenceModel::ReadType) );
- fKeyword.insert( pair<string, EQM_MthPtr>( "k_fissil", & EquivalenceModel::ReadFissil) );
- fKeyword.insert( pair<string, EQM_MthPtr>( "k_fertil", & EquivalenceModel::ReadFertil) );
- fKeyword.insert( pair<string, EQM_MthPtr>( "k_specpower", & EquivalenceModel::ReadSpecificPower));
+ fKeyword.insert( pair<string, EQM_MthPtr>( "k_zail", & EquivalenceModel::ReadZAIlimits) );
+ fKeyword.insert( pair<string, EQM_MthPtr>( "k_reactor", & EquivalenceModel::ReadType) );
+ fKeyword.insert( pair<string, EQM_MthPtr>( "k_fuel", & EquivalenceModel::ReadType) );
+ fKeyword.insert( pair<string, EQM_MthPtr>( "k_firstguesscontent", & EquivalenceModel::ReadFirstGuessContent) );
+ fKeyword.insert( pair<string, EQM_MthPtr>( "k_list", & EquivalenceModel::ReadList) );
+ fKeyword.insert( pair<string, EQM_MthPtr>( "k_specpower", & EquivalenceModel::ReadSpecificPower) );
DBGL
}
//________________________________________________________________________
@@ -91,19 +90,23 @@ void EquivalenceModel::PrintInfo()
INFO << "Reactor Type : "<< fDBRType << endl;
INFO << "Fuel Type : "<< fDBFType << endl;
INFO << "Specific Power [W/g]: "<< fSpecificPower << endl;
+
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
- INFO << "Fissile Liste (Z A I) :" << endl;
- map<ZAI ,double >::iterator it1;
- map<ZAI ,double > fMap1 = fFissileList.GetIsotopicQuantity();
- for(it1 = fMap1.begin() ; it1 != fMap1.end() ; it1++)
- INFO << (*it1).first.Z() <<" "<< (*it1).first.A() <<" "<< (*it1).first.I() << endl;
-
- INFO<<"Fertile Liste (Z A I Quantity) :"<<endl;
- map<ZAI ,double >::iterator it2;
- map<ZAI ,double > fMap2 = fFertileList.GetIsotopicQuantity();
- for(it2 = fMap2.begin() ; it2 != fMap2.end() ; it2++)
- INFO << (*it2).first.Z() <<" "<< (*it2).first.A() << " "<< (*it2).first.I() <<" "<< (*it2).second << endl;
-
+ for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++)
+ {
+ INFO <<(* it_s_IV).first<<" (Z A I) :" << endl;
+ map<ZAI ,double >::iterator it1;
+ map<ZAI ,double > fMap1 = fStreamList[(* it_s_IV).first].GetIsotopicQuantity();
+ for(it1 = fMap1.begin() ; it1 != fMap1.end() ; it1++)
+ INFO << (*it1).first.Z() <<" "<< (*it1).first.A() <<" "<< (*it1).first.I() << endl;
+ }
+ INFO<<"First guess content in fuel : "<<endl;
+ for( it_s_D = fFirstGuessContent.begin(); it_s_D != fFirstGuessContent.end(); it_s_D++)
+ {
+ INFO <<(* it_s_D).first<<" "<<fFirstGuessContent[(* it_s_D).first]<<endl;
+ }
INFO<<"ZAI limits (validity domain)[prop in fresh fuel] (Z A I min max) :"<<endl;
for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++)
@@ -148,61 +151,56 @@ void EquivalenceModel::ReadZAIlimits(const string &line)
exit(1);
}
- int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- double downLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
- double upLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
+ double downLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
+ double upLimit = atof(StringLine::NextWord(line, pos, ' ').c_str());
if (upLimit < downLimit)
{
- double tmp = upLimit;
- upLimit = downLimit;
- downLimit = tmp;
+ double tmp = upLimit;
+ upLimit = downLimit;
+ downLimit = tmp;
}
fZAILimits.insert(pair<ZAI, pair<double, double> >(ZAI(Z,A,I), pair<double,double>(downLimit, upLimit)));
DBGL
}
//________________________________________________________________________
-void EquivalenceModel::ReadFissil(const string &line)
+void EquivalenceModel::ReadList(const string &line)
{
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if( keyword != "k_fissil" ) // Check the keyword
+ if( keyword != "k_list" ) // Check the keyword
{
- ERROR << " Bad keyword : \"k_fissil\" not found !" << endl;
+ ERROR << " Bad keyword : \"k_list\" not found !" << endl;
exit(1);
}
-
- int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
-
- fFissileList.Add(Z, A, I, 1.0);
+ string ListName = StringLine::NextWord(line, pos, ' ');
+ int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
+ double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
+ fStreamList[ListName].Add(Z, A, I, Q);
DBGL
}
//________________________________________________________________________
-void EquivalenceModel::ReadFertil(const string &line)
+void EquivalenceModel::ReadFirstGuessContent(const string &line)
{
DBGL
int pos = 0;
string keyword = tlc(StringLine::NextWord(line, pos, ' '));
- if( keyword != "k_fertil" ) // Check the keyword
+ if( keyword != "k_firstguesscontent" ) // Check the keyword
{
- ERROR << " Bad keyword : \"k_fertil\" not found !" << endl;
+ ERROR << " Bad keyword : \"k_firstguesscontent\" not found !" << endl;
exit(1);
}
-
- int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
- double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
-
- fFertileList.Add(Z, A, I, Q);
-
+ string ListName = StringLine::NextWord(line, pos, ' ');
+ double Q = atof(StringLine::NextWord(line, pos, ' ').c_str());
+ fFirstGuessContent[ListName] = Q;
DBGL
}
@@ -222,144 +220,206 @@ void EquivalenceModel::ReadSpecificPower(const string &line)
DBGL
}
+
//________________________________________________________________________
-double EquivalenceModel::LAMBDA_TOT_FOR(double MassNeeded, vector<IsotopicVector> Stocks, string FisOrFer)
+void EquivalenceModel::StocksTotalMassCalculation(map < string , vector <IsotopicVector> > Stocks)
{
- double Lambda_tot = 0;
-
// Calculating total mass of stock once and for all
- if( fTotalFissileMassInStocks == 0 || fTotalFertileMassInStocks == 0 )
- {
- double TotalMassInStocks = 0;
- for( int i = 0 ; i<(int)Stocks.size() ; i++ )
- TotalMassInStocks += Stocks[i].GetTotalMass() ;
-
- if(FisOrFer == "Fis")
- fTotalFissileMassInStocks = TotalMassInStocks * 1e6; // in grams
- else
- fTotalFertileMassInStocks = TotalMassInStocks * 1e6; // in grams
- }
- double TotalMassInStocks = 0;
-
- if(FisOrFer == "Fis")
- TotalMassInStocks = fTotalFissileMassInStocks;
- else
- TotalMassInStocks = fTotalFertileMassInStocks;
- // If there is not enought matter in stocks construction fails
- if( MassNeeded > TotalMassInStocks )
+ double TotalMassInStocks = 0;
+ map < string , vector <IsotopicVector> >::iterator it_s_vIV;
+
+ for( it_s_vIV = Stocks.begin(); it_s_vIV != Stocks.end(); it_s_vIV++)
{
- WARNING << "Not enought " << FisOrFer << " material to build fuel" << endl;
- WARNING << TotalMassInStocks << endl;
- return -1;
+ fTotalMassInStocks[(*it_s_vIV).first] = 0;
+ fLambdaMax[(*it_s_vIV).first] = 0;
+
}
-
- for( int i = 0 ; i<(int)Stocks.size() ; i++ )
- {
-
- if( MassNeeded >= (Stocks[i].GetTotalMass()*1e6) )
+
+ for( it_s_vIV = Stocks.begin(); it_s_vIV != Stocks.end(); it_s_vIV++)
+ {
+ TotalMassInStocks = 0;
+ for(int i=0; i<Stocks[(* it_s_vIV).first].size(); i++)
{
- Lambda_tot+= 1;
- MassNeeded -= (Stocks[i].GetTotalMass()*1e6);
+ TotalMassInStocks += Stocks[(* it_s_vIV).first][i].GetTotalMass();
}
- else
+ fLambdaMax[(*it_s_vIV).first] = Stocks[(* it_s_vIV).first].size();
+ fTotalMassInStocks[(* it_s_vIV).first] = TotalMassInStocks * 1e6; // in grams
+ }
+}
+
+//________________________________________________________________________
+double EquivalenceModel::LambdaCalculation(string MaterialDenomination, double LambdaPreviousStep, double MaterialMassNeeded, double DeltaMass, vector <IsotopicVector> Stocks)
+{
+ double Lambda_tot = 0;
+
+ // If there is not enough matter in stocks construction fails
+ if( MaterialMassNeeded > fTotalMassInStocks[MaterialDenomination] )
+ {
+ if(DeltaMass > 0)
+ Lambda_tot = LambdaPreviousStep - (fLambdaMax[MaterialDenomination] - LambdaPreviousStep)/2 ;
+
+ if(DeltaMass < 0)
+ Lambda_tot = LambdaPreviousStep + (fLambdaMax[MaterialDenomination] - LambdaPreviousStep)/2 ;
+
+ if( (fLambdaMax[MaterialDenomination] - Lambda_tot)<1 && (fLambdaMax[MaterialDenomination]-Lambda_tot)*Stocks.back().GetTotalMass()*1e6<MaterialMassNeeded*fRelativMassPrecision/2.)
{
- Lambda_tot+= MassNeeded/(Stocks[i].GetTotalMass()*1e6);
- break;
+ WARNING << "Not enough " << MaterialDenomination << " material to build fuel" << endl;
+ WARNING << "Mass available "<<fTotalMassInStocks[MaterialDenomination] << endl;
+ WARNING << "Mass needed "<<MaterialMassNeeded<< endl;
+ Lambda_tot = -1;
}
}
+ // Lambda calculation
+ else
+ {
+ for( int i=0; i<Stocks.size(); i++)
+ {
+ if( MaterialMassNeeded >= (Stocks[i].GetTotalMass()*1e6))
+ {
+ Lambda_tot += 1;
+ MaterialMassNeeded -= (Stocks[i].GetTotalMass()*1e6);
+ }
+ else
+ {
+ Lambda_tot += MaterialMassNeeded/(Stocks[i].GetTotalMass()*1e6);
+ break;
+ }
+ }
+ }
+
return Lambda_tot;
}
//________________________________________________________________________
-bool EquivalenceModel::Build_Fuel_According_Lambda(vector<double> &lambda,vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray, double HMMass, IsotopicVector &Fissile, IsotopicVector &Fertile)
+void EquivalenceModel::SetLambda(vector<double>& lambda, double Lambda_tot)
{
- //Build the Plutonium vector from stocks
- Fissile.Clear();
- for( int i = 0; i < (int)FissilArray.size(); i++ )
- Fissile += lambda[i] * FissilArray[i];
-
-
- double AvailablePuMass = Fissile.GetTotalMass() * 1e6; // in grams
-
- // Building complementary Fertile from stocks
- double FertilMassNeeded = HMMass - AvailablePuMass;
- double LAMBDA_FERTILE = LAMBDA_TOT_FOR( FertilMassNeeded , FertilArray , "Fer" );
-
- 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++)
+ if(Lambda_tot > (int)lambda.size() )
{
- j++;
- Fertile += lambda[i] * FertilArray[j];
+ cout<<Lambda_tot<<" "<<lambda.size()<<endl;
+ ERROR << " FATAL ERROR " <<endl;
+ exit(0);
}
-
- if( fabs(Fertile.GetTotalMass()*1e6 - FertilMassNeeded) > FertilMassNeeded * 1e-6) // Not enought fertile in stocks
- {
- WARNING << "Not enought fertile material to build fuel" << endl;
- return false;
+
+ for(int i = 0 ; i < (int)lambda.size() ; i++) //set to 0 all non touched value (to be sure)
+ lambda[i] = 0 ;
+
+ int IntegerPart = floor( Lambda_tot );
+ double DecimalPart = Lambda_tot - IntegerPart;
+
+ for(int i=0 ; i < IntegerPart; i++ )
+ lambda[i]=1;
+
+ lambda[IntegerPart] = DecimalPart;
+}
+
+//________________________________________________________________________
+void EquivalenceModel::SetLambdaToErrorCode(vector<double>& lambda)
+{
+ for( int i=0; i<lambda.size(); i++)
+ {
+ lambda[i] = -1;
}
-
- return true;
}
+
//________________________________________________________________________
-vector<double> EquivalenceModel::BuildFuel(double BurnUp, double HMMass, vector<IsotopicVector> FissilArray, vector<IsotopicVector> FertilArray)
+map <string , vector<double> > EquivalenceModel::BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray)
{
+DBGL
+
+ map <string , vector<double> > lambda ; // map containing name of the list and associated vector of proportions taken from stocks
- DBGL
- vector<double> lambda ; // vector of portion of stocks taken (fissile & fertil)
- for(int i = 0 ; i < (int)FissilArray.size() + (int)FertilArray.size() ; i++ )
- lambda.push_back(0);
+ //Iterators declaration
+ map < string , vector <IsotopicVector> >::iterator it_s_vIV;
+ map < string , vector <double> >::iterator it_s_vD;
+ map < string , IsotopicVector >::iterator it_s_IV;
+ map < string , double >::iterator it_s_D;
+ map < string , bool >::iterator it_s_B;
- /*** Test if there is a stock **/
- if( (int)FissilArray.size() == 0 )
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i<StreamArray[(*it_s_vIV).first].size(); i++)
+ {
+ lambda[(*it_s_vIV).first].push_back(0);
+ }
+ }
+
+ /*** Test if there is stocks **/
+ bool BreakReturnLambda = false;
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
{
- WARNING << " No fissile stocks available ! Fuel not build" << endl;
- SetLambdaToErrorCode(lambda);
- return lambda;
+ if(StreamArray[(*it_s_vIV).first].size() == 0)
+ {
+ WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<". Fuel not build." << endl;
+ SetLambdaToErrorCode(lambda[(*it_s_vIV).first]);
+ BreakReturnLambda = true;
+ }
}
-
- HMMass *= 1e6; // Unit onversion : tons to gram
+ if(BreakReturnLambda) { return lambda;}
+ HMMass *= 1e6; //Unit onversion : tons to gram
/**** Some initializations **/
- fTotalFissileMassInStocks = 0;
- fTotalFertileMassInStocks = 0;
-
- fActualFissileContent = GetBuildFuelFirstGuess();
-
- IsotopicVector Fertile;
- IsotopicVector Fissile;
-
- double AvailablePuMass = 0;
- double PuMassNeeded = HMMass * fActualFissileContent;
- double WeightPuContent = 0;
- int loopCount = 0;
+
+ StocksTotalMassCalculation(StreamArray);
+
+ fActualMassContentInFuel = GetBuildFuelFirstGuess();
+ fActualMolarContentInFuel = GetBuildFuelFirstGuess();
+ int loopCount = 0;
+ bool FuelBuiltCorrectly = false ;
+
+ map <string , IsotopicVector > IVStream;
+ map < string , double > MaterialMassNeeded ;
+ map < string , double > LambdaNeeded;
+ map < string , double > DeltaMass;
- do
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ LambdaNeeded[(*it_s_vIV).first] = 0;
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
{
- double LAMBDA_NEEDED = LAMBDA_TOT_FOR(PuMassNeeded,FissilArray,"Fis");
- if( LAMBDA_NEEDED == -1 ) // Check if previous lambda was well calculated
- {
- SetLambdaToErrorCode(lambda);
- WARNING << "Not enought fissile material to build fuel" << endl;
- return lambda;
+ MaterialMassNeeded[(*it_s_vIV).first] = HMMass*fActualMassContentInFuel[(*it_s_vIV).first];
+ DeltaMass[(*it_s_vIV).first] = - MaterialMassNeeded[(*it_s_vIV).first];
+ }
+ do
+ {
+ map < string , double > LambdaPreviousStep;
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ LambdaPreviousStep[(*it_s_D).first] =LambdaNeeded[(*it_s_D).first];
+
+ BreakReturnLambda = false;
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ LambdaNeeded[(*it_s_vIV).first] = LambdaCalculation((*it_s_vIV).first, LambdaPreviousStep[(*it_s_vIV).first], MaterialMassNeeded[(*it_s_vIV).first], DeltaMass[(*it_s_vIV).first], StreamArray[(*it_s_vIV).first]);
+
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ {
+ if( LambdaNeeded[(*it_s_D).first] == -1 )
+ {
+ SetLambdaToErrorCode(lambda[(*it_s_D).first]);
+ WARNING << "Not enough : "<< (*it_s_D).first <<" to build fuel." << endl;
+ BreakReturnLambda = true;
+ }
}
- SetLambda(lambda, 0, FissilArray.size()-1, LAMBDA_NEEDED );
+ if(BreakReturnLambda) { return lambda;}
- bool succeed = Build_Fuel_According_Lambda(lambda, FissilArray, FertilArray, HMMass, Fissile, Fertile);
+ BreakReturnLambda = false;
- if(!succeed)
- {
- SetLambdaToErrorCode(lambda);
- return lambda;
+ for( it_s_D = LambdaNeeded.begin(); it_s_D != LambdaNeeded.end(); it_s_D++)
+ SetLambda(lambda[(*it_s_D).first], LambdaNeeded[(*it_s_D).first]);
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ IVStream[(*it_s_vIV).first].Clear();
+
+ for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++)
+ {
+ for(int i=0; i<StreamArray[(*it_s_vIV).first].size(); i++)
+ {
+ IVStream[(*it_s_vIV).first] += lambda[(*it_s_vIV).first][i] * StreamArray[(*it_s_vIV).first][i];
+ }
}
- AvailablePuMass = Fissile.GetTotalMass() * 1e6; //in grams
-
if (loopCount > fMaxInterration)
{
ERROR << "Too much iterration in BuildFuel Method !";
@@ -368,67 +428,83 @@ vector<double> EquivalenceModel::BuildFuel(double BurnUp, double HMMass, vector<
}
/* Calcul the quantity of this composition needed to reach the burnup */
- double MolarPuContent = GetFissileMolarFraction(Fissile, Fertile, BurnUp);
+
+ map < string , double> MaterialMolarContent = GetMolarFraction(IVStream, BurnUp);
- if( MolarPuContent < 0 || MolarPuContent > 1 )
- { SetLambdaToErrorCode(lambda);
- WARNING << "GetFissileMolarFraction return negative or greater than one value";
- return lambda;
+ for( it_s_D = MaterialMolarContent.begin(); it_s_D != MaterialMolarContent.end(); it_s_D++)
+ {
+ if( MaterialMolarContent[(*it_s_D).first] < 0 || MaterialMolarContent[(*it_s_D).first] > 1 )
+ {
+ SetLambdaToErrorCode(lambda[(*it_s_D).first]);
+ WARNING << "GetFissileMolarFraction return negative or greater than one value, at least for one material : "<<(*it_s_D).first;
+ BreakReturnLambda = true;
+ }
}
+ if(BreakReturnLambda) { return lambda;}
- double MeanMolarPu = Fissile.GetMeanMolarMass();
- double MeanMolarDepletedU = Fertile.GetMeanMolarMass();
-
- double MeanMolar = MeanMolarPu * MolarPuContent + (1-MolarPuContent) * MeanMolarDepletedU;
+ map < string , double > MaterialMassContent;
+ map < string , double > MaterialMeanMolar;
+ map < string , double > MaterialMassAvailable;
+ map < string , bool > CheckOnMass;
+
+ double FuelMeanMolar = 0;
+
+ for( it_s_IV = IVStream.begin(); it_s_IV != IVStream.end(); it_s_IV++)
+ MaterialMassAvailable[(*it_s_IV).first] = IVStream[(*it_s_IV).first].GetTotalMass()*1e06;
+ for( it_s_D = MaterialMolarContent.begin(); it_s_D != MaterialMolarContent.end(); it_s_D++)
+ {
+ MaterialMeanMolar[(*it_s_D).first] = IVStream[(*it_s_D).first].GetMeanMolarMass();
+ FuelMeanMolar += MaterialMolarContent[(*it_s_D).first] * MaterialMeanMolar[(*it_s_D).first];
+ }
- WeightPuContent = MolarPuContent * MeanMolarPu / MeanMolar;
- fActualFissileContent = MolarPuContent; //fActualFissileContent can be accessed by a derivated EquivalenModel to accelerate GetFissileMolarFraction function (exemple in EQM_MLP_Kinf)
- PuMassNeeded = WeightPuContent * HMMass ;
+ for( it_s_D = MaterialMolarContent.begin(); it_s_D != MaterialMolarContent.end(); it_s_D++)
+ MaterialMassContent [(*it_s_D).first] = MaterialMolarContent[(*it_s_D).first] * MaterialMeanMolar[(*it_s_D).first] / FuelMeanMolar;
+
+ fActualMolarContentInFuel = MaterialMolarContent; //fActualContent can be accessed by a derivated EquivalenModel to accelerate GetFissileMolarFraction function (exemple in EQM_MLP_Kinf)
- DBGV( "MolarPuContent " << MolarPuContent << " DeltaM " << PuMassNeeded - AvailablePuMass << " g" );
+ for( it_s_D = MaterialMassContent.begin(); it_s_D != MaterialMassContent.end(); it_s_D++)
+ {
+ MaterialMassNeeded[(*it_s_D).first] = HMMass*MaterialMassContent[(*it_s_D).first];
+ DeltaMass[(*it_s_D).first] = MaterialMassAvailable[(*it_s_D).first] - MaterialMassNeeded[(*it_s_D).first];
+ }
+ for( it_s_D = MaterialMassNeeded.begin(); it_s_D != MaterialMassNeeded.end(); it_s_D++)
+ {
+ double DeltaM = fabs(MaterialMassNeeded[(*it_s_D).first] - MaterialMassAvailable[(*it_s_D).first]) / HMMass ;
+
+ if(DeltaM<fRelativMassPrecision) {CheckOnMass[(*it_s_D).first] = true;}
+ else{CheckOnMass[(*it_s_D).first] = false;}
+ }
+
+ FuelBuiltCorrectly = true;
+ for( it_s_B = CheckOnMass.begin(); it_s_B != CheckOnMass.end(); it_s_B++)
+ {
+ if(!CheckOnMass[(*it_s_B).first]) {FuelBuiltCorrectly = false;}
+ }
+
loopCount++;
- }while( fabs( PuMassNeeded - AvailablePuMass )/HMMass > fRelativMassPrecision );
-
- (*this).isIVInDomain(Fissile);
-
- DBGV( "Weight percent fissile : " << PuMassNeeded/HMMass );
- DBGV( "Lambda vector : " );
+ }while(!FuelBuiltCorrectly);
+
+ for( it_s_IV = IVStream.begin(); it_s_IV != IVStream.end(); it_s_IV++)
+ (*this).isIVInDomain(IVStream[(*it_s_IV).first]);
- for(int i = 0; i < (int)FissilArray.size() + (int)FertilArray.size(); i++ )
- DBGV(lambda[i]);
+ for( it_s_D = MaterialMassNeeded.begin(); it_s_D != MaterialMassNeeded.end(); it_s_D++)
+ DBGV( "Weight percent : "<<(*it_s_D).first<<" "<< MaterialMassNeeded[(*it_s_D).first]/HMMass);
+ for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++)
+ {
+ DBGV( "Lambda vector : "<<(*it_s_vD).first );
+
+ for(int i=0; i<lambda[(*it_s_vD).first].size(); i++)
+ {
+ DBGV(lambda[(*it_s_vD).first][i]);
+ }
+ }
return lambda;
}
-//________________________________________________________________________
-void EquivalenceModel::SetLambda(vector<double>& lambda ,int FirstStockID, int LastStockID, double LAMBDA_TOT)
-{
- if( LAMBDA_TOT > LastStockID - FirstStockID + 1 )
- {
- ERROR << " FATAL ERROR " << endl;
- exit(0);
- }
-
- for( int i = FirstStockID; i <= LastStockID; i++ ) //set to 0 all non touched value (to be sure)
- lambda[i] = 0 ;
-
- int IntegerPart = floor( LAMBDA_TOT );
- double DecimalPart = LAMBDA_TOT - IntegerPart;
-
- for( int i = FirstStockID; i < FirstStockID +IntegerPart; i++ )
- lambda[i] = 1;
-
- lambda[FirstStockID + IntegerPart] = DecimalPart;
-}
-//________________________________________________________________________
-void EquivalenceModel::SetLambdaToErrorCode(vector<double>& lambda)
-{
- for(int i = 0 ; i < (int)lambda.size() ;i++ )
- lambda[i] = -1;
-}
//________________________________________________________________________
bool EquivalenceModel::isIVInDomain(IsotopicVector IV)
{
@@ -448,7 +524,7 @@ bool EquivalenceModel::isIVInDomain(IsotopicVector IV)
IsotopicVector IVNorm = IV /IV.GetSumOfAll();
for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++)
{
- double ThatZAIProp = IVNorm.GetIsotopicQuantity()[Domain_it->first] ;
+ double ThatZAIProp = IVNorm.GetIsotopicQuantity()[Domain_it->first];
double ThatZAIMin = Domain_it->second.first;
double ThatZAIMax = Domain_it->second.second;
if( (ThatZAIProp > ThatZAIMax) || (ThatZAIProp < ThatZAIMin) )
diff --git a/source/branches/SL3/src/EvolutionData.cxx b/source/branches/SL3/src/EvolutionData.cxx
index 853166f40..2d7d08327 100755
--- a/source/branches/SL3/src/EvolutionData.cxx
+++ b/source/branches/SL3/src/EvolutionData.cxx
@@ -661,7 +661,6 @@ void EvolutionData::ReadDB(string DBfile, bool oldread)
{
ERROR << " Bad Database file : " << DBfile << endl;
ERROR << " The first Line MUST be the time line !!!" << endl;
- ERROR << " last line red : \n" <<line<< endl;
exit (1);
}
@@ -1229,86 +1228,6 @@ void EvolutionData::OldReadDB(string DBfile)
}
-//________________________________________________________________________
-void EvolutionData::Print(string filename)
-{
-
- map<ZAI ,TGraph* >::iterator iterator;
- ofstream out(filename.c_str());
-
-
- out<<"time ";
- for(int t = 0 ; t < fInventoryEvolution.begin()->second->GetN() ; t++ )
- out<<fInventoryEvolution.begin()->second->GetX()[t]<< " ";
- out<<endl;
-
- if(fFlux)
- {
- out<<"flux ";
- for(int t = 0 ; t < fFlux->GetN() ; t++ )
- out<<fFlux->GetY()[t]<< " ";
- out<<endl;
- }
-
- if(fKeff)
- { out<<"keff ";
- for(int t = 0 ; t < fKeff->GetN() ; t++ )
- out<<fKeff->GetY()[t]<< " ";
- out<<endl;
- }
- for( iterator = fInventoryEvolution.begin(); iterator != fInventoryEvolution.end(); iterator++)
- {
- int N = (*iterator).second->GetN();
- if((*iterator).second->GetY()[N-1] != 0)
- { out<<"Inv "<<(*iterator).first.Z()<<" "<<(*iterator).first.A()<<" "<<(*iterator).first.I()<<" ";
- for(int t = 0 ; t < N ; t++ )
- out<<(*iterator).second->GetY()[t]<< " ";
- out<<endl;
- }
- }
-
- for( iterator = fFissionXS.begin(); iterator != fFissionXS.end(); iterator++)
- {
- int N = (*iterator).second->GetN();
- if((*iterator).second->GetY()[N-1] != 0)
- { out<<"XSFis "<<(*iterator).first.Z()<<" "<<(*iterator).first.A()<<" "<<(*iterator).first.I()<<" ";
-
- for(int t = 0 ; t < (*iterator).second->GetN() ; t++ )
- out<<-(*iterator).second->GetY()[t]*1e24<< " ";
-
- out<<endl;
- }
- }
-
- for( iterator = fCaptureXS.begin(); iterator != fCaptureXS.end(); iterator++)
- {
- int N = (*iterator).second->GetN();
- if((*iterator).second->GetY()[N-1] != 0)
- { out<<"XSCap "<<(*iterator).first.Z()<<" "<<(*iterator).first.A()<<" "<<(*iterator).first.I()<<" ";
-
- for(int t = 0 ; t < (*iterator).second->GetN() ; t++ )
- out<<-(*iterator).second->GetY()[t]*1e24<< " ";
-
- out<<endl;
- }
- }
-
- for( iterator = fn2nXS.begin(); iterator != fn2nXS.end(); iterator++)
- {
- int N = (*iterator).second->GetN();
- if((*iterator).second->GetY()[N-1] != 0)
- { out<<"XSn2n "<<(*iterator).first.Z()<<" "<<(*iterator).first.A()<<" "<<(*iterator).first.I()<<" ";
-
- for(int t = 0 ; t < (*iterator).second->GetN() ; t++ )
- out<<-(*iterator).second->GetY()[t]*1e24<< " ";
-
- out<<endl;
- }
- }
-
-
-
-}
diff --git a/source/branches/SL3/src/FabricationPlant.cxx b/source/branches/SL3/src/FabricationPlant.cxx
index 8bdbb0201..55295a1f2 100644
--- a/source/branches/SL3/src/FabricationPlant.cxx
+++ b/source/branches/SL3/src/FabricationPlant.cxx
@@ -3,15 +3,13 @@
#include "Storage.hxx"
#include "Reactor.hxx"
#include "EvolutionData.hxx"
+#include "DecayDataBank.hxx"
#include "PhysicsModels.hxx"
#include "IsotopicVector.hxx"
#include "Scenario.hxx"
#include "CLASSLogger.hxx"
#include "CLASSConstante.hxx"
-
-
-
#include "TMatrixT.h"
#include <sstream>
@@ -19,24 +17,18 @@
#include <iostream>
#include <cmath>
#include <algorithm>
-#include <ctime>
-#include <cstdlib>
+#include <stdarg.h>
ClassImp(FabricationPlant)
- //________________________________________________________________________
//________________________________________________________________________
//________________________________________________________________________
//
// FabricationPlant
- //
- //
- //
- //
//________________________________________________________________________
//________________________________________________________________________
-//________________________________________________________________________
+
FabricationPlant::FabricationPlant():CLASSFacility(16)
{
SetName("F_FabricationPLant.");
@@ -45,37 +37,34 @@ FabricationPlant::FabricationPlant():CLASSFacility(16)
fIsReusable = false;
}
-//________________________________________________________________________
+
FabricationPlant::FabricationPlant(CLASSLogger* log, double fabricationtime):CLASSFacility(log, fabricationtime, 16)
{
DBGL
SetName("F_FabricationPLant.");
- fStorageManagement = kpLiFo;
- fIsSeparationManagement = true;
+ fFiFo = false;
fSubstitutionFuel = false;
- fSubstitutionFissile = false;
- fIsReplaceFissileStock = false;
fReUsable = 0;
fIsReusable = false;
-
- INFO << " A FabricationPlant has been define :" << endl;
- INFO << "\t Chronological Stock Priority has been set! " << endl;
- INFO << "\t Fabrication time set to \t " << (double)(GetCycleTime()/cYear) << " year" << endl << endl;
+
+ INFO << " A FabricationPlant has been define :" << endl;
+ INFO << "\t Chronological Stock Priority has been set! "<< endl;
+ INFO << "\t Fabrication time set to \t " << (double)(GetCycleTime()/cYear) << " year" << endl << endl;
DBGL
}
-//________________________________________________________________________
+
+ //________________________________________________________________________
FabricationPlant::~FabricationPlant()
{
}
-
-//________________________________________________________________________
-void FabricationPlant::SetSeparartionEfficiencyIV(ZAI zai, double factor)
+ //________________________________________________________________________
+void FabricationPlant::SetSeparartionEfficiencyIV(ZAI zai, double factor)
{
pair<map<ZAI, double>::iterator, bool> IResult;
@@ -83,23 +72,21 @@ void FabricationPlant::SetSeparartionEfficiencyIV(ZAI zai, double factor)
if(factor > 0)
{
- IResult = fSeparationLostFraction.GetIsotopicQuantity().insert( pair<ZAI ,double>(zai, 1 - factor));
+ IResult = fSeparationLostFraction.GetIsotopicQuantity().insert( pair<ZAI ,double>(zai, 1 - factor));
if(!IResult.second)
IResult.first->second = 1 - factor;
}
}
-//________________________________________________________________________
-//_______________________________ Evolution ______________________________
-//________________________________________________________________________
-
-//________________________________________________________________________
+ //________________________________________________________________________
+ //_______________________________ Evolution ______________________________
+ //________________________________________________________________________
void FabricationPlant::Evolution(cSecond t)
{
// Check if the FabricationPlant has been created ...
- if(t == fInternalTime && t != 0) return;
+ if(t == fInternalTime && t != 0) return;
// Make the evolution for the FabricationPlant ...
FabricationPlantEvolution(t);
//Update Inside IsotopicVector
@@ -109,40 +96,39 @@ void FabricationPlant::Evolution(cSecond t)
}
-//________________________________________________________________________
+ //________________________________________________________________________
void FabricationPlant::FabricationPlantEvolution(cSecond t)
{
DBGL
+
map<int ,cSecond >::iterator it;
- for( it = fReactorNextStep.begin(); it != fReactorNextStep.end(); it++ )
+ for( it = fReactorNextStep.begin(); it!= fReactorNextStep.end(); it++ )
{
double R_CreactionTime = GetParc()->GetReactor()[ (*it).first ]->GetCreationTime();
- double R_LifeTime = GetParc()->GetReactor()[ (*it).first ]->GetLifeTime();
+ double R_LifeTime = GetParc()->GetReactor()[ (*it).first ]->GetLifeTime();
int ReactorId = (*it).first;
pair<CLASSFuel, double> R_Fuel = GetParc()->GetReactor()[ReactorId]->GetFuelPlan()->GetFuelAt( t + GetCycleTime() );
- double R_BU = R_Fuel.second;
- double R_Power = GetParc()->GetReactor()[ReactorId]->GetPower();
- double R_HMMass = GetParc()->GetReactor()[ReactorId]->GetHeavyMetalMass();
- cSecond R_CycleTime = (cSecond) (R_BU*1e9 / (R_Power) * R_HMMass * 3600 * 24);
-
-
+ double R_BU = R_Fuel.second;
+ double R_Power = GetParc()->GetReactor()[ReactorId]->GetPower();
+ double R_HMMass = GetParc()->GetReactor()[ReactorId]->GetHeavyMetalMass();
+ cSecond R_CycleTime = (cSecond) (R_BU*1e9 / (R_Power) * R_HMMass * 3600 * 24);
if( R_CycleTime < GetCycleTime())
{
- ERROR << "Reactor Cycle Time is shorter than Fabrication Time of the fuel, we cannot deal it upto now!!!" << endl;
+ ERROR << "Reactor Cycle Time is shorter than Fabrication Time of the fuel, we cannot deal it upto now!!!"<< endl;
exit(1);
}
if( t + GetCycleTime() >= R_CreactionTime
&& t + GetCycleTime() < R_CreactionTime + R_LifeTime)
{
- if( (*it).second == t )
+ if( (*it).second == t )
{
#pragma omp critical(FuelBuild)
{
if( R_Fuel.first.GetPhysicsModels() )
{
- BuildFuelForReactor( (*it).first, t );
+ BuildFuelForReactor( (*it).first, t );
}
(*it).second += R_CycleTime;
}
@@ -161,7 +147,6 @@ DBGL
DBGL
}
-//________________________________________________________________________
void FabricationPlant::UpdateInsideIV()
{
DBGL
@@ -174,74 +159,75 @@ void FabricationPlant::UpdateInsideIV()
DBGL
}
-//________________________________________________________________________
+
+ //________________________________________________________________________
void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
{
DBGL
- if(fFissileStorage.size() == 0)
- {
- ERROR << " One need at least one Fissile storage to build fuel " << endl;
- ERROR << " Use AddFissileStorage to add a stock to provide fissil material... " << endl;
- ERROR << " One need at least one Fissile storage to build fuel " << endl;
- exit(1);
- }
-
-
- double R_HM_Mass = GetParc()->GetReactor()[ ReactorId ]->GetHeavyMetalMass();
- double R_CycleTime = GetParc()->GetReactor()[ ReactorId ]->GetCycleTime();
- double R_Power = GetParc()->GetReactor()[ ReactorId ]->GetPower();
-
- pair<CLASSFuel, double > FuelBU = GetParc()->GetReactor()[ReactorId]->GetFuelPlan()->GetFuelAt(t+GetCycleTime()) ;
- PhysicsModels FuelType = *FuelBU.first.GetPhysicsModels();
- double R_BU = FuelBU.second;
+ double R_HM_Mass = GetParc()->GetReactor()[ ReactorId ]->GetHeavyMetalMass();
+ double R_CycleTime = GetParc()->GetReactor()[ ReactorId ]->GetCycleTime();
+ double R_Power = GetParc()->GetReactor()[ ReactorId ]->GetPower();
+
+ 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();
-
- // If there is not enough Fissile its possible to take from an infinite fissile stock
- if( !fIsReplaceFissileStock && fSubstitutionFissile )//if it is defined and wanted
- { IsotopicVector CooledSeparatedIV = GetDecay( fSubstitutionFissileIV , GetCycleTime());
- fFissileArray.push_back( CooledSeparatedIV/ CooledSeparatedIV.GetTotalMass() * R_HM_Mass );
+ map < string , vector <IsotopicVector> >::iterator it_s_vIV;
+ map < string , vector <double> >::iterator it_s_vD;
+ map < string , bool >::iterator it_s_B;
+
+ fStreamList = FuelType.GetEquivalenceModel()->GetAllStreamList();
+
+ BuildArray(ReactorId); // Checker chez les stocks si les StreamList sont présentes
+ // Grosse map qui contient tous les IV
+ // séparation + refroidissement virtuel
+ // Construit les stocks de matière infini (=taille du réacteur)
+
+ // string="MA, .." LambdaArray = tableau sur les IV
+ map < string , vector <double> > LambdaArray = FuelType.GetEquivalenceModel()->BuildFuel(R_BU, R_HM_Mass, fStreamArray);
+
+ bool FuelCanBeBuilt = true;
+ double LambdaSum = 0;
+
+ for( it_s_vD = LambdaArray.begin(); it_s_vD != LambdaArray.end(); it_s_vD++)
+ {
+ for(int i = 0; i < (int)LambdaArray[(*it_s_vD).first].size();i++) {fErrorOnLambda[(*it_s_vD).first] = false;}
}
-
- fFertileList = FuelType.GetEquivalenceModel()->GetFertileList();
-
-
- if(fFertileStorage.size() != 0) // If the fertile need to be taken in stock
- BuildFertileArray();
- else // if their is not stock and the fertile come from outside of the park
- {
- fFertileArray.push_back( fFertileList / fFertileList.GetTotalMass() * R_HM_Mass );
- DBGV("Fertile Array size : " << fFertileArray.size())
+ for( it_s_vD = LambdaArray.begin(); it_s_vD != LambdaArray.end(); it_s_vD++)
+ {
+ for(int i = 0; i < (int)LambdaArray[(*it_s_vD).first].size();i++)
+ {
+ if(LambdaArray[(*it_s_vD).first][i] == -1){fErrorOnLambda[(*it_s_vD).first] = true;}
+ LambdaSum += LambdaArray[(*it_s_vD).first][i];
+ }
}
-
- vector<double> LambdaArray = FuelType.GetEquivalenceModel()->BuildFuel(R_BU, R_HM_Mass, fFissileArray, fFertileArray);
-
- double LambdaSum = 0;
- for(int i = 0; i < (int)fFissileArray.size();i++)
- LambdaSum += LambdaArray[i];
+ for( it_s_B = fErrorOnLambda.begin(); it_s_B != fErrorOnLambda.end(); it_s_B++)
+ {
+ if(fErrorOnLambda[(*it_s_B).first]){FuelCanBeBuilt = false;}
+ }
- if(LambdaArray[0] != -1 && LambdaSum > 0 )
- {
- DBGV("Building process suceed: ")
+ if(FuelCanBeBuilt && LambdaSum > 0 )
+ {
+ DBGV("Building process from initial stocks has succeeded : ")
+ IsotopicVector IV = BuildFuelFromEqModel(LambdaArray);
+ IsotopicVector LoadedIV = GetDecay(IV,fCycleTime);
- IsotopicVector IV = BuildFuelFromEqModel(LambdaArray);
- EvolutionData EvolDB = FuelType.GenerateEvolutionData( GetDecay(IV,fCycleTime), R_CycleTime, R_Power);
+ EvolutionData EvolDB = FuelType.GenerateEvolutionData(GetDecay(IV,fCycleTime), R_CycleTime, R_Power);
{
pair<map<int, IsotopicVector>::iterator, bool> IResult;
IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
if(!IResult.second)
- IResult.first->second = IV;
+ IResult.first->second = IV;
}
{
pair<map<int, EvolutionData>::iterator, bool> IResult;
IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,EvolDB) );
if(!IResult.second)
- IResult.first->second = EvolDB;
+ IResult.first->second = EvolDB;
}
fInsideIV += IV;
AddCumulativeIVIn(IV);
@@ -250,10 +236,90 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
}
else
{
- DBGV("Building process failed: ")
- if (!fSubstitutionFuel)
+ DBGV("Building process from initial stocks failed: ")
+ //If the building fails : possibility to take a materials from an infinite stock composed of a substitution IV defined by user
+ bool IsSubstitutionMaterials = false;
+
+ for( it_s_B = fErrorOnLambda.begin(); it_s_B != fErrorOnLambda.end(); it_s_B++)
+ {
+ if(fErrorOnLambda[(*it_s_B).first] && fSubstitutionMaterialFromIV[(*it_s_B).first])
+ IsSubstitutionMaterials = true;
+ else if(fErrorOnLambda[(*it_s_B).first] && !fSubstitutionMaterialFromIV[(*it_s_B).first])
+ IsSubstitutionMaterials = false;
+ }
+
+ if(IsSubstitutionMaterials)
+ {
+ DBGV("Using substitute : -> From infinite substitutions IV ")
+
+ //Make the user specified composition to decay the fabrication time
+
+ map < string , IsotopicVector> CooledSeparatedIV;
+
+ for( it_s_B = fSubstitutionMaterialFromIV.begin(); it_s_B != fSubstitutionMaterialFromIV.end(); it_s_B++)
+ CooledSeparatedIV[(*it_s_B).first] = GetDecay(fSubstitutionIV[(*it_s_B).first], GetCycleTime());
+
+ for( it_s_vIV = fStreamArray.begin(); it_s_vIV != fStreamArray.end(); it_s_vIV++)
+ {
+ if(fErrorOnLambda[(*it_s_vIV).first])
+ {
+ fStreamArray[(*it_s_vIV).first].clear();
+ fStreamArray[(*it_s_vIV).first].push_back(CooledSeparatedIV[(*it_s_vIV).first]);
+ }
+ }
+
+ //Building the fuel :
+ for( it_s_vD = LambdaArray.begin(); it_s_vD != LambdaArray.end(); it_s_vD++)
+ LambdaArray[(*it_s_vD).first].clear();
+
+ LambdaArray = FuelType.GetEquivalenceModel()->BuildFuel(R_BU, R_HM_Mass, fStreamArray);
+ IsotopicVector IV = BuildFuelFromEqModel(LambdaArray);
+
+ //Generating the EvolutionData
+ EvolutionData EvolDB = FuelType.GenerateEvolutionData(GetDecay(IV,fCycleTime), R_CycleTime, R_Power);
+ {
+ pair<map<int, IsotopicVector>::iterator, bool> IResult;
+ IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
+ if(!IResult.second)
+ IResult.first->second = IV;
+ }
+ {
+ pair<map<int, EvolutionData>::iterator, bool> IResult;
+ IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,EvolDB) );
+ if(!IResult.second)
+ IResult.first->second = EvolDB;
+ }
+ GetParc()->AddOutIncome(IV);
+ fInsideIV += IV;
+ AddCumulativeIVIn(IV);
+ DBGL
+ }
+ else if(fSubstitutionFuel)
{
- DBGV("Reactor not loaded ")
+ DBGV("Using substitute : -> From a fixed data base")
+ IsotopicVector IV = fSubstitutionEvolutionData.GetIsotopicVectorAt(0);
+ EvolutionData evolutiondb = fSubstitutionEvolutionData * R_HM_Mass / IV.GetTotalMass();
+
+ IV = IV* R_HM_Mass / IV.GetTotalMass();
+ {
+ pair<map<int, IsotopicVector>::iterator, bool> IResult;
+ IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
+ if(!IResult.second)
+ IResult.first->second = IV;
+ }
+ {
+ pair<map<int, EvolutionData>::iterator, bool> IResult;
+ IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,evolutiondb) );
+ if(!IResult.second)
+ IResult.first->second = evolutiondb;
+ }
+ GetParc()->AddOutIncome( IV );
+ fInsideIV += IV;
+ AddCumulativeIVIn(IV);
+ }
+ else
+ {
+ DBGV("No Alternative Solution. Reactor not loaded. ")
{
EvolutionData EmptyDB;
pair<map<int, EvolutionData>::iterator, bool> IResult;
@@ -269,66 +335,7 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
IResult.first->second = EmptyIV;
}
}
- else
- {
- DBGV("Using substitute : ")
- if(fSubstitutionFissile)//if the build fail possibility to take a fissile material from an infinite fissile stock (if infinite stock defined)
- {
- DBGV("->From an infinite stock")
- //make the user specified fissil composition to decay the fabrication time
- IsotopicVector CooledSeparatedIV = GetDecay(fSubstitutionFissileIV, GetCycleTime());
- //Building the fuel :
- double MolarFissileContent = FuelType.GetEquivalenceModel()->GetFissileMolarFraction(CooledSeparatedIV,fFertileList, R_BU);
- IsotopicVector BuiltFuel = MolarFissileContent*fSubstitutionFissileIV/fSubstitutionFissileIV.GetSumOfAll() + (1-MolarFissileContent)*fFertileList/fFertileList.GetSumOfAll();
- IsotopicVector IV = BuiltFuel/ BuiltFuel.GetTotalMass() * R_HM_Mass;
-
- //Generating the EvolutionData
- EvolutionData EvolDB = FuelType.GenerateEvolutionData( GetDecay(IV,fCycleTime), R_CycleTime, R_Power);
-
- {
- pair<map<int, IsotopicVector>::iterator, bool> IResult;
- IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
- if(!IResult.second)
- IResult.first->second = IV;
- }
- {
- pair<map<int, EvolutionData>::iterator, bool> IResult;
- IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,EvolDB) );
- if(!IResult.second)
- IResult.first->second = EvolDB;
- }
- GetParc()->AddOutIncome(IV);
- fInsideIV += IV;
- AddCumulativeIVIn(IV);
-
- DBGL
- }
- else
- {
- DBGV("->From a fixed data base")
- IsotopicVector IV = fSubstitutionEvolutionData.GetIsotopicVectorAt(0);
- EvolutionData evolutiondb = fSubstitutionEvolutionData * R_HM_Mass / IV.GetTotalMass();
-
- IV = IV* R_HM_Mass / IV.GetTotalMass();
- {
- pair<map<int, IsotopicVector>::iterator, bool> IResult;
- IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
- if(!IResult.second)
- IResult.first->second = IV;
- }
- {
- pair<map<int, EvolutionData>::iterator, bool> IResult;
- IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,evolutiondb) );
- if(!IResult.second)
- IResult.first->second = evolutiondb;
- }
- GetParc()->AddOutIncome( IV );
- fInsideIV += IV;
- AddCumulativeIVIn(IV);
- }
-
- }
DBGL
ResetArrays();
@@ -337,312 +344,141 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
DBGL
}
-//________________________________________________________________________
-void FabricationPlant::BuildFissileArray()
+void FabricationPlant::BuildArray(int ReactorId)
{
DBGL
-
- for(int i = 0; i < (int)fFissileStorage.size(); i++)
- {
-
- vector<IsotopicVector> IVArray = fFissileStorage[i]->GetIVArray();
+ double R_HM_Mass = GetParc()->GetReactor()[ ReactorId ]->GetHeavyMetalMass();
- for(int j = 0; j < (int)IVArray.size(); j++)
+ vector <IsotopicVector> StreamArray;
+ vector <cSecond> StreamArrayTime;
+ vector < pair<int,int> > StreamArrayAdress;
+
+ map < string , IsotopicVector>::iterator it;
+
+ for( it = fStreamList.begin(); it != fStreamList.end(); it++)
+ {
+ if(fInfiniteMaterialFromList[(*it).first])
{
-
- IsotopicVector SeparatedIV = Separation(IVArray[j], fFissileList).first;
-
- if(Norme(SeparatedIV) != 0)
- {
- IsotopicVector CooledSeparatedIV = GetDecay( SeparatedIV , GetCycleTime());
-
- fFissileArray.push_back( CooledSeparatedIV );
- fFissileArrayAdress.push_back( pair<int,int>(i,j) );
- fFissileArrayTime.push_back(fFissileStorage[i]->GetIVArrayArrivalTime()[j]);
- }
-
+ IsotopicVector IV = (*it).second / (*it).second.GetTotalMass() * R_HM_Mass;
+ StreamArray.push_back(IV);
+ StreamArrayAdress.push_back(pair<int,int>(0,0));
+ StreamArrayTime.push_back(0);
}
-
- }
-
- SortArray(0);
-DBGL
-}
-
-//________________________________________________________________________
-void FabricationPlant::BuildFertileArray()
-{
-DBGL
-
- for(int i = 0; i < (int)fFertileStorage.size(); i++)
- {
- vector<IsotopicVector> IVArray = fFertileStorage[i]->GetIVArray();
- for(int j = 0; j < (int)IVArray.size(); j++)
+ else
{
-
- IsotopicVector SeparatedIV = Separation(IVArray[j], fFertileList).first;
- if(Norme(SeparatedIV) != 0)
+ for(int j = 0; j < (int)fStorage[(*it).first].size(); j++)
{
- IsotopicVector CooledSeparatedIV = GetDecay( SeparatedIV , GetCycleTime());
-
- fFertileArray.push_back( CooledSeparatedIV );
- fFertileArrayAdress.push_back( pair<int,int>(i,j) );
- fFertileArrayTime.push_back(fFertileStorage[i]->GetIVArrayArrivalTime()[j]);
+ vector<IsotopicVector> IVArray = fStorage[(*it).first][j]->GetIVArray();
+ for(int k = 0; k < (int)IVArray.size(); k++)
+ {
+ IsotopicVector SeparatedIV = Separation(IVArray[k], (*it).second).first;
+ if(Norme(SeparatedIV) != 0)
+ {
+ IsotopicVector CooledSeparatedIV = GetDecay( SeparatedIV , GetCycleTime());
+ StreamArray.push_back( CooledSeparatedIV );
+ StreamArrayAdress.push_back( pair<int,int>(j,k) );
+ StreamArrayTime.push_back(fStorage[(*it).first][j]->GetIVArrayArrivalTime()[k]);
+ }
+ }
}
}
+ fStreamArray[(*it).first] = StreamArray; StreamArray.clear();
+ fStreamArrayAdress[(*it).first] = StreamArrayAdress; StreamArrayAdress.clear();
+ fStreamArrayTime[(*it).first] = StreamArrayTime; StreamArrayTime.clear();
}
- SortArray(1);
+ SortArray();
DBGL
}
-//________________________________________________________________________
-void FabricationPlant::SortArray(int i)
+void FabricationPlant::SortArray()
{
+ map < string , vector <IsotopicVector> >::iterator it_s_vIV;
+
+ for( it_s_vIV = fStreamArray.begin(); it_s_vIV != fStreamArray.end(); it_s_vIV++)
+ {
+ vector<IsotopicVector> IVArray;
+ vector<cSecond> TimeArray;
+ vector< pair<int,int> > AdressArray;
- vector<IsotopicVector> IVArray;
- vector<cSecond> TimeArray;
- vector< pair<int,int> > AdressArray;
-
- if(i == 0) //Fissile
- {
- IVArray = fFissileArray;
- TimeArray = fFissileArrayTime;
- AdressArray = fFissileArrayAdress;
- }
- else if (i == 1) //Fertile
- {
- IVArray = fFertileArray;
- TimeArray = fFertileArrayTime;
- AdressArray = fFertileArrayAdress;
-
- }
-
-
- switch(fStorageManagement)
- {
- case kpFiFo: SortFiFo(IVArray, TimeArray, AdressArray);
- break;
-
- case kpLiFo: SortLiFo(IVArray, TimeArray, AdressArray);
- break;
-
- case kpMix: SortMix(IVArray, TimeArray, AdressArray);
- break;
-
- case kpRand: SortRandom(IVArray, TimeArray, AdressArray);
- break;
-
- default:
- ERROR<<" Posibble Storage Management are"<<endl;
- ERROR<<" YourFabPlant->SetStorageManagement(key); //with key can be either"<<endl;
- ERROR<<"\tkFiFo : First In First Out (i.e the older storage first)"<<endl;
- ERROR<<"\tkLiFo : Last In First Out (i.e the youger storage first)"<<endl;
- ERROR<<"\tkMix : IVs are sorted that way : "<<"\n"<<"First: The younger , Second: The older, Third: The second younger ,4th : the second older ...."<<endl;
- ERROR<<"\tkRand : IVs order in storage is random"<<endl;
-
- exit(1);
- }
-
-
- if(i == 0) //Fissile
- {
- fFissileArray = IVArray;
- fFissileArrayTime = TimeArray;
- fFissileArrayAdress = AdressArray;
- }
- else if (i == 1) //Fertile
- {
- fFertileArray = IVArray;
- fFertileArrayTime = TimeArray;
- fFertileArrayAdress = AdressArray;
-
- }
-
-
-}
+ IVArray = fStreamArray[(*it_s_vIV).first] ;
+ TimeArray = fStreamArrayTime[(*it_s_vIV).first] ;
+ AdressArray = fStreamArrayAdress[(*it_s_vIV).first] ;
-//________________________________________________________________________
-void FabricationPlant::SortFiFo(vector<IsotopicVector> &IVArray, vector<cSecond> &TimeArray, vector< pair<int,int> > &AdressArray)
-{
- for(int j = 0; j < (int)TimeArray.size(); j++)
+ if(fFiFo)
{
- for (int k = j+1; k < (int)TimeArray.size(); k++)
+ for(int j = 0; j < (int)TimeArray.size(); j++)
{
- cSecond time_tmp = TimeArray[j];
- pair<int,int> Adress_tmp = AdressArray[j];
- IsotopicVector IV_tmp = IVArray[j];
-
- if(time_tmp > TimeArray[k])
+ for (int k = j+1; k < (int)TimeArray.size(); k++)
{
- TimeArray[j] = TimeArray[k];
- TimeArray[k] = time_tmp;
+ cSecond time_tmp = TimeArray[j];
+ pair<int,int> Adress_tmp = AdressArray[j];
+ IsotopicVector IV_tmp = IVArray[j];
- AdressArray[j] = AdressArray[k];
- AdressArray[k] = Adress_tmp;
+ if(time_tmp > TimeArray[k])
+ {
+ TimeArray[j] = TimeArray[k];
+ TimeArray[k] = time_tmp;
- IVArray[j] = IVArray[k];
- IVArray[k] = IV_tmp;
- }
+ AdressArray[j] = AdressArray[k];
+ AdressArray[k] = Adress_tmp;
+ IVArray[j] = IVArray[k];
+ IVArray[k] = IV_tmp;
+ }
+ }
}
}
-
-}
-
-//________________________________________________________________________
-void FabricationPlant::SortLiFo(vector<IsotopicVector> &IVArray, vector<cSecond> &TimeArray, vector< pair<int,int> > &AdressArray)
-{
- for(int j = 0; j < (int)TimeArray.size(); j++)
- {
- for (int k = j+1; k < (int)TimeArray.size(); k++)
+ else
+ {
+ for(int j = 0; j < (int)TimeArray.size(); j++)
{
- cSecond time_tmp = TimeArray[j];
- pair<int,int> Adress_tmp = AdressArray[j];
- IsotopicVector IV_tmp = IVArray[j];
-
- if(time_tmp < TimeArray[k])
+
+ for (int k = j+1; k < (int)TimeArray.size(); k++)
{
- TimeArray[j] = TimeArray[k];
- TimeArray[k] = time_tmp;
-
- AdressArray[j] = AdressArray[k];
- AdressArray[k] = Adress_tmp;
-
- IVArray[j] = IVArray[k];
- IVArray[k] = IV_tmp;
- }
+ cSecond time_tmp = TimeArray[j];
+ pair<int,int> Adress_tmp = AdressArray[j];
+ IsotopicVector IV_tmp = IVArray[j];
+ if(time_tmp < TimeArray[k])
+ {
+ TimeArray[j] = TimeArray[k];
+ TimeArray[k] = time_tmp;
+
+ AdressArray[j] = AdressArray[k];
+ AdressArray[k] = Adress_tmp;
+
+ IVArray[j] = IVArray[k];
+ IVArray[k] = IV_tmp;
+ }
+ }
}
}
-
+ fStreamArray[(*it_s_vIV).first] = IVArray;
+ fStreamArrayTime[(*it_s_vIV).first] = TimeArray;
+ fStreamArrayAdress[(*it_s_vIV).first] = AdressArray;
+ }
}
//________________________________________________________________________
-void FabricationPlant::SortMix(vector<IsotopicVector> &IVArray, vector<cSecond> &TimeArray, vector< pair<int,int> > &AdressArray)
-{
-
- //Sort by anti-chronoligical order (youger first)
- SortLiFo(IVArray, TimeArray, AdressArray);
- /*******Store it ******/
- vector<IsotopicVector> Saved_IVArray = IVArray ;
- vector<cSecond> Saved_TimeArray = TimeArray ;
- vector< pair<int,int> > Saved_AdressArray = AdressArray;
-
- int IVsize = (int)IVArray.size();
- /*******Then reset the vectors ******/
- IVArray.clear();
- TimeArray.clear();
- AdressArray.clear();
-
-
- int HalfSize = floor( (double)IVsize/2. );
-
- int old = IVsize;
-
- bool isYoung=true;//change to false to begin with an old isotopicvector
-
- int RemainingIV = -1;
-
-
- for(int young = 0 ; young < HalfSize ; young++)
- {
-
- if(isYoung)
- {
- IVArray.push_back(Saved_IVArray[young]);
- TimeArray.push_back(Saved_TimeArray[young]);
- AdressArray.push_back(Saved_AdressArray[young]);
-
- isYoung=!isYoung;
- RemainingIV = young+1; // +1 ? -> The next young will be the +1
- }
- if(!isYoung)
- {
- old--;
-
- IVArray.push_back(Saved_IVArray[old]);
- TimeArray.push_back(Saved_TimeArray[old]);
- AdressArray.push_back(Saved_AdressArray[old]);
-
- isYoung=!isYoung;
- RemainingIV = old-1; // -1 ? -> The next old will be the -1
-
- }
-
- }
-
- if( (double)IVsize/2. - (double)HalfSize > 0.0 ) //if odd number of isotopic vector : one isotopic vector is still missing add it @ the end
- {
-
- IVArray.push_back(Saved_IVArray[RemainingIV]);
- TimeArray.push_back(Saved_TimeArray[RemainingIV]);
- AdressArray.push_back(Saved_AdressArray[RemainingIV]);
-
- }
-
-
-}
-
+// Substitution Fuel
//________________________________________________________________________
-void FabricationPlant::SortRandom(vector<IsotopicVector> &IVArray, vector<cSecond> &TimeArray, vector< pair<int,int> > &AdressArray)
-{
- int SizeOfIVArray = IVArray.size();
-
-/*********Create a Random list of vector position**********/
- srand ( unsigned ( std::time(0) ) );
- vector<int> RandomPosition;
- for (int i=0 ; i < (int) SizeOfIVArray ; ++i)
- RandomPosition.push_back(i);
-
- random_shuffle(RandomPosition.begin(), RandomPosition.end());
-
-/*******Store old vectors ******/
- vector<IsotopicVector> Saved_IVArray = IVArray ;
- vector<cSecond> Saved_TimeArray = TimeArray ;
- vector< pair<int,int> > Saved_AdressArray = AdressArray;
-
-/*******Asign values ******/
- for (int i=0 ; i < (int) SizeOfIVArray ; ++i)
- {
- IVArray[i] = Saved_IVArray[RandomPosition[i]];
- TimeArray[i] = Saved_TimeArray[RandomPosition[i]];
- AdressArray[i] = Saved_AdressArray[RandomPosition[i]];
- }
-
-}
-
-//________________________________________________________________________
-void FabricationPlant::SetSubstitutionFuel(EvolutionData fuel, bool ReplaceTheStock)
+void FabricationPlant::SetSubstitutionFuel(EvolutionData fuel)
{
fSubstitutionFuel = true;
- fIsReplaceFissileStock = ReplaceTheStock;
double M0 = cZAIMass.GetMass( fuel.GetIsotopicVectorAt(0.).GetActinidesComposition() );
fSubstitutionEvolutionData = fuel / M0;
}
-//________________________________________________________________________
-void FabricationPlant::SetSubstitutionFissile(IsotopicVector IV)
-{
-
- fSubstitutionFuel = true;
- fSubstitutionFissile = true;
-
- fSubstitutionFissileIV = IV / IV.GetSumOfAll();
-
-}
//________________________________________________________________________
//_____________________________ Reactor & DB _____________________________
-//________________________________________________________________________
-
-//________________________________________________________________________
void FabricationPlant::TakeReactorFuel(int Id)
{
DBGL
@@ -655,7 +491,7 @@ DBGL
(*it2).second = IV;
map< int,EvolutionData >::iterator it = fReactorFuturDB.find(Id);
- (*it).second.DeleteEvolutionDataCopy();
+ (*it).second = EvolutionData();
UpdateInsideIV();
DBGL
@@ -668,56 +504,42 @@ EvolutionData FabricationPlant::GetReactorEvolutionDB(int ReactorId)
map< int,EvolutionData >::iterator it = fReactorFuturDB.find(ReactorId);
return (*it).second;
}
-
//________________________________________________________________________
//_______________________________ Storage ________________________________
//________________________________________________________________________
-IsotopicVector FabricationPlant::BuildFuelFromEqModel(vector<double> LambdaArray)
+IsotopicVector FabricationPlant::BuildFuelFromEqModel(map <string , vector<double> > LambdaArray)
{
DBGL
IsotopicVector BuildedFuel;
IsotopicVector Lost;
+
+ map < string , IsotopicVector>::iterator it;
- int StockCorrection = 0;
- if( !fIsReplaceFissileStock && fSubstitutionFissile)
+ for( it = fStreamList.begin(); it != fStreamList.end(); it++)
{
- StockCorrection = 1;
- BuildedFuel += fFissileArray.back()*LambdaArray[fFissileArray.size()-1];
- }
-
- for(int i = 0; i < (int)fFissileArray.size() - StockCorrection ; i++)
- {
- if(LambdaArray[i] != 0)
+ for(int i = 0; i < (int)fStreamArray[(*it).first].size(); i++)
{
- int Stor_N = fFissileArrayAdress[i].first;
- int IV_N = fFissileArrayAdress[i].second;
-
- pair<IsotopicVector, IsotopicVector> Separated_Lost;
- Separated_Lost = Separation( fFissileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFissileList );
- BuildedFuel += Separated_Lost.first;
- Lost += Separated_Lost.second;
-
- }
- }
+ if(fInfiniteMaterialFromList[(*it).first])
+ BuildedFuel += fStreamArray[(*it).first][i]*LambdaArray[(*it).first][i];
- if(fFertileStorage.size() != 0)
- {
- for(int i = fFissileArray.size(); i < (int)(fFertileArray.size()+fFissileArray.size()); i++)
- {
- if(LambdaArray[i] != 0)
+ else if (fSubstitutionMaterialFromIV[(*it).first] && fErrorOnLambda[(*it).first])
+ BuildedFuel += fStreamArray[(*it).first][i]*LambdaArray[(*it).first][i];
+ else
{
- int Stor_N = fFertileArrayAdress[i].first;
- int IV_N = fFertileArrayAdress[i].second;
-
- pair<IsotopicVector, IsotopicVector> Separated_Lost;
- Separated_Lost = Separation( fFertileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFertileList);
- BuildedFuel += Separated_Lost.first;
- Lost += Separated_Lost.second;
+ if(LambdaArray[(*it).first][i] != 0)
+ {
+ int Stor_N = fStreamArrayAdress[(*it).first][i].first;
+ int IV_N = fStreamArrayAdress[(*it).first][i].second;
+
+ pair<IsotopicVector, IsotopicVector> Separated_Lost;
+ Separated_Lost = Separation( fStorage[(*it).first][Stor_N]->GetIVArray()[IV_N]*LambdaArray[(*it).first][i], (*it).second );
+ BuildedFuel += Separated_Lost.first;
+ Lost += Separated_Lost.second;
+
+ }
}
}
}
- else
- BuildedFuel += fFertileArray[0]*LambdaArray.back();
if(fIsReusable)
fReUsable->AddIV(Lost);
@@ -731,84 +553,67 @@ DBGL
}
//________________________________________________________________________
-void FabricationPlant::DumpStock(vector<double> LambdaArray)
+void FabricationPlant::DumpStock(map <string , vector<double> > LambdaArray)
{
DBGL
- int StockCorrection = 0;
- if( !fIsReplaceFissileStock && fSubstitutionFissile)
- { StockCorrection = 1;
- GetParc()->AddOutIncome( fFissileArray.back()*LambdaArray[fFissileArray.size()-1] );
- }
- for(int i = 0; i < (int)fFissileArray.size() - StockCorrection; i++)
- {
- if(LambdaArray[i] != 0)
- {
- int Stor_N = fFissileArrayAdress[i].first;
- int IV_N = fFissileArrayAdress[i].second;
- fFissileStorage[Stor_N]->TakeFractionFromStock( IV_N, LambdaArray[i] );
- }
- }
- if(fFertileStorage.size() != 0)
- {
- for(int i = fFissileArray.size(); i < (int)(fFertileArray.size()+fFissileArray.size()); i++)
+
+ map < string , IsotopicVector>::iterator it;
+
+ for( it = fStreamList.begin(); it != fStreamList.end(); it++)
+ {
+ for(int i = 0; i < (int)fStreamArray[(*it).first].size(); i++)
{
- if(LambdaArray[i] != 0)
+ if(fInfiniteMaterialFromList[(*it).first])
{
- int Stor_N = fFertileArrayAdress[i].first;
- int IV_N = fFertileArrayAdress[i].second;
-
- fFertileStorage[Stor_N]->TakeFractionFromStock( IV_N, LambdaArray[i] );
+ GetParc()->AddOutIncome( fStreamArray[(*it).first][0]*LambdaArray[(*it).first][i] );
+ }
+
+ else
+ {
+ if(LambdaArray[(*it).first][i] != 0)
+ {
+ int Stor_N = fStreamArrayAdress[(*it).first][i].first;
+ int IV_N = fStreamArrayAdress[(*it).first][i].second;
+ fStorage[(*it).first][Stor_N]->TakeFractionFromStock( IV_N, LambdaArray[(*it).first][i] );
+ }
}
}
}
- else
- GetParc()->AddOutIncome( fFertileArray[0]*LambdaArray.back() );
ResetArrays();
-
DBGL
}
-
//________________________________________________________________________
void FabricationPlant::ResetArrays()
{
- //Clear the Building Array (Fissile and Fertile)
- fFissileArray.clear();
- fFissileArrayTime.clear();
- fFissileArrayAdress.clear();
- fFertileArray.clear();
- fFertileArrayTime.clear();
- fFertileArrayAdress.clear();
-
- fFertileList = fFissileList = IsotopicVector();
-}
+ //Clear the Building Array
+ map < string , IsotopicVector>::iterator it;
+
+ for( it = fStreamList.begin(); it != fStreamList.end(); it++)
+ {
+ fStreamList[(*it).first]= IsotopicVector();
+ fStreamArray[(*it).first].clear();
+ fStreamArrayTime[(*it).first].clear();
+ fStreamArrayAdress[(*it).first].clear();
+ }
+
+
+}
//________________________________________________________________________
pair<IsotopicVector, IsotopicVector> FabricationPlant::Separation(IsotopicVector isotopicvector, IsotopicVector ExtractedList)
{
DBGL
- IsotopicVector SeparatedPart;
- IsotopicVector LostPart;
-
- if(fIsSeparationManagement)
- {
//[0] = re-use ; [1] = waste
- IsotopicVector LostInReprocessing = isotopicvector.GetThisComposition(ExtractedList) * fSeparationLostFraction;
- SeparatedPart = isotopicvector.GetThisComposition(ExtractedList) - LostInReprocessing;
- LostPart = isotopicvector - SeparatedPart;
- }
- else
- {
- //[0] = re-use ; [1] = waste
- //IsotopicVector LostInReprocessing = isotopicvector.GetThisComposition(ExtractedList) * fSeparationLostFraction;
- SeparatedPart = isotopicvector;
- LostPart = isotopicvector - SeparatedPart;
- }
+ IsotopicVector LostInReprocessing = isotopicvector.GetThisComposition(ExtractedList) * fSeparationLostFraction;
+ IsotopicVector SeparatedPart = isotopicvector.GetThisComposition(ExtractedList) - LostInReprocessing;
+ IsotopicVector LostPart = isotopicvector - SeparatedPart;
DBGL
return pair<IsotopicVector, IsotopicVector> (SeparatedPart, LostPart);
}
+
//________________________________________________________________________
// Get Decay
//________________________________________________________________________
@@ -832,3 +637,9 @@ IsotopicVector FabricationPlant::GetDecay(IsotopicVector isotopicvector, cSecond
}
+void FabricationPlant::AddInfiniteStorage(string keyword)
+{
+ Storage* Stock;
+ fStorage[keyword].push_back(Stock);
+ fInfiniteMaterialFromList[keyword] = true;
+}
diff --git a/source/branches/SL3/src/Reactor.cxx b/source/branches/SL3/src/Reactor.cxx
index a37839c42..7a6d40dfa 100755
--- a/source/branches/SL3/src/Reactor.cxx
+++ b/source/branches/SL3/src/Reactor.cxx
@@ -6,10 +6,11 @@
#include "Storage.hxx"
#include "Scenario.hxx"
#include "CLASSConstante.hxx"
+#include "../../config/config.hxx"
#include <iostream>
#include <cmath>
-#include "../../config/config.hxx"
+
#include <typeinfo>
//________________________________________________________________________
@@ -174,7 +175,7 @@ Reactor::Reactor(CLASSLogger* log, PhysicsModels* fueltypeDB, FabricationPlant*
fFuelPlan = new CLASSFuelPlan(log);
fFuelPlan->AddFuel(creationtime, CLASSFuel(fueltypeDB), fBurnUp);
-
+ CheckListConsistency(fueltypeDB, fabricationplant);
INFO << " A Reactor has been define :" << endl;
INFO << "\t Fuel Composition is not fixed ! " << endl;
@@ -218,8 +219,9 @@ Reactor::Reactor(CLASSLogger* log, PhysicsModels* fueltypeDB,
fFuelPlan = new CLASSFuelPlan(log);
fFuelPlan->AddFuel(creationtime, CLASSFuel(fueltypeDB), fBurnUp);
-
-
+
+ CheckListConsistency(fueltypeDB, fabricationplant);
+
INFO << " A Reactor has been define :" << endl;
INFO << "\t Fuel Composition is not fixed ! " << endl;
INFO << "\t Creation time set at \t " << ((double)GetCreationTime())/((double)cYear) << " year" << endl;
@@ -460,7 +462,7 @@ void Reactor::Evolution(cSecond t)
cSecond EvolutionTime = t - fInternalTime; // Calculation of the evolution time (relativ)
- if( EvolutionTime + fInCycleTime == fCycleTime ) //End of Cycle
+ if( EvolutionTime + fInCycleTime == fCycleTime ) //End of Cycle
{
fIsAtEndOfCycle = true;
fInternalTime += EvolutionTime; // Update Internal Time
@@ -638,3 +640,74 @@ cSecond Reactor::GetNextCycleTime(cSecond time)
return LastCycle;
}
+void Reactor::CheckListConsistency(PhysicsModels* fueltypeDB, FabricationPlant* fabricationplant)
+{
+ //Get Lists containers defined in EqM and FP
+ map < string , IsotopicVector> StreamList = fueltypeDB->GetEquivalenceModel()->GetAllStreamList();
+ map < string , vector <Storage*> > Stocks = fabricationplant->GetAllStorage();
+
+ vector <string> StreamListEqM;
+ vector <string> StreamListFP;
+
+ vector <bool> CheckOnLists;
+ bool AreListsConsistent = true;
+
+ //Iterators
+ map < string , vector <Storage*> >::iterator it_s_vS;
+ map < string , IsotopicVector>::iterator it_s_IV;
+
+ for( it_s_vS = Stocks.begin(); it_s_vS != Stocks.end(); it_s_vS++)
+ StreamListFP.push_back((*it_s_vS).first);
+
+ for( it_s_IV = StreamList.begin(); it_s_IV != StreamList.end(); it_s_IV++)
+ StreamListEqM.push_back((*it_s_IV).first);
+
+ if((int)StreamListEqM.size() != (int)StreamListFP.size())
+ {
+ WARNING <<" Not the same number of lists in FP and associated EqM"<< endl;
+ AreListsConsistent = false;
+ }
+
+ else
+ {
+ for(int i=0; i<StreamListEqM.size(); i++)
+ {
+ CheckOnLists.push_back(false);
+ }
+
+ for(int i=0; i<StreamListEqM.size(); i++)
+ {
+ for(int j=0; j<StreamListFP.size(); j++)
+ {
+ if (StreamListEqM[i] == StreamListFP[j]) {CheckOnLists[i] = true;}
+ }
+ }
+
+ for(int i=0; i<CheckOnLists.size(); i++)
+ {
+ if (!CheckOnLists[i]) {AreListsConsistent = false;}
+ }
+ }
+
+ if (!AreListsConsistent)
+ {
+ ERROR <<"Lists defined in Fabrication plant and in EqM are not the same."<< endl;
+ ERROR <<"Lists defined in Fabrication plant :"<< endl;
+ for(int j=0; j<StreamListFP.size(); j++)
+ {
+ ERROR <<"-> "<<StreamListFP[j]<<" "<< endl;
+ }
+ ERROR <<"Lists defined in EqM :"<< endl;
+
+ for(int i=0; i<StreamListEqM.size(); i++)
+ {
+ ERROR <<"-> "<<StreamListEqM[i]<<" "<< endl;
+ }
+ ERROR <<"Check in your scenario."<< endl;
+
+ exit(1);
+
+ }
+
+
+}
diff --git a/source/branches/SL3/src/Scenario.cxx b/source/branches/SL3/src/Scenario.cxx
index 780e07c34..da917c4b1 100755
--- a/source/branches/SL3/src/Scenario.cxx
+++ b/source/branches/SL3/src/Scenario.cxx
@@ -774,7 +774,7 @@ void Scenario::PrintCLASSPresentation()
cout<<"├───────────────────────────────────────────────┤"<<endl;
cout<<"│ Authors : │"<<endl;
cout<<"│ B. MOUGINOT (@BaM) B. LENIAU (@BLG) │"<<endl;
- cout<<"│ F. COURTIN (@FaC) N. THIOLIERE (@NT) │"<<endl;
+ cout<<"│ F. COURTIN (@FaC) N. THIOLLIERE (@NT) │"<<endl;
cout<<"╰───────────────────────────────────────────────╯"<<endl;
cout<<" "<<endl;
cout<<" BEGINING FUEL CYCLE EVOLUTION "<<endl;
@@ -803,7 +803,7 @@ void Scenario::PrintClover(int i)
cout<<"├───────────────────────────────────────────────┤ @@@@@@@ "<<endl;
cout<<"│ Authors : │ @@@@@@@ "<<endl;
cout<<"│ B. MOUGINOT (@BaM) B. LENIAU (@BLG) │ @@@@@@ "<<endl;
- cout<<"│ F. COURTIN (@FaC) N. THIOLIERE (@NT) │ @@@@@ "<<endl;
+ cout<<"│ F. COURTIN (@FaC) N. THIOLLIERE (@NT) │ @@@@@ "<<endl;
cout<<"╰───────────────────────────────────────────────╯ @@@ "<<endl;
cout<<" "<<endl;
}
@@ -825,7 +825,7 @@ void Scenario::PrintClover(int i)
cout<<"├───────────────────────────────────────────────┤ @@@@@ @@@@@ "<<endl;
cout<<"│ Authors : │ @@@@ @@@@ "<<endl;
cout<<"│ B. MOUGINOT (@BaM) B. LENIAU (@BLG) │ @@ @@@ "<<endl;
- cout<<"│ F. COURTIN (@FaC) N. THIOLIERE (@NT) │ @ @ "<<endl;
+ cout<<"│ F. COURTIN (@FaC) N. THIOLLIERE (@NT) │ @ @ "<<endl;
cout<<"╰───────────────────────────────────────────────╯ "<<endl;
cout<<" "<<endl;
}
@@ -848,7 +848,7 @@ void Scenario::PrintClover(int i)
cout<<"├───────────────────────────────────────────────┤ @@@@@@@ "<<endl;
cout<<"│ Authors : │ @@@@@@@ "<<endl;
cout<<"│ B. MOUGINOT (@BaM) B. LENIAU (@BLG) │ @@@@@@@ "<<endl;
- cout<<"│ F. COURTIN (@FaC) N. THIOLIERE (@NT) │ @@@@@ "<<endl;
+ cout<<"│ F. COURTIN (@FaC) N. THIOLLIERE (@NT) │ @@@@@ "<<endl;
cout<<"╰───────────────────────────────────────────────╯ @@@@ "<<endl;
cout<<" "<<endl;
@@ -872,7 +872,7 @@ void Scenario::PrintClover(int i)
cout<<"├───────────────────────────────────────────────┤ @@@@@@ "<<endl;
cout<<"│ Authors : │ @@@@@@@ "<<endl;
cout<<"│ B. MOUGINOT (@BaM) B. LENIAU (@BLG) │ @@@@@@@@ "<<endl;
- cout<<"│ F. COURTIN (@FaC) N. THIOLIERE (@NT) │ "<<endl;
+ cout<<"│ F. COURTIN (@FaC) N. THIOLLIERE (@NT) │ "<<endl;
cout<<"╰───────────────────────────────────────────────╯ "<<endl;
cout<<" "<<endl;
diff --git a/source/branches/SL3/src/XSModel.cxx b/source/branches/SL3/src/XSModel.cxx
index cb6707c07..0caac5ebc 100644
--- a/source/branches/SL3/src/XSModel.cxx
+++ b/source/branches/SL3/src/XSModel.cxx
@@ -37,13 +37,12 @@ void XSModel::ReadNFO()
{
DBGL
ifstream NFO(fInformationFile.c_str());
-
+
if(!NFO)
{
ERROR << "Can't find/open file " << fInformationFile << endl;
exit(0);
}
-
do
{
string line;
@@ -70,7 +69,7 @@ void XSModel::ReadLine(string line)
if(it != fKeyword.end())
(this->*(it->second))( line );
-
+
freaded = true;
ReadLine(line);
@@ -86,7 +85,7 @@ void XSModel::LoadKeyword()
{
DBGL
fKeyword.insert( pair<string, XSM_MthPtr>( "k_zail", & XSModel::ReadZAIlimits));
- fKeyword.insert( pair<string, XSM_MthPtr>( "k_reactor", & XSModel::ReadType) );
+ fKeyword.insert( pair<string, XSM_MthPtr>( "k_reactor",& XSModel::ReadType) );
fKeyword.insert( pair<string, XSM_MthPtr>( "k_fuel", & XSModel::ReadType) );
fKeyword.insert( pair<string, XSM_MthPtr>( "k_mass", & XSModel::ReadRParam) );
fKeyword.insert( pair<string, XSM_MthPtr>( "k_power", & XSModel::ReadRParam) );
--
GitLab