1 //============================================================================
2 // Name : Problème CoreFlows
5 // Copyright : CEA Saclay 2014
6 // Description : Classe générique des problèmes de Thermohydraulique
7 // Il s'agit d'une classe virtuelle pure (non instanciable qui regroupe
8 // les fonctionalités communes à tous les codes thermohydraulique.
9 // Le but est d'avoir une interface commune rendant
10 // plus lisibles et plus evolutifs les codes dévelopés avec CDMATH
11 //============================================================================
13 #include "ProblemCoreFlows.hxx"
14 #include "SparseMatrixPetsc.hxx"
21 ProblemCoreFlows::ProblemCoreFlows()
23 PetscBool petscInitialized;
24 PetscInitialized(&petscInitialized);
26 PetscInitialize(NULL,NULL,0,0);
34 _precision_Newton=_precision;
40 _fileName = "myCoreFlowsProblem";
42 _initialDataSet=false;
43 _initializedMemory=false;
44 _restartWithNewTimeScheme=false;
45 _restartWithNewFileName=false;
47 _wellBalancedCorrection=false;
50 _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
53 int _PetscIts=0;//the number of iterations of the linear solver
54 _ksptype = (char*)&KSPGMRES;
55 _pctype = (char*)&PCLU;
56 _heatPowerFieldSet=false;
57 _heatTransfertCoeff=0;
58 _rodTemperatureFieldSet=false;
62 _conditionNumber=false;
64 _runLogFile=new ofstream;
66 //extracting current directory
67 char result[ PATH_MAX ];
68 getcwd(result, PATH_MAX );
69 _path=string( result );
73 TimeScheme ProblemCoreFlows::getTimeScheme()
78 void ProblemCoreFlows::setTimeScheme(TimeScheme timeScheme)
80 if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
81 _restartWithNewTimeScheme=true;
82 _timeScheme = timeScheme;
85 bool ProblemCoreFlows::isStationary() const
90 double ProblemCoreFlows::presentTime() const
94 void ProblemCoreFlows::setTimeMax(double timeMax){
97 void ProblemCoreFlows::setPresentTime (double time)
101 void ProblemCoreFlows::setMaxNbOfTimeStep(int maxNbOfTimeStep){
102 _maxNbOfTimeStep = maxNbOfTimeStep;
104 void ProblemCoreFlows::setCFL(double cfl)
108 void ProblemCoreFlows::setPrecision(double precision)
110 _precision=precision;
112 void ProblemCoreFlows::setInitialField(const Field &VV)
115 if(_Ndim != VV.getSpaceDimension()){
116 *_runLogFile<<"ProblemCoreFlows::setInitialField: mesh has incorrect space dimension"<<endl;
117 _runLogFile->close();
118 throw CdmathException("ProblemCoreFlows::setInitialField: mesh has incorrect space dimension");
120 if(_nVar!=VV.getNumberOfComponents())
122 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect number of components"<<endl;
123 _runLogFile->close();
124 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect number of components");
130 _Nmailles = _mesh.getNumberOfCells();
131 _Nnodes = _mesh.getNumberOfNodes();
132 _Nfaces = _mesh.getNumberOfFaces();
133 _perimeters=Field("Perimeters", CELLS, _mesh,1);
135 // find _minl and maximum nb of neibourghs
137 int nbNeib,indexFace;
142 cout<<"Computing cell perimeters and mesh minimal diameter"<<endl;
144 if(VV.getTypeOfField()==NODES)
146 _minl = _mesh.getMaxNbNeighbours(NODES);
147 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
150 for (int i=0; i<_mesh.getNumberOfCells(); i++){
151 Ci = _mesh.getCell(i);
152 //Detect mesh with junction
154 for(int j=0; j<Ci.getNumberOfFaces(); j++){
155 Fk=_mesh.getFace(Ci.getFacesId()[j]);
156 nbNeib+=Fk.getNumberOfCells()-1;
158 if(nbNeib>_neibMaxNb)
163 for (int k=0 ; k<Ci.getNumberOfFaces() ; k++){
164 indexFace=Ci.getFacesId()[k];
165 Fk = _mesh.getFace(indexFace);
166 _minl = min(_minl,Ci.getMeasure()/Fk.getMeasure());
167 _perimeters(i)+=Fk.getMeasure();
170 _minl = min(_minl,Ci.getMeasure());
171 _perimeters(i)=Ci.getNumberOfFaces();
174 _initialDataSet=true;
177 cout<<_perimeters<<endl;
179 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber)
181 Field VV(fileName,CELLS,fieldName,timeStepNumber,0);
184 void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant)
187 Field VV("Primitive", CELLS, M, Vconstant.size());
189 for (int j = 0; j < M.getNumberOfCells(); j++) {
190 for (int i=0; i< VV.getNumberOfComponents(); i++)
191 VV(j,i) = Vconstant[i];
196 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant)
198 Field VV("Primitive", CELLS, M, Vconstant.getNumberOfRows());
200 for (int j = 0; j < M.getNumberOfCells(); j++) {
201 for (int i=0; i< VV.getNumberOfComponents(); i++)
202 VV(j,i) = Vconstant(i);
206 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant)
208 Field VV("Primitive", CELLS, M, Vconstant.size());
210 for (int j = 0; j < M.getNumberOfCells(); j++) {
211 for (int i=0; i< VV.getNumberOfComponents(); i++)
212 VV(j,i) = Vconstant[i];
216 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
217 double ymin, double ymax, int ny, string backSide, string frontSide,
218 double zmin, double zmax, int nz, string bottomSide, string topSide)
222 //cout<<"coucou1 xmin="<<xmin<<", xmax= "<<xmax<< ", nx= "<<nx<<endl;
223 M=Mesh(xmin,xmax,nx);
224 //cout<<"coucou2"<<endl;
227 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
229 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
231 cout<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
232 *_runLogFile<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
233 _runLogFile->close();
234 throw CdmathException("Space dimension nDim should be between 1 and 3");
237 M.setGroupAtPlan(xmax,0,_precision,rightSide);
238 M.setGroupAtPlan(xmin,0,_precision,leftSide);
240 M.setGroupAtPlan(ymax,1,_precision,frontSide);
241 M.setGroupAtPlan(ymin,1,_precision,backSide);
244 M.setGroupAtPlan(zmax,2,_precision,topSide);
245 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
248 setInitialFieldConstant(M, Vconstant);
250 void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction)
252 if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
254 *_runLogFile<<"ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes"<<endl;
255 _runLogFile->close();
256 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes");
258 Field VV("Primitive", CELLS, M, VV_Left.getNumberOfRows());
260 double component_value;
262 for (int j = 0; j < M.getNumberOfCells(); j++) {
264 component_value=M.getCell(j).x();
265 else if(direction==1)
266 component_value=M.getCell(j).y();
267 else if(direction==2)
268 component_value=M.getCell(j).z();
270 _runLogFile->close();
271 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: direction should be an integer between 0 and 2");
274 for (int i=0; i< VV.getNumberOfComponents(); i++)
275 if (component_value< disc_pos )
276 VV(j,i) = VV_Left[i];
278 VV(j,i) = VV_Right[i];
282 void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
283 double xmin, double xmax, int nx, string leftSide, string rightSide,
284 double ymin, double ymax, int ny, string backSide, string frontSide,
285 double zmin, double zmax, int nz, string bottomSide, string topSide)
289 M=Mesh(xmin,xmax,nx);
291 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
293 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
296 _runLogFile->close();
297 throw CdmathException("Space dimension nDim should be between 1 and 3");
300 M.setGroupAtPlan(xmax,0,_precision,rightSide);
301 M.setGroupAtPlan(xmin,0,_precision,leftSide);
303 M.setGroupAtPlan(ymax,1,_precision,frontSide);
304 M.setGroupAtPlan(ymin,1,_precision,backSide);
307 M.setGroupAtPlan(zmax,2,_precision,topSide);
308 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
310 Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
311 for(int i=0;i<VV_Left.size(); i++){
312 V_Left(i)=VV_Left[i];
313 V_Right(i)=VV_Right[i];
315 setInitialFieldStepFunction(M, V_Left, V_Right, xstep);
318 void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center)
320 if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
322 cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
323 throw CdmathException("ProblemCoreFlows::setInitialFieldSphericalStepFunction : Vector size error");
325 int nVar=Vout.size();
326 int spaceDim=M.getSpaceDimension();
327 Field VV("Primitive variables for spherical step function", CELLS, M, nVar);
329 Vector currentPoint(spaceDim);
330 for(int i=0;i<M.getNumberOfCells();i++)
332 currentPoint(0)=M.getCell(i).x();
335 currentPoint(1)=M.getCell(i).y();
337 currentPoint(2)=M.getCell(i).z();
339 if((currentPoint-Center).norm()<radius)
340 for(int j=0;j<nVar;j++)
343 for(int j=0;j<nVar;j++)
349 double ProblemCoreFlows::getTime()
353 unsigned ProblemCoreFlows::getNbTimeStep()
357 double ProblemCoreFlows::getCFL()
361 double ProblemCoreFlows::getPrecision()
365 Mesh ProblemCoreFlows::getMesh()
370 void ProblemCoreFlows::setLinearSolver(linearSolver kspType, preconditioner pcType)
372 //_maxPetscIts=maxIterationsPetsc;
373 // set linear solver algorithm
375 _ksptype = (char*)&KSPGMRES;
376 else if (kspType==CG)
377 _ksptype = (char*)&KSPCG;
378 else if (kspType==BCGS)
379 _ksptype = (char*)&KSPBCGS;
381 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
382 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
383 _runLogFile->close();
384 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
386 // set preconditioner
388 _pctype = (char*)&PCNONE;
389 else if (pcType ==LU)
390 _pctype = (char*)&PCLU;
391 else if (pcType == ILU)
392 _pctype = (char*)&PCILU;
393 else if (pcType ==CHOLESKY)
394 _pctype = (char*)&PCCHOLESKY;
395 else if (pcType == ICC)
396 _pctype = (char*)&PCICC;
398 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
399 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
400 _runLogFile->close();
401 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
406 // Cette methode lance une execution du ProblemCoreFlows
407 // Elle peut etre utilisee si le probleme n'est couple a aucun autre.
408 // (s'il n'a besoin d'aucun champ d'entree).
409 // Precondition: initialize
410 // Seule la methode terminate peut etre appelée apres
411 bool ProblemCoreFlows::run()
413 if(!_initializedMemory)
415 _runLogFile->close();
416 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
418 bool stop=false; // Does the Problem want to stop (error) ?
419 bool ok; // Is the time interval successfully solved ?
420 _isStationary=false;//in case of a second run with a different physics or cfl
422 cout<< "Running test case "<< _fileName<<endl;
424 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
425 *_runLogFile<< "Running test case "<< _fileName<<endl;
428 while(!stop && !_isStationary &&_time<_timeMax && _nbTimeStep<_maxNbOfTimeStep)
430 ok=false; // Is the time interval successfully solved ?
432 // Guess the next time step length
433 _dt=computeTimeStep(stop);
435 cout << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
436 *_runLogFile << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
439 // Loop on the time interval tries
440 while (!ok && !stop )
442 stop=!initTimeStep(_dt);
443 // Prepare the next time step
445 cout << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
446 *_runLogFile << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
449 // Solve the next time step
452 if (!ok) // The resolution failed, try with a new time interval.
455 cout<<"ComputeTimeStep returned _dt="<<_dt<<endl;
456 cout << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
457 *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
458 double dt=_dt/2;//We chose to divide the time step by 2
459 abortTimeStep();//Cancel the initTimeStep
460 _dt=dt;//new value of time step is previous time step divided by 2 (we do not call computeTimeStep
461 //_cfl*=0.5;//If we change the cfl, we must compute the new time step with computeTimeStep
462 //_dt=computeTimeStep(stop);
465 cout << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
466 *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
467 stop=true; // Impossible to solve the next time step, the Problem has given up
471 else // The resolution was successful, validate and go to the next time step.
474 if (_nbTimeStep%_freqSave ==0){
475 cout << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl;
476 *_runLogFile << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl;
482 cout << "Stationary state reached" <<endl;
483 *_runLogFile << "Stationary state reached" <<endl;
485 else if(_time>=_timeMax){
486 cout<<"Maximum time "<<_timeMax<<" reached"<<endl;
487 *_runLogFile<<"Maximum time "<<_timeMax<<" reached"<<endl;
489 else if(_nbTimeStep>=_maxNbOfTimeStep){
490 cout<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
491 *_runLogFile<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
494 cout<<"Error problem wants to stop!"<<endl;
495 *_runLogFile<<"Error problem wants to stop!"<<endl;
497 cout << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
498 *_runLogFile << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
500 _runLogFile->close();
504 void ProblemCoreFlows::displayMatrix(double *matrix, int size, string name)
507 for(int p=0; p<size; p++)
509 for(int q=0; q<size; q++)
511 cout << matrix[p*size+q] << "\t";
518 void ProblemCoreFlows::displayVector(double *vector, int size, string name)
521 for(int q=0; q<size; q++)
523 cout << vector[q] << "\t";
527 void ProblemCoreFlows::setFileName(string fileName){
528 if( _nbTimeStep>0 && fileName!=_fileName)//This is a change of file name during a simulation
529 _restartWithNewFileName=true;
530 _fileName = fileName;
533 void ProblemCoreFlows::setFreqSave(int freqSave){
534 _freqSave = freqSave;
536 bool ProblemCoreFlows::solveTimeStep(){
538 bool converged=false, ok=true;
539 while(!converged && ok && _NEWTON_its < _maxNewtonIts){
540 ok=iterateTimeStep(converged);//resolution du systeme lineaire si schema implicite
542 if(_timeScheme == Implicit && _nbTimeStep%_freqSave ==0)//To monitor the convergence of the newton scheme
544 cout << "\n Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << " : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
545 *_runLogFile<< "\n Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << " : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
549 PetscReal sv_max, sv_min;
550 KSPComputeExtremeSingularValues(_ksp, &sv_max, &sv_min);
551 cout<<" singular value max = " << sv_max <<" singular value min = " << sv_min <<" condition number = " << sv_max/sv_min <<endl;
552 *_runLogFile<<" singular value max = " << sv_max <<" singular value min = " << sv_min <<" condition number = " << sv_max/sv_min <<endl;
558 if(_NEWTON_its >= _maxNewtonIts){
559 cout << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
560 *_runLogFile << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
563 cout<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
564 *_runLogFile<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
567 else if(_timeScheme == Implicit && _nbTimeStep%_freqSave ==0)
570 cout << "Nombre d'iterations de Newton "<< _NEWTON_its << ", Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl;
572 *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its << " Variation ||Un+1-Un||= "<< _erreur_rel<<endl;
573 *_runLogFile << "Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl;
574 _MaxIterLinearSolver = 0;
579 ProblemCoreFlows::~ProblemCoreFlows()
582 PetscBool petscInitialized;
583 PetscInitialized(&petscInitialized);
591 ProblemCoreFlows::getConditionNumber(bool isSingular, double tol) const
593 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
594 return A.getConditionNumber( isSingular, tol);
596 std::vector< double >
597 ProblemCoreFlows::getEigenvalues(int nev, EPSWhich which, double tol) const
599 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
600 return A.getEigenvalues( nev, which, tol);
602 std::vector< Vector >
603 ProblemCoreFlows::getEigenvectors(int nev, EPSWhich which, double tol) const
605 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
606 return A.getEigenvectors( nev, which, tol);
609 ProblemCoreFlows::getEigenvectorsField(int nev, EPSWhich which, double tol) const
611 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
612 MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
616 my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
618 my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
620 my_eigenfield.setFieldByDataArrayDouble(d);
622 return my_eigenfield;