CoreFlows/Models/src/TransportEquation.cxx

   1 #include "TransportEquation.hxx"
   2 #include "math.h"
   3 #include <fstream>
   4 #include <sstream>
   5
   6 using namespace std;
   7
   8 TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport, MPI_Comm comm):ProblemCoreFlows(comm)
   9 {
  10         if(pEstimate==around1bar300KTransport){
  11                 _Tref=300;
  12                 if(fluid==GasPhase){//Nitrogen pressure 1 bar and temperature 27°C
  13                         _href=3.11e5; //nitrogen enthalpy at 1 bar and 300K
  14                         _cpref=1041;//nitrogen specific heat at constant pressure 1 bar and 300K
  15                         //saturation data for nitrogen at 1 bar and 77K
  16                         _hsatv=0.77e5;//nitrogen vapour enthalpy at saturation at 1 bar
  17                         _hsatl=-1.22e5;//nitrogen liquid enthalpy at saturation at 1 bar
  18                         _rhosatv=4.556;//nitrogen vapour density at saturation at 1 bar
  19                         _rhosatl=806.6;//nitrogen liquid density at saturation at 1 bar
  20                 }
  21                 else{
  22                         //Water at pressure 1 bar and temperature 27°C
  23                         _href=1.127e5; //water enthalpy at 1 bar and 300K
  24                         _cpref=4181;//water specific heat at 1 bar and 300K
  25                         //saturation data for water at 1 bar and 373K
  26                         _hsatv=2.675e6;//Gas enthalpy at saturation at 1 bar
  27                         _hsatl=4.175e5;//water enthalpy at saturation at 1 bar
  28                         _rhosatv=0.6;//Gas density at saturation at 1 bar
  29                         _rhosatl=958;//water density at saturation at 1 bar
  30                 }
  31         }
  32         else{//around155bars600K
  33                 _Tref=618;//=Tsat
  34                 if(fluid==GasPhase){
  35                         _href=2.675e6; //Gas enthalpy at 155 bars and 618K
  36                         _cpref=14001;//Gas specific heat at 155 bar and 618K
  37                 }
  38                 else{//Liquid
  39                         _href=4.175e5;//water enthalpy at 155 bars and 618K
  40                         _cpref=8950;//water specific heat at 155 bar and 618K
  41                 }
  42                 //saturation data for water at 155 bars and 618K
  43                 _hsatv=2.6e6;//Gas enthalpy at saturation at 155 bars
  44                 _hsatl=1.63e6;//water enthalpy at saturation at 155 bars
  45                 _rhosatv=101.9;//Gas density at saturation at 155 bars
  46                 _rhosatl=594.4;//water density at saturation at 155 bars
  47         }
  48         _Ndim=vitesseTransport.size();
  49         _vitesseTransport=Vector(_Ndim);
  50         for(int i=0;i<_Ndim;i++)
  51                 _vitesseTransport[i]=vitesseTransport[i];
  52         _nVar=1;
  53         _dt_transport=0;
  54         _dt_src=0;
  55         _transportMatrixSet=false;
  56         _FECalculation=false;//Only finite volumes available
  57         _rodTemperatureFieldSet=false;
  58         _rodTemperature=0;
  59
  60         _fileName = "SolverlabTransportProblem";
  61     PetscPrintf(PETSC_COMM_WORLD,"\n Transport problem of fluid enthalpy with constant velocity\n");
  62 }
  63
  64 void TransportEquation::initialize()
  65 {
  66         if(_mpi_rank==0)
  67         {
  68                 if(!_initialDataSet)
  69                         throw CdmathException("!!!!!!!!TransportEquation::initialize() set initial data first");
  70                 else if (_VV.getTypeOfField() != CELLS)
  71                         throw CdmathException("!!!!!!!!TransportEquation::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
  72                 else
  73                         PetscPrintf(PETSC_COMM_SELF,"\n Initialising the transport of a fluid enthalpy\n");
  74
  75                 /**************** Field creation *********************/
  76
  77                 //post processing fields used only for saving results
  78                 _TT=Field ("Temperature", CELLS, _mesh, 1);
  79                 _Alpha=Field ("Void fraction", CELLS, _mesh, 1);
  80                 _Rho=Field ("Mixture density", CELLS, _mesh, 1);
  81                 //Construction des champs de post-traitement
  82                 for(int i =0; i<_Nmailles;i++){
  83                         _TT(i)=temperature(_VV(i));
  84                         _Alpha(i)=voidFraction(_VV(i));
  85                         _Rho(i)=density(_Alpha(i));
  86                 }
  87                 if(!_heatPowerFieldSet){
  88                         _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
  89                         for(int i =0; i<_Nmailles; i++)
  90                                 _heatPowerField(i) = _heatSource;
  91                 }
  92                 if(!_rodTemperatureFieldSet){
  93                         _rodTemperatureField=Field("Rod temperature",CELLS,_mesh,1);
  94                         for(int i =0; i<_Nmailles; i++)
  95                                 _rodTemperatureField(i) = _rodTemperature;
  96                 }
  97         }
  98
  99         _globalNbUnknowns = _Nmailles;
 100
 101         /* Vectors creations */
 102         VecCreate(PETSC_COMM_WORLD, &_Hn);
 103         VecSetSizes(_Hn,PETSC_DECIDE,_Nmailles);
 104         VecSetFromOptions(_Hn);
 105         VecGetLocalSize(_Hn, &_localNbUnknowns);
 106
 107         VecDuplicate(_Hn, &_Hk);
 108         VecDuplicate(_Hn, &_Hkm1);
 109         VecDuplicate(_Hn, &_deltaH);
 110         VecDuplicate(_Hn, &_b);//RHS of the linear system: _b=Hn/dt + _b0 + puisance
 111         VecDuplicate(_Hn, &_b0);//part of the RHS that comes from the boundary conditions
 112
 113         if(_mpi_rank == 0)//Process 0 reads and distributes initial data
 114                 for(int i =0; i<_Nmailles;i++)
 115                         VecSetValue(_Hn,i,_VV(i), INSERT_VALUES);
 116         VecAssemblyBegin(_Hn);
 117         VecAssemblyEnd(_Hn);
 118
 119         //creation de la matrice
 120         MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
 121
 122         /* Local sequential vector creation */
 123         if(_mpi_size>1 && _mpi_rank == 0)
 124                 VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Hn_seq);//For saving results on proc 0
 125         VecScatterCreateToZero(_Hn,&_scat,&_Hn_seq);
 126
 127         //Linear solver
 128         KSPCreate(PETSC_COMM_SELF, &_ksp);
 129         KSPSetType(_ksp, _ksptype);
 130         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 131         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 132         KSPGetPC(_ksp, &_pc);
 133         PCSetType(_pc, _pctype);
 134
 135         _initializedMemory=true;
 136         save();//save initial data
 137 }
 138
 139 double TransportEquation::computeTimeStep(bool & stop){
 140         if(!_transportMatrixSet)
 141                 _dt_transport=computeTransportMatrix();
 142
 143         _dt_src=computeRHS();
 144
 145     if(_verbose or _system)
 146         {
 147                 PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
 148         VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
 149         }
 150
 151         stop=false;
 152         return min(_dt_transport,_dt_src);
 153 }
 154 double TransportEquation::computeTransportMatrix(){
 155         MatZeroEntries(_A);
 156         VecZeroEntries(_b0);
 157
 158         if(_mpi_rank == 0)
 159         {
 160                 long nbFaces = _mesh.getNumberOfFaces();
 161                 Face Fj;
 162                 Cell Cell1,Cell2;
 163                 string nameOfGroup;
 164                 double inv_dxi, inv_dxj;
 165                 Vector normale(_Ndim);
 166                 double un, hk;
 167                 PetscInt idm, idn;
 168                 std::vector< int > idCells;
 169                 for (int j=0; j<nbFaces;j++){
 170                         Fj = _mesh.getFace(j);
 171
 172                         // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 173                         idCells = Fj.getCellsId();
 174                         Cell1 = _mesh.getCell(idCells[0]);
 175                         idm = idCells[0];
 176                         if (_Ndim >1){
 177                                 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 178                                         if (j == Cell1.getFacesId()[l]){
 179                                                 for (int idim = 0; idim < _Ndim; ++idim)
 180                                                         normale[idim] = Cell1.getNormalVector(l,idim);
 181                                                 break;
 182                                         }
 183                                 }
 184                         }else{ // _Ndim = 1 : assume that this is normal mesh : the face index increases in positive direction
 185                                 if (Fj.getNumberOfCells()<2) {
 186                                         if (j==0)
 187                                                 normale[0] = -1;
 188                                         else if (j==nbFaces-1)
 189                                                 normale[0] = 1;
 190                                         else
 191                                                 throw CdmathException("TransportEquation::ComputeTimeStep(): computation of normal vector failed");
 192                                 } else if(Fj.getNumberOfCells()==2){
 193                                         if (idCells[0] < idCells[1])
 194                                                 normale[0] = 1;
 195                                         else
 196                                                 normale[0] = -1;
 197                                 }
 198                         }
 199                         //Compute velocity at the face Fj
 200                         un=normale*_vitesseTransport;
 201                         if(abs(un)>_maxvp)
 202                                 _maxvp=abs(un);
 203
 204                         // compute 1/dxi = volume of Ci/area of Fj
 205                         if (_Ndim > 1)
 206                                 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 207                         else
 208                                 inv_dxi = 1/Cell1.getMeasure();
 209
 210                         // If Fj is on the boundary
 211                         if (Fj.getNumberOfCells()==1) {
 212                                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 213                                 {
 214                                         cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 215                                         for(int p=0; p<_Ndim; p++)
 216                                                 cout << normale[p] << ",";
 217                                         cout << ") "<<endl;
 218                                 }
 219                                 nameOfGroup = Fj.getGroupName();
 220
 221                                 if     (_limitField[nameOfGroup].bcType==NeumannTransport || _limitField[nameOfGroup].bcType==OutletTransport ){
 222                                         MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
 223                                 }
 224                                 else if(_limitField[nameOfGroup].bcType==InletTransport   || _limitField[nameOfGroup].bcType==DirichletTransport){
 225                                         if(un>0){
 226                                                 MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
 227                                         }
 228                                         else{
 229                                                 hk=_limitField[nameOfGroup].h;
 230                                                 VecSetValue(_b0,idm,-inv_dxi*un*hk, ADD_VALUES);
 231                                         }
 232                                 }
 233                                 else {//bcType=NoneBCTransport
 234                                         cout<<"!!!!!!!!!!!!!!! Error TransportEquation::computeTransportMatrix() !!!!!!!!!!"<<endl;
 235                                         cout<<"!!!!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<" !!!!!!!!!!!!!! "<<endl;
 236                                         cout<<"Accepted boundary conditions are NeumannTransport "<<NeumannTransport<< " and InletTransport "<< InletTransport <<endl;
 237                                         throw CdmathException("Boundary condition not accepted");
 238                                 }
 239                                 // if Fj is inside the domain
 240                         } else  if (Fj.getNumberOfCells()==2 ){
 241                                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 242                                 {
 243                                         cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 244                                         cout << " ; vecteur normal=(";
 245                                         for(int p=0; p<_Ndim; p++)
 246                                                 cout << normale[p] << ",";
 247                                         cout << "). "<<endl;
 248                                 }
 249                                 Cell2 = _mesh.getCell(idCells[1]);
 250                                 idn = idCells[1];
 251                                 if (_Ndim > 1)
 252                                         inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 253                                 else
 254                                         inv_dxj = 1/Cell2.getMeasure();
 255
 256                                 if(un>0){
 257                                         MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
 258                                         MatSetValue(_A,idn,idm,-inv_dxj*un, ADD_VALUES);
 259                                 }
 260                                 else{
 261                                         MatSetValue(_A,idm,idn,inv_dxi*un, ADD_VALUES);
 262                                         MatSetValue(_A,idn,idn,-inv_dxj*un, ADD_VALUES);
 263                                 }
 264                         }
 265                         else
 266                                 throw CdmathException("TransportEquation::ComputeTimeStep(): incompatible number of cells around a face");
 267                 }
 268         }
 269
 270         _transportMatrixSet=true;
 271
 272         MPI_Bcast(&_maxvp, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
 273         PetscPrintf(PETSC_COMM_WORLD, "Maximum conductivity is %.2e, CFL = %.2f, Delta x = %.2e\n",_maxvp,_cfl,_minl);
 274
 275     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 276         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 277         VecAssemblyBegin(_b0);
 278         VecAssemblyEnd(  _b0);
 279
 280         if(abs(_maxvp)<_precision)
 281                 throw CdmathException("TransportEquation::computeTransportMatrix(): maximum eigenvalue for time step is zero");
 282         else
 283                 return _cfl*_minl/_maxvp;
 284 }
 285 double TransportEquation::computeRHS(){
 286         double rhomin=INFINITY;
 287
 288         VecCopy(_b0,_b);
 289
 290         if(_mpi_rank == 0)
 291         {
 292                 if(_system)
 293                         cout<<"Second membre of transport problem"<<endl;
 294                 for (int i=0; i<_Nmailles;i++){
 295                         VecSetValue(_b,i,_heatTransfertCoeff*(_rodTemperatureField(i)-_TT(i))/_Rho(i),ADD_VALUES);
 296                         VecSetValue(_b,i,_heatPowerField(i)/_Rho(i),ADD_VALUES);
 297                         if(_system)
 298                                 cout<<_heatPowerField(i)/_Rho(i)<<endl;
 299                         if(_Rho(i)<rhomin)
 300                                 rhomin=_Rho(i);
 301                 }
 302         }
 303         VecAssemblyBegin(_b);
 304         VecAssemblyEnd(  _b);
 305
 306     if(_verbose or _system)
 307         {
 308                 PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
 309         VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
 310         }
 311
 312         return rhomin*_cpref/_heatTransfertCoeff;
 313 }
 314 void TransportEquation::updatePrimitives()
 315 {
 316         if(_mpi_size>1){
 317                 VecScatterBegin(_scat,_Hn,_Hn_seq,INSERT_VALUES,SCATTER_FORWARD);
 318                 VecScatterEnd(  _scat,_Hn,_Hn_seq,INSERT_VALUES,SCATTER_FORWARD);
 319         }
 320
 321     if(_verbose or _system)
 322         {
 323                 PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n");
 324         VecView(_Hn,PETSC_VIEWER_STDOUT_WORLD);
 325         }
 326
 327         if(_mpi_rank == 0)
 328         {
 329                 double hi;
 330                 for(int i=0; i<_Nmailles; i++)
 331                 {
 332                         if(_mpi_size>1)
 333                                 VecGetValues(_Hn_seq, 1, &i, &hi);
 334                         else
 335                                 VecGetValues(_Hn    , 1, &i, &hi);
 336                         _VV(i)=hi;
 337                         _TT(i)=temperature(hi);
 338                         _Alpha(i)=voidFraction(hi);
 339                         _Rho(i)=density(_Alpha(i));
 340                 }
 341         }
 342 }
 343
 344 bool TransportEquation::initTimeStep(double dt){
 345         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 346         {
 347                 PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n");
 348                 MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
 349         }
 350
 351     if(_dt>0 and dt>0)
 352     {
 353         //Remove the contribution from dt to prepare for new time step. The diffusion matrix is not recomputed
 354         if(_timeScheme == Implicit)
 355             MatShift(_A,-1/_dt+1/dt);
 356         //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
 357     }
 358     else if(dt>0)//_dt==0, first time step
 359     {
 360                 if(_timeScheme == Implicit)
 361                         MatShift(_A,1/dt);
 362         }
 363     else//dt<=0
 364     {
 365         PetscPrintf(PETSC_COMM_WORLD,"TransportEquation::initTimeStep %.2e = \n",dt);
 366         throw CdmathException("Error TransportEquation::initTimeStep : cannot set time step to zero");
 367     }
 368     //At this stage _b contains _b0 + power + heat exchange
 369     VecAXPY(_b, 1/dt, _Hn);
 370
 371         _dt=dt;
 372
 373         return _dt>0;
 374 }
 375
 376 void TransportEquation::abortTimeStep(){
 377     //Remove contribution of dt to the RHS
 378         VecAXPY(_b,  -1/_dt, _Hn);
 379
 380     //Remove contribution of dt to the matrix
 381         if(_timeScheme == Implicit)
 382                 MatShift(_A,-1/_dt);
 383
 384         _dt = 0;
 385 }
 386
 387 bool TransportEquation::iterateTimeStep(bool &converged)
 388 {
 389         bool stop=false;
 390
 391         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 392                 _maxvp=0;
 393                 computeTimeStep(stop);
 394         }
 395         if(stop){
 396                 converged=false;
 397                 return false;
 398         }
 399         VecAXPY(_b, 1/_dt, _Hn);
 400         if(_system)
 401         {
 402                 PetscPrintf(PETSC_COMM_WORLD,"Vecteur Hn : \n");
 403                 VecView(_Hn,  PETSC_VIEWER_STDOUT_WORLD);
 404                 cout << endl;
 405                 PetscPrintf(PETSC_COMM_WORLD,"Vecteur _b : \n");
 406                 VecView(_b,  PETSC_VIEWER_STDOUT_WORLD);
 407                 PetscPrintf(PETSC_COMM_WORLD,"Matrice A : \n");
 408                 MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
 409         }
 410
 411         if(_timeScheme == Explicit)
 412         {
 413                 MatMult(_A, _Hn, _Hk);
 414                 if(_system)
 415                 {
 416                         PetscPrintf(PETSC_COMM_WORLD,"Nouveau vecteur Hk: \n");
 417                         VecView(_Hk,  PETSC_VIEWER_STDOUT_WORLD);
 418                         cout << endl;
 419                 }
 420                 VecAXPY(_Hk, -1, _b);
 421                 if(_system)
 422                 {
 423                         PetscPrintf(PETSC_COMM_WORLD,"Nouveau vecteur Hk-b: \n");
 424                         VecView(_Hk,  PETSC_VIEWER_STDOUT_WORLD);
 425                         cout << endl;
 426                 }
 427                 VecScale(_Hk, -_dt);
 428                 if(_system)
 429                 {
 430                         PetscPrintf(PETSC_COMM_WORLD,"Nouveau vecteur dt*(Hk-b): \n");
 431                         VecView(_Hk,  PETSC_VIEWER_STDOUT_WORLD);
 432                         cout << endl;
 433                 }
 434                 converged = true;
 435         }
 436         else
 437         {
 438
 439 #if PETSC_VERSION_GREATER_3_5
 440                 KSPSetOperators(_ksp, _A, _A);
 441 #else
 442                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 443 #endif
 444
 445                 if(_conditionNumber)
 446                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 447
 448                 KSPSolve(_ksp, _b, _Hk);
 449                 MatShift(_A,-1/_dt);
 450
 451                 KSPGetIterationNumber(_ksp, &_PetscIts);
 452                 if( _MaxIterLinearSolver < _PetscIts)
 453                         _MaxIterLinearSolver = _PetscIts;
 454                 if(_PetscIts>=_maxPetscIts)
 455                 {
 456                         PetscPrintf(PETSC_COMM_WORLD,"Systeme lineaire : pas de convergence de Petsc. Itérations maximales %d atteintes", _maxPetscIts);
 457                         converged=false;
 458                         return false;
 459                 }
 460                 else
 461                 {
 462                         VecCopy(_Hk, _deltaH);//ici on a deltaH=Hk
 463                         VecAXPY(_deltaH,  -1, _Hkm1);//On obtient deltaH=Hk-Hkm1
 464                         VecNorm(_deltaH,NORM_INFINITY,&_erreur_rel);
 465                         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 466                 }
 467         }
 468
 469         VecCopy(_Hk, _Hkm1);
 470
 471         return true;
 472 }
 473 void TransportEquation::validateTimeStep()
 474 {
 475         VecCopy(_Hk, _deltaH);//ici Hk=Hnp1 donc on a deltaH=Hnp1
 476         VecAXPY(_deltaH,  -1, _Hn);//On obtient deltaH=Hnp1-Hn
 477         VecNorm(_deltaH,NORM_INFINITY,&_erreur_rel);
 478
 479         _isStationary =(_erreur_rel <_precision);
 480
 481         VecCopy(_Hk, _Hn);
 482         VecCopy(_Hk, _Hkm1);
 483
 484         updatePrimitives();
 485
 486         if(_mpi_rank == 0)
 487                 if((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 488                 {
 489                         cout <<"Valeur propre locale max: " << _maxvp << endl;
 490                         //Find minimum and maximum void fractions
 491                         double alphamin=INFINITY;
 492                         double alphamax=-INFINITY;
 493                         for(int i=0; i<_Nmailles; i++)
 494                         {
 495                                 if(_Alpha(i)>alphamax)
 496                                         alphamax=_Alpha(i);
 497                                 if(_Alpha(i)<alphamin)
 498                                         alphamin=_Alpha(i);
 499                         }
 500                         cout<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<endl;
 501                 }
 502
 503         _time+=_dt;
 504         _nbTimeStep++;
 505         save();
 506 }
 507
 508 void TransportEquation::terminate(){
 509         VecDestroy(&_Hn);
 510         VecDestroy(&_Hk);
 511         VecDestroy(&_Hkm1);
 512         VecDestroy(&_deltaH);
 513         VecDestroy(&_b0);
 514         VecDestroy(&_b);
 515         MatDestroy(&_A);
 516         if(_mpi_size>1 && _mpi_rank == 0)
 517                 VecDestroy(&_Hn_seq);
 518 }
 519
 520 void TransportEquation::save(){
 521     PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results at time step number %d \n\n", _nbTimeStep);
 522     *_runLogFile<< "Saving numerical results at time step number "<< _nbTimeStep << endl<<endl;
 523
 524         string resultFile(_path+"/TransportEquation_");///Results
 525         resultFile+=_fileName;
 526
 527         if(_mpi_rank==0){
 528                 _VV.setTime(_time,_nbTimeStep);
 529                 _TT.setTime(_time,_nbTimeStep);
 530                 _Alpha.setTime(_time,_nbTimeStep);
 531                 _Rho.setTime(_time,_nbTimeStep);
 532
 533                 // create mesh and component info
 534                 if (_nbTimeStep ==0 || _restartWithNewFileName){
 535                         if (_restartWithNewFileName)
 536                                 _restartWithNewFileName=false;
 537                         string suppress ="rm -rf "+resultFile+"_*";
 538                         system(suppress.c_str());//Nettoyage des précédents calculs identiques
 539
 540                         switch(_saveFormat)
 541                         {
 542                         case VTK :
 543                                 _VV.writeVTK(resultFile+"Enthalpy");
 544                                 _TT.writeVTK(resultFile+"Temperature");
 545                                 _Alpha.writeVTK(resultFile+"GasFraction");
 546                                 _Rho.writeVTK(resultFile+"MixtureDensity");
 547                                 break;
 548                         case MED :
 549                                 _VV.writeMED(resultFile+"Enthalpy");
 550                                 _TT.writeMED(resultFile+"Temperature");
 551                                 _Alpha.writeMED(resultFile+"GasFraction");
 552                                 _Rho.writeMED(resultFile+"MixtureDensity");
 553                                 break;
 554                         case CSV :
 555                                 _VV.writeCSV(resultFile+"Enthalpy");
 556                                 _TT.writeCSV(resultFile+"Temperature");
 557                                 _Alpha.writeCSV(resultFile+"GasFraction");
 558                                 _Rho.writeCSV(resultFile+"MixtureDensity");
 559                                 break;
 560                         }
 561                 }
 562                 // do not create mesh
 563                 else{
 564                         switch(_saveFormat)
 565                         {
 566                         case VTK :
 567                                 _VV.writeVTK(resultFile+"Enthalpy",false);
 568                                 _TT.writeVTK(resultFile+"Temperature",false);
 569                                 _Alpha.writeVTK(resultFile+"GasFraction",false);
 570                                 _Rho.writeVTK(resultFile+"MixtureDensity",false);
 571                                 break;
 572                         case MED :
 573                                 _VV.writeMED(resultFile+"Enthalpy",false);
 574                                 _TT.writeMED(resultFile+"Temperature",false);
 575                                 _Alpha.writeMED(resultFile+"GasFraction",false);
 576                                 _Rho.writeMED(resultFile+"MixtureDensity",false);
 577                                 break;
 578                         case CSV :
 579                                 _VV.writeCSV(resultFile+"Enthalpy");
 580                                 _TT.writeCSV(resultFile+"Temperature");
 581                                 _Alpha.writeCSV(resultFile+"GasFraction");
 582                                 _Rho.writeCSV(resultFile+"MixtureDensity");
 583                                 break;
 584                         }
 585                 }
 586         }
 587 }
 588
 589 vector<string> TransportEquation::getInputFieldsNames()
 590 {
 591         vector<string> result(2);
 592
 593         result[0]="HeatPower";
 594         result[1]="RodTemperature";
 595
 596         return result;
 597 }
 598
 599 vector<string> TransportEquation::getOutputFieldsNames()
 600 {
 601         vector<string> result(4);
 602
 603         result[0]="Enthalpy";
 604         result[1]="FluidTemperature";
 605         result[2]="VoidFraction";
 606         result[3]="Density";
 607
 608         return result;
 609 }
 610
 611 Field& TransportEquation::getOutputField(const string& nameField )
 612 {
 613         if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
 614                 return getFluidTemperatureField();
 615         else if(nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" )
 616                 return getEnthalpyField();
 617     else        if(nameField=="VoidFraction" || nameField=="VOIDFRACTION" || nameField=="TauxDeVide" || nameField=="TAUXDEVIDE")
 618                 return getVoidFractionField();
 619         else if(nameField=="Density" || nameField=="DENSITY" || nameField=="Densité" || nameField=="DENSITE" )
 620                 return getDensityField();
 621         else
 622     {
 623         cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
 624         throw CdmathException("TransportEquation::getOutputField error : Unknown Field name");
 625     }
 626 }
 627
 628 void
 629 TransportEquation::setInputField(const string& nameField, Field& inputField )
 630 {
 631         if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TemperatureCombustible" || nameField=="TEMPERATURECOMBUSTIBLE")
 632                 return setRodTemperatureField( inputField) ;
 633         else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
 634                 return setHeatPowerField( inputField );
 635         else
 636     {
 637         cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
 638         throw CdmathException("TransportEquation::setInputField error : Unknown Field name");
 639     }
 640 }
 641
 642 void
 643 TransportEquation::setRodTemperatureField(Field rodTemperature){
 644         if(!_initialDataSet)
 645                 throw CdmathException("!!!!!!!! TransportEquation::setRodTemperatureField() set initial field first");
 646
 647         rodTemperature.getMesh().checkFastEquivalWith(_mesh);
 648         _rodTemperatureField=rodTemperature;
 649         _rodTemperatureFieldSet=true;
 650         _isStationary=false;//Source term may be changed after previously reaching a stationary state
 651 }