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