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 "SparseMatrixPetsc.hxx"
  15
  16 #include <limits.h>
  17 #include <unistd.h>
  18
  19 using namespace std;
  20
  21 ProblemCoreFlows::ProblemCoreFlows()
  22 {
  23         PetscBool petscInitialized;
  24         PetscInitialized(&petscInitialized);
  25         if(!petscInitialized)
  26                 PetscInitialize(NULL,NULL,0,0);
  27         _dt = 0;
  28         _time = 0;
  29         _nbTimeStep=0;
  30         _cfl=0;
  31         _timeMax =0.;
  32         _maxNbOfTimeStep = 0;
  33         _precision=1.e-6;
  34         _precision_Newton=_precision;
  35         _erreur_rel= 0;
  36         _isStationary=false;
  37         _Ndim=0;
  38         _minl=0;
  39         _neibMaxNb=0;
  40         _fileName = "myCoreFlowsProblem";
  41         _freqSave = 1;
  42         _initialDataSet=false;
  43         _initializedMemory=false;
  44         _restartWithNewTimeScheme=false;
  45         _restartWithNewFileName=false;
  46         _timeScheme=Explicit;
  47         _wellBalancedCorrection=false;
  48     _FECalculation=false;
  49         _maxPetscIts=50;
  50         _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
  51         _maxNewtonIts=50;
  52         _NEWTON_its=0;
  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;
  59         _heatSource=0;
  60         _verbose = false;
  61         _system = false;
  62         _conditionNumber=false;
  63         _maxvp=0;
  64         _runLogFile=new ofstream;
  65
  66         //extracting current directory
  67         char result[ PATH_MAX ];
  68         getcwd(result, PATH_MAX );
  69         _path=string( result );
  70         _saveFormat=VTK;
  71 }
  72
  73 TimeScheme ProblemCoreFlows::getTimeScheme()
  74 {
  75         return _timeScheme;
  76 }
  77
  78 void ProblemCoreFlows::setTimeScheme(TimeScheme timeScheme)
  79 {
  80         if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
  81                 _restartWithNewTimeScheme=true;
  82         _timeScheme = timeScheme;
  83 }
  84
  85 bool ProblemCoreFlows::isStationary() const
  86 {
  87         return _isStationary;
  88 }
  89
  90 double ProblemCoreFlows::presentTime() const
  91 {
  92         return _time;
  93 }
  94 void ProblemCoreFlows::setTimeMax(double timeMax){
  95         _timeMax = timeMax;
  96 }
  97 void ProblemCoreFlows::setPresentTime (double time)
  98 {
  99         _time=time;
 100 }
 101 void ProblemCoreFlows::setMaxNbOfTimeStep(int maxNbOfTimeStep){
 102         _maxNbOfTimeStep = maxNbOfTimeStep;
 103 }
 104 void ProblemCoreFlows::setCFL(double cfl)
 105 {
 106         _cfl=cfl;
 107 }
 108 void ProblemCoreFlows::setPrecision(double precision)
 109 {
 110         _precision=precision;
 111 }
 112 void ProblemCoreFlows::setInitialField(const Field &VV)
 113 {
 114
 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");
 119         }
 120         if(_nVar!=VV.getNumberOfComponents())
 121         {
 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");
 125         }
 126
 127         _VV=VV;
 128         _time=_VV.getTime();
 129         _mesh=_VV.getMesh();
 130         _Nmailles = _mesh.getNumberOfCells();
 131         _Nnodes =   _mesh.getNumberOfNodes();
 132         _Nfaces =   _mesh.getNumberOfFaces();
 133         _perimeters=Field("Perimeters", CELLS, _mesh,1);
 134
 135         // find _minl and maximum nb of neibourghs
 136         _minl  = INFINITY;
 137         int nbNeib,indexFace;
 138         Cell Ci;
 139         Face Fk;
 140
 141         if(_verbose)
 142                 cout<<"Computing cell perimeters and mesh minimal diameter"<<endl;
 143
 144     if(VV.getTypeOfField()==NODES)
 145     {
 146         _minl = _mesh.getMaxNbNeighbours(NODES);
 147         _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 148     }
 149     else
 150         for (int i=0; i<_mesh.getNumberOfCells(); i++){
 151             Ci = _mesh.getCell(i);
 152             //Detect mesh with junction
 153             nbNeib=0;
 154             for(int j=0; j<Ci.getNumberOfFaces(); j++){
 155                 Fk=_mesh.getFace(Ci.getFacesId()[j]);
 156                 nbNeib+=Fk.getNumberOfCells()-1;
 157             }
 158             if(nbNeib>_neibMaxNb)
 159                 _neibMaxNb=nbNeib;
 160             //Compute mesh data
 161             if (_Ndim > 1){
 162                 _perimeters(i)=0;
 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();
 168                 }
 169             }else{
 170                 _minl = min(_minl,Ci.getMeasure());
 171                 _perimeters(i)=Ci.getNumberOfFaces();
 172             }
 173         }
 174         _initialDataSet=true;
 175
 176         if(_verbose)
 177                 cout<<_perimeters<<endl;
 178 }
 179 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber)
 180 {
 181         Field VV(fileName,CELLS,fieldName,timeStepNumber,0);
 182         setInitialField(VV);
 183 }
 184 void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant)
 185 {
 186         Mesh M(fileName);
 187         Field VV("Primitive", CELLS, M, Vconstant.size());
 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
 194         setInitialField(VV);
 195 }
 196 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant)
 197 {
 198         Field VV("Primitive", CELLS, M, Vconstant.getNumberOfRows());
 199
 200         for (int j = 0; j < M.getNumberOfCells(); j++) {
 201                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 202                         VV(j,i) = Vconstant(i);
 203         }
 204         setInitialField(VV);
 205 }
 206 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant)
 207 {
 208         Field VV("Primitive", CELLS, M, Vconstant.size());
 209
 210         for (int j = 0; j < M.getNumberOfCells(); j++) {
 211                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 212                         VV(j,i) = Vconstant[i];
 213         }
 214         setInitialField(VV);
 215 }
 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)
 219 {
 220         Mesh M;
 221         if(nDim==1){
 222                 //cout<<"coucou1 xmin="<<xmin<<", xmax= "<<xmax<< ", nx= "<<nx<<endl;
 223                 M=Mesh(xmin,xmax,nx);
 224                 //cout<<"coucou2"<<endl;
 225         }
 226         else if(nDim==2)
 227                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 228         else if(nDim==3)
 229                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 230         else{
 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");
 235         }
 236
 237         M.setGroupAtPlan(xmax,0,_precision,rightSide);
 238         M.setGroupAtPlan(xmin,0,_precision,leftSide);
 239         if(nDim>=2){
 240                 M.setGroupAtPlan(ymax,1,_precision,frontSide);
 241                 M.setGroupAtPlan(ymin,1,_precision,backSide);
 242         }
 243         if(nDim==3){
 244                 M.setGroupAtPlan(zmax,2,_precision,topSide);
 245                 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
 246         }
 247
 248         setInitialFieldConstant(M, Vconstant);
 249 }
 250 void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction)
 251 {
 252         if  (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
 253         {
 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");
 257         }
 258         Field VV("Primitive", CELLS, M, VV_Left.getNumberOfRows());
 259
 260         double component_value;
 261
 262         for (int j = 0; j < M.getNumberOfCells(); j++) {
 263                 if(direction==0)
 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();
 269                 else{
 270                         _runLogFile->close();
 271                         throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: direction should be an integer between 0 and 2");
 272                 }
 273
 274                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 275                         if (component_value< disc_pos )
 276                                 VV(j,i) = VV_Left[i];
 277                         else
 278                                 VV(j,i) = VV_Right[i];
 279         }
 280         setInitialField(VV);
 281 }
 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)
 286 {
 287         Mesh M;
 288         if(nDim==1)
 289                 M=Mesh(xmin,xmax,nx);
 290         else if(nDim==2)
 291                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 292         else if(nDim==3)
 293                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 294         else
 295         {
 296                 _runLogFile->close();
 297                 throw CdmathException("Space dimension nDim should be between 1 and 3");
 298         }
 299
 300         M.setGroupAtPlan(xmax,0,_precision,rightSide);
 301         M.setGroupAtPlan(xmin,0,_precision,leftSide);
 302         if(nDim>=2){
 303                 M.setGroupAtPlan(ymax,1,_precision,frontSide);
 304                 M.setGroupAtPlan(ymin,1,_precision,backSide);
 305         }
 306         if(nDim==3){
 307                 M.setGroupAtPlan(zmax,2,_precision,topSide);
 308                 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
 309         }
 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];
 314         }
 315         setInitialFieldStepFunction(M, V_Left, V_Right, xstep);
 316 }
 317
 318 void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center)
 319 {
 320         if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
 321         {
 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");
 324         }
 325         int nVar=Vout.size();
 326         int spaceDim=M.getSpaceDimension();
 327         Field VV("Primitive variables for spherical step function", CELLS, M, nVar);
 328
 329         Vector currentPoint(spaceDim);
 330         for(int i=0;i<M.getNumberOfCells();i++)
 331         {
 332                 currentPoint(0)=M.getCell(i).x();
 333                 if(spaceDim>1)
 334                 {
 335                         currentPoint(1)=M.getCell(i).y();
 336                         if(spaceDim>2)
 337                                 currentPoint(2)=M.getCell(i).z();
 338                 }
 339                 if((currentPoint-Center).norm()<radius)
 340                         for(int j=0;j<nVar;j++)
 341                                 VV(i,j)=Vin[j];
 342                 else
 343                         for(int j=0;j<nVar;j++)
 344                                 VV(i,j)=Vout[j];
 345         }
 346         setInitialField(VV);
 347 }
 348
 349 double ProblemCoreFlows::getTime()
 350 {
 351         return _time;
 352 }
 353 unsigned ProblemCoreFlows::getNbTimeStep()
 354 {
 355         return _nbTimeStep;
 356 }
 357 double ProblemCoreFlows::getCFL()
 358 {
 359         return _cfl;
 360 }
 361 double ProblemCoreFlows::getPrecision()
 362 {
 363         return _precision;
 364 }
 365 Mesh ProblemCoreFlows::getMesh()
 366 {
 367         return _mesh;
 368 }
 369
 370 void ProblemCoreFlows::setLinearSolver(linearSolver kspType, preconditioner pcType)
 371 {
 372         //_maxPetscIts=maxIterationsPetsc;
 373         // set linear solver algorithm
 374         if (kspType==GMRES)
 375                 _ksptype = (char*)&KSPGMRES;
 376         else if (kspType==CG)
 377                 _ksptype = (char*)&KSPCG;
 378         else if (kspType==BCGS)
 379                 _ksptype = (char*)&KSPBCGS;
 380         else {
 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 !!!");
 385         }
 386         // set preconditioner
 387         if (pcType == NONE)
 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;
 397         else {
 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 !!!" );
 402         }
 403 }
 404
 405 // Description:
 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()
 412 {
 413         if(!_initializedMemory)
 414         {
 415                 _runLogFile->close();
 416                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 417         }
 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
 421
 422         cout<< "Running test case "<< _fileName<<endl;
 423
 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;
 426
 427         // Time step loop
 428         while(!stop && !_isStationary &&_time<_timeMax && _nbTimeStep<_maxNbOfTimeStep)
 429         {
 430                 ok=false; // Is the time interval successfully solved ?
 431
 432                 // Guess the next time step length
 433                 _dt=computeTimeStep(stop);
 434                 if (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;
 437                         break;
 438                 }
 439                 // Loop on the time interval tries
 440                 while (!ok && !stop )
 441                 {
 442                         stop=!initTimeStep(_dt);
 443                         // Prepare the next time step
 444                         if (stop){
 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;
 447                                 break;
 448                         }
 449                         // Solve the next time step
 450                         ok=solveTimeStep();
 451
 452                         if (!ok)   // The resolution failed, try with a new time interval.
 453                         {
 454                                 if(_dt>_precision){
 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);
 463                                 }
 464                                 else{
 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
 468                                         break;
 469                                 }
 470                         }
 471                         else // The resolution was successful, validate and go to the next time step.
 472                         {
 473                                 validateTimeStep();
 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;
 477                                 }
 478                         }
 479                 }
 480         }
 481         if(_isStationary){
 482                 cout << "Stationary state reached" <<endl;
 483                 *_runLogFile << "Stationary state reached" <<endl;
 484         }
 485         else if(_time>=_timeMax){
 486                 cout<<"Maximum time "<<_timeMax<<" reached"<<endl;
 487                 *_runLogFile<<"Maximum time "<<_timeMax<<" reached"<<endl;
 488         }
 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;
 492         }
 493         else{
 494                 cout<<"Error problem wants to stop!"<<endl;
 495                 *_runLogFile<<"Error problem wants to stop!"<<endl;
 496         }
 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;
 499
 500         _runLogFile->close();
 501         return !stop;
 502 }
 503
 504 void ProblemCoreFlows::displayMatrix(double *matrix, int size, string name)
 505 {
 506         cout<<name<<endl;
 507         for(int p=0; p<size; p++)
 508         {
 509                 for(int q=0; q<size; q++)
 510                 {
 511                         cout << matrix[p*size+q] << "\t";
 512                 }
 513                 cout << endl;
 514         }
 515         cout << endl;
 516 }
 517
 518 void ProblemCoreFlows::displayVector(double *vector, int size, string name)
 519 {
 520         cout<<name<<endl;
 521         for(int q=0; q<size; q++)
 522         {
 523                 cout << vector[q] << "\t";
 524         }
 525         cout << endl;
 526 }
 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;
 531 }
 532
 533 void ProblemCoreFlows::setFreqSave(int freqSave){
 534         _freqSave = freqSave;
 535 }
 536 bool ProblemCoreFlows::solveTimeStep(){
 537         _NEWTON_its=0;
 538         bool converged=false, ok=true;
 539         while(!converged && ok && _NEWTON_its < _maxNewtonIts){
 540                 ok=iterateTimeStep(converged);//resolution du systeme lineaire si schema implicite
 541
 542                 if(_timeScheme == Implicit && _nbTimeStep%_freqSave ==0)//To monitor the convergence of the newton scheme
 543                 {
 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;
 546
 547                         if(_conditionNumber)
 548                         {
 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;
 553                         }
 554                 }
 555                 _NEWTON_its++;
 556         }
 557         if(!converged){
 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;
 561                 }
 562                 else if(!ok){
 563                         cout<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
 564                         *_runLogFile<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
 565                 }
 566         }
 567         else if(_timeScheme == Implicit && _nbTimeStep%_freqSave ==0)
 568         {
 569                 cout<<endl;
 570                 cout << "Nombre d'iterations de Newton "<< _NEWTON_its << ", Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl;
 571                 *_runLogFile <<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;
 575         }
 576
 577         return converged;
 578 }
 579 ProblemCoreFlows::~ProblemCoreFlows()
 580 {
 581         /*
 582         PetscBool petscInitialized;
 583         PetscInitialized(&petscInitialized);
 584         if(petscInitialized)
 585                 PetscFinalize();
 586          */
 587         delete _runLogFile;
 588 }
 589
 590 double
 591 ProblemCoreFlows::getConditionNumber(bool isSingular, double tol) const
 592 {
 593   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 594   return A.getConditionNumber( isSingular, tol);
 595 }
 596 std::vector< double >
 597 ProblemCoreFlows::getEigenvalues(int nev, EPSWhich which, double tol) const
 598 {
 599   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 600   return A.getEigenvalues( nev, which, tol);
 601 }
 602 std::vector< Vector >
 603 ProblemCoreFlows::getEigenvectors(int nev, EPSWhich which, double tol) const
 604 {
 605   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 606   return A.getEigenvectors( nev, which, tol);
 607 }
 608 Field
 609 ProblemCoreFlows::getEigenvectorsField(int nev, EPSWhich which, double tol) const
 610 {
 611   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 612   MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
 613   Field my_eigenfield;
 614
 615   if(_FECalculation)
 616     my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
 617   else
 618     my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
 619
 620   my_eigenfield.setFieldByDataArrayDouble(d);
 621
 622   return my_eigenfield;
 623 }
 624
 625 Field
 626 ProblemCoreFlows::getUnknownField() const
 627 {
 628         if(!_initializedMemory)
 629         {
 630                 _runLogFile->close();
 631                 throw CdmathException("ProblemCoreFlows::getUnknownField() call initialize() first");
 632         }
 633         return _VV;
 634 }