CoreFlows/Models/src/ProblemFluid.cxx

   1 #include "ProblemFluid.hxx"
   2 #include "math.h"
   3
   4 using namespace std;
   5
   6 ProblemFluid::ProblemFluid(MPI_Comm comm):ProblemCoreFlows(comm)
   7 {
   8         _latentHeat=1e30;
   9         _Tsat=1e30;
  10         _Psat=-1e30;
  11         _dHsatl_over_dp=0;
  12         _porosityFieldSet=false;
  13         _pressureLossFieldSet=false;
  14         _sectionFieldSet=false;
  15         _GravityField3d=vector<double>(3,0);
  16         _gravityReferencePoint=vector<double>(3,0);
  17         _Uroe=NULL;_Udiff=NULL;_temp=NULL;_l=NULL;_vec_normal=NULL;;
  18         _idm=NULL;_idn=NULL;
  19         _saveVelocity=false;
  20         _saveConservativeField=false;
  21         _usePrimitiveVarsInNewton=false;
  22         _saveInterfacialField=false;
  23         _err_press_max=0; _part_imag_max=0; _nbMaillesNeg=0; _nbVpCplx=0;_minm1=1e30;_minm2=1e30;//Pour affichage paramètres diphasiques dans IterateTimeStep
  24         _isScaling=false;
  25         _entropicCorrection=false;
  26         _pressureCorrectionOrder=2;
  27         _nonLinearFormulation=Roe;
  28         _maxvploc=0.;
  29         _spaceScheme=upwind;
  30 }
  31
  32 SpaceScheme ProblemFluid::getSpaceScheme() const
  33 {
  34         return _spaceScheme;
  35 }
  36 void ProblemFluid::setNumericalScheme(SpaceScheme spaceScheme, TimeScheme timeScheme)
  37 {
  38         if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
  39                 _restartWithNewTimeScheme=true;
  40         _timeScheme = timeScheme;
  41         _spaceScheme = spaceScheme;
  42 }
  43
  44 void ProblemFluid::initialize()
  45 {
  46         if(!_initialDataSet)
  47         {
  48                 *_runLogFile<<"!!!!!!!!ProblemFluid::initialize() set initial data first"<<endl;
  49                 _runLogFile->close();
  50                 throw CdmathException("!!!!!!!!ProblemFluid::initialize() set initial data first");
  51         }
  52         else if (_VV.getTypeOfField() != CELLS)
  53         {
  54                 *_runLogFile<<"!!!!!!!!Initial data should be a field on CELLS, not NODES, neither FACES"<<endl;
  55                 _runLogFile->close();
  56                 throw CdmathException("!!!!!!!!ProblemFluid::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
  57         }
  58         cout << "Number of Phases = " << _nbPhases << " mesh dimension = "<<_Ndim<<" number of variables = "<<_nVar<<endl;
  59         *_runLogFile << "Number of Phases = " << _nbPhases << " spaceDim= "<<_Ndim<<" number of variables= "<<_nVar<<endl;
  60
  61         /********* local arrays ****************/
  62         _AroePlus = new PetscScalar[_nVar*_nVar];
  63         _Diffusion = new PetscScalar[_nVar*_nVar];
  64         _AroeMinus = new PetscScalar[_nVar*_nVar];
  65         _Aroe = new PetscScalar[_nVar*_nVar];
  66         _absAroe = new PetscScalar[_nVar*_nVar];
  67         _signAroe = new PetscScalar[_nVar*_nVar];
  68         _invAroe = new PetscScalar[_nVar*_nVar];
  69         _AroeImplicit = new PetscScalar[_nVar*_nVar];
  70         _AroeMinusImplicit = new PetscScalar[_nVar*_nVar];
  71         _AroePlusImplicit = new PetscScalar[_nVar*_nVar];
  72         _absAroeImplicit = new PetscScalar[_nVar*_nVar];
  73         _primToConsJacoMat = new PetscScalar[_nVar*_nVar];
  74         _phi = new PetscScalar[_nVar];
  75         _Jcb = new PetscScalar[_nVar*_nVar];
  76         _JcbDiff = new PetscScalar[_nVar*_nVar];
  77         _a = new PetscScalar[_nVar*_nVar];
  78         _externalStates = new PetscScalar[_nVar];
  79         _Ui = new PetscScalar[_nVar];
  80         _Uj = new PetscScalar[_nVar];
  81         _Vi = new PetscScalar[_nVar];
  82         _Vj = new PetscScalar[_nVar];
  83         _Uij = new PetscScalar[_nVar];//used for VFRoe scheme
  84         _Vij = new PetscScalar[_nVar];//used for VFRoe scheme
  85         _Si = new PetscScalar[_nVar];
  86         _Sj = new PetscScalar[_nVar];
  87         _pressureLossVector= new PetscScalar[_nVar];
  88         _porosityGradientSourceVector= new PetscScalar[_nVar];
  89
  90         _l = new double[_nVar];
  91         _r = new double[_nVar];
  92         _Udiff = new double[_nVar];
  93         _temp = new double[_nVar*_nVar];
  94
  95         _idm = new int[_nVar];
  96         _idn = new int[_nVar];
  97         _vec_normal = new double[_Ndim];
  98
  99         for(int k=0;k<_nVar*_nVar;k++)
 100         {
 101                 _Diffusion[k]=0;
 102                 _JcbDiff[k]=0;
 103         }
 104         for(int k=0; k<_nVar; k++){
 105                 _idm[k] = k;
 106                 _idn[k] = k;
 107         }
 108
 109         /**************** Field creation *********************/
 110         if(!_heatPowerFieldSet){
 111                 _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
 112                 for(int i =0; i<_Nmailles; i++)
 113                         _heatPowerField(i) = _heatSource;
 114         }
 115         if(_Psat>-1e30)
 116                 _dp_over_dt=Field("dP/dt",CELLS,_mesh,1);
 117         if(!_porosityFieldSet){
 118                 _porosityField=Field("Porosity field",CELLS,_mesh,1);
 119                 for(int i =0; i<_Nmailles; i++)
 120                         _porosityField(i) = 1;
 121                 _porosityFieldSet=true;
 122         }
 123
 124         //conservative field used only for saving results
 125         _UU=Field ("Conservative vec", CELLS, _mesh, _nVar);
 126         if(_saveInterfacialField && _nonLinearFormulation==VFRoe)
 127         {
 128                 _UUstar=Field ("Interfacial U", CELLS, _mesh, _nVar);
 129                 _VVstar=Field ("Interfacial V", CELLS, _mesh, _nVar);
 130         }
 131
 132         //Construction des champs primitifs et conservatifs initiaux comme avant dans ParaFlow
 133         double * initialFieldPrim = new double[_nVar*_Nmailles];
 134         for(int i =0; i<_Nmailles;i++)
 135                 for(int j =0; j<_nVar;j++)
 136                         initialFieldPrim[i*_nVar+j]=_VV(i,j);
 137         double *initialFieldCons = new double[_nVar*_Nmailles];
 138         for(int i=0; i<_Nmailles; i++)
 139                 primToCons(initialFieldPrim, i, initialFieldCons, i);
 140         for(int i =0; i<_Nmailles;i++)
 141                 for(int j =0; j<_nVar;j++)
 142                         _UU(i,j)=initialFieldCons[i*_nVar+j];
 143
 144         /**********Petsc structures:  ****************/
 145
 146         //creation de la matrice
 147         if(_timeScheme == Implicit)
 148                 MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
 149
 150         //creation des vecteurs
 151         VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uext);
 152         VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Vext);
 153         VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uextdiff);
 154         //        VecCreateSeq(PETSC_COMM_SELF, _nVar*_Nmailles, &_conservativeVars);
 155         VecCreate(PETSC_COMM_SELF, &_conservativeVars);//Current conservative variables at Newton iteration k between time steps n and n+1
 156         VecSetSizes(_conservativeVars,PETSC_DECIDE,_nVar*_Nmailles);
 157         VecSetBlockSize(_conservativeVars,_nVar);
 158         VecSetFromOptions(_conservativeVars);
 159         VecDuplicate(_conservativeVars, &_old);//Old conservative variables at time step n
 160         VecDuplicate(_conservativeVars, &_newtonVariation);//Newton variation Uk+1-Uk to be computed between time steps n and n+1
 161         VecDuplicate(_conservativeVars, &_b);//Right hand side of Newton method at iteration k between time steps n and n+1
 162         VecDuplicate(_conservativeVars, &_primitiveVars);//Current primitive variables at Newton iteration k between time steps n and n+1
 163
 164         if(_isScaling)
 165         {
 166                 VecDuplicate(_conservativeVars, &_vecScaling);
 167                 VecDuplicate(_conservativeVars, &_invVecScaling);
 168                 VecDuplicate(_conservativeVars, &_bScaling);
 169
 170                 _blockDiag = new PetscScalar[_nVar];
 171                 _invBlockDiag = new PetscScalar[_nVar];
 172         }
 173
 174
 175         int *indices = new int[_Nmailles];
 176         for(int i=0; i<_Nmailles; i++)
 177                 indices[i] = i;
 178         VecSetValuesBlocked(_conservativeVars, _Nmailles, indices, initialFieldCons, INSERT_VALUES);
 179         VecAssemblyBegin(_conservativeVars);
 180         VecAssemblyEnd(_conservativeVars);
 181         VecCopy(_conservativeVars, _old);
 182         VecAssemblyBegin(_old);
 183         VecAssemblyEnd(_old);
 184         VecSetValuesBlocked(_primitiveVars, _Nmailles, indices, initialFieldPrim, INSERT_VALUES);
 185         VecAssemblyBegin(_primitiveVars);
 186         VecAssemblyEnd(_primitiveVars);
 187         if(_system)
 188         {
 189                 cout << "Variables primitives initiales : " << endl;
 190                 VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_WORLD);
 191                 cout << endl;
 192                 cout<<"Variables conservatives initiales : "<<endl;
 193                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_SELF);
 194         }
 195
 196         delete[] initialFieldPrim;
 197         delete[] initialFieldCons;
 198         delete[] indices;
 199
 200         //Linear solver
 201         KSPCreate(PETSC_COMM_SELF, &_ksp);
 202         KSPSetType(_ksp, _ksptype);
 203         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 204         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 205         KSPGetPC(_ksp, &_pc);
 206         PCSetType(_pc, _pctype);
 207
 208         // Creation du solveur de Newton de PETSc
 209         if( _timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
 210         {
 211                 SNESType snestype;
 212
 213                 // set nonlinear solver
 214                 if (_nonLinearSolver == Newton_PETSC_LINESEARCH || _nonLinearSolver == Newton_PETSC_LINESEARCH_BASIC || _nonLinearSolver == Newton_PETSC_LINESEARCH_BT || _nonLinearSolver == Newton_PETSC_LINESEARCH_SECANT || _nonLinearSolver == Newton_PETSC_LINESEARCH_NLEQERR)
 215                         snestype = (char*)&SNESNEWTONLS;
 216                 else if (_nonLinearSolver == Newton_PETSC_TRUSTREGION)
 217                         snestype = (char*)&SNESNEWTONTR;
 218                 else if (_nonLinearSolver == Newton_PETSC_NGMRES)
 219                         snestype = (char*)&SNESNGMRES;
 220                 else if (_nonLinearSolver ==Newton_PETSC_ASPIN)
 221                         snestype = (char*)&SNESASPIN;
 222                 else if(_nonLinearSolver != Newton_SOLVERLAB)
 223                 {
 224                         cout << "!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" << endl;
 225                         *_runLogFile << "!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" << endl;
 226                         _runLogFile->close();
 227                         throw CdmathException("!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" );
 228                 }
 229
 230                 SNESCreate(PETSC_COMM_WORLD, &_snes);
 231                 SNESSetType( _snes, snestype);
 232                 SNESGetLineSearch( _snes, &_linesearch);
 233                 if(_nonLinearSolver == Newton_PETSC_LINESEARCH_BASIC)
 234                         SNESLineSearchSetType( _linesearch,     SNESLINESEARCHBASIC );
 235                 else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_BT)
 236                         SNESLineSearchSetType( _linesearch,     SNESLINESEARCHBT );
 237                 else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_SECANT)
 238                         SNESLineSearchSetType( _linesearch,     SNESLINESEARCHL2 );
 239                 else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_NLEQERR)
 240                         SNESLineSearchSetType( _linesearch,     SNESLINESEARCHNLEQERR );
 241
 242                 PetscViewerCreate(PETSC_COMM_WORLD,&_monitorLineSearch);
 243                 PetscViewerSetType(_monitorLineSearch, PETSCVIEWERASCII);
 244
 245                 SNESSetTolerances(_snes,_precision_Newton,_precision_Newton,_precision_Newton,_maxNewtonIts,-1);
 246
 247                 SNESSetFunction(_snes,_newtonVariation,computeSnesRHS,this);
 248                 SNESSetJacobian(_snes,_A,_A,computeSnesJacobian,this);
 249         }
 250
 251         _initializedMemory=true;
 252         save();//save initial data
 253 }
 254
 255 bool ProblemFluid::initTimeStep(double dt){
 256         _dt = dt;
 257         return _dt>0;//No need to call MatShift as the linear system matrix is filled at each Newton iteration (unlike linear problem)
 258 }
 259
 260 bool ProblemFluid::solveTimeStep(){
 261         if(_timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
 262                 return solveNewtonPETSc();
 263         else
 264                 return ProblemCoreFlows::solveTimeStep();
 265 }
 266
 267 bool ProblemFluid::solveNewtonPETSc()
 268 {
 269         if( (_nbTimeStep-1)%_freqSave ==0)
 270                 SNESLineSearchSetDefaultMonitor(_linesearch, _monitorLineSearch);
 271         else
 272                 SNESLineSearchSetDefaultMonitor(_linesearch, NULL);
 273
 274     SNESSolve(_snes,NULL,_conservativeVars);
 275
 276     SNESConvergedReason reason;
 277         SNESGetConvergedReason(_snes,&reason);
 278
 279         if( (_nbTimeStep-1)%_freqSave ==0)
 280         {
 281                 if(reason == SNES_CONVERGED_FNORM_ABS  )
 282                         cout<<"Converged with absolute norm (absolute tolerance) less than "<<_precision_Newton<<", (||F|| < atol)"<<endl;
 283                 else if(reason == SNES_CONVERGED_FNORM_RELATIVE  )
 284                         cout<<"Converged because residual has decreased by a factor less than "<<_precision_Newton<<", (||F|| < rtol*||F_initial||)"<<endl;
 285                 else if(reason == SNES_CONVERGED_SNORM_RELATIVE  )
 286                         cout<<"Converged with  relative norm (relative tolerance) less than "<<_precision_Newton<<", (|| delta x || < stol || x ||)"<<endl;
 287                 else if(reason == SNES_CONVERGED_ITS  )
 288                         cout<<"SNESConvergedSkip() was chosen as the convergence test; thus the usual convergence criteria have not been checked and may or may not be satisfied"<<endl;
 289                 else if(reason == SNES_DIVERGED_LINEAR_SOLVE  )
 290                         cout<<"Solving linear system failed"<<endl;
 291                 else if(reason ==  SNES_DIVERGED_LINE_SEARCH   )
 292                         cout<<"Line search failed"<<endl;
 293                 else if(reason ==   SNES_DIVERGED_MAX_IT    )
 294                         cout<<"Reached the maximum number of iterations"<<endl;
 295
 296             SNESGetIterationNumber(_snes, &_NEWTON_its);
 297             PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n\n", _NEWTON_its);
 298                 *_runLogFile <<endl;
 299                 *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its <<endl;
 300         }
 301
 302         return reason>0;
 303 }
 304
 305 bool ProblemFluid::iterateTimeStep(bool &converged)
 306 {
 307         bool stop=false;
 308
 309         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 310                 _maxvp=0.;
 311                 computeTimeStep(stop);//This compute timestep is just to update the linear system. The time step was imposed before starting the Newton iterations
 312         }
 313         if(stop){//Le compute time step ne s'est pas bien passé
 314                 cout<<"ComputeTimeStep failed"<<endl;
 315                 converged=false;
 316                 return false;
 317         }
 318
 319         computeNewtonVariation();
 320
 321         //converged=convergence des iterations
 322         if(_timeScheme == Explicit)
 323                 converged=true;
 324         else{//Implicit scheme
 325
 326                 KSPGetIterationNumber(_ksp, &_PetscIts);
 327                 if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
 328                         _MaxIterLinearSolver = _PetscIts;
 329
 330                 KSPConvergedReason reason;
 331                 KSPGetConvergedReason(_ksp,&reason);
 332
 333                 if(reason<0)//solving the linear system failed
 334                 {
 335                         cout<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
 336                         cout<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
 337                         *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
 338                         *_runLogFile<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
 339                         converged=false;
 340                         return false;
 341                 }
 342                 else{//solving the linear system succeeded
 343                         //Calcul de la variation relative Uk+1-Uk
 344                         _erreur_rel = 0.;
 345                         double x, dx;
 346                         int I;
 347                         for(int j=1; j<=_Nmailles; j++)
 348                         {
 349                                 for(int k=0; k<_nVar; k++)
 350                                 {
 351                                         I = (j-1)*_nVar + k;
 352                                         VecGetValues(_newtonVariation, 1, &I, &dx);
 353                                         VecGetValues(_conservativeVars, 1, &I, &x);
 354                                         if (fabs(x)*fabs(x)< _precision)
 355                                         {
 356                                                 if(_erreur_rel < fabs(dx))
 357                                                         _erreur_rel = fabs(dx);
 358                                         }
 359                                         else if(_erreur_rel < fabs(dx/x))
 360                                                 _erreur_rel = fabs(dx/x);
 361                                 }
 362                         }
 363                 }
 364                 converged = _erreur_rel <= _precision_Newton;
 365         }
 366
 367         double relaxation=1;//Uk+1=Uk+relaxation*deltaU
 368
 369         VecAXPY(_conservativeVars,  relaxation, _newtonVariation);
 370
 371         //mise a jour du champ primitif
 372         updatePrimitives();
 373
 374         if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
 375         {
 376                 cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 377                 *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 378                 converged=false;
 379                 return false;
 380         }
 381         if(_system)
 382         {
 383                 cout<<"Vecteur Ukp1-Uk "<<endl;
 384                 VecView(_newtonVariation,  PETSC_VIEWER_STDOUT_SELF);
 385                 cout << "Nouvel etat courant Uk de l'iteration Newton: " << endl;
 386                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_SELF);
 387         }
 388
 389         if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
 390                 if(_minm1<-_precision || _minm2<-_precision)
 391                 {
 392                         cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 393                         *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 394                 }
 395
 396                 if (_nbVpCplx>0){
 397                         cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 398                         *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 399                 }
 400         }
 401
 402         return true;
 403 }
 404
 405 double ProblemFluid::computeTimeStep(bool & stop){//dt is not known and will not contribute to the Newton scheme
 406
 407         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 408         {
 409                 cout << "ProblemFluid::computeTimeStep : Début calcul matrice implicite et second membre"<<endl;
 410                 cout << endl;
 411         }
 412         if(_restartWithNewTimeScheme)//This is a change of time scheme during a simulation
 413         {
 414                 if(_timeScheme == Implicit)
 415                         MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
 416                 else
 417                         MatDestroy(&_A);
 418                 _restartWithNewTimeScheme=false;
 419         }
 420         if(_timeScheme == Implicit)
 421                 MatZeroEntries(_A);
 422
 423         VecAssemblyBegin(_b);
 424         VecZeroEntries(_b);
 425
 426         std::vector< int > idCells(2);
 427         PetscInt idm, idn, size = 1;
 428
 429         long nbFaces = _mesh.getNumberOfFaces();
 430         Face Fj;
 431         Cell Ctemp1,Ctemp2;
 432         string nameOfGroup;
 433
 434         for (int j=0; j<nbFaces;j++){
 435                 Fj = _mesh.getFace(j);
 436                 _isBoundary=Fj.isBorder();
 437                 idCells = Fj.getCellsId();
 438
 439                 // If Fj is on the boundary
 440                 if (_isBoundary)
 441                 {
 442                         for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
 443                         {
 444                                 // compute the normal vector corresponding to face j : from Ctemp1 outward
 445                                 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
 446                                 if (_Ndim >1){
 447                                         for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
 448                                         {//we look for l the index of the face Fj for the cell Ctemp1
 449                                                 if (j == Ctemp1.getFacesId()[l])
 450                                                 {
 451                                                         for (int idim = 0; idim < _Ndim; ++idim)
 452                                                                 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
 453                                                         break;
 454                                                 }
 455                                         }
 456                                 }else{ // _Ndim = 1, build normal vector (bug cdmath)
 457                                         if(!_sectionFieldSet)
 458                                         {
 459                                                 if (Fj.x()<Ctemp1.x())
 460                                                         _vec_normal[0] = -1;
 461                                                 else
 462                                                         _vec_normal[0] = 1;
 463                                         }
 464                                         else
 465                                         {
 466                                                 if(idCells[0]==0)
 467                                                         _vec_normal[0] = -1;
 468                                                 else//idCells[0]==31
 469                                                         _vec_normal[0] = 1;
 470                                         }
 471                                 }
 472                                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 473                                 {
 474                                         cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
 475                                         for(int p=0; p<_Ndim; p++)
 476                                                 cout << _vec_normal[p] << ",";
 477                                         cout << "). "<<endl;
 478                                 }
 479                                 nameOfGroup = Fj.getGroupName();
 480                                 _porosityi=_porosityField(idCells[k]);
 481                                 _porosityj=_porosityi;
 482                                 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
 483                                 convectionState(idCells[k],0,true);
 484                                 convectionMatrices();
 485                                 diffusionStateAndMatrices(idCells[k], 0, true);
 486                                 // compute 1/dxi
 487                                 if (_Ndim > 1)
 488                                         _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
 489                                 else
 490                                         _inv_dxi = 1/Ctemp1.getMeasure();
 491
 492                                 addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
 493                                 addDiffusionToSecondMember(idCells[k],-1,true);
 494                                 addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
 495
 496                                 if(_timeScheme == Implicit){
 497                                         for(int l=0; l<_nVar*_nVar;l++){
 498                                                 _AroeMinusImplicit[l] *= _inv_dxi;
 499                                                 _Diffusion[l] *=_inv_dxi*_inv_dxi;
 500                                         }
 501
 502                                         jacobian(idCells[k],nameOfGroup,_vec_normal);
 503                                         jacobianDiff(idCells[k],nameOfGroup);
 504                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
 505                                                 cout << "Matrice Jacobienne CL convection:" << endl;
 506                                                 for(int p=0; p<_nVar; p++){
 507                                                         for(int q=0; q<_nVar; q++)
 508                                                                 cout << _Jcb[p*_nVar+q] << "\t";
 509                                                         cout << endl;
 510                                                 }
 511                                                 cout << endl;
 512                                                 cout << "Matrice Jacobienne CL diffusion:" << endl;
 513                                                 for(int p=0; p<_nVar; p++){
 514                                                         for(int q=0; q<_nVar; q++)
 515                                                                 cout << _JcbDiff[p*_nVar+q] << "\t";
 516                                                         cout << endl;
 517                                                 }
 518                                                 cout << endl;
 519                                         }
 520                                         idm = idCells[k];
 521                                         //calcul et insertion de A^-*Jcb
 522                                         Polynoms::matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
 523                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
 524
 525                                         if(_verbose)
 526                                                 displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
 527
 528                                         //insertion de -A^-
 529                                         for(int k=0; k<_nVar*_nVar;k++){
 530                                                 _AroeMinusImplicit[k] *= -1;
 531                                         }
 532                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 533                                         if(_verbose)
 534                                                 displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
 535
 536                                         //calcul et insertion de D*JcbDiff
 537                                         Polynoms::matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
 538                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
 539                                         for(int k=0; k<_nVar*_nVar;k++)
 540                                                 _Diffusion[k] *= -1;
 541                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 542                                 }
 543                         }
 544                 } else  if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
 545                         // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
 546                         Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
 547                         Ctemp2 = _mesh.getCell(idCells[1]);
 548                         if (_Ndim >1){
 549                                 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
 550                                         if (j == Ctemp1.getFacesId()[l]){
 551                                                 for (int idim = 0; idim < _Ndim; ++idim)
 552                                                         _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
 553                                                 break;
 554                                         }
 555                                 }
 556                         }else{ // _Ndim = 1, build normal vector (bug cdmath)
 557                                 if(!_sectionFieldSet)
 558                                 {
 559                                         if (Fj.x()<Ctemp1.x())
 560                                                 _vec_normal[0] = -1;
 561                                         else
 562                                                 _vec_normal[0] = 1;
 563                                 }
 564                                 else
 565                                 {
 566                                         if(idCells[0]>idCells[1])
 567                                                 _vec_normal[0] = -1;
 568                                         else
 569                                                 _vec_normal[0] = 1;
 570                                 }
 571                         }
 572                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 573                         {
 574                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 575                                 cout<<" Normal vector= ";
 576                                 for (int idim = 0; idim < _Ndim; ++idim)
 577                                         cout<<_vec_normal[idim]<<", ";
 578                                 cout<<endl;
 579                         }
 580                         _porosityi=_porosityField(idCells[0]);
 581                         _porosityj=_porosityField(idCells[1]);
 582                         convectionState(idCells[0],idCells[1],false);
 583                         convectionMatrices();
 584                         diffusionStateAndMatrices(idCells[0], idCells[1], false);
 585
 586                         // compute 1/dxi and 1/dxj
 587                         if (_Ndim > 1)
 588                         {
 589                                 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
 590                                 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
 591                         }
 592                         else
 593                         {
 594                                 _inv_dxi = 1/Ctemp1.getMeasure();
 595                                 _inv_dxj = 1/Ctemp2.getMeasure();
 596                         }
 597
 598                         addConvectionToSecondMember(idCells[0],idCells[1], false);
 599                         addDiffusionToSecondMember( idCells[0],idCells[1], false);
 600                         addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
 601
 602                         if(_timeScheme == Implicit){
 603                                 for(int k=0; k<_nVar*_nVar;k++)
 604                                 {
 605                                         _AroeMinusImplicit[k] *= _inv_dxi;
 606                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
 607                                 }
 608                                 idm = idCells[0];
 609                                 idn = idCells[1];
 610                                 //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNbCells<<endl;
 611                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
 612                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
 613
 614                                 if(_verbose){
 615                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
 616                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
 617                                 }
 618                                 for(int k=0;k<_nVar*_nVar;k++){
 619                                         _AroeMinusImplicit[k] *= -1;
 620                                         _Diffusion[k] *= -1;
 621                                 }
 622                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 623                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 624                                 if(_verbose){
 625                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
 626                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
 627                                 }
 628                                 for(int k=0; k<_nVar*_nVar;k++)
 629                                 {
 630                                         _AroePlusImplicit[k]  *= _inv_dxj;
 631                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;
 632                                 }
 633                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
 634                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
 635                                 if(_verbose)
 636                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
 637
 638                                 for(int k=0;k<_nVar*_nVar;k++){
 639                                         _AroePlusImplicit[k] *= -1;
 640                                         _Diffusion[k] *= -1;
 641                                 }
 642                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
 643                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
 644
 645                                 if(_verbose)
 646                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
 647                         }
 648                 }
 649                 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
 650                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 651                                 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
 652                         *_runLogFile<<"Warning: treatment of a junction node"<<endl;
 653
 654                         if(!_sectionFieldSet)
 655                         {
 656                                 _runLogFile->close();
 657                                 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
 658                         }
 659                         int largestSectionCellIndex=0;
 660                         for(int i=1;i<Fj.getNumberOfCells();i++){
 661                                 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
 662                                         largestSectionCellIndex=i;
 663                         }
 664                         idm = idCells[largestSectionCellIndex];
 665                         Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
 666                         _porosityi=_porosityField(idm);
 667
 668                         if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
 669                                 _vec_normal[0] = 1;
 670                         else//j==16
 671                                 _vec_normal[0] = -1;
 672                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 673                         {
 674                                 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
 675                                 cout << " ; vecteur normal=(";
 676                                 for(int p=0; p<_Ndim; p++)
 677                                         cout << _vec_normal[p] << ",";
 678                                 cout << "). "<<endl;
 679                         }
 680                         for(int i=0;i<Fj.getNumberOfCells();i++){
 681                                 if(i != largestSectionCellIndex){
 682                                         idn = idCells[i];
 683                                         Ctemp2 = _mesh.getCell(idn);
 684                                         _porosityj=_porosityField(idn);
 685                                         convectionState(idm,idn,false);
 686                                         convectionMatrices();
 687                                         diffusionStateAndMatrices(idm, idn,false);
 688
 689                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 690                                                 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
 691
 692                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
 693                                         _inv_dxj = 1/Ctemp2.getMeasure();
 694
 695                                         addConvectionToSecondMember(idm,idn, false);
 696                                         _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
 697                                         addDiffusionToSecondMember(idm,idn, false);
 698                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
 699                                         addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
 700
 701                                         if(_timeScheme == Implicit){
 702                                                 for(int k=0; k<_nVar*_nVar;k++)
 703                                                 {
 704                                                         _AroeMinusImplicit[k] *= _inv_dxi;
 705                                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
 706                                                 }
 707                                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
 708                                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
 709
 710                                                 if(_verbose){
 711                                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
 712                                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
 713                                                 }
 714                                                 for(int k=0;k<_nVar*_nVar;k++){
 715                                                         _AroeMinusImplicit[k] *= -1;
 716                                                         _Diffusion[k] *= -1;
 717                                                 }
 718                                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 719                                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 720                                                 if(_verbose){
 721                                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
 722                                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
 723                                                 }
 724                                                 for(int k=0; k<_nVar*_nVar;k++)
 725                                                 {
 726                                                         _AroePlusImplicit[k] *= _inv_dxj;
 727                                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
 728                                                 }
 729                                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
 730                                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
 731                                                 if(_verbose)
 732                                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
 733
 734                                                 for(int k=0;k<_nVar*_nVar;k++){
 735                                                         _AroePlusImplicit[k] *= -1;
 736                                                         _Diffusion[k] *= -1;
 737                                                 }
 738                                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
 739                                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
 740
 741                                                 if(_verbose)
 742                                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
 743                                         }
 744                                 }
 745                         }
 746                 }
 747                 else
 748                 {
 749                         cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
 750                         _runLogFile->close();
 751                         throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
 752                 }
 753
 754         }
 755         VecAssemblyEnd(_b);
 756
 757         if(_timeScheme == Implicit){
 758                 for(int imaille = 0; imaille<_Nmailles; imaille++)
 759                         MatSetValuesBlocked(_A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
 760
 761                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 762                         displayMatrix(_GravityImplicitationMatrix,_nVar,"Gravity matrix:");
 763
 764                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 765                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 766                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
 767                         cout << "ProblemFluid::computeTimeStep : Fin calcul matrice implicite et second membre"<<endl;
 768                         cout << "ProblemFluid::computeTimeStep : Matrice implicite :"<<endl;
 769                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 770                         cout << "ProblemFluid::computeTimeStep : Second membre :"<<endl;
 771                         VecView(_b,  PETSC_VIEWER_STDOUT_WORLD);
 772                         cout << endl;
 773                 }
 774         }
 775
 776         stop=false;
 777
 778         /*
 779         if(_nbTimeStep+1<_cfl)
 780                 return (_nbTimeStep+1)*_minl/_maxvp;
 781         else
 782          */
 783         return _cfl*_minl/_maxvp;
 784 }
 785
 786 void ProblemFluid::computeNewtonVariation()
 787 {
 788         if(_system)
 789         {
 790                 cout<<"Vecteur courant Uk "<<endl;
 791                 VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
 792                 cout << endl;
 793         }
 794         if(_timeScheme == Explicit)
 795         {
 796                 VecCopy(_b,_newtonVariation);
 797                 VecScale(_newtonVariation, _dt);
 798                 if(_system && (_nbTimeStep-1)%_freqSave ==0)
 799                 {
 800                         cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
 801                         VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
 802                         cout << endl;
 803                 }
 804         }
 805         else
 806         {
 807                 if(_system)
 808                 {
 809                         cout << "Matrice du système linéaire avant contribution delta t" << endl;
 810                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 811                         cout << endl;
 812                         cout << "Second membre du système linéaire avant contribution delta t" << endl;
 813                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 814                         cout << endl;
 815                 }
 816                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 817                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 818
 819                 VecAXPY(_b, 1/_dt, _old);
 820                 VecAXPY(_b, -1/_dt, _conservativeVars);
 821                 MatShift(_A, 1/_dt);
 822
 823 #if PETSC_VERSION_GREATER_3_5
 824                 KSPSetOperators(_ksp, _A, _A);
 825 #else
 826                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 827 #endif
 828
 829                 if(_system)
 830                 {
 831                         cout << "Matrice du système linéaire après contribution delta t" << endl;
 832                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 833                         cout << endl;
 834                         cout << "Second membre du système linéaire après contribution delta t" << endl;
 835                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 836                         cout << endl;
 837                 }
 838
 839                 if(_conditionNumber)
 840                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 841                 if(!_isScaling)
 842                 {
 843                         KSPSolve(_ksp, _b, _newtonVariation);
 844                 }
 845                 else
 846                 {
 847                         computeScaling(_maxvp);
 848                         int indice;
 849                         VecAssemblyBegin(_vecScaling);
 850                         VecAssemblyBegin(_invVecScaling);
 851                         for(int imaille = 0; imaille<_Nmailles; imaille++)
 852                         {
 853                                 indice = imaille;
 854                                 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
 855                                 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
 856                         }
 857                         VecAssemblyEnd(_vecScaling);
 858                         VecAssemblyEnd(_invVecScaling);
 859
 860                         if(_system)
 861                         {
 862                                 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
 863                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 864                                 cout << endl;
 865                                 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
 866                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 867                                 cout << endl;
 868                         }
 869                         MatDiagonalScale(_A,_vecScaling,_invVecScaling);
 870                         if(_system)
 871                         {
 872                                 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
 873                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 874                                 cout << endl;
 875                         }
 876                         VecCopy(_b,_bScaling);
 877                         VecPointwiseMult(_b,_vecScaling,_bScaling);
 878                         if(_system)
 879                         {
 880                                 cout << "Produit du second membre par le preconditionneur bloc diagonal  a gauche" << endl;
 881                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 882                                 cout << endl;
 883                         }
 884
 885                         KSPSolve(_ksp,_b, _bScaling);
 886                         VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
 887                 }
 888                 if(_system)
 889                 {
 890                         cout << "solution du systeme lineaire local:" << endl;
 891                         VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
 892                         cout << endl;
 893                 }
 894         }
 895 }
 896
 897 void ProblemFluid::computeNewtonRHS( Vec X, Vec F_X){//dt is known and will contribute to the right hand side of the Newton scheme
 898
 899         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 900         {
 901                 cout << "ProblemFluid::computeNewtonRHS : Début calcul second membre pour PETSc, _dt="<<_dt<<endl;
 902                 cout << endl;
 903                 cout<<"Vecteur courant Uk "<<endl;
 904                 VecView(X,PETSC_VIEWER_STDOUT_SELF);
 905                 cout << endl;
 906         }
 907
 908         VecCopy(X,_conservativeVars);
 909         updatePrimitives();
 910
 911         VecAssemblyBegin(_b);
 912         VecZeroEntries(_b);
 913
 914         std::vector< int > idCells(2);
 915         PetscInt idm, idn, size = 1;
 916
 917         long nbFaces = _mesh.getNumberOfFaces();
 918         Face Fj;
 919         Cell Ctemp1,Ctemp2;
 920         string nameOfGroup;
 921
 922         for (int j=0; j<nbFaces;j++){
 923                 Fj = _mesh.getFace(j);
 924                 _isBoundary=Fj.isBorder();
 925                 idCells = Fj.getCellsId();
 926
 927                 // If Fj is on the boundary
 928                 if (_isBoundary)
 929                 {
 930                         for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
 931                         {
 932                                 // compute the normal vector corresponding to face j : from Ctemp1 outward
 933                                 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
 934                                 if (_Ndim >1){
 935                                         for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
 936                                         {//we look for l the index of the face Fj for the cell Ctemp1
 937                                                 if (j == Ctemp1.getFacesId()[l])
 938                                                 {
 939                                                         for (int idim = 0; idim < _Ndim; ++idim)
 940                                                                 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
 941                                                         break;
 942                                                 }
 943                                         }
 944                                 }else{ // _Ndim = 1, build normal vector (bug cdmath)
 945                                         if(!_sectionFieldSet)
 946                                         {
 947                                                 if (Fj.x()<Ctemp1.x())
 948                                                         _vec_normal[0] = -1;
 949                                                 else
 950                                                         _vec_normal[0] = 1;
 951                                         }
 952                                         else
 953                                         {
 954                                                 if(idCells[0]==0)
 955                                                         _vec_normal[0] = -1;
 956                                                 else//idCells[0]==31
 957                                                         _vec_normal[0] = 1;
 958                                         }
 959                                 }
 960                                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 961                                 {
 962                                         cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
 963                                         for(int p=0; p<_Ndim; p++)
 964                                                 cout << _vec_normal[p] << ",";
 965                                         cout << "). "<<endl;
 966                                 }
 967                                 nameOfGroup = Fj.getGroupName();
 968                                 _porosityi=_porosityField(idCells[k]);
 969                                 _porosityj=_porosityi;
 970                                 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
 971                                 convectionState(idCells[k],0,true);
 972                                 convectionMatrices();
 973                                 diffusionStateAndMatrices(idCells[k], 0, true);
 974                                 // compute 1/dxi
 975                                 if (_Ndim > 1)
 976                                         _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
 977                                 else
 978                                         _inv_dxi = 1/Ctemp1.getMeasure();
 979
 980                                 addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
 981                                 addDiffusionToSecondMember(idCells[k],-1,true);
 982                                 addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
 983                         }
 984                 } else  if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
 985                         // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
 986                         Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
 987                         Ctemp2 = _mesh.getCell(idCells[1]);
 988                         if (_Ndim >1){
 989                                 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
 990                                         if (j == Ctemp1.getFacesId()[l]){
 991                                                 for (int idim = 0; idim < _Ndim; ++idim)
 992                                                         _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
 993                                                 break;
 994                                         }
 995                                 }
 996                         }else{ // _Ndim = 1, build normal vector (bug cdmath)
 997                                 if(!_sectionFieldSet)
 998                                 {
 999                                         if (Fj.x()<Ctemp1.x())
1000                                                 _vec_normal[0] = -1;
1001                                         else
1002                                                 _vec_normal[0] = 1;
1003                                 }
1004                                 else
1005                                 {
1006                                         if(idCells[0]>idCells[1])
1007                                                 _vec_normal[0] = -1;
1008                                         else
1009                                                 _vec_normal[0] = 1;
1010                                 }
1011                         }
1012                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1013                         {
1014                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
1015                                 cout<<" Normal vector= ";
1016                                 for (int idim = 0; idim < _Ndim; ++idim)
1017                                         cout<<_vec_normal[idim]<<", ";
1018                                 cout<<endl;
1019                         }
1020                         _porosityi=_porosityField(idCells[0]);
1021                         _porosityj=_porosityField(idCells[1]);
1022                         convectionState(idCells[0],idCells[1],false);
1023                         convectionMatrices();
1024                         diffusionStateAndMatrices(idCells[0], idCells[1], false);
1025
1026                         // compute 1/dxi and 1/dxj
1027                         if (_Ndim > 1)
1028                         {
1029                                 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
1030                                 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
1031                         }
1032                         else
1033                         {
1034                                 _inv_dxi = 1/Ctemp1.getMeasure();
1035                                 _inv_dxj = 1/Ctemp2.getMeasure();
1036                         }
1037
1038                         addConvectionToSecondMember(idCells[0],idCells[1], false);
1039                         addDiffusionToSecondMember( idCells[0],idCells[1], false);
1040                         addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
1041
1042                 }
1043                 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
1044                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1045                                 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
1046                         *_runLogFile<<"Warning: treatment of a junction node"<<endl;
1047
1048                         if(!_sectionFieldSet)
1049                         {
1050                                 _runLogFile->close();
1051                                 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
1052                         }
1053                         int largestSectionCellIndex=0;
1054                         for(int i=1;i<Fj.getNumberOfCells();i++){
1055                                 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
1056                                         largestSectionCellIndex=i;
1057                         }
1058                         idm = idCells[largestSectionCellIndex];
1059                         Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
1060                         _porosityi=_porosityField(idm);
1061
1062                         if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
1063                                 _vec_normal[0] = 1;
1064                         else//j==16
1065                                 _vec_normal[0] = -1;
1066                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1067                         {
1068                                 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
1069                                 cout << " ; vecteur normal=(";
1070                                 for(int p=0; p<_Ndim; p++)
1071                                         cout << _vec_normal[p] << ",";
1072                                 cout << "). "<<endl;
1073                         }
1074                         for(int i=0;i<Fj.getNumberOfCells();i++){
1075                                 if(i != largestSectionCellIndex){
1076                                         idn = idCells[i];
1077                                         Ctemp2 = _mesh.getCell(idn);
1078                                         _porosityj=_porosityField(idn);
1079                                         convectionState(idm,idn,false);
1080                                         convectionMatrices();
1081                                         diffusionStateAndMatrices(idm, idn,false);
1082
1083                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1084                                                 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
1085
1086                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
1087                                         _inv_dxj = 1/Ctemp2.getMeasure();
1088
1089                                         addConvectionToSecondMember(idm,idn, false);
1090                                         _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
1091                                         addDiffusionToSecondMember(idm,idn, false);
1092                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
1093                                         addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
1094                                 }
1095                         }
1096                 }
1097                 else
1098                 {
1099                         cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
1100                         _runLogFile->close();
1101                         throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
1102                 }
1103
1104         }
1105
1106         //Contribution from delta t
1107         VecAXPY(_b, 1/_dt, _old);
1108         VecAXPY(_b, -1/_dt, _conservativeVars);
1109
1110         VecAssemblyEnd(_b);
1111         VecCopy(_b,F_X);
1112         VecScale(F_X,-1.);
1113
1114         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1115         {
1116                 cout << "ProblemFluid::computeNewtonRHS : Fin calcul second membre pour PETSc"<<endl;
1117                 VecView(F_X,  PETSC_VIEWER_STDOUT_WORLD);
1118                 cout << endl;
1119         }
1120 }
1121
1122 int ProblemFluid::computeSnesRHS(SNES snes, Vec X, Vec F_X, void *ctx)//Prototype imposé par PETSc
1123 {
1124         ProblemFluid * myProblem = (ProblemFluid *) ctx;
1125         myProblem->computeNewtonRHS( X, F_X);
1126
1127         return 0;
1128 }
1129
1130 void ProblemFluid::computeNewtonJacobian( Vec X, Mat A){//dt is known and will contribute to the jacobian matrix of the Newton scheme
1131
1132         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1133         {
1134                 cout << "ProblemFluid::computeNewtonJacobian : Début calcul Jacobienne schéma Newton pour PETSc, _dt="<<_dt<<endl;
1135                 cout << endl;
1136         }
1137
1138         if(_timeScheme == Explicit){
1139                 MatCreateConstantDiagonal(PETSC_COMM_SELF, _nVar, _nVar, _nVar*_Nmailles, _nVar*_Nmailles,1./_dt, &A);
1140                 return ;
1141         }
1142
1143         MatZeroEntries(A);
1144         VecCopy(X,_conservativeVars);
1145
1146         std::vector< int > idCells(2);
1147         PetscInt idm, idn, size = 1;
1148
1149         long nbFaces = _mesh.getNumberOfFaces();
1150         Face Fj;
1151         Cell Ctemp1,Ctemp2;
1152         string nameOfGroup;
1153
1154         for (int j=0; j<nbFaces;j++){
1155                 Fj = _mesh.getFace(j);
1156                 _isBoundary=Fj.isBorder();
1157                 idCells = Fj.getCellsId();
1158
1159                 // If Fj is on the boundary
1160                 if (_isBoundary)
1161                 {
1162                         for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
1163                         {
1164                                 // compute the normal vector corresponding to face j : from Ctemp1 outward
1165                                 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
1166                                 if (_Ndim >1){
1167                                         for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
1168                                         {//we look for l the index of the face Fj for the cell Ctemp1
1169                                                 if (j == Ctemp1.getFacesId()[l])
1170                                                 {
1171                                                         for (int idim = 0; idim < _Ndim; ++idim)
1172                                                                 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
1173                                                         break;
1174                                                 }
1175                                         }
1176                                 }else{ // _Ndim = 1, build normal vector (bug cdmath)
1177                                         if(!_sectionFieldSet)
1178                                         {
1179                                                 if (Fj.x()<Ctemp1.x())
1180                                                         _vec_normal[0] = -1;
1181                                                 else
1182                                                         _vec_normal[0] = 1;
1183                                         }
1184                                         else
1185                                         {
1186                                                 if(idCells[0]==0)
1187                                                         _vec_normal[0] = -1;
1188                                                 else//idCells[0]==31
1189                                                         _vec_normal[0] = 1;
1190                                         }
1191                                 }
1192                                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1193                                 {
1194                                         cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
1195                                         for(int p=0; p<_Ndim; p++)
1196                                                 cout << _vec_normal[p] << ",";
1197                                         cout << "). "<<endl;
1198                                 }
1199                                 nameOfGroup = Fj.getGroupName();
1200                                 _porosityi=_porosityField(idCells[k]);
1201                                 _porosityj=_porosityi;
1202                                 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
1203                                 convectionState(idCells[k],0,true);
1204                                 convectionMatrices();
1205                                 diffusionStateAndMatrices(idCells[k], 0, true);
1206                                 // compute 1/dxi
1207                                 if (_Ndim > 1)
1208                                         _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
1209                                 else
1210                                         _inv_dxi = 1/Ctemp1.getMeasure();
1211
1212                                         for(int l=0; l<_nVar*_nVar;l++){
1213                                                 _AroeMinusImplicit[l] *= _inv_dxi;
1214                                                 _Diffusion[l] *=_inv_dxi*_inv_dxi;
1215                                         }
1216
1217                                         jacobian(idCells[k],nameOfGroup,_vec_normal);
1218                                         jacobianDiff(idCells[k],nameOfGroup);
1219                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
1220                                                 cout << "Matrice Jacobienne CL convection:" << endl;
1221                                                 for(int p=0; p<_nVar; p++){
1222                                                         for(int q=0; q<_nVar; q++)
1223                                                                 cout << _Jcb[p*_nVar+q] << "\t";
1224                                                         cout << endl;
1225                                                 }
1226                                                 cout << endl;
1227                                                 cout << "Matrice Jacobienne CL diffusion:" << endl;
1228                                                 for(int p=0; p<_nVar; p++){
1229                                                         for(int q=0; q<_nVar; q++)
1230                                                                 cout << _JcbDiff[p*_nVar+q] << "\t";
1231                                                         cout << endl;
1232                                                 }
1233                                                 cout << endl;
1234                                         }
1235                                         idm = idCells[k];
1236                                         //calcul et insertion de A^-*Jcb
1237                                         Polynoms::matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
1238                                         MatSetValuesBlocked(A, size, &idm, size, &idm, _a, ADD_VALUES);
1239
1240                                         if(_verbose)
1241                                                 displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
1242
1243                                         //insertion de -A^-
1244                                         for(int k=0; k<_nVar*_nVar;k++){
1245                                                 _AroeMinusImplicit[k] *= -1;
1246                                         }
1247                                         MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
1248                                         if(_verbose)
1249                                                 displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
1250
1251                                         //calcul et insertion de D*JcbDiff
1252                                         Polynoms::matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
1253                                         MatSetValuesBlocked(A, size, &idm, size, &idm, _a, ADD_VALUES);
1254                                         for(int k=0; k<_nVar*_nVar;k++)
1255                                                 _Diffusion[k] *= -1;
1256                                         MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
1257                                 }
1258                 } else  if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
1259                         // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
1260                         Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
1261                         Ctemp2 = _mesh.getCell(idCells[1]);
1262                         if (_Ndim >1){
1263                                 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
1264                                         if (j == Ctemp1.getFacesId()[l]){
1265                                                 for (int idim = 0; idim < _Ndim; ++idim)
1266                                                         _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
1267                                                 break;
1268                                         }
1269                                 }
1270                         }else{ // _Ndim = 1, build normal vector (bug cdmath)
1271                                 if(!_sectionFieldSet)
1272                                 {
1273                                         if (Fj.x()<Ctemp1.x())
1274                                                 _vec_normal[0] = -1;
1275                                         else
1276                                                 _vec_normal[0] = 1;
1277                                 }
1278                                 else
1279                                 {
1280                                         if(idCells[0]>idCells[1])
1281                                                 _vec_normal[0] = -1;
1282                                         else
1283                                                 _vec_normal[0] = 1;
1284                                 }
1285                         }
1286                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1287                         {
1288                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
1289                                 cout<<" Normal vector= ";
1290                                 for (int idim = 0; idim < _Ndim; ++idim)
1291                                         cout<<_vec_normal[idim]<<", ";
1292                                 cout<<endl;
1293                         }
1294                         _porosityi=_porosityField(idCells[0]);
1295                         _porosityj=_porosityField(idCells[1]);
1296                         convectionState(idCells[0],idCells[1],false);
1297                         convectionMatrices();
1298                         diffusionStateAndMatrices(idCells[0], idCells[1], false);
1299
1300                         // compute 1/dxi and 1/dxj
1301                         if (_Ndim > 1)
1302                         {
1303                                 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
1304                                 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
1305                         }
1306                         else
1307                         {
1308                                 _inv_dxi = 1/Ctemp1.getMeasure();
1309                                 _inv_dxj = 1/Ctemp2.getMeasure();
1310                         }
1311
1312                                 for(int k=0; k<_nVar*_nVar;k++)
1313                                 {
1314                                         _AroeMinusImplicit[k] *= _inv_dxi;
1315                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
1316                                 }
1317                                 idm = idCells[0];
1318                                 idn = idCells[1];
1319                                 //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNbCells<<endl;
1320                                 MatSetValuesBlocked(A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
1321                                 MatSetValuesBlocked(A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
1322
1323                                 if(_verbose){
1324                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
1325                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
1326                                 }
1327                                 for(int k=0;k<_nVar*_nVar;k++){
1328                                         _AroeMinusImplicit[k] *= -1;
1329                                         _Diffusion[k] *= -1;
1330                                 }
1331                                 MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
1332                                 MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
1333                                 if(_verbose){
1334                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
1335                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
1336                                 }
1337                                 for(int k=0; k<_nVar*_nVar;k++)
1338                                 {
1339                                         _AroePlusImplicit[k]  *= _inv_dxj;
1340                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;
1341                                 }
1342                                 MatSetValuesBlocked(A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
1343                                 MatSetValuesBlocked(A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
1344                                 if(_verbose)
1345                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
1346
1347                                 for(int k=0;k<_nVar*_nVar;k++){
1348                                         _AroePlusImplicit[k] *= -1;
1349                                         _Diffusion[k] *= -1;
1350                                 }
1351                                 MatSetValuesBlocked(A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
1352                                 MatSetValuesBlocked(A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
1353
1354                                 if(_verbose)
1355                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
1356                 }
1357                 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
1358                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1359                                 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
1360                         *_runLogFile<<"Warning: treatment of a junction node"<<endl;
1361
1362                         if(!_sectionFieldSet)
1363                         {
1364                                 _runLogFile->close();
1365                                 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
1366                         }
1367                         int largestSectionCellIndex=0;
1368                         for(int i=1;i<Fj.getNumberOfCells();i++){
1369                                 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
1370                                         largestSectionCellIndex=i;
1371                         }
1372                         idm = idCells[largestSectionCellIndex];
1373                         Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
1374                         _porosityi=_porosityField(idm);
1375
1376                         if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
1377                                 _vec_normal[0] = 1;
1378                         else//j==16
1379                                 _vec_normal[0] = -1;
1380                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1381                         {
1382                                 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
1383                                 cout << " ; vecteur normal=(";
1384                                 for(int p=0; p<_Ndim; p++)
1385                                         cout << _vec_normal[p] << ",";
1386                                 cout << "). "<<endl;
1387                         }
1388                         for(int i=0;i<Fj.getNumberOfCells();i++){
1389                                 if(i != largestSectionCellIndex){
1390                                         idn = idCells[i];
1391                                         Ctemp2 = _mesh.getCell(idn);
1392                                         _porosityj=_porosityField(idn);
1393                                         convectionState(idm,idn,false);
1394                                         convectionMatrices();
1395                                         diffusionStateAndMatrices(idm, idn,false);
1396
1397                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1398                                                 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
1399
1400                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
1401                                         _inv_dxj = 1/Ctemp2.getMeasure();
1402
1403                                                 for(int k=0; k<_nVar*_nVar;k++)
1404                                                 {
1405                                                         _AroeMinusImplicit[k] *= _inv_dxi;
1406                                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
1407                                                 }
1408                                                 MatSetValuesBlocked(A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
1409                                                 MatSetValuesBlocked(A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
1410
1411                                                 if(_verbose){
1412                                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
1413                                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
1414                                                 }
1415                                                 for(int k=0;k<_nVar*_nVar;k++){
1416                                                         _AroeMinusImplicit[k] *= -1;
1417                                                         _Diffusion[k] *= -1;
1418                                                 }
1419                                                 MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
1420                                                 MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
1421                                                 if(_verbose){
1422                                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
1423                                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
1424                                                 }
1425                                                 for(int k=0; k<_nVar*_nVar;k++)
1426                                                 {
1427                                                         _AroePlusImplicit[k] *= _inv_dxj;
1428                                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
1429                                                 }
1430                                                 MatSetValuesBlocked(A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
1431                                                 MatSetValuesBlocked(A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
1432                                                 if(_verbose)
1433                                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
1434
1435                                                 for(int k=0;k<_nVar*_nVar;k++){
1436                                                         _AroePlusImplicit[k] *= -1;
1437                                                         _Diffusion[k] *= -1;
1438                                                 }
1439                                                 MatSetValuesBlocked(A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
1440                                                 MatSetValuesBlocked(A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
1441
1442                                                 if(_verbose)
1443                                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
1444                                         }
1445                                 }
1446                 }
1447                 else
1448                 {
1449                         cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
1450                         _runLogFile->close();
1451                         throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
1452                 }
1453
1454         }
1455
1456         for(int imaille = 0; imaille<_Nmailles; imaille++)
1457                 MatSetValuesBlocked(A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
1458
1459         MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1460         MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1461
1462         MatShift(A, 1/_dt);
1463
1464         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1465         {
1466                 cout << "ProblemFluid::computeNewtonJacobian : Fin calcul Jacobienne schéma Newton pour PETSc"<<endl;
1467                 MatView(A,PETSC_VIEWER_STDOUT_SELF);
1468                 cout << endl;
1469         }
1470 }
1471
1472 int ProblemFluid::computeSnesJacobian(SNES snes, Vec X, Mat A, Mat Aapprox, void *ctx)//Propotype imposé par PETSc
1473 {
1474         ProblemFluid * myProblem = (ProblemFluid *) ctx;
1475         myProblem->computeNewtonJacobian( X, A);
1476
1477         Aapprox = A;
1478
1479         return 0;
1480 }
1481
1482 void ProblemFluid::validateTimeStep()
1483 {
1484         if(_system)
1485         {
1486                 cout<<" Vecteur Un"<<endl;
1487                 VecView(_old,  PETSC_VIEWER_STDOUT_WORLD);
1488                 cout<<" Vecteur Un+1"<<endl;
1489                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
1490         }
1491         VecAXPY(_old,  -1, _conservativeVars);//old contient old-courant
1492
1493         //Calcul de la variation Un+1-Un
1494         _erreur_rel= 0;
1495         double x, dx;
1496         int I;
1497
1498         for(int j=1; j<=_Nmailles; j++)
1499         {
1500                 for(int k=0; k<_nVar; k++)
1501                 {
1502                         I = (j-1)*_nVar + k;
1503                         VecGetValues(_old, 1, &I, &dx);
1504                         VecGetValues(_conservativeVars, 1, &I, &x);
1505                         if (fabs(x)< _precision)
1506                         {
1507                                 if(_erreur_rel < fabs(dx))
1508                                         _erreur_rel = fabs(dx);
1509                         }
1510                         else if(_erreur_rel < fabs(dx/x))
1511                                 _erreur_rel = fabs(dx/x);
1512                 }
1513         }
1514
1515         _isStationary =_erreur_rel <_precision;
1516
1517         VecCopy(_conservativeVars, _old);
1518
1519         if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
1520                 if(!_usePrimitiveVarsInNewton)
1521                         testConservation();
1522                 cout <<"Valeur propre locale max: " << _maxvp << endl;
1523         }
1524
1525         if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
1526                 //Find minimum and maximum void fractions
1527                 double alphamin=1e30;
1528                 double alphamax=-1e30;
1529                 double T, Tmax=-1e30;
1530                 int J=_nVar-1;
1531                 I = 0;
1532                 for(int j=0; j<_Nmailles; j++)
1533                 {
1534                         VecGetValues(_primitiveVars, 1, &I, &x);//extract void fraction
1535                         if(x>alphamax)
1536                                 alphamax=x;
1537                         if(x<alphamin)
1538                                 alphamin=x;
1539                         I += _nVar;
1540                         VecGetValues(_primitiveVars, 1, &J, &T);//extract void fraction
1541                         if(T>Tmax)
1542                                 Tmax=T;
1543                         J += _nVar;
1544                 }
1545                 cout<<"Alpha min = " << alphamin << ", Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
1546                 cout<<endl;
1547                 *_runLogFile<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
1548         }
1549
1550         _time+=_dt;
1551         _nbTimeStep++;
1552         if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
1553                 save();
1554 }
1555
1556 void ProblemFluid::abortTimeStep(){
1557         _dt = 0;
1558 }
1559
1560 void ProblemFluid::addConvectionToSecondMember
1561 (               const int &i,
1562                 const int &j, bool isBord, string groupname
1563 )
1564 {
1565         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1566                 cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
1567
1568         //extraction des valeurs
1569         for(int k=0; k<_nVar; k++)
1570                 _idm[k] = _nVar*i + k;
1571         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1572         VecGetValues(_primitiveVars,    _nVar, _idm, _Vi);
1573
1574         if(!isBord){
1575                 for(int k=0; k<_nVar; k++)
1576                         _idn[k] = _nVar*j + k;
1577                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1578                 VecGetValues(_primitiveVars,    _nVar, _idn, _Vj);
1579         }
1580         else{
1581                 for(int k=0; k<_nVar; k++)
1582                         _idn[k] = k;
1583                 VecGetValues(_Uext, _nVar, _idn, _Uj);
1584                 consToPrim(_Uj, _Vj,_porosityj);
1585         }
1586         _idm[0] = i;
1587         _idn[0] = j;
1588
1589         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1590         {
1591                 cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
1592                 for(int q=0; q<_nVar; q++)
1593                         cout << _Ui[q] << endl;
1594                 cout << endl;
1595                 cout << "addConvectionToSecondMember : état j= " << j << ", _Uj=" << endl;
1596                 for(int q=0; q<_nVar; q++)
1597                         cout << _Uj[q] << endl;
1598                 cout << endl;
1599         }
1600         if(_nonLinearFormulation!=reducedRoe){//VFRoe, Roe or VFFC
1601                 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar),Fij(_nVar);
1602                 for(int i1=0;i1<_nVar;i1++)
1603                 {
1604                         Ui(i1)=_Ui[i1];
1605                         Uj(i1)=_Uj[i1];
1606                         Vi(i1)=_Vi[i1];
1607                         Vj(i1)=_Vj[i1];
1608                 }
1609                 Vector normale(_Ndim);
1610                 for(int i1=0;i1<_Ndim;i1++)
1611                         normale(i1)=_vec_normal[i1];
1612
1613                 Matrix signAroe(_nVar);
1614                 for(int i1=0;i1<_nVar;i1++)
1615                         for(int i2=0;i2<_nVar;i2++)
1616                                 signAroe(i1,i2)=_signAroe[i1*_nVar+i2];
1617
1618                 Matrix absAroe(_nVar);
1619                 for(int i1=0;i1<_nVar;i1++)
1620                         for(int i2=0;i2<_nVar;i2++)
1621                                 absAroe(i1,i2)=_absAroe[i1*_nVar+i2];
1622
1623                 if(_nonLinearFormulation==VFRoe)//VFRoe
1624                 {
1625                         Vector Uij(_nVar),Vij(_nVar);
1626                         double porosityij=sqrt(_porosityi*_porosityj);
1627
1628                         Uij=(Ui+Uj)/2+signAroe*(Ui-Uj)/2;
1629
1630                         for(int i1=0;i1<_nVar;i1++)
1631                                 _Uij[i1]=Uij(i1);
1632
1633                         consToPrim(_Uij, _Vij,porosityij);
1634
1635                         applyVFRoeLowMachCorrections(isBord, groupname);
1636
1637                         for(int i1=0;i1<_nVar;i1++)
1638                         {
1639                                 Uij(i1)=_Uij[i1];
1640                                 Vij(i1)=_Vij[i1];
1641                         }
1642
1643                         Fij=convectionFlux(Uij,Vij,normale,porosityij);
1644
1645                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1646                         {
1647                                 cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
1648                                 cout<<Uij<<endl;
1649                                 cout<<"Etat interfacial primitif "<<i<<", "<<j<< endl;
1650                                 cout<<Vij<<endl;
1651                         }
1652                 }
1653                 else //Roe or VFFC
1654                 {
1655                         if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
1656                                 Fij=staggeredVFFCFlux();
1657                         else
1658                         {
1659                                 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
1660                                 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
1661                                 if(_nonLinearFormulation==VFFC)//VFFC
1662                                         Fij=(Fi+Fj)/2+signAroe*(Fi-Fj)/2;
1663                                 else if(_nonLinearFormulation==Roe)//Roe
1664                                         Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
1665
1666                                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1667                                 {
1668                                         cout<<"Flux cellule "<<i<<" = "<< endl;
1669                                         cout<<Fi<<endl;
1670                                         cout<<"Flux cellule "<<j<<" = "<< endl;
1671                                         cout<<Fj<<endl;
1672                                 }
1673                         }
1674                 }
1675                 for(int i1=0;i1<_nVar;i1++)
1676                         _phi[i1]=-Fij(i1)*_inv_dxi;
1677                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1678                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1679                 {
1680                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1681                         cout<<"Flux interfacial "<<i<<", "<<j<< endl;
1682                         cout<<Fij<<endl;
1683                         cout << "Contribution convection à " << i << ", -Fij(i1)*_inv_dxi= "<<endl;
1684                         for(int q=0; q<_nVar; q++)
1685                                 cout << _phi[q] << endl;
1686                         cout << endl;
1687                 }
1688                 if(!isBord){
1689                         for(int i1=0;i1<_nVar;i1++)
1690                                 _phi[i1]*=-_inv_dxj/_inv_dxi;
1691                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1692                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1693                         {
1694                                 cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
1695                                 for(int q=0; q<_nVar; q++)
1696                                         cout << _phi[q] << endl;
1697                                 cout << endl;
1698                         }
1699                 }
1700         }else //_nonLinearFormulation==reducedRoe)
1701         {
1702                 for(int k=0; k<_nVar; k++)
1703                         _temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
1704                 Polynoms::matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
1705                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1706
1707                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1708                 {
1709                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1710                         cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
1711                         for(int q=0; q<_nVar; q++)
1712                                 cout << _temp[q] << endl;
1713                         cout << endl;
1714                         cout << "Contribution convection à " << i << ", A^-*(Ui - Uj)*_inv_dxi= "<<endl;
1715                         for(int q=0; q<_nVar; q++)
1716                                 cout << _phi[q] << endl;
1717                         cout << endl;
1718                 }
1719
1720                 if(!isBord)
1721                 {
1722                         for(int k=0; k<_nVar; k++)
1723                                 _temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
1724                         Polynoms::matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
1725                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1726
1727                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1728                         {
1729                                 cout << "Contribution convection à  " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
1730                                 for(int q=0; q<_nVar; q++)
1731                                         cout << _phi[q] << endl;
1732                                 cout << endl;
1733                         }
1734                 }
1735         }
1736         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1737         {
1738                 cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
1739                 displayMatrix(_absAroe,   _nVar,"Valeur absolue matrice de Roe");
1740                 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
1741                 displayMatrix(_AroePlus,  _nVar,"Matrice _AroePlus");
1742                 displayMatrix(_signAroe,  _nVar,"Signe de la matrice de Roe");
1743         }
1744 }
1745
1746 void ProblemFluid::addSourceTermToSecondMember
1747 (       const int i, int nbVoisinsi,
1748                 const int j, int nbVoisinsj,
1749                 bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
1750 {
1751         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1752                 cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
1753
1754         _idm[0] = i*_nVar;
1755         for(int k=1; k<_nVar;k++)
1756                 _idm[k] = _idm[k-1] + 1;
1757         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1758         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1759         sourceVector(_Si,_Ui,_Vi,i);
1760
1761         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1762         {
1763                 cout << "Terme source S(Ui), i= " << i<<endl;
1764                 for(int q=0; q<_nVar; q++)
1765                         cout << _Si[q] << endl;
1766                 cout << endl;
1767         }
1768         if(!isBord){
1769                 for(int k=0; k<_nVar; k++)
1770                         _idn[k] = _nVar*j + k;
1771                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1772                 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1773                 sourceVector(_Sj,_Uj,_Vj,j);
1774         }else
1775         {
1776                 for(int k=0; k<_nVar; k++)
1777                         _idn[k] = k;
1778                 VecGetValues(_Uext, _nVar, _idn, _Uj);
1779                 consToPrim(_Uj, _Vj,_porosityj);
1780                 sourceVector(_Sj,_Uj,_Vj,i);
1781         }
1782         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1783         {
1784                 if(!isBord)
1785                         cout << "Terme source S(Uj), j= " << j<<endl;
1786                 else
1787                         cout << "Terme source S(Uj), cellule fantôme "<<endl;
1788                 for(int q=0; q<_nVar; q++)
1789                         cout << _Sj[q] << endl;
1790                 cout << endl;
1791         }
1792
1793         //Compute pressure loss vector
1794         double K;
1795         if(_pressureLossFieldSet){
1796                 K=_pressureLossField(ij);
1797                 pressureLossVector(_pressureLossVector, K,_Ui, _Vi, _Uj, _Vj);
1798         }
1799         else{
1800                 K=0;
1801                 for(int k=0; k<_nVar;k++)
1802                         _pressureLossVector[k]=0;
1803         }
1804         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1805         {
1806                 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
1807                 for(int k=0; k<_nVar;k++)
1808                         cout<< _pressureLossVector[k]<<", ";
1809                 cout<<endl;
1810         }
1811         //Contribution of the porosityField gradient:
1812         if(!isBord)
1813                 porosityGradientSourceVector();
1814         else{
1815                 for(int k=0; k<_nVar;k++)
1816                         _porosityGradientSourceVector[k]=0;
1817         }
1818
1819         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1820         {
1821                 if(!isBord)
1822                         cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
1823                 else
1824                         cout<<"interface frontière i= "<<i<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<endl;
1825                 cout<<"Gradient de porosite à l'interface"<<endl;
1826                 for(int k=0; k<_nVar;k++)
1827                         cout<< _porosityGradientSourceVector[k]<<", ";
1828                 cout<<endl;
1829         }
1830
1831         if(!isBord){
1832                 if(_wellBalancedCorrection){
1833                         for(int k=0; k<_nVar;k++)
1834                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1835                         Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1836                         for(int k=0; k<_nVar;k++){
1837                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1838                                 _Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
1839                         }
1840                         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1841                         {
1842                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,j), j="<<j<<endl;
1843                                 for(int q=0; q<_nVar; q++)
1844                                         cout << _Si[q] << endl;
1845                                 cout << "Contribution au terme source Sj de la cellule j= " << j<<" venant  (après décentrement) de la face (i,j), i="<<i<<endl;
1846                                 for(int q=0; q<_nVar; q++)
1847                                         cout << _Sj[q] << endl;
1848                                 cout << endl;
1849                                 cout<<"ratio surface sur volume i = "<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1850                                 cout<<"ratio surface sur volume j = "<<mesureFace/_perimeters(j)<<" perimeter = "<< _perimeters(j)<<endl;
1851                                 cout << endl;
1852                         }
1853                 }
1854                 else{
1855                         for(int k=0; k<_nVar;k++){
1856                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
1857                                 _Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
1858                         }
1859                         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1860                         {
1861                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant  de la face (i,j), j="<<j<<endl;
1862                                 for(int q=0; q<_nVar; q++)
1863                                         cout << _Si[q] << endl;
1864                                 cout << "Contribution au terme source Sj de la cellule j = " << j<<" venant  de la face (i,j), i="<<i <<endl;
1865                                 for(int q=0; q<_nVar; q++)
1866                                         cout << _Sj[q] << endl;
1867                                 cout << endl;
1868                         }
1869                 }
1870                 _idn[0] = j;
1871                 VecSetValuesBlocked(_b, 1, _idn, _Sj, ADD_VALUES);
1872         }else{
1873                 if(_wellBalancedCorrection){
1874                         for(int k=0; k<_nVar;k++)
1875                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1876                         Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1877                         for(int k=0; k<_nVar;k++)
1878                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1879                         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1880                         {
1881                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,bord)"<<endl;
1882                                 for(int q=0; q<_nVar; q++)
1883                                         cout << _Si[q] << endl;
1884                                 cout<<"ratio surface sur volume i ="<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1885                                 cout << endl;
1886                         }
1887                 }
1888                 else
1889                 {
1890                         for(int k=0; k<_nVar;k++)
1891                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
1892                         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1893                         {
1894                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
1895                                 for(int q=0; q<_nVar; q++)
1896                                         cout << _Si[q] << endl;
1897                                 cout << endl;
1898                         }
1899                 }
1900         }
1901         _idm[0] = i;
1902         VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
1903
1904         if(_verbose && (_nbTimeStep-1)%_freqSave ==0 && _wellBalancedCorrection)
1905                 displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
1906 }
1907
1908 void ProblemFluid::updatePrimitives()
1909 {
1910         VecAssemblyBegin(_primitiveVars);
1911         for(int i=1; i<=_Nmailles; i++)
1912         {
1913                 _idm[0] = _nVar*( (i-1));
1914                 for(int k=1; k<_nVar; k++)
1915                         _idm[k] = _idm[k-1] + 1;
1916
1917                 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1918                 consToPrim(_Ui,_Vi,_porosityField(i-1));
1919                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1920                 {
1921                         cout << "ProblemFluid::updatePrimitives() cell " << i-1 << endl;
1922                         cout << "Ui = ";
1923                         for(int q=0; q<_nVar; q++)
1924                                 cout << _Ui[q]  << "\t";
1925                         cout << endl;
1926                         cout << "Vi = ";
1927                         for(int q=0; q<_nVar; q++)
1928                                 cout << _Vi[q]  << "\t";
1929                         cout << endl;
1930                 }
1931
1932                 if(_nbPhases==2 && _Psat>-1e30){//Cas simulation flashing
1933                         double pressure;
1934                         VecGetValues(_primitiveVars, 1, _idm+1, &pressure);
1935                         _dp_over_dt(i-1)=(_Vi[1]-pressure)/_dt;//pn+1-pn
1936                 }
1937                 _idm[0] = i-1;
1938                 VecSetValuesBlocked(_primitiveVars, 1, _idm, _Vi, INSERT_VALUES);
1939         }
1940         VecAssemblyEnd(_primitiveVars);
1941
1942         if(_system)
1943         {
1944                 cout << "Nouvelles variables primitives : " << endl;
1945                 VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_WORLD);
1946                 cout << endl;
1947         }
1948 }
1949
1950 void ProblemFluid::updateConservatives()
1951 {
1952         VecAssemblyBegin(_conservativeVars);
1953         for(int i=1; i<=_Nmailles; i++)
1954         {
1955                 _idm[0] = _nVar*( (i-1));
1956                 for(int k=1; k<_nVar; k++)
1957                         _idm[k] = _idm[k-1] + 1;
1958
1959                 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1960
1961                 primToCons(_Vi,0,_Ui,0);
1962                 _idm[0] = i-1;
1963                 VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
1964
1965                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1966                 {
1967                         cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
1968                         cout << "Vi = ";
1969                         for(int q=0; q<_nVar; q++)
1970                                 cout << _Vi[q]  << "\t";
1971                         cout << endl;
1972                         cout << "Ui = ";
1973                         for(int q=0; q<_nVar; q++)
1974                                 cout << _Ui[q]  << "\t";
1975                         cout << endl;
1976                 }
1977         }
1978         VecAssemblyEnd(_conservativeVars);
1979
1980         if(_system)
1981         {
1982                 cout << "Nouvelles variables conservatives : " << endl;
1983                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
1984                 cout << endl;
1985         }
1986 }
1987
1988 vector< complex<double> > ProblemFluid::getRacines(vector< double > pol_car){
1989         int degre_polynom = pol_car.size()-1;
1990         double tmp;
1991         vector< complex< double > > valeurs_propres (degre_polynom);
1992         vector< double > zeror(degre_polynom);
1993         vector< double > zeroi(degre_polynom);
1994         for(int j=0; j<(degre_polynom+1)/2; j++){//coefficients in order of decreasing powers for rpoly
1995                 tmp=pol_car[j];
1996                 pol_car[j] =pol_car[degre_polynom-j];
1997                 pol_car[degre_polynom-j]=tmp;
1998         }
1999
2000         //Calcul des racines du polynome
2001         roots_polynoms roots;
2002         int size_vp = roots.rpoly(&pol_car[0],degre_polynom,&zeror[0],&zeroi[0]);
2003
2004         //On ordonne les valeurs propres
2005         if(zeror[1]<zeror[0])
2006         {
2007                 tmp=zeror[0];
2008                 zeror[0]=zeror[1];
2009                 zeror[1]=tmp;
2010                 tmp=zeroi[0];
2011                 zeroi[0]=zeroi[1];
2012                 zeroi[1]=tmp;
2013         }
2014
2015         if(size_vp == degre_polynom) {
2016                 for(int ct =0; ct<degre_polynom; ct++)
2017                         valeurs_propres[ct] = complex<double> (zeror[ct],zeroi[ct]);  //vp non triviales
2018         }
2019         else {
2020                 cout << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
2021                 *_runLogFile << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
2022                 cout<<"Coefficients polynome caracteristique: "<<endl;
2023                 for(int ct =0; ct<degre_polynom+1; ct++)
2024                         cout<<pol_car[ct]<< " , " <<endl;
2025
2026                 *_runLogFile<<"getRacines computation failed"<<endl;
2027                 _runLogFile->close();
2028                 throw CdmathException("getRacines computation failed");
2029         }
2030
2031         return valeurs_propres;
2032 }
2033
2034 void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
2035 {
2036         int nbVp_dist=valeurs_propres_dist.size();
2037         vector< complex< double > > y (nbVp_dist,0);
2038         for( int i=0 ; i<nbVp_dist ; i++)
2039                 y[i] = Polynoms::abs_generalise(valeurs_propres_dist[i]);
2040         Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
2041         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2042         {
2043                 cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
2044                 for(int ct =0; ct<nbVp_dist; ct++)
2045                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
2046                 cout<< endl;
2047                 cout<<"ProblemFluid::AbsMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
2048                 for(int ct =0; ct<nbVp_dist; ct++)
2049                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
2050                 cout<< endl;
2051                 displayMatrix(_absAroe,_nVar,"Valeur absolue de la matrice de Roe");
2052         }
2053 }
2054
2055 void ProblemFluid::SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist)
2056 {
2057         int nbVp_dist=valeurs_propres_dist.size();
2058         vector< complex< double > > y (nbVp_dist,0);
2059         for( int i=0 ; i<nbVp_dist ; i++)
2060         {
2061                 if(valeurs_propres_dist[i].real()>0)
2062                         y[i] = 1;
2063                 else if(valeurs_propres_dist[i].real()<0)
2064                         y[i] = -1;
2065                 else
2066                         y[i] = 0;
2067         }
2068
2069         Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
2070         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2071         {
2072                 cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
2073                 for(int ct =0; ct<nbVp_dist; ct++)
2074                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
2075                 cout<< endl;
2076                 cout<<"ProblemFluid::SigneMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
2077                 for(int ct =0; ct<nbVp_dist; ct++)
2078                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
2079                 cout<< endl;
2080                 displayMatrix(_signAroe,_nVar,"Signe matrice de Roe");
2081         }
2082 }
2083 void ProblemFluid::InvMatriceRoe(vector< complex<double> > valeurs_propres_dist)
2084 {
2085         int nbVp_dist=valeurs_propres_dist.size();
2086         vector< complex< double > > y (nbVp_dist,0);
2087
2088         for( int i=0 ; i<nbVp_dist ; i++)
2089         {
2090                 if(Polynoms::module(valeurs_propres_dist[i])>_precision)
2091                         y[i] = 1./valeurs_propres_dist[i];
2092                 else
2093                         y[i] = 1./_precision;
2094         }
2095         Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
2096         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2097         {
2098                 cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
2099                 for(int ct =0; ct<nbVp_dist; ct++)
2100                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
2101                 cout<< endl;
2102                 cout<<"ProblemFluid::InvMatriceRoe : Valeurs à atteindre pour l'interpolation"<<endl;
2103                 for(int ct =0; ct<nbVp_dist; ct++)
2104                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
2105                 cout<< endl;
2106                 displayMatrix(_invAroe,_nVar,"Inverse matrice de Roe");
2107         }
2108 }
2109
2110 Field ProblemFluid::getConservativeField() const
2111 {
2112         if(!_initializedMemory)
2113         {
2114                 _runLogFile->close();
2115                 throw CdmathException("ProblemFluid::getConservativeField call initialize() first");
2116         }
2117         return _UU;
2118 }
2119
2120 Field ProblemFluid::getPrimitiveField() const
2121 {
2122         if(!_initializedMemory)
2123         {
2124                 _runLogFile->close();
2125                 throw CdmathException("ProblemFluid::getPrimitiveField call initialize() first");
2126         }
2127         return _VV;
2128 }
2129
2130 void ProblemFluid::terminate(){
2131
2132         delete[] _AroePlus;
2133         delete[] _Diffusion;
2134         delete[] _GravityImplicitationMatrix;
2135         delete[] _AroeMinus;
2136         delete[] _Aroe;
2137         delete[] _absAroe;
2138         delete[] _signAroe;
2139         delete[] _invAroe;
2140         delete[] _AroeImplicit;
2141         delete[] _AroeMinusImplicit;
2142         delete[] _AroePlusImplicit;
2143         delete[] _absAroeImplicit;
2144         delete[] _phi;
2145         delete[] _Jcb;
2146         delete[] _JcbDiff;
2147         delete[] _a;
2148         delete[] _primToConsJacoMat;
2149
2150         delete[] _l;
2151         delete[] _r;
2152         delete[] _Uroe;
2153         delete[] _Udiff;
2154         delete[] _temp;
2155         delete[] _externalStates;
2156         delete[] _idm;
2157         delete[] _idn;
2158         delete[] _vec_normal;
2159         delete[] _Ui;
2160         delete[] _Uj;
2161         delete[] _Vi;
2162         delete[] _Vj;
2163         if(_nonLinearFormulation==VFRoe){
2164                 delete[] _Uij;
2165                 delete[] _Vij;
2166         }
2167         delete[] _Si;
2168         delete[] _Sj;
2169         delete[] _pressureLossVector;
2170         delete[] _porosityGradientSourceVector;
2171         if(_isScaling)
2172         {
2173                 delete[] _blockDiag;
2174                 delete[] _invBlockDiag;
2175
2176                 VecDestroy(&_vecScaling);
2177                 VecDestroy(&_invVecScaling);
2178                 VecDestroy(&_bScaling);
2179         }
2180
2181         VecDestroy(&_conservativeVars);
2182         VecDestroy(&_newtonVariation);
2183         VecDestroy(&_b);
2184         VecDestroy(&_primitiveVars);
2185         VecDestroy(&_Uext);
2186         VecDestroy(&_Vext);
2187         VecDestroy(&_Uextdiff);
2188
2189         //      PCDestroy(_pc);
2190         KSPDestroy(&_ksp);
2191         for(int i=0;i<_nbPhases;i++)
2192                 delete _fluides[i];
2193
2194         // Destruction du solveur de Newton de PETSc
2195         if(_timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
2196                 SNESDestroy(&_snes);
2197 }
2198
2199 vector<string>
2200 ProblemFluid::getInputFieldsNames()
2201 {
2202         vector<string> result(1);
2203
2204         result[0]="HeatPower";
2205         result[1]="Porosity";
2206         result[2]="PressureLoss";
2207         result[3]="Section";
2208
2209         return result;
2210 }
2211
2212 void
2213 ProblemFluid::setInputField(const string& nameField, Field& inputField )
2214 {
2215         if(nameField=="Porosity" || nameField=="POROSITY" || nameField=="Porosité" || nameField=="POROSITE")
2216                 return setPorosityField( inputField) ;
2217         else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
2218                 return setHeatPowerField( inputField );
2219         else    if(nameField=="PressureLoss" || nameField=="PRESSURELOSS" || nameField=="PerteDeCharge" || nameField=="PPERTEDECHARGE")
2220                 return setPressureLossField( inputField) ;
2221         else if(nameField=="Section" || nameField=="SECTION" )
2222                 return setSectionField( inputField );
2223         else
2224     {
2225         cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
2226         throw CdmathException("SinglePhase::setInputField error : Unknown Field name");
2227     }
2228 }
2229
2230 void
2231 ProblemFluid::setPorosityField(string fileName, string fieldName){
2232         if(!_initialDataSet)
2233                 throw CdmathException("!!!!!!!! ProblemFluid::setPorosityField set initial field first");
2234
2235         _porosityField=Field(fileName, CELLS,fieldName);
2236         _porosityField.getMesh().checkFastEquivalWith(_mesh);
2237         _porosityFieldSet=true;
2238 }
2239
2240 void
2241 ProblemFluid::setPressureLossField(Field PressureLoss){
2242         if(!_initialDataSet)
2243                 throw CdmathException("!!!!!!!! ProblemFluid::setPressureLossField set initial field first");
2244
2245         PressureLoss.getMesh().checkFastEquivalWith(_mesh);
2246         _pressureLossField=PressureLoss;
2247         _pressureLossFieldSet=true;
2248 }
2249
2250 void
2251 ProblemFluid::setPressureLossField(string fileName, string fieldName){
2252         if(!_initialDataSet)
2253                 throw CdmathException("!!!!!!!! ProblemFluid::setPressureLossField set initial field first");
2254
2255         _pressureLossField=Field(fileName, FACES,fieldName);
2256         _pressureLossField.getMesh().checkFastEquivalWith(_mesh);
2257         _pressureLossFieldSet=true;
2258 }
2259
2260 void
2261 ProblemFluid::setSectionField(Field sectionField){
2262         if(!_initialDataSet)
2263                 throw CdmathException("!!!!!!!! ProblemFluid::setSectionField set initial field first");
2264
2265         sectionField.getMesh().checkFastEquivalWith(_mesh);
2266         _sectionField=sectionField;
2267         _sectionFieldSet=true;
2268 }
2269
2270 void
2271 ProblemFluid::setSectionField(string fileName, string fieldName){
2272         if(!_initialDataSet)
2273                 throw CdmathException("!!!!!!!! ProblemFluid::setSectionField set initial field first");
2274
2275         _sectionField=Field(fileName, CELLS,fieldName);
2276         _sectionField.getMesh().checkFastEquivalWith(_mesh);
2277         _sectionFieldSet=true;
2278 }
2279
2280 void
2281 ProblemFluid::setPorosityField(Field Porosity){
2282         if(!_initialDataSet)
2283                 throw CdmathException("!!!!!!!! ProblemFluid::setPorosityField set initial field first");
2284
2285         Porosity.getMesh().checkFastEquivalWith(_mesh);
2286         _porosityField=Porosity;
2287         _porosityFieldSet=true;
2288 }