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