1 #include "StationaryDiffusionEquation.hxx"
2 #include "DiffusionEquation.hxx"
4 #include "SparseMatrixPetsc.hxx"
12 StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda,MPI_Comm comm){
13 /* Initialisation of PETSC */
14 //check if PETSC is already initialised
15 PetscBool petscInitialized;
16 PetscInitialized(&petscInitialized);
18 {//check if MPI is already initialised
20 MPI_Initialized(&mpiInitialized);
22 PETSC_COMM_WORLD = comm;
23 PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
25 MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
26 MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
27 PetscPrintf(PETSC_COMM_WORLD,"\n Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
28 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
29 PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
33 std::cout<<"Conductivity="<<lambda<<endl;
34 throw CdmathException("!!!!!!!!Error : conductivity parameter lambda cannot be negative");
38 std::cout<<"Space dimension="<<dim<<endl;
39 throw CdmathException("!!!!!!!!Error : parameter dim cannot be negative");
42 _FECalculation=FECalculation;
46 _nVar=1;//scalar prolem
48 _diffusionMatrixSet=false;
49 _initializedMemory=false;
57 _boundaryNodeIds=std::vector< int >(0);
58 _dirichletNodeIds=std::vector< int >(0);
62 _dirichletValuesSet=false;
63 _neumannValuesSet=false;
67 _precision_Newton=_precision;
68 _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
72 int _PetscIts=0;//the number of iterations of the linear solver
73 _ksptype = (char*)&KSPGMRES;
75 _pctype = (char*)&PCNONE;
77 _pctype = (char*)&PCLU;
79 _conditionNumber=false;
82 //parameters for monitoring simulation
85 _runLogFile=new ofstream;
87 //result save parameters
88 _fileName = "StationaryDiffusionProblem";
89 char result[ PATH_MAX ];//extracting current directory
90 getcwd(result, PATH_MAX );
91 _path=string( result );
93 _computationCompletedSuccessfully=false;
95 //heat transfer parameters
97 _fluidTemperatureFieldSet=false;
99 _heatPowerFieldSet=false;
100 _heatTransfertCoeff=0;
103 /* Default diffusion tensor is diagonal */
104 _DiffusionTensor=Matrix(_Ndim);
105 for(int idim=0;idim<_Ndim;idim++)
106 _DiffusionTensor(idim,idim)=_conductivity;
108 PetscPrintf(PETSC_COMM_WORLD,"\n Stationary diffusion problem with conductivity %.2f", lambda);
110 PetscPrintf(PETSC_COMM_WORLD," and finite elements method\n\n");
112 PetscPrintf(PETSC_COMM_WORLD," and finite volumes method\n\n");
115 void StationaryDiffusionEquation::initialize()
119 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
122 throw CdmathException("!!!!!!!!StationaryDiffusionEquation::initialize() set mesh first");
125 cout<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
126 *_runLogFile<<"\n Initialisation of the computation of the temperature diffusion in a solid using ";
129 cout<< "Finite volumes method"<<endl<<endl;
130 *_runLogFile<< "Finite volumes method"<<endl<<endl;
134 cout<< "Finite elements method"<<endl<<endl;
135 *_runLogFile<< "Finite elements method"<<endl<<endl;
139 /**************** Field creation *********************/
141 if(!_heatPowerFieldSet){
142 _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
143 for(int i =0; i<_VV.getNumberOfElements(); i++)
144 _heatPowerField(i) = _heatSource;
145 _heatPowerFieldSet=true;
147 if(!_fluidTemperatureFieldSet){
148 _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
149 for(int i =0; i<_VV.getNumberOfElements(); i++)
150 _fluidTemperatureField(i) = _fluidTemperature;
151 _fluidTemperatureFieldSet=true;
154 /* Détection des noeuds frontière avec une condition limite de Dirichlet */
157 if(_NboundaryNodes==_Nnodes)
158 cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
160 for(int i=0; i<_NboundaryNodes; i++)
162 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
163 if( it != _dirichletBoundaryValues.end() )
164 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
165 else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
167 cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
168 *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
169 _runLogFile->close();
170 throw CdmathException("Missing boundary group");
172 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
174 cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
175 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
176 *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
177 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
178 _runLogFile->close();
179 throw CdmathException("Missing boundary condition");
181 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
182 _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
183 else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
185 cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
186 cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
187 *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
188 *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
189 _runLogFile->close();
190 throw CdmathException("Wrong boundary condition");
193 _NdirichletNodes=_dirichletNodeIds.size();
194 _NunknownNodes=_Nnodes - _NdirichletNodes;
195 cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
196 *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
201 _globalNbUnknowns = _Nmailles*_nVar;
203 MPI_Bcast(&_NunknownNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
204 _globalNbUnknowns = _NunknownNodes*_nVar;
207 /* Vectors creations */
208 VecCreate(PETSC_COMM_WORLD, &_Tk);//main unknown
209 VecSetSizes(_Tk,PETSC_DECIDE,_globalNbUnknowns);
210 VecSetFromOptions(_Tk);
211 VecGetLocalSize(_Tk, &_localNbUnknowns);
213 VecDuplicate(_Tk, &_Tkm1);
214 VecDuplicate(_Tk, &_deltaT);
215 VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
217 /* Matrix creation */
218 MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
220 /* Local sequential vector creation */
221 if(_mpi_size>1 && _mpi_rank == 0)
222 VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Tk_seq);//For saving results on proc 0
225 VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq);
228 KSPCreate(PETSC_COMM_WORLD, &_ksp);
229 KSPSetType(_ksp, _ksptype);
230 // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
231 KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
232 KSPGetPC(_ksp, &_pc);
233 PCSetType(_pc, _pctype);
235 //Checking whether all boundary conditions are Neumann boundary condition
236 //if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
237 if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
239 map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
240 while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
242 _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
246 _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
248 _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
250 //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
253 std::cout<<"### Warning : all boundary conditions are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
254 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;
255 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;
256 *_runLogFile<<"### Warning : all boundary condition are Neumann. System matrix is not invertible since constant vectors are in the kernel."<<std::endl;
257 *_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;
258 *_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;
261 MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
262 MatSetNullSpace(_A, nullsp);
263 MatSetTransposeNullSpace(_A, nullsp);
264 MatNullSpaceDestroy(&nullsp);
265 //PCFactorSetShiftType(_pc,MAT_SHIFT_NONZERO);
266 //PCFactorSetShiftAmount(_pc,1e-10);
268 _initializedMemory=true;
271 double StationaryDiffusionEquation::computeTimeStep(bool & stop){
272 if(!_diffusionMatrixSet)//The diffusion matrix is computed once and for all time steps
273 computeDiffusionMatrix(stop);
275 _dt_src=computeRHS(stop);
280 double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
288 result=computeDiffusionMatrixFE(stop);
290 result=computeDiffusionMatrixFV(stop);
292 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
293 MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
295 //Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
296 //update value here if variable heat transfer coefficient
297 if(_heatTransfertCoeff>_precision)
298 MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
300 if(_verbose or _system)
301 MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
306 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
313 Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
314 std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
315 std::vector< int > nodeIds(_Ndim+1);//cell node Ids
316 std::vector< Node > nodes(_Ndim+1);//cell nodes
317 int i_int, j_int; //index of nodes j and k considered as unknown nodes
318 bool dirichletCell_treated;
320 std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
321 for (int idim=0; idim<_Ndim+1;idim++)
322 values[idim][idim]=1;
324 /* parameters for boundary treatment */
325 vector< double > valuesBorder(_Ndim+1);
326 Vector GradShapeFuncBorder(_Ndim+1);
328 for (int j=0; j<_Nmailles;j++)
330 Cj = _mesh.getCell(j);
332 for (int idim=0; idim<_Ndim+1;idim++){
333 nodeIds[idim]=Cj.getNodeId(idim);
334 nodes[idim]=_mesh.getNode(nodeIds[idim]);
335 for (int jdim=0; jdim<_Ndim;jdim++)
336 M(idim,jdim)=nodes[idim].getPoint()[jdim];
339 for (int idim=0; idim<_Ndim+1;idim++)
340 GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim);
342 /* Loop on the edges of the cell */
343 for (int idim=0; idim<_Ndim+1;idim++)
345 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
346 {//First node of the edge is not Dirichlet node
347 i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
348 dirichletCell_treated=false;
349 for (int jdim=0; jdim<_Ndim+1;jdim++)
351 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
352 {//Second node of the edge is not Dirichlet node
353 j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
354 MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
356 else if (!dirichletCell_treated)
357 {//Second node of the edge is a Dirichlet node
358 dirichletCell_treated=true;
359 for (int kdim=0; kdim<_Ndim+1;kdim++)
361 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
362 if( it != _dirichletBoundaryValues.end() )
364 if( _dirichletValuesSet )
365 valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
367 valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
370 valuesBorder[kdim]=0;
372 GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
373 coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
374 VecSetValue(_b,i_int,coeff, ADD_VALUES);
381 //Calcul de la contribution de la condition limite de Neumann au second membre
382 if( _NdirichletNodes !=_NboundaryNodes)
384 vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
385 int NboundaryFaces=boundaryFaces.size();
386 for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
388 Face Fi = _mesh.getFace(boundaryFaces[i]);
389 for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
391 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
393 j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
394 if( _neumannValuesSet )
395 coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
397 coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
398 VecSetValue(_b, j_int, coeff, ADD_VALUES);
405 if(_onlyNeumannBC) //Check that the matrix is symmetric
407 PetscBool isSymetric;
408 MatIsSymmetric(_A,_precision,&isSymetric);
411 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
412 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
416 _diffusionMatrixSet=true;
422 double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
425 long nbFaces = _mesh.getNumberOfFaces();
429 double inv_dxi, inv_dxj;
430 double barycenterDistance;
431 Vector normale(_Ndim);
434 std::vector< int > idCells;
436 for (int j=0; j<nbFaces;j++){
437 Fj = _mesh.getFace(j);
439 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
440 idCells = Fj.getCellsId();
441 Cell1 = _mesh.getCell(idCells[0]);
443 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
444 if (j == Cell1.getFacesId()[l]){
445 for (int idim = 0; idim < _Ndim; ++idim)
446 normale[idim] = Cell1.getNormalVector(l,idim);
451 //Compute velocity at the face Fj
452 dn=(_DiffusionTensor*normale)*normale;
454 // compute 1/dxi = volume of Ci/area of Fj
455 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
457 // If Fj is on the boundary
458 if (Fj.getNumberOfCells()==1) {
461 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
462 for(int p=0; p<_Ndim; p++)
463 cout << normale[p] << ",";
467 std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
468 if( it != _dirichletBoundaryValues.end() )
470 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
471 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
472 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
476 nameOfGroup = Fj.getGroupName();
478 if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
479 VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
481 else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
482 barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
483 MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
484 VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
488 cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
489 cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
490 cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
491 *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
492 _runLogFile->close();
493 throw CdmathException("Boundary condition not set");
496 // if Fj is inside the domain
497 } else if (Fj.getNumberOfCells()==2 ){
500 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
501 cout << " ; vecteur normal=(";
502 for(int p=0; p<_Ndim; p++)
503 cout << normale[p] << ",";
506 Cell2 = _mesh.getCell(idCells[1]);
509 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
511 inv_dxj = 1/Cell2.getMeasure();
513 barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
515 MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
516 MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
517 MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
518 MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
522 *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
523 throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
528 if(_onlyNeumannBC) //Check that the matrix is symmetric
530 PetscBool isSymetric;
531 MatIsSymmetric(_A,_precision,&isSymetric);
534 cout<<"Singular matrix is not symmetric, tolerance= "<< _precision<<endl;
535 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
539 _diffusionMatrixSet=true;
545 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
551 for (int i=0; i<_Nmailles;i++)
552 VecSetValue(_b,i, _heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i) ,ADD_VALUES);
556 std::vector< int > nodesId;
557 for (int i=0; i<_Nmailles;i++)
560 nodesId=Ci.getNodesId();
561 for (int j=0; j<nodesId.size();j++)
562 if(!_mesh.isBorderNode(nodesId[j]))
564 double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
565 VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
570 VecAssemblyBegin(_b);
573 if(_verbose or _system)
574 VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
580 bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
584 //Only implicit scheme considered
586 #if PETSC_VERSION_GREATER_3_5
587 KSPSetOperators(_ksp, _A, _A);
589 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
593 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
594 KSPSolve(_ksp, _b, _Tk);
596 KSPConvergedReason reason;
597 KSPGetConvergedReason(_ksp,&reason);
598 KSPGetIterationNumber(_ksp, &_PetscIts);
600 KSPGetResidualNorm(_ksp,&residu);
601 if (reason!=2 and reason!=3)
603 PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc.");
604 PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Itérations maximales %d atteintes, résidu = %1.2e, précision demandée= %1.2e",_maxPetscIts,residu,_precision);
605 PetscPrintf(PETSC_COMM_WORLD,"Solver used %s, preconditioner %s, Final number of iteration = %d",_ksptype,_pctype,_PetscIts);
606 *_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
607 *_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
608 *_runLogFile<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
609 _runLogFile->close();
614 if( _MaxIterLinearSolver < _PetscIts)
615 _MaxIterLinearSolver = _PetscIts;
616 PetscPrintf(PETSC_COMM_WORLD,"## Système linéaire résolu en %d itérations par le solveur %s et le preconditioneur %s, précision demandée = %1.2e",_PetscIts,_ksptype,_pctype,_precision);
617 *_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;
618 VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
619 VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
621 VecAssemblyBegin(_Tk);
622 VecAssemblyEnd( _Tk);
623 VecAssemblyBegin(_deltaT);
624 VecAssemblyEnd( _deltaT);
627 PetscPrintf(PETSC_COMM_WORLD,"Début calcul de l'erreur maximale");
629 VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
632 PetscPrintf(PETSC_COMM_WORLD,"Fin calcul de la variation relative, erreur maximale : %1.2e", _erreur_rel );
634 converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
642 void StationaryDiffusionEquation::setMesh(const Mesh &M)
646 if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
648 cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
649 *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
650 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
651 _runLogFile->close();
652 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
656 _Nmailles = _mesh.getNumberOfCells();
657 _Nnodes = _mesh.getNumberOfNodes();
659 cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
660 *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
662 // find maximum nb of neibourghs
665 _VV=Field ("Temperature", CELLS, _mesh, 1);
666 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
670 if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
671 cout<<"1D Finite element method on segments"<<endl;
674 if( _mesh.isTriangular() )//Mesh dim=2
675 cout<<"2D Finite element method on triangles"<<endl;
676 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
677 cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
680 cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
681 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
682 _runLogFile->close();
683 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
688 if( _mesh.isTetrahedral() )//Mesh dim=3
689 cout<<"3D Finite element method on tetrahedra"<<endl;
690 else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
691 cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
692 else if (_mesh.getMeshDimension()==1)//Mesh dim=1
693 cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
696 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;
697 *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
698 _runLogFile->close();
699 throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
703 _VV=Field ("Temperature", NODES, _mesh, 1);
705 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
706 _boundaryNodeIds = _mesh.getBoundaryNodeIds();
707 _NboundaryNodes=_boundaryNodeIds.size();
711 /* MPI distribution parameters */
712 int nbVoisinsMax;//Mettre en attribut ?
714 MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD);
715 MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
716 nbVoisinsMax = _neibMaxNbCells;
719 MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD);
720 MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
721 nbVoisinsMax = _neibMaxNbNodes;
723 _d_nnz = (nbVoisinsMax+1)*_nVar;
724 _o_nnz = nbVoisinsMax *_nVar;
729 void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, preconditioner pcType)
731 //_maxPetscIts=maxIterationsPetsc;
732 // set linear solver algorithm
734 _ksptype = (char*)&KSPGMRES;
735 else if (kspType==CG)
736 _ksptype = (char*)&KSPCG;
737 else if (kspType==BCGS)
738 _ksptype = (char*)&KSPBCGS;
740 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
741 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
742 _runLogFile->close();
743 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
745 // set preconditioner
747 _pctype = (char*)&PCNONE;
748 else if (pcType ==LU)
749 _pctype = (char*)&PCLU;
750 else if (pcType == ILU)
751 _pctype = (char*)&PCILU;
752 else if (pcType ==CHOLESKY)
753 _pctype = (char*)&PCCHOLESKY;
754 else if (pcType == ICC)
755 _pctype = (char*)&PCICC;
757 cout << "!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
758 *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
759 _runLogFile->close();
760 throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
764 bool StationaryDiffusionEquation::solveStationaryProblem()
766 if(!_initializedMemory)
768 *_runLogFile<< "ProblemCoreFlows::run() call initialize() first"<< _fileName<<endl;
769 _runLogFile->close();
770 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
772 bool stop=false; // Does the Problem want to stop (error) ?
773 bool converged=false; // has the newton scheme converged (end) ?
775 PetscPrintf(PETSC_COMM_WORLD,"!!! Running test case %s using ",_fileName);
776 *_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
780 PetscPrintf(PETSC_COMM_WORLD,"Finite volumes method\n\n");
781 *_runLogFile<< "Finite volumes method"<<endl<<endl;
785 PetscPrintf(PETSC_COMM_WORLD,"Finite elements method\n\n");
786 *_runLogFile<< "Finite elements method"<< endl<<endl;
789 computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
791 PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing diffusion matrix, stopping calculation\n");
792 *_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
793 _runLogFile->close();
794 throw CdmathException("Failed computing diffusion matrix");
796 computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
798 PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing right hand side, stopping calculation\n");
799 *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
800 throw CdmathException("Failed computing right hand side");
802 stop = iterateNewtonStep(converged);
804 PetscPrintf(PETSC_COMM_WORLD,"Error : failed solving linear system, stopping calculation\n");
805 *_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
806 _runLogFile->close();
807 throw CdmathException("Failed solving linear system");
810 _computationCompletedSuccessfully=true;
813 // Newton iteration loop for non linear problems
815 while(!stop and !converged and _NEWTON_its<_maxNewtonIts)
817 computeDiffusionMatrix( stop);//case when the conductivity depends on the temperature (nonlinear LHS)
818 computeRHS(stop);//case the heat power depends on the unknown temperature (nonlinear RHS)
819 stop = iterateNewtonStep(converged);
823 cout << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
824 *_runLogFile << "Error : failed solving Newton iteration "<<_NEWTON_its<<", stopping calculation"<< endl;
825 throw CdmathException("Failed solving a Newton iteration");
827 else if(_NEWTON_its==_maxNewtonIts){
828 cout << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
829 *_runLogFile << "Error : no convergence of Newton scheme. Maximum Newton iterations "<<_maxNewtonIts<<" reached, stopping calculation"<< endl;
830 throw CdmathException("No convergence of Newton scheme");
833 cout << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
834 *_runLogFile << "Convergence of Newton scheme at iteration "<<_NEWTON_its<<", end of calculation"<< endl;
839 *_runLogFile<< "!!!!!! Computation successful !!!!!!"<< endl;
840 _runLogFile->close();
845 void StationaryDiffusionEquation::save(){
846 PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n");
847 *_runLogFile<< "Saving numerical results"<< endl<<endl;
849 string resultFile(_path+"/StationaryDiffusionEquation");//Results
852 resultFile+=_fileName;
854 // create mesh and component info
855 string suppress ="rm -rf "+resultFile+"_*";
856 system(suppress.c_str());//Nettoyage des précédents calculs identiques
859 VecScatterBegin(_scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
860 VecScatterEnd( _scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
863 if(_verbose or _system)
864 VecView(_Tk,PETSC_VIEWER_STDOUT_WORLD);
869 for(int i=0; i<_Nmailles; i++)
872 VecGetValues(_Tk_seq, 1, &i, &Ti);
874 VecGetValues(_Tk , 1, &i, &Ti);
880 for(int i=0; i<_NunknownNodes; i++)
883 VecGetValues(_Tk_seq, 1, &i, &Ti);
885 VecGetValues(_Tk , 1, &i, &Ti);
886 globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
892 for(int i=0; i<_NdirichletNodes; i++)
894 Ni=_mesh.getNode(_dirichletNodeIds[i]);
895 nameOfGroup = Ni.getGroupName();
896 _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
900 _VV.setInfoOnComponent(0,"Temperature_(K)");
904 _VV.writeVTK(resultFile);
907 _VV.writeMED(resultFile);
910 _VV.writeCSV(resultFile);
916 void StationaryDiffusionEquation::terminate()
920 VecDestroy(&_deltaT);
923 if(_mpi_size>1 && _mpi_rank == 0)
924 VecDestroy(&_Tk_seq);
926 PetscBool petscInitialized;
927 PetscInitialized(&petscInitialized);
934 StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
936 _dirichletValuesSet=true;
937 _dirichletBoundaryValues=dirichletBoundaryValues;
941 StationaryDiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues)
943 _neumannValuesSet=true;
944 _neumannBoundaryValues=neumannBoundaryValues;
948 StationaryDiffusionEquation::getConditionNumber(bool isSingular, double tol) const
950 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
951 return A.getConditionNumber( isSingular, tol);
953 std::vector< double >
954 StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) const
956 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
958 if(_FECalculation)//We need to scale the FE matrix, otherwise the eigenvalues go to zero as the mesh is refined
960 Vector nodal_volumes(_NunknownNodes);
962 for(int i = 0; i< _Nmailles ; i++)//On parcourt les cellules du maillage
964 Cell Ci = _mesh.getCell(i);
965 for(int j = 0 ; j<_Ndim+1 ; j++)//On parcourt les noeuds de la cellule
967 if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
969 j_int=DiffusionEquation::unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
970 nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
974 for( j_int = 0; j_int< _NunknownNodes ; j_int++)
975 nodal_volumes[j_int]=1/nodal_volumes[j_int];
976 A.leftDiagonalScale(nodal_volumes);
979 return A.getEigenvalues( nev, which, tol);
981 std::vector< Vector >
982 StationaryDiffusionEquation::getEigenvectors(int nev, EPSWhich which, double tol) const
984 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
985 return A.getEigenvectors( nev, which, tol);
988 StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, double tol) const
990 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
991 MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
994 my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
996 my_eigenfield.setFieldByDataArrayDouble(d);
998 return my_eigenfield;
1002 StationaryDiffusionEquation::getOutputTemperatureField()
1004 if(!_computationCompletedSuccessfully)
1005 throw("Computation not performed yet or failed. No temperature field available");
1011 StationaryDiffusionEquation::getRodTemperatureField()
1013 return getOutputTemperatureField();
1017 StationaryDiffusionEquation::getInputFieldsNames()
1019 vector<string> result(2);
1021 result[0]="FluidTemperature";
1022 result[1]="HeatPower";
1027 StationaryDiffusionEquation::getOutputFieldsNames()
1029 vector<string> result(1);
1031 result[0]="RodTemperature";
1037 StationaryDiffusionEquation::getOutputField(const string& nameField )
1039 if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
1040 return getRodTemperatureField();
1043 cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
1044 throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
1049 StationaryDiffusionEquation::setInputField(const string& nameField, Field& inputField )
1052 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setInputField set the mesh first");
1054 if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
1055 return setFluidTemperatureField( inputField) ;
1056 else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
1057 return setHeatPowerField( inputField );
1060 cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
1061 throw CdmathException("StationaryDiffusionEquation::setInputField error : Unknown Field name");
1066 StationaryDiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){
1068 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setFluidTemperatureField set initial field first");
1070 coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh);
1071 _fluidTemperatureField=coupledTemperatureField;
1072 _fluidTemperatureFieldSet=true;
1076 StationaryDiffusionEquation::setHeatPowerField(Field heatPower){
1078 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
1080 heatPower.getMesh().checkFastEquivalWith(_mesh);
1081 _heatPowerField=heatPower;
1082 _heatPowerFieldSet=true;
1086 StationaryDiffusionEquation::setHeatPowerField(string fileName, string fieldName, int iteration, int order, int meshLevel){
1088 throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first");
1090 _heatPowerField=Field(fileName, CELLS,fieldName, iteration, order, meshLevel);
1091 _heatPowerField.getMesh().checkFastEquivalWith(_mesh);
1092 _heatPowerFieldSet=true;