CoreFlows/Models/src/ProblemFluid.cxx

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