1 #include "StationaryDiffusionEquation.hxx"
3 #include "SparseMatrixPetsc.hxx"
11 int StationaryDiffusionEquation::fact(int n)
13 return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
15 int StationaryDiffusionEquation::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes) const
16 {//assumes Dirichlet node numbering is strictly increasing
17 int j=0;//indice de parcours des noeuds frontière avec CL Dirichlet
18 int boundarySize=dirichletNodes.size();
19 while(j<boundarySize and dirichletNodes[j]<globalIndex)
22 return globalIndex-boundarySize;
23 else if (dirichletNodes[j]>globalIndex)
26 throw CdmathException("StationaryDiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
29 int StationaryDiffusionEquation::globalNodeIndex(int unknownNodeIndex, std::vector< int > dirichletNodes) const
30 {//assumes Dirichlet boundary node numbering is strictly increasing
31 int boundarySize=dirichletNodes.size();
32 /* trivial case where all boundary nodes are Neumann BC */
34 return unknownNodeIndex;
36 double unknownNodeMax=-1;//max unknown node number in the interval between jth and (j+1)th Dirichlet boundary nodes
37 int j=0;//indice de parcours des noeuds frontière
38 //On cherche l'intervale [j,j+1] qui contient le noeud de numéro interieur unknownNodeIndex
39 while(j+1<boundarySize and unknownNodeMax<unknownNodeIndex)
41 unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
46 return unknownNodeIndex+boundarySize;
47 else //unknownNodeMax>=unknownNodeIndex, hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j]
48 return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1;
51 StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda){
52 PetscBool petscInitialized;
53 PetscInitialized(&petscInitialized);
55 PetscInitialize(NULL,NULL,0,0);
59 std::cout<<"Conductivity="<<lambda<<endl;
60 throw CdmathException("Error : conductivity parameter lambda cannot be negative");
64 std::cout<<"Space dimension="<<dim<<endl;
65 throw CdmathException("Error : parameter dim cannot be negative");
68 _FECalculation=FECalculation;
72 _nVar=1;//scalar prolem
74 _diffusionMatrixSet=false;
75 _initializedMemory=false;
83 _boundaryNodeIds=std::vector< int >(0);
84 _dirichletNodeIds=std::vector< int >(0);
88 _dirichletValuesSet=false;
89 _neumannValuesSet=false;
93 _precision_Newton=_precision;
94 _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
98 int _PetscIts=0;//the number of iterations of the linear solver
99 _ksptype = (char*)&KSPGMRES;
100 _pctype = (char*)&PCLU;
101 _conditionNumber=false;
104 //parameters for monitoring simulation
107 _runLogFile=new ofstream;
109 //result save parameters
110 _fileName = "StationaryDiffusionProblem";
111 char result[ PATH_MAX ];//extracting current directory
112 getcwd(result, PATH_MAX );
113 _path=string( result );
115 _computationCompletedSuccessfully=false;
117 //heat transfer parameters
118 _conductivity=lambda;
119 _fluidTemperatureFieldSet=false;
121 _heatPowerFieldSet=false;
122 _heatTransfertCoeff=0;
126 void StationaryDiffusionEquation::initialize()
128 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
131 throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first");
134 cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
135 *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
138 cout<< "Finite volumes method"<<endl<<endl;
139 *_runLogFile<< "Finite volumes method"<<endl<<endl;
143 cout<< "Finite elements method"<<endl<<endl;
144 *_runLogFile<< "Finite elements method"<<endl<<endl;
148 _DiffusionTensor=Matrix(_Ndim);
149 for(int idim=0;idim<_Ndim;idim++)
150 _DiffusionTensor(idim,idim)=1;
151 /**************** Field creation *********************/
153 if(!_heatPowerFieldSet){
155 _heatPowerField=Field("Heat power",NODES,_mesh,1);
156 for(int i =0; i<_Nnodes; i++)
157 _heatPowerField(i) = _heatSource;
160 _heatPowerField=Field("Heat power",CELLS,_mesh,1);
161 for(int i =0; i<_Nmailles; i++)
162 _heatPowerField(i) = _heatSource;
164 _heatPowerFieldSet=true;
166 if(!_fluidTemperatureFieldSet){
168 _fluidTemperatureField=Field("Fluid temperature",NODES,_mesh,1);
169 for(int i =0; i<_Nnodes; i++)
170 _fluidTemperatureField(i) = _fluidTemperature;
173 _fluidTemperatureField=Field("Fluid temperature",CELLS,_mesh,1);
174 for(int i =0; i<_Nmailles; i++)
175 _fluidTemperatureField(i) = _fluidTemperature;
177 _fluidTemperatureFieldSet=true;
180 /* Détection des noeuds frontière avec une condition limite de Dirichlet */
183 if(_NboundaryNodes==_Nnodes)
184 cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
186 for(int i=0; i<_NboundaryNodes; i++)
188 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
189 if( it != _dirichletBoundaryValues.end() )
190 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
191 else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
193 cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
194 *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
195 _runLogFile->close();
196 throw CdmathException("Missing boundary group");
198 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
200 cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
201 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
202 *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
203 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
204 _runLogFile->close();
205 throw CdmathException("Missing boundary condition");
207 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
208 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
209 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
211 cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
212 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
213 *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
214 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
215 _runLogFile->close();
216 throw CdmathException("Wrong boundary condition");
219 _NdirichletNodes=_dirichletNodeIds.size();
220 _NunknownNodes=_Nnodes - _NdirichletNodes;
221 cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
222 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
225 //creation de la matrice
227 MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
229 MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
231 VecCreate(PETSC_COMM_SELF, &_Tk);
234 VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
236 VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
238 VecSetFromOptions(_Tk);
239 VecDuplicate(_Tk, &_Tkm1);
240 VecDuplicate(_Tk, &_deltaT);
241 VecDuplicate(_Tk, &_b);//RHS of the linear system
244 KSPCreate(PETSC_COMM_SELF, &_ksp);
245 KSPSetType(_ksp, _ksptype);
246 // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
247 KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
248 KSPGetPC(_ksp, &_pc);
249 PCSetType(_pc, _pctype);
251 //Checking whether all boundary conditions are Neumann boundary condition
252 //if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
253 if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
255 map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
256 while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
258 _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
262 _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
264 _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
266 //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
269 std::cout<<"### Warning : all boundary conditions are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
270 std::cout<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
271 std::cout<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
272 *_runLogFile<<"### Warning : all boundary condition are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
273 *_runLogFile<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."<<std::endl<<endl;
274 *_runLogFile<<"### Check the compatibility condition between the right hand side and the boundary data. For homogeneous Neumann BCs, the right hand side must have integral equal to zero."<<std::endl;
276 //Check that the matrix is symmetric
277 PetscBool isSymetric;
278 MatIsSymmetric(_A,_precision,&isSymetric);
281 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
282 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
285 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
286 MatSetNullSpace(_A, nullsp);
287 MatSetTransposeNullSpace(_A, nullsp);
288 MatNullSpaceDestroy(&nullsp);
289 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
290 //PCFactorSetShiftAmount(_pc,1e-10);
293 _initializedMemory=true;
297 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
298 if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
299 computeDiffusionMatrix(stop);
301 _dt_src=computeRHS(stop);
306 Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values){
307 vector< Matrix > matrices(_Ndim);
309 for (int idim=0; idim<_Ndim;idim++){
310 matrices[idim]=M.deepCopy();
311 for (int jdim=0; jdim<_Ndim+1;jdim++)
312 matrices[idim](jdim,idim) = values[jdim] ;
315 Vector result(_Ndim);
316 for (int idim=0; idim<_Ndim;idim++)
317 result[idim] = matrices[idim].determinant();
322 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
327 result=computeDiffusionMatrixFE(stop);
329 result=computeDiffusionMatrixFV(stop);
331 //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
332 //update value here if variable heat transfer coefficient
333 if(_heatTransfertCoeff>_precision)
334 MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
336 if(_verbose or _system)
337 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
342 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
349 Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
350 std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
351 std::vector< int > nodeIds(_Ndim+1);//cell node Ids
352 std::vector< Node > nodes(_Ndim+1);//cell nodes
353 int i_int, j_int; //index of nodes j and k considered as unknown nodes
354 bool dirichletCell_treated;
356 std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
357 for (int idim=0; idim<_Ndim+1;idim++)
358 values[idim][idim]=1;
360 /* parameters for boundary treatment */
361 vector< double > valuesBorder(_Ndim+1);
362 Vector GradShapeFuncBorder(_Ndim+1);
364 for (int j=0; j<_Nmailles;j++)
366 Cj = _mesh.getCell(j);
368 for (int idim=0; idim<_Ndim+1;idim++){
369 nodeIds[idim]=Cj.getNodeId(idim);
370 nodes[idim]=_mesh.getNode(nodeIds[idim]);
371 for (int jdim=0; jdim<_Ndim;jdim++)
372 M(idim,jdim)=nodes[idim].getPoint()[jdim];
375 for (int idim=0; idim<_Ndim+1;idim++)
376 GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
378 /* Loop on the edges of the cell */
379 for (int idim=0; idim<_Ndim+1;idim++)
381 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
382 {//First node of the edge is not Dirichlet node
383 i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
384 dirichletCell_treated=false;
385 for (int jdim=0; jdim<_Ndim+1;jdim++)
387 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
388 {//Second node of the edge is not Dirichlet node
389 j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
390 MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
392 else if (!dirichletCell_treated)
393 {//Second node of the edge is a Dirichlet node
394 dirichletCell_treated=true;
395 for (int kdim=0; kdim<_Ndim+1;kdim++)
397 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
398 if( it != _dirichletBoundaryValues.end() )
400 if( _dirichletValuesSet )
401 valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
403 valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
406 valuesBorder[kdim]=0;
408 GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
409 coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
410 VecSetValue(_b,i_int,coeff, ADD_VALUES);
417 //Calcul de la contribution de la condition limite de Neumann au second membre
418 if( _NdirichletNodes !=_NboundaryNodes)
420 vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
421 int NboundaryFaces=boundaryFaces.size();
422 for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
424 Face Fi = _mesh.getFace(i);
425 for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
427 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
429 j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
430 if( _neumannValuesSet )
431 coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
433 coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
434 VecSetValue(_b, j_int, coeff, ADD_VALUES);
439 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
440 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
441 VecAssemblyBegin(_b);
444 _diffusionMatrixSet=true;
450 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
451 long nbFaces = _mesh.getNumberOfFaces();
455 double inv_dxi, inv_dxj;
456 double barycenterDistance;
457 Vector normale(_Ndim);
460 std::vector< int > idCells;
464 for (int j=0; j<nbFaces;j++){
465 Fj = _mesh.getFace(j);
467 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
468 idCells = Fj.getCellsId();
469 Cell1 = _mesh.getCell(idCells[0]);
471 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
472 if (j == Cell1.getFacesId()[l]){
473 for (int idim = 0; idim < _Ndim; ++idim)
474 normale[idim] = Cell1.getNormalVector(l,idim);
479 //Compute velocity at the face Fj
480 dn=_conductivity*(_DiffusionTensor*normale)*normale;
482 // compute 1/dxi = volume of Ci/area of Fj
483 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
485 // If Fj is on the boundary
486 if (Fj.getNumberOfCells()==1) {
489 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
490 for(int p=0; p<_Ndim; p++)
491 cout << normale[p] << ",";
495 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
496 if( it != _dirichletBoundaryValues.end() )
498 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
499 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
500 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
504 nameOfGroup = Fj.getGroupName();
506 if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
507 VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
509 else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
510 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
511 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
512 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
516 cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
517 cout<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
518 cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
519 *_runLogFile<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
520 _runLogFile->close();
521 throw CdmathException("Boundary condition not accepted");
524 // if Fj is inside the domain
525 } else if (Fj.getNumberOfCells()==2 ){
528 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
529 cout << " ; vecteur normal=(";
530 for(int p=0; p<_Ndim; p++)
531 cout << normale[p] << ",";
534 Cell2 = _mesh.getCell(idCells[1]);
537 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
539 inv_dxj = 1/Cell2.getMeasure();
541 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
543 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
544 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
545 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
546 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
550 *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
551 throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
555 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
556 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
557 VecAssemblyBegin(_b);
560 _diffusionMatrixSet=true;
566 double StationaryDiffusionEquation::computeRHS(bool & stop)//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeDiffusionMatrix)
568 VecAssemblyBegin(_b);
571 for (int i=0; i<_Nmailles;i++)
572 VecSetValue(_b,i,_heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i),ADD_VALUES);
576 std::vector< int > nodesId;
577 for (int i=0; i<_Nmailles;i++)
580 nodesId=Ci.getNodesId();
581 for (int j=0; j<nodesId.size();j++)
582 if(!_mesh.isBorderNode(nodesId[j]))
584 double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
585 VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
592 if(_verbose or _system)
593 VecView(_b,PETSC_VIEWER_STDOUT_SELF);
599 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
603 //Only implicit scheme considered
604 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
605 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
607 #if PETSC_VERSION_GREATER_3_5
608 KSPSetOperators(_ksp, _A, _A);
610 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
614 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
615 KSPSolve(_ksp, _b, _Tk);
617 KSPConvergedReason reason;
618 KSPGetConvergedReason(_ksp,&reason);
619 KSPGetIterationNumber(_ksp, &_PetscIts);
621 KSPGetResidualNorm(_ksp,&residu);
622 if (reason!=2 and reason!=3)
624 cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
625 cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
626 cout<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
627 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
628 *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
629 *_runLogFile<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
630 _runLogFile->close();
635 if( _MaxIterLinearSolver < _PetscIts)
636 _MaxIterLinearSolver = _PetscIts;
637 cout<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<< _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
638 *_runLogFile<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<< _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
639 VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
640 VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
645 VecAssemblyBegin(_Tk);
646 VecAssemblyEnd( _Tk);
647 VecAssemblyBegin(_deltaT);
648 VecAssemblyEnd( _deltaT);
651 for(int i=0; i<_Nmailles; i++)
653 VecGetValues(_deltaT, 1, &i, &dTi);
654 VecGetValues(_Tk, 1, &i, &Ti);
655 if(_erreur_rel < fabs(dTi/Ti))
656 _erreur_rel = fabs(dTi/Ti);
659 for(int i=0; i<_NunknownNodes; i++)
661 VecGetValues(_deltaT, 1, &i, &dTi);
662 VecGetValues(_Tk, 1, &i, &Ti);
663 if(_erreur_rel < fabs(dTi/Ti))
664 _erreur_rel = fabs(dTi/Ti);
667 converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
675 void StationaryDiffusionEquation::setMesh(const Mesh &M)
677 if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
679 cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
680 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
681 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
682 _runLogFile->close();
683 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
687 _Nmailles = _mesh.getNumberOfCells();
688 _Nnodes = _mesh.getNumberOfNodes();
690 cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
691 *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
693 // find maximum nb of neibourghs
696 _VV=Field ("Temperature", CELLS, _mesh, 1);
697 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
701 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
702 cout<<"1D Finite element method on segments"<<endl;
705 if( _mesh.isTriangular() )//Mesh dim=2
706 cout<<"2D Finite element method on triangles"<<endl;
707 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
708 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
711 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
712 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
713 _runLogFile->close();
714 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
719 if( _mesh.isTetrahedral() )//Mesh dim=3
720 cout<<"3D Finite element method on tetrahedra"<<endl;
721 else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
722 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
723 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
724 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
727 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;
728 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
729 _runLogFile->close();
730 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
734 _VV=Field ("Temperature", NODES, _mesh, 1);
736 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
737 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
738 _NboundaryNodes=_boundaryNodeIds.size();
744 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
746 //_maxPetscIts=maxIterationsPetsc;
747 // set linear solver algorithm
749 _ksptype = (char*)&KSPGMRES;
750 else if (kspType==CG)
751 _ksptype = (char*)&KSPCG;
752 else if (kspType==BCGS)
753 _ksptype = (char*)&KSPBCGS;
755 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
756 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
757 _runLogFile->close();
758 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
760 // set preconditioner
762 _pctype = (char*)&PCNONE;
763 else if (pcType ==LU)
764 _pctype = (char*)&PCLU;
765 else if (pcType == ILU)
766 _pctype = (char*)&PCILU;
767 else if (pcType ==CHOLESKY)
768 _pctype = (char*)&PCCHOLESKY;
769 else if (pcType == ICC)
770 _pctype = (char*)&PCICC;
772 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
773 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
774 _runLogFile->close();
775 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
779 bool StationaryDiffusionEquation::solveStationaryProblem()
781 if(!_initializedMemory)
783 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
784 _runLogFile->close();
785 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
787 bool stop=false; // Does the Problem want to stop (error) ?
788 bool converged=false; // has the newton scheme converged (end) ?
790 cout<< "!!! Running test case "<< _fileName << " using ";
791 *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
795 cout<< "Finite volumes method"<<endl<<endl;
796 *_runLogFile<< "Finite volumes method"<<endl<<endl;
800 cout<< "Finite elements method"<<endl<<endl;
801 *_runLogFile<< "Finite elements method"<< endl<<endl;
804 computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
806 cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
807 *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
808 _runLogFile->close();
809 throw CdmathException("Failed computing diffusion matrix");
811 computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
813 cout << "Error : failed computing right hand side, stopping calculation"<< endl;
814 *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
815 throw CdmathException("Failed computing right hand side");
817 stop = iterateNewtonStep(converged);
819 cout << "Error : failed solving linear system, stopping calculation"<< endl;
820 *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
821 _runLogFile->close();
822 throw CdmathException("Failed solving linear system");
825 _computationCompletedSuccessfully=true;
828 // Newton iteration loop for non linear problems
830 while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
832 computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
833 computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
834 stop = iterateNewtonStep(converged);
838 cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
839 *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
840 throw CdmathException("Failed solving a Newton iteration");
842 else if(_NEWTON_its==_maxNewtonIts){
843 cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
844 *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
845 throw CdmathException("No convergence of Newton scheme");
848 cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
849 *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
854 *_runLogFile<< "!!!!!! Computation successful !!!!!!"<< endl;
855 _runLogFile->close();
860 void StationaryDiffusionEquation::save(){
861 cout<< "Saving numerical results"<<endl<<endl;
862 *_runLogFile<< "Saving numerical results"<< endl<<endl;
864 string resultFile(_path+"/StationaryDiffusionEquation");//Results
867 resultFile+=_fileName;
869 // create mesh and component info
870 string suppress ="rm -rf "+resultFile+"_*";
871 system(suppress.c_str());//Nettoyage des précédents calculs identiques
873 if(_verbose or _system)
874 VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
878 for(int i=0; i<_Nmailles; i++)
880 VecGetValues(_Tk, 1, &i, &Ti);
886 for(int i=0; i<_NunknownNodes; i++)
888 VecGetValues(_Tk, 1, &i, &Ti);
889 globalIndex = globalNodeIndex(i, _dirichletNodeIds);
895 for(int i=0; i<_NdirichletNodes; i++)
897 Ni=_mesh.getNode(_dirichletNodeIds[i]);
898 nameOfGroup = Ni.getGroupName();
899 _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
903 _VV.setInfoOnComponent(0,"Temperature_(K)");
907 _VV.writeVTK(resultFile);
910 _VV.writeMED(resultFile);
913 _VV.writeCSV(resultFile);
918 StationaryDiffusionEquation::getOutputTemperatureField()
920 if(!_computationCompletedSuccessfully)
921 throw("Computation not performed yet or failed. No temperature field available");
926 void StationaryDiffusionEquation::terminate()
930 VecDestroy(&_deltaT);
935 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
937 _dirichletValuesSet=true;
938 _dirichletBoundaryValues=dirichletBoundaryValues;
942 StationaryDiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
944 _neumannValuesSet=true;
945 _neumannBoundaryValues=neumannBoundaryValues;
949 StationaryDiffusionEquation::getConditionNumber(bool isSingular, double tol) const
951 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
952 return A.getConditionNumber( isSingular, tol);
954 std::vector< double >
955 StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) const
957 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
959 if(_FECalculation)//We need to scale the FE matrix, otherwise the eigenvalues go to zero as the mesh is refined
961 Vector nodal_volumes(_NunknownNodes);
963 for(int i = 0; i< _Nmailles ; i++)//On parcourt les cellules du maillage
965 Cell Ci = _mesh.getCell(i);
966 for(int j = 0 ; j<_Ndim+1 ; j++)//On parcourt les noeuds de la cellule
968 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
970 j_int=unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
971 nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
975 for( j_int = 0; j_int< _NunknownNodes ; j_int++)
976 nodal_volumes[j_int]=1/nodal_volumes[j_int];
977 A.leftDiagonalScale(nodal_volumes);
980 return A.getEigenvalues( nev, which, tol);
982 std::vector< Vector >
983 StationaryDiffusionEquation::getEigenvectors(int nev, EPSWhich which, double tol) const
985 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
986 return A.getEigenvectors( nev, which, tol);
989 StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, double tol) const
991 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
992 MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
996 my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
998 my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
1000 my_eigenfield.setFieldByDataArrayDouble(d);
1002 return my_eigenfield;