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==DirichletDiffusion)
 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==NeumannDiffusion){
 443                 VecSetValue(_b,idm,   -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
 444                         }
 445                         else if(_limitField[nameOfGroup].bcType==DirichletDiffusion){
 446                                 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 447                                 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 448                                 VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 449                         }
 450                         else {
 451                 stop=true ;
 452                                 cout<<"!!!!!!!!!!!!!!!!! Error DiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 453                 cout<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 454                                 cout<<"Accepted boundary conditions are NeumannDiffusion "<<NeumannDiffusion<< " and DirichletDiffusion "<<DirichletDiffusion<<endl;
 455                 *_runLogFile<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 456                                 throw CdmathException("Boundary condition not accepted");
 457                         }
 458
 459                         // if Fj is inside the domain
 460                 } else  if (Fj.getNumberOfCells()==2 ){
 461                         if(_verbose )
 462                         {
 463                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 464                                 cout << " ; vecteur normal=(";
 465                                 for(int p=0; p<_Ndim; p++)
 466                                         cout << normale[p] << ",";
 467                                 cout << ") "<<endl;
 468                         }
 469                         Cell2 = _mesh.getCell(idCells[1]);
 470                         idn = idCells[1];
 471                         if (_Ndim > 1)
 472                                 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 473                         else
 474                                 inv_dxj = 1/Cell2.getMeasure();
 475
 476                         barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 477
 478                         MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 479                         MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 480                         MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 481                         MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 482                 }
 483                 else
 484                         throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 485         }
 486
 487         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 488         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 489         VecAssemblyBegin(_b0);
 490         VecAssemblyEnd(_b0);
 491
 492         _diffusionMatrixSet=true;
 493     stop=false ;
 494
 495         cout<<"_maxvp= "<<_maxvp<< " cfl= "<<_cfl<<" minl= "<<_minl<<endl;
 496         if(fabs(_maxvp)<_precision)
 497                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): maximum eigenvalue for time step is zero");
 498         else
 499                 return _cfl*_minl*_minl/_maxvp;
 500 }
 501
 502 double DiffusionEquation::computeRHS(bool & stop){
 503         VecAssemblyBegin(_b);
 504     double Ti;
 505     if(!_FECalculation)
 506         for (int i=0; i<_Nmailles;i++)
 507         {
 508             VecSetValue(_b,i,_heatPowerField(i)/(_rho*_cp),ADD_VALUES);//Contribution of the volumic heat power
 509             //Contribution due to fluid/solide heat exchange
 510             if(_timeScheme == Explicit)
 511             {
 512                 VecGetValues(_Tn, 1, &i, &Ti);
 513                 VecSetValue(_b,i,_heatTransfertCoeff/(_rho*_cp)*(_fluidTemperatureField(i)-Ti),ADD_VALUES);
 514             }
 515             else//Implicit scheme
 516                 VecSetValue(_b,i,_heatTransfertCoeff/(_rho*_cp)* _fluidTemperatureField(i)    ,ADD_VALUES);
 517         }
 518     else
 519         {
 520             Cell Ci;
 521             std::vector< int > nodesId;
 522             for (int i=0; i<_Nmailles;i++)
 523             {
 524                 Ci=_mesh.getCell(i);
 525                 nodesId=Ci.getNodesId();
 526                 for (int j=0; j<nodesId.size();j++)
 527                     if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
 528                     {
 529                         double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
 530                         VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 531                     }
 532             }
 533         }
 534         VecAssemblyEnd(_b);
 535
 536     if(_verbose or _system)
 537         VecView(_b,PETSC_VIEWER_STDOUT_SELF);
 538
 539     stop=false ;
 540     if(_heatTransfertCoeff>_precision)
 541         return _rho*_cp/_heatTransfertCoeff;
 542     else
 543         return INFINITY;
 544 }
 545
 546 bool DiffusionEquation::initTimeStep(double dt){
 547
 548     if(_dt>0 and dt>0)
 549     {
 550         //Remove the contribution from dt to prepare for new initTimeStep. The diffusion matrix is not recomputed
 551         if(_timeScheme == Implicit)
 552             MatShift(_A,-1/_dt+1/dt);
 553         //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
 554     }
 555     else if(dt>0)//_dt==0, first time step
 556     {
 557         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 558         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 559         if(_timeScheme == Implicit)
 560             MatShift(_A,1/_dt);
 561     }
 562     else//dt<=0
 563     {
 564         cout<<"DiffusionEquation::initTimeStep dt= "<<dt<<endl;
 565         throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero");
 566     }
 567     //At this stage _b contains _b0 + power + heat exchange
 568     VecAXPY(_b, 1/_dt, _Tn);
 569
 570         _dt = dt;
 571
 572         if(_verbose && _nbTimeStep%_freqSave ==0)
 573                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 574
 575         return _dt>0;
 576 }
 577
 578 void DiffusionEquation::abortTimeStep(){
 579     //Remove contribution od dt to the RHS
 580         VecAXPY(_b,  -1/_dt, _Tn);
 581         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 582         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 583     //Remove contribution od dt to the matrix
 584         if(_timeScheme == Implicit)
 585                 MatShift(_A,-1/_dt);
 586         _dt = 0;
 587 }
 588
 589 bool DiffusionEquation::iterateTimeStep(bool &converged)
 590 {
 591         bool stop=false;
 592
 593         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 594                 _maxvp=0;
 595                 computeTimeStep(stop);
 596         }
 597         if(stop){
 598                 converged=false;
 599                 return false;
 600         }
 601
 602         if(_timeScheme == Explicit)
 603         {
 604                 MatMult(_A, _Tn, _Tk);
 605                 VecAXPY(_Tk, -1, _b);
 606                 VecScale(_Tk, -_dt);
 607
 608                 converged = true;
 609         }
 610         else
 611         {
 612                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 613                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 614
 615 #if PETSC_VERSION_GREATER_3_5
 616                 KSPSetOperators(_ksp, _A, _A);
 617 #else
 618                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 619 #endif
 620
 621                 if(_conditionNumber)
 622                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 623                 KSPSolve(_ksp, _b, _Tk);
 624
 625                 KSPGetIterationNumber(_ksp, &_PetscIts);
 626                 if( _MaxIterLinearSolver < _PetscIts)
 627                         _MaxIterLinearSolver = _PetscIts;
 628                 if(_PetscIts>=_maxPetscIts)
 629                 {
 630                         cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 631                         *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 632                         converged=false;
 633                         return false;
 634                 }
 635                 else{
 636                         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 637                         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 638                         _erreur_rel= 0;
 639                         double Ti, dTi;
 640
 641             if(!_FECalculation)
 642                 for(int i=0; i<_Nmailles; i++)
 643                 {
 644                     VecGetValues(_deltaT, 1, &i, &dTi);
 645                     VecGetValues(_Tk, 1, &i, &Ti);
 646                     if(_erreur_rel < fabs(dTi/Ti))
 647                         _erreur_rel = fabs(dTi/Ti);
 648                 }
 649             else
 650                 for(int i=0; i<_Nnodes; i++)
 651                 {
 652                     VecGetValues(_deltaT, 1, &i, &dTi);
 653                     VecGetValues(_Tk, 1, &i, &Ti);
 654                     if(_erreur_rel < fabs(dTi/Ti))
 655                         _erreur_rel = fabs(dTi/Ti);
 656                 }
 657                         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 658                 }
 659         }
 660
 661         VecCopy(_Tk, _Tkm1);
 662
 663         return true;
 664 }
 665 void DiffusionEquation::validateTimeStep()
 666 {
 667         VecCopy(_Tk, _deltaT);//ici Tk=Tnp1 donc on a deltaT=Tnp1
 668         VecAXPY(_deltaT,  -1, _Tn);//On obtient deltaT=Tnp1-Tn
 669
 670         _erreur_rel= 0;
 671         double Ti, dTi;
 672
 673     if(!_FECalculation)
 674         for(int i=0; i<_Nmailles; i++)
 675         {
 676             VecGetValues(_deltaT, 1, &i, &dTi);
 677             VecGetValues(_Tk, 1, &i, &Ti);
 678             if(_erreur_rel < fabs(dTi/Ti))
 679                 _erreur_rel = fabs(dTi/Ti);
 680         }
 681     else
 682         for(int i=0; i<_Nnodes; i++)
 683         {
 684             VecGetValues(_deltaT, 1, &i, &dTi);
 685             VecGetValues(_Tk, 1, &i, &Ti);
 686             if(_erreur_rel < fabs(dTi/Ti))
 687                 _erreur_rel = fabs(dTi/Ti);
 688         }
 689
 690         _isStationary =(_erreur_rel <_precision);
 691
 692         VecCopy(_Tk, _Tn);
 693         VecCopy(_Tk, _Tkm1);
 694
 695         if(_verbose && _nbTimeStep%_freqSave ==0)
 696                 cout <<"Valeur propre locale max: " << _maxvp << endl;
 697
 698         _time+=_dt;
 699         _nbTimeStep++;
 700         if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 701         save();
 702 }
 703
 704 void DiffusionEquation::save(){
 705     cout<< "Saving numerical results"<<endl<<endl;
 706     *_runLogFile<< "Saving numerical results"<< endl<<endl;
 707
 708         string resultFile(_path+"/DiffusionEquation");//Results
 709
 710         resultFile+="_";
 711         resultFile+=_fileName;
 712
 713     if(_verbose or _system)
 714         VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
 715
 716     //On remplit le champ
 717     double Ti;
 718     if(!_FECalculation)
 719         for(int i =0; i<_Nmailles;i++)
 720         {
 721             VecGetValues(_Tn, 1, &i, &Ti);
 722             _VV(i)=Ti;
 723         }
 724     else
 725     {
 726         int globalIndex;
 727         for(int i=0; i<_NunknownNodes; i++)
 728         {
 729             VecGetValues(_Tk, 1, &i, &Ti);
 730             globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 731             _VV(globalIndex)=Ti;//Assumes node numbering starts with border nodes
 732         }
 733
 734         Node Ni;
 735         string nameOfGroup;
 736         for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
 737         {
 738             Ni=_mesh.getNode(_dirichletNodeIds[i]);
 739             nameOfGroup = Ni.getGroupName();
 740             _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 741         }
 742     }
 743         _VV.setTime(_time,_nbTimeStep);
 744
 745         // create mesh and component info
 746         if (_nbTimeStep ==0 || _restartWithNewFileName){
 747                 if (_restartWithNewFileName)
 748                         _restartWithNewFileName=false;
 749                 string suppress ="rm -rf "+resultFile+"_*";
 750                 system(suppress.c_str());//Nettoyage des précédents calculs identiques
 751
 752         _VV.setInfoOnComponent(0,"Temperature_(K)");
 753                 switch(_saveFormat)
 754                 {
 755                 case VTK :
 756                         _VV.writeVTK(resultFile);
 757                         break;
 758                 case MED :
 759                         _VV.writeMED(resultFile);
 760                         break;
 761                 case CSV :
 762                         _VV.writeCSV(resultFile);
 763                         break;
 764                 }
 765         }
 766         else{   // do not create mesh
 767                 switch(_saveFormat)
 768                 {
 769                 case VTK :
 770                         _VV.writeVTK(resultFile,false);
 771                         break;
 772                 case MED :
 773                         _VV.writeMED(resultFile,false);
 774                         break;
 775                 case CSV :
 776                         _VV.writeCSV(resultFile);
 777                         break;
 778                 }
 779         }
 780
 781     if(_isStationary)
 782         {
 783         resultFile+="_Stat";
 784         switch(_saveFormat)
 785         {
 786         case VTK :
 787             _VV.writeVTK(resultFile);
 788             break;
 789         case MED :
 790             _VV.writeMED(resultFile);
 791             break;
 792         case CSV :
 793             _VV.writeCSV(resultFile);
 794             break;
 795         }
 796     }
 797 }
 798
 799 void DiffusionEquation::terminate(){
 800         VecDestroy(&_Tn);
 801         VecDestroy(&_Tk);
 802         VecDestroy(&_Tkm1);
 803         VecDestroy(&_deltaT);
 804         VecDestroy(&_b0);
 805         VecDestroy(&_b);
 806         MatDestroy(&_A);
 807 }
 808
 809 vector<string> DiffusionEquation::getOutputFieldsNames()
 810 {
 811         vector<string> result(2);
 812
 813         result[0]="FluidTemperature";
 814         result[1]="RodTemperature";
 815
 816         return result;
 817 }
 818
 819 Field& DiffusionEquation::getOutputField(const string& nameField )
 820 {
 821         if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE" )
 822                 return getFluidTemperatureField();
 823         else if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
 824                 return getRodTemperatureField();
 825     else
 826     {
 827         cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
 828         throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
 829     }
 830 }
 831