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){
   9         if(pEstimate==around1bar300KTransport){
  10                 _Tref=300;
  11                 if(fluid==GasPhase){//Nitrogen pressure 1 bar and temperature 27°C
  12                         _href=3.11e5; //nitrogen enthalpy at 1 bar and 300K
  13                         _cpref=1041;//nitrogen specific heat at constant pressure 1 bar and 300K
  14                         //saturation data for nitrogen at 1 bar and 77K
  15                         _hsatv=0.77e5;//nitrogen vapour enthalpy at saturation at 1 bar
  16                         _hsatl=-1.22e5;//nitrogen liquid enthalpy at saturation at 1 bar
  17                         _rhosatv=4.556;//nitrogen vapour density at saturation at 1 bar
  18                         _rhosatl=806.6;//nitrogen liquid density at saturation at 1 bar
  19                 }
  20                 else{
  21                         //Water at pressure 1 bar and temperature 27°C
  22                         _href=1.127e5; //water enthalpy at 1 bar and 300K
  23                         _cpref=4181;//water specific heat at 1 bar and 300K
  24                         //saturation data for water at 1 bar and 373K
  25                         _hsatv=2.675e6;//Gas enthalpy at saturation at 1 bar
  26                         _hsatl=4.175e5;//water enthalpy at saturation at 1 bar
  27                         _rhosatv=0.6;//Gas density at saturation at 1 bar
  28                         _rhosatl=958;//water density at saturation at 1 bar
  29                 }
  30         }
  31         else{//around155bars600K
  32                 _Tref=618;//=Tsat
  33                 if(fluid==GasPhase){
  34                         _href=2.675e6; //Gas enthalpy at 155 bars and 618K
  35                         _cpref=14001;//Gas specific heat at 155 bar and 618K
  36                 }
  37                 else{//Liquid
  38                         _href=4.175e5;//water enthalpy at 155 bars and 618K
  39                         _cpref=8950;//water specific heat at 155 bar and 618K
  40                 }
  41                 //saturation data for water at 155 bars and 618K
  42                 _hsatv=2.6e6;//Gas enthalpy at saturation at 155 bars
  43                 _hsatl=1.63e6;//water enthalpy at saturation at 155 bars
  44                 _rhosatv=101.9;//Gas density at saturation at 155 bars
  45                 _rhosatl=594.4;//water density at saturation at 155 bars
  46         }
  47         _Ndim=vitesseTransport.size();
  48         _vitesseTransport=Vector(_Ndim);
  49         for(int i=0;i<_Ndim;i++)
  50                 _vitesseTransport[i]=vitesseTransport[i];
  51         _nVar=1;
  52         _dt_transport=0;
  53         _dt_src=0;
  54         _transportMatrixSet=false;
  55 }
  56
  57 void TransportEquation::initialize()
  58 {
  59         if(!_initialDataSet)
  60                 throw CdmathException("TransportEquation::initialize() set initial data first");
  61         else
  62                 cout<<"Initialising the transport of a fluid enthalpy"<<endl;
  63         /**************** Field creation *********************/
  64
  65         //post processing fields used only for saving results
  66         _TT=Field ("Temperature", CELLS, _mesh, 1);
  67         _Alpha=Field ("Void fraction", CELLS, _mesh, 1);
  68         _Rho=Field ("Mixture density", CELLS, _mesh, 1);
  69         //Construction des champs de post-traitement
  70         VecCreate(PETSC_COMM_SELF, &_Hn);
  71         VecSetSizes(_Hn,PETSC_DECIDE,_Nmailles);
  72         VecSetFromOptions(_Hn);
  73         for(int i =0; i<_Nmailles;i++){
  74                 _TT(i)=temperature(_VV(i));
  75                 _Alpha(i)=voidFraction(_VV(i));
  76                 _Rho(i)=density(_Alpha(i));
  77                 VecSetValue(_Hn,i,_VV(i), INSERT_VALUES);
  78         }
  79         if(!_heatPowerFieldSet){
  80                 _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
  81                 for(int i =0; i<_Nmailles; i++)
  82                         _heatPowerField(i) = _heatSource;
  83         }
  84         if(!_rodTemperatureFieldSet){
  85                 _rodTemperatureField=Field("Rod temperature",CELLS,_mesh,1);
  86                 for(int i =0; i<_Nmailles; i++)
  87                         _rodTemperatureField(i) = _rodTemperature;
  88         }
  89
  90         //creation de la matrice
  91         MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
  92         VecDuplicate(_Hn, &_Hk);
  93         VecDuplicate(_Hn, &_Hkm1);
  94         VecDuplicate(_Hn, &_deltaH);
  95         VecDuplicate(_Hn, &_b);//RHS of the linear system: _b=Hn/dt + _b0 + puisance
  96         VecDuplicate(_Hn, &_b0);//part of the RHS that comes from the boundary conditions
  97
  98         //Linear solver
  99         KSPCreate(PETSC_COMM_SELF, &_ksp);
 100         KSPSetType(_ksp, _ksptype);
 101         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 102         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 103         KSPGetPC(_ksp, &_pc);
 104         PCSetType(_pc, _pctype);
 105
 106         _initializedMemory=true;
 107         save();//save initial data
 108 }
 109
 110 double TransportEquation::computeTimeStep(bool & stop){
 111         if(!_transportMatrixSet)
 112                 _dt_transport=computeTransportMatrix();
 113         _dt_src=computeRHS();
 114
 115         stop=false;
 116         return min(_dt_transport,_dt_src);
 117 }
 118 double TransportEquation::computeTransportMatrix(){
 119         long nbFaces = _mesh.getNumberOfFaces();
 120         Face Fj;
 121         Cell Cell1,Cell2;
 122         string nameOfGroup;
 123         double inv_dxi, inv_dxj;
 124         Vector normale(_Ndim);
 125         double un, hk;
 126         PetscInt idm, idn;
 127         std::vector< int > idCells;
 128         MatZeroEntries(_A);
 129         VecZeroEntries(_b0);
 130         for (int j=0; j<nbFaces;j++){
 131                 Fj = _mesh.getFace(j);
 132
 133                 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 134                 idCells = Fj.getCellsId();
 135                 Cell1 = _mesh.getCell(idCells[0]);
 136                 idm = idCells[0];
 137                 if (_Ndim >1){
 138                         for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 139                                 if (j == Cell1.getFacesId()[l]){
 140                                         for (int idim = 0; idim < _Ndim; ++idim)
 141                                                 normale[idim] = Cell1.getNormalVector(l,idim);
 142                                         break;
 143                                 }
 144                         }
 145                 }else{ // _Ndim = 1 : assume that this is normal mesh : the face index increases in positive direction
 146                         if (Fj.getNumberOfCells()<2) {
 147                                 if (j==0)
 148                                         normale[0] = -1;
 149                                 else if (j==nbFaces-1)
 150                                         normale[0] = 1;
 151                                 else
 152                                         throw CdmathException("TransportEquation::ComputeTimeStep(): computation of normal vector failed");
 153                         } else if(Fj.getNumberOfCells()==2){
 154                                 if (idCells[0] < idCells[1])
 155                                         normale[0] = 1;
 156                                 else
 157                                         normale[0] = -1;
 158                         }
 159                 }
 160                 //Compute velocity at the face Fj
 161                 un=normale*_vitesseTransport;
 162                 if(abs(un)>_maxvp)
 163                         _maxvp=abs(un);
 164
 165                 // compute 1/dxi = volume of Ci/area of Fj
 166                 if (_Ndim > 1)
 167                         inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 168                 else
 169                         inv_dxi = 1/Cell1.getMeasure();
 170
 171                 // If Fj is on the boundary
 172                 if (Fj.getNumberOfCells()==1) {
 173                         if(_verbose && _nbTimeStep%_freqSave ==0)
 174                         {
 175                                 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 176                                 for(int p=0; p<_Ndim; p++)
 177                                         cout << normale[p] << ",";
 178                                 cout << ") "<<endl;
 179                         }
 180                         nameOfGroup = Fj.getGroupName();
 181
 182                         if     (_limitField[nameOfGroup].bcType==NeumannTransport || _limitField[nameOfGroup].bcType==OutletTransport ){
 183                                 MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
 184                         }
 185                         else if(_limitField[nameOfGroup].bcType==InletTransport   || _limitField[nameOfGroup].bcType==DirichletTransport){
 186                                 if(un>0){
 187                                         MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
 188                                 }
 189                                 else{
 190                                         hk=_limitField[nameOfGroup].h;
 191                                         VecSetValue(_b0,idm,-inv_dxi*un*hk, ADD_VALUES);
 192                                 }
 193                         }
 194                         else {//bcType=NoneBCTransport
 195                                 cout<<"!!!!!!!!!!!!!!! Error TransportEquation::computeTransportMatrix() !!!!!!!!!!"<<endl;
 196                                 cout<<"!!!!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<" !!!!!!!!!!!!!! "<<endl;
 197                                 cout<<"Accepted boundary conditions are NeumannTransport "<<NeumannTransport<< " and InletTransport "<< InletTransport <<endl;
 198                                 throw CdmathException("Boundary condition not accepted");
 199                         }
 200                         // if Fj is inside the domain
 201                 } else  if (Fj.getNumberOfCells()==2 ){
 202                         if(_verbose && _nbTimeStep%_freqSave ==0)
 203                         {
 204                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 205                                 cout << " ; vecteur normal=(";
 206                                 for(int p=0; p<_Ndim; p++)
 207                                         cout << normale[p] << ",";
 208                                 cout << "). "<<endl;
 209                         }
 210                         Cell2 = _mesh.getCell(idCells[1]);
 211                         idn = idCells[1];
 212                         if (_Ndim > 1)
 213                                 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 214                         else
 215                                 inv_dxj = 1/Cell2.getMeasure();
 216
 217                         if(un>0){
 218                                 MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
 219                                 MatSetValue(_A,idn,idm,-inv_dxj*un, ADD_VALUES);
 220                         }
 221                         else{
 222                                 MatSetValue(_A,idm,idn,inv_dxi*un, ADD_VALUES);
 223                                 MatSetValue(_A,idn,idn,-inv_dxj*un, ADD_VALUES);
 224                         }
 225                 }
 226                 else
 227                         throw CdmathException("TransportEquation::ComputeTimeStep(): incompatible number of cells around a face");
 228         }
 229         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 230         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 231         VecAssemblyBegin(_b0);
 232         VecAssemblyEnd(_b0);
 233         _transportMatrixSet=true;
 234         if(abs(_maxvp)<_precision)
 235                 throw CdmathException("TransportEquation::computeTransportMatrix(): maximum eigenvalue for time step is zero");
 236         else
 237                 return _cfl*_minl/_maxvp;
 238 }
 239 double TransportEquation::computeRHS(){
 240         double rhomin=INFINITY;
 241         VecCopy(_b0,_b);
 242         VecAssemblyBegin(_b);
 243         if(_system)
 244                 cout<<"second membre of transport problem"<<endl;
 245         for (int i=0; i<_Nmailles;i++){
 246                 VecSetValue(_b,i,_heatTransfertCoeff*(_rodTemperatureField(i)-_TT(i))/_Rho(i),ADD_VALUES);
 247                 VecSetValue(_b,i,_heatPowerField(i)/_Rho(i),ADD_VALUES);
 248                 if(_system)
 249                         cout<<_heatPowerField(i)/_Rho(i)<<endl;
 250                 if(_Rho(i)<rhomin)
 251                         rhomin=_Rho(i);
 252         }
 253         VecAssemblyEnd(_b);
 254         if(_system)
 255                 VecView(_b,  PETSC_VIEWER_STDOUT_WORLD);
 256
 257         return rhomin*_cpref/_heatTransfertCoeff;
 258 }
 259 void TransportEquation::updatePrimitives()
 260 {
 261         double hi;
 262         for(int i=0; i<_Nmailles; i++)
 263         {
 264                 VecGetValues(_Hk, 1, &i, &hi);
 265                 _VV(i)=hi;
 266                 _TT(i)=temperature(hi);
 267                 _Alpha(i)=voidFraction(hi);
 268                 _Rho(i)=density(_Alpha(i));
 269         }
 270 }
 271
 272 bool TransportEquation::initTimeStep(double dt){
 273         _dt = dt;
 274         if(_verbose && _nbTimeStep%_freqSave ==0)
 275                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 276         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 277         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 278
 279         if(_timeScheme == Implicit)
 280                 MatShift(_A,1/_dt);
 281
 282         return _dt>0;
 283 }
 284
 285 void TransportEquation::abortTimeStep(){
 286         VecAXPY(_b,  -1/_dt, _Hn);
 287         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 288         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 289
 290         if(_timeScheme == Implicit)
 291                 MatShift(_A,-1/_dt);
 292         _dt = 0;
 293 }
 294
 295 bool TransportEquation::iterateTimeStep(bool &converged)
 296 {
 297         bool stop=false;
 298
 299         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 300                 _maxvp=0;
 301                 computeTimeStep(stop);
 302         }
 303         if(stop){
 304                 converged=false;
 305                 return false;
 306         }
 307         VecAXPY(_b, 1/_dt, _Hn);
 308         if(_system)
 309         {
 310                 cout << "Vecteur Hn: " << endl;
 311                 VecView(_Hn,  PETSC_VIEWER_STDOUT_WORLD);
 312                 cout << endl;
 313                 cout<<"Vecteur _b "<<endl;
 314                 VecView(_b,  PETSC_VIEWER_STDOUT_SELF);
 315                 cout << "Matrice A "<<endl;
 316                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 317         }
 318
 319         if(_timeScheme == Explicit)
 320         {
 321                 MatMult(_A, _Hn, _Hk);
 322                 if(_system)
 323                 {
 324                         cout << "Nouveau vecteur Hk: " << endl;
 325                         VecView(_Hk,  PETSC_VIEWER_STDOUT_WORLD);
 326                         cout << endl;
 327                 }
 328                 VecAXPY(_Hk, -1, _b);
 329                 if(_system)
 330                 {
 331                         cout << "Nouveau vecteur Hk-b: " << endl;
 332                         VecView(_Hk,  PETSC_VIEWER_STDOUT_WORLD);
 333                         cout << endl;
 334                 }
 335                 VecScale(_Hk, -_dt);
 336                 if(_system)
 337                 {
 338                         cout << "Nouveau vecteur dt*(Hk-b): " << endl;
 339                         VecView(_Hk,  PETSC_VIEWER_STDOUT_WORLD);
 340                         cout << endl;
 341                 }
 342                 converged = true;
 343         }
 344         else
 345         {
 346                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 347                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 348
 349 #if PETSC_VERSION_GREATER_3_5
 350                 KSPSetOperators(_ksp, _A, _A);
 351 #else
 352                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 353 #endif
 354
 355                 if(_conditionNumber)
 356                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 357
 358                 KSPSolve(_ksp, _b, _Hk);
 359                 MatShift(_A,-1/_dt);
 360
 361                 KSPGetIterationNumber(_ksp, &_PetscIts);
 362                 if( _MaxIterLinearSolver < _PetscIts)
 363                         _MaxIterLinearSolver = _PetscIts;
 364                 if(_PetscIts>=_maxPetscIts)
 365                 {
 366                         cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 367                         converged=false;
 368                         return false;
 369                 }
 370                 else
 371                 {
 372                         VecCopy(_Hk, _deltaH);//ici on a deltaH=Hk
 373                         VecAXPY(_deltaH,  -1, _Hkm1);//On obtient deltaH=Hk-Hkm1
 374                         _erreur_rel= 0;
 375                         double hi, dhi;
 376
 377                         for(int i=0; i<_Nmailles; i++)
 378                         {
 379                                 VecGetValues(_deltaH, 1, &i, &dhi);
 380                                 VecGetValues(_Hk, 1, &i, &hi);
 381                                 if(_erreur_rel < fabs(dhi/hi))
 382                                         _erreur_rel = fabs(dhi/hi);
 383                         }
 384                 }
 385
 386                 converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 387         }
 388
 389         updatePrimitives();
 390
 391         VecCopy(_Hk, _Hkm1);
 392
 393
 394         return true;
 395 }
 396 void TransportEquation::validateTimeStep()
 397 {
 398         VecCopy(_Hk, _deltaH);//ici Hk=Hnp1 donc on a deltaH=Hnp1
 399         VecAXPY(_deltaH,  -1, _Hn);//On obtient deltaH=Hnp1-Hn
 400
 401         _erreur_rel= 0;
 402         double hi, dhi;
 403
 404         for(int i=0; i<_Nmailles; i++)
 405         {
 406                 VecGetValues(_deltaH, 1, &i, &dhi);
 407                 VecGetValues(_Hk, 1, &i, &hi);
 408                 if(_erreur_rel < fabs(dhi/hi))
 409                         _erreur_rel = fabs(dhi/hi);
 410         }
 411         _isStationary =(_erreur_rel <_precision);
 412
 413         VecCopy(_Hk, _Hn);
 414         VecCopy(_Hk, _Hkm1);
 415
 416         if(_nbTimeStep%_freqSave ==0)
 417         {
 418                 cout <<"Valeur propre locale max: " << _maxvp << endl;
 419                 //Find minimum and maximum void fractions
 420                 double alphamin=INFINITY;
 421                 double alphamax=-INFINITY;
 422                 for(int i=0; i<_Nmailles; i++)
 423                 {
 424                         if(_Alpha(i)>alphamax)
 425                                 alphamax=_Alpha(i);
 426                         if(_Alpha(i)<alphamin)
 427                                 alphamin=_Alpha(i);
 428                 }
 429                 cout<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<endl;
 430         }
 431         _time+=_dt;
 432         _nbTimeStep++;
 433         save();
 434
 435 }
 436
 437 void TransportEquation::terminate(){
 438         VecDestroy(&_Hn);
 439         VecDestroy(&_Hk);
 440         VecDestroy(&_Hkm1);
 441         VecDestroy(&_deltaH);
 442         VecDestroy(&_b0);
 443         VecDestroy(&_b);
 444         MatDestroy(&_A);
 445 }
 446
 447 void TransportEquation::save(){
 448         string resultFile(_path+"/TransportEquation_");///Results
 449         resultFile+=_fileName;
 450
 451         _VV.setTime(_time,_nbTimeStep);
 452         _TT.setTime(_time,_nbTimeStep);
 453         _Alpha.setTime(_time,_nbTimeStep);
 454         _Rho.setTime(_time,_nbTimeStep);
 455
 456         // create mesh and component info
 457         if (_nbTimeStep ==0 || _restartWithNewFileName){
 458                 if (_restartWithNewFileName)
 459                         _restartWithNewFileName=false;
 460                 string suppress ="rm -rf "+resultFile+"_*";
 461                 system(suppress.c_str());//Nettoyage des précédents calculs identiques
 462
 463                 switch(_saveFormat)
 464                 {
 465                 case VTK :
 466                         _VV.writeVTK(resultFile+"Enthalpy");
 467                         _TT.writeVTK(resultFile+"Temperature");
 468                         _Alpha.writeVTK(resultFile+"GasFraction");
 469                         _Rho.writeVTK(resultFile+"MixtureDensity");
 470                         break;
 471                 case MED :
 472                         _VV.writeMED(resultFile+"Enthalpy");
 473                         _TT.writeMED(resultFile+"Temperature");
 474                         _Alpha.writeMED(resultFile+"GasFraction");
 475                         _Rho.writeMED(resultFile+"MixtureDensity");
 476                         break;
 477                 case CSV :
 478                         _VV.writeCSV(resultFile+"Enthalpy");
 479                         _TT.writeCSV(resultFile+"Temperature");
 480                         _Alpha.writeCSV(resultFile+"GasFraction");
 481                         _Rho.writeCSV(resultFile+"MixtureDensity");
 482                         break;
 483                 }
 484         }
 485         // do not create mesh
 486         else{
 487                 switch(_saveFormat)
 488                 {
 489                 case VTK :
 490                         _VV.writeVTK(resultFile+"Enthalpy",false);
 491                         _TT.writeVTK(resultFile+"Temperature",false);
 492                         _Alpha.writeVTK(resultFile+"GasFraction",false);
 493                         _Rho.writeVTK(resultFile+"MixtureDensity",false);
 494                         break;
 495                 case MED :
 496                         _VV.writeMED(resultFile+"Enthalpy",false);
 497                         _TT.writeMED(resultFile+"Temperature",false);
 498                         _Alpha.writeMED(resultFile+"GasFraction",false);
 499                         _Rho.writeMED(resultFile+"MixtureDensity",false);
 500                         break;
 501                 case CSV :
 502                         _VV.writeCSV(resultFile+"Enthalpy");
 503                         _TT.writeCSV(resultFile+"Temperature");
 504                         _Alpha.writeCSV(resultFile+"GasFraction");
 505                         _Rho.writeCSV(resultFile+"MixtureDensity");
 506                         break;
 507                 }
 508         }
 509 }
 510
 511 vector<string> TransportEquation::getOutputFieldsNames()
 512 {
 513         vector<string> result(2);
 514
 515         result[0]="Enthalpy";
 516         result[1]="FluidTemperature";
 517
 518         return result;
 519 }
 520
 521 Field& TransportEquation::getOutputField(const string& nameField )
 522 {
 523         if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" )
 524                 return getFluidTemperatureField();
 525         else if(nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPY" )
 526                 return getEnthalpyField();
 527     else
 528     {
 529         cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
 530         throw CdmathException("TransportEquation::getOutputField error : Unknown Field name");
 531     }
 532 }
 533