CoreFlows/Models/src/DiffusionEquation.cxx

   1 #include "DiffusionEquation.hxx"
   2 #include "math.h"
   3 #include <algorithm>
   4 #include <fstream>
   5 #include <sstream>
   6
   7 using namespace std;
   8
   9 int DiffusionEquation::fact(int n)
  10 {
  11   return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
  12 }
  13 int DiffusionEquation::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes)
  14 {//assumes Dirichlet node numbering is strictly increasing
  15     int j=0;//indice de parcours des noeuds frontière avec CL Dirichlet
  16     int boundarySize=dirichletNodes.size();
  17     while(j<boundarySize and dirichletNodes[j]<globalIndex)
  18         j++;
  19     if(j==boundarySize)
  20         return globalIndex-boundarySize;
  21     else if (dirichletNodes[j]>globalIndex)
  22         return globalIndex-j;
  23     else
  24         throw CdmathException("DiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
  25 }
  26
  27 int DiffusionEquation::globalNodeIndex(int unknownNodeIndex, std::vector< int > dirichletNodes)
  28 {//assumes Dirichlet boundary node numbering is strictly increasing
  29     int boundarySize=dirichletNodes.size();
  30     /* trivial case where all boundary nodes are Neumann BC */
  31     if(boundarySize==0)
  32         return unknownNodeIndex;
  33
  34     double unknownNodeMax=-1;//max unknown node number in the interval between jth and (j+1)th Dirichlet boundary nodes
  35     int j=0;//indice de parcours des noeuds frontière
  36     //On cherche l'intervale [j,j+1] qui contient le noeud de numéro interieur unknownNodeIndex
  37     while(j+1<boundarySize and unknownNodeMax<unknownNodeIndex)
  38     {
  39         unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
  40         j++;
  41     }
  42
  43     if(j+1==boundarySize)
  44         return unknownNodeIndex+boundarySize;
  45     else //unknownNodeMax>=unknownNodeIndex) hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j]
  46         return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1;
  47 }
  48
  49 Vector DiffusionEquation::gradientNodal(Matrix M, vector< double > values){
  50     vector< Matrix > matrices(_Ndim);
  51
  52     for (int idim=0; idim<_Ndim;idim++){
  53         matrices[idim]=M.deepCopy();
  54         for (int jdim=0; jdim<_Ndim+1;jdim++)
  55                         matrices[idim](jdim,idim) = values[jdim] ;
  56     }
  57
  58         Vector result(_Ndim);
  59     for (int idim=0; idim<_Ndim;idim++)
  60         result[idim] = matrices[idim].determinant();
  61
  62         return result;
  63 }
  64
  65 DiffusionEquation::DiffusionEquation(int dim, bool FECalculation,double rho,double cp, double lambda){
  66     /* Control input value are acceptable */
  67     if(rho<_precision or cp<_precision)
  68     {
  69         std::cout<<"rho="<<rho<<", cp= "<<cp<< ", precision= "<<_precision;
  70         throw CdmathException("Error : parameters rho and cp should be strictly positive");
  71     }
  72     if(lambda < 0.)
  73     {
  74         std::cout<<"conductivity="<<lambda<<endl;
  75         throw CdmathException("Error : conductivity parameter lambda cannot  be negative");
  76     }
  77     if(dim<=0)
  78     {
  79         std::cout<<"space dimension="<<dim<<endl;
  80         throw CdmathException("Error : parameter dim cannot  be negative");
  81     }
  82
  83     cout<<"Diffusion problem with density "<<rho<<", specific heat "<< cp<<", conductivity "<< lambda;
  84     if(FECalculation)
  85         cout<<" and finite elements method"<<endl;
  86     else
  87         cout<<" and finite volumes method"<<endl;
  88
  89     _FECalculation=FECalculation;
  90
  91     /* Finite element data */
  92     _neibMaxNbNodes=0;
  93     _boundaryNodeIds=std::vector< int >(0);
  94     _dirichletNodeIds=std::vector< int >(0);
  95     _NboundaryNodes=0;
  96     _NdirichletNodes=0;
  97     _NunknownNodes=0;
  98
  99     /* Physical parameters */
 100         _conductivity=lambda;
 101         _cp=cp;
 102         _rho=rho;
 103         _diffusivity=_conductivity/(_rho*_cp);
 104         _fluidTemperatureFieldSet=false;
 105         _fluidTemperature=0;
 106
 107     /* Numerical parameters */
 108         _Ndim=dim;
 109         _nVar=1;
 110         _dt_diffusion=0;
 111         _dt_src=0;
 112         _diffusionMatrixSet=false;
 113
 114         _fileName = "CoreFlowsDiffusionProblem";
 115
 116         _runLogFile=new ofstream;
 117
 118     /* Default diffusion tensor is identity matrix */
 119         _DiffusionTensor=Matrix(_Ndim);
 120         for(int idim=0;idim<_Ndim;idim++)
 121                 _DiffusionTensor(idim,idim)=1;
 122 }
 123
 124 void DiffusionEquation::initialize()
 125 {
 126         _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 127
 128         if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim
 129         {
 130         cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<_mesh.getMeshDimension()<<", and space dimension is "<<_mesh.getSpaceDimension()<<endl;
 131                 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<_mesh.getMeshDimension()<<", mesh space dim= "<<_mesh.getSpaceDimension()<<endl;
 132                 *_runLogFile<<"DiffusionEquation::initialize: mesh has incorrect dimension"<<endl;
 133                 _runLogFile->close();
 134                 throw CdmathException("DiffusionEquation::initialize: mesh has incorrect  dimension");
 135         }
 136
 137         if(!_initialDataSet)
 138                 throw CdmathException("DiffusionEquation::initialize() set initial data first");
 139         else
 140         {
 141             cout<<"Initialising the diffusion of a solid temperature using ";
 142             *_runLogFile<<"Initialising the diffusion of a solid temperature using ";
 143             if(!_FECalculation)
 144             {
 145                 cout<< "Finite volumes method"<<endl<<endl;
 146                 *_runLogFile<< "Finite volumes method"<<endl<<endl;
 147             }
 148             else
 149             {
 150                 cout<< "Finite elements method"<<endl<<endl;
 151                 *_runLogFile<< "Finite elements method"<<endl<<endl;
 152             }
 153         }
 154
 155         /**************** Field creation *********************/
 156
 157         if(!_heatPowerFieldSet){
 158         if(_FECalculation){
 159             _heatPowerField=Field("Heat power",NODES,_mesh,1);
 160             for(int i =0; i<_Nnodes; i++)
 161                 _heatPowerField(i) = _heatSource;
 162         }
 163         else{
 164             _heatPowerField=Field("Heat power",CELLS,_mesh,1);
 165             for(int i =0; i<_Nmailles; i++)
 166                 _heatPowerField(i) = _heatSource;
 167         }
 168         _heatPowerFieldSet=true;
 169     }
 170         if(!_fluidTemperatureFieldSet){
 171         if(_FECalculation){
 172             _fluidTemperatureField=Field("Fluid temperature",NODES,_mesh,1);
 173             for(int i =0; i<_Nnodes; i++)
 174                 _fluidTemperatureField(i) = _fluidTemperature;
 175         }
 176         else{
 177             _fluidTemperatureField=Field("Fluid temperature",CELLS,_mesh,1);
 178             for(int i =0; i<_Nmailles; i++)
 179                 _fluidTemperatureField(i) = _fluidTemperature;
 180         }
 181         _fluidTemperatureFieldSet=true;
 182         }
 183
 184     /* Détection des noeuds frontière avec une condition limite de Dirichlet */
 185     if(_FECalculation)
 186     {
 187         if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
 188                         cout<<"1D Finite element method on segments"<<endl;
 189         else if(_Ndim==2)
 190         {
 191                         if( _mesh.isTriangular() )//Mesh dim=2
 192                                 cout<<"2D Finite element method on triangles"<<endl;
 193                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 194                                 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 195                         else
 196                         {
 197                                 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
 198                                 *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
 199                                 _runLogFile->close();
 200                                 throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
 201                         }
 202         }
 203         else if(_Ndim==3)
 204         {
 205                         if( _mesh.isTetrahedral() )//Mesh dim=3
 206                                 cout<<"3D Finite element method on tetrahedra"<<endl;
 207                         else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
 208                                 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 209                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 210                                 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 211                         else
 212                         {
 213                                 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either tetrahedral, either a triangularised surface or 1D network"<<endl;
 214                                 *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
 215                                 _runLogFile->close();
 216                                 throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
 217                         }
 218         }
 219
 220         _boundaryNodeIds = _mesh.getBoundaryNodeIds();
 221         _NboundaryNodes=_boundaryNodeIds.size();
 222
 223         if(_NboundaryNodes==_Nnodes)
 224             cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl;
 225
 226         for(int i=0; i<_NboundaryNodes; i++)
 227             if(_limitField[(_mesh.getNode(_boundaryNodeIds[i])).getGroupName()].bcType==Dirichlet)
 228                 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 229         _NdirichletNodes=_dirichletNodeIds.size();
 230         _NunknownNodes=_Nnodes - _NdirichletNodes;
 231         cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 232         *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 233     }
 234
 235         //creation de la matrice
 236     if(!_FECalculation)
 237         MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
 238     else
 239         MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
 240
 241         VecCreate(PETSC_COMM_SELF, &_Tk);
 242
 243     if(!_FECalculation)
 244         VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
 245     else
 246         VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
 247
 248         VecSetFromOptions(_Tk);
 249         VecDuplicate(_Tk, &_Tn);
 250         VecDuplicate(_Tk, &_Tkm1);
 251         VecDuplicate(_Tk, &_deltaT);
 252         VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
 253         VecDuplicate(_Tk, &_b0);//part of the RHS that comes from the boundary conditions. Computed only once at the first time step
 254
 255     if(!_FECalculation)
 256         for(int i =0; i<_Nmailles;i++)
 257             VecSetValue(_Tn,i,_VV(i), INSERT_VALUES);
 258     else
 259         for(int i =0; i<_Nnodes;i++)
 260             VecSetValue(_Tn,i,_VV(i), INSERT_VALUES);
 261
 262         //Linear solver
 263         KSPCreate(PETSC_COMM_SELF, &_ksp);
 264         KSPSetType(_ksp, _ksptype);
 265         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 266         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 267         KSPGetPC(_ksp, &_pc);
 268         PCSetType(_pc, _pctype);
 269
 270         _initializedMemory=true;
 271         save();//save initial data
 272 }
 273
 274 double DiffusionEquation::computeTimeStep(bool & stop){
 275         if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
 276                 _dt_diffusion=computeDiffusionMatrix(stop);
 277
 278     //reset right hand side
 279         VecCopy(_b0,_b);
 280
 281         _dt_src=computeRHS(stop);
 282
 283         stop=false;
 284         return min(_dt_diffusion,_dt_src);
 285 }
 286
 287 double DiffusionEquation::computeDiffusionMatrix(bool & stop)
 288 {
 289     double result;
 290
 291     if(_FECalculation)
 292         result=computeDiffusionMatrixFE(stop);
 293     else
 294         result=computeDiffusionMatrixFV(stop);
 295
 296     //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
 297     //update value here if variable  heat transfer coefficient
 298     if(_timeScheme == Implicit and _heatTransfertCoeff/(_rho*_cp)>_precision)
 299         MatShift(_A,_heatTransfertCoeff/(_rho*_cp));//Contribution from the liquit/solid heat transfer
 300
 301     if(_verbose or _system)
 302         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 303
 304     return  result;
 305 }
 306
 307 double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){
 308         Cell Cj;
 309         string nameOfGroup;
 310         double dij;//Diffusion coefficients between nodes i and j
 311         MatZeroEntries(_A);
 312         VecZeroEntries(_b);
 313
 314     Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
 315     std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
 316     std::vector< int > nodeIds(_Ndim+1);//cell node Ids
 317     std::vector< Node > nodes(_Ndim+1);//cell nodes
 318     int i_int, j_int; //index of nodes j and k considered as unknown nodes
 319     bool dirichletCell_treated;
 320
 321     std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
 322     for (int idim=0; idim<_Ndim+1;idim++)
 323         values[idim][idim]=1;
 324
 325     /* parameters for boundary treatment */
 326     vector< double > valuesBorder(_Ndim+1);
 327     Vector GradShapeFuncBorder(_Ndim+1);
 328
 329         for (int j=0; j<_Nmailles;j++)
 330     {
 331                 Cj = _mesh.getCell(j);
 332
 333         for (int idim=0; idim<_Ndim+1;idim++){
 334             nodeIds[idim]=Cj.getNodeId(idim);
 335             nodes[idim]=_mesh.getNode(nodeIds[idim]);
 336             for (int jdim=0; jdim<_Ndim;jdim++)
 337                 M(idim,jdim)=nodes[idim].getPoint()[jdim];
 338             M(idim,_Ndim)=1;
 339         }
 340         for (int idim=0; idim<_Ndim+1;idim++)
 341             GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
 342
 343         /* Loop on the edges of the cell */
 344         for (int idim=0; idim<_Ndim+1;idim++)
 345         {
 346             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
 347             {//First node of the edge is not Dirichlet node
 348                 i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 349                 dirichletCell_treated=false;
 350                 for (int jdim=0; jdim<_Ndim+1;jdim++)
 351                 {
 352                     if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
 353                     {//Second node of the edge is not Dirichlet node
 354                         j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 355                         dij=_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure();
 356                         MatSetValue(_A,i_int,j_int,dij, ADD_VALUES);
 357                         if(fabs(dij)>_maxvp)
 358                             _maxvp=fabs(dij);
 359                     }
 360                     else if (!dirichletCell_treated)
 361                     {//Second node of the edge is a Dirichlet node
 362                         dirichletCell_treated=true;
 363                         for (int kdim=0; kdim<_Ndim+1;kdim++)
 364                         {
 365                             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[kdim])!=_dirichletNodeIds.end())
 366                                 valuesBorder[kdim]=_limitField[nameOfGroup].T;
 367                             else
 368                                 valuesBorder[kdim]=0;
 369                         }
 370                         GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
 371                         dij =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 372                         VecSetValue(_b,i_int,dij, ADD_VALUES);
 373                     }
 374                 }
 375             }
 376         }
 377         }
 378
 379     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 380         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 381         VecAssemblyBegin(_b);
 382         VecAssemblyEnd(_b);
 383
 384         _diffusionMatrixSet=true;
 385     stop=false ;
 386
 387         cout<<"_maxvp= "<<_maxvp<< " cfl= "<<_cfl<<" minl= "<<_minl<<endl;
 388         if(fabs(_maxvp)<_precision)
 389                 throw CdmathException("DiffusionEquation::computeDiffusionMatrix(): maximum eigenvalue for time step is zero");
 390         else
 391                 return _cfl/_maxvp;
 392 }
 393
 394 double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){
 395         long nbFaces = _mesh.getNumberOfFaces();
 396         Face Fj;
 397         Cell Cell1,Cell2;
 398         string nameOfGroup;
 399         double inv_dxi, inv_dxj;
 400         double barycenterDistance;
 401         Vector normale(_Ndim);
 402         double dn;
 403         PetscInt idm, idn;
 404         std::vector< int > idCells;
 405         MatZeroEntries(_A);
 406         VecZeroEntries(_b0);
 407
 408         for (int j=0; j<nbFaces;j++){
 409                 Fj = _mesh.getFace(j);
 410
 411                 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 412                 idCells = Fj.getCellsId();
 413                 Cell1 = _mesh.getCell(idCells[0]);
 414                 idm = idCells[0];
 415         for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 416             if (j == Cell1.getFacesId()[l]){
 417                 for (int idim = 0; idim < _Ndim; ++idim)
 418                     normale[idim] = Cell1.getNormalVector(l,idim);
 419                 break;
 420             }
 421         }
 422
 423                 //Compute velocity at the face Fj
 424                 dn=_conductivity*(_DiffusionTensor*normale)*normale;
 425                 if(fabs(dn)>_maxvp)
 426                         _maxvp=fabs(dn);
 427
 428                 // compute 1/dxi = volume of Ci/area of Fj
 429         inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 430
 431                 // If Fj is on the boundary
 432                 if (Fj.getNumberOfCells()==1) {
 433                         if(_verbose )
 434                         {
 435                                 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 436                                 for(int p=0; p<_Ndim; p++)
 437                                         cout << normale[p] << ",";
 438                                 cout << ") "<<endl;
 439                         }
 440                         nameOfGroup = Fj.getGroupName();
 441
 442                         if (_limitField[nameOfGroup].bcType==Neumann){//Nothing to do
 443                         }
 444                         else if(_limitField[nameOfGroup].bcType==Dirichlet){
 445                                 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 446                                 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 447                                 VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 448                         }
 449                         else {
 450                 stop=true ;
 451                                 cout<<"!!!!!!!!!!!!!!!!! Error DiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 452                 cout<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 453                                 cout<<"Accepted boundary conditions are Neumann "<<Neumann<< " and Dirichlet "<<Dirichlet<<endl;
 454                 *_runLogFile<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 455                                 throw CdmathException("Boundary condition not accepted");
 456                         }
 457
 458                         // if Fj is inside the domain
 459                 } else  if (Fj.getNumberOfCells()==2 ){
 460                         if(_verbose )
 461                         {
 462                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 463                                 cout << " ; vecteur normal=(";
 464                                 for(int p=0; p<_Ndim; p++)
 465                                         cout << normale[p] << ",";
 466                                 cout << ") "<<endl;
 467                         }
 468                         Cell2 = _mesh.getCell(idCells[1]);
 469                         idn = idCells[1];
 470                         if (_Ndim > 1)
 471                                 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 472                         else
 473                                 inv_dxj = 1/Cell2.getMeasure();
 474
 475                         barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 476
 477                         MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 478                         MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 479                         MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 480                         MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 481                 }
 482                 else
 483                         throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 484         }
 485
 486         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 487         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 488         VecAssemblyBegin(_b0);
 489         VecAssemblyEnd(_b0);
 490
 491         _diffusionMatrixSet=true;
 492     stop=false ;
 493
 494         cout<<"_maxvp= "<<_maxvp<< " cfl= "<<_cfl<<" minl= "<<_minl<<endl;
 495         if(fabs(_maxvp)<_precision)
 496                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): maximum eigenvalue for time step is zero");
 497         else
 498                 return _cfl*_minl*_minl/_maxvp;
 499 }
 500
 501 double DiffusionEquation::computeRHS(bool & stop){
 502         VecAssemblyBegin(_b);
 503     double Ti;
 504     if(!_FECalculation)
 505         for (int i=0; i<_Nmailles;i++)
 506         {
 507             VecSetValue(_b,i,_heatPowerField(i)/(_rho*_cp),ADD_VALUES);//Contribution of the volumic heat power
 508             //Contribution due to fluid/solide heat exchange
 509             if(_timeScheme == Explicit)
 510             {
 511                 VecGetValues(_Tn, 1, &i, &Ti);
 512                 VecSetValue(_b,i,_heatTransfertCoeff/(_rho*_cp)*(_fluidTemperatureField(i)-Ti),ADD_VALUES);
 513             }
 514             else//Implicit scheme
 515                 VecSetValue(_b,i,_heatTransfertCoeff/(_rho*_cp)* _fluidTemperatureField(i)    ,ADD_VALUES);
 516         }
 517     else
 518         {
 519             Cell Ci;
 520             std::vector< int > nodesId;
 521             for (int i=0; i<_Nmailles;i++)
 522             {
 523                 Ci=_mesh.getCell(i);
 524                 nodesId=Ci.getNodesId();
 525                 for (int j=0; j<nodesId.size();j++)
 526                     if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
 527                     {
 528                         double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
 529                         VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 530                     }
 531             }
 532         }
 533         VecAssemblyEnd(_b);
 534
 535     if(_verbose or _system)
 536         VecView(_b,PETSC_VIEWER_STDOUT_SELF);
 537
 538     stop=false ;
 539     if(_heatTransfertCoeff>_precision)
 540         return _rho*_cp/_heatTransfertCoeff;
 541     else
 542         return INFINITY;
 543 }
 544
 545 bool DiffusionEquation::initTimeStep(double dt){
 546
 547     if(_dt>0 and dt>0)
 548     {
 549         //Remove the contribution from dt to prepare for new initTimeStep. The diffusion matrix is not recomputed
 550         if(_timeScheme == Implicit)
 551             MatShift(_A,-1/_dt+1/dt);
 552         //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
 553     }
 554     else if(dt>0)//_dt==0, first time step
 555     {
 556         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 557         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 558         if(_timeScheme == Implicit)
 559             MatShift(_A,1/_dt);
 560     }
 561     else//dt<=0
 562     {
 563         cout<<"DiffusionEquation::initTimeStep dt= "<<dt<<endl;
 564         throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero");
 565     }
 566     //At this stage _b contains _b0 + power + heat exchange
 567     VecAXPY(_b, 1/_dt, _Tn);
 568
 569         _dt = dt;
 570
 571         if(_verbose && _nbTimeStep%_freqSave ==0)
 572                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 573
 574         return _dt>0;
 575 }
 576
 577 void DiffusionEquation::abortTimeStep(){
 578     //Remove contribution od dt to the RHS
 579         VecAXPY(_b,  -1/_dt, _Tn);
 580         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 581         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 582     //Remove contribution od dt to the matrix
 583         if(_timeScheme == Implicit)
 584                 MatShift(_A,-1/_dt);
 585         _dt = 0;
 586 }
 587
 588 bool DiffusionEquation::iterateTimeStep(bool &converged)
 589 {
 590         bool stop=false;
 591
 592         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 593                 _maxvp=0;
 594                 computeTimeStep(stop);
 595         }
 596         if(stop){
 597                 converged=false;
 598                 return false;
 599         }
 600
 601         if(_timeScheme == Explicit)
 602         {
 603                 MatMult(_A, _Tn, _Tk);
 604                 VecAXPY(_Tk, -1, _b);
 605                 VecScale(_Tk, -_dt);
 606
 607                 converged = true;
 608         }
 609         else
 610         {
 611                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 612                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 613
 614 #if PETSC_VERSION_GREATER_3_5
 615                 KSPSetOperators(_ksp, _A, _A);
 616 #else
 617                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 618 #endif
 619
 620                 if(_conditionNumber)
 621                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 622                 KSPSolve(_ksp, _b, _Tk);
 623
 624                 KSPGetIterationNumber(_ksp, &_PetscIts);
 625                 if( _MaxIterLinearSolver < _PetscIts)
 626                         _MaxIterLinearSolver = _PetscIts;
 627                 if(_PetscIts>=_maxPetscIts)
 628                 {
 629                         cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 630                         *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 631                         converged=false;
 632                         return false;
 633                 }
 634                 else{
 635                         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 636                         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 637                         _erreur_rel= 0;
 638                         double Ti, dTi;
 639
 640             if(!_FECalculation)
 641                 for(int i=0; i<_Nmailles; i++)
 642                 {
 643                     VecGetValues(_deltaT, 1, &i, &dTi);
 644                     VecGetValues(_Tk, 1, &i, &Ti);
 645                     if(_erreur_rel < fabs(dTi/Ti))
 646                         _erreur_rel = fabs(dTi/Ti);
 647                 }
 648             else
 649                 for(int i=0; i<_Nnodes; i++)
 650                 {
 651                     VecGetValues(_deltaT, 1, &i, &dTi);
 652                     VecGetValues(_Tk, 1, &i, &Ti);
 653                     if(_erreur_rel < fabs(dTi/Ti))
 654                         _erreur_rel = fabs(dTi/Ti);
 655                 }
 656                         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 657                 }
 658         }
 659
 660         VecCopy(_Tk, _Tkm1);
 661
 662         return true;
 663 }
 664 void DiffusionEquation::validateTimeStep()
 665 {
 666         VecCopy(_Tk, _deltaT);//ici Tk=Tnp1 donc on a deltaT=Tnp1
 667         VecAXPY(_deltaT,  -1, _Tn);//On obtient deltaT=Tnp1-Tn
 668
 669         _erreur_rel= 0;
 670         double Ti, dTi;
 671
 672     if(!_FECalculation)
 673         for(int i=0; i<_Nmailles; i++)
 674         {
 675             VecGetValues(_deltaT, 1, &i, &dTi);
 676             VecGetValues(_Tk, 1, &i, &Ti);
 677             if(_erreur_rel < fabs(dTi/Ti))
 678                 _erreur_rel = fabs(dTi/Ti);
 679         }
 680     else
 681         for(int i=0; i<_Nnodes; i++)
 682         {
 683             VecGetValues(_deltaT, 1, &i, &dTi);
 684             VecGetValues(_Tk, 1, &i, &Ti);
 685             if(_erreur_rel < fabs(dTi/Ti))
 686                 _erreur_rel = fabs(dTi/Ti);
 687         }
 688
 689         _isStationary =(_erreur_rel <_precision);
 690
 691         VecCopy(_Tk, _Tn);
 692         VecCopy(_Tk, _Tkm1);
 693
 694         if(_verbose && _nbTimeStep%_freqSave ==0)
 695                 cout <<"Valeur propre locale max: " << _maxvp << endl;
 696
 697         _time+=_dt;
 698         _nbTimeStep++;
 699         if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 700         save();
 701 }
 702
 703 void DiffusionEquation::save(){
 704     cout<< "Saving numerical results"<<endl<<endl;
 705     *_runLogFile<< "Saving numerical results"<< endl<<endl;
 706
 707         string resultFile(_path+"/DiffusionEquation");//Results
 708
 709         resultFile+="_";
 710         resultFile+=_fileName;
 711
 712     if(_verbose or _system)
 713         VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
 714
 715     //On remplit le champ
 716     double Ti;
 717     if(!_FECalculation)
 718         for(int i =0; i<_Nmailles;i++)
 719         {
 720             VecGetValues(_Tn, 1, &i, &Ti);
 721             _VV(i)=Ti;
 722         }
 723     else
 724     {
 725         int globalIndex;
 726         for(int i=0; i<_NunknownNodes; i++)
 727         {
 728             VecGetValues(_Tk, 1, &i, &Ti);
 729             globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 730             _VV(globalIndex)=Ti;//Assumes node numbering starts with border nodes
 731         }
 732
 733         Node Ni;
 734         string nameOfGroup;
 735         for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
 736         {
 737             Ni=_mesh.getNode(_dirichletNodeIds[i]);
 738             nameOfGroup = Ni.getGroupName();
 739             _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 740         }
 741     }
 742         _VV.setTime(_time,_nbTimeStep);
 743
 744         // create mesh and component info
 745         if (_nbTimeStep ==0){
 746                 string suppress ="rm -rf "+resultFile+"_*";
 747                 system(suppress.c_str());//Nettoyage des précédents calculs identiques
 748
 749         _VV.setInfoOnComponent(0,"Temperature_(K)");
 750                 switch(_saveFormat)
 751                 {
 752                 case VTK :
 753                         _VV.writeVTK(resultFile);
 754                         break;
 755                 case MED :
 756                         _VV.writeMED(resultFile);
 757                         break;
 758                 case CSV :
 759                         _VV.writeCSV(resultFile);
 760                         break;
 761                 }
 762         }
 763         else{   // do not create mesh
 764                 switch(_saveFormat)
 765                 {
 766                 case VTK :
 767                         _VV.writeVTK(resultFile,false);
 768                         break;
 769                 case MED :
 770                         _VV.writeMED(resultFile,false);
 771                         break;
 772                 case CSV :
 773                         _VV.writeCSV(resultFile);
 774                         break;
 775                 }
 776         }
 777
 778     if(_isStationary)
 779         {
 780         resultFile+="_Stat";
 781         switch(_saveFormat)
 782         {
 783         case VTK :
 784             _VV.writeVTK(resultFile);
 785             break;
 786         case MED :
 787             _VV.writeMED(resultFile);
 788             break;
 789         case CSV :
 790             _VV.writeCSV(resultFile);
 791             break;
 792         }
 793     }
 794 }
 795
 796 void DiffusionEquation::terminate(){
 797         VecDestroy(&_Tn);
 798         VecDestroy(&_Tk);
 799         VecDestroy(&_Tkm1);
 800         VecDestroy(&_deltaT);
 801         VecDestroy(&_b0);
 802         VecDestroy(&_b);
 803         MatDestroy(&_A);
 804 }
 805