From 858e6d4c805b1ce0beca09babfe7291547117a10 Mon Sep 17 00:00:00 2001 From: Baptiste LENIAU <baptiste.leniau@subatech.in2p3.fr> Date: Fri, 22 Apr 2016 10:38:05 +0000 Subject: [PATCH] git-svn-id: svn+ssh://svn.in2p3.fr/class@830 0e7d625b-0364-4367-a6be-d5be4a48d228 --- .../To_Merge/src/EquivalenceModel.cxx | 486 ++++++++++-------- .../branches/To_Merge/src/EvolutionData.cxx | 1 - .../To_Merge/src/FabricationPlant.cxx | 5 +- source/branches/To_Merge/src/Makefile | 6 +- source/branches/To_Merge/src/Reactor.cxx | 92 +++- source/branches/To_Merge/src/Scenario.cxx | 20 +- 6 files changed, 378 insertions(+), 232 deletions(-) diff --git a/source/branches/To_Merge/src/EquivalenceModel.cxx b/source/branches/To_Merge/src/EquivalenceModel.cxx index 460e5f285..26c268d4d 100644 --- a/source/branches/To_Merge/src/EquivalenceModel.cxx +++ b/source/branches/To_Merge/src/EquivalenceModel.cxx @@ -5,26 +5,25 @@ //________________________________________________________________________ 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() @@ -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 < (int)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 < (int)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 < (int)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 < (int)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 < (int)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 < (int)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/To_Merge/src/EvolutionData.cxx b/source/branches/To_Merge/src/EvolutionData.cxx index 853166f40..1864265c5 100755 --- a/source/branches/To_Merge/src/EvolutionData.cxx +++ b/source/branches/To_Merge/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); } diff --git a/source/branches/To_Merge/src/FabricationPlant.cxx b/source/branches/To_Merge/src/FabricationPlant.cxx index 4572a1706..7451941d6 100644 --- a/source/branches/To_Merge/src/FabricationPlant.cxx +++ b/source/branches/To_Merge/src/FabricationPlant.cxx @@ -359,7 +359,6 @@ void FabricationPlant::BuildArray(int ReactorId) vector < pair<int,int> > StreamArrayAdress; map < string , IsotopicVector>::iterator it; - for( it = fStreamList.begin(); it != fStreamList.end(); it++) { if(fInfiniteMaterialFromList[(*it).first]) @@ -370,7 +369,7 @@ void FabricationPlant::BuildArray(int ReactorId) StreamArrayTime.push_back(0); } else - { + { for(int j = 0; j < (int)fStorage[(*it).first].size(); j++) { vector<IsotopicVector> IVArray = fStorage[(*it).first][j]->GetIVArray(); @@ -387,7 +386,7 @@ void FabricationPlant::BuildArray(int ReactorId) } } } - + fStreamArray[(*it).first] = StreamArray; StreamArray.clear(); fStreamArrayAdress[(*it).first] = StreamArrayAdress; StreamArrayAdress.clear(); fStreamArrayTime[(*it).first] = StreamArrayTime; StreamArrayTime.clear(); diff --git a/source/branches/To_Merge/src/Makefile b/source/branches/To_Merge/src/Makefile index 277e7ee8e..dec557d1f 100755 --- a/source/branches/To_Merge/src/Makefile +++ b/source/branches/To_Merge/src/Makefile @@ -38,9 +38,9 @@ OBJS = CLASSLogger.o \ Scenario.o OBJMODEL = $(EQMINC)/EQM_PWR_MLP_MOX.o $(EQMINC)/EQM_PWR_MLP_MOX_Am.o \ - $(EQMINC)/EQM_PWR_QUAD_MOX.o $(EQMINC)/EQM_PWR_POL_UO2.o\ - $(EQMINC)/EQM_PWR_LIN_MOX.o $(EQMINC)/EQM_FBR_BakerRoss_MOX.o $(EQMINC)/EQM_MLP_Kinf.o\ - $(EQMINC)/EQM_FBR_MLP_Keff.o $(EQMINC)/EQM_FBR_MLP_Keff_BOUND.o\ + $(EQMINC)/EQM_PWR_POL_UO2.o\ + $(EQMINC)/EQM_FBR_BakerRoss_MOX.o $(EQMINC)/EQM_MLP_Kinf.o\ + $(EQMINC)/EQM_FBR_MLP_Keff.o \ $(XSMINC)/XSM_MLP.o $(XSMINC)/XSM_CLOSEST.o \ $(IMINC)/IM_RK4.o $(IMINC)/IM_Matrix.o diff --git a/source/branches/To_Merge/src/Reactor.cxx b/source/branches/To_Merge/src/Reactor.cxx index a37839c42..811ec1ffe 100755 --- a/source/branches/To_Merge/src/Reactor.cxx +++ b/source/branches/To_Merge/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> //________________________________________________________________________ @@ -37,7 +38,7 @@ Reactor::Reactor():CLASSFacility(4) fFuelPlan = 0; } - +//________________________________________________________________________ Reactor::Reactor(CLASSLogger* log):CLASSFacility(log, 4) { DBGL @@ -50,7 +51,7 @@ Reactor::Reactor(CLASSLogger* log):CLASSFacility(log, 4) DBGL } - +//________________________________________________________________________ Reactor::Reactor(CLASSLogger* log, CLASSBackEnd* Pool, cSecond creationtime, @@ -100,7 +101,7 @@ Reactor::Reactor(CLASSLogger* log, DBGL } - +//________________________________________________________________________ Reactor::Reactor(CLASSLogger* log, FabricationPlant* fabricationplant, CLASSBackEnd* Pool, cSecond creationtime, cSecond lifetime, @@ -143,8 +144,7 @@ Reactor::Reactor(CLASSLogger* log, DBGL } - - +//________________________________________________________________________ Reactor::Reactor(CLASSLogger* log, PhysicsModels* fueltypeDB, FabricationPlant* fabricationplant, CLASSBackEnd* Pool, cSecond creationtime, cSecond lifetime, double Power, double HMMass, double BurnUp, double CapacityFactor):CLASSFacility(log, creationtime, lifetime, 4) @@ -174,7 +174,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 +218,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 +461,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 +639,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 < (int)StreamListEqM.size(); i++) + { + CheckOnLists.push_back(false); + } + + for(int i=0; i < (int)StreamListEqM.size(); i++) + { + for(int j=0; j < (int)StreamListFP.size(); j++) + { + if (StreamListEqM[i] == StreamListFP[j]) {CheckOnLists[i] = true;} + } + } + + for(int i=0; i < (int)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 < (int)StreamListFP.size(); j++) + { + ERROR <<"-> "<<StreamListFP[j]<<" "<< endl; + } + ERROR <<"Lists defined in EqM :"<< endl; + + for(int i=0; i < (int)StreamListEqM.size(); i++) + { + ERROR <<"-> "<<StreamListEqM[i]<<" "<< endl; + } + ERROR <<"Check in your scenario."<< endl; + + exit(1); + + } + + +} diff --git a/source/branches/To_Merge/src/Scenario.cxx b/source/branches/To_Merge/src/Scenario.cxx index 780e07c34..c6812a580 100755 --- a/source/branches/To_Merge/src/Scenario.cxx +++ b/source/branches/To_Merge/src/Scenario.cxx @@ -764,7 +764,7 @@ void Scenario::PrintCLASSPresentation() cout<<"│ ██║ ██║ ██╔â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘ │"<<endl; cout<<"│ ╚██████╗███████╗██║ ██║███████║███████║ │"<<endl; cout<<"│ â•šâ•â•â•â•â•â•â•šâ•â•â•â•â•â•â•â•šâ•â• â•šâ•â•â•šâ•â•â•â•â•â•â•â•šâ•â•â•â•â•â•â• │"<<endl; - cout<<"│ Version 4.1 │"<<endl; + cout<<"│ Version dev │"<<endl; cout<<"├───────────────────────────────────────────────┤"<<endl; cout<<"│ Core Lybrary for Advance Scenario Simulation │"<<endl; cout<<"│ │"<<endl; @@ -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; @@ -793,7 +793,7 @@ void Scenario::PrintClover(int i) cout<<"│ ██║ ██║ ██╔â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘ │ @@@@@@@ "<<endl; cout<<"│ ╚██████╗███████╗██║ ██║███████║███████║ │ @@@@@@@ "<<endl; cout<<"│ â•šâ•â•â•â•â•â•â•šâ•â•â•â•â•â•â•â•šâ•â• â•šâ•â•â•šâ•â•â•â•â•â•â•â•šâ•â•â•â•â•â•â• │ @@@@@@ @ "<<endl; - cout<<"│ Version 4.1 │ @@@@ @@@@ "<<endl; + cout<<"│ Version dev │ @@@@ @@@@ "<<endl; cout<<"├───────────────────────────────────────────────┤ @@ @@@@@@ "<<endl; cout<<"│ Core Lybrary for Advance Scenario Simulation │ @ @@@@@@@ "<<endl; cout<<"│ │ @@ @@@@@@@ "<<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; } @@ -815,7 +815,7 @@ void Scenario::PrintClover(int i) cout<<"│ ██║ ██║ ██╔â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘ │ @@@@@@ "<<endl; cout<<"│ ╚██████╗███████╗██║ ██║███████║███████║ │ @@@@@ "<<endl; cout<<"│ â•šâ•â•â•â•â•â•â•šâ•â•â•â•â•â•â•â•šâ•â• â•šâ•â•â•šâ•â•â•â•â•â•â•â•šâ•â•â•â•â•â•â• │ @@@@ "<<endl; - cout<<"│ Version 4.1 │ @@ "<<endl; + cout<<"│ Version dev │ @@ "<<endl; cout<<"├───────────────────────────────────────────────┤ @ "<<endl; cout<<"│ Core Lybrary for Advance Scenario Simulation │ @ "<<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; } @@ -838,7 +838,7 @@ void Scenario::PrintClover(int i) cout<<"│ ██║ ██║ ██╔â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘ │ @@@@@@@ "<<endl; cout<<"│ ╚██████╗███████╗██║ ██║███████║███████║ │ @@@@@@@@ "<<endl; cout<<"│ â•šâ•â•â•â•â•â•â•šâ•â•â•â•â•â•â•â•šâ•â• â•šâ•â•â•šâ•â•â•â•â•â•â•â•šâ•â•â•â•â•â•â• │ @@ @@@@@@ "<<endl; - cout<<"│ Version 4.1 │ @@@ @@@@ "<<endl; + cout<<"│ Version dev │ @@@ @@@@ "<<endl; cout<<"├───────────────────────────────────────────────┤ @@@@@@ @@ "<<endl; cout<<"│ Core Lybrary for Advance Scenario Simulation │ @@@@@@@@ @ "<<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; @@ -862,7 +862,7 @@ void Scenario::PrintClover(int i) cout<<"│ ██║ ██║ ██╔â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘â•šâ•â•â•â•â–ˆâ–ˆâ•‘ │ @@@@ @@@ "<<endl; cout<<"│ ╚██████╗███████╗██║ ██║███████║███████║ │ @@@@@ @@@@@ "<<endl; cout<<"│ â•šâ•â•â•â•â•â•â•šâ•â•â•â•â•â•â•â•šâ•â• â•šâ•â•â•šâ•â•â•â•â•â•â•â•šâ•â•â•â•â•â•â• │ @@@@@@ @@@@@@ "<<endl; - cout<<"│ Version 4.1 │ @@@@@@ @@@@@@@ "<<endl; + cout<<"│ Version dev │ @@@@@@ @@@@@@@ "<<endl; cout<<"├───────────────────────────────────────────────┤ @@@@@@@@ @@@@@@@@ "<<endl; cout<<"│ Core Lybrary for Advance Scenario Simulation │ @@@@@@@@ @ @@@@@@@ "<<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; -- GitLab