From fc811030085f58d6db775cf4b6332c7d80388117 Mon Sep 17 00:00:00 2001
From: Baptiste Mouginot <mouginot.baptiste@gmail.com>
Date: Mon, 16 Feb 2015 10:47:21 +0000
Subject: [PATCH] main algorythm for IM decay Matrix
git-svn-id: svn+ssh://svn.in2p3.fr/class@604 0e7d625b-0364-4367-a6be-d5be4a48d228
---
.../Model/Irradiation/IM_Matrix_Decay.cxx | 272 ++++--------------
.../Model/Irradiation/IM_Matrix_Decay.hxx | 15 +-
2 files changed, 64 insertions(+), 223 deletions(-)
diff --git a/source/trunk/Model/Irradiation/IM_Matrix_Decay.cxx b/source/trunk/Model/Irradiation/IM_Matrix_Decay.cxx
index 27e3f2ca6..a74c953ae 100644
--- a/source/trunk/Model/Irradiation/IM_Matrix_Decay.cxx
+++ b/source/trunk/Model/Irradiation/IM_Matrix_Decay.cxx
@@ -32,254 +32,101 @@ using namespace std;
//________________________________________________________________________
IM_Matrix_Decay::IM_Matrix_Decay():IrradiationModel(new CLASSLogger("IM_Matrix_Decay.log"))
{
- fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
+ fShorstestHalflife = 0;
+ fZAIThreshold = 0;
+ NuclearDataInitialization();
+ fExponentialDecayMatrix = ExponentialCalculation(fDecayMatrix);
+
}
IM_Matrix_Decay::IM_Matrix_Decay(CLASSLogger* log):IrradiationModel(log)
{
- fShorstestHalflife = 3600.*24*160.; //cut by default all nuclei with a shorter liftime than the Cm242 -> remain 33 actinides
-}
-
+ fShorstestHalflife = 0;
+ fZAIThreshold = 0;
+ NuclearDataInitialization();
+ fExponentialDecayMatrix = ExponentialCalculation(fDecayMatrix);
+}
-//________________________________________________________________________
-/* Evolution Calculation */
-//________________________________________________________________________
-EvolutionData IM_Matrix_Decay::GenerateEvolutionData(IsotopicVector isotopicvector, EvolutionData XSSet, double Power, double cycletime)
+TMatrixT<double> IM_Matrix_Decay::ExponentialCalculation(TMatrixT<double> myMatrix)
{
- DBGL
- if(fFastDecay.size() == 0)
- NuclearDataInitialization();
-
+ TMatrixT<double> IdMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());
+ for(int j = 0; j < (int)fReverseMatrixIndex.size(); j++)
+ for(int k = 0; k < (int)fReverseMatrixIndex.size(); k++)
+ {
+ if(k == j) IdMatrix[j][k] = 1;
+ else IdMatrix[j][k] = 0;
+ }
- string ReactorType;
+ TMatrixT<double> MatrixDL = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size()); // Order 0 Term from the DL : Id
+ TMatrixT<double> MatrixDLTermN = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size()); // Addind it;
- vector< TMatrixT<double> > NMatrix ;// TMatrixT<double>(decayindex.size(),1))
- { // Filling the t=0 State;
- map<ZAI, double > isotopicquantity = isotopicvector.GetIsotopicQuantity();
- TMatrixT<double> N_0Matrix = TMatrixT<double>( fReverseMatrixIndex.size(),1) ;
- for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
- N_0Matrix[i] = 0;
-
- map<ZAI, double >::iterator it ;
- for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
- N_0Matrix[i] = 0;
+ {
+ MatrixDLTermN = IdMatrix;
+ MatrixDL = MatrixDLTermN;
+ int j = 1;
+ double NormN;
- for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
+ do
{
- /// Need TO change with FP managment
- map<ZAI, int >::iterator it2;
+ TMatrixT<double> MatrixDLTermtmp = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size()); // Adding it;
+ MatrixDLTermtmp = MatrixDLTermN;
- if( (*it).first.Z() < fZAIThreshold )
- it2 = fMatrixIndex.find( ZAI(-2,-2,-2) );
- else it2 = fMatrixIndex.find( (*it).first );
+ MatrixDLTermN.Mult(MatrixDLTermtmp, myMatrix );
- if(it2 == fMatrixIndex.end() ) //If not in index should be TMP, can't be fast decay for new Fuel !!!
- it2 = fMatrixIndex.find( ZAI(-3,-3,-3) );
- N_0Matrix[ (*it2).second ][0] = (*it).second ;
+ MatrixDLTermN *= 1./j;
+ MatrixDL += MatrixDLTermN;
+ NormN = 0;
+ for(int m = 0; m < (int)fReverseMatrixIndex.size(); m++)
+ for(int n = 0; n < (int)fReverseMatrixIndex.size(); n++)
+ NormN += MatrixDLTermN[m][n]*MatrixDLTermN[m][n];
+ j++;
- }
-
- isotopicquantity.clear();
-
- NMatrix.push_back(N_0Matrix);
- N_0Matrix.Clear();
-
+ } while ( NormN != 0 );
}
+ return MatrixDL;
+}
+
+
+IsotopicVector IM_Matrix_Decay::GetDecay(IsotopicVector Mother_IV, time)
+{
- //-------------------------//
- //--- Perform Evolution ---//
- //-------------------------//
- ReactorType = XSSet.GetReactorType();
-
- double M_ref = XSSet.GetHeavyMetalMass();
- double M = cZAIMass.GetMass(isotopicvector.GetActinidesComposition());
- double Power_ref = XSSet.GetPower();
-
- int NStep = XSSet.GetFissionXS().begin()->second->GetN();
- double* DBTimeStep = XSSet.GetFissionXS().begin()->second->GetX();
-
- int InsideStep = 10;
-
- double timevector[NStep];
- timevector[0] = 0;
-
- double Flux[NStep];
- TMatrixT<double> FissionEnergy = TMatrixT<double>(fReverseMatrixIndex.size(),1);
+ TMatrixT<double> NMatrix = TMatrixT<double>(decayindex.size(),1))
for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
- FissionEnergy[i] = 0;
-
- {
- map< ZAI, int >::iterator it;
- for(it = fMatrixIndex.begin(); it != fMatrixIndex.end(); it++)
- {
- map< ZAI, double >::iterator it2 = fFissionEnergy.find(it->first);
- if(it2 == fFissionEnergy.end())
- {
- if(it->first.Z() > fZAIThreshold)
- FissionEnergy[it->second][0] = 1.9679e6*it->first.A()-2.601e8; // //simple linear fit to known values ;extrapolation to unknown isotopes
- else FissionEnergy[it->second][0] = 0;
- }
- else
- FissionEnergy[it->second][0] = it2->second;
-
- }
- }
-
- 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
- DBGL
- for(int i = 0; i < NStep-1; i++)
- {
- double TStepMax = ( (DBTimeStep[i+1]-DBTimeStep[i] ) ) * Power_ref/M_ref / Power*M ; // Get the next Time step
-
-
- TMatrixT<double> BatemanMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());
- TMatrixT<double> BatemanReactionMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());
-
- TMatrixT<double> NEvolutionMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),1);
- NEvolutionMatrix = NMatrix.back();
-
-
-
- FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[i])); //Feel the fission reaction Matrix
- CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[i])); //Feel the capture reaction Matrix
- n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[i])); //Feel the (n,2n) reaction Matrix
-
- // ---------------- Evolution
+ NMatrix[i] = 0;
+ { // Filling the t=0 State;
+ map<ZAI, double > isotopicquantity = Mother_IV.GetIsotopicQuantity();
- BatemanReactionMatrix = FissionXSMatrix[i];
- BatemanReactionMatrix += CaptureXSMatrix[i];
- BatemanReactionMatrix += n2nXSMatrix[i];
+ map<ZAI, double >::iterator it ;
+
- for(int k=0; k < InsideStep; k++)
+ for(it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
{
- double ESigmaN = 0;
- for (int j = 0; j < (int)fReverseMatrixIndex.size() ; j++)
- ESigmaN -= FissionXSMatrix[i][j][j]*NEvolutionMatrix[j][0]*1.6e-19*FissionEnergy[j][0];
- // Update Flux
- double Flux_k = Power/ESigmaN;
-
- if(k==0)
- Flux[i]=Flux_k;
+ /// Need TO change with FP managment
+ map<ZAI, int >::iterator it2;
- BatemanMatrix = BatemanReactionMatrix;
- BatemanMatrix *= Flux_k;
- BatemanMatrix += fDecayMatrix ;
- BatemanMatrix *= TStepMax/InsideStep ;
+ if( (*it).first.Z() < fZAIThreshold )
+ it2 = fMatrixIndex.find( ZAI(-2,-2,-2) );
+ else it2 = fMatrixIndex.find( (*it).first );
+ if(it2 == fMatrixIndex.end() ) //If not in index should be TMP, can't be fast decay for new Fuel !!!
+ it2 = fMatrixIndex.find( ZAI(-3,-3,-3) );
- TMatrixT<double> IdMatrix = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size());
- for(int j = 0; j < (int)fReverseMatrixIndex.size(); j++)
- for(int k = 0; k < (int)fReverseMatrixIndex.size(); k++)
- {
- if(k == j) IdMatrix[j][k] = 1;
- else IdMatrix[j][k] = 0;
- }
+ NMatrix[ (*it2).second ][0] = (*it).second ;
- TMatrixT<double> BatemanMatrixDL = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size()); // Order 0 Term from the DL : Id
- TMatrixT<double> BatemanMatrixDLTermN = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size()); // Addind it;
- {
- BatemanMatrix *= TStepMax ;
- BatemanMatrixDLTermN = IdMatrix;
- BatemanMatrixDL = BatemanMatrixDLTermN;
- int j = 1;
- double NormN;
-
- do
- {
- TMatrixT<double> BatemanMatrixDLTermtmp = TMatrixT<double>(fReverseMatrixIndex.size(),fReverseMatrixIndex.size()); // Adding it;
- BatemanMatrixDLTermtmp = BatemanMatrixDLTermN;
-
- BatemanMatrixDLTermN.Mult(BatemanMatrixDLTermtmp, BatemanMatrix );
-
- BatemanMatrixDLTermN *= 1./j;
- BatemanMatrixDL += BatemanMatrixDLTermN;
-
- NormN = 0;
- for(int m = 0; m < (int)fReverseMatrixIndex.size(); m++)
- for(int n = 0; n < (int)fReverseMatrixIndex.size(); n++)
- NormN += BatemanMatrixDLTermN[m][n]*BatemanMatrixDLTermN[m][n];
- j++;
-
- } while ( NormN != 0 );
- }
- NEvolutionMatrix = BatemanMatrixDL * NEvolutionMatrix ;
- }
- NMatrix.push_back(NEvolutionMatrix);
-
-
- timevector[i+1] = timevector[i] + TStepMax;
-
- BatemanMatrix.Clear();
- BatemanReactionMatrix.Clear();
- NEvolutionMatrix.Clear();
-
-
- }
- DBGL
- FissionXSMatrix.push_back(GetFissionXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
- CaptureXSMatrix.push_back(GetCaptureXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
- n2nXSMatrix.push_back(Getn2nXsMatrix(XSSet, DBTimeStep[NStep-1])); //Feel the reaction Matrix
-
-
- EvolutionData GeneratedDB = EvolutionData(GetLog());
-
- double ESigmaN = 0;
- for (int j = 0; j < (int)fReverseMatrixIndex.size() ; j++)
- ESigmaN -= FissionXSMatrix.back()[j][j]*NMatrix.back()[j][0]*1.6e-19*FissionEnergy[j][0];
-
- Flux[NStep-1] = Power/ESigmaN;
-
- GeneratedDB.SetFlux( new TGraph(NStep, timevector, Flux) );
-
- for(int i = 0; i < (int)fReverseMatrixIndex.size(); i++)
- {
- double ZAIQuantity[NMatrix.size()];
- double FissionXS[NStep];
- double CaptureXS[NStep];
- double n2nXS[NStep];
- for(int j = 0; j < (int)NMatrix.size(); j++)
- ZAIQuantity[j] = (NMatrix[j])[i][0];
-
- for(int j = 0; j < NStep; j++)
- {
- FissionXS[j] = FissionXSMatrix[j][i][i];
- CaptureXS[j] = CaptureXSMatrix[j][i][i];
- n2nXS[j] = n2nXSMatrix[j][i][i];
}
- 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.CaptureXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, CaptureXS)));
- GeneratedDB.n2nXSInsert(pair<ZAI, TGraph*> (fReverseMatrixIndex[i], new TGraph(NStep, timevector, n2nXS)));
+ isotopicquantity.clear();
}
- GeneratedDB.SetPower(Power );
- GeneratedDB.SetHeavyMetalMass(M);
- GeneratedDB.SetReactorType(ReactorType );
- GeneratedDB.SetCycleTime(cycletime);
-
- for (int i = 0; i < (int) FissionXSMatrix.size(); i++)
- {
- FissionXSMatrix[i].Clear();
- CaptureXSMatrix[i].Clear();
- n2nXSMatrix[i].Clear();
- }
- FissionXSMatrix.clear();
- CaptureXSMatrix.clear();
- n2nXSMatrix.clear();
- DBGL
- return GeneratedDB;
}
@@ -298,4 +145,3 @@ EvolutionData IM_Matrix_Decay::GenerateEvolutionData(IsotopicVector isotopicvect
-
diff --git a/source/trunk/Model/Irradiation/IM_Matrix_Decay.hxx b/source/trunk/Model/Irradiation/IM_Matrix_Decay.hxx
index e7d321647..f38ee6694 100644
--- a/source/trunk/Model/Irradiation/IM_Matrix_Decay.hxx
+++ b/source/trunk/Model/Irradiation/IM_Matrix_Decay.hxx
@@ -11,6 +11,7 @@
@version 2.0
*/
#include "IrradiationModel.hxx"
+#include "TMatrix.h"
using namespace std;
@@ -61,21 +62,15 @@ class IM_Matrix_Decay : public IrradiationModel
//@}
- /// virtual method called to perform the irradiation calculation using a set of cross section.
- /*!
- Perform the Irradiation Calcultion using the XSSet data
- \param IsotopicVector IV isotopic vector to irradiate
- \param EvolutionData XSSet set of corss section to use to perform the evolution calculation
- */
- virtual EvolutionData GenerateEvolutionData(IsotopicVector IV, EvolutionData XSSet, double Power, double cycletime);
- //}
-
+ IsotopicVector GetDecay(IsotopicVector Mother_IV, time);
+ TMatrixT<double> ExponentialCalculation(TMatrixT<double> myMatrix);
private :
-
+ TMatrixT<double> fExponentialDecayMatrix; //!< Matrix with half life for each nuclei
+
};
--
GitLab