From edc6f6631919de077e7a2d026d65a4f39dfd9dce Mon Sep 17 00:00:00 2001
From: Baptiste LENIAU <baptiste.leniau@subatech.in2p3.fr>
Date: Wed, 29 Oct 2014 09:35:07 +0000
Subject: [PATCH] =?UTF-8?q?Correct=20negative=20value=20when=20using=20?=
=?UTF-8?q?=C2=AB=C2=A0IV-=3D=C2=AB=C2=A0=20due=20to=20double=20precision.?=
=?UTF-8?q?=20Correct=20a=20bug=20due=20to=20internal=20time=20of=20some?=
=?UTF-8?q?=20factories=20not=20defined=20when=20created.=20This=20was=20c?=
=?UTF-8?q?ausing=20the=20code=20to=20make=20an=20initial=20stock=20to=20d?=
=?UTF-8?q?ecay=20while=20we=20don=E2=80=99t=20want=20it=20to=20do=20so?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
git-svn-id: svn+ssh://svn.in2p3.fr/class@412 0e7d625b-0364-4367-a6be-d5be4a48d228
---
source/branches/CLASSV3/src/CLASSBackEnd.cxx | 10 +++
source/branches/CLASSV3/src/CLASSFuelPlan.cxx | 36 +++++----
.../branches/CLASSV3/src/FabricationPlant.cxx | 74 +++++++++++--------
.../branches/CLASSV3/src/IsotopicVector.cxx | 10 ++-
source/branches/CLASSV3/src/Pool.cxx | 10 +--
source/branches/CLASSV3/src/Scenario.cxx | 8 +-
.../branches/CLASSV3/src/SeparationPlant.cxx | 1 -
source/branches/CLASSV3/src/Storage.cxx | 19 ++---
8 files changed, 94 insertions(+), 74 deletions(-)
diff --git a/source/branches/CLASSV3/src/CLASSBackEnd.cxx b/source/branches/CLASSV3/src/CLASSBackEnd.cxx
index 03870e12e..faf588e85 100644
--- a/source/branches/CLASSV3/src/CLASSBackEnd.cxx
+++ b/source/branches/CLASSV3/src/CLASSBackEnd.cxx
@@ -58,6 +58,16 @@ void CLASSBackEnd::AddIV(IsotopicVector isotopicvector)
}
+//________________________________________________________________________
+void CLASSBackEnd::UpdateInsideIV()
+{
+ DBGL
+ fInsideIV = IsotopicVector();
+ for(int i = 0; i < (int)fIVArray.size(); i++)
+ fInsideIV += fIVArray[i];
+ DBGL
+}
+
//________________________________________________________________________
// Get Decay
diff --git a/source/branches/CLASSV3/src/CLASSFuelPlan.cxx b/source/branches/CLASSV3/src/CLASSFuelPlan.cxx
index 42cc1e0d9..a57a24d78 100644
--- a/source/branches/CLASSV3/src/CLASSFuelPlan.cxx
+++ b/source/branches/CLASSV3/src/CLASSFuelPlan.cxx
@@ -38,31 +38,29 @@ pair< CLASSFuel, double > CLASSFuelPlan::GetFuelAt(cSecond t)
pair< CLASSFuel, double > FuelAtT = fLoadingPlan.begin()->second;
map< cSecond, pair< CLASSFuel, double > >::iterator it;
- bool AfterPrevious = false;
- bool AfterNext = true;
+
+ cSecond ThisFuelTime;
+
for (it = fLoadingPlan.begin(); it != fLoadingPlan.end(); it++ )
{
- if (t < (*it).first )
- AfterNext = false;
- else
- AfterNext = true;
-
- if (AfterPrevious && !AfterNext)
- {
- DBGL
- return FuelAtT;
- }
- else if (!AfterPrevious && !AfterNext)
+ if( it == fLoadingPlan.begin())
{
- WARNING << "The time asked is before the first laoding time..."<< endl;
- WARNING << "The first Fuel will be loaded... Check your FuelPLan!!!!!" << endl;
- DBGL
- return FuelAtT;
+ FuelAtT = (*it).second;
}
else
{
- FuelAtT = (*it).second;
- AfterPrevious = true;
+ ThisFuelTime = (*it).first;
+
+ if( t < ThisFuelTime )
+ {
+ DBGL
+ return FuelAtT;
+ }
+ else
+ {
+ FuelAtT = (*it).second;
+ }
+
}
}
diff --git a/source/branches/CLASSV3/src/FabricationPlant.cxx b/source/branches/CLASSV3/src/FabricationPlant.cxx
index a033d338a..8cc636737 100644
--- a/source/branches/CLASSV3/src/FabricationPlant.cxx
+++ b/source/branches/CLASSV3/src/FabricationPlant.cxx
@@ -100,6 +100,8 @@ void FabricationPlant::Evolution(cSecond t)
if(t == fInternalTime && t != 0) return;
// Make the evolution for the FabricationPlant ...
FabricationPlantEvolution(t);
+ //Update Inside IsotopicVector
+ UpdateInsideIV();
// ... And Finaly update the AbsoluteInternalTime
fInternalTime = t;
@@ -109,41 +111,44 @@ void FabricationPlant::Evolution(cSecond t)
void FabricationPlant::FabricationPlantEvolution(cSecond t)
{
DBGL
- IsotopicVector fInsideIV;
-
-
map<int ,cSecond >::iterator 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();
+
+ 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 / R_Power * R_HMMass * 1e9 *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;
+ exit(1);
+ }
+
if( t + GetCycleTime() >= R_CreactionTime
&& t + GetCycleTime() < R_CreactionTime + R_LifeTime)
{
if( (*it).second == t )
{
- int ReactorId = (*it).first;
- pair<CLASSFuel, double> R_Fuel = GetParc()->GetReactor()[ReactorId]->GetFuelPlan()->GetFuelAt( t + GetCycleTime() );
#pragma omp critical(FuelBuild)
{
if( R_Fuel.first.GetPhysicsModels() )
{
BuildFuelForReactor( (*it).first, t );
}
+ (*it).second += R_CycleTime;
}
- double R_BU = R_Fuel.second;
- double R_Power = GetParc()->GetReactor()[ReactorId]->GetPower();
- double R_HMMass = GetParc()->GetReactor()[ReactorId]->GetHeavyMetalMass();
- (*it).second += (cSecond) (R_BU / R_Power * R_HMMass * 1e9 *3600*24);
}
- else if ( (*it).second - GetParc()->GetReactor()[ (*it).first ]->GetCycleTime() + GetCycleTime() > t )
+ else if ( (*it).second - R_CycleTime + GetCycleTime() >= t && (*it).second - R_CycleTime < t )
{
map<int ,IsotopicVector >::iterator it2 = fReactorFuturIV.find( (*it).first );
if (it2 != fReactorFuturIV.end())
- (*it2).second = GetDecay((*it2).second, t - fInternalTime );
-
- fInsideIV += (*it2).second;
+ (*it2).second = GetDecay((*it2).second, t - fInternalTime );
}
}
}
@@ -152,10 +157,23 @@ DBGL
DBGL
}
+void FabricationPlant::UpdateInsideIV()
+{
+ DBGL
+ fInsideIV = IsotopicVector();
+
+ map< int,IsotopicVector >::iterator it;
+ for( it = fReactorFuturIV.begin(); it != fReactorFuturIV.end(); it++ )
+ fInsideIV += (*it).second;
+
+ 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;
@@ -213,7 +231,7 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
}
fInsideIV += IV;
AddCumulativeIVIn(IV);
-
+ DBGL
return;
}
else
@@ -261,7 +279,7 @@ void FabricationPlant::BuildFuelForReactor(int ReactorId, cSecond t)
- }
+ }DBGL
return;
}
DBGL
@@ -440,21 +458,20 @@ void FabricationPlant::SetSubstitutionFuel(EvolutionData fuel)
//________________________________________________________________________
void FabricationPlant::TakeReactorFuel(int Id)
{
-
-
+DBGL
IsotopicVector IV;
map<int ,IsotopicVector >::iterator it2 = fReactorFuturIV.find( Id );
- AddCumulativeIVOut(it2->second);
- fInsideIV -= (*it2).second;
+ AddCumulativeIVOut(it2->second);
if (it2 != fReactorFuturIV.end())
(*it2).second = IV;
-
map< int,EvolutionData >::iterator it = fReactorFuturDB.find(Id);
(*it).second = EvolutionData();
+ UpdateInsideIV();
+DBGL
}
//________________________________________________________________________
@@ -481,9 +498,10 @@ DBGL
int IV_N = fFissileArrayAdress[i].second;
pair<IsotopicVector, IsotopicVector> Separated_Lost;
- Separated_Lost = Separation( fFissileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFertileList);
+ Separated_Lost = Separation( fFissileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFissileList );
BuildedFuel += Separated_Lost.first;
Lost += Separated_Lost.second;
+
}
}
@@ -497,7 +515,7 @@ DBGL
int IV_N = fFertileArrayAdress[i].second;
pair<IsotopicVector, IsotopicVector> Separated_Lost;
- Separated_Lost = Separation( fFertileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFissileList);
+ Separated_Lost = Separation( fFertileStorage[Stor_N]->GetIVArray()[IV_N]*LambdaArray[i], fFertileList);
BuildedFuel += Separated_Lost.first;
Lost += Separated_Lost.second;
}
@@ -565,14 +583,12 @@ DBGL
//________________________________________________________________________
pair<IsotopicVector, IsotopicVector> FabricationPlant::Separation(IsotopicVector isotopicvector, IsotopicVector ExtractedList)
{
-
+DBGL
//[0] = re-use ; [1] = waste
- IsotopicVector LostPart = isotopicvector.GetThisComposition(ExtractedList) * fSeparationLostFraction;
-
- IsotopicVector SeparatedPart = isotopicvector.GetThisComposition(ExtractedList) - LostPart;
- 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);
}
diff --git a/source/branches/CLASSV3/src/IsotopicVector.cxx b/source/branches/CLASSV3/src/IsotopicVector.cxx
index 177562c47..982fa5056 100755
--- a/source/branches/CLASSV3/src/IsotopicVector.cxx
+++ b/source/branches/CLASSV3/src/IsotopicVector.cxx
@@ -392,10 +392,11 @@ void IsotopicVector::Remove(const ZAI& zai, double quantity)
if ( it != fIsotopicQuantity.end() )
{
if (it->second > quantity)
- it->second = it->second - quantity;
+ it->second = it->second - quantity;
else
{
- Need(zai, quantity - it->second );
+ if( (it->second - quantity)/it->second > 1e-16 ) // to fit with double precision : 16 digits
+ Need(zai, quantity - it->second );
it->second = 0;
}
}
@@ -422,11 +423,12 @@ void IsotopicVector::Remove(const IsotopicVector& isotopicvector)
//________________________________________________________________________
void IsotopicVector::Need(const ZAI& zai, double quantity)
{
- cout << "Negative quantity for ZAI " << zai.Z() << " " << zai.A() << " " << zai.I() << " in this IsotopicVector" << endl;
+ if(quantity < 0.5) quantity = 0;
pair<map<ZAI, double>::iterator, bool> IResult;
if(quantity > 0)
- {
+ { cout << "Negative quantity : "<<quantity <<" for ZAI " << zai.Z() << " " << zai.A() << " " << zai.I() << " in this IsotopicVector" << endl;
+ exit(0);
IResult = fIsotopicQuantityNeeded.insert( pair<ZAI ,double>(zai, quantity));
if(!IResult.second)
IResult.first->second += quantity;
diff --git a/source/branches/CLASSV3/src/Pool.cxx b/source/branches/CLASSV3/src/Pool.cxx
index 60f97d763..cbbb3cc0b 100755
--- a/source/branches/CLASSV3/src/Pool.cxx
+++ b/source/branches/CLASSV3/src/Pool.cxx
@@ -43,7 +43,6 @@ Pool::Pool(CLASSLogger* log, cSecond coolingtime):CLASSBackEnd(log, coolingtime,
INFO << " A new Pool has been define :" << endl;
- INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
INFO << "\t The Cooling Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
WARNING << " All Cooled Fuel goes directly to WASTE after cooling !! " << endl;
@@ -66,7 +65,6 @@ Pool::Pool(CLASSLogger* log, CLASSBackEnd* storage, cSecond coolingtime):CLASSBa
INFO << " A new Pool has been define :" << endl;
- INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
INFO << "\t The Cooling Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
DBGL
@@ -123,7 +121,7 @@ void Pool::RemoveIVCooling(int i) //!< Remove a Cooling IsotopicVector
fIVArray.erase(fIVArray.begin()+i);
fIVArrayArrivalTime.erase( fIVArrayArrivalTime.begin()+i);
fCoolingIndex.erase(fCoolingIndex.begin()+i);
-
+ UpdateInsideIV();
}
@@ -148,8 +146,6 @@ DBGL
int RemainingCoolingTime;
cSecond EvolutionTime = t - fInternalTime;
- fInsideIV = IsotopicVector();
-
#pragma omp parallel for
for ( int i = 0 ; i < (int)fIVArray.size() ; i++)
{
@@ -173,7 +169,6 @@ DBGL
else if ( fIVArrayArrivalTime[i] != t )
{
fIVArray[i] = GetDecay( fIVArray[i] , EvolutionTime);
- fInsideIV += fIVArray[i];
}
}
#pragma omp critical(DeleteCoolingIVPB)
@@ -191,7 +186,8 @@ void Pool::Evolution(cSecond t)
if(t == fInternalTime && t!=0) return;
// Make the evolution for the Cooling IV ...
CoolingEvolution(t);
-
+ // Update Inside IV
+ UpdateInsideIV();
// ... And Finaly update the AbsoluteInternalTime
fInternalTime = t;
diff --git a/source/branches/CLASSV3/src/Scenario.cxx b/source/branches/CLASSV3/src/Scenario.cxx
index 36d80e9eb..0feee545d 100755
--- a/source/branches/CLASSV3/src/Scenario.cxx
+++ b/source/branches/CLASSV3/src/Scenario.cxx
@@ -184,7 +184,7 @@ void Scenario::AddPool(Pool* Pool)
fPool.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
fPool.back()->SetLog(GetLog());
fPool.back()->SetId((int)fPool.size()-1);
- fPool.back()->SetCreationTime(fAbsoluteTime);
+ fPool.back()->SetInternalTime(fAbsoluteTime);
string Pool_name = fPool.back()->GetName();
@@ -218,6 +218,7 @@ void Scenario::AddReactor(Reactor* reactor)
fReactor.back()->SetParc(this);
fReactor.back()->SetLog(GetLog());
fReactor.back()->SetId((int)fReactor.size()-1);
+ fReactor.back()->SetInternalTime(fAbsoluteTime);
string Reactor_name = fReactor.back()->GetName();
@@ -252,7 +253,7 @@ void Scenario::AddStorage(Storage* storage)
fStorage.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
fStorage.back()->SetLog(GetLog());
fStorage.back()->SetId((int)fStorage.size()-1);
- fStorage.back()->SetCreationTime(fAbsoluteTime);
+ fStorage.back()->SetInternalTime(fAbsoluteTime);
string Storage_name = fStorage.back()->GetName();
@@ -286,6 +287,7 @@ void Scenario::AddFabricationPlant(FabricationPlant* fabricationplant)
fFabricationPlant.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
fFabricationPlant.back()->SetLog(GetLog());
fFabricationPlant.back()->SetId((int)fStorage.size()-1);
+ fFabricationPlant.back()->SetInternalTime(fAbsoluteTime);
string FP_name = fFabricationPlant.back()->GetName();
@@ -318,7 +320,7 @@ void Scenario::AddSeparationPlant(SeparationPlant* SeparationPlant)
fSeparationPlant.back()->SetDecayDataBank( (*this).GetDecayDataBase() );
fSeparationPlant.back()->SetLog(GetLog());
fSeparationPlant.back()->SetId((int)fSeparationPlant.size()-1);
- fSeparationPlant.back()->SetCreationTime(fAbsoluteTime);
+ fSeparationPlant.back()->SetInternalTime(fAbsoluteTime);
string SeparationPlant_name = fSeparationPlant.back()->GetName();
diff --git a/source/branches/CLASSV3/src/SeparationPlant.cxx b/source/branches/CLASSV3/src/SeparationPlant.cxx
index cf2568af6..065f66a41 100644
--- a/source/branches/CLASSV3/src/SeparationPlant.cxx
+++ b/source/branches/CLASSV3/src/SeparationPlant.cxx
@@ -45,7 +45,6 @@ SeparationPlant::SeparationPlant(CLASSLogger* log):CLASSBackEnd(log, -2)
SetIsStorageType();
INFO << " A new SeparationPlant has been define :" << endl;
- INFO << "\t Creation time set at \t " << (double)(GetCreationTime()/3600/24/365.25) << " year" << endl;
INFO << "\t The Separation Time set at\t " << (double)(fCycleTime/3600/24/365.25) << " year" << endl;
WARNING << " All Separated Fuel go directly to WASTE after cooling !! " << endl;
diff --git a/source/branches/CLASSV3/src/Storage.cxx b/source/branches/CLASSV3/src/Storage.cxx
index 43baac497..5bbff10fe 100644
--- a/source/branches/CLASSV3/src/Storage.cxx
+++ b/source/branches/CLASSV3/src/Storage.cxx
@@ -91,7 +91,7 @@ void Storage::AddIV(IsotopicVector isotopicvector)
//________________________________________________________________________
void Storage::TakeFractionFromStock(int IVId,double fraction)
{
-
+DBGL
if(GetParc())
{
if(GetParc()->GetStockManagement() )
@@ -103,8 +103,6 @@ void Storage::TakeFractionFromStock(int IVId,double fraction)
else
{
AddCumulativeIVOut(fIVArray[IVId]*fraction);
-
- fInsideIV -= fIVArray[IVId] * fraction;
fIVArray[IVId] -= fIVArray[IVId] * fraction;
}
@@ -125,15 +123,12 @@ void Storage::TakeFractionFromStock(int IVId,double fraction)
else
{
AddCumulativeIVOut(fIVArray[IVId]*fraction);
-
- fInsideIV -= fIVArray[IVId] * fraction;
fIVArray[IVId] -= fIVArray[IVId] * fraction;
}
}
-
-
-
+ UpdateInsideIV();
+ DBGL
}
void Storage::TakeFromStock(IsotopicVector isotopicvector)
@@ -171,7 +166,7 @@ DBGL
RemoveEmptyStocks();
- int EvolutionTime = t - fInternalTime;
+ cSecond EvolutionTime = t - fInternalTime;
fInsideIV = GetDecay(fInsideIV , EvolutionTime);
@@ -191,11 +186,13 @@ void Storage::Evolution(cSecond t)
if(t == fInternalTime && t!=0) return;
// Make the evolution for the Storage ...
StorageEvolution(t);
+
+ // Update Inside IV;
+ UpdateInsideIV();
+
// ... And Finaly update the AbsoluteInternalTime
fInternalTime = t;
-
-
}
void Storage::Write(string filename, cSecond date)
--
GitLab