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(MPI_Comm comm)
  22 {
  23         /* Initialisation of PETSC */
  24         //check if PETSC is already initialised
  25         PetscBool petscInitialized;
  26         PetscInitialized(&petscInitialized);
  27         if(!petscInitialized)
  28         {//check if MPI is already initialised
  29                 int mpiInitialized;
  30                 MPI_Initialized(&mpiInitialized);
  31                 if(mpiInitialized)
  32                         PETSC_COMM_WORLD = comm;
  33                 PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
  34         }
  35         MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
  36         MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
  37         PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
  38         PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
  39         PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
  40
  41         if(_mpi_size>1)
  42         {
  43                 PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n");
  44                 PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output fields remain sequential ----\n");
  45                 PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix operations are done in parallel thanks to PETSc----\n");
  46                 PetscPrintf(PETSC_COMM_WORLD,"---- Processor %d is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n",_mpi_rank);
  47         }
  48
  49         /* Numerical parameter */
  50         _dt = 0;
  51         _time = 0;
  52         _nbTimeStep=0;
  53         _cfl=0;
  54         _timeMax =0.;
  55         _maxNbOfTimeStep = 0;
  56         _precision=1.e-6;
  57         _precision_Newton=_precision;
  58         _erreur_rel= 0;
  59         _isStationary=false;
  60         _stationaryMode=false;
  61         _timeScheme=Explicit;
  62         _wellBalancedCorrection=false;
  63     _FECalculation=false;
  64         _nonLinearSolver = Newton_SOLVERLAB;
  65         _conditionNumber=false;
  66         _maxvp=0;
  67         _runLogFile=new ofstream;
  68
  69         /* Monitoring of simulation */
  70         _restartWithNewTimeScheme=false;
  71         _restartWithNewFileName=false;
  72         _fileName = "myCoreFlowsProblem";
  73         _freqSave = 1;
  74         _verbose = false;
  75         _system = false;
  76
  77         /* Mesh parameters */
  78         _Ndim=0;
  79         _minl=0;
  80         _neibMaxNbCells=0;
  81         _neibMaxNbNodes=0;
  82
  83         /* Memory and restart */
  84         _initialDataSet=false;
  85         _initializedMemory=false;
  86
  87         /* PETSc and linear solver parameters */
  88         _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
  89         _NEWTON_its=0;
  90         _maxPetscIts=50;
  91         _maxNewtonIts=50;
  92         _PetscIts=0;//the number of iterations of the linear solver
  93         _ksptype = (char*)&KSPGMRES;
  94         if( _mpi_size>1)
  95                 _pctype = (char*)&PCNONE;
  96         else
  97                 _pctype = (char*)&PCLU;
  98
  99         /* Physical parameters */
 100         _heatPowerFieldSet=false;
 101         _heatTransfertCoeff=0;
 102         _heatSource=0;
 103
 104         //extracting current directory
 105         char result[ PATH_MAX ];
 106         getcwd(result, PATH_MAX );
 107         _path=string( result );
 108         _saveFormat=VTK;
 109 }
 110
 111 TimeScheme ProblemCoreFlows::getTimeScheme()
 112 {
 113         return _timeScheme;
 114 }
 115
 116 void ProblemCoreFlows::setTimeScheme(TimeScheme timeScheme)
 117 {
 118         if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
 119                 _restartWithNewTimeScheme=true;
 120         _timeScheme = timeScheme;
 121 }
 122
 123 bool ProblemCoreFlows::isStationary() const
 124 {
 125         return _isStationary;
 126 }
 127
 128 double ProblemCoreFlows::presentTime() const
 129 {
 130         return _time;
 131 }
 132 void ProblemCoreFlows::setTimeMax(double timeMax){
 133         _timeMax = timeMax;
 134 }
 135 void ProblemCoreFlows::resetTime (double time)
 136 {
 137         _time=time;
 138 }
 139 void ProblemCoreFlows::setMaxNbOfTimeStep(int maxNbOfTimeStep){
 140         _maxNbOfTimeStep = maxNbOfTimeStep;
 141 }
 142 void ProblemCoreFlows::setCFL(double cfl)
 143 {
 144         _cfl=cfl;
 145 }
 146 void ProblemCoreFlows::setPrecision(double precision)
 147 {
 148         _precision=precision;
 149 }
 150 void ProblemCoreFlows::setInitialField(const Field &VV)
 151 {
 152         if(_Ndim != VV.getSpaceDimension()){
 153                 *_runLogFile<<"ProblemCoreFlows::setInitialField: mesh has incorrect space dimension"<<endl;
 154                 _runLogFile->close();
 155                 throw CdmathException("ProblemCoreFlows::setInitialField: mesh has incorrect space dimension");
 156         }
 157         if(_nVar!=VV.getNumberOfComponents())
 158         {
 159                 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect number of components"<<endl;
 160                 _runLogFile->close();
 161                 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect number of components");
 162         }
 163         if(_FECalculation && VV.getTypeOfField()!=NODES)
 164         {
 165                 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method"<<endl;
 166                 _runLogFile->close();
 167                 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method");
 168         }
 169         else if(!_FECalculation && VV.getTypeOfField()==NODES)
 170         {
 171                 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method"<<endl;
 172                 _runLogFile->close();
 173                 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method");
 174         }
 175
 176         _VV=VV;
 177         _VV.setName("SOLVERLAB results");
 178         _time=_VV.getTime();
 179         _mesh=_VV.getMesh();
 180
 181         _initialDataSet=true;
 182
 183         //Mesh data
 184         _Nmailles = _mesh.getNumberOfCells();
 185         _Nnodes =   _mesh.getNumberOfNodes();
 186         _Nfaces =   _mesh.getNumberOfFaces();
 187         _perimeters=Field("Perimeters", CELLS, _mesh,1);
 188
 189         // find _minl (delta x) and maximum nb of neibourghs
 190         _minl  = INFINITY;
 191         int nbNeib,indexFace;
 192         Cell Ci;
 193         Face Fk;
 194
 195         if(_verbose)
 196                 PetscPrintf(PETSC_COMM_SELF,"Processor %d Computing cell perimeters and mesh minimal diameter\n", _mpi_rank);
 197
 198         //Compute the maximum number of neighbours for nodes or cells
 199     if(VV.getTypeOfField()==NODES)
 200         _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
 201     else
 202     {
 203         _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
 204
 205                 //Compute Delta x and the cell perimeters
 206                 for (int i=0; i<_mesh.getNumberOfCells(); i++){
 207                         Ci = _mesh.getCell(i);
 208                         if (_Ndim > 1){
 209                                 _perimeters(i)=0;
 210                                 for (int k=0 ; k<Ci.getNumberOfFaces() ; k++){
 211                                         indexFace=Ci.getFacesId()[k];
 212                                         Fk = _mesh.getFace(indexFace);
 213                                         _minl = min(_minl,Ci.getMeasure()/Fk.getMeasure());
 214                                         _perimeters(i)+=Fk.getMeasure();
 215                                 }
 216                         }else{
 217                                 _minl = min(_minl,Ci.getMeasure());
 218                                 _perimeters(i)=Ci.getNumberOfFaces();
 219                         }
 220                 }
 221                 if(_verbose)
 222                         cout<<_perimeters<<endl;
 223         }
 224
 225         /*** MPI distribution of parameters ***/
 226         MPI_Allreduce(&_initialDataSet, &_initialDataSet, 1, MPIU_BOOL, MPI_LOR, PETSC_COMM_WORLD);
 227
 228         int nbVoisinsMax;//Mettre en attribut ?
 229         if(!_FECalculation){
 230                 MPI_Bcast(&_Nmailles      , 1, MPI_INT, 0, PETSC_COMM_WORLD);
 231                 MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 232                 nbVoisinsMax = _neibMaxNbCells;
 233         }
 234         else{
 235                 MPI_Bcast(&_Nnodes        , 1, MPI_INT, 0, PETSC_COMM_WORLD);
 236                 MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
 237                 nbVoisinsMax = _neibMaxNbNodes;
 238         }
 239     _d_nnz = (nbVoisinsMax+1)*_nVar;
 240     _o_nnz =  nbVoisinsMax   *_nVar;
 241 }
 242
 243 //Function needed because swig of enum EntityType fails
 244 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type)
 245 {
 246         if(_FECalculation && field_support_type!= NODES)
 247                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 248
 249         Field VV;
 250
 251         switch(field_support_type)
 252         {
 253         case CELLS:
 254                 VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
 255                 break;
 256         case NODES:
 257                 VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
 258                 break;
 259         case FACES:
 260                 VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
 261                 break;
 262         default:
 263                 std::ostringstream message;
 264                 message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
 265                 throw CdmathException(message.str().c_str());
 266         }
 267
 268         setInitialField(VV);
 269 }
 270 //Function needed because swig of enum EntityType fails
 271 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
 272                 double ymin, double ymax, int ny, string backSide, string frontSide,
 273                 double zmin, double zmax, int nz, string bottomSide, string topSide, int field_support_type)
 274 {
 275         if(_FECalculation && field_support_type!= NODES)
 276                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 277
 278         EntityType typeField;
 279
 280         switch(field_support_type)
 281         {
 282         case CELLS:
 283                 typeField=CELLS;
 284                 break;
 285         case NODES:
 286                 typeField=NODES;
 287                 break;
 288         case FACES:
 289                 typeField=FACES;
 290                 break;
 291         default:
 292                 std::ostringstream message;
 293                 message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
 294                 throw CdmathException(message.str().c_str());
 295         }
 296
 297         setInitialFieldConstant( nDim, Vconstant, xmin, xmax, nx, leftSide, rightSide, ymin, ymax, ny, backSide, frontSide, zmin, zmax, nz, bottomSide, topSide, typeField);
 298 }
 299 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, EntityType typeField)
 300 {
 301         if(_FECalculation && typeField!= NODES)
 302                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 303
 304         Field VV(fileName, typeField, fieldName, timeStepNumber, order, meshLevel);
 305
 306         setInitialField(VV);
 307 }
 308 void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant, EntityType typeField)
 309 {
 310         if(_FECalculation && typeField!= NODES)
 311                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 312         Mesh M(fileName);
 313         Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
 314
 315         for (int j = 0; j < VV.getNumberOfElements(); j++) {
 316                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 317                         VV(j,i) = Vconstant[i];
 318         }
 319
 320         setInitialField(VV);
 321 }
 322 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant, EntityType typeField)
 323 {
 324         if(_FECalculation && typeField!= NODES)
 325                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 326
 327         Field VV("SOLVERLAB results", typeField, M, Vconstant.getNumberOfRows());
 328
 329         for (int j = 0; j < VV.getNumberOfElements(); j++) {
 330                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 331                         VV(j,i) = Vconstant(i);
 332         }
 333         setInitialField(VV);
 334 }
 335 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant, EntityType typeField)
 336 {
 337         if(_FECalculation && typeField!= NODES)
 338                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 339
 340         Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
 341
 342         for (int j = 0; j < VV.getNumberOfElements(); j++) {
 343                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 344                         VV(j,i) = Vconstant[i];
 345         }
 346         setInitialField(VV);
 347 }
 348 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
 349                 double ymin, double ymax, int ny, string backSide, string frontSide,
 350                 double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
 351 {
 352         Mesh M;
 353         if(nDim==1)
 354                 M=Mesh(xmin,xmax,nx);
 355         else if(nDim==2)
 356                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 357         else if(nDim==3)
 358                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 359         else{
 360                 PetscPrintf(PETSC_COMM_WORLD,"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3\n");
 361                 *_runLogFile<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
 362                 _runLogFile->close();
 363                 throw CdmathException("Space dimension nDim should be between 1 and 3");
 364         }
 365
 366         M.setGroupAtPlan(xmax,0,_precision,rightSide);
 367         M.setGroupAtPlan(xmin,0,_precision,leftSide);
 368         if(nDim>=2){
 369                 M.setGroupAtPlan(ymax,1,_precision,frontSide);
 370                 M.setGroupAtPlan(ymin,1,_precision,backSide);
 371         }
 372         if(nDim==3){
 373                 M.setGroupAtPlan(zmax,2,_precision,topSide);
 374                 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
 375         }
 376
 377         setInitialFieldConstant(M, Vconstant, typeField);
 378 }
 379 void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction, EntityType typeField)
 380 {
 381         if(_FECalculation && typeField!= NODES)
 382                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 383
 384         if  (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
 385         {
 386                 *_runLogFile<<"ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes"<<endl;
 387                 _runLogFile->close();
 388                 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes");
 389         }
 390         Field VV("SOLVERLAB results", typeField, M, VV_Left.getNumberOfRows());
 391
 392         double component_value;
 393
 394         for (int j = 0; j < VV.getNumberOfElements(); j++)
 395         {
 396                 if(direction<=2)
 397                         component_value=VV.getElementComponent(j, direction);
 398                 else
 399                 {
 400                         PetscPrintf(PETSC_COMM_WORLD,"Error : space dimension is %d,  direction asked is \%d \n",M.getSpaceDimension(),direction);
 401                         _runLogFile->close();
 402                         throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: direction should be an integer between 0 and 2");
 403                 }
 404
 405                 for (int i=0; i< VV.getNumberOfComponents(); i++)
 406                         if (component_value< disc_pos )
 407                                 VV(j,i) = VV_Left[i];
 408                         else
 409                                 VV(j,i) = VV_Right[i];
 410         }
 411         setInitialField(VV);
 412 }
 413 void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
 414                 double xmin, double xmax, int nx, string leftSide, string rightSide,
 415                 double ymin, double ymax, int ny, string backSide, string frontSide,
 416                 double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
 417 {
 418         Mesh M;
 419         if(nDim==1)
 420                 M=Mesh(xmin,xmax,nx);
 421         else if(nDim==2)
 422                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 423         else if(nDim==3)
 424                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 425         else
 426         {
 427                 _runLogFile->close();
 428                 throw CdmathException("Space dimension nDim should be between 1 and 3");
 429         }
 430
 431         M.setGroupAtPlan(xmax,0,_precision,rightSide);
 432         M.setGroupAtPlan(xmin,0,_precision,leftSide);
 433         if(nDim>=2){
 434                 M.setGroupAtPlan(ymax,1,_precision,frontSide);
 435                 M.setGroupAtPlan(ymin,1,_precision,backSide);
 436         }
 437         if(nDim==3){
 438                 M.setGroupAtPlan(zmax,2,_precision,topSide);
 439                 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
 440         }
 441         Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
 442         for(int i=0;i<VV_Left.size(); i++){
 443                 V_Left(i)=VV_Left[i];
 444                 V_Right(i)=VV_Right[i];
 445         }
 446         setInitialFieldStepFunction(M, V_Left, V_Right, xstep, typeField);
 447 }
 448
 449 void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center, EntityType typeField)
 450 {
 451         if(_FECalculation && typeField!= NODES)
 452                 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
 453
 454         if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
 455         {
 456                 PetscPrintf(PETSC_COMM_WORLD,"Vout.size() = %d, Vin.size()= %d, Center.size() = %d, M.getSpaceDim = %d \n",Vout.size(),Vin.size(),Center.size(), M.getSpaceDimension());
 457                 throw CdmathException("ProblemCoreFlows::setInitialFieldSphericalStepFunction : Vector size error");
 458         }
 459         int nVar=Vout.size();
 460         int spaceDim=M.getSpaceDimension();
 461         Field VV("Primitive variables for spherical step function", typeField, M, nVar);
 462
 463         Vector currentPoint(spaceDim);
 464         for(int i=0;i<VV.getNumberOfElements();i++)
 465         {
 466                 currentPoint(0)=VV.getElementComponent(i,0);
 467                 if(spaceDim>1)
 468                 {
 469                         currentPoint(1)=VV.getElementComponent(i,1);
 470                         if(spaceDim>2)
 471                                 currentPoint(2)=VV.getElementComponent(i,2);
 472                 }
 473                 if((currentPoint-Center).norm()<radius)
 474                         for(int j=0;j<nVar;j++)
 475                                 VV(i,j)=Vin[j];
 476                 else
 477                         for(int j=0;j<nVar;j++)
 478                                 VV(i,j)=Vout[j];
 479         }
 480         setInitialField(VV);
 481 }
 482
 483 double ProblemCoreFlows::getTime()
 484 {
 485         return _time;
 486 }
 487 unsigned ProblemCoreFlows::getNbTimeStep()
 488 {
 489         return _nbTimeStep;
 490 }
 491 double ProblemCoreFlows::getCFL()
 492 {
 493         return _cfl;
 494 }
 495 double ProblemCoreFlows::getPrecision()
 496 {
 497         return _precision;
 498 }
 499 Mesh ProblemCoreFlows::getMesh()
 500 {
 501         return _mesh;
 502 }
 503
 504 void ProblemCoreFlows::setLinearSolver(linearSolver kspType, preconditioner pcType, double maxIts)
 505 {
 506         _maxPetscIts=maxIts;
 507         // set linear solver algorithm
 508         if (kspType==GMRES)
 509                 _ksptype = (char*)&KSPGMRES;
 510         else if (kspType==CG)
 511                 _ksptype = (char*)&KSPCG;
 512         else if (kspType==BCGS)
 513                 _ksptype = (char*)&KSPBCGS;
 514         else {
 515                 PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!\n");
 516                 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
 517                 _runLogFile->close();
 518                 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
 519         }
 520         // set preconditioner
 521         if (pcType == NONE)
 522                 _pctype = (char*)&PCNONE;
 523         else if (pcType ==LU)
 524                 _pctype = (char*)&PCLU;
 525         else if (pcType == ILU)
 526                 _pctype = (char*)&PCILU;
 527         else if (pcType ==CHOLESKY)
 528                 _pctype = (char*)&PCCHOLESKY;
 529         else if (pcType == ICC)
 530                 _pctype = (char*)&PCICC;
 531         else {
 532                 PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!\n");
 533                 *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
 534                 _runLogFile->close();
 535                 throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
 536         }
 537 }
 538
 539 void ProblemCoreFlows::setNewtonSolver(double precision, int iterations, nonLinearSolver solverName)
 540 {
 541         _maxNewtonIts=iterations;
 542         _precision_Newton=precision;
 543         _nonLinearSolver=solverName;
 544 }
 545
 546
 547 // Description:
 548 // Cette methode lance une execution du ProblemCoreFlows
 549 // Elle peut etre utilisee si le probleme n'est couple a aucun autre.
 550 // (s'il n'a besoin d'aucun champ d'entree).
 551 // Precondition: initialize
 552 // Seule la methode terminate peut etre appelée apres
 553 bool ProblemCoreFlows::run()
 554 {
 555         if(!_initializedMemory)
 556         {
 557                 _runLogFile->close();
 558                 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
 559         }
 560         bool stop=false; // Does the Problem want to stop (error) ?
 561         bool ok; // Is the time interval successfully solved ?
 562         _isStationary=false;//in case of a second run with a different physics or cfl
 563
 564         PetscPrintf(PETSC_COMM_WORLD,"Running test case %s\n",_fileName);
 565
 566         _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
 567         *_runLogFile<< "Running test case "<< _fileName<<endl;
 568
 569         // Time step loop
 570         while(!stop && !_isStationary &&_time<_timeMax && _nbTimeStep<_maxNbOfTimeStep)
 571         {
 572                 ok=false; // Is the time interval successfully solved ?
 573
 574                 // Guess the next time step length
 575                 _dt=computeTimeStep(stop);
 576                 if (stop){
 577                         PetscPrintf(PETSC_COMM_WORLD,"Failed computing time step %d, time = %.2f, dt= %.2f, stopping calculation",_nbTimeStep,_time,_dt);
 578                         *_runLogFile << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
 579                         break;
 580                 }
 581                 // Loop on the time interval tries
 582                 while (!ok && !stop )
 583                 {
 584                         stop=!initTimeStep(_dt);
 585                         // Prepare the next time step
 586                         if (stop){
 587                                 PetscPrintf(PETSC_COMM_WORLD,"Failed initializing time step %d, time = %.2f, dt= %.2f, stopping calculation",_nbTimeStep,_time,_dt);
 588                                 *_runLogFile << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
 589                                 break;
 590                         }
 591                         // Solve the next time step
 592                         ok=solveTimeStep();
 593
 594                         if (!ok)   // The resolution failed, try with a new time interval.
 595                         {
 596                                 /* if(_dt>_precision){
 597                                         cout<<"ComputeTimeStep returned _dt="<<_dt<<endl;
 598                                         cout << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
 599                                         *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
 600                                         double dt=_dt/2;//We chose to divide the time step by 2
 601                                         abortTimeStep();//Cancel the initTimeStep
 602                                         _dt=dt;//new value of time step is previous time step divided by 2 (we do not call computeTimeStep
 603                                         //_cfl*=0.5;//If we change the cfl, we must compute the new time step with computeTimeStep
 604                                         //_dt=computeTimeStep(stop);
 605                                 }
 606                                 else{*/
 607                                         PetscPrintf(PETSC_COMM_WORLD,"Failed solving time step %d, time = %.2f, dt= %.2f, cfl = %.2f, stopping calculation \n",_nbTimeStep,_time,_dt,_cfl);
 608                                         *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
 609                                         stop=true; // Impossible to solve the next time step, the Problem has given up
 610                                         break;
 611                                 //}
 612                         }
 613                         else // The resolution was successful, validate and go to the next time step.
 614                         {
 615                                 validateTimeStep();
 616                                 if ((_nbTimeStep-1)%_freqSave ==0){
 617                                         PetscPrintf(PETSC_COMM_WORLD,"Time step = %d, dt = %.2f, time = %.2f, ||Un+1-Un||= %.2f\n\n",_nbTimeStep,_dt,_time,_erreur_rel);
 618                                         *_runLogFile << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl<<endl;
 619                                 }
 620                         }
 621                 }
 622         }
 623         if(_isStationary){
 624                 PetscPrintf(PETSC_COMM_WORLD,"Stationary state reached\n");
 625                 *_runLogFile << "Stationary state reached" <<endl;
 626         }
 627         else if(_time>=_timeMax){
 628                 PetscPrintf(PETSC_COMM_WORLD,"Maximum time %.2f reached\n",_timeMax);
 629                 *_runLogFile<<"Maximum time "<<_timeMax<<" reached"<<endl;
 630         }
 631         else if(_nbTimeStep>=_maxNbOfTimeStep){
 632                 PetscPrintf(PETSC_COMM_WORLD,"Maximum number of time steps %d reached\n",_maxNbOfTimeStep);
 633                 *_runLogFile<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
 634         }
 635         else{
 636                 PetscPrintf(PETSC_COMM_WORLD,"Error problem wants to stop!\n");
 637                 *_runLogFile<<"Error problem wants to stop!"<<endl;
 638         }
 639         PetscPrintf(PETSC_COMM_WORLD,"End of calculation at time t = %.2f and time step number %d\n",_time,_nbTimeStep);
 640         *_runLogFile << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
 641
 642         _runLogFile->close();
 643         return !stop;
 644 }
 645
 646 void ProblemCoreFlows::displayMatrix(double *matrix, int size, string name)
 647 {
 648         cout<<name<<endl;
 649         for(int p=0; p<size; p++)
 650         {
 651                 for(int q=0; q<size; q++)
 652                 {
 653                         cout << matrix[p*size+q] << "\t";
 654                 }
 655                 cout << endl;
 656         }
 657         cout << endl;
 658 }
 659
 660 void ProblemCoreFlows::displayVector(double *vector, int size, string name)
 661 {
 662         cout<<name<<endl;
 663         for(int q=0; q<size; q++)
 664         {
 665                 cout << vector[q] << "\t";
 666         }
 667         cout << endl;
 668 }
 669 void ProblemCoreFlows::setFileName(string fileName){
 670         if( _nbTimeStep>0 && fileName!=_fileName)//This is a change of file name during a simulation
 671                 _restartWithNewFileName=true;
 672         _fileName = fileName;
 673 }
 674
 675 void ProblemCoreFlows::setFreqSave(int freqSave){
 676         _freqSave = freqSave;
 677 }
 678
 679 bool ProblemCoreFlows::solveTimeStep(){
 680         _NEWTON_its=0;
 681         bool converged=false, ok=true;
 682         while(!converged && ok && _NEWTON_its < _maxNewtonIts){
 683                 ok=iterateTimeStep(converged);//resolution du systeme lineaire si schema implicite
 684
 685                 if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)//To monitor the convergence of the newton scheme
 686                 {
 687                         PetscPrintf(PETSC_COMM_WORLD," Newton iteration %d, %s iterations : %d maximum variation ||Uk+1-Uk||: %.2f\n",_NEWTON_its,_ksptype,_PetscIts,_erreur_rel);
 688                         *_runLogFile<< " Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
 689
 690                         if(_conditionNumber)
 691                         {
 692                                 PetscReal sv_max, sv_min;
 693                                 KSPComputeExtremeSingularValues(_ksp, &sv_max, &sv_min);
 694                                 PetscPrintf(PETSC_COMM_WORLD," Singular value max = %.2f, singular value min = %.2f, condition number = %.2f\n",sv_max,sv_min,sv_max/sv_min);
 695                                 *_runLogFile<<" Singular value max = " << sv_max <<", singular value min = " << sv_min <<", condition number = " << sv_max/sv_min <<endl;
 696                         }
 697                 }
 698                 _NEWTON_its++;
 699         }
 700         if(!converged){
 701                 if(_NEWTON_its >= _maxNewtonIts){
 702                         PetscPrintf(PETSC_COMM_WORLD,"Maximum number of Newton iterations %d reached\n",_maxNewtonIts);
 703                         *_runLogFile << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
 704                 }
 705                 else if(!ok){
 706                         PetscPrintf(PETSC_COMM_WORLD,"iterateTimeStep: solving Newton iteration %d Failed\n",_NEWTON_its);
 707                         *_runLogFile<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
 708                 }
 709         }
 710         else if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)
 711         {
 712                 PetscPrintf(PETSC_COMM_WORLD,"Nombre d'iterations de Newton %d, Nombre max d'iterations %s : %d\n\n",_NEWTON_its, _ksptype, _MaxIterLinearSolver);
 713                 *_runLogFile <<endl;
 714                 *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its << "Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl << endl;
 715                 _MaxIterLinearSolver = 0;
 716         }
 717
 718         return converged;
 719 }
 720
 721 ProblemCoreFlows::~ProblemCoreFlows()
 722 {
 723         PetscBool petscInitialized;
 724         PetscInitialized(&petscInitialized);
 725         if(petscInitialized)
 726                 PetscFinalize();
 727
 728         delete _runLogFile;
 729 }
 730
 731 double
 732 ProblemCoreFlows::getConditionNumber(bool isSingular, double tol) const
 733 {
 734   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 735   return A.getConditionNumber( isSingular, tol);
 736 }
 737 std::vector< double >
 738 ProblemCoreFlows::getEigenvalues(int nev, EPSWhich which, double tol) const
 739 {
 740   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 741   return A.getEigenvalues( nev, which, tol);
 742 }
 743 std::vector< Vector >
 744 ProblemCoreFlows::getEigenvectors(int nev, EPSWhich which, double tol) const
 745 {
 746   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 747   return A.getEigenvectors( nev, which, tol);
 748 }
 749 Field
 750 ProblemCoreFlows::getEigenvectorsField(int nev, EPSWhich which, double tol) const
 751 {
 752   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 753   MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
 754   Field my_eigenfield;
 755
 756   my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
 757
 758   my_eigenfield.setFieldByDataArrayDouble(d);
 759
 760   return my_eigenfield;
 761 }
 762
 763 std::vector< double >
 764 ProblemCoreFlows::getSingularValues( int nsv, SVDWhich which, double tol) const
 765 {
 766   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 767   return A.getSingularValues( nsv, which, tol);
 768 }
 769 std::vector< Vector >
 770 ProblemCoreFlows::getSingularVectors(int nsv, SVDWhich which, double tol) const
 771 {
 772   SparseMatrixPetsc A = SparseMatrixPetsc(_A);
 773   return A.getSingularVectors( nsv, which, tol);
 774 }
 775
 776 Field
 777 ProblemCoreFlows::getUnknownField() const
 778 {
 779         if(!_initializedMemory)
 780         {
 781                 _runLogFile->close();
 782                 throw CdmathException("ProblemCoreFlows::getUnknownField() call initialize() first");
 783         }
 784         return _VV;
 785 }
 786
 787 bool
 788 ProblemCoreFlows::isMEDCoupling64Bits() const
 789 {
 790 #ifdef MEDCOUPLING_USE_64BIT_IDS
 791         return  true;
 792 #else
 793         return false;
 794 #endif
 795 };