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