CoreFlows/Models/src/ProblemCoreFlows.cxx

   1 //============================================================================
   2 // Name        : Problème CoreFlows
   3 // Author      : M. Ndjinga
   4 // Version     :
   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 //============================================================================
  12
  13 #include "ProblemCoreFlows.hxx"
  14 #include <limits.h>
  15 #include <unistd.h>
  16
  17 using namespace std;
  18
  19 ProblemCoreFlows::ProblemCoreFlows()
  20 {
  21         PetscBool petscInitialized;
  22         PetscInitialized(&petscInitialized);
  23         if(!petscInitialized)
  24                 PetscInitialize(NULL,NULL,0,0);
  25         _dt = 0;
  26         _time = 0;
  27         _nbTimeStep=0;
  28         _cfl=0;
  29         _timeMax =0.;
  30         _maxNbOfTimeStep = 0;
  31         _precision=1.e-6;
  32         _precision_Newton=_precision;
  33         _erreur_rel= 0;
  34         _isStationary=false;
  35         _Ndim=0;
  36         _minl=0;
  37         _neibMaxNb=0;
  38         _fileName = "myCoreFlowsProblem";
  39         _freqSave = 1;
  40         _initialDataSet=false;
  41         _initializedMemory=false;
  42         _spaceScheme=upwind;
  43         _timeScheme=Explicit;
  44         _wellBalancedCorrection=false;
  45         _maxPetscIts=50;
  46         _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
  47         _maxNewtonIts=50;
  48         _NEWTON_its=0;
  49         int _PetscIts=0;//the number of iterations of the linear solver
  50         _ksptype = (char*)&KSPGMRES;
  51         _pctype = (char*)&PCLU;
  52         _heatPowerFieldSet=false;
  53         _heatTransfertCoeff=0;
  54         _rodTemperatureFieldSet=false;
  55         _heatSource=0;
  56         _verbose = false;
  57         _system = false;
  58         _conditionNumber=false;
  59         _maxvp=0;
  60         _runLogFile=new ofstream;
  61
  62         //extracting current directory
  63         char result[ PATH_MAX ];
  64         getcwd(result, PATH_MAX );
  65         _path=string( result );
  66         _saveFormat=VTK;
  67 }
  68
  69 bool ProblemCoreFlows::isStationary() const
  70 {
  71         return _isStationary;
  72 }
  73
  74 double ProblemCoreFlows::presentTime() const
  75 {
  76         return _time;
  77 }
  78 void ProblemCoreFlows::setTimeMax(double timeMax){
  79         _timeMax = timeMax;
  80 }
  81 void ProblemCoreFlows::setPresentTime (double time)
  82 {
  83         _time=time;
  84 }
  85 void ProblemCoreFlows::setMaxNbOfTimeStep(int maxNbOfTimeStep){
  86         _maxNbOfTimeStep = maxNbOfTimeStep;
  87 }
  88 void ProblemCoreFlows::setCFL(double cfl)
  89 {
  90         _cfl=cfl;
  91 }
  92 void ProblemCoreFlows::setPrecision(double precision)
  93 {
  94         _precision=precision;
  95 }
  96 void ProblemCoreFlows::setNumericalScheme(SpaceScheme spaceScheme, TimeScheme timeScheme){
  97         _timeScheme = timeScheme;
  98         _spaceScheme = spaceScheme;
  99 }
 100
 101 void ProblemCoreFlows::setInitialField(const Field &VV)
 102 {
 103
 104         if(_Ndim != VV.getSpaceDimension()){
 105                 *_runLogFile<<"ProblemCoreFlows::setInitialField: mesh has incorrect space dimension"<<endl;
 106                 _runLogFile->close();
 107                 throw CdmathException("ProblemCoreFlows::setInitialField: mesh has incorrect space dimension");
 108         }
 109         if(_nVar!=VV.getNumberOfComponents())
 110         {
 111                 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect number of components"<<endl;
 112                 _runLogFile->close();
 113                 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect number of components");
 114         }
 115
 116         _VV=VV;
 117         _time=_VV.getTime();
 118         _mesh=_VV.getMesh();
 119         _Nmailles = _mesh.getNumberOfCells();
 120         _Nnodes =   _mesh.getNumberOfNodes();
 121         _Nfaces =   _mesh.getNumberOfFaces();
 122         _perimeters=Field("Perimeters", CELLS, _mesh,1);
 123
 124         // find _minl and maximum nb of neibourghs
 125         _minl  = INFINITY;
 126         int nbNeib,indexFace;
 127         Cell Ci;
 128         Face Fk;
 129
 130         if(_verbose)
 131                 cout<<"Computing cell perimeters and mesh minimal diameter"<<endl;
 132
 133     if(VV.getTypeOfField()==NODES)
 134     {
 135         _minl = _mesh.getMaxNbNeighbours(NODES);
 136         _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 137     }
 138     else
 139         for (int i=0; i<_mesh.getNumberOfCells(); i++){
 140             Ci = _mesh.getCell(i);
 141             //Detect mesh with junction
 142             nbNeib=0;
 143             for(int j=0; j<Ci.getNumberOfFaces(); j++){
 144                 Fk=_mesh.getFace(Ci.getFacesId()[j]);
 145                 nbNeib+=Fk.getNumberOfCells()-1;
 146             }
 147             if(nbNeib>_neibMaxNb)
 148                 _neibMaxNb=nbNeib;
 149             //Compute mesh data
 150             if (_Ndim > 1){
 151                 _perimeters(i)=0;
 152                 for (int k=0 ; k<Ci.getNumberOfFaces() ; k++){
 153                     indexFace=Ci.getFacesId()[k];
 154                     Fk = _mesh.getFace(indexFace);
 155                     _minl = min(_minl,Ci.getMeasure()/Fk.getMeasure());
 156                     _perimeters(i)+=Fk.getMeasure();
 157                 }
 158             }else{
 159                 _minl = min(_minl,Ci.getMeasure());
 160                 _perimeters(i)=Ci.getNumberOfFaces();
 161             }
 162         }
 163         _initialDataSet=true;
 164
 165         if(_verbose)
 166                 cout<<_perimeters<<endl;
 167 }
 168 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber)
 169 {
 170         Field VV(fileName,CELLS,fieldName,timeStepNumber,0);
 171         setInitialField(VV);
 172 }
 173 void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant)
 174 {
 175         Mesh M(fileName);
 176         Field VV("Primitive", CELLS, M, Vconstant.size());
 177
 178         for (int j = 0; j < M.getNumberOfCells(); j++) {
 179                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 180                         VV(j,i) = Vconstant[i];
 181         }
 182
 183         setInitialField(VV);
 184 }
 185 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant)
 186 {
 187         Field VV("Primitive", CELLS, M, Vconstant.getNumberOfRows());
 188
 189         for (int j = 0; j < M.getNumberOfCells(); j++) {
 190                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 191                         VV(j,i) = Vconstant(i);
 192         }
 193         setInitialField(VV);
 194 }
 195 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant)
 196 {
 197         Field VV("Primitive", CELLS, M, Vconstant.size());
 198
 199         for (int j = 0; j < M.getNumberOfCells(); j++) {
 200                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 201                         VV(j,i) = Vconstant[i];
 202         }
 203         setInitialField(VV);
 204 }
 205 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
 206                 double ymin, double ymax, int ny, string backSide, string frontSide,
 207                 double zmin, double zmax, int nz, string bottomSide, string topSide)
 208 {
 209         Mesh M;
 210         if(nDim==1){
 211                 //cout<<"coucou1 xmin="<<xmin<<", xmax= "<<xmax<< ", nx= "<<nx<<endl;
 212                 M=Mesh(xmin,xmax,nx);
 213                 //cout<<"coucou2"<<endl;
 214         }
 215         else if(nDim==2)
 216                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 217         else if(nDim==3)
 218                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 219         else{
 220                 cout<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
 221                 *_runLogFile<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
 222                 _runLogFile->close();
 223                 throw CdmathException("Space dimension nDim should be between 1 and 3");
 224         }
 225
 226         M.setGroupAtPlan(xmax,0,_precision,rightSide);
 227         M.setGroupAtPlan(xmin,0,_precision,leftSide);
 228         if(nDim>=2){
 229                 M.setGroupAtPlan(ymax,1,_precision,frontSide);
 230                 M.setGroupAtPlan(ymin,1,_precision,backSide);
 231         }
 232         if(nDim==3){
 233                 M.setGroupAtPlan(zmax,2,_precision,topSide);
 234                 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
 235         }
 236
 237         setInitialFieldConstant(M, Vconstant);
 238 }
 239 void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction)
 240 {
 241         if  (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
 242         {
 243                 *_runLogFile<<"ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes"<<endl;
 244                 _runLogFile->close();
 245                 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes");
 246         }
 247         Field VV("Primitive", CELLS, M, VV_Left.getNumberOfRows());
 248
 249         double component_value;
 250
 251         for (int j = 0; j < M.getNumberOfCells(); j++) {
 252                 if(direction==0)
 253                         component_value=M.getCell(j).x();
 254                 else if(direction==1)
 255                         component_value=M.getCell(j).y();
 256                 else if(direction==2)
 257                         component_value=M.getCell(j).z();
 258                 else{
 259                         _runLogFile->close();
 260                         throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: direction should be an integer between 0 and 2");
 261                 }
 262
 263                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 264                         if (component_value< disc_pos )
 265                                 VV(j,i) = VV_Left[i];
 266                         else
 267                                 VV(j,i) = VV_Right[i];
 268         }
 269         setInitialField(VV);
 270 }
 271 void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
 272                 double xmin, double xmax, int nx, string leftSide, string rightSide,
 273                 double ymin, double ymax, int ny, string backSide, string frontSide,
 274                 double zmin, double zmax, int nz, string bottomSide, string topSide)
 275 {
 276         Mesh M;
 277         if(nDim==1)
 278                 M=Mesh(xmin,xmax,nx);
 279         else if(nDim==2)
 280                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 281         else if(nDim==3)
 282                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 283         else
 284         {
 285                 _runLogFile->close();
 286                 throw CdmathException("Space dimension nDim should be between 1 and 3");
 287         }
 288
 289         M.setGroupAtPlan(xmax,0,_precision,rightSide);
 290         M.setGroupAtPlan(xmin,0,_precision,leftSide);
 291         if(nDim>=2){
 292                 M.setGroupAtPlan(ymax,1,_precision,frontSide);
 293                 M.setGroupAtPlan(ymin,1,_precision,backSide);
 294         }
 295         if(nDim==3){
 296                 M.setGroupAtPlan(zmax,2,_precision,topSide);
 297                 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
 298         }
 299         Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
 300         for(int i=0;i<VV_Left.size(); i++){
 301                 V_Left(i)=VV_Left[i];
 302                 V_Right(i)=VV_Right[i];
 303         }
 304         setInitialFieldStepFunction(M, V_Left, V_Right, xstep);
 305 }
 306
 307 void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center)
 308 {
 309         if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
 310         {
 311                 cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
 312                 throw CdmathException("ProblemCoreFlows::setInitialFieldSphericalStepFunction : Vector size error");
 313         }
 314         int nVar=Vout.size();
 315         int spaceDim=M.getSpaceDimension();
 316         Field VV("Primitive variables for spherical step function", CELLS, M, nVar);
 317
 318         Vector currentPoint(spaceDim);
 319         for(int i=0;i<M.getNumberOfCells();i++)
 320         {
 321                 currentPoint(0)=M.getCell(i).x();
 322                 if(spaceDim>1)
 323                 {
 324                         currentPoint(1)=M.getCell(i).y();
 325                         if(spaceDim>2)
 326                                 currentPoint(2)=M.getCell(i).z();
 327                 }
 328                 if((currentPoint-Center).norm()<radius)
 329                         for(int j=0;j<nVar;j++)
 330                                 VV(i,j)=Vin[j];
 331                 else
 332                         for(int j=0;j<nVar;j++)
 333                                 VV(i,j)=Vout[j];
 334         }
 335         setInitialField(VV);
 336 }
 337
 338 double ProblemCoreFlows::getTime()
 339 {
 340         return _time;
 341 }
 342 unsigned ProblemCoreFlows::getNbTimeStep()
 343 {
 344         return _nbTimeStep;
 345 }
 346 double ProblemCoreFlows::getCFL()
 347 {
 348         return _cfl;
 349 }
 350 double ProblemCoreFlows::getPrecision()
 351 {
 352         return _precision;
 353 }
 354 SpaceScheme ProblemCoreFlows::getSpaceScheme()
 355 {
 356         return _spaceScheme;
 357 }
 358 TimeScheme ProblemCoreFlows::getTimeScheme()
 359 {
 360         return _timeScheme;
 361 }
 362 Mesh ProblemCoreFlows::getMesh()
 363 {
 364         return _mesh;
 365 }
 366
 367 void ProblemCoreFlows::setLinearSolver(linearSolver kspType, preconditioner pcType)
 368 {
 369         //_maxPetscIts=maxIterationsPetsc;
 370         // set linear solver algorithm
 371         if (kspType==GMRES)
 372                 _ksptype = (char*)&KSPGMRES;
 373         else if (kspType==CG)
 374                 _ksptype = (char*)&KSPCG;
 375         else if (kspType==BCGS)
 376                 _ksptype = (char*)&KSPBCGS;
 377         else {
 378                 cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 379                 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 380                 _runLogFile->close();
 381                 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
 382         }
 383         // set preconditioner
 384         if (pcType == NONE)
 385                 _pctype = (char*)&PCNONE;
 386         else if (pcType ==LU)
 387                 _pctype = (char*)&PCLU;
 388         else if (pcType == ILU)
 389                 _pctype = (char*)&PCILU;
 390         else if (pcType ==CHOLESKY)
 391                 _pctype = (char*)&PCCHOLESKY;
 392         else if (pcType == ICC)
 393                 _pctype = (char*)&PCICC;
 394         else {
 395                 cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
 396                 *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
 397                 _runLogFile->close();
 398                 throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
 399         }
 400 }
 401
 402 // Description:
 403 // Cette methode lance une execution du ProblemCoreFlows
 404 // Elle peut etre utilisee si le probleme n'est couple a aucun autre.
 405 // (s'il n'a besoin d'aucun champ d'entree).
 406 // Precondition: initialize
 407 // Seule la methode terminate peut etre appelee apres
 408 bool ProblemCoreFlows::run()
 409 {
 410         if(!_initializedMemory)
 411         {
 412                 _runLogFile->close();
 413                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 414         }
 415         bool stop=false; // Does the Problem want to stop (error) ?
 416         bool ok; // Is the time interval successfully solved ?
 417         _isStationary=false;//in case of a second run with a different physics or cfl
 418
 419         cout<< "Running test case "<< _fileName<<endl;
 420
 421         _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 422         *_runLogFile<< "Running test case "<< _fileName<<endl;
 423
 424         // Time step loop
 425         while(!stop && !_isStationary &&_time<_timeMax && _nbTimeStep<_maxNbOfTimeStep)
 426         {
 427                 ok=false; // Is the time interval successfully solved ?
 428
 429                 // Guess the next time step length
 430                 _dt=computeTimeStep(stop);
 431                 if (stop){
 432                         cout << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
 433                         *_runLogFile << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
 434                         break;
 435                 }
 436                 // Loop on the time interval tries
 437                 while (!ok && !stop )
 438                 {
 439                         stop=!initTimeStep(_dt);
 440                         // Prepare the next time step
 441                         if (stop){
 442                                 cout << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
 443                                 *_runLogFile << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
 444                                 break;
 445                         }
 446                         // Solve the next time step
 447                         ok=solveTimeStep();
 448
 449                         if (!ok)   // The resolution failed, try with a new time interval.
 450                         {
 451                                 abortTimeStep();
 452                                 if(_dt>_precision){
 453                                         cout << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with cfl/2"<< endl;
 454                                         *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with cfl/2"<< endl;
 455                                         _dt*=0.5;
 456                                         _cfl*=0.5;
 457                                 }
 458                                 else{
 459                                         cout << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
 460                                         *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
 461                                         stop=true; // Impossible to solve the next time step, the Problem has given up
 462                                         break;
 463                                 }
 464                         }
 465                         else // The resolution was successful, validate and go to the next time step.
 466                         {
 467                                 validateTimeStep();
 468                                 if (_nbTimeStep%_freqSave ==0){
 469                                         cout << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl;
 470                                         *_runLogFile << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl;
 471                                 }
 472                         }
 473                 }
 474         }
 475         if(_isStationary){
 476                 cout << "Stationary state reached" <<endl;
 477                 *_runLogFile << "Stationary state reached" <<endl;
 478         }
 479         else if(_time>=_timeMax){
 480                 cout<<"Maximum time "<<_timeMax<<" reached"<<endl;
 481                 *_runLogFile<<"Maximum time "<<_timeMax<<" reached"<<endl;
 482         }
 483         else if(_nbTimeStep>=_maxNbOfTimeStep){
 484                 cout<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
 485                 *_runLogFile<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
 486         }
 487         else{
 488                 cout<<"Error problem wants to stop!"<<endl;
 489                 *_runLogFile<<"Error problem wants to stop!"<<endl;
 490         }
 491         cout << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
 492         *_runLogFile << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
 493
 494         _runLogFile->close();
 495         return !stop;
 496 }
 497
 498 void ProblemCoreFlows::displayMatrix(double *matrix, int size, string name)
 499 {
 500         cout<<name<<endl;
 501         for(int p=0; p<size; p++)
 502         {
 503                 for(int q=0; q<size; q++)
 504                 {
 505                         cout << matrix[p*size+q] << "\t";
 506                 }
 507                 cout << endl;
 508         }
 509         cout << endl;
 510 }
 511
 512 void ProblemCoreFlows::displayVector(double *vector, int size, string name)
 513 {
 514         cout<<name<<endl;
 515         for(int q=0; q<size; q++)
 516         {
 517                 cout << vector[q] << "\t";
 518         }
 519         cout << endl;
 520 }
 521 void ProblemCoreFlows::setFileName(string fileName){
 522         _fileName = fileName;
 523 }
 524
 525 void ProblemCoreFlows::setFreqSave(int freqSave){
 526         _freqSave = freqSave;
 527 }
 528 bool ProblemCoreFlows::solveTimeStep(){
 529         _NEWTON_its=0;
 530         bool converged=false, ok=true;
 531         while(!converged && ok && _NEWTON_its < _maxNewtonIts){
 532                 ok=iterateTimeStep(converged);//resolution du systeme lineaire si schema implicite
 533
 534                 if(_timeScheme == Implicit && _nbTimeStep%_freqSave ==0)//To monitor the convergence of the newton scheme
 535                 {
 536                         cout << "\n Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << " : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
 537                         *_runLogFile<< "\n Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << " : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
 538
 539                         if(_conditionNumber)
 540                         {
 541                                 PetscReal sv_max, sv_min;
 542                                 KSPComputeExtremeSingularValues(_ksp, &sv_max, &sv_min);
 543                                 cout<<" singular value max = " << sv_max <<" singular value min = " << sv_min <<" condition number = " << sv_max/sv_min <<endl;
 544                                 *_runLogFile<<" singular value max = " << sv_max <<" singular value min = " << sv_min <<" condition number = " << sv_max/sv_min <<endl;
 545                         }
 546                 }
 547                 _NEWTON_its++;
 548         }
 549         if(!converged){
 550                 if(_NEWTON_its >= _maxNewtonIts){
 551                         cout << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
 552                         *_runLogFile << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
 553                 }
 554                 else if(!ok){
 555                         cout<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
 556                         *_runLogFile<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
 557                 }
 558         }
 559         else if(_timeScheme == Implicit && _nbTimeStep%_freqSave ==0)
 560         {
 561                 cout<<endl;
 562                 cout << "Nombre d'iterations de Newton "<< _NEWTON_its << ", Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl;
 563                 *_runLogFile <<endl;
 564                 *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its << " Variation ||Un+1-Un||= "<< _erreur_rel<<endl;
 565                 *_runLogFile << "Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl;
 566                 _MaxIterLinearSolver = 0;
 567         }
 568
 569         return converged;
 570 }
 571 ProblemCoreFlows::~ProblemCoreFlows()
 572 {
 573         /*
 574         PetscBool petscInitialized;
 575         PetscInitialized(&petscInitialized);
 576         if(petscInitialized)
 577                 PetscFinalize();
 578          */
 579         delete _runLogFile;
 580 }