CoreFlows/Models/src/StationaryDiffusionEquation.cxx

   1 #include "StationaryDiffusionEquation.hxx"
   2 #include "DiffusionEquation.hxx"
   3 #include "Node.hxx"
   4 #include "SparseMatrixPetsc.hxx"
   5 #include "math.h"
   6 #include <algorithm>
   7 #include <fstream>
   8 #include <sstream>
   9
  10 using namespace std;
  11
  12 StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda,MPI_Comm comm){
  13         /* Initialisation of PETSC */
  14         //check if PETSC is already initialised
  15         PetscBool petscInitialized;
  16         PetscInitialized(&petscInitialized);
  17         if(!petscInitialized)
  18         {//check if MPI is already initialised
  19                 int mpiInitialized;
  20                 MPI_Initialized(&mpiInitialized);
  21                 if(mpiInitialized)
  22                         PETSC_COMM_WORLD = comm;
  23                 PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
  24         }
  25         MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
  26         MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
  27         PetscPrintf(PETSC_COMM_WORLD,"\n Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
  28         PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
  29         PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
  30
  31     if(lambda < 0.)
  32     {
  33         std::cout<<"Conductivity="<<lambda<<endl;
  34         throw CdmathException("!!!!!!!!Error : conductivity parameter lambda cannot  be negative");
  35     }
  36     if(dim<=0)
  37     {
  38         std::cout<<"Space dimension="<<dim<<endl;
  39         throw CdmathException("!!!!!!!!Error : parameter dim cannot  be negative");
  40     }
  41
  42     _FECalculation=FECalculation;
  43     _onlyNeumannBC=false;
  44
  45         _Ndim=dim;
  46         _nVar=1;//scalar prolem
  47         _dt_src=0;
  48         _diffusionMatrixSet=false;
  49         _initializedMemory=false;
  50
  51     //Mesh data
  52     _neibMaxNbCells=0;
  53     _meshSet=false;
  54     _neibMaxNbNodes=0;
  55
  56     //Boundary conditions
  57     _boundaryNodeIds=std::vector< int >(0);
  58     _dirichletNodeIds=std::vector< int >(0);
  59     _NboundaryNodes=0;
  60     _NdirichletNodes=0;
  61     _NunknownNodes=0;
  62     _dirichletValuesSet=false;
  63     _neumannValuesSet=false;
  64
  65     //Linear solver data
  66         _precision=1.e-6;
  67         _precision_Newton=_precision;
  68         _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
  69         _maxPetscIts=50;
  70         _maxNewtonIts=50;
  71         _NEWTON_its=0;
  72         int _PetscIts=0;//the number of iterations of the linear solver
  73         _ksptype = (char*)&KSPGMRES;
  74         if( _mpi_size>1)
  75                 _pctype = (char*)&PCNONE;
  76         else
  77                 _pctype = (char*)&PCLU;
  78
  79         _conditionNumber=false;
  80         _erreur_rel= 0;
  81
  82     //parameters for monitoring simulation
  83         _verbose = false;
  84         _system = false;
  85         _runLogFile=new ofstream;
  86
  87         //result save parameters
  88         _fileName = "StationaryDiffusionProblem";
  89         char result[ PATH_MAX ];//extracting current directory
  90         getcwd(result, PATH_MAX );
  91         _path=string( result );
  92         _saveFormat=VTK;
  93     _computationCompletedSuccessfully=false;
  94
  95     //heat transfer parameters
  96         _conductivity=lambda;
  97         _fluidTemperatureFieldSet=false;
  98         _fluidTemperature=0;
  99         _heatPowerFieldSet=false;
 100         _heatTransfertCoeff=0;
 101         _heatSource=0;
 102
 103     /* Default diffusion tensor is diagonal */
 104         _DiffusionTensor=Matrix(_Ndim);
 105         for(int idim=0;idim<_Ndim;idim++)
 106                 _DiffusionTensor(idim,idim)=_conductivity;
 107
 108     PetscPrintf(PETSC_COMM_WORLD,"\n Stationary diffusion problem with conductivity %.2f", lambda);
 109     if(FECalculation)
 110         PetscPrintf(PETSC_COMM_WORLD," and finite elements method\n\n");
 111     else
 112         PetscPrintf(PETSC_COMM_WORLD," and finite volumes method\n\n");
 113 }
 114
 115 void StationaryDiffusionEquation::initialize()
 116 {
 117         if(_mpi_rank==0)
 118         {
 119                 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 120
 121                 if(!_meshSet)
 122                         throw CdmathException("!!!!!!!!StationaryDiffusionEquation::initialize() set mesh first");
 123                 else
 124             {
 125                         cout<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
 126                 *_runLogFile<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
 127                 if(!_FECalculation)
 128                 {
 129                     cout<< "Finite volumes method"<<endl<<endl;
 130                     *_runLogFile<< "Finite volumes method"<<endl<<endl;
 131                 }
 132                 else
 133                 {
 134                     cout<< "Finite elements method"<<endl<<endl;
 135                     *_runLogFile<< "Finite elements method"<<endl<<endl;
 136                 }
 137             }
 138
 139                 /**************** Field creation *********************/
 140
 141                 if(!_heatPowerFieldSet){
 142                         _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
 143                         for(int i =0; i<_VV.getNumberOfElements(); i++)
 144                                 _heatPowerField(i) = _heatSource;
 145                 _heatPowerFieldSet=true;
 146             }
 147                 if(!_fluidTemperatureFieldSet){
 148                         _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
 149                         for(int i =0; i<_VV.getNumberOfElements(); i++)
 150                                 _fluidTemperatureField(i) = _fluidTemperature;
 151                 _fluidTemperatureFieldSet=true;
 152                 }
 153
 154             /* Détection des noeuds frontière avec une condition limite de Dirichlet */
 155             if(_FECalculation)
 156             {
 157                 if(_NboundaryNodes==_Nnodes)
 158                     cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
 159
 160                 for(int i=0; i<_NboundaryNodes; i++)
 161                 {
 162                     std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
 163                     if( it != _dirichletBoundaryValues.end() )
 164                         _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 165                     else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
 166                     {
 167                         cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
 168                         *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
 169                         _runLogFile->close();
 170                         throw CdmathException("Missing boundary group");
 171                     }
 172                     else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
 173                     {
 174                         cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
 175                         cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
 176                         *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
 177                         *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
 178                         _runLogFile->close();
 179                         throw CdmathException("Missing boundary condition");
 180                     }
 181                     else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
 182                         _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 183                     else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
 184                     {
 185                         cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
 186                         cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
 187                         *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
 188                         *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
 189                         _runLogFile->close();
 190                         throw CdmathException("Wrong boundary condition");
 191                     }
 192                 }
 193                 _NdirichletNodes=_dirichletNodeIds.size();
 194                 _NunknownNodes=_Nnodes - _NdirichletNodes;
 195                 cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 196                         *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 197             }
 198         }
 199
 200     if(!_FECalculation)
 201                 _globalNbUnknowns = _Nmailles*_nVar;
 202     else{
 203         MPI_Bcast(&_NunknownNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 204                 _globalNbUnknowns = _NunknownNodes*_nVar;
 205         }
 206
 207         /* Vectors creations */
 208         VecCreate(PETSC_COMM_WORLD, &_Tk);//main unknown
 209     VecSetSizes(_Tk,PETSC_DECIDE,_globalNbUnknowns);
 210         VecSetFromOptions(_Tk);
 211         VecGetLocalSize(_Tk, &_localNbUnknowns);
 212
 213         VecDuplicate(_Tk, &_Tkm1);
 214         VecDuplicate(_Tk, &_deltaT);
 215         VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
 216
 217         /* Matrix creation */
 218         MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
 219
 220         /* Local sequential vector creation */
 221         if(_mpi_size>1 && _mpi_rank == 0)
 222                 VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Tk_seq);//For saving results on proc 0
 223         if(_mpi_size==0)
 224                 _Tk_seq=_Tk;
 225         VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq);
 226
 227         //Linear solver
 228         KSPCreate(PETSC_COMM_WORLD, &_ksp);
 229         KSPSetType(_ksp, _ksptype);
 230         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 231         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 232         KSPGetPC(_ksp, &_pc);
 233         PCSetType(_pc, _pctype);
 234
 235     //Checking whether all boundary conditions are Neumann boundary condition
 236     //if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
 237     if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
 238     {
 239         map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
 240         while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
 241             it++;
 242         _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
 243     }
 244     else
 245         if(_FECalculation)
 246             _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
 247         else
 248             _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
 249
 250     //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
 251     if(_onlyNeumannBC)
 252     {
 253         std::cout<<"### Warning : all boundary conditions are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 254         std::cout<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
 255         std::cout<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 256         *_runLogFile<<"### Warning : all boundary condition are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
 257         *_runLogFile<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
 258         *_runLogFile<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
 259
 260                 MatNullSpace nullsp;
 261                 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
 262                 MatSetNullSpace(_A, nullsp);
 263                 MatSetTransposeNullSpace(_A, nullsp);
 264                 MatNullSpaceDestroy(&nullsp);
 265                 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
 266                 //PCFactorSetShiftAmount(_pc,1e-10);
 267     }
 268         _initializedMemory=true;
 269 }
 270
 271 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
 272         if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
 273                 computeDiffusionMatrix(stop);
 274
 275         _dt_src=computeRHS(stop);
 276         stop=false;
 277         return _dt_src;
 278 }
 279
 280 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
 281 {
 282     double result;
 283
 284         MatZeroEntries(_A);
 285         VecZeroEntries(_b);
 286
 287     if(_FECalculation)
 288         result=computeDiffusionMatrixFE(stop);
 289     else
 290         result=computeDiffusionMatrixFV(stop);
 291
 292     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 293         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 294
 295     //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
 296     //update value here if variable  heat transfer coefficient
 297     if(_heatTransfertCoeff>_precision)
 298         MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
 299
 300     if(_verbose or _system)
 301         MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
 302
 303     return  result;
 304 }
 305
 306 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
 307         if(_mpi_rank == 0)
 308                 {
 309                 Cell Cj;
 310                 string nameOfGroup;
 311                 double coeff;
 312
 313             Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
 314             std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
 315             std::vector< int > nodeIds(_Ndim+1);//cell node Ids
 316             std::vector< Node > nodes(_Ndim+1);//cell nodes
 317             int i_int, j_int; //index of nodes j and k considered as unknown nodes
 318             bool dirichletCell_treated;
 319
 320             std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
 321             for (int idim=0; idim<_Ndim+1;idim++)
 322                 values[idim][idim]=1;
 323
 324             /* parameters for boundary treatment */
 325             vector< double > valuesBorder(_Ndim+1);
 326             Vector GradShapeFuncBorder(_Ndim+1);
 327
 328                 for (int j=0; j<_Nmailles;j++)
 329             {
 330                         Cj = _mesh.getCell(j);
 331
 332                 for (int idim=0; idim<_Ndim+1;idim++){
 333                     nodeIds[idim]=Cj.getNodeId(idim);
 334                     nodes[idim]=_mesh.getNode(nodeIds[idim]);
 335                     for (int jdim=0; jdim<_Ndim;jdim++)
 336                         M(idim,jdim)=nodes[idim].getPoint()[jdim];
 337                     M(idim,_Ndim)=1;
 338                 }
 339                 for (int idim=0; idim<_Ndim+1;idim++)
 340                     GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim);
 341
 342                 /* Loop on the edges of the cell */
 343                 for (int idim=0; idim<_Ndim+1;idim++)
 344                 {
 345                     if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
 346                     {//First node of the edge is not Dirichlet node
 347                         i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 348                         dirichletCell_treated=false;
 349                         for (int jdim=0; jdim<_Ndim+1;jdim++)
 350                         {
 351                             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
 352                             {//Second node of the edge is not Dirichlet node
 353                                 j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 354                                 MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
 355                             }
 356                             else if (!dirichletCell_treated)
 357                             {//Second node of the edge is a Dirichlet node
 358                                 dirichletCell_treated=true;
 359                                 for (int kdim=0; kdim<_Ndim+1;kdim++)
 360                                 {
 361                                                                 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
 362                                                                 if( it != _dirichletBoundaryValues.end() )
 363                                     {
 364                                         if( _dirichletValuesSet )
 365                                             valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
 366                                         else
 367                                             valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
 368                                     }
 369                                     else
 370                                         valuesBorder[kdim]=0;
 371                                 }
 372                                 GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
 373                                 coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 374                                 VecSetValue(_b,i_int,coeff, ADD_VALUES);
 375                             }
 376                         }
 377                     }
 378                 }
 379                 }
 380
 381             //Calcul de la contribution de la condition limite de Neumann au second membre
 382             if( _NdirichletNodes !=_NboundaryNodes)
 383             {
 384                 vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
 385                 int NboundaryFaces=boundaryFaces.size();
 386                 for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
 387                 {
 388                     Face Fi = _mesh.getFace(boundaryFaces[i]);
 389                     for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
 390                     {
 391                         if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
 392                         {
 393                             j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 394                             if( _neumannValuesSet )
 395                                 coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
 396                             else
 397                                 coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
 398                             VecSetValue(_b, j_int, coeff, ADD_VALUES);
 399                         }
 400                     }
 401                 }
 402             }
 403         }
 404
 405         if(_onlyNeumannBC)      //Check that the matrix is symmetric
 406         {
 407                 PetscBool isSymetric;
 408                 MatIsSymmetric(_A,_precision,&isSymetric);
 409                 if(!isSymetric)
 410                 {
 411                         cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
 412                         throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 413                 }
 414         }
 415
 416         _diffusionMatrixSet=true;
 417     stop=false ;
 418
 419         return INFINITY;
 420 }
 421
 422 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
 423         if(_mpi_rank == 0)
 424         {
 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
 436                 for (int j=0; j<nbFaces;j++){
 437                         Fj = _mesh.getFace(j);
 438
 439                         // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 440                         idCells = Fj.getCellsId();
 441                         Cell1 = _mesh.getCell(idCells[0]);
 442                         idm = idCells[0];
 443                 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 444                     if (j == Cell1.getFacesId()[l]){
 445                         for (int idim = 0; idim < _Ndim; ++idim)
 446                             normale[idim] = Cell1.getNormalVector(l,idim);
 447                         break;
 448                     }
 449                 }
 450
 451                         //Compute velocity at the face Fj
 452                         dn=(_DiffusionTensor*normale)*normale;
 453
 454                         // compute 1/dxi = volume of Ci/area of Fj
 455                 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 456
 457                         // If Fj is on the boundary
 458                         if (Fj.getNumberOfCells()==1) {
 459                                 if(_verbose )
 460                                 {
 461                                         cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 462                                         for(int p=0; p<_Ndim; p++)
 463                                                 cout << normale[p] << ",";
 464                                         cout << ") "<<endl;
 465                                 }
 466
 467                     std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
 468                     if( it != _dirichletBoundaryValues.end() )
 469                     {
 470                         barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 471                         MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                                     , ADD_VALUES);
 472                         VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
 473                     }
 474                     else
 475                     {
 476                         nameOfGroup = Fj.getGroupName();
 477
 478                         if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
 479                             VecSetValue(_b,idm,    -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
 480                         }
 481                         else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
 482                             barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 483                             MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 484                             VecSetValue(_b,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 485                         }
 486                         else {
 487                             stop=true ;
 488                             cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 489                             cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 490                             cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
 491                             *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 492                             _runLogFile->close();
 493                             throw CdmathException("Boundary condition not set");
 494                         }
 495                     }
 496                                 // if Fj is inside the domain
 497                         } else  if (Fj.getNumberOfCells()==2 ){
 498                                 if(_verbose )
 499                                 {
 500                                         cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 501                                         cout << " ; vecteur normal=(";
 502                                         for(int p=0; p<_Ndim; p++)
 503                                                 cout << normale[p] << ",";
 504                                         cout << ") "<<endl;
 505                                 }
 506                                 Cell2 = _mesh.getCell(idCells[1]);
 507                                 idn = idCells[1];
 508                                 if (_Ndim > 1)
 509                                         inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 510                                 else
 511                                         inv_dxj = 1/Cell2.getMeasure();
 512
 513                                 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 514
 515                                 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 516                                 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 517                                 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 518                                 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 519                         }
 520                         else
 521                 {
 522                     *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
 523                                 throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 524                 }
 525                 }
 526         }
 527
 528         if(_onlyNeumannBC)      //Check that the matrix is symmetric
 529         {
 530                 PetscBool isSymetric;
 531                 MatIsSymmetric(_A,_precision,&isSymetric);
 532                 if(!isSymetric)
 533                 {
 534                         cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
 535                         throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 536                 }
 537         }
 538
 539         _diffusionMatrixSet=true;
 540     stop=false ;
 541
 542     return INFINITY;
 543 }
 544
 545 double StationaryDiffusionEquation::computeRHS(bool & stop)//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeDiffusionMatrix
 546 {
 547
 548         if(_mpi_rank == 0)
 549         {
 550             if(!_FECalculation)
 551                 for (int i=0; i<_Nmailles;i++)
 552                     VecSetValue(_b,i, _heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i) ,ADD_VALUES);
 553             else
 554                 {
 555                     Cell Ci;
 556                     std::vector< int > nodesId;
 557                     for (int i=0; i<_Nmailles;i++)
 558                     {
 559                         Ci=_mesh.getCell(i);
 560                         nodesId=Ci.getNodesId();
 561                         for (int j=0; j<nodesId.size();j++)
 562                             if(!_mesh.isBorderNode(nodesId[j]))
 563                             {
 564                                 double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
 565                                 VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 566                             }
 567                     }
 568                 }
 569         }
 570         VecAssemblyBegin(_b);
 571         VecAssemblyEnd(_b);
 572
 573     if(_verbose or _system)
 574         VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
 575
 576     stop=false ;
 577         return INFINITY;
 578 }
 579
 580 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
 581 {
 582         bool stop;
 583
 584     //Only implicit scheme considered
 585
 586 #if PETSC_VERSION_GREATER_3_5
 587     KSPSetOperators(_ksp, _A, _A);
 588 #else
 589     KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 590 #endif
 591
 592     if(_conditionNumber)
 593         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 594     KSPSolve(_ksp, _b, _Tk);
 595
 596         KSPConvergedReason reason;
 597         KSPGetConvergedReason(_ksp,&reason);
 598     KSPGetIterationNumber(_ksp, &_PetscIts);
 599     double residu;
 600     KSPGetResidualNorm(_ksp,&residu);
 601         if (reason!=2 and reason!=3)
 602     {
 603         PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc.");
 604         PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Itérations maximales %d atteintes, résidu = %1.2e, précision demandée= %1.2e",_maxPetscIts,residu,_precision);
 605         PetscPrintf(PETSC_COMM_WORLD,"Solver used %s, preconditioner %s, Final number of iteration = %d",_ksptype,_pctype,_PetscIts);
 606                 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
 607         *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
 608         *_runLogFile<<"Solver used "<<  _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
 609                 _runLogFile->close();
 610         converged = false;
 611         stop = true;
 612     }
 613     else{
 614         if( _MaxIterLinearSolver < _PetscIts)
 615             _MaxIterLinearSolver = _PetscIts;
 616         PetscPrintf(PETSC_COMM_WORLD,"## Système linéaire résolu en %d itérations par le solveur %s et le preconditioneur %s, précision demandée = %1.2e",_PetscIts,_ksptype,_pctype,_precision);
 617                 *_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;
 618         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 619         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 620
 621         VecAssemblyBegin(_Tk);
 622         VecAssemblyEnd(  _Tk);
 623         VecAssemblyBegin(_deltaT);
 624         VecAssemblyEnd(  _deltaT);
 625
 626                 if(_verbose)
 627                         PetscPrintf(PETSC_COMM_WORLD,"Début calcul de l'erreur maximale");
 628
 629                 VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
 630
 631                 if(_verbose)
 632                         PetscPrintf(PETSC_COMM_WORLD,"Fin calcul de la variation relative, erreur maximale : %1.2e", _erreur_rel );
 633         stop=false;
 634         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 635     }
 636
 637         VecCopy(_Tk, _Tkm1);
 638
 639         return stop;
 640 }
 641
 642 void StationaryDiffusionEquation::setMesh(const Mesh &M)
 643 {
 644         if(_mpi_rank==0)
 645         {
 646                 if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
 647                 {
 648                 cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
 649                         *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
 650                         *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 651                         _runLogFile->close();
 652                         throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect  dimension");
 653                 }
 654
 655                 _mesh=M;
 656                 _Nmailles = _mesh.getNumberOfCells();
 657                 _Nnodes =   _mesh.getNumberOfNodes();
 658
 659             cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
 660                 *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
 661
 662                 // find  maximum nb of neibourghs
 663             if(!_FECalculation)
 664             {
 665                 _VV=Field ("Temperature", CELLS, _mesh, 1);
 666                 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
 667             }
 668             else
 669             {
 670                 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
 671                                 cout<<"1D Finite element method on segments"<<endl;
 672                 else if(_Ndim==2)
 673                 {
 674                                 if( _mesh.isTriangular() )//Mesh dim=2
 675                                         cout<<"2D Finite element method on triangles"<<endl;
 676                                 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 677                                         cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 678                                 else
 679                                 {
 680                                         cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension  "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
 681                                         *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 682                                         _runLogFile->close();
 683                                         throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 684                                 }
 685                 }
 686                 else if(_Ndim==3)
 687                 {
 688                                 if( _mesh.isTetrahedral() )//Mesh dim=3
 689                                         cout<<"3D Finite element method on tetrahedra"<<endl;
 690                                 else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
 691                                         cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 692                                 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 693                                         cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
 694                                 else
 695                                 {
 696                                         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;
 697                                         *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
 698                                         _runLogFile->close();
 699                                         throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
 700                                 }
 701                 }
 702
 703                         _VV=Field ("Temperature", NODES, _mesh, 1);
 704
 705                 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 706                 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
 707                 _NboundaryNodes=_boundaryNodeIds.size();
 708             }
 709         }
 710
 711         /* MPI distribution parameters */
 712         int nbVoisinsMax;//Mettre en attribut ?
 713         if(!_FECalculation){
 714                 MPI_Bcast(&_Nmailles      , 1, MPI_INT, 0, PETSC_COMM_WORLD);
 715                 MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 716                 nbVoisinsMax = _neibMaxNbCells;
 717         }
 718         else{
 719                 MPI_Bcast(&_Nnodes        , 1, MPI_INT, 0, PETSC_COMM_WORLD);
 720                 MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 721                 nbVoisinsMax = _neibMaxNbNodes;
 722         }
 723     _d_nnz = (nbVoisinsMax+1)*_nVar;
 724     _o_nnz =  nbVoisinsMax   *_nVar;
 725
 726         _meshSet=true;
 727 }
 728
 729 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
 730 {
 731         //_maxPetscIts=maxIterationsPetsc;
 732         // set linear solver algorithm
 733         if (kspType==GMRES)
 734                 _ksptype = (char*)&KSPGMRES;
 735         else if (kspType==CG)
 736                 _ksptype = (char*)&KSPCG;
 737         else if (kspType==BCGS)
 738                 _ksptype = (char*)&KSPBCGS;
 739         else {
 740                 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 741                 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 742                 _runLogFile->close();
 743                 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
 744         }
 745         // set preconditioner
 746         if (pcType == NONE)
 747                 _pctype = (char*)&PCNONE;
 748         else if (pcType ==LU)
 749                 _pctype = (char*)&PCLU;
 750         else if (pcType == ILU)
 751                 _pctype = (char*)&PCILU;
 752         else if (pcType ==CHOLESKY)
 753                 _pctype = (char*)&PCCHOLESKY;
 754         else if (pcType == ICC)
 755                 _pctype = (char*)&PCICC;
 756         else {
 757                 cout << "!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
 758                 *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
 759                 _runLogFile->close();
 760                 throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
 761         }
 762 }
 763
 764 bool StationaryDiffusionEquation::solveStationaryProblem()
 765 {
 766         if(!_initializedMemory)
 767         {
 768                 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
 769                 _runLogFile->close();
 770                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 771         }
 772         bool stop=false; // Does the Problem want to stop (error) ?
 773         bool converged=false; // has the newton scheme converged (end) ?
 774
 775         PetscPrintf(PETSC_COMM_WORLD,"!!! Running test case %s using ",_fileName);
 776         *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
 777
 778     if(!_FECalculation)
 779     {
 780         PetscPrintf(PETSC_COMM_WORLD,"Finite volumes method\n\n");
 781                 *_runLogFile<< "Finite volumes method"<<endl<<endl;
 782         }
 783     else
 784         {
 785         PetscPrintf(PETSC_COMM_WORLD,"Finite elements method\n\n");
 786                 *_runLogFile<< "Finite elements method"<< endl<<endl;
 787         }
 788
 789     computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
 790     if (stop){
 791         PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing diffusion matrix, stopping calculation\n");
 792         *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
 793                 _runLogFile->close();
 794        throw CdmathException("Failed computing diffusion matrix");
 795     }
 796     computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
 797     if (stop){
 798         PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing right hand side, stopping calculation\n");
 799         *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
 800         throw CdmathException("Failed computing right hand side");
 801     }
 802     stop = iterateNewtonStep(converged);
 803     if (stop){
 804         PetscPrintf(PETSC_COMM_WORLD,"Error : failed solving linear system, stopping calculation\n");
 805         *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
 806                 _runLogFile->close();
 807         throw CdmathException("Failed solving linear system");
 808     }
 809
 810     _computationCompletedSuccessfully=true;
 811     save();
 812
 813         // Newton iteration loop for non linear problems
 814     /*
 815         while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
 816         {
 817         computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
 818         computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
 819         stop = iterateNewtonStep(converged);
 820         _NEWTON_its++;
 821         }
 822     if (stop){
 823         cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 824         *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
 825         throw CdmathException("Failed solving a Newton iteration");
 826     }
 827     else if(_NEWTON_its==_maxNewtonIts){
 828         cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 829         *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
 830         throw CdmathException("No convergence of Newton scheme");
 831     }
 832     else{
 833         cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 834         *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
 835         save();
 836     }
 837     */
 838
 839     *_runLogFile<< "!!!!!! Computation successful !!!!!!"<< endl;
 840         _runLogFile->close();
 841
 842         return !stop;
 843 }
 844
 845 void StationaryDiffusionEquation::save(){
 846     PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n");
 847     *_runLogFile<< "Saving numerical results"<< endl<<endl;
 848
 849         string resultFile(_path+"/StationaryDiffusionEquation");//Results
 850
 851         resultFile+="_";
 852         resultFile+=_fileName;
 853
 854         // create mesh and component info
 855     string suppress ="rm -rf "+resultFile+"_*";
 856     system(suppress.c_str());//Nettoyage des précédents calculs identiques
 857
 858         if(_mpi_size>1){
 859                 VecScatterBegin(_scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
 860                 VecScatterEnd(  _scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
 861         }
 862
 863     if(_verbose or _system)
 864         VecView(_Tk,PETSC_VIEWER_STDOUT_WORLD);
 865
 866         if(_mpi_rank==0){
 867             double Ti;
 868             if(!_FECalculation)
 869                 for(int i=0; i<_Nmailles; i++)
 870                     {
 871                                         if(_mpi_size>1)
 872                                                 VecGetValues(_Tk_seq, 1, &i, &Ti);
 873                                         else
 874                                                 VecGetValues(_Tk    , 1, &i, &Ti);
 875                         _VV(i)=Ti;
 876                     }
 877             else
 878             {
 879                 int globalIndex;
 880                 for(int i=0; i<_NunknownNodes; i++)
 881                 {
 882                                 if(_mpi_size>1)
 883                                         VecGetValues(_Tk_seq, 1, &i, &Ti);
 884                                 else
 885                                         VecGetValues(_Tk    , 1, &i, &Ti);
 886                     globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
 887                     _VV(globalIndex)=Ti;
 888                 }
 889
 890                 Node Ni;
 891                 string nameOfGroup;
 892                 for(int i=0; i<_NdirichletNodes; i++)
 893                 {
 894                     Ni=_mesh.getNode(_dirichletNodeIds[i]);
 895                     nameOfGroup = Ni.getGroupName();
 896                     _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 897                 }
 898             }
 899
 900             _VV.setInfoOnComponent(0,"Temperature_(K)");
 901             switch(_saveFormat)
 902             {
 903                 case VTK :
 904                     _VV.writeVTK(resultFile);
 905                     break;
 906                 case MED :
 907                     _VV.writeMED(resultFile);
 908                     break;
 909                 case CSV :
 910                     _VV.writeCSV(resultFile);
 911                     break;
 912             }
 913         }
 914 }
 915
 916 void StationaryDiffusionEquation::terminate()
 917 {
 918         VecDestroy(&_Tk);
 919         VecDestroy(&_Tkm1);
 920         VecDestroy(&_deltaT);
 921         VecDestroy(&_b);
 922         MatDestroy(&_A);
 923         if(_mpi_size>1 && _mpi_rank == 0)
 924                 VecDestroy(&_Tk_seq);
 925
 926         PetscBool petscInitialized;
 927         PetscInitialized(&petscInitialized);
 928         if(petscInitialized)
 929                 PetscFinalize();
 930
 931         delete _runLogFile;
 932 }
 933 void
 934 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
 935 {
 936     _dirichletValuesSet=true;
 937     _dirichletBoundaryValues=dirichletBoundaryValues;
 938 }
 939
 940 void
 941 StationaryDiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
 942 {
 943     _neumannValuesSet=true;
 944     _neumannBoundaryValues=neumannBoundaryValues;
 945 }
 946
 947 double
 948 StationaryDiffusionEquation::getConditionNumber(bool isSingular, double tol) const
 949 {
 950   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 951   return A.getConditionNumber( isSingular, tol);
 952 }
 953 std::vector< double >
 954 StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) const
 955 {
 956   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 957
 958   if(_FECalculation)//We need to scale the FE matrix, otherwise the eigenvalues go to zero as the mesh is refined
 959   {
 960       Vector nodal_volumes(_NunknownNodes);
 961       int j_int;
 962       for(int i = 0; i< _Nmailles ; i++)//On parcourt les cellules du maillage
 963       {
 964         Cell Ci = _mesh.getCell(i);
 965         for(int j = 0 ; j<_Ndim+1 ; j++)//On parcourt les noeuds de la cellule
 966         {
 967             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
 968                         {
 969                 j_int=DiffusionEquation::unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 970                 nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
 971             }
 972         }
 973        }
 974       for( j_int = 0; j_int< _NunknownNodes ; j_int++)
 975         nodal_volumes[j_int]=1/nodal_volumes[j_int];
 976       A.leftDiagonalScale(nodal_volumes);
 977   }
 978
 979   return A.getEigenvalues( nev, which, tol);
 980 }
 981 std::vector< Vector >
 982 StationaryDiffusionEquation::getEigenvectors(int nev, EPSWhich which, double tol) const
 983 {
 984   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 985   return A.getEigenvectors( nev, which, tol);
 986 }
 987 Field
 988 StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, double tol) const
 989 {
 990   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 991   MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
 992   Field my_eigenfield;
 993
 994     my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
 995
 996   my_eigenfield.setFieldByDataArrayDouble(d);
 997
 998   return my_eigenfield;
 999 }
1000
1001 Field&
1002 StationaryDiffusionEquation::getOutputTemperatureField()
1003 {
1004     if(!_computationCompletedSuccessfully)
1005         throw("Computation not performed yet or failed. No temperature field available");
1006     else
1007         return _VV;
1008 }
1009
1010 Field&
1011 StationaryDiffusionEquation::getRodTemperatureField()
1012 {
1013    return getOutputTemperatureField();
1014 }
1015
1016 vector<string>
1017 StationaryDiffusionEquation::getInputFieldsNames()
1018 {
1019         vector<string> result(2);
1020
1021         result[0]="FluidTemperature";
1022         result[1]="HeatPower";
1023
1024         return result;
1025 }
1026 vector<string>
1027 StationaryDiffusionEquation::getOutputFieldsNames()
1028 {
1029         vector<string> result(1);
1030
1031         result[0]="RodTemperature";
1032
1033         return result;
1034 }
1035
1036 Field&
1037 StationaryDiffusionEquation::getOutputField(const string& nameField )
1038 {
1039         if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
1040                 return getRodTemperatureField();
1041     else
1042     {
1043         cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
1044         throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
1045     }
1046 }
1047
1048 void
1049 StationaryDiffusionEquation::setInputField(const string& nameField, Field& inputField )
1050 {
1051         if(!_meshSet)
1052                 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setInputField set the mesh first");
1053
1054         if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
1055                 return setFluidTemperatureField( inputField) ;
1056         else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
1057                 return setHeatPowerField( inputField );
1058         else
1059     {
1060         cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
1061         throw CdmathException("StationaryDiffusionEquation::setInputField error : Unknown Field name");
1062     }
1063 }
1064
1065 void
1066 StationaryDiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){
1067         if(!_meshSet)
1068                 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setFluidTemperatureField set initial field first");
1069
1070         coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
1071         _fluidTemperatureField=coupledTemperatureField;
1072         _fluidTemperatureFieldSet=true;
1073 };
1074
1075 void
1076 StationaryDiffusionEquation::setHeatPowerField(Field heatPower){
1077         if(!_meshSet)
1078                 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
1079
1080         heatPower.getMesh().checkFastEquivalWith(_mesh);
1081         _heatPowerField=heatPower;
1082         _heatPowerFieldSet=true;
1083 }
1084
1085 void
1086 StationaryDiffusionEquation::setHeatPowerField(string fileName, string fieldName, int iteration, int order, int meshLevel){
1087         if(!_meshSet)
1088                 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
1089
1090         _heatPowerField=Field(fileName, CELLS,fieldName, iteration, order, meshLevel);
1091         _heatPowerField.getMesh().checkFastEquivalWith(_mesh);
1092         _heatPowerFieldSet=true;
1093 }