CoreFlows/Models/src/StationaryDiffusionEquation.cxx

   1 #include "StationaryDiffusionEquation.hxx"
   2 #include "math.h"
   3 #include <algorithm>
   4 #include <fstream>
   5 #include <sstream>
   6
   7 using namespace std;
   8
   9 int StationaryDiffusionEquation::fact(int n)
  10 {
  11   return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
  12 }
  13 int StationaryDiffusionEquation::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("StationaryDiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
  25 }
  26
  27 int StationaryDiffusionEquation::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 StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda){
  50         PetscBool petscInitialized;
  51         PetscInitialized(&petscInitialized);
  52         if(!petscInitialized)
  53                 PetscInitialize(NULL,NULL,0,0);
  54
  55     if(lambda < 0.)
  56     {
  57         std::cout<<"conductivity="<<lambda<<endl;
  58         throw CdmathException("Error : conductivity parameter lambda cannot  be negative");
  59     }
  60     if(dim<=0)
  61     {
  62         std::cout<<"space dimension="<<dim<<endl;
  63         throw CdmathException("Error : parameter dim cannot  be negative");
  64     }
  65
  66     _FECalculation=FECalculation;
  67     _onlyNeumannBC=false;
  68
  69         _Ndim=dim;
  70         _nVar=1;//scalar prolem
  71         _dt_src=0;
  72         _diffusionMatrixSet=false;
  73         _initializedMemory=false;
  74
  75     //Mesh data
  76     _neibMaxNbCells=0;
  77     _meshSet=false;
  78     _neibMaxNbNodes=0;
  79
  80     //Boundary conditions
  81     _boundaryNodeIds=std::vector< int >(0);
  82     _dirichletNodeIds=std::vector< int >(0);
  83     _NboundaryNodes=0;
  84     _NdirichletNodes=0;
  85     _NunknownNodes=0;
  86     _dirichletValuesSet=false;
  87
  88     //Linear solver data
  89         _precision=1.e-6;
  90         _precision_Newton=_precision;
  91         _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
  92         _maxPetscIts=50;
  93         _maxNewtonIts=50;
  94         _NEWTON_its=0;
  95         int _PetscIts=0;//the number of iterations of the linear solver
  96         _ksptype = (char*)&KSPGMRES;
  97         _pctype = (char*)&PCLU;
  98         _conditionNumber=false;
  99         _erreur_rel= 0;
 100
 101     //parameters for monitoring simulation
 102         _verbose = false;
 103         _system = false;
 104         _runLogFile=new ofstream;
 105
 106         //result save parameters
 107         _fileName = "StationaryDiffusionProblem";
 108         char result[ PATH_MAX ];//extracting current directory
 109         getcwd(result, PATH_MAX );
 110         _path=string( result );
 111         _saveFormat=VTK;
 112     _computationCompletedSuccessfully=false;
 113
 114     //heat transfer parameters
 115         _conductivity=lambda;
 116         _fluidTemperatureFieldSet=false;
 117         _fluidTemperature=0;
 118         _heatPowerFieldSet=false;
 119         _heatTransfertCoeff=0;
 120         _heatSource=0;
 121 }
 122
 123 void StationaryDiffusionEquation::initialize()
 124 {
 125         _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 126
 127         if(!_meshSet)
 128                 throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first");
 129         else
 130     {
 131                 cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
 132         *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
 133         if(!_FECalculation)
 134         {
 135             cout<< "Finite volumes method"<<endl<<endl;
 136             *_runLogFile<< "Finite volumes method"<<endl<<endl;
 137         }
 138         else
 139         {
 140             cout<< "Finite elements method"<<endl<<endl;
 141             *_runLogFile<< "Finite elements method"<<endl<<endl;
 142         }
 143     }
 144
 145         _DiffusionTensor=Matrix(_Ndim);
 146         for(int idim=0;idim<_Ndim;idim++)
 147                 _DiffusionTensor(idim,idim)=1;
 148         /**************** Field creation *********************/
 149
 150         if(!_heatPowerFieldSet){
 151         if(_FECalculation){
 152             _heatPowerField=Field("Heat power",NODES,_mesh,1);
 153             for(int i =0; i<_Nnodes; i++)
 154                 _heatPowerField(i) = _heatSource;
 155         }
 156         else{
 157             _heatPowerField=Field("Heat power",CELLS,_mesh,1);
 158             for(int i =0; i<_Nmailles; i++)
 159                 _heatPowerField(i) = _heatSource;
 160         }
 161         _heatPowerFieldSet=true;
 162     }
 163         if(!_fluidTemperatureFieldSet){
 164         if(_FECalculation){
 165             _fluidTemperatureField=Field("Fluid temperature",NODES,_mesh,1);
 166             for(int i =0; i<_Nnodes; i++)
 167                 _fluidTemperatureField(i) = _fluidTemperature;
 168         }
 169         else{
 170             _fluidTemperatureField=Field("Fluid temperature",CELLS,_mesh,1);
 171             for(int i =0; i<_Nmailles; i++)
 172                 _fluidTemperatureField(i) = _fluidTemperature;
 173         }
 174         _fluidTemperatureFieldSet=true;
 175         }
 176
 177     /* Détection des noeuds frontière avec une condition limite de Dirichlet */
 178     if(_FECalculation)
 179     {
 180         /*
 181         vector<int> _boundaryFaceIds = _mesh.getBoundaryFaceIds();
 182
 183         cout <<"Total number of faces " <<_mesh.getNumberOfFaces()<<", Number of boundary faces " << _boundaryFaceIds.size()<<endl;
 184         for(int i=0; i<_boundaryFaceIds.size(); i++)
 185             cout<<", "<<_boundaryFaceIds[i];
 186         cout<<endl;
 187
 188         cout <<"Total number of nodes " <<_mesh.getNumberOfNodes()<<", Number of boundary nodes " << _NboundaryNodes<<endl;
 189         for(int i=0; i<_NboundaryNodes; i++)
 190             cout<<", "<<_boundaryNodeIds[i];
 191         cout<<endl;
 192         */
 193         if(_NboundaryNodes==_Nnodes)
 194             cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
 195
 196         for(int i=0; i<_NboundaryNodes; i++)
 197         {
 198             std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
 199             if( it != _dirichletBoundaryValues.end() )
 200                 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 201             else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
 202             {
 203                 cout<<"!!! No boundary value set for boundary node" << _boundaryNodeIds[i]<< endl;
 204                 *_runLogFile<< "!!! No boundary value set for boundary node" << _boundaryNodeIds[i]<<endl;
 205                 _runLogFile->close();
 206                 throw CdmathException("Missing boundary value");
 207             }
 208             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoTypeSpecified)
 209             {
 210                 cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
 211                 *_runLogFile<< "!!!No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
 212                 _runLogFile->close();
 213                 throw CdmathException("Missing boundary condition");
 214             }
 215             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==Dirichlet)
 216                 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 217             else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=Neumann)
 218             {
 219                 cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
 220                 cout<<"!!! Accepted boundary conditions are Dirichlet "<< Dirichlet <<" and Neumann "<< Neumann << endl;
 221                 *_runLogFile<< "Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
 222                 *_runLogFile<< "Accepted boundary conditions are Dirichlet "<< Dirichlet <<" and Neumann "<< Neumann <<endl;
 223                 _runLogFile->close();
 224                 throw CdmathException("Wrong boundary condition");
 225             }
 226         }
 227         _NdirichletNodes=_dirichletNodeIds.size();
 228         _NunknownNodes=_Nnodes - _NdirichletNodes;
 229         cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 230                 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 231     }
 232
 233         //creation de la matrice
 234     if(!_FECalculation)
 235         MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
 236     else
 237         MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
 238
 239         VecCreate(PETSC_COMM_SELF, &_Tk);
 240
 241     if(!_FECalculation)
 242         VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
 243     else
 244         VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
 245
 246         VecSetFromOptions(_Tk);
 247         VecDuplicate(_Tk, &_Tkm1);
 248         VecDuplicate(_Tk, &_deltaT);
 249         VecDuplicate(_Tk, &_b);//RHS of the linear system
 250
 251         //Linear solver
 252         KSPCreate(PETSC_COMM_SELF, &_ksp);
 253         KSPSetType(_ksp, _ksptype);
 254         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 255         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 256         KSPGetPC(_ksp, &_pc);
 257         PCSetType(_pc, _pctype);
 258
 259     //Checking whether all boundaries are Neumann boundaries
 260     map<string, LimitField>::iterator it = _limitField.begin();
 261     while(it != _limitField.end() and (it->second).bcType == Neumann)
 262         it++;
 263     _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);
 264     //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
 265     if(_onlyNeumannBC)
 266     {
 267         std::cout<<"## Warning all boundary conditions are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 268         std::cout<<"## As a consequence we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 269         *_runLogFile<<"## Warning all boundary condition are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 270         *_runLogFile<<"## As a consequence we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 271
 272                 //Check that the matrix is symmetric
 273                 PetscBool isSymetric;
 274                 MatIsSymmetric(_mat,_precision,&isSymetric);
 275                 if(!isSymetric)
 276                         {
 277                                 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
 278                                 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 279                         }
 280                 MatNullSpace nullsp;
 281                 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
 282                 MatSetNullSpace(_A, nullsp);
 283                 MatSetTransposeNullSpace(_A, nullsp);
 284                 MatNullSpaceDestroy(&nullsp);
 285                 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
 286                 //PCFactorSetShiftAmount(_pc,1e-10);
 287     }
 288
 289         _initializedMemory=true;
 290
 291 }
 292
 293 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
 294         if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
 295                 computeDiffusionMatrix(stop);
 296
 297         _dt_src=computeRHS(stop);
 298         stop=false;
 299         return _dt_src;
 300 }
 301
 302 Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values){
 303     vector< Matrix > matrices(_Ndim);
 304
 305     for (int idim=0; idim<_Ndim;idim++){
 306         matrices[idim]=M.deepCopy();
 307         for (int jdim=0; jdim<_Ndim+1;jdim++)
 308                         matrices[idim](jdim,idim) = values[jdim] ;
 309     }
 310
 311         Vector result(_Ndim);
 312     for (int idim=0; idim<_Ndim;idim++)
 313         result[idim] = matrices[idim].determinant();
 314
 315         return result;
 316 }
 317
 318 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
 319 {
 320     double result;
 321
 322     if(_FECalculation)
 323         result=computeDiffusionMatrixFE(stop);
 324     else
 325         result=computeDiffusionMatrixFV(stop);
 326
 327     //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
 328     //update value here if variable  heat transfer coefficient
 329     if(_heatTransfertCoeff>_precision)
 330         MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
 331
 332     if(_verbose or _system)
 333         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 334
 335     return  result;
 336 }
 337
 338 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
 339         Cell Cj;
 340         string nameOfGroup;
 341         double dn;
 342         MatZeroEntries(_A);
 343         VecZeroEntries(_b);
 344
 345     Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
 346     std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
 347     std::vector< int > nodeIds(_Ndim+1);//cell node Ids
 348     std::vector< Node > nodes(_Ndim+1);//cell nodes
 349     int i_int, j_int; //index of nodes j and k considered as unknown nodes
 350     bool dirichletCell_treated;
 351
 352     std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
 353     for (int idim=0; idim<_Ndim+1;idim++)
 354         values[idim][idim]=1;
 355
 356     /* parameters for boundary treatment */
 357     vector< double > valuesBorder(_Ndim+1);
 358     Vector GradShapeFuncBorder(_Ndim+1);
 359
 360         for (int j=0; j<_Nmailles;j++)
 361     {
 362                 Cj = _mesh.getCell(j);
 363
 364         for (int idim=0; idim<_Ndim+1;idim++){
 365             nodeIds[idim]=Cj.getNodeId(idim);
 366             nodes[idim]=_mesh.getNode(nodeIds[idim]);
 367             for (int jdim=0; jdim<_Ndim;jdim++)
 368                 M(idim,jdim)=nodes[idim].getPoint()[jdim];
 369             M(idim,_Ndim)=1;
 370         }
 371         for (int idim=0; idim<_Ndim+1;idim++)
 372             GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
 373
 374         /* Loop on the edges of the cell */
 375         for (int idim=0; idim<_Ndim+1;idim++)
 376         {
 377             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
 378             {//First node of the edge is not Dirichlet node
 379                 i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 380                 dirichletCell_treated=false;
 381                 for (int jdim=0; jdim<_Ndim+1;jdim++)
 382                 {
 383                     if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
 384                     {//Second node of the edge is not Dirichlet node
 385                         j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 386                         MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
 387                     }
 388                     else if (!dirichletCell_treated)
 389                     {//Second node of the edge is a Dirichlet node
 390                         dirichletCell_treated=true;
 391                         for (int kdim=0; kdim<_Ndim+1;kdim++)
 392                         {
 393                                                         std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
 394                                                         if( it != _dirichletBoundaryValues.end() )
 395                             {
 396                                 if( _dirichletValuesSet )
 397                                     valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
 398                                 else
 399                                     valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
 400                             }
 401                             else
 402                                 valuesBorder[kdim]=0;
 403                         }
 404                         GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
 405                         double coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 406                         VecSetValue(_b,i_int,coeff, ADD_VALUES);
 407                     }
 408                 }
 409             }
 410         }
 411         }
 412
 413     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 414         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 415         VecAssemblyBegin(_b);
 416         VecAssemblyEnd(_b);
 417
 418         _diffusionMatrixSet=true;
 419     stop=false ;
 420
 421         return INFINITY;
 422 }
 423
 424 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
 425         long nbFaces = _mesh.getNumberOfFaces();
 426         Face Fj;
 427         Cell Cell1,Cell2;
 428         string nameOfGroup;
 429         double inv_dxi, inv_dxj;
 430         double barycenterDistance;
 431         Vector normale(_Ndim);
 432         double dn;
 433         PetscInt idm, idn;
 434         std::vector< int > idCells;
 435         MatZeroEntries(_A);
 436         VecZeroEntries(_b);
 437
 438         for (int j=0; j<nbFaces;j++){
 439                 Fj = _mesh.getFace(j);
 440
 441                 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 442                 idCells = Fj.getCellsId();
 443                 Cell1 = _mesh.getCell(idCells[0]);
 444                 idm = idCells[0];
 445         for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 446             if (j == Cell1.getFacesId()[l]){
 447                 for (int idim = 0; idim < _Ndim; ++idim)
 448                     normale[idim] = Cell1.getNormalVector(l,idim);
 449                 break;
 450             }
 451         }
 452
 453                 //Compute velocity at the face Fj
 454                 dn=_conductivity*(_DiffusionTensor*normale)*normale;
 455
 456                 // compute 1/dxi = volume of Ci/area of Fj
 457         inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 458
 459                 // If Fj is on the boundary
 460                 if (Fj.getNumberOfCells()==1) {
 461                         if(_verbose )
 462                         {
 463                                 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 464                                 for(int p=0; p<_Ndim; p++)
 465                                         cout << normale[p] << ",";
 466                                 cout << ") "<<endl;
 467                         }
 468
 469             std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
 470             if( it != _dirichletBoundaryValues.end() )
 471             {
 472                 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 473                 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                                     , ADD_VALUES);
 474                 VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
 475             }
 476             else
 477             {
 478                 nameOfGroup = Fj.getGroupName();
 479
 480                 if (_limitField[nameOfGroup].bcType==Neumann){//Nothing to do
 481                 }
 482                 else if(_limitField[nameOfGroup].bcType==Dirichlet){
 483                     barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 484                     MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 485                     VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 486                 }
 487                 else {
 488                     stop=true ;
 489                     cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 490                     cout<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 491                     cout<<"Accepted boundary conditions are Neumann "<<Neumann<< " and Dirichlet "<<Dirichlet<<endl;
 492                     *_runLogFile<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 493                     _runLogFile->close();
 494                     throw CdmathException("Boundary condition not accepted");
 495                 }
 496             }
 497                         // if Fj is inside the domain
 498                 } else  if (Fj.getNumberOfCells()==2 ){
 499                         if(_verbose )
 500                         {
 501                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 502                                 cout << " ; vecteur normal=(";
 503                                 for(int p=0; p<_Ndim; p++)
 504                                         cout << normale[p] << ",";
 505                                 cout << ") "<<endl;
 506                         }
 507                         Cell2 = _mesh.getCell(idCells[1]);
 508                         idn = idCells[1];
 509                         if (_Ndim > 1)
 510                                 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 511                         else
 512                                 inv_dxj = 1/Cell2.getMeasure();
 513
 514                         barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 515
 516                         MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 517                         MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 518                         MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 519                         MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 520                 }
 521                 else
 522         {
 523             *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
 524                         throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 525         }
 526         }
 527
 528         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 529         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 530         VecAssemblyBegin(_b);
 531         VecAssemblyEnd(_b);
 532
 533         _diffusionMatrixSet=true;
 534     stop=false ;
 535
 536     return INFINITY;
 537 }
 538
 539 double StationaryDiffusionEquation::computeRHS(bool & stop)//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS)
 540 {
 541         VecAssemblyBegin(_b);
 542
 543     if(!_FECalculation)
 544         for (int i=0; i<_Nmailles;i++)
 545             VecSetValue(_b,i,_heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i),ADD_VALUES);
 546     else
 547         {
 548             Cell Ci;
 549             std::vector< int > nodesId;
 550             for (int i=0; i<_Nmailles;i++)
 551             {
 552                 Ci=_mesh.getCell(i);
 553                 nodesId=Ci.getNodesId();
 554                 for (int j=0; j<nodesId.size();j++)
 555                     if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
 556                     {
 557                         double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
 558                         VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);//assumes node numbering starts with unknown nodes. otherwise unknownNodes.index(j)
 559                     }
 560             }
 561         }
 562
 563         VecAssemblyEnd(_b);
 564
 565     if(_verbose or _system)
 566         VecView(_b,PETSC_VIEWER_STDOUT_SELF);
 567
 568     stop=false ;
 569         return INFINITY;
 570 }
 571
 572 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
 573 {
 574         bool stop;
 575
 576     //Only implicit scheme considered
 577     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 578     MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 579
 580 #if PETSC_VERSION_GREATER_3_5
 581     KSPSetOperators(_ksp, _A, _A);
 582 #else
 583     KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 584 #endif
 585
 586     if(_conditionNumber)
 587         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 588     KSPSolve(_ksp, _b, _Tk);
 589
 590         KSPConvergedReason reason;
 591         KSPGetConvergedReason(_ksp,&reason);
 592     KSPGetIterationNumber(_ksp, &_PetscIts);
 593     double residu;
 594     KSPGetResidualNorm(_ksp,&residu);
 595         if (reason!=2 and reason!=3)
 596     {
 597         cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 598         cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 599         cout<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 600                 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 601         *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 602         *_runLogFile<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 603                 _runLogFile->close();
 604         converged = false;
 605         stop = true;
 606     }
 607     else{
 608         if( _MaxIterLinearSolver < _PetscIts)
 609             _MaxIterLinearSolver = _PetscIts;
 610         cout<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<<  _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
 611                 *_runLogFile<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<<  _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
 612         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 613         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 614
 615         _erreur_rel= 0;
 616         double Ti, dTi;
 617
 618         VecAssemblyBegin(_Tk);
 619         VecAssemblyEnd(  _Tk);
 620         VecAssemblyBegin(_deltaT);
 621         VecAssemblyEnd(  _deltaT);
 622
 623         if(!_FECalculation)
 624             for(int i=0; i<_Nmailles; i++)
 625             {
 626                 VecGetValues(_deltaT, 1, &i, &dTi);
 627                 VecGetValues(_Tk, 1, &i, &Ti);
 628                 if(_erreur_rel < fabs(dTi/Ti))
 629                     _erreur_rel = fabs(dTi/Ti);
 630             }
 631         else
 632             for(int i=0; i<_NunknownNodes; i++)
 633             {
 634                 VecGetValues(_deltaT, 1, &i, &dTi);
 635                 VecGetValues(_Tk, 1, &i, &Ti);
 636                 if(_erreur_rel < fabs(dTi/Ti))
 637                     _erreur_rel = fabs(dTi/Ti);
 638             }
 639         stop=false;
 640         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 641     }
 642
 643         VecCopy(_Tk, _Tkm1);
 644
 645         return stop;
 646 }
 647
 648 void StationaryDiffusionEquation::setMesh(const Mesh &M)
 649 {
 650         if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
 651         {
 652         cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
 653                 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
 654                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 655                 _runLogFile->close();
 656                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect  dimension");
 657         }
 658
 659         _mesh=M;
 660         _Nmailles = _mesh.getNumberOfCells();
 661         _Nnodes =   _mesh.getNumberOfNodes();
 662
 663     cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
 664         *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
 665
 666         // find  maximum nb of neibourghs
 667     if(!_FECalculation)
 668     {
 669         _VV=Field ("Temperature", CELLS, _mesh, 1);
 670         _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
 671     }
 672     else
 673     {
 674         if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
 675                         cout<<"1D Finite element method on segments"<<endl;
 676         else if(_Ndim==2)
 677         {
 678                         if( _mesh.isTriangular() )//Mesh dim=2
 679                                 cout<<"2D Finite element method on triangles"<<endl;
 680                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 681                                 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 682                         else
 683                         {
 684                                 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
 685                                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 686                                 _runLogFile->close();
 687                                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 688                         }
 689         }
 690         else if(_Ndim==3)
 691         {
 692                         if( _mesh.isTetrahedral() )//Mesh dim=3
 693                                 cout<<"3D Finite element method on tetrahedra"<<endl;
 694                         else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
 695                                 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 696                         else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 697                                 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 698                         else
 699                         {
 700                                 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;
 701                                 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 702                                 _runLogFile->close();
 703                                 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 704                         }
 705         }
 706
 707                 _VV=Field ("Temperature", NODES, _mesh, 1);
 708
 709         _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 710         _boundaryNodeIds = _mesh.getBoundaryNodeIds();
 711         _NboundaryNodes=_boundaryNodeIds.size();
 712     }
 713
 714         _meshSet=true;
 715 }
 716
 717 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
 718 {
 719         //_maxPetscIts=maxIterationsPetsc;
 720         // set linear solver algorithm
 721         if (kspType==GMRES)
 722                 _ksptype = (char*)&KSPGMRES;
 723         else if (kspType==CG)
 724                 _ksptype = (char*)&KSPCG;
 725         else if (kspType==BCGS)
 726                 _ksptype = (char*)&KSPBCGS;
 727         else {
 728                 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 729                 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 730                 _runLogFile->close();
 731                 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
 732         }
 733         // set preconditioner
 734         if (pcType == NONE)
 735                 _pctype = (char*)&PCNONE;
 736         else if (pcType ==LU)
 737                 _pctype = (char*)&PCLU;
 738         else if (pcType == ILU)
 739                 _pctype = (char*)&PCILU;
 740         else if (pcType ==CHOLESKY)
 741                 _pctype = (char*)&PCCHOLESKY;
 742         else if (pcType == ICC)
 743                 _pctype = (char*)&PCICC;
 744         else {
 745                 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
 746                 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
 747                 _runLogFile->close();
 748                 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
 749         }
 750 }
 751
 752 bool StationaryDiffusionEquation::solveStationaryProblem()
 753 {
 754         if(!_initializedMemory)
 755         {
 756                 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
 757                 _runLogFile->close();
 758                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 759         }
 760         bool stop=false; // Does the Problem want to stop (error) ?
 761         bool converged=false; // has the newton scheme converged (end) ?
 762
 763         cout<< "!!! Running test case "<< _fileName << " using ";
 764         *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
 765
 766     if(!_FECalculation)
 767     {
 768         cout<< "Finite volumes method"<<endl<<endl;
 769                 *_runLogFile<< "Finite volumes method"<<endl<<endl;
 770         }
 771     else
 772         {
 773         cout<< "Finite elements method"<<endl<<endl;
 774                 *_runLogFile<< "Finite elements method"<< endl<<endl;
 775         }
 776
 777     computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
 778     if (stop){
 779         cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 780         *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 781                 _runLogFile->close();
 782        throw CdmathException("Failed computing diffusion matrix");
 783     }
 784     computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
 785     if (stop){
 786         cout << "Error : failed computing right hand side, stopping calculation"<< endl;
 787         *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
 788         throw CdmathException("Failed computing right hand side");
 789     }
 790     stop = iterateNewtonStep(converged);
 791     if (stop){
 792         cout << "Error : failed solving linear system, stopping calculation"<< endl;
 793         *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
 794                 _runLogFile->close();
 795         throw CdmathException("Failed solving linear system");
 796     }
 797
 798     _computationCompletedSuccessfully=true;
 799     save();
 800
 801         // Newton iteration loop for non linear problems
 802     /*
 803         while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
 804         {
 805         computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
 806         computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
 807         stop = iterateNewtonStep(converged);
 808         _NEWTON_its++;
 809         }
 810     if (stop){
 811         cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 812         *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 813         throw CdmathException("Failed solving a Newton iteration");
 814     }
 815     else if(_NEWTON_its==_maxNewtonIts){
 816         cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 817         *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 818         throw CdmathException("No convergence of Newton scheme");
 819     }
 820     else{
 821         cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 822         *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 823         save();
 824     }
 825     */
 826
 827     *_runLogFile<< "!!!!!! Computation successful"<< endl;
 828         _runLogFile->close();
 829
 830         return !stop;
 831 }
 832
 833 void StationaryDiffusionEquation::save(){
 834     cout<< "Saving numerical results"<<endl<<endl;
 835     *_runLogFile<< "Saving numerical results"<< endl<<endl;
 836
 837         string resultFile(_path+"/StationaryDiffusionEquation");//Results
 838
 839         resultFile+="_";
 840         resultFile+=_fileName;
 841
 842         // create mesh and component info
 843     string suppress ="rm -rf "+resultFile+"_*";
 844     system(suppress.c_str());//Nettoyage des précédents calculs identiques
 845
 846     if(_verbose or _system)
 847         VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
 848
 849     double Ti;
 850     if(!_FECalculation)
 851         for(int i=0; i<_Nmailles; i++)
 852             {
 853                 VecGetValues(_Tk, 1, &i, &Ti);
 854                 _VV(i)=Ti;
 855             }
 856     else
 857     {
 858         int globalIndex;
 859         for(int i=0; i<_NunknownNodes; i++)
 860         {
 861             VecGetValues(_Tk, 1, &i, &Ti);
 862             globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 863             _VV(globalIndex)=Ti;//Assumes node numbering starts with border nodes
 864         }
 865
 866         Node Ni;
 867         string nameOfGroup;
 868         for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
 869         {
 870             Ni=_mesh.getNode(_dirichletNodeIds[i]);
 871             nameOfGroup = Ni.getGroupName();
 872             _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 873         }
 874     }
 875
 876     _VV.setInfoOnComponent(0,"Temperature_(K)");
 877     switch(_saveFormat)
 878     {
 879         case VTK :
 880             _VV.writeVTK(resultFile);
 881             break;
 882         case MED :
 883             _VV.writeMED(resultFile);
 884             break;
 885         case CSV :
 886             _VV.writeCSV(resultFile);
 887             break;
 888     }
 889 }
 890 Field
 891 StationaryDiffusionEquation::getOutputTemperatureField()
 892 {
 893     if(!_computationCompletedSuccessfully)
 894         throw("Computation not performed yet or failed. No temperature field available");
 895     else
 896         return _VV;
 897 }
 898
 899 void StationaryDiffusionEquation::terminate()
 900 {
 901         VecDestroy(&_Tk);
 902         VecDestroy(&_Tkm1);
 903         VecDestroy(&_deltaT);
 904         VecDestroy(&_b);
 905         MatDestroy(&_A);
 906 }
 907 void
 908 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
 909 {
 910     _dirichletValuesSet=true;
 911     _dirichletBoundaryValues=dirichletBoundaryValues;
 912 }
 913