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         //Linear solver
 270         KSPCreate(PETSC_COMM_WORLD, &_ksp);
 271         KSPSetType(_ksp, _ksptype);
 272         KSPGetPC(_ksp, &_pc);
 273         if(_mpi_size==1 )
 274                 PCSetType(_pc, _pctype);
 275         else
 276         {
 277                 PCSetType(_pc, PCBJACOBI);//Global preconditioner is block jacobi
 278                 if(_pctype != (char*)&PCILU)//Default pc type is ilu
 279                 {
 280                         PetscOptionsSetValue(NULL,"-sub_pc_type ",_pctype);
 281                         PetscOptionsSetValue(NULL,"-sub_ksp_type ","preonly");
 282                         //If the above setvalue does not work, try the following
 283                         /*
 284                         KSPSetUp(_ksp);//to set the block Jacobi data structures (including creation of an internal KSP context for each block)
 285                         KSP * subKSP;
 286                         PC subpc;
 287                         int nlocal;//nb local blocs (should equal 1)
 288                         PCBJacobiGetSubKSP(_pc,&nlocal,NULL,&subKSP);
 289                         if(nlocal==1)
 290                         {
 291                                 KSPSetType(subKSP[0], KSPPREONLY);//local block solver is same as global
 292                                 KSPGetPC(subKSP[0],&subpc);
 293                                 PCSetType(subpc,_pctype);
 294                         }
 295                         else
 296                                 throw CdmathException("PC Block Jacobi, more than one block in this processor!!");
 297                         */
 298                 }
 299         }
 300         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 301
 302         _initializedMemory=true;
 303         save();//save initial data
 304 }
 305
 306 double DiffusionEquation::computeTimeStep(bool & stop){
 307         if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
 308                 _dt_diffusion=computeDiffusionMatrix(stop);
 309
 310     //reset right hand side
 311         VecCopy(_b0,_b);
 312
 313         _dt_src=computeRHS(stop);
 314
 315         VecAssemblyBegin(_b);
 316         VecAssemblyEnd(  _b);
 317
 318     if(_verbose or _system)
 319         {
 320                 PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
 321         VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
 322         }
 323
 324         stop=false;
 325         return min(_dt_diffusion,_dt_src);
 326 }
 327
 328 double DiffusionEquation::computeDiffusionMatrix(bool & stop)
 329 {
 330     double result;
 331
 332         MatZeroEntries(_A);
 333         VecZeroEntries(_b0);
 334
 335     if(_FECalculation)
 336         result=computeDiffusionMatrixFE(stop);
 337     else
 338         result=computeDiffusionMatrixFV(stop);
 339
 340         PetscPrintf(PETSC_COMM_WORLD,"Maximum diffusivity is %.2e, CFL = %.2f, Delta x = %.2e\n",_maxvp,_cfl,_minl);
 341
 342     MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 343         MatAssemblyEnd(  _A, MAT_FINAL_ASSEMBLY);
 344         VecAssemblyBegin(_b0);
 345         VecAssemblyEnd(  _b0);
 346
 347
 348     //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
 349     //update value here if variable  heat transfer coefficient
 350     if(_timeScheme == Implicit and _heatTransfertCoeff/(_rho*_cp)>_precision)
 351         MatShift(_A,_heatTransfertCoeff/(_rho*_cp));//Contribution from the liquit/solid heat transfer
 352
 353     if(_verbose or _system)
 354         MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
 355
 356     return  result;
 357 }
 358
 359 double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){
 360
 361         if(_mpi_rank == 0)
 362         {
 363                 Cell Cj;
 364                 string nameOfGroup;
 365                 double coeff;//Diffusion coefficients between nodes i and j
 366
 367             Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
 368             std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
 369             std::vector< int > nodeIds(_Ndim+1);//cell node Ids
 370             std::vector< Node > nodes(_Ndim+1);//cell nodes
 371             int i_int, j_int; //index of nodes j and k considered as unknown nodes
 372             bool dirichletCell_treated;
 373
 374             std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
 375             for (int idim=0; idim<_Ndim+1;idim++)
 376                 values[idim][idim]=1;
 377
 378             /* parameters for boundary treatment */
 379             vector< double > valuesBorder(_Ndim+1);
 380             Vector GradShapeFuncBorder(_Ndim+1);
 381
 382                 for (int j=0; j<_Nmailles;j++)
 383             {
 384                         Cj = _mesh.getCell(j);
 385
 386                 for (int idim=0; idim<_Ndim+1;idim++){
 387                     nodeIds[idim]=Cj.getNodeId(idim);
 388                     nodes[idim]=_mesh.getNode(nodeIds[idim]);
 389                     for (int jdim=0; jdim<_Ndim;jdim++)
 390                         M(idim,jdim)=nodes[idim].getPoint()[jdim];
 391                     M(idim,_Ndim)=1;
 392                 }
 393                 for (int idim=0; idim<_Ndim+1;idim++)
 394                     GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
 395
 396                 /* Loop on the edges of the cell */
 397                 for (int idim=0; idim<_Ndim+1;idim++)
 398                 {
 399                     if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
 400                     {//First node of the edge is not Dirichlet node
 401                         i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 402                         dirichletCell_treated=false;
 403                         for (int jdim=0; jdim<_Ndim+1;jdim++)
 404                         {
 405                             if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
 406                             {//Second node of the edge is not Dirichlet node
 407                                 j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
 408                                 MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
 409                             }
 410                             else if (!dirichletCell_treated)
 411                             {//Second node of the edge is a Dirichlet node
 412                                 dirichletCell_treated=true;
 413                                 for (int kdim=0; kdim<_Ndim+1;kdim++)
 414                                 {
 415                                                                 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
 416                                                                 if( it != _dirichletBoundaryValues.end() )
 417                                     {
 418                                         if( _dirichletValuesSet )
 419                                             valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
 420                                         else
 421                                             valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
 422                                     }
 423                                     else
 424                                         valuesBorder[kdim]=0;
 425                                 }
 426                                 GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
 427                                 coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
 428                                 VecSetValue(_b0,i_int,coeff, ADD_VALUES);
 429                             }
 430                         }
 431                     }
 432                 }
 433                 }
 434
 435             //Calcul de la contribution de la condition limite de Neumann au second membre
 436             if( _NdirichletNodes !=_NboundaryNodes)
 437             {
 438                 vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
 439                 int NboundaryFaces=boundaryFaces.size();
 440                 for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
 441                 {
 442                     Face Fi = _mesh.getFace(boundaryFaces[i]);
 443                     for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
 444                     {
 445                         if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
 446                         {
 447                             j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
 448                             if( _neumannValuesSet )
 449                                 coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
 450                             else
 451                                 coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
 452                             VecSetValue(_b0, j_int, coeff, ADD_VALUES);
 453                         }
 454                     }
 455                 }
 456             }
 457         }
 458
 459         _diffusionMatrixSet=true;
 460
 461         _maxvp=_diffusivity;//To do : optimise value with the mesh while respecting stability
 462         if(fabs(_maxvp)<_precision)
 463                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFE(): Error computing time step ! Maximum diffusivity is zero => division by zero");
 464         else
 465                 return _cfl*_minl*_minl/_maxvp;
 466 }
 467
 468 double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){
 469         if(_mpi_rank == 0)
 470         {
 471                 long nbFaces = _mesh.getNumberOfFaces();
 472                 Face Fj;
 473                 Cell Cell1,Cell2;
 474                 string nameOfGroup;
 475                 double inv_dxi, inv_dxj;
 476                 double barycenterDistance;
 477                 Vector normale(_Ndim);
 478                 double dn;
 479                 PetscInt idm, idn;
 480                 std::vector< int > idCells;
 481
 482                 for (int j=0; j<nbFaces;j++){
 483                         Fj = _mesh.getFace(j);
 484
 485                         // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
 486                         idCells = Fj.getCellsId();
 487                         Cell1 = _mesh.getCell(idCells[0]);
 488                         idm = idCells[0];
 489                 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
 490                     if (j == Cell1.getFacesId()[l]){
 491                         for (int idim = 0; idim < _Ndim; ++idim)
 492                             normale[idim] = Cell1.getNormalVector(l,idim);
 493                         break;
 494                     }
 495                 }
 496
 497                         //Compute velocity at the face Fj
 498                         dn=(_DiffusionTensor*normale)*normale;
 499                         if(fabs(dn)>_maxvp)
 500                                 _maxvp=fabs(dn);
 501
 502                         // compute 1/dxi = volume of Ci/area of Fj
 503                 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
 504
 505                         // If Fj is on the boundary
 506                         if (Fj.getNumberOfCells()==1) {
 507                                 if(_verbose )
 508                                 {
 509                                         cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
 510                                         for(int p=0; p<_Ndim; p++)
 511                                                 cout << normale[p] << ",";
 512                                         cout << ") "<<endl;
 513                                 }
 514                                 nameOfGroup = Fj.getGroupName();
 515
 516                                 if (_limitField[nameOfGroup].bcType==NeumannDiffusion){
 517                         VecSetValue(_b0,idm,   -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
 518                                 }
 519                                 else if(_limitField[nameOfGroup].bcType==DirichletDiffusion){
 520                                         barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
 521                                         MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance                           , ADD_VALUES);
 522                                         VecSetValue(_b0,idm,    dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
 523                                 }
 524                                 else {
 525                         stop=true ;
 526                                         cout<<"!!!!!!!!!!!!!!!!! Error DiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
 527                         cout<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 528                                         cout<<"Accepted boundary conditions are NeumannDiffusion "<<NeumannDiffusion<< " and DirichletDiffusion "<<DirichletDiffusion<<endl;
 529                         *_runLogFile<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 530                                         throw CdmathException("Boundary condition not accepted");
 531                                 }
 532
 533                                 // if Fj is inside the domain
 534                         } else  if (Fj.getNumberOfCells()==2 ){
 535                                 if(_verbose )
 536                                 {
 537                                         cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 538                                         cout << " ; vecteur normal=(";
 539                                         for(int p=0; p<_Ndim; p++)
 540                                                 cout << normale[p] << ",";
 541                                         cout << ") "<<endl;
 542                                 }
 543                                 Cell2 = _mesh.getCell(idCells[1]);
 544                                 idn = idCells[1];
 545                                 if (_Ndim > 1)
 546                                         inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
 547                                 else
 548                                         inv_dxj = 1/Cell2.getMeasure();
 549
 550                                 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
 551
 552                                 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
 553                                 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
 554                                 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
 555                                 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
 556                         }
 557                         else
 558                                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
 559                 }
 560         }
 561
 562         _diffusionMatrixSet=true;
 563
 564         MPI_Bcast(&_maxvp, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);//The determination of _maxvp is optimal, unlike in the FE case
 565
 566         if(fabs(_maxvp)<_precision)
 567                 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): Error computing time step ! Maximum diffusivity is zero => division by zero");
 568         else
 569                 return _cfl*_minl*_minl/_maxvp;
 570 }
 571
 572 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
 573
 574         if(_mpi_rank == 0)
 575         {
 576             double Ti;
 577             if(!_FECalculation)
 578                 for (int i=0; i<_Nmailles;i++)
 579                     //Contribution due to fluid/solide heat exchange + Contribution of the volumic heat power
 580                     if(_timeScheme == Explicit)
 581                     {
 582                         VecGetValues(_Tn, 1, &i, &Ti);
 583                         VecSetValue(_b,i,(_heatTransfertCoeff*(_fluidTemperatureField(i)-Ti)+_heatPowerField(i))/(_rho*_cp),ADD_VALUES);
 584                     }
 585                     else//Implicit scheme
 586                         VecSetValue(_b,i,(_heatTransfertCoeff* _fluidTemperatureField(i)    +_heatPowerField(i))/(_rho*_cp)    ,ADD_VALUES);
 587             else
 588                 {
 589                     Cell Ci;
 590                     std::vector< int > nodesId;
 591                     for (int i=0; i<_Nmailles;i++)
 592                     {
 593                         Ci=_mesh.getCell(i);
 594                         nodesId=Ci.getNodesId();
 595                         for (int j=0; j<nodesId.size();j++)
 596                             if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
 597                             {
 598                                 double coeff = (_heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]))/(_rho*_cp);
 599                                 VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
 600                             }
 601                     }
 602                 }
 603         }
 604
 605     if(_heatTransfertCoeff>_precision)
 606         return _rho*_cp/_heatTransfertCoeff;
 607     else
 608         return INFINITY;
 609 }
 610
 611 bool DiffusionEquation::initTimeStep(double dt){
 612
 613         /* tricky because of code coupling */
 614     if(_dt>0 and dt>0)//Previous time step was set and used
 615     {
 616         //Remove the contribution from dt to prepare for new time step. The diffusion matrix is not recomputed
 617         if(_timeScheme == Implicit)
 618             MatShift(_A,-1/_dt+1/dt);
 619         //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
 620     }
 621     else if(dt>0)//_dt==0, first time step
 622     {
 623         if(_timeScheme == Implicit)
 624             MatShift(_A,1/dt);
 625     }
 626     else//dt<=0
 627     {
 628         PetscPrintf(PETSC_COMM_WORLD,"DiffusionEquation::initTimeStep %.2e = \n",dt);
 629         throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero");
 630     }
 631     //At this stage _b contains _b0 + power + heat exchange
 632     VecAXPY(_b, 1/dt, _Tn);
 633
 634         _dt = dt;
 635
 636         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 637         {
 638                 PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n");
 639                 MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
 640         }
 641
 642         return _dt>0;
 643 }
 644
 645 void DiffusionEquation::abortTimeStep(){
 646     //Remove contribution of dt to the RHS
 647         VecAXPY(_b,  -1/_dt, _Tn);
 648     //Remove contribution of dt to the matrix
 649         if(_timeScheme == Implicit)
 650                 MatShift(_A,-1/_dt);
 651         _dt = 0;
 652 }
 653
 654 bool DiffusionEquation::iterateTimeStep(bool &converged)
 655 {
 656         bool stop=false;
 657
 658         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 659                 _maxvp=0;
 660                 computeTimeStep(stop);
 661         }
 662         if(stop){
 663                 converged=false;
 664                 return false;
 665         }
 666
 667         if(_timeScheme == Explicit)
 668         {
 669                 MatMult(_A, _Tn, _Tk);
 670                 VecAXPY(_Tk, -1, _b);
 671                 VecScale(_Tk, -_dt);
 672
 673                 converged = true;
 674         }
 675         else
 676         {
 677
 678 #if PETSC_VERSION_GREATER_3_5
 679                 KSPSetOperators(_ksp, _A, _A);
 680 #else
 681                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 682 #endif
 683
 684                 if(_conditionNumber)
 685                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 686                 KSPSolve(_ksp, _b, _Tk);
 687
 688                 KSPGetIterationNumber(_ksp, &_PetscIts);
 689                 if( _MaxIterLinearSolver < _PetscIts)
 690                         _MaxIterLinearSolver = _PetscIts;
 691                 if(_PetscIts>=_maxPetscIts)
 692                 {
 693                         PetscPrintf(PETSC_COMM_WORLD,"Systeme lineaire : pas de convergence de Petsc. Itérations maximales %d atteintes \n",_maxPetscIts);
 694                         *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 695                         converged=false;
 696                         return false;
 697                 }
 698                 else
 699                 {
 700                         VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
 701                         VecAXPY(_deltaT,  -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
 702                         VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
 703                         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
 704                 }
 705         }
 706
 707         VecCopy(_Tk, _Tkm1);
 708
 709         return true;
 710 }
 711 void DiffusionEquation::validateTimeStep()
 712 {
 713         VecCopy(_Tk, _deltaT);//ici Tk=Tnp1 donc on a deltaT=Tnp1
 714         VecAXPY(_deltaT,  -1, _Tn);//On obtient deltaT=Tnp1-Tn
 715         VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
 716
 717         _isStationary =(_erreur_rel <_precision);
 718
 719         VecCopy(_Tk, _Tn);
 720         VecCopy(_Tk, _Tkm1);
 721
 722         _time+=_dt;
 723         _nbTimeStep++;
 724
 725         if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 726                 save();
 727 }
 728
 729 void DiffusionEquation::save(){
 730     PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results at time step number %d \n\n", _nbTimeStep);
 731     *_runLogFile<< "Saving numerical results at time step number "<< _nbTimeStep << endl<<endl;
 732
 733         string resultFile(_path+"/DiffusionEquation");//Results
 734
 735         resultFile+="_";
 736         resultFile+=_fileName;
 737
 738         if(_mpi_size>1){
 739                 VecScatterBegin(_scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD);
 740                 VecScatterEnd(  _scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD);
 741         }
 742
 743     if(_verbose or _system)
 744         {
 745                 PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n");
 746         VecView(_Tn,PETSC_VIEWER_STDOUT_WORLD);
 747         }
 748
 749         if(_mpi_rank==0){
 750             //On remplit le champ
 751             double Ti;
 752             if(!_FECalculation)
 753                 for(int i =0; i<_Nmailles;i++)
 754                 {
 755                                 if(_mpi_size>1)
 756                                         VecGetValues(_Tn_seq, 1, &i, &Ti);
 757                                 else
 758                                         VecGetValues(_Tn    , 1, &i, &Ti);
 759
 760                     _VV(i)=Ti;
 761                 }
 762             else
 763             {
 764                 int globalIndex;
 765                 for(int i=0; i<_NunknownNodes; i++)
 766                 {
 767                                 if(_mpi_size>1)
 768                                         VecGetValues(_Tn_seq, 1, &i, &Ti);
 769                                 else
 770                                         VecGetValues(_Tk    , 1, &i, &Ti);
 771                     globalIndex = globalNodeIndex(i, _dirichletNodeIds);
 772                     _VV(globalIndex)=Ti;
 773                 }
 774
 775                 Node Ni;
 776                 string nameOfGroup;
 777                 for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
 778                 {
 779                     Ni=_mesh.getNode(_dirichletNodeIds[i]);
 780                     nameOfGroup = Ni.getGroupName();
 781                     _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
 782                 }
 783             }
 784                 _VV.setTime(_time,_nbTimeStep);
 785
 786                 // create mesh and component info
 787                 if (_nbTimeStep ==0 || _restartWithNewFileName){
 788                         if (_restartWithNewFileName)
 789                                 _restartWithNewFileName=false;
 790                         string suppress ="rm -rf "+resultFile+"_*";
 791                         system(suppress.c_str());//Nettoyage des précédents calculs identiques
 792
 793                 _VV.setInfoOnComponent(0,"Temperature_(K)");
 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                 else{   // do not create mesh
 808                         switch(_saveFormat)
 809                         {
 810                         case VTK :
 811                                 _VV.writeVTK(resultFile,false);
 812                                 break;
 813                         case MED :
 814                                 _VV.writeMED(resultFile,false);
 815                                 break;
 816                         case CSV :
 817                                 _VV.writeCSV(resultFile);
 818                                 break;
 819                         }
 820                 }
 821
 822             if(_isStationary)
 823                 {
 824                 resultFile+="_Stat";
 825                 switch(_saveFormat)
 826                 {
 827                 case VTK :
 828                     _VV.writeVTK(resultFile);
 829                     break;
 830                 case MED :
 831                     _VV.writeMED(resultFile);
 832                     break;
 833                 case CSV :
 834                     _VV.writeCSV(resultFile);
 835                     break;
 836                 }
 837             }
 838         }
 839 }
 840
 841 void DiffusionEquation::terminate(){
 842         VecDestroy(&_Tn);
 843         VecDestroy(&_Tk);
 844         VecDestroy(&_Tkm1);
 845         VecDestroy(&_deltaT);
 846         VecDestroy(&_b0);
 847         VecDestroy(&_b);
 848         MatDestroy(&_A);
 849         if(_mpi_size>1 && _mpi_rank == 0)
 850                 VecDestroy(&_Tn_seq);
 851 }
 852
 853 void
 854 DiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
 855 {
 856     _dirichletValuesSet=true;
 857     _dirichletBoundaryValues=dirichletBoundaryValues;
 858 }
 859
 860 void
 861 DiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
 862 {
 863     _neumannValuesSet=true;
 864     _neumannBoundaryValues=neumannBoundaryValues;
 865 }
 866
 867
 868 Field&
 869 DiffusionEquation::getOutputTemperatureField()
 870 {
 871     if(!_initializedMemory)
 872         throw("Computation not initialized. No temperature field available");
 873     else
 874         return _VV;
 875 }
 876
 877 Field&
 878 DiffusionEquation::getRodTemperatureField()
 879 {
 880    return getOutputTemperatureField();
 881 }
 882
 883 vector<string>
 884 DiffusionEquation::getInputFieldsNames()
 885 {
 886         vector<string> result(2);
 887
 888         result[0]="FluidTemperature";
 889         result[1]="HeatPower";
 890
 891         return result;
 892 }
 893 vector<string>
 894 DiffusionEquation::getOutputFieldsNames()
 895 {
 896         vector<string> result(1);
 897
 898         result[0]="RodTemperature";
 899
 900         return result;
 901 }
 902
 903 Field&
 904 DiffusionEquation::getOutputField(const string& nameField )
 905 {
 906         if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
 907                 return getRodTemperatureField();
 908     else
 909     {
 910         cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
 911         throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
 912     }
 913 }
 914
 915 void
 916 DiffusionEquation::setInputField(const string& nameField, Field& inputField )
 917 {
 918         if(!_initialDataSet)
 919                 throw CdmathException("!!!!!!!! DiffusionEquation::setInputField set initial field first");
 920
 921         if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
 922                 return setFluidTemperatureField( inputField) ;
 923         else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
 924                 return setHeatPowerField( inputField );
 925         else
 926     {
 927         cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
 928         throw CdmathException("DiffusionEquation::setInputField error : Unknown Field name");
 929     }
 930 }
 931
 932 void
 933 DiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){
 934         if(!_initialDataSet)
 935                 throw CdmathException("!!!!!!!! DiffusionEquation::setFluidTemperatureField set initial field first");
 936
 937         coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
 938         _fluidTemperatureField=coupledTemperatureField;
 939         _fluidTemperatureFieldSet=true;
 940 };