CLASS
1.1
|
00001 #include "CLASS.hxx" 00002 #include "Reactor.hxx" 00003 #include "TreatmentFactory.hxx" 00004 #include "Defines.hxx" 00005 #include "IsotopicVector.hxx" 00006 #include "LogFile.hxx" 00007 00008 #include <ctime> 00009 #include "time.h" 00010 #include <cmath> 00011 #include <iomanip> 00012 #include <fstream> 00013 #include <iostream> 00014 #include <algorithm> 00015 #include <omp.h> 00016 #include <TFile.h> 00017 #include <TTree.h> 00018 //#include <gsl/gsl_vector.h> 00019 //#include <gsl/gsl_multimin.h> 00020 00021 00022 //________________________________________________________________________ 00023 // 00024 // CLASS 00025 // 00026 // 00027 // 00028 // 00029 //________________________________________________________________________ 00030 00031 00032 float random(float a, float b) //peak random numebr between a and b 00033 { 00034 float range = pow(2., 31); 00035 srand(time(NULL)); //initialize the srand 00036 return (float)a + (float)(b-a)*rand()/range; 00037 } 00038 00039 //________________________________________________________________________ 00040 CLASS::CLASS() 00041 { 00042 DBGL; 00043 fPrintStep = (long int)(3600*24*365.4); // One Step per Month 00044 fAbsoluteTime = 0; 00045 fTreatmentFactoryLastIndex = 0; 00046 fReactorLastIndex = 0; 00047 fStockManagement = true; 00048 fBuildingMethod = 0; 00049 00050 fOutputName = "CLASS_Default.root"; 00051 string logname = "CLASS.log"; 00052 fLog = new LogFile("CLASS.log"); 00053 00054 DBGL; 00055 } 00056 00057 //________________________________________________________________________ 00058 CLASS::CLASS(long int abstime) 00059 { 00060 DBGL; 00061 fPrintStep = (long int)(3600*24*365.4); // One Step per Month 00062 fAbsoluteTime = abstime; 00063 fTreatmentFactoryLastIndex = 0; 00064 fReactorLastIndex = 0; 00065 fStockManagement = true; 00066 fBuildingMethod = 0; 00067 00068 fOutputName = "CLASS_Default.root"; 00069 string logname = "CLASS.log"; 00070 fLog = new LogFile("CLASS.log"); 00071 00072 DBGL; 00073 } 00074 00075 //________________________________________________________________________ 00076 CLASS::CLASS(string name, long int abstime) 00077 { 00078 DBGL; 00079 fPrintStep = (long int)(3600*24*365.4); // One Step per Month 00080 fAbsoluteTime = abstime; 00081 fTreatmentFactoryLastIndex = 0; 00082 fReactorLastIndex = 0; 00083 fStockManagement = true; 00084 fBuildingMethod = 0; 00085 fOutputName = name; 00086 string logname = "CLASS.log"; 00087 fLog = new LogFile("CLASS.log"); 00088 00089 DBGL; 00090 } 00091 //________________________________________________________________________ 00092 CLASS::~CLASS() 00093 { 00094 DBGL; 00095 DBGL; 00096 } 00097 00098 //________________________________________________________________________ 00099 void CLASS::AddTreatmentFactory(TreatmentFactory* treatmentfactory) 00100 { 00101 DBGL; 00102 00103 fTreatmentFactory.push_back(treatmentfactory); 00104 fTreatmentFactory.back()->SetParc(this); 00105 fTreatmentFactory.back()->SetLog(fLog); 00106 fTreatmentFactory.back()->SetStockManagement(fStockManagement); 00107 fTreatmentFactoryLastIndex++; 00108 fTreatmentFactoryIndex.push_back(fTreatmentFactoryLastIndex); 00109 DBGL; 00110 } 00111 00112 //________________________________________________________________________ 00113 void CLASS::AddReactor(Reactor* reactor) 00114 { 00115 DBGL; 00116 fReactor.push_back(reactor); 00117 fReactor.back()->SetParc(this); 00118 fReactor.back()->SetLog(fLog); 00119 fReactorLastIndex++; 00120 fReactorIndex.push_back(fReactorLastIndex); 00121 DBGL; 00122 } 00123 00124 //________________________________________________________________________ 00125 void CLASS::RemoveReactor() 00126 { 00127 DBGL; 00128 for(int i = (int)fReactor.size()-1; i >= 0; i--) 00129 { 00130 if(fAbsoluteTime == fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime()) 00131 { 00132 fReactor[i]->GetAssociedTreatmentFactory()->AddIVCooling(fReactor[i]->GetIVReactor()); 00133 delete (fReactor[i]); 00134 fReactor.erase(fReactor.begin()+i); 00135 fReactorIndex.erase(fReactorIndex.begin()+i); 00136 } 00137 } 00138 DBGL; 00139 } 00140 00141 //________________________________________________________________________ 00142 //____________________________ Building Method ___________________________ 00143 //________________________________________________________________________ 00144 00145 void CLASS::UpdateParcStock() 00146 { 00147 DBGL; 00148 00149 fParcStock.clear(); 00150 00151 for(int i=0; i< (int)fTreatmentFactory.size(); i++ ) 00152 { 00153 for(int j = 0; j < (int)fTreatmentFactory[i]->GetIVStock().size(); j++ ) 00154 { 00155 fParcStock.insert( pair< pair< int ,int > , IsotopicVector > 00156 ( pair< int ,int >(i, j), fTreatmentFactory[i]->GetIVStock()[j]) ) ; 00157 } 00158 } 00159 00160 DBGL; 00161 } 00162 00163 //________________________________________________________________________ 00164 IsotopicVector CLASS::BuildIsotopicVector(IsotopicVector isotopicvector) 00165 { 00166 DBGL; 00167 UpdateParcStock(); 00168 00169 IsotopicVector BuildIV; 00170 if (fStockManagement == false ) 00171 return BuildIV; 00172 00173 if(fBuildingMethod == 0) 00174 return MonteCarloBuild(isotopicvector); 00175 00176 // else if(fBuildingMethod == 1) 00177 // return MinimizationBuild(isotopicvector); 00178 else 00179 { 00180 cout << "!!Warning!! !!!CLASS!!! Wrong BUilding Method !\n" << endl; 00181 fLog->fLog << "!!Warning!! !!!CLASS!!! Wrong BUilding Method !\n" << endl; 00182 exit(1); 00183 } 00184 00185 DBGL; 00186 } 00187 00188 00189 //________________________________________________________________________ 00190 void CLASS::DumpParcStock() 00191 { 00192 DBGL; 00193 if(fStockManagement == false) return; 00194 for(int i = 0; i < (int)fTreatmentFactory.size(); i++) 00195 fTreatmentFactory[i]->ClearIVStock(); 00196 00197 00198 map< pair<int,int>, IsotopicVector >::iterator it; 00199 for(it = fParcStock.begin(); it != fParcStock.end(); it++) 00200 fTreatmentFactory[(*it).first.first]->AddIVStock((*it).second); 00201 00202 DBGL; 00203 } 00204 00205 00206 00207 00208 00209 //________________________________________________________________________ 00210 00211 //________________________________________________________________________ 00212 00213 /*struct MinimizationParameter { 00214 int z; 00215 vector<int> FactoryIndex; 00216 vector<int> StockIndex; 00217 map< pair<int,int>, IsotopicVector > ParcStock; 00218 IsotopicVector isotopicvector; 00219 IsotopicVector TotalStock; 00220 00221 }; 00222 00223 00224 double my_f (const gsl_vector *v, void *params) 00225 { 00226 00227 MinimizationParameter *par = (MinimizationParameter *)params; 00228 int z = (*par).z; 00229 vector<int> FactoryIndex = (*par).FactoryIndex; 00230 vector<int> StockIndex = (*par).StockIndex; 00231 map< pair<int,int>, IsotopicVector > ParcStock = (*par).ParcStock; 00232 IsotopicVector isotopicvector = (*par).isotopicvector; 00233 IsotopicVector TotalStock = (*par).TotalStock; 00234 00235 double x[FactoryIndex.size()]; 00236 00237 IsotopicVector IVSum; 00238 for (int i =0; i < (int)FactoryIndex.size(); i++) 00239 { 00240 x[i] = gsl_vector_get(v, i); 00241 if(x[i]<0) x[i] = -x[i]; 00242 } 00243 00244 00245 for (int i = 0; i < (int)FactoryIndex.size(); i++) 00246 { 00247 map< pair<int,int>, IsotopicVector >::iterator it; 00248 it = ParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] )); 00249 IVSum += (*it).second.GetAtomicComposition(z)*x[i]/Norme((*it).second.GetAtomicComposition(z)); 00250 } 00251 return Distance(IVSum,isotopicvector.GetAtomicComposition(z)); 00252 } 00253 00254 // The gradient of f, df = (df/dx, df/dy). 00255 void my_df (const gsl_vector *v, void *params, gsl_vector *df) 00256 { 00257 MinimizationParameter *par = (MinimizationParameter *)params; 00258 int z = (*par).z; 00259 vector<int> FactoryIndex = (*par).FactoryIndex; 00260 vector<int> StockIndex = (*par).StockIndex; 00261 map< pair<int,int>, IsotopicVector > ParcStock = (*par).ParcStock; 00262 IsotopicVector isotopicvector = (*par).isotopicvector; 00263 IsotopicVector TotalStock = (*par).TotalStock; 00264 00265 00266 double x[FactoryIndex.size()]; 00267 00268 for (int i =0; i < (int)FactoryIndex.size(); i++) 00269 { 00270 x[i] = gsl_vector_get(v, i); 00271 if(x[i]<0) x[i] = -x[i]; 00272 } 00273 for (int i = 0; i < (int)FactoryIndex.size(); i++) 00274 { 00275 double Sum = 0; 00276 map< pair<int,int>, IsotopicVector >::iterator it_i; 00277 it_i = ParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] )); 00278 map<ZAI ,double > isotopicquantity = TotalStock.GetAtomicComposition(z).GetIsotopicQuantity(); 00279 for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++) 00280 { 00281 Sum += 2*(*it_i).second.GetAtomicComposition(z).GetZAIIsotopicQuantity((*it).first) 00282 *isotopicvector.GetAtomicComposition(z).GetZAIIsotopicQuantity((*it).first); 00283 00284 for (int j = 0; j < (int)FactoryIndex.size(); j++) 00285 { 00286 //if(i == j) j++; 00287 00288 map< pair<int,int>, IsotopicVector >::iterator it_j; 00289 it_j = ParcStock.find(pair<int,int>(FactoryIndex[j], StockIndex[j] )); 00290 00291 Sum -= 2*x[j]*(*it_i).second.GetZAIIsotopicQuantity((*it).first) 00292 *(*it_j).second.GetZAIIsotopicQuantity((*it).first); 00293 } 00294 } 00295 gsl_vector_set(df, i, Sum); 00296 00297 } 00298 } 00299 00300 // Compute both f and df together. 00301 void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df) 00302 { 00303 *f = my_f(x, params); 00304 my_df(x, params, df); 00305 } 00306 00307 00308 00309 IsotopicVector CLASS::MinimizationBuild(IsotopicVector isotopicvector ) 00310 { 00311 DBGL; 00312 IsotopicVector BuildIV; 00313 vector<int> FactoryIndex; 00314 vector<int> StockIndex; 00315 00316 map< pair<int,int>, IsotopicVector >::iterator it; 00317 00318 for(it = fParcStock.begin(); it != fParcStock.end(); it++ ) 00319 { 00320 FactoryIndex.push_back((*it).first.first); 00321 StockIndex.push_back((*it).first.second); 00322 } 00323 00324 00325 if( (int)FactoryIndex.size() <= 3) return BuildIV; 00326 00327 for(int k = 0; k < (int)isotopicvector.GetAtomicSpecies().size(); k++ ) // Loop on the Atomic Species 00328 { 00329 int z = isotopicvector.GetAtomicSpecies()[k]; // Get the Atomic Species Number 00330 IsotopicVector BuildIVz; 00331 IsotopicVector IVAtomicCompo = isotopicvector.GetAtomicComposition(z); 00332 //double delta = 1e-3*Norme(IVAtomicCompo); 00333 00334 00335 00336 00337 00338 size_t iter = 0; 00339 int status; 00340 00341 const gsl_multimin_fdfminimizer_type *T; 00342 gsl_multimin_fdfminimizer *s; 00343 00344 // Index of IV 00345 MinimizationParameter par ; 00346 par.z = z; 00347 par.FactoryIndex = FactoryIndex; 00348 par.StockIndex = StockIndex; 00349 par.ParcStock = fParcStock; 00350 par.isotopicvector = isotopicvector; 00351 par.TotalStock = fTotalStock; 00352 //double par[5] = { 1.0, 2.0, 10.0, 20.0, 30.0 }; 00353 00354 gsl_vector *x; 00355 gsl_multimin_function_fdf my_func; 00356 00357 my_func.n = (int)FactoryIndex.size(); // number of function components 00358 my_func.f = my_f; 00359 my_func.df = my_df; 00360 my_func.fdf = my_fdf; 00361 my_func.params = ∥ 00362 00363 // Starting point, x = (5,7) 00364 x = gsl_vector_alloc (FactoryIndex.size()); 00365 for (int i = 0; i < (int)FactoryIndex.size(); i++) 00366 { 00367 map< pair<int,int>, IsotopicVector >::iterator it; 00368 it = fParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] )); 00369 gsl_vector_set (x, i, Norme((*it).second.GetAtomicComposition(z))/2. ); 00370 } 00371 T = gsl_multimin_fdfminimizer_conjugate_fr; 00372 s = gsl_multimin_fdfminimizer_alloc (T, FactoryIndex.size()); 00373 00374 gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-1); 00375 00376 do 00377 { 00378 iter++; 00379 status = gsl_multimin_fdfminimizer_iterate (s); 00380 00381 if (status) 00382 break; 00383 00384 status = gsl_multimin_test_gradient (s->gradient, 1e-3); 00385 //status = gsl_multimin_test_size (0, 1e-3); 00386 if (status == GSL_SUCCESS) 00387 printf ("Minimum found !\n"); 00388 00389 } 00390 while (status == GSL_CONTINUE && iter < 100); 00391 00392 00393 for (int i = 0; i < (int)FactoryIndex.size(); i++) 00394 { 00395 map< pair<int,int>, IsotopicVector >::iterator it; 00396 it = fParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] )); 00397 double Contribution_i = gsl_vector_get (s->x, i); 00398 BuildIVz += (*it).second.GetAtomicComposition(z)*Contribution_i/Norme((*it).second.GetAtomicComposition(z)); 00399 (*it).second -= (*it).second.GetAtomicComposition(z)*Contribution_i/Norme((*it).second.GetAtomicComposition(z)); 00400 } 00401 00402 00403 gsl_multimin_fdfminimizer_free (s); 00404 gsl_vector_free (x); 00405 00406 BuildIV += BuildIVz; 00407 } 00408 return BuildIV; 00409 DBGL; 00410 } 00411 */ 00412 00413 //________________________________________________________________________ 00414 IsotopicVector CLASS::MonteCarloBuild(IsotopicVector isotopicvector ) 00415 { 00416 00417 IsotopicVector BuildIV; 00418 00419 00420 int z = 0; 00421 00422 vector<int> FactoryIndex; 00423 vector<int> StockIndex; 00424 map< pair<int,int>, IsotopicVector >::iterator it; 00425 00426 for(it = fParcStock.begin(); it != fParcStock.end(); it++ ) 00427 { 00428 FactoryIndex.push_back((*it).first.first); 00429 StockIndex.push_back((*it).first.second); 00430 } 00431 00432 if( (int)FactoryIndex.size() == 0) return BuildIV; 00433 00434 for(int k = 0; k < (int)isotopicvector.GetAtomicSpecies().size(); k++ ) // Loop on the Atomic Species 00435 { 00436 z = isotopicvector.GetAtomicSpecies()[k]; // Get the Atomic Species Number 00437 00438 IsotopicVector IVAtomicCompo = isotopicvector.GetAtomicComposition(z); 00439 IsotopicVector BuildIVz; 00440 double delta = 1e-3*Norme(IVAtomicCompo); 00441 00442 00443 00444 vector<int> FactoryIndextmp = FactoryIndex; 00445 vector<int> StockIndextmp = StockIndex; 00446 00447 double DistanceN = Norme(IVAtomicCompo); // Distance at the Nst iteration 00448 double DistanceN_1 = Norme(IVAtomicCompo); // Distance at the N-1 st iteration 00449 00450 while(FactoryIndextmp.size() !=0) 00451 { 00452 int IVIndex = FactoryIndextmp.size(); 00453 00454 // Find a random IV contribution 00455 while(IVIndex == (int)FactoryIndextmp.size()) IVIndex = (int)random(0,FactoryIndextmp.size()); 00456 00457 map< pair<int,int>, IsotopicVector >::iterator it; 00458 it = fParcStock.find(pair<int,int>(FactoryIndextmp[IVIndex], StockIndextmp[IVIndex] )); 00459 IsotopicVector IVToAdd = (*it).second.GetAtomicComposition(z) ; 00460 00461 00462 IsotopicVector BuildIVzTmp = BuildIVz + IVToAdd/Norme(IVToAdd) * delta; // Build the new IVBuild 00463 00464 DistanceN = Distance( BuildIVzTmp, isotopicvector.GetAtomicComposition(z) ); //Progress determination 00465 00466 if(DistanceN > DistanceN_1 || delta > Norme(IVToAdd) ) // Not Good Peak remove it from Index 00467 { 00468 FactoryIndextmp.erase(FactoryIndextmp.begin()+IVIndex); 00469 StockIndextmp.erase(StockIndextmp.begin()+IVIndex); 00470 } 00471 else // Good Peak, Add the contribution to IVBuild and remove from stock... 00472 { 00473 // Validation of Values.... 00474 DistanceN_1 = DistanceN; 00475 BuildIVz = BuildIVzTmp; 00476 00477 // remove contribution from stock 00478 map< pair<int,int>, IsotopicVector >::iterator it; 00479 it = fParcStock.find(pair<int,int>(FactoryIndextmp[IVIndex], StockIndextmp[IVIndex] )); 00480 00481 (*it).second -= IVToAdd/Norme(IVToAdd) * delta; 00482 } 00483 } 00484 BuildIV += BuildIVz; 00485 } 00486 00487 DBGL; 00488 return BuildIV; 00489 } 00490 00491 00492 //________________________________________________________________________ 00493 void CLASS::BuildTimeVector(long int t) 00494 { 00495 DBGL; 00496 fTimeStep.clear(); 00497 fTimeStep.insert( pair<long int ,int>(t,1) ); 00498 00499 //********* Printing Step *********// 00500 { 00501 long int step = 0; 00502 if(step >= fAbsoluteTime ) 00503 fTimeStep.insert( pair<long int ,int>(step,1) ); 00504 step += fPrintStep; 00505 do 00506 { 00507 00508 if(step >= fAbsoluteTime ) 00509 fTimeStep.insert( pair<long int ,int>(step,1) ); 00510 step += fPrintStep; 00511 } 00512 while( step < t ); 00513 } 00514 00515 for(int i = 0; i < (int)fReactor.size();i++) 00516 { 00517 long int step = fReactor[i]->GetCreationTime(); 00518 long int coolingstep = fReactor[i]->GetAssociedTreatmentFactory()->GetCoolingTime(); 00519 long int seprationstep = fReactor[i]->GetAssociedTreatmentFactory()->GetSeparationTime(); 00520 00521 //********* Reactor Evolution Step *********// 00522 // set destruction of a reactor 00523 if( (fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime() > fAbsoluteTime) && 00524 (fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime()<t) ) 00525 { 00526 pair< map<long int, int>::iterator, bool > IResult ; 00527 IResult = fTimeStep.insert( pair<long int ,int>(fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime(),2) ); 00528 if( IResult.second == false ) IResult.first->second |= 2; 00529 } 00530 00531 // Set End of reactor cycle 00532 if(step >= fAbsoluteTime && step <= t && step < fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime()) 00533 { 00534 pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(step,4) ); 00535 if( IResult.second == false ) IResult.first->second |= 4; 00536 } 00537 00538 00539 //********* TF Evolution Step *********// 00540 step += fReactor[i]->GetCycleTime(); 00541 do 00542 { 00543 00544 if(step > fAbsoluteTime && step <= t && step < fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime()) 00545 { // Set End of reactor cycle 00546 pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(step,4) ); 00547 if( IResult.second == false ) IResult.first->second |= 4; 00548 } 00549 00550 if(step >= fAbsoluteTime && step + coolingstep <= t) // Set End of Cooling 00551 { 00552 pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(step+coolingstep,8) ); 00553 if( IResult.second == false ) IResult.first->second |= 8; 00554 } 00555 00556 if(step >= fAbsoluteTime && step + coolingstep + seprationstep <= t) // Set end of Separation 00557 { 00558 pair< map<long int, int>::iterator, bool > IResult; 00559 IResult = fTimeStep.insert( pair<long int ,int>(step+coolingstep+seprationstep,16) ); 00560 if( IResult.second == false ) IResult.first->second |= 16; 00561 } 00562 step += fReactor[i]->GetCycleTime(); 00563 } 00564 while(step <= t && step <= fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime() ); 00565 } 00566 00567 //********* TF Evolution Step *********// 00568 //*** In Case of Evolution Restart ****// 00569 for(int i =0; i < (int)fTreatmentFactory.size(); i++) 00570 { 00571 00572 for(int j = 0; j<(int)fTreatmentFactory[i]->GetIVCooling().size(); j++ )// Set End of Cooling 00573 { 00574 if(fTreatmentFactory[i]->GetCoolingStartingTime()[j] + fTreatmentFactory[i]->GetCoolingTime() > fAbsoluteTime ) 00575 { 00576 pair< map<long int, int>::iterator, bool > IResult; 00577 IResult = fTimeStep.insert( pair<long int ,int>(fTreatmentFactory[i]->GetCoolingStartingTime()[j] + fTreatmentFactory[i]->GetCoolingTime(),8) ); 00578 if( IResult.second == false ) IResult.first->second |= 8; 00579 } 00580 } 00581 for(int j = 0; j< (int)fTreatmentFactory[i]->GetIVSeparating().size(); j++ )// Set end of Separation 00582 { 00583 if(fTreatmentFactory[i]->GetSeparatingStartingTime()[j] + fTreatmentFactory[i]->GetSeparationTime() > fAbsoluteTime ) 00584 { 00585 pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(fTreatmentFactory[i]->GetSeparatingStartingTime()[j] + fTreatmentFactory[i]->GetSeparationTime(),16) ); 00586 if( IResult.second == false ) IResult.first->second |= 16; 00587 } 00588 } 00589 } 00590 00591 00592 00593 //****** Print the Time Index ******// 00594 ofstream TimeStepfile("CLASS_TimeStep", ios_base::app); // Open the File 00595 00596 if(!TimeStepfile) 00597 { 00598 cout << "!!Warning!! !!!CLASS!!! Can't open \" CLASS_TimeStep \"\n" << endl; 00599 fLog->fLog << "!!Warning!! !!!CLASS!!! Can't open \" CLASS_TimeStep \"\n" << endl; 00600 } 00601 for(map<long int ,int >::iterator it = fTimeStep.begin(); it != fTimeStep.end(); it++) 00602 TimeStepfile << (*it).first/365.4/3600/24 << " " << (*it).second << endl; 00603 00604 DBGL; 00605 } 00606 00607 00608 00609 //________________________________________________________________________ 00610 //___________________________ Evolution Method ___________________________ 00611 //________________________________________________________________________ 00612 00613 void CLASS::TreatmentEvolution() 00614 { 00615 DBGL; 00616 for(int i = 0; i < (int) fTreatmentFactory.size();i++) 00617 fTreatmentFactory[i]->Evolution(fAbsoluteTime); 00618 DBGL; 00619 } 00620 00621 //________________________________________________________________________ 00622 void CLASS::ReactorEvolution() 00623 { 00624 DBGL; 00625 00626 00627 #pragma omp sections 00628 { 00629 #pragma omp section 00630 {TreatmentEvolution();} 00631 00632 #pragma omp section 00633 { 00634 #pragma omp parallel for 00635 for(int i = 0; i < (int)fReactor.size(); i++) 00636 fReactor[i]->Evolution(fAbsoluteTime); 00637 } 00638 } 00639 00640 00641 for(int i = 0; i < (int) fTreatmentFactory.size();i++) 00642 fTreatmentFactory[i]->Dump(); 00643 00644 UpdateParcStock(); 00645 for(int i = 0; i < (int)fReactor.size(); i++) 00646 fReactor[i]->Dump(); 00647 00648 DumpParcStock(); 00649 00650 DBGL; 00651 } 00652 00653 //________________________________________________________________________ 00654 void CLASS::Evolution(long int t) 00655 { 00656 DBGL; 00657 OpenOutputTree(); 00658 OutAttach(); 00659 BuildTimeVector(t); 00660 00661 int Start = time(NULL); 00662 for(map<long int ,int >::iterator it = fTimeStep.begin(); it != fTimeStep.end(); it++) 00663 { 00664 00665 ResetQuantity(); 00666 fAbsoluteTime = (*it).first; 00667 //if(fAbsoluteTime > t) return; 00668 00669 if( (*it).second & 2 || (*it).second & 1 ) 00670 { 00671 ReactorEvolution(); 00672 if((*it).second & 2 ) 00673 RemoveReactor(); 00674 00675 if((*it).second & 2 ) 00676 (*it).second ^= 2; 00677 if((*it).second & 4 ) 00678 (*it).second ^= 4; 00679 if((*it).second & 8 ) 00680 (*it).second ^= 8; 00681 if((*it).second & 16 ) 00682 (*it).second ^= 16; 00683 } 00684 00685 if( (*it).second & 4 ) 00686 { 00687 ReactorEvolution(); 00688 if((*it).second & 4 ) 00689 (*it).second ^= 4; 00690 if((*it).second & 8 ) 00691 (*it).second ^= 8; 00692 if((*it).second & 16 ) 00693 (*it).second ^= 16; 00694 } 00695 00696 if( (*it).second & 8 || (*it).second & 16 ) 00697 { 00698 TreatmentEvolution(); 00699 for(int i = 0; i < (int) fTreatmentFactory.size();i++) 00700 fTreatmentFactory[i]->Dump(); 00701 00702 if((*it).second & 8 ) 00703 (*it).second ^= 8; 00704 if((*it).second & 16 ) 00705 (*it).second ^= 16; 00706 } 00707 00708 00709 00710 if( (*it).second & 1 ) 00711 #pragma omp single 00712 { 00713 // Print(); 00714 // Write(); 00715 fOutT->Fill(); 00716 //fTotalStock.Print(); 00717 ProgressPrintout(Start, t); 00718 } 00719 } 00720 #pragma omp single 00721 {CloseOutputTree();} 00722 DBGL; 00723 } 00724 void CLASS::ProgressPrintout(int starttime, long int t) 00725 { 00726 DBGL; 00727 int TimeNow = time(NULL); 00728 int Spent = difftime(TimeNow, starttime); 00729 double Time = fAbsoluteTime/3600/24/365.4; 00730 double Total = t/3600/24/365.4; 00731 double Remain = (Spent/Time * Total - Spent); 00732 cout << " " << flush ; 00733 cout << "\rProcessed " << setprecision(4) << Time << " / " << setprecision(4) << Total << " STEP " 00734 << setprecision(4) << (Time/Total*100.0) << "%) --- I still need : " 00735 << (long int)Remain/60 << " min "<< (long int)(Remain)%60 << " sec to finish ! \r" << flush; 00736 DBGL; 00737 } 00738 00739 //________________________________________________________________________ 00740 //______________________________ Out Method ______________________________ 00741 //________________________________________________________________________ 00742 00743 void CLASS::ResetQuantity() 00744 { 00745 DBGL; 00746 IsotopicVector EmptyIV; 00747 fTotalInReactor = EmptyIV; 00748 fTotalWaste = EmptyIV; 00749 fTotalStock = EmptyIV; 00750 // fTotalGodIncome = EmptyIV; 00751 fTotalCooling = EmptyIV; 00752 fTotalSeparating = EmptyIV; 00753 fIVInCycleTotal = EmptyIV; 00754 fIVTotal = EmptyIV; 00755 DBGL; 00756 } 00757 00758 //________________________________________________________________________ 00759 void CLASS::OpenOutputTree() 00760 { 00761 DBGL; 00762 00763 cout << "Opening : " << fOutputName << " ..."; 00764 fLog->fLog << "Opening : " << fOutputName << " ..."; 00765 fOutTree = new TFile(fOutputName.c_str(),"UPDATE"); 00766 00767 if(!fOutTree) 00768 { 00769 cout << "\nCould not open " << fOutputName <<endl; 00770 fLog->fLog << "\nCould not open " << fOutputName <<endl; 00771 exit(-1); 00772 } 00773 cout << "\t ...O.K." << endl; 00774 00775 fOutT = new TTree("Data","Data Tree"); 00776 cout << "Creating Data Tree..."; 00777 fLog->fLog << "Creating Data Tree..."; 00778 if(!fOutTree) 00779 { 00780 cout << "\nCould not create Data Tree in " << fOutputName << endl; 00781 fLog->fLog << "\nCould not create Data Tree in " << fOutputName << endl; 00782 exit(-1); 00783 } 00784 cout << "\t ...O.K!" << endl; 00785 fLog->fLog << "\t ...O.K!" << endl; 00786 // fOutT->SetAutoFlush(1000000); 00787 // fOutT->SetAutoSave(100000); 00788 } 00789 void CLASS::CloseOutputTree() 00790 { 00791 DBGL; 00792 00793 fOutTree->ls(); 00794 cout << "Writing outTree " << fOutputName << endl; 00795 fLog->fLog << "Writing outTree " << fOutputName << endl; 00796 fOutTree->Write(); 00797 00798 if(fOutTree->IsOpen()) { 00799 cout << "Deleting outTree : " << endl; 00800 fLog->fLog << "Deleting outTree : " << endl; 00801 delete fOutT; 00802 cout << "Closing file : " << fOutputName <<endl; 00803 fLog->fLog << "Closing file : " << fOutputName <<endl; 00804 fOutTree-> Close(); 00805 delete fOutTree; 00806 } else { 00807 cout << "File was not opened " << fOutputName << endl; 00808 fLog->fLog << "File was not opened " << fOutputName << endl; 00809 exit(-1); 00810 } 00811 } 00812 //________________________________________________________________________ 00813 void CLASS::OutAttach() 00814 { 00815 DBGL; 00816 ResetQuantity(); 00817 //Branch Absolut Time 00818 fOutT->Branch("AbsTime",&fAbsoluteTime,"AbsoluteTime/L"); 00819 00820 // Branch the Sum IV 00821 00822 00823 fOutT->Branch("TF_WASTE.", "IsotopicVector", &fTotalWaste); 00824 fOutT->Branch("TF_STOCK.", "IsotopicVector", &fTotalStock); 00825 fOutT->Branch("TF_SEPARATING.", "IsotopicVector", &fTotalSeparating); 00826 fOutT->Branch("TF_COOLING.", "IsotopicVector", &fTotalCooling); 00827 fOutT->Branch("GODINCOME.", "IsotopicVector", &fTotalGodIncome); 00828 fOutT->Branch("INTOTALINREACTOR.", "IsotopicVector", &fTotalInReactor); 00829 fOutT->Branch("INCYCLE.", "IsotopicVector", &fIVInCycleTotal); 00830 fOutT->Branch("TOTAL.", "IsotopicVector", &fIVTotal); 00831 00832 fOutT->Branch("TreatmentFactory", "TreatmentFactory", &fTreatmentFactory.at(0)); 00833 // fOutT->Branch("Reactor", "Reactor", &fReactor.at(0)); 00834 00835 DBGL; 00836 } 00837 00838 //________________________________________________________________________ 00839 void CLASS::Write() 00840 { 00841 DBGL; 00842 string basename = "CLASS"; 00843 00844 //****** Total *****// 00845 IsotopicVector TotalInReactor; 00846 for(int i = 0; i < (int)fReactor.size(); i++) 00847 { 00848 TotalInReactor += fReactor[i]->GetIVReactor(); 00849 } 00850 string RTotal = basename + "_R_TOTAL"; 00851 TotalInReactor.Write(RTotal, fAbsoluteTime); 00852 00853 //****** Total *****// 00854 IsotopicVector TotalWaste; 00855 IsotopicVector TotalGodIncome; 00856 IsotopicVector TotalStock; 00857 IsotopicVector TotalCooling; 00858 IsotopicVector TotalSeparating; 00859 00860 00861 for(int i = 0; i < (int) fTreatmentFactory.size();i++) 00862 { 00863 TotalWaste += fTreatmentFactory[i]->GetIVWaste(); 00864 TotalStock += fTreatmentFactory[i]->GetIVFullStock(); 00865 TotalGodIncome += fTreatmentFactory[i]->GetIVGodIncome(); 00866 00867 for(int j=0; j < (int)fTreatmentFactory[i]->GetIVCooling().size(); j++) 00868 TotalCooling += fTreatmentFactory[i]->GetIVCooling()[j]; 00869 for(int j=0; j < (int)fTreatmentFactory[i]->GetIVSeparating().size(); j++) 00870 TotalSeparating += fTreatmentFactory[i]->GetIVSeparating()[j]; 00871 } 00872 00873 string TFWaste = basename + "_TF_WASTE"; 00874 string TFStock = basename + "_TF_STOCK"; 00875 string TFGodIncome = basename + "_TF_GODINCOME"; 00876 string TFSeparating = basename + "_TF_SPERATING"; 00877 00878 TotalWaste.Write(TFWaste, fAbsoluteTime); 00879 TotalStock.Write(TFStock, fAbsoluteTime); 00880 TotalGodIncome.Write(TFGodIncome, fAbsoluteTime); 00881 TotalSeparating.Write(TFSeparating, fAbsoluteTime); 00882 00883 //---- Total Writing ----/ 00884 00885 string InCycleTotal = basename + "_IN_CYCLE_TOTAL"; 00886 string Total = basename + "_TOTAL"; 00887 00888 IsotopicVector IVInCycleTotal; 00889 IsotopicVector IVTotal; 00890 00891 IVInCycleTotal = TotalCooling + TotalInReactor + TotalSeparating + TotalStock; 00892 IVTotal = IVInCycleTotal + TotalWaste; 00893 00894 IVInCycleTotal.Write(InCycleTotal, fAbsoluteTime); 00895 IVTotal.Write(Total, fAbsoluteTime); 00896 00897 DBGL; 00898 } 00899 00900 //________________________________________________________________________ 00901 void CLASS::Print() 00902 { 00903 DBGL; 00904 for(int i = 0; i < (int) fTreatmentFactory.size();i++) 00905 { 00906 cout << "!!!!!!!!!STEP : " << fAbsoluteTime/(int)(3600*24*365.4) << endl; 00907 cout << "TreatmentFactory : " << endl; 00908 cout << "Cooling "; 00909 cout << fTreatmentFactory[i]->GetIVCooling().size()<< endl; 00910 cout << "Waste " << endl; 00911 fTreatmentFactory[i]->GetIVWaste().Print(); 00912 cout << "Stock :" << endl; 00913 for(int j=0; j < (int)fTreatmentFactory[i]->GetIVStock().size(); j++) 00914 { 00915 cout << j << endl; 00916 fTreatmentFactory[i]->GetIVStock()[j].Print(); 00917 } 00918 } 00919 00920 for(int i = 0; i < (int)fReactor.size(); i++) 00921 { 00922 cout << "Reactor" << endl; 00923 fReactor[i]->GetIVReactor().Print(); 00924 } 00925 DBGL; 00926 }