|
#include "FabricationPlant.hxx"
|
|
|
|
#include "Storage.hxx"
|
|
#include "FabricationPlant.hxx"
|
|
#include "Reactor.hxx"
|
|
#include "EvolutionData.hxx"
|
|
#include "DataBank.hxx"
|
|
#include "IsotopicVector.hxx"
|
|
#include "CLASS.hxx"
|
|
#include "LogFile.hxx"
|
|
|
|
|
|
|
|
|
|
#include "TMatrixT.h"
|
|
|
|
#include <sstream>
|
|
#include <string>
|
|
#include <iostream>
|
|
#include <cmath>
|
|
#include <algorithm>
|
|
|
|
//________________________________________________________________________
|
|
//
|
|
// FabricationPlant
|
|
//
|
|
//
|
|
//
|
|
//
|
|
//________________________________________________________________________
|
|
ClassImp(FabricationPlant)
|
|
|
|
template <class T> T random(T a, T b) //peak random numebr between a and b
|
|
{
|
|
double range = pow(2., 31);
|
|
srand(time(NULL)); //initialize the srand
|
|
return (T)a + (T)(b-a)*rand()/range;
|
|
}
|
|
|
|
|
|
FabricationPlant::FabricationPlant()
|
|
{
|
|
|
|
|
|
|
|
}
|
|
|
|
FabricationPlant::FabricationPlant(LogFile* log)
|
|
{
|
|
|
|
fLog = log;
|
|
fChronologicalTimePriority = false;
|
|
fFabricationTime = -1;
|
|
fUpdateReferenceDBatEachStep = false;
|
|
fSubstitutionFuel = false;
|
|
|
|
// Warning
|
|
|
|
cout << "!!INFO!! !!!FabricationPlant!!! A FabricationPlant has been define :" << endl;
|
|
cout << "\t Chronological Stock Priority set! "<< endl << endl;
|
|
cout << "!!WARNING!! !!!FabricationPlant!!! You need to set the different stock manually as well as the Fabrication Time Manualy !! " << endl;
|
|
|
|
fLog->fLog << "!!INFO!! !!!FabricationPlant!!! A FabricationPlant has been define :" << endl;
|
|
fLog->fLog << "\t Chronological Stock Priority set! "<< endl << endl;
|
|
fLog->fLog << "!!WARNING!! !!!FabricationPlant!!! You need to set the different stock manually as well as the Fabrication Time Manualy !! " << endl;
|
|
|
|
|
|
}
|
|
|
|
FabricationPlant::FabricationPlant(LogFile* log, Storage* storage, Storage* reusable, double fabircationtime)
|
|
{
|
|
|
|
fLog = log;
|
|
|
|
fChronologicalTimePriority = false;
|
|
fUpdateReferenceDBatEachStep = false;
|
|
fSubstitutionFuel = false;
|
|
|
|
fFabricationTime = (cSecond)fabircationtime;
|
|
fStorage = storage;
|
|
fReUsable = reusable;
|
|
|
|
|
|
cout << "!!INFO!! !!!FabricationPlant!!! A FabricationPlant has been define :" << endl;
|
|
cout << "\t Chronological Stock Priority has been set! "<< endl;
|
|
cout << "\t Fabrication time set to \t " << (double)(fFabricationTime/3600/24/365.25) << " year" << endl << endl;
|
|
|
|
fLog->fLog << "!!INFO!! !!!FabricationPlant!!! A FabricationPlant has been define :" << endl;
|
|
fLog->fLog << "\t Chronological Stock Priority has been set! "<< endl;
|
|
fLog->fLog << "\t Fabrication time set to \t " << (double)(fFabricationTime/3600/24/365.25) << " year" << endl << endl;
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
//________________________________________________________________________
|
|
FabricationPlant::~FabricationPlant()
|
|
{
|
|
|
|
|
|
}
|
|
|
|
|
|
//________________________________________________________________________
|
|
void FabricationPlant::AddValorisableIV(ZAI zai, double factor)
|
|
{
|
|
|
|
pair<map<ZAI, double>::iterator, bool> IResult;
|
|
if(factor > 1) factor = 1;
|
|
|
|
if(factor > 0)
|
|
{
|
|
IResult = fValorisableIV.insert( pair<ZAI ,double>(zai, factor));
|
|
if(IResult.second == false)
|
|
IResult.first->second = factor;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
//________________________________________________________________________
|
|
//_______________________________ Evolution ______________________________
|
|
//________________________________________________________________________
|
|
void FabricationPlant::Evolution(cSecond t)
|
|
{
|
|
|
|
// Check if the FabricationPlant has been created ...
|
|
if(t == fInternalTime && t != 0) return;
|
|
// Make the evolution for the FabricationPlant ...
|
|
FabricationPlantEvolution(t);
|
|
// ... And Finaly update the AbsoluteInternalTime
|
|
fInternalTime = t;
|
|
|
|
}
|
|
|
|
//________________________________________________________________________
|
|
void FabricationPlant::FabricationPlantEvolution(cSecond t)
|
|
{
|
|
|
|
|
|
|
|
map<int ,cSecond >::iterator it;
|
|
for( it = fReactorNextStep.begin(); it!= fReactorNextStep.end(); it++ )
|
|
{
|
|
if( t + fFabricationTime >= fParc->GetReactor()[ (*it).first ]->GetCreationTime()
|
|
&& t + fFabricationTime < fParc->GetReactor()[ (*it).first ]->GetCreationTime() + fParc->GetReactor()[ (*it).first ]->GetLifeTime())
|
|
{
|
|
if( (*it).second == t )
|
|
{
|
|
BuildFuelForReactor( (*it).first );
|
|
(*it).second += fParc->GetReactor()[ (*it).first ]->GetCycleTime();
|
|
}
|
|
else if ( (*it).second - fParc->GetReactor()[ (*it).first ]->GetCycleTime() + fFabricationTime > t )
|
|
{
|
|
map<int ,IsotopicVector >::iterator it2 = fReactorFuturIV.find( (*it).first );
|
|
if (it2 != fReactorFuturIV.end())
|
|
(*it2).second = GetDecay((*it2).second, t - fInternalTime );
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
//________________________________________________________________________
|
|
void FabricationPlant::BuildFuelForReactor(int ReactorId)
|
|
{
|
|
|
|
|
|
DataBank<IsotopicVector>* FuelType = fParc->GetReactor()[ReactorId]->GetFuelType();
|
|
|
|
|
|
if(FuelType->GetFuelType() != "ThPu")
|
|
{
|
|
cout << "!!Bad Trouble!! !!!FabricationPlant!!! Try to do ThPu with a not ThPued DB "<< endl;
|
|
fLog->fLog << "!!Bad Trouble!! !!!FabricationPlant!!! Try to do ThPu with a not ThPued DB" << endl;
|
|
exit (1);
|
|
}
|
|
map<ZAI, double> ZAImass;
|
|
////!!!!
|
|
ZAImass.insert( pair< ZAI,double >( ZAI(90,232,0), 232038060.000e-6 ) );
|
|
|
|
ZAIMass.insert( pair< ZAI,double >( ZAI(92,232,0), 232000000e-6 ) );
|
|
ZAIMass.insert( pair< ZAI,double >( ZAI(92,233,0), 233000000e-6 ) );
|
|
ZAIMass.insert( pair< ZAI,double >( ZAI(92,234,0), 234000000e-6 ) );
|
|
ZAIMass.insert( pair< ZAI,double >( ZAI(92,235,0), 235043929.918e-6 ) );
|
|
ZAIMass.insert( pair< ZAI,double >( ZAI(92,236,0), 236000000e-6 ) );
|
|
ZAIMass.insert( pair< ZAI,double >( ZAI(92,238,0), 238050788.247e-6 ) );
|
|
|
|
ZAImass.insert( pair< ZAI,double >( ZAI(94,238,0), 238049559.894e-6 ) );
|
|
ZAImass.insert( pair< ZAI,double >( ZAI(94,239,0), 239052163.381e-6 ) );
|
|
ZAImass.insert( pair< ZAI,double >( ZAI(94,240,0), 240053813.545e-6 ) );
|
|
ZAImass.insert( pair< ZAI,double >( ZAI(94,241,0), 241056851.456e-6 ) );
|
|
ZAImass.insert( pair< ZAI,double >( ZAI(94,242,0), 242058742.611e-6 ) );
|
|
|
|
ZAImass.insert( pair< ZAI,double >( ZAI(95,241,0), 241056829.144e-6 ) );
|
|
|
|
|
|
double Na = 6.02214129e23; //N Avogadro
|
|
|
|
double HMmass = fParc->GetReactor()[ReactorId]->GetHeavyMetalMass();
|
|
double BU = fParc->GetReactor()[ReactorId]->GetBurnUp();
|
|
IsotopicVector FullUsedStock;
|
|
IsotopicVector stock;
|
|
|
|
bool FuelBuild = false;
|
|
while(!FuelBuild)
|
|
{
|
|
double nPu_0 = 0;
|
|
double MPu_0 = 0;
|
|
{
|
|
map<ZAI ,double>::iterator it;
|
|
|
|
map<ZAI ,double> isotopicquantity = GetDecay( FullUsedStock , fFabricationTime).GetSpeciesComposition(94).GetIsotopicQuantity();
|
|
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
|
|
nPu_0 += (*it).second;
|
|
|
|
isotopicquantity = FullUsedStock.GetSpeciesComposition(94).GetIsotopicQuantity();
|
|
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
|
|
MPu_0 += (*it).second*ZAImass.find( (*it).first )->second/Na*1e-6;
|
|
}
|
|
|
|
stock = GetStockToRecycle();
|
|
|
|
if( stock.GetZAIIsotopicQuantity(ZAI(-1,-1,-1)) == 1 ) // Not enought stock to build the needed fuel
|
|
{
|
|
if (!fSubstitutionFuel)
|
|
{
|
|
{
|
|
EvolutionData evolutiondb;
|
|
pair<map<int, EvolutionData>::iterator, bool> IResult;
|
|
IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,evolutiondb) );
|
|
if(IResult.second == false)
|
|
IResult.first->second = evolutiondb;
|
|
}
|
|
{
|
|
IsotopicVector EmptyIV;
|
|
pair<map<int, IsotopicVector>::iterator, bool> IResult;
|
|
IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId,EmptyIV) );
|
|
if(IResult.second == false)
|
|
IResult.first->second = EmptyIV;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
{
|
|
EvolutionData evolutiondb = fSubstitutionEvolutionData* HMmass;
|
|
pair<map<int, EvolutionData>::iterator, bool> IResult;
|
|
IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,evolutiondb) );
|
|
if(IResult.second == false)
|
|
IResult.first->second = evolutiondb;
|
|
}
|
|
{
|
|
IsotopicVector IV = fSubstitutionEvolutionData.GetIsotopicVectorAt(0)* HMmass;
|
|
pair<map<int, IsotopicVector>::iterator, bool> IResult;
|
|
IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId, IV) );
|
|
if(IResult.second == false)
|
|
IResult.first->second = IV;
|
|
}
|
|
|
|
}
|
|
FuelBuild = true;
|
|
fFractionToTake.clear();
|
|
}
|
|
else
|
|
{
|
|
double nPu_1 = 0;
|
|
double MPu_1 = 0;
|
|
double Sum_AlphaI_nPuI = 0;
|
|
double Sum_AlphaI_nPuI0 = 0;
|
|
{
|
|
map<ZAI ,double>::iterator it;
|
|
map<ZAI ,double> isotopicquantity = GetDecay( stock , fFabricationTime).GetSpeciesComposition(94).GetIsotopicQuantity();
|
|
|
|
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
|
|
{
|
|
if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
|
|
{
|
|
nPu_1 += (*it).second;
|
|
Sum_AlphaI_nPuI += FuelType->GetFuelParameter()[(*it).first.A() -237]*(*it).second;
|
|
}
|
|
}
|
|
|
|
isotopicquantity = stock.GetSpeciesComposition(94).GetIsotopicQuantity();
|
|
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
|
|
if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
|
|
{
|
|
MPu_1 += (*it).second * (ZAImass.find( (*it).first )->second)/Na*1e-6;
|
|
}
|
|
|
|
isotopicquantity = GetDecay( FullUsedStock , fFabricationTime).GetSpeciesComposition(94).GetIsotopicQuantity();
|
|
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
|
|
if ((*it).first.A() >= 238 && (*it).first.A() <= 242)
|
|
{
|
|
Sum_AlphaI_nPuI0 += FuelType->GetFuelParameter()[(*it).first.A() -237]*(*it).second;
|
|
}
|
|
}
|
|
|
|
double StockFactionToUse = 0;
|
|
////!!!
|
|
double NT = HMmass*1e6 * Na / (ZAImass.find( ZAI(90,232,0) )->second); //I changed the the depleted uranium into Th
|
|
|
|
double N1 = (BU - FuelType->GetFuelParameter()[6]) * NT;
|
|
double N2 = -Sum_AlphaI_nPuI0;
|
|
////!!!
|
|
double N3 = -FuelType->GetFuelParameter()[0] * Na / (ZAImass.find( ZAI(90,232,0) )->second) * (HMmass*1e6 - MPu_0*1e6); //The same as above
|
|
|
|
double D1 = Sum_AlphaI_nPuI;
|
|
////!!!
|
|
double D2 = -FuelType->GetFuelParameter()[0] * MPu_1*1e6 * Na / (ZAImass.find( ZAI(90,232,0) )->second ) ; //The same as above
|
|
|
|
StockFactionToUse = (N1 + N2 + N3) / (D1 + D2);
|
|
|
|
if(StockFactionToUse < 0)
|
|
{
|
|
cout << "!!Bad Trouble!! !!!FabricationPlant!!! Oups Bug in calculating stock fraction to use "<< endl;
|
|
fLog->fLog << "!!Bad Trouble!! !!!FabricationPlant!!! Oups Bug in calculating stock fraction to use" << endl;
|
|
exit (1);
|
|
}
|
|
|
|
if( StockFactionToUse > 1 )
|
|
{
|
|
|
|
FullUsedStock += stock;
|
|
RecycleStock(1);
|
|
FuelBuild = false;
|
|
}
|
|
else
|
|
{
|
|
RecycleStock(StockFactionToUse);
|
|
|
|
IsotopicVector IVBeginCycle;
|
|
FuelBuild = true;
|
|
////!!!!
|
|
ZAI Th = ZAI(90,232,0); //I changed the U8 = ZAI(92,238,0) here
|
|
// ZAI U5 = ZAI(92,235,0);
|
|
////!!!!
|
|
double Th_Quantity = (HMmass - (MPu_0+StockFactionToUse*MPu_1 ))/(ZAImass.find( ZAI(90,232,0) )->second ) *Na/1e-6; //I changed the the depleted uranium into Th
|
|
////!!!
|
|
fParc->AddGodIncome( Th, Th_Quantity );
|
|
//fParc->AddGodIncome( U5, U8_Quantity*0.003 );
|
|
|
|
for(int i = (int)fFractionToTake.size()-1; i >= 0; i--)
|
|
{
|
|
IVBeginCycle += fStorage->GetStock()[fFractionToTake[i].first].GetSpeciesComposition(94)*( fFractionToTake[i].second );
|
|
fReUsable->AddToStock(fStorage->GetStock()[fFractionToTake[i].first]*(fFractionToTake[i].second)
|
|
- fStorage->GetStock()[fFractionToTake[i].first].GetSpeciesComposition(94)*(fFractionToTake[i].second));
|
|
|
|
fStorage->TakeFractionFromStock(fFractionToTake[i].first,fFractionToTake[i].second);
|
|
|
|
}
|
|
fFractionToTake.clear();
|
|
|
|
////!!!!
|
|
IVBeginCycle += Th_Quantity*Th;
|
|
EvolutionData evolutiondb = BuildEvolutiveDB(ReactorId, IVBeginCycle);
|
|
|
|
{
|
|
pair<map<int, EvolutionData>::iterator, bool> IResult;
|
|
IResult = fReactorFuturDB.insert( pair<int, EvolutionData>(ReactorId,evolutiondb) );
|
|
if(IResult.second == false)
|
|
IResult.first->second = evolutiondb;
|
|
}
|
|
{
|
|
pair<map<int, IsotopicVector>::iterator, bool> IResult;
|
|
IResult = fReactorFuturIV.insert( pair<int, IsotopicVector>(ReactorId,IVBeginCycle) );
|
|
if(IResult.second == false)
|
|
IResult.first->second = IVBeginCycle;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
}
|
|
|
|
void FabricationPlant::SetSubstitutionFuel(EvolutionData fuel)
|
|
{
|
|
|
|
fSubstitutionFuel = true;
|
|
fSubstitutionEvolutionData = fuel / fuel.GetHMMass();
|
|
|
|
}
|
|
|
|
|
|
//________________________________________________________________________
|
|
//_________________________________ Decay ________________________________
|
|
//________________________________________________________________________
|
|
IsotopicVector FabricationPlant::GetDecay(IsotopicVector isotopicvector, cSecond t)
|
|
{
|
|
|
|
IsotopicVector IV;
|
|
|
|
map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
|
|
|
|
map<ZAI ,double >::iterator it;
|
|
for( it = isotopicquantity.begin(); it != isotopicquantity.end(); it++ )
|
|
{
|
|
if((*it).second > 0)
|
|
{
|
|
IsotopicVector ivtmp = fDecayDataBase->Evolution(it->first, t) * (*it).second ;
|
|
IV += ivtmp;
|
|
}
|
|
}
|
|
|
|
return IV;
|
|
}
|
|
|
|
|
|
//________________________________________________________________________
|
|
//_____________________________ Reactor & DB _____________________________
|
|
//________________________________________________________________________
|
|
EvolutionData FabricationPlant::BuildEvolutiveDB(int ReactorId,IsotopicVector isotopicvector)
|
|
{
|
|
|
|
DataBank<IsotopicVector>* evolutiondb = fParc->GetReactor()[ReactorId]->GetFuelType();
|
|
|
|
isotopicvector = GetDecay(isotopicvector, fFabricationTime);
|
|
|
|
EvolutionData EvolBuild;
|
|
/*
|
|
if( fUpdateReferenceDBatEachStep == true )
|
|
{
|
|
EvolutionData EvolBuild = evolutiondb->GenerateDB(isotopicvector,
|
|
fParc->GetReactor()[ReactorId]->GetCycleTime(),
|
|
fParc->GetReactor()[ReactorId]->GetPower());
|
|
}
|
|
else
|
|
{
|
|
map<double, EvolutionData> distances = evolutiondb->GetDistancesTo(isotopicvector);
|
|
EvolBuild = distances.begin()->second.GenerateDBFor(isotopicvector);
|
|
|
|
}
|
|
*/
|
|
EvolBuild = evolutiondb->GenerateEvolutionData(isotopicvector,
|
|
fParc->GetReactor()[ReactorId]->GetCycleTime(),
|
|
fParc->GetReactor()[ReactorId]->GetPower());
|
|
|
|
return EvolBuild;
|
|
|
|
}
|
|
|
|
//________________________________________________________________________
|
|
void FabricationPlant::TakeReactorFuel(int Id)
|
|
{
|
|
|
|
|
|
IsotopicVector IV;
|
|
map<int ,IsotopicVector >::iterator it2 = fReactorFuturIV.find( Id );
|
|
if (it2 != fReactorFuturIV.end())
|
|
(*it2).second = IV;
|
|
|
|
}
|
|
|
|
//________________________________________________________________________
|
|
EvolutionData FabricationPlant::GetReactorEvolutionDB(int ReactorId)
|
|
{
|
|
|
|
map< int,EvolutionData >::iterator it = fReactorFuturDB.find(ReactorId);
|
|
return (*it).second;
|
|
|
|
}
|
|
|
|
IsotopicVector FabricationPlant::GetFullFabrication()
|
|
{
|
|
|
|
IsotopicVector tmp = 0*ZAI(0,0,0);
|
|
|
|
map<int, IsotopicVector > reactorNextStep = fReactorFuturIV;
|
|
map<int, IsotopicVector >::iterator it;
|
|
for ( it = reactorNextStep.begin(); it != reactorNextStep.end(); it++)
|
|
tmp += (*it).second;
|
|
|
|
return tmp;
|
|
|
|
}
|
|
|
|
//________________________________________________________________________
|
|
//_______________________________ Storage ________________________________
|
|
//________________________________________________________________________
|
|
|
|
IsotopicVector FabricationPlant::GetStockToRecycle()
|
|
{
|
|
|
|
IsotopicVector NextStock;
|
|
int IdToTake = -1;
|
|
|
|
if(fChronologicalTimePriority == true)
|
|
IdToTake = (int)( fFractionToTake.size() );
|
|
else
|
|
IdToTake = (int)( fStorage->GetStock().size() -1 - fFractionToTake.size() );
|
|
if(0 <= IdToTake && IdToTake <= (int)fStorage->GetStock().size()-1)
|
|
{
|
|
NextStock = fStorage->GetStock()[IdToTake];
|
|
fFractionToTake.push_back( pair<int,double>(IdToTake,0.) );
|
|
}
|
|
else NextStock += ZAI(-1,-1,-1) *1;
|
|
|
|
return NextStock;
|
|
|
|
}
|
|
|
|
//________________________________________________________________________
|
|
void FabricationPlant::RecycleStock(double fraction)
|
|
{
|
|
|
|
fFractionToTake.back().second = fraction;
|
|
|
|
}
|
|
|
|
|
|
|
|
//________________________________________________________________________
|
|
void FabricationPlant::DumpStock()
|
|
{
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
//________________________________________________________________________
|
|
pair<IsotopicVector, IsotopicVector> FabricationPlant::Separation(IsotopicVector isotopicvector)
|
|
{
|
|
|
|
//[0] = re-use ; [1] = waste
|
|
pair<IsotopicVector, IsotopicVector> IVTmp;
|
|
|
|
map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
|
|
map<ZAI ,double >::iterator it;
|
|
for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
|
|
{
|
|
map<ZAI ,double>::iterator it2;
|
|
it2 = fValorisableIV.find((*it).first);
|
|
if ( it2 != fValorisableIV.end() )
|
|
{
|
|
IVTmp.first.Add( (*it).first, (*it).second * (*it2).second ); //re-use
|
|
IVTmp.second.Add( (*it).first, (*it).second * (1-(*it2).second) ); //waste
|
|
}
|
|
else IVTmp.second.Add( (*it).first, (*it).second ); //waste
|
|
}
|
|
|
|
return IVTmp;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|