1 #include "DiffusionEquation.hxx"
9 int DiffusionEquation::fact(int n)
11 return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
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)
20 return globalIndex-boundarySize;
21 else if (dirichletNodes[j]>globalIndex)
24 throw CdmathException("DiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
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 */
32 return unknownNodeIndex;
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)
39 unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
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;
49 Vector DiffusionEquation::gradientNodal(Matrix M, vector< double > values){
50 vector< Matrix > matrices(_Ndim);
52 for (int idim=0; idim<_Ndim;idim++){
53 matrices[idim]=M.deepCopy();
54 for (int jdim=0; jdim<_Ndim+1;jdim++)
55 matrices[idim](jdim,idim) = values[jdim] ;
59 for (int idim=0; idim<_Ndim;idim++)
60 result[idim] = matrices[idim].determinant();
65 DiffusionEquation::DiffusionEquation(int dim, bool FECalculation,double rho,double cp, double lambda){
66 /* Control input value are acceptable */
67 if(rho<_precision or cp<_precision)
69 std::cout<<"rho="<<rho<<", cp= "<<cp<< ", precision= "<<_precision;
70 throw CdmathException("Error : parameters rho and cp should be strictly positive");
74 std::cout<<"conductivity="<<lambda<<endl;
75 throw CdmathException("Error : conductivity parameter lambda cannot be negative");
79 std::cout<<"space dimension="<<dim<<endl;
80 throw CdmathException("Error : parameter dim cannot be negative");
83 cout<<"Diffusion problem with density "<<rho<<", specific heat "<< cp<<", conductivity "<< lambda;
85 cout<<" and finite elements method"<<endl;
87 cout<<" and finite volumes method"<<endl;
89 _FECalculation=FECalculation;
91 /* Finite element data */
93 _boundaryNodeIds=std::vector< int >(0);
94 _dirichletNodeIds=std::vector< int >(0);
99 /* Physical parameters */
100 _conductivity=lambda;
103 _diffusivity=_conductivity/(_rho*_cp);
104 _fluidTemperatureFieldSet=false;
107 /* Numerical parameters */
112 _diffusionMatrixSet=false;
114 _fileName = "CoreFlowsDiffusionProblem";
116 _runLogFile=new ofstream;
118 /* Default diffusion tensor is identity matrix */
119 _DiffusionTensor=Matrix(_Ndim);
120 for(int idim=0;idim<_Ndim;idim++)
121 _DiffusionTensor(idim,idim)=1;
124 void DiffusionEquation::initialize()
126 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
128 if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim
130 cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<_mesh.getMeshDimension()<<", and space dimension is "<<_mesh.getSpaceDimension()<<endl;
131 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<_mesh.getMeshDimension()<<", mesh space dim= "<<_mesh.getSpaceDimension()<<endl;
132 *_runLogFile<<"DiffusionEquation::initialize: mesh has incorrect dimension"<<endl;
133 _runLogFile->close();
134 throw CdmathException("DiffusionEquation::initialize: mesh has incorrect dimension");
138 throw CdmathException("DiffusionEquation::initialize() set initial data first");
141 cout<<"Initialising the diffusion of a solid temperature using ";
142 *_runLogFile<<"Initialising the diffusion of a solid temperature using ";
145 cout<< "Finite volumes method"<<endl<<endl;
146 *_runLogFile<< "Finite volumes method"<<endl<<endl;
150 cout<< "Finite elements method"<<endl<<endl;
151 *_runLogFile<< "Finite elements method"<<endl<<endl;
155 /**************** Field creation *********************/
157 if(!_heatPowerFieldSet){
159 _heatPowerField=Field("Heat power",NODES,_mesh,1);
160 for(int i =0; i<_Nnodes; i++)
161 _heatPowerField(i) = _heatSource;
164 _heatPowerField=Field("Heat power",CELLS,_mesh,1);
165 for(int i =0; i<_Nmailles; i++)
166 _heatPowerField(i) = _heatSource;
168 _heatPowerFieldSet=true;
170 if(!_fluidTemperatureFieldSet){
172 _fluidTemperatureField=Field("Fluid temperature",NODES,_mesh,1);
173 for(int i =0; i<_Nnodes; i++)
174 _fluidTemperatureField(i) = _fluidTemperature;
177 _fluidTemperatureField=Field("Fluid temperature",CELLS,_mesh,1);
178 for(int i =0; i<_Nmailles; i++)
179 _fluidTemperatureField(i) = _fluidTemperature;
181 _fluidTemperatureFieldSet=true;
184 /* Détection des noeuds frontière avec une condition limite de Dirichlet */
187 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
188 cout<<"1D Finite element method on segments"<<endl;
191 if( _mesh.isTriangular() )//Mesh dim=2
192 cout<<"2D Finite element method on triangles"<<endl;
193 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
194 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
197 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
198 *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
199 _runLogFile->close();
200 throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
205 if( _mesh.isTetrahedral() )//Mesh dim=3
206 cout<<"3D Finite element method on tetrahedra"<<endl;
207 else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
208 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
209 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
210 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
213 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either tetrahedral, either a triangularised surface or 1D network"<<endl;
214 *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
215 _runLogFile->close();
216 throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
220 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
221 _NboundaryNodes=_boundaryNodeIds.size();
223 if(_NboundaryNodes==_Nnodes)
224 cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl;
226 for(int i=0; i<_NboundaryNodes; i++)
227 if(_limitField[(_mesh.getNode(_boundaryNodeIds[i])).getGroupName()].bcType==Dirichlet)
228 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
229 _NdirichletNodes=_dirichletNodeIds.size();
230 _NunknownNodes=_Nnodes - _NdirichletNodes;
231 cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
232 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
235 //creation de la matrice
237 MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
239 MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
241 VecCreate(PETSC_COMM_SELF, &_Tk);
244 VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
246 VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
248 VecSetFromOptions(_Tk);
249 VecDuplicate(_Tk, &_Tn);
250 VecDuplicate(_Tk, &_Tkm1);
251 VecDuplicate(_Tk, &_deltaT);
252 VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
253 VecDuplicate(_Tk, &_b0);//part of the RHS that comes from the boundary conditions. Computed only once at the first time step
256 for(int i =0; i<_Nmailles;i++)
257 VecSetValue(_Tn,i,_VV(i), INSERT_VALUES);
259 for(int i =0; i<_Nnodes;i++)
260 VecSetValue(_Tn,i,_VV(i), INSERT_VALUES);
263 KSPCreate(PETSC_COMM_SELF, &_ksp);
264 KSPSetType(_ksp, _ksptype);
265 // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
266 KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
267 KSPGetPC(_ksp, &_pc);
268 PCSetType(_pc, _pctype);
270 _initializedMemory=true;
271 save();//save initial data
274 double DiffusionEquation::computeTimeStep(bool & stop){
275 if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
276 _dt_diffusion=computeDiffusionMatrix(stop);
278 //reset right hand side
281 _dt_src=computeRHS(stop);
284 return min(_dt_diffusion,_dt_src);
287 double DiffusionEquation::computeDiffusionMatrix(bool & stop)
292 result=computeDiffusionMatrixFE(stop);
294 result=computeDiffusionMatrixFV(stop);
296 //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
297 //update value here if variable heat transfer coefficient
298 if(_timeScheme == Implicit and _heatTransfertCoeff/(_rho*_cp)>_precision)
299 MatShift(_A,_heatTransfertCoeff/(_rho*_cp));//Contribution from the liquit/solid heat transfer
301 if(_verbose or _system)
302 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
307 double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){
310 double dij;//Diffusion coefficients between nodes i and j
314 Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
315 std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
316 std::vector< int > nodeIds(_Ndim+1);//cell node Ids
317 std::vector< Node > nodes(_Ndim+1);//cell nodes
318 int i_int, j_int; //index of nodes j and k considered as unknown nodes
319 bool dirichletCell_treated;
321 std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
322 for (int idim=0; idim<_Ndim+1;idim++)
323 values[idim][idim]=1;
325 /* parameters for boundary treatment */
326 vector< double > valuesBorder(_Ndim+1);
327 Vector GradShapeFuncBorder(_Ndim+1);
329 for (int j=0; j<_Nmailles;j++)
331 Cj = _mesh.getCell(j);
333 for (int idim=0; idim<_Ndim+1;idim++){
334 nodeIds[idim]=Cj.getNodeId(idim);
335 nodes[idim]=_mesh.getNode(nodeIds[idim]);
336 for (int jdim=0; jdim<_Ndim;jdim++)
337 M(idim,jdim)=nodes[idim].getPoint()[jdim];
340 for (int idim=0; idim<_Ndim+1;idim++)
341 GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
343 /* Loop on the edges of the cell */
344 for (int idim=0; idim<_Ndim+1;idim++)
346 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
347 {//First node of the edge is not Dirichlet node
348 i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
349 dirichletCell_treated=false;
350 for (int jdim=0; jdim<_Ndim+1;jdim++)
352 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
353 {//Second node of the edge is not Dirichlet node
354 j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
355 dij=_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure();
356 MatSetValue(_A,i_int,j_int,dij, ADD_VALUES);
360 else if (!dirichletCell_treated)
361 {//Second node of the edge is a Dirichlet node
362 dirichletCell_treated=true;
363 for (int kdim=0; kdim<_Ndim+1;kdim++)
365 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[kdim])!=_dirichletNodeIds.end())
366 valuesBorder[kdim]=_limitField[nameOfGroup].T;
368 valuesBorder[kdim]=0;
370 GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
371 dij =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
372 VecSetValue(_b,i_int,dij, ADD_VALUES);
379 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
380 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
381 VecAssemblyBegin(_b);
384 _diffusionMatrixSet=true;
387 cout<<"_maxvp= "<<_maxvp<< " cfl= "<<_cfl<<" minl= "<<_minl<<endl;
388 if(fabs(_maxvp)<_precision)
389 throw CdmathException("DiffusionEquation::computeDiffusionMatrix(): maximum eigenvalue for time step is zero");
394 double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){
395 long nbFaces = _mesh.getNumberOfFaces();
399 double inv_dxi, inv_dxj;
400 double barycenterDistance;
401 Vector normale(_Ndim);
404 std::vector< int > idCells;
408 for (int j=0; j<nbFaces;j++){
409 Fj = _mesh.getFace(j);
411 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
412 idCells = Fj.getCellsId();
413 Cell1 = _mesh.getCell(idCells[0]);
415 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
416 if (j == Cell1.getFacesId()[l]){
417 for (int idim = 0; idim < _Ndim; ++idim)
418 normale[idim] = Cell1.getNormalVector(l,idim);
423 //Compute velocity at the face Fj
424 dn=_conductivity*(_DiffusionTensor*normale)*normale;
428 // compute 1/dxi = volume of Ci/area of Fj
429 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
431 // If Fj is on the boundary
432 if (Fj.getNumberOfCells()==1) {
435 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
436 for(int p=0; p<_Ndim; p++)
437 cout << normale[p] << ",";
440 nameOfGroup = Fj.getGroupName();
442 if (_limitField[nameOfGroup].bcType==Neumann){//Nothing to do
444 else if(_limitField[nameOfGroup].bcType==Dirichlet){
445 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
446 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
447 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
451 cout<<"!!!!!!!!!!!!!!!!! Error DiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
452 cout<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
453 cout<<"Accepted boundary conditions are Neumann "<<Neumann<< " and Dirichlet "<<Dirichlet<<endl;
454 *_runLogFile<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
455 throw CdmathException("Boundary condition not accepted");
458 // if Fj is inside the domain
459 } else if (Fj.getNumberOfCells()==2 ){
462 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
463 cout << " ; vecteur normal=(";
464 for(int p=0; p<_Ndim; p++)
465 cout << normale[p] << ",";
468 Cell2 = _mesh.getCell(idCells[1]);
471 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
473 inv_dxj = 1/Cell2.getMeasure();
475 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
477 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
478 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
479 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
480 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
483 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
486 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
487 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
488 VecAssemblyBegin(_b0);
491 _diffusionMatrixSet=true;
494 cout<<"_maxvp= "<<_maxvp<< " cfl= "<<_cfl<<" minl= "<<_minl<<endl;
495 if(fabs(_maxvp)<_precision)
496 throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): maximum eigenvalue for time step is zero");
498 return _cfl*_minl*_minl/_maxvp;
501 double DiffusionEquation::computeRHS(bool & stop){
502 VecAssemblyBegin(_b);
505 for (int i=0; i<_Nmailles;i++)
507 VecSetValue(_b,i,_heatPowerField(i)/(_rho*_cp),ADD_VALUES);//Contribution of the volumic heat power
508 //Contribution due to fluid/solide heat exchange
509 if(_timeScheme == Explicit)
511 VecGetValues(_Tn, 1, &i, &Ti);
512 VecSetValue(_b,i,_heatTransfertCoeff/(_rho*_cp)*(_fluidTemperatureField(i)-Ti),ADD_VALUES);
514 else//Implicit scheme
515 VecSetValue(_b,i,_heatTransfertCoeff/(_rho*_cp)* _fluidTemperatureField(i) ,ADD_VALUES);
520 std::vector< int > nodesId;
521 for (int i=0; i<_Nmailles;i++)
524 nodesId=Ci.getNodesId();
525 for (int j=0; j<nodesId.size();j++)
526 if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
528 double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
529 VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
535 if(_verbose or _system)
536 VecView(_b,PETSC_VIEWER_STDOUT_SELF);
539 if(_heatTransfertCoeff>_precision)
540 return _rho*_cp/_heatTransfertCoeff;
545 bool DiffusionEquation::initTimeStep(double dt){
549 //Remove the contribution from dt to prepare for new initTimeStep. The diffusion matrix is not recomputed
550 if(_timeScheme == Implicit)
551 MatShift(_A,-1/_dt+1/dt);
552 //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
554 else if(dt>0)//_dt==0, first time step
556 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
557 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
558 if(_timeScheme == Implicit)
563 cout<<"DiffusionEquation::initTimeStep dt= "<<dt<<endl;
564 throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero");
566 //At this stage _b contains _b0 + power + heat exchange
567 VecAXPY(_b, 1/_dt, _Tn);
571 if(_verbose && _nbTimeStep%_freqSave ==0)
572 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
577 void DiffusionEquation::abortTimeStep(){
578 //Remove contribution od dt to the RHS
579 VecAXPY(_b, -1/_dt, _Tn);
580 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
581 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
582 //Remove contribution od dt to the matrix
583 if(_timeScheme == Implicit)
588 bool DiffusionEquation::iterateTimeStep(bool &converged)
592 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
594 computeTimeStep(stop);
601 if(_timeScheme == Explicit)
603 MatMult(_A, _Tn, _Tk);
604 VecAXPY(_Tk, -1, _b);
611 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
612 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
614 #if PETSC_VERSION_GREATER_3_5
615 KSPSetOperators(_ksp, _A, _A);
617 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
621 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
622 KSPSolve(_ksp, _b, _Tk);
624 KSPGetIterationNumber(_ksp, &_PetscIts);
625 if( _MaxIterLinearSolver < _PetscIts)
626 _MaxIterLinearSolver = _PetscIts;
627 if(_PetscIts>=_maxPetscIts)
629 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
630 *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
635 VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
636 VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
641 for(int i=0; i<_Nmailles; i++)
643 VecGetValues(_deltaT, 1, &i, &dTi);
644 VecGetValues(_Tk, 1, &i, &Ti);
645 if(_erreur_rel < fabs(dTi/Ti))
646 _erreur_rel = fabs(dTi/Ti);
649 for(int i=0; i<_Nnodes; i++)
651 VecGetValues(_deltaT, 1, &i, &dTi);
652 VecGetValues(_Tk, 1, &i, &Ti);
653 if(_erreur_rel < fabs(dTi/Ti))
654 _erreur_rel = fabs(dTi/Ti);
656 converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
664 void DiffusionEquation::validateTimeStep()
666 VecCopy(_Tk, _deltaT);//ici Tk=Tnp1 donc on a deltaT=Tnp1
667 VecAXPY(_deltaT, -1, _Tn);//On obtient deltaT=Tnp1-Tn
673 for(int i=0; i<_Nmailles; i++)
675 VecGetValues(_deltaT, 1, &i, &dTi);
676 VecGetValues(_Tk, 1, &i, &Ti);
677 if(_erreur_rel < fabs(dTi/Ti))
678 _erreur_rel = fabs(dTi/Ti);
681 for(int i=0; i<_Nnodes; i++)
683 VecGetValues(_deltaT, 1, &i, &dTi);
684 VecGetValues(_Tk, 1, &i, &Ti);
685 if(_erreur_rel < fabs(dTi/Ti))
686 _erreur_rel = fabs(dTi/Ti);
689 _isStationary =(_erreur_rel <_precision);
694 if(_verbose && _nbTimeStep%_freqSave ==0)
695 cout <<"Valeur propre locale max: " << _maxvp << endl;
699 if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
703 void DiffusionEquation::save(){
704 cout<< "Saving numerical results"<<endl<<endl;
705 *_runLogFile<< "Saving numerical results"<< endl<<endl;
707 string resultFile(_path+"/DiffusionEquation");//Results
710 resultFile+=_fileName;
712 if(_verbose or _system)
713 VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
715 //On remplit le champ
718 for(int i =0; i<_Nmailles;i++)
720 VecGetValues(_Tn, 1, &i, &Ti);
726 for(int i=0; i<_NunknownNodes; i++)
728 VecGetValues(_Tk, 1, &i, &Ti);
729 globalIndex = globalNodeIndex(i, _dirichletNodeIds);
730 _VV(globalIndex)=Ti;//Assumes node numbering starts with border nodes
735 for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
737 Ni=_mesh.getNode(_dirichletNodeIds[i]);
738 nameOfGroup = Ni.getGroupName();
739 _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
742 _VV.setTime(_time,_nbTimeStep);
744 // create mesh and component info
745 if (_nbTimeStep ==0){
746 string suppress ="rm -rf "+resultFile+"_*";
747 system(suppress.c_str());//Nettoyage des précédents calculs identiques
749 _VV.setInfoOnComponent(0,"Temperature_(K)");
753 _VV.writeVTK(resultFile);
756 _VV.writeMED(resultFile);
759 _VV.writeCSV(resultFile);
763 else{ // do not create mesh
767 _VV.writeVTK(resultFile,false);
770 _VV.writeMED(resultFile,false);
773 _VV.writeCSV(resultFile);
784 _VV.writeVTK(resultFile);
787 _VV.writeMED(resultFile);
790 _VV.writeCSV(resultFile);
796 void DiffusionEquation::terminate(){
800 VecDestroy(&_deltaT);