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