CoreFlows/Models/src/DiffusionEquation.cxx

   1 #include "DiffusionEquation.hxx"
   2 #include "math.h"
   3 #include <algorithm>
   4 #include <fstream>
   5 #include <sstream>
   6
   7 using namespace std;
   8
   9 int DiffusionEquation::fact(int n)
  10 {
  11   return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
  12 }
  13 int DiffusionEquation::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes)
  14 {//assumes Dirichlet node numbering is strictly increasing
  15     int j=0;//indice de parcours des noeuds frontière avec CL Dirichlet
  16     int boundarySize=dirichletNodes.size();
  17     while(j<boundarySize and dirichletNodes[j]<globalIndex)
  18         j++;
  19     if(j==boundarySize)
  20         return globalIndex-boundarySize;
  21     else if (dirichletNodes[j]>globalIndex)
  22         return globalIndex-j;
  23     else
  24         throw CdmathException("DiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
  25 }
  26
  27 int DiffusionEquation::globalNodeIndex(int unknownNodeIndex, std::vector< int > dirichletNodes)
  28 {//assumes Dirichlet boundary node numbering is strictly increasing
  29     int boundarySize=dirichletNodes.size();
  30     /* trivial case where all boundary nodes are Neumann BC */
  31     if(boundarySize==0)
  32         return unknownNodeIndex;
  33
  34     double unknownNodeMax=-1;//max unknown node number in the interval between jth and (j+1)th Dirichlet boundary nodes
  35     int j=0;//indice de parcours des noeuds frontière
  36     //On cherche l'intervale [j,j+1] qui contient le noeud de numéro interieur unknownNodeIndex
  37     while(j+1<boundarySize and unknownNodeMax<unknownNodeIndex)
  38     {
  39         unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
  40         j++;
  41     }
  42
  43     if(j+1==boundarySize)
  44         return unknownNodeIndex+boundarySize;
  45     else //unknownNodeMax>=unknownNodeIndex, hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j]
  46         return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1;
  47 }
  48
  49 Vector DiffusionEquation::gradientNodal(Matrix M, vector< double > values)
  50 {
  51         if(! M.isSquare() )
  52                 throw CdmathException("DiffusionEquation::gradientNodal Matrix M should be square !!!");
  53
  54         int Ndim = M.getNumberOfRows()-1;
  55     vector< Matrix > matrices(Ndim);
  56
  57     for (int idim=0; idim<Ndim;idim++){
  58         matrices[idim]=M.deepCopy();
  59         for (int jdim=0; jdim<Ndim+1;jdim++)
  60                         matrices[idim](jdim,idim) = values[jdim] ;
  61     }
  62
  63         Vector result(Ndim);
  64     for (int idim=0; idim<Ndim;idim++)
  65         result[idim] = matrices[idim].determinant();
  66
  67         return result;
  68 }
  69
  70 DiffusionEquation::DiffusionEquation(int dim, bool FECalculation,double rho,double cp, double lambda, MPI_Comm comm ):ProblemCoreFlows(comm)
  71 {
  72     /* Control input value are acceptable */
  73     if(rho<_precision or cp<_precision)
  74     {
  75         PetscPrintf(PETSC_COMM_WORLD,"rho = %.2f, cp = %.2f, precision = %.2e\n",rho,cp,_precision);
  76         throw CdmathException("Error : parameters rho and cp should be strictly positive");
  77     }
  78     if(lambda < 0.)
  79     {
  80         PetscPrintf(PETSC_COMM_WORLD,"Conductivity = %.2f\n",lambda);
  81         throw CdmathException("Error : conductivity parameter lambda cannot  be negative");
  82     }
  83     if(dim<=0)
  84     {
  85         PetscPrintf(PETSC_COMM_WORLD,"Space dimension = %.2f\n",dim);
  86         throw CdmathException("Error : parameter dim cannot  be negative");
  87     }
  88
  89     PetscPrintf(PETSC_COMM_WORLD,"\n Diffusion problem with density %.2e, specific heat %.2e, conductivity %.2e", rho,cp,lambda);
  90     if(FECalculation)
  91         PetscPrintf(PETSC_COMM_WORLD," and finite elements method\n\n");
  92     else
  93         PetscPrintf(PETSC_COMM_WORLD," and finite volumes method\n\n");
  94
  95     _FECalculation=FECalculation;
  96
  97     /* Finite element data */
  98     _boundaryNodeIds=std::vector< int >(0);
  99     _dirichletNodeIds=std::vector< int >(0);
 100     _NboundaryNodes=0;
 101     _NdirichletNodes=0;
 102     _NunknownNodes=0;
 103     _dirichletValuesSet=false;
 104     _neumannValuesSet=false;
 105
 106     /* Physical parameters */
 107         _conductivity=lambda;
 108         _cp=cp;
 109         _rho=rho;
 110         _diffusivity=_conductivity/(_rho*_cp);
 111         _fluidTemperatureFieldSet=false;
 112         _fluidTemperature=0;
 113
 114     /* Numerical parameters */
 115         _Ndim=dim;
 116         _nVar=1;
 117         _dt_diffusion=0;
 118         _dt_src=0;
 119         _diffusionMatrixSet=false;
 120
 121         _fileName = "SolverlabDiffusionProblem";
 122
 123     /* Default diffusion tensor is diagonal */
 124         _DiffusionTensor=Matrix(_Ndim);
 125         for(int idim=0;idim<_Ndim;idim++)
 126                 _DiffusionTensor(idim,idim)=_diffusivity;
 127 }
 128
 129 void DiffusionEquation::initialize()
 130 {
 131         if(_mpi_rank==0)
 132         {
 133                 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 134
 135                 if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim
 136                 {
 137                 PetscPrintf(PETSC_COMM_SELF,"Problem : dimension defined is %d but mesh dimension= %d, and space dimension is %d",_Ndim,_mesh.getMeshDimension(),_mesh.getSpaceDimension());
 138                         *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<_mesh.getMeshDimension()<<", mesh space dim= "<<_mesh.getSpaceDimension()<<endl;
 139                         *_runLogFile<<"DiffusionEquation::initialize: mesh has incorrect dimension"<<endl;
 140                         _runLogFile->close();
 141                         throw CdmathException("!!!!!!!!DiffusionEquation::initialize: mesh has incorrect  dimension");
 142                 }
 143
 144                 if(!_initialDataSet)
 145                         throw CdmathException("!!!!!!!!DiffusionEquation::initialize() set initial data first");
 146                 else
 147                 {
 148                     PetscPrintf(PETSC_COMM_SELF,"\n Initialising the diffusion of a solid temperature using ");
 149                     *_runLogFile<<"Initialising the diffusion of a solid temperature using ";
 150                     if(!_FECalculation)
 151                     {
 152                         PetscPrintf(PETSC_COMM_SELF,"Finite volumes method\n\n");
 153                         *_runLogFile<< "Finite volumes method"<<endl<<endl;
 154                     }
 155                     else
 156                     {
 157                         PetscPrintf(PETSC_COMM_SELF,"Finite elements method\n\n");
 158                         *_runLogFile<< "Finite elements method"<<endl<<endl;
 159                     }
 160                 }
 161
 162                 /**************** Field creation *********************/
 163
 164                 if(!_heatPowerFieldSet){
 165                 _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
 166                 for(int i =0; i<_VV.getNumberOfElements(); i++)
 167                     _heatPowerField(i) = _heatSource;
 168                 _heatPowerFieldSet=true;
 169             }
 170                 if(!_fluidTemperatureFieldSet){
 171                         _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
 172                         for(int i =0; i<_VV.getNumberOfElements(); i++)
 173                                 _fluidTemperatureField(i) = _fluidTemperature;
 174                 _fluidTemperatureFieldSet=true;
 175                 }
 176
 177             /* Détection des noeuds frontière avec une condition limite de Dirichlet */
 178             if(_FECalculation)
 179             {
 180                 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
 181                                 PetscPrintf(PETSC_COMM_SELF,"1D Finite element method on segments\n");
 182                 else if(_Ndim==2)
 183                 {
 184                                 if( _mesh.isTriangular() )//Mesh dim=2
 185                                         PetscPrintf(PETSC_COMM_SELF,"2D Finite element method on triangles\n");
 186                                 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 187                                         PetscPrintf(PETSC_COMM_SELF,"1D Finite element method on a 2D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
 188                                 else
 189                                 {
 190                                         PetscPrintf(PETSC_COMM_SELF,"Error Finite element with space dimension %, and mesh dimension  %d, mesh should be either triangular either 1D network\n",_Ndim,_mesh.getMeshDimension());
 191                                         *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
 192                                         _runLogFile->close();
 193                                         throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
 194                                 }
 195                 }
 196                 else if(_Ndim==3)
 197                 {
 198                                 if( _mesh.isTetrahedral() )//Mesh dim=3
 199                                         PetscPrintf(PETSC_COMM_SELF,"3D Finite element method on tetrahedra\n");
 200                                 else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
 201                                         PetscPrintf(PETSC_COMM_SELF,"2D Finite element method on a 3D surface : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
 202                                 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
 203                                         PetscPrintf(PETSC_COMM_SELF,"1D Finite element method on a 3D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
 204                                 else
 205                                 {
 206                                         PetscPrintf(PETSC_COMM_SELF,"Error Finite element with space dimension %d, and mesh dimension  %d, mesh should be either tetrahedral, either a triangularised surface or 1D network",_Ndim,_mesh.getMeshDimension());
 207                                         *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
 208                                         _runLogFile->close();
 209                                         throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
 210                                 }
 211                 }
 212
 213                 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
 214                 _NboundaryNodes=_boundaryNodeIds.size();
 215
 216                 if(_NboundaryNodes==_Nnodes)
 217                     PetscPrintf(PETSC_COMM_SELF,"!!!!! Warning : all nodes are boundary nodes !!!!!");
 218
 219                 for(int i=0; i<_NboundaryNodes; i++)
 220                     if(_limitField[(_mesh.getNode(_boundaryNodeIds[i])).getGroupName()].bcType==DirichletDiffusion)
 221                         _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
 222                 _NdirichletNodes=_dirichletNodeIds.size();
 223                 _NunknownNodes=_Nnodes - _NdirichletNodes;
 224                 PetscPrintf(PETSC_COMM_SELF,"Number of unknown nodes %d, Number of boundary nodes %d, Number of Dirichlet boundary nodes %d\n\n", _NunknownNodes,_NboundaryNodes, _NdirichletNodes);
 225                 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
 226             }
 227         }
 228
 229     if(!_FECalculation)
 230                 _globalNbUnknowns = _Nmailles*_nVar;
 231     else{
 232         MPI_Bcast(&_NunknownNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 233                 _globalNbUnknowns = _NunknownNodes*_nVar;
 234         }
 235
 236         /* Vectors creations */
 237         VecCreate(PETSC_COMM_WORLD, &_Tk);//main unknown
 238     VecSetSizes(_Tk,PETSC_DECIDE,_globalNbUnknowns);
 239         VecSetFromOptions(_Tk);
 240         VecGetLocalSize(_Tk, &_localNbUnknowns);
 241
 242         VecDuplicate(_Tk, &_Tn);
 243         VecDuplicate(_Tk, &_Tkm1);
 244         VecDuplicate(_Tk, &_deltaT);
 245         VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
 246         VecDuplicate(_Tk, &_b0);//part of the RHS that comes from the boundary conditions. Computed only once at the first time step
 247
 248         if(_mpi_rank == 0)//Process 0 reads and distributes initial data
 249                 if(_FECalculation)
 250                         for(int i = 0; i<_NunknownNodes; i++)
 251                 {
 252                                 int globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 253                                 VecSetValue(_Tn,i,_VV(globalIndex), INSERT_VALUES);
 254                         }
 255                 else
 256                         for(int i = 0; i<_Nmailles; i++)
 257                                 VecSetValue( _Tn, i, _VV(i), INSERT_VALUES);
 258         VecAssemblyBegin(_Tn);
 259         VecAssemblyEnd(_Tn);
 260
 261         /* Matrix creation */
 262         MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
 263
 264         /* Local sequential vector creation */
 265         if(_mpi_size>1 && _mpi_rank == 0)
 266                 VecCreateSeq(PETSC_COMM_SELF, _globalNbUnknowns, &_Tn_seq);//For saving results on proc 0
 267         VecScatterCreateToZero(_Tn,&_scat,&_Tn_seq);
 268
 269         createKSP();
 270
 271         _initializedMemory=true;
 272         save();//save initial data
 273 }
 274
 275 double DiffusionEquation::computeTimeStep(bool & stop){
 276         if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
 277                 _dt_diffusion=computeDiffusionMatrix(stop);
 278
 279     //reset right hand side
 280         VecCopy(_b0,_b);
 281
 282         _dt_src=computeRHS(stop);
 283
 284         VecAssemblyBegin(_b);
 285         VecAssemblyEnd(  _b);
 286
 287     if(_verbose or _system)
 288         {
 289                 PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
 290         VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
 291         }
 292
 293         stop=false;
 294         return min(_dt_diffusion,_dt_src);
 295 }
 296
 297 double DiffusionEquation::computeDiffusionMatrix(bool & stop)
 298 {
 299     double result;
 300
 301         MatZeroEntries(_A);
 302         VecZeroEntries(_b0);
 303
 304     if(_FECalculation)
 305         result=computeDiffusionMatrixFE(stop);
 306     else
 307         result=computeDiffusionMatrixFV(stop);
 308
 309         PetscPrintf(PETSC_COMM_WORLD,"Maximum diffusivity is %.2e, CFL = %.2f, Delta x = %.2e\n",_maxvp,_cfl,_minl);
 310
 311     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 312         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 313         VecAssemblyBegin(_b0);
 314         VecAssemblyEnd(  _b0);
 315
 316
 317     //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
 318     //update value here if variable  heat transfer coefficient
 319     if(_timeScheme == Implicit and _heatTransfertCoeff/(_rho*_cp)>_precision)
 320         MatShift(_A,_heatTransfertCoeff/(_rho*_cp));//Contribution from the liquit/solid heat transfer
 321
 322     if(_verbose or _system)
 323         MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
 324
 325     return  result;
 326 }
 327
 328 double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){
 329
 330         if(_mpi_rank == 0)
 331         {
 332                 Cell Cj;
 333                 string nameOfGroup;
 334                 double coeff;//Diffusion coefficients between nodes i and j
 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]=gradientNodal(M,values[idim])/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=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= 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=gradientNodal(M,valuesBorder)/fact(_Ndim);
 396                                 coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 397                                 VecSetValue(_b0,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=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(_b0, j_int, coeff, ADD_VALUES);
 422                         }
 423                     }
 424                 }
 425             }
 426         }
 427
 428         _diffusionMatrixSet=true;
 429
 430         _maxvp=_diffusivity;//To do : optimise value with the mesh while respecting stability
 431         if(fabs(_maxvp)<_precision)
 432                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFE(): Error computing time step ! Maximum diffusivity is zero => division by zero");
 433         else
 434                 return _cfl*_minl*_minl/_maxvp;
 435 }
 436
 437 double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){
 438         if(_mpi_rank == 0)
 439         {
 440                 long nbFaces = _mesh.getNumberOfFaces();
 441                 Face Fj;
 442                 Cell Cell1,Cell2;
 443                 string nameOfGroup;
 444                 double inv_dxi, inv_dxj;
 445                 double barycenterDistance;
 446                 Vector normale(_Ndim);
 447                 double dn;
 448                 PetscInt idm, idn;
 449                 std::vector< int > idCells;
 450
 451                 for (int j=0; j<nbFaces;j++){
 452                         Fj = _mesh.getFace(j);
 453
 454                         // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 455                         idCells = Fj.getCellsId();
 456                         Cell1 = _mesh.getCell(idCells[0]);
 457                         idm = idCells[0];
 458                 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 459                     if (j == Cell1.getFacesId()[l]){
 460                         for (int idim = 0; idim < _Ndim; ++idim)
 461                             normale[idim] = Cell1.getNormalVector(l,idim);
 462                         break;
 463                     }
 464                 }
 465
 466                         //Compute velocity at the face Fj
 467                         dn=(_DiffusionTensor*normale)*normale;
 468                         if(fabs(dn)>_maxvp)
 469                                 _maxvp=fabs(dn);
 470
 471                         // compute 1/dxi = volume of Ci/area of Fj
 472                 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 473
 474                         // If Fj is on the boundary
 475                         if (Fj.getNumberOfCells()==1) {
 476                                 if(_verbose )
 477                                 {
 478                                         cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 479                                         for(int p=0; p<_Ndim; p++)
 480                                                 cout << normale[p] << ",";
 481                                         cout << ") "<<endl;
 482                                 }
 483                                 nameOfGroup = Fj.getGroupName();
 484
 485                                 if (_limitField[nameOfGroup].bcType==NeumannDiffusion){
 486                         VecSetValue(_b0,idm,   -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
 487                                 }
 488                                 else if(_limitField[nameOfGroup].bcType==DirichletDiffusion){
 489                                         barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 490                                         MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 491                                         VecSetValue(_b0,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 492                                 }
 493                                 else {
 494                         stop=true ;
 495                                         cout<<"!!!!!!!!!!!!!!!!! Error DiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 496                         cout<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 497                                         cout<<"Accepted boundary conditions are NeumannDiffusion "<<NeumannDiffusion<< " and DirichletDiffusion "<<DirichletDiffusion<<endl;
 498                         *_runLogFile<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 499                                         throw CdmathException("Boundary condition not accepted");
 500                                 }
 501
 502                                 // if Fj is inside the domain
 503                         } else  if (Fj.getNumberOfCells()==2 ){
 504                                 if(_verbose )
 505                                 {
 506                                         cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 507                                         cout << " ; vecteur normal=(";
 508                                         for(int p=0; p<_Ndim; p++)
 509                                                 cout << normale[p] << ",";
 510                                         cout << ") "<<endl;
 511                                 }
 512                                 Cell2 = _mesh.getCell(idCells[1]);
 513                                 idn = idCells[1];
 514                                 if (_Ndim > 1)
 515                                         inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 516                                 else
 517                                         inv_dxj = 1/Cell2.getMeasure();
 518
 519                                 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 520
 521                                 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 522                                 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 523                                 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 524                                 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 525                         }
 526                         else
 527                                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 528                 }
 529         }
 530
 531         _diffusionMatrixSet=true;
 532
 533         MPI_Bcast(&_maxvp, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);//The determination of _maxvp is optimal, unlike in the FE case
 534
 535         if(fabs(_maxvp)<_precision)
 536                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): Error computing time step ! Maximum diffusivity is zero => division by zero");
 537         else
 538                 return _cfl*_minl*_minl/_maxvp;
 539 }
 540
 541 double DiffusionEquation::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
 542
 543         if(_mpi_rank == 0)
 544         {
 545             double Ti;
 546             if(!_FECalculation)
 547                 for (int i=0; i<_Nmailles;i++)
 548                     //Contribution due to fluid/solide heat exchange + Contribution of the volumic heat power
 549                     if(_timeScheme == Explicit)
 550                     {
 551                         VecGetValues(_Tn, 1, &i, &Ti);
 552                         VecSetValue(_b,i,(_heatTransfertCoeff*(_fluidTemperatureField(i)-Ti)+_heatPowerField(i))/(_rho*_cp),ADD_VALUES);
 553                     }
 554                     else//Implicit scheme
 555                         VecSetValue(_b,i,(_heatTransfertCoeff* _fluidTemperatureField(i)    +_heatPowerField(i))/(_rho*_cp)    ,ADD_VALUES);
 556             else
 557                 {
 558                     Cell Ci;
 559                     std::vector< int > nodesId;
 560                     for (int i=0; i<_Nmailles;i++)
 561                     {
 562                         Ci=_mesh.getCell(i);
 563                         nodesId=Ci.getNodesId();
 564                         for (int j=0; j<nodesId.size();j++)
 565                             if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
 566                             {
 567                                 double coeff = (_heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]))/(_rho*_cp);
 568                                 VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 569                             }
 570                     }
 571                 }
 572         }
 573
 574     if(_heatTransfertCoeff>_precision)
 575         return _rho*_cp/_heatTransfertCoeff;
 576     else
 577         return INFINITY;
 578 }
 579
 580 bool DiffusionEquation::initTimeStep(double dt){
 581
 582         /* tricky because of code coupling */
 583     if(_dt>0 and dt>0)//Previous time step was set and used
 584     {
 585         //Remove the contribution from dt to prepare for new time step. The diffusion matrix is not recomputed
 586         if(_timeScheme == Implicit)
 587             MatShift(_A,-1/_dt+1/dt);
 588         //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
 589     }
 590     else if(dt>0)//_dt==0, first time step
 591     {
 592         if(_timeScheme == Implicit)
 593             MatShift(_A,1/dt);
 594     }
 595     else//dt<=0
 596     {
 597         PetscPrintf(PETSC_COMM_WORLD,"DiffusionEquation::initTimeStep %.2e = \n",dt);
 598         throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero");
 599     }
 600     //At this stage _b contains _b0 + power + heat exchange
 601     VecAXPY(_b, 1/dt, _Tn);
 602
 603         _dt = dt;
 604
 605         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 606         {
 607                 PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n");
 608                 MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
 609         }
 610
 611         return _dt>0;
 612 }
 613
 614 void DiffusionEquation::abortTimeStep(){
 615     //Remove contribution of dt to the RHS
 616         VecAXPY(_b,  -1/_dt, _Tn);
 617     //Remove contribution of dt to the matrix
 618         if(_timeScheme == Implicit)
 619                 MatShift(_A,-1/_dt);
 620         _dt = 0;
 621 }
 622
 623 bool DiffusionEquation::iterateTimeStep(bool &converged)
 624 {
 625         bool stop=false;
 626
 627         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 628                 _maxvp=0;
 629                 computeTimeStep(stop);
 630         }
 631         if(stop){
 632                 converged=false;
 633                 return false;
 634         }
 635
 636         if(_timeScheme == Explicit)
 637         {
 638                 MatMult(_A, _Tn, _Tk);
 639                 VecAXPY(_Tk, -1, _b);
 640                 VecScale(_Tk, -_dt);
 641
 642                 converged = true;
 643         }
 644         else
 645         {
 646
 647 #if PETSC_VERSION_GREATER_3_5
 648                 KSPSetOperators(_ksp, _A, _A);
 649 #else
 650                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 651 #endif
 652
 653                 if(_conditionNumber)
 654                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 655                 KSPSolve(_ksp, _b, _Tk);
 656
 657                 KSPGetIterationNumber(_ksp, &_PetscIts);
 658                 if( _MaxIterLinearSolver < _PetscIts)
 659                         _MaxIterLinearSolver = _PetscIts;
 660                 if(_PetscIts>=_maxPetscIts)
 661                 {
 662                         PetscPrintf(PETSC_COMM_WORLD,"Systeme lineaire : pas de convergence de Petsc. Itérations maximales %d atteintes \n",_maxPetscIts);
 663                         *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 664                         converged=false;
 665                         return false;
 666                 }
 667                 else
 668                 {
 669                         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 670                         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 671                         VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
 672                         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 673                 }
 674         }
 675
 676         VecCopy(_Tk, _Tkm1);
 677
 678         return true;
 679 }
 680 void DiffusionEquation::validateTimeStep()
 681 {
 682         VecCopy(_Tk, _deltaT);//ici Tk=Tnp1 donc on a deltaT=Tnp1
 683         VecAXPY(_deltaT,  -1, _Tn);//On obtient deltaT=Tnp1-Tn
 684         VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
 685
 686         _isStationary =(_erreur_rel <_precision);
 687
 688         VecCopy(_Tk, _Tn);
 689         VecCopy(_Tk, _Tkm1);
 690
 691         _time+=_dt;
 692         _nbTimeStep++;
 693
 694         if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 695                 save();
 696 }
 697
 698 void DiffusionEquation::save(){
 699     PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results at time step number %d \n\n", _nbTimeStep);
 700     *_runLogFile<< "Saving numerical results at time step number "<< _nbTimeStep << endl<<endl;
 701
 702         string resultFile(_path+"/DiffusionEquation");//Results
 703
 704         resultFile+="_";
 705         resultFile+=_fileName;
 706
 707         if(_mpi_size>1){
 708                 VecScatterBegin(_scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD);
 709                 VecScatterEnd(  _scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD);
 710         }
 711
 712     if(_verbose or _system)
 713         {
 714                 PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n");
 715         VecView(_Tn,PETSC_VIEWER_STDOUT_WORLD);
 716         }
 717
 718         if(_mpi_rank==0){
 719             //On remplit le champ
 720             double Ti;
 721             if(!_FECalculation)
 722                 for(int i =0; i<_Nmailles;i++)
 723                 {
 724                                 if(_mpi_size>1)
 725                                         VecGetValues(_Tn_seq, 1, &i, &Ti);
 726                                 else
 727                                         VecGetValues(_Tn    , 1, &i, &Ti);
 728
 729                     _VV(i)=Ti;
 730                 }
 731             else
 732             {
 733                 int globalIndex;
 734                 for(int i=0; i<_NunknownNodes; i++)
 735                 {
 736                                 if(_mpi_size>1)
 737                                         VecGetValues(_Tn_seq, 1, &i, &Ti);
 738                                 else
 739                                         VecGetValues(_Tk    , 1, &i, &Ti);
 740                     globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 741                     _VV(globalIndex)=Ti;
 742                 }
 743
 744                 Node Ni;
 745                 string nameOfGroup;
 746                 for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
 747                 {
 748                     Ni=_mesh.getNode(_dirichletNodeIds[i]);
 749                     nameOfGroup = Ni.getGroupName();
 750                     _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 751                 }
 752             }
 753                 _VV.setTime(_time,_nbTimeStep);
 754
 755                 // create mesh and component info
 756                 if (_nbTimeStep ==0 || _restartWithNewFileName){
 757                         if (_restartWithNewFileName)
 758                                 _restartWithNewFileName=false;
 759                         string suppress ="rm -rf "+resultFile+"_*";
 760                         system(suppress.c_str());//Nettoyage des précédents calculs identiques
 761
 762                 _VV.setInfoOnComponent(0,"Temperature_(K)");
 763                         switch(_saveFormat)
 764                         {
 765                         case VTK :
 766                                 _VV.writeVTK(resultFile);
 767                                 break;
 768                         case MED :
 769                                 _VV.writeMED(resultFile);
 770                                 break;
 771                         case CSV :
 772                                 _VV.writeCSV(resultFile);
 773                                 break;
 774                         }
 775                 }
 776                 else{   // do not create mesh
 777                         switch(_saveFormat)
 778                         {
 779                         case VTK :
 780                                 _VV.writeVTK(resultFile,false);
 781                                 break;
 782                         case MED :
 783                                 _VV.writeMED(resultFile,false);
 784                                 break;
 785                         case CSV :
 786                                 _VV.writeCSV(resultFile);
 787                                 break;
 788                         }
 789                 }
 790
 791             if(_isStationary)
 792                 {
 793                 resultFile+="_Stat";
 794                 switch(_saveFormat)
 795                 {
 796                 case VTK :
 797                     _VV.writeVTK(resultFile);
 798                     break;
 799                 case MED :
 800                     _VV.writeMED(resultFile);
 801                     break;
 802                 case CSV :
 803                     _VV.writeCSV(resultFile);
 804                     break;
 805                 }
 806             }
 807         }
 808 }
 809
 810 void DiffusionEquation::terminate(){
 811         VecDestroy(&_Tn);
 812         VecDestroy(&_Tk);
 813         VecDestroy(&_Tkm1);
 814         VecDestroy(&_deltaT);
 815         VecDestroy(&_b0);
 816         VecDestroy(&_b);
 817         MatDestroy(&_A);
 818         if(_mpi_size>1 && _mpi_rank == 0)
 819                 VecDestroy(&_Tn_seq);
 820 }
 821
 822 void
 823 DiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
 824 {
 825     _dirichletValuesSet=true;
 826     _dirichletBoundaryValues=dirichletBoundaryValues;
 827 }
 828
 829 void
 830 DiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
 831 {
 832     _neumannValuesSet=true;
 833     _neumannBoundaryValues=neumannBoundaryValues;
 834 }
 835
 836
 837 Field&
 838 DiffusionEquation::getOutputTemperatureField()
 839 {
 840     if(!_initializedMemory)
 841         throw("Computation not initialized. No temperature field available");
 842     else
 843         return _VV;
 844 }
 845
 846 Field&
 847 DiffusionEquation::getRodTemperatureField()
 848 {
 849    return getOutputTemperatureField();
 850 }
 851
 852 vector<string>
 853 DiffusionEquation::getInputFieldsNames()
 854 {
 855         vector<string> result(2);
 856
 857         result[0]="FluidTemperature";
 858         result[1]="HeatPower";
 859
 860         return result;
 861 }
 862 vector<string>
 863 DiffusionEquation::getOutputFieldsNames()
 864 {
 865         vector<string> result(1);
 866
 867         result[0]="RodTemperature";
 868
 869         return result;
 870 }
 871
 872 Field&
 873 DiffusionEquation::getOutputField(const string& nameField )
 874 {
 875         if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
 876                 return getRodTemperatureField();
 877     else
 878     {
 879         cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
 880         throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
 881     }
 882 }
 883
 884 void
 885 DiffusionEquation::setInputField(const string& nameField, Field& inputField )
 886 {
 887         if(!_initialDataSet)
 888                 throw CdmathException("!!!!!!!! DiffusionEquation::setInputField set initial field first");
 889
 890         if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
 891                 return setFluidTemperatureField( inputField) ;
 892         else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
 893                 return setHeatPowerField( inputField );
 894         else
 895     {
 896         cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
 897         throw CdmathException("DiffusionEquation::setInputField error : Unknown Field name");
 898     }
 899 }
 900
 901 void
 902 DiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){
 903         if(!_initialDataSet)
 904                 throw CdmathException("!!!!!!!! DiffusionEquation::setFluidTemperatureField set initial field first");
 905
 906         coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
 907         _fluidTemperatureField=coupledTemperatureField;
 908         _fluidTemperatureFieldSet=true;
 909 };
 910
 911 void
 912 DiffusionEquation::setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
 913         if(_FECalculation && field_support_type != NODES)
 914                 cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
 915         else if(!_FECalculation && field_support_type == NODES)
 916                 cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
 917
 918         Field VV;
 919
 920         switch(field_support_type)
 921         {
 922         case CELLS:
 923                 VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
 924                 break;
 925         case NODES:
 926                 VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
 927                 break;
 928         case FACES:
 929                 VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
 930                 break;
 931         default:
 932                 std::ostringstream message;
 933                 message << "Error DiffusionEquation::setDirichletBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
 934                 throw CdmathException(message.str().c_str());
 935         }
 936         /* For the moment the boundary value is taken constant equal to zero */
 937         _limitField[groupName]=LimitFieldDiffusion(DirichletDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
 938 }
 939
 940 void
 941 DiffusionEquation::setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
 942         if(_FECalculation && field_support_type != NODES)
 943                 cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
 944         else if(!_FECalculation && field_support_type == NODES)
 945                 cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
 946
 947         Field VV;
 948
 949         switch(field_support_type)
 950         {
 951         case CELLS:
 952                 VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
 953                 break;
 954         case NODES:
 955                 VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
 956                 break;
 957         case FACES:
 958                 VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
 959                 break;
 960         default:
 961                 std::ostringstream message;
 962                 message << "Error DiffusionEquation::setNeumannBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
 963                 throw CdmathException(message.str().c_str());
 964         }
 965         /* For the moment the boundary value is taken constant equal to zero */
 966         _limitField[groupName]=LimitFieldDiffusion(NeumannDiffusion,-1,0);//This line will be deleted when variable BC are properly treated in solverlab
 967 }