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)
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)
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;
277 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
278 MatSetNullSpace(_A, nullsp);
279 MatSetTransposeNullSpace(_A, nullsp);
280 MatNullSpaceDestroy(&nullsp);
281 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
282 //PCFactorSetShiftAmount(_pc,1e-10);
284 _initializedMemory=true;
287 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
288 if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
289 computeDiffusionMatrix(stop);
291 _dt_src=computeRHS(stop);
296 Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values){
297 vector< Matrix > matrices(_Ndim);
299 for (int idim=0; idim<_Ndim;idim++){
300 matrices[idim]=M.deepCopy();
301 for (int jdim=0; jdim<_Ndim+1;jdim++)
302 matrices[idim](jdim,idim) = values[jdim] ;
305 Vector result(_Ndim);
306 for (int idim=0; idim<_Ndim;idim++)
307 result[idim] = matrices[idim].determinant();
312 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
317 result=computeDiffusionMatrixFE(stop);
319 result=computeDiffusionMatrixFV(stop);
321 //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
322 //update value here if variable heat transfer coefficient
323 if(_heatTransfertCoeff>_precision)
324 MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
326 if(_verbose or _system)
327 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
332 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
339 Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
340 std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
341 std::vector< int > nodeIds(_Ndim+1);//cell node Ids
342 std::vector< Node > nodes(_Ndim+1);//cell nodes
343 int i_int, j_int; //index of nodes j and k considered as unknown nodes
344 bool dirichletCell_treated;
346 std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
347 for (int idim=0; idim<_Ndim+1;idim++)
348 values[idim][idim]=1;
350 /* parameters for boundary treatment */
351 vector< double > valuesBorder(_Ndim+1);
352 Vector GradShapeFuncBorder(_Ndim+1);
354 for (int j=0; j<_Nmailles;j++)
356 Cj = _mesh.getCell(j);
358 for (int idim=0; idim<_Ndim+1;idim++){
359 nodeIds[idim]=Cj.getNodeId(idim);
360 nodes[idim]=_mesh.getNode(nodeIds[idim]);
361 for (int jdim=0; jdim<_Ndim;jdim++)
362 M(idim,jdim)=nodes[idim].getPoint()[jdim];
365 for (int idim=0; idim<_Ndim+1;idim++)
366 GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
368 /* Loop on the edges of the cell */
369 for (int idim=0; idim<_Ndim+1;idim++)
371 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
372 {//First node of the edge is not Dirichlet node
373 i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
374 dirichletCell_treated=false;
375 for (int jdim=0; jdim<_Ndim+1;jdim++)
377 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
378 {//Second node of the edge is not Dirichlet node
379 j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
380 MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
382 else if (!dirichletCell_treated)
383 {//Second node of the edge is a Dirichlet node
384 dirichletCell_treated=true;
385 for (int kdim=0; kdim<_Ndim+1;kdim++)
387 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
388 if( it != _dirichletBoundaryValues.end() )
390 if( _dirichletValuesSet )
391 valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
393 valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
396 valuesBorder[kdim]=0;
398 GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
399 coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
400 VecSetValue(_b,i_int,coeff, ADD_VALUES);
407 //Calcul de la contribution de la condition limite de Neumann au second membre
408 if( _NdirichletNodes !=_NboundaryNodes)
410 vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
411 int NboundaryFaces=boundaryFaces.size();
412 for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
414 Face Fi = _mesh.getFace(i);
415 for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
417 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
419 j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
420 if( _neumannValuesSet )
421 coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
423 coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
424 VecSetValue(_b, j_int, coeff, ADD_VALUES);
429 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
430 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
431 VecAssemblyBegin(_b);
434 if(_onlyNeumannBC) //Check that the matrix is symmetric
436 PetscBool isSymetric;
437 MatIsSymmetric(_A,_precision,&isSymetric);
440 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
441 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
445 _diffusionMatrixSet=true;
451 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
452 long nbFaces = _mesh.getNumberOfFaces();
456 double inv_dxi, inv_dxj;
457 double barycenterDistance;
458 Vector normale(_Ndim);
461 std::vector< int > idCells;
465 for (int j=0; j<nbFaces;j++){
466 Fj = _mesh.getFace(j);
468 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
469 idCells = Fj.getCellsId();
470 Cell1 = _mesh.getCell(idCells[0]);
472 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
473 if (j == Cell1.getFacesId()[l]){
474 for (int idim = 0; idim < _Ndim; ++idim)
475 normale[idim] = Cell1.getNormalVector(l,idim);
480 //Compute velocity at the face Fj
481 dn=_conductivity*(_DiffusionTensor*normale)*normale;
483 // compute 1/dxi = volume of Ci/area of Fj
484 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
486 // If Fj is on the boundary
487 if (Fj.getNumberOfCells()==1) {
490 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
491 for(int p=0; p<_Ndim; p++)
492 cout << normale[p] << ",";
496 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
497 if( it != _dirichletBoundaryValues.end() )
499 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
500 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
501 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
505 nameOfGroup = Fj.getGroupName();
507 if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
508 VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
510 else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
511 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
512 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
513 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
515 else {//Problem for 3D tetrahera
516 //cout<< "Fj.getGroupNames().size()= "<<Fj.getGroupNames().size()<<", Fj.x()= "<<Fj.x()<<", Fj.y()= "<<Fj.y()<<", Fj.z()= "<<Fj.z()<<endl;
518 cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
519 cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
520 cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
521 *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
522 _runLogFile->close();
523 throw CdmathException("Boundary condition not set");
526 // if Fj is inside the domain
527 } else if (Fj.getNumberOfCells()==2 ){
530 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
531 cout << " ; vecteur normal=(";
532 for(int p=0; p<_Ndim; p++)
533 cout << normale[p] << ",";
536 Cell2 = _mesh.getCell(idCells[1]);
539 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
541 inv_dxj = 1/Cell2.getMeasure();
543 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
545 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
546 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
547 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
548 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
552 *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
553 throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
557 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
558 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
559 VecAssemblyBegin(_b);
562 if(_onlyNeumannBC) //Check that the matrix is symmetric
564 PetscBool isSymetric;
565 MatIsSymmetric(_A,_precision,&isSymetric);
568 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
569 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
573 _diffusionMatrixSet=true;
579 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)
581 VecAssemblyBegin(_b);
584 for (int i=0; i<_Nmailles;i++)
585 VecSetValue(_b,i,_heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i),ADD_VALUES);
589 std::vector< int > nodesId;
590 for (int i=0; i<_Nmailles;i++)
593 nodesId=Ci.getNodesId();
594 for (int j=0; j<nodesId.size();j++)
595 if(!_mesh.isBorderNode(nodesId[j]))
597 double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
598 VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
605 if(_verbose or _system)
606 VecView(_b,PETSC_VIEWER_STDOUT_SELF);
612 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
616 //Only implicit scheme considered
617 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
618 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
620 #if PETSC_VERSION_GREATER_3_5
621 KSPSetOperators(_ksp, _A, _A);
623 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
627 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
628 KSPSolve(_ksp, _b, _Tk);
630 KSPConvergedReason reason;
631 KSPGetConvergedReason(_ksp,&reason);
632 KSPGetIterationNumber(_ksp, &_PetscIts);
634 KSPGetResidualNorm(_ksp,&residu);
635 if (reason!=2 and reason!=3)
637 cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
638 cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
639 cout<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
640 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
641 *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
642 *_runLogFile<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
643 _runLogFile->close();
648 if( _MaxIterLinearSolver < _PetscIts)
649 _MaxIterLinearSolver = _PetscIts;
650 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;
651 *_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;
652 VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
653 VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
658 VecAssemblyBegin(_Tk);
659 VecAssemblyEnd( _Tk);
660 VecAssemblyBegin(_deltaT);
661 VecAssemblyEnd( _deltaT);
664 for(int i=0; i<_Nmailles; i++)
666 VecGetValues(_deltaT, 1, &i, &dTi);
667 VecGetValues(_Tk, 1, &i, &Ti);
668 if(_erreur_rel < fabs(dTi/Ti))
669 _erreur_rel = fabs(dTi/Ti);
672 for(int i=0; i<_NunknownNodes; i++)
674 VecGetValues(_deltaT, 1, &i, &dTi);
675 VecGetValues(_Tk, 1, &i, &Ti);
676 if(_erreur_rel < fabs(dTi/Ti))
677 _erreur_rel = fabs(dTi/Ti);
680 converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
688 void StationaryDiffusionEquation::setMesh(const Mesh &M)
690 if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
692 cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
693 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
694 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
695 _runLogFile->close();
696 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
700 _Nmailles = _mesh.getNumberOfCells();
701 _Nnodes = _mesh.getNumberOfNodes();
703 cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
704 *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
706 // find maximum nb of neibourghs
709 _VV=Field ("Temperature", CELLS, _mesh, 1);
710 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
714 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
715 cout<<"1D Finite element method on segments"<<endl;
718 if( _mesh.isTriangular() )//Mesh dim=2
719 cout<<"2D Finite element method on triangles"<<endl;
720 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
721 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
724 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
725 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
726 _runLogFile->close();
727 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
732 if( _mesh.isTetrahedral() )//Mesh dim=3
733 cout<<"3D Finite element method on tetrahedra"<<endl;
734 else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
735 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
736 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
737 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
740 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;
741 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
742 _runLogFile->close();
743 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
747 _VV=Field ("Temperature", NODES, _mesh, 1);
749 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
750 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
751 _NboundaryNodes=_boundaryNodeIds.size();
757 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
759 //_maxPetscIts=maxIterationsPetsc;
760 // set linear solver algorithm
762 _ksptype = (char*)&KSPGMRES;
763 else if (kspType==CG)
764 _ksptype = (char*)&KSPCG;
765 else if (kspType==BCGS)
766 _ksptype = (char*)&KSPBCGS;
768 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
769 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
770 _runLogFile->close();
771 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
773 // set preconditioner
775 _pctype = (char*)&PCNONE;
776 else if (pcType ==LU)
777 _pctype = (char*)&PCLU;
778 else if (pcType == ILU)
779 _pctype = (char*)&PCILU;
780 else if (pcType ==CHOLESKY)
781 _pctype = (char*)&PCCHOLESKY;
782 else if (pcType == ICC)
783 _pctype = (char*)&PCICC;
785 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
786 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
787 _runLogFile->close();
788 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
792 bool StationaryDiffusionEquation::solveStationaryProblem()
794 if(!_initializedMemory)
796 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
797 _runLogFile->close();
798 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
800 bool stop=false; // Does the Problem want to stop (error) ?
801 bool converged=false; // has the newton scheme converged (end) ?
803 cout<< "!!! Running test case "<< _fileName << " using ";
804 *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
808 cout<< "Finite volumes method"<<endl<<endl;
809 *_runLogFile<< "Finite volumes method"<<endl<<endl;
813 cout<< "Finite elements method"<<endl<<endl;
814 *_runLogFile<< "Finite elements method"<< endl<<endl;
817 computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
819 cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
820 *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
821 _runLogFile->close();
822 throw CdmathException("Failed computing diffusion matrix");
824 computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
826 cout << "Error : failed computing right hand side, stopping calculation"<< endl;
827 *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
828 throw CdmathException("Failed computing right hand side");
830 stop = iterateNewtonStep(converged);
832 cout << "Error : failed solving linear system, stopping calculation"<< endl;
833 *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
834 _runLogFile->close();
835 throw CdmathException("Failed solving linear system");
838 _computationCompletedSuccessfully=true;
841 // Newton iteration loop for non linear problems
843 while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
845 computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
846 computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
847 stop = iterateNewtonStep(converged);
851 cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
852 *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
853 throw CdmathException("Failed solving a Newton iteration");
855 else if(_NEWTON_its==_maxNewtonIts){
856 cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
857 *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
858 throw CdmathException("No convergence of Newton scheme");
861 cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
862 *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
867 *_runLogFile<< "!!!!!! Computation successful !!!!!!"<< endl;
868 _runLogFile->close();
873 void StationaryDiffusionEquation::save(){
874 cout<< "Saving numerical results"<<endl<<endl;
875 *_runLogFile<< "Saving numerical results"<< endl<<endl;
877 string resultFile(_path+"/StationaryDiffusionEquation");//Results
880 resultFile+=_fileName;
882 // create mesh and component info
883 string suppress ="rm -rf "+resultFile+"_*";
884 system(suppress.c_str());//Nettoyage des précédents calculs identiques
886 if(_verbose or _system)
887 VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
891 for(int i=0; i<_Nmailles; i++)
893 VecGetValues(_Tk, 1, &i, &Ti);
899 for(int i=0; i<_NunknownNodes; i++)
901 VecGetValues(_Tk, 1, &i, &Ti);
902 globalIndex = globalNodeIndex(i, _dirichletNodeIds);
908 for(int i=0; i<_NdirichletNodes; i++)
910 Ni=_mesh.getNode(_dirichletNodeIds[i]);
911 nameOfGroup = Ni.getGroupName();
912 _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
916 _VV.setInfoOnComponent(0,"Temperature_(K)");
920 _VV.writeVTK(resultFile);
923 _VV.writeMED(resultFile);
926 _VV.writeCSV(resultFile);
931 StationaryDiffusionEquation::getOutputTemperatureField()
933 if(!_computationCompletedSuccessfully)
934 throw("Computation not performed yet or failed. No temperature field available");
939 void StationaryDiffusionEquation::terminate()
943 VecDestroy(&_deltaT);
948 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
950 _dirichletValuesSet=true;
951 _dirichletBoundaryValues=dirichletBoundaryValues;
955 StationaryDiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
957 _neumannValuesSet=true;
958 _neumannBoundaryValues=neumannBoundaryValues;
962 StationaryDiffusionEquation::getConditionNumber(bool isSingular, double tol) const
964 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
965 return A.getConditionNumber( isSingular, tol);
967 std::vector< double >
968 StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) const
970 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
972 if(_FECalculation)//We need to scale the FE matrix, otherwise the eigenvalues go to zero as the mesh is refined
974 Vector nodal_volumes(_NunknownNodes);
976 for(int i = 0; i< _Nmailles ; i++)//On parcourt les cellules du maillage
978 Cell Ci = _mesh.getCell(i);
979 for(int j = 0 ; j<_Ndim+1 ; j++)//On parcourt les noeuds de la cellule
981 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
983 j_int=unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
984 nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
988 for( j_int = 0; j_int< _NunknownNodes ; j_int++)
989 nodal_volumes[j_int]=1/nodal_volumes[j_int];
990 A.leftDiagonalScale(nodal_volumes);
993 return A.getEigenvalues( nev, which, tol);
995 std::vector< Vector >
996 StationaryDiffusionEquation::getEigenvectors(int nev, EPSWhich which, double tol) const
998 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
999 return A.getEigenvectors( nev, which, tol);
1002 StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, double tol) const
1004 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
1005 MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
1006 Field my_eigenfield;
1009 my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
1011 my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
1013 my_eigenfield.setFieldByDataArrayDouble(d);
1015 return my_eigenfield;