CoreFlows/Models/src/ProblemFluid.cxx

   1 #include "ProblemFluid.hxx"
   2 #include "Fluide.h"
   3 #include "math.h"
   4
   5 using namespace std;
   6
   7 ProblemFluid::ProblemFluid(void){
   8         _latentHeat=1e30;
   9         _Tsat=1e30;
  10         _Psat=-1e30;
  11         _dHsatl_over_dp=0;
  12         _porosityFieldSet=false;
  13         _pressureLossFieldSet=false;
  14         _sectionFieldSet=false;
  15         _GravityField3d=vector<double>(3,0);
  16         _gravityReferencePoint=vector<double>(3,0);
  17         _Uroe=NULL;_Udiff=NULL;_temp=NULL;_l=NULL;_vec_normal=NULL;;
  18         _idm=NULL;_idn=NULL;
  19         _saveVelocity=false;
  20         _saveConservativeField=false;
  21         _usePrimitiveVarsInNewton=false;
  22         _saveInterfacialField=false;
  23         _err_press_max=0; _part_imag_max=0; _nbMaillesNeg=0; _nbVpCplx=0;_minm1=1e30;_minm2=1e30;//Pour affichage paramètres diphasiques dans IterateTimeStep
  24         _isScaling=false;
  25         _entropicCorrection=false;
  26         _pressureCorrectionOrder=2;
  27         _nonLinearFormulation=Roe;
  28         _maxvploc=0.;
  29         _spaceScheme=upwind;
  30 }
  31
  32 SpaceScheme ProblemFluid::getSpaceScheme() const
  33 {
  34         return _spaceScheme;
  35 }
  36 void ProblemFluid::setNumericalScheme(SpaceScheme spaceScheme, TimeScheme timeScheme)
  37 {
  38         if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
  39                 _restartWithNewTimeScheme=true;
  40         _timeScheme = timeScheme;
  41         _spaceScheme = spaceScheme;
  42 }
  43
  44 void ProblemFluid::initialize()
  45 {
  46         if(!_initialDataSet){
  47                 *_runLogFile<<"ProblemFluid::initialize() set initial data first"<<endl;
  48                 _runLogFile->close();
  49                 throw CdmathException("ProblemFluid::initialize() set initial data first");
  50         }
  51         cout << "Number of Phases = " << _nbPhases << " mesh dimension = "<<_Ndim<<" number of variables = "<<_nVar<<endl;
  52         *_runLogFile << "Number of Phases = " << _nbPhases << " spaceDim= "<<_Ndim<<" number of variables= "<<_nVar<<endl;
  53
  54         /********* local arrays ****************/
  55         _AroePlus = new PetscScalar[_nVar*_nVar];
  56         _Diffusion = new PetscScalar[_nVar*_nVar];
  57         _AroeMinus = new PetscScalar[_nVar*_nVar];
  58         _Aroe = new PetscScalar[_nVar*_nVar];
  59         _absAroe = new PetscScalar[_nVar*_nVar];
  60         _signAroe = new PetscScalar[_nVar*_nVar];
  61         _invAroe = new PetscScalar[_nVar*_nVar];
  62         _AroeImplicit = new PetscScalar[_nVar*_nVar];
  63         _AroeMinusImplicit = new PetscScalar[_nVar*_nVar];
  64         _AroePlusImplicit = new PetscScalar[_nVar*_nVar];
  65         _absAroeImplicit = new PetscScalar[_nVar*_nVar];
  66         _primToConsJacoMat = new PetscScalar[_nVar*_nVar];
  67         _phi = new PetscScalar[_nVar];
  68         _Jcb = new PetscScalar[_nVar*_nVar];
  69         _JcbDiff = new PetscScalar[_nVar*_nVar];
  70         _a = new PetscScalar[_nVar*_nVar];
  71         _externalStates = new PetscScalar[_nVar];
  72         _Ui = new PetscScalar[_nVar];
  73         _Uj = new PetscScalar[_nVar];
  74         _Vi = new PetscScalar[_nVar];
  75         _Vj = new PetscScalar[_nVar];
  76         _Uij = new PetscScalar[_nVar];//used for VFRoe scheme
  77         _Vij = new PetscScalar[_nVar];//used for VFRoe scheme
  78         _Si = new PetscScalar[_nVar];
  79         _Sj = new PetscScalar[_nVar];
  80         _pressureLossVector= new PetscScalar[_nVar];
  81         _porosityGradientSourceVector= new PetscScalar[_nVar];
  82
  83         _l = new double[_nVar];
  84         _r = new double[_nVar];
  85         _Udiff = new double[_nVar];
  86         _temp = new double[_nVar*_nVar];
  87
  88         _idm = new int[_nVar];
  89         _idn = new int[_nVar];
  90         _vec_normal = new double[_Ndim];
  91
  92         for(int k=0;k<_nVar*_nVar;k++)
  93         {
  94                 _Diffusion[k]=0;
  95                 _JcbDiff[k]=0;
  96         }
  97         for(int k=0; k<_nVar; k++){
  98                 _idm[k] = k;
  99                 _idn[k] = k;
 100         }
 101
 102         /**************** Field creation *********************/
 103         if(!_heatPowerFieldSet){
 104                 _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
 105                 for(int i =0; i<_Nmailles; i++)
 106                         _heatPowerField(i) = _heatSource;
 107         }
 108         if(_Psat>-1e30)
 109                 _dp_over_dt=Field("dP/dt",CELLS,_mesh,1);
 110         if(!_porosityFieldSet){
 111                 _porosityField=Field("Porosity field",CELLS,_mesh,1);
 112                 for(int i =0; i<_Nmailles; i++)
 113                         _porosityField(i) = 1;
 114                 _porosityFieldSet=true;
 115         }
 116
 117         //conservative field used only for saving results
 118         _UU=Field ("Conservative vec", CELLS, _mesh, _nVar);
 119         if(_saveInterfacialField && _nonLinearFormulation==VFRoe)
 120         {
 121                 _UUstar=Field ("Interfacial U", CELLS, _mesh, _nVar);
 122                 _VVstar=Field ("Interfacial V", CELLS, _mesh, _nVar);
 123         }
 124
 125         //Construction des champs primitifs et conservatifs initiaux comme avant dans ParaFlow
 126         double * initialFieldPrim = new double[_nVar*_Nmailles];
 127         for(int i =0; i<_Nmailles;i++)
 128                 for(int j =0; j<_nVar;j++)
 129                         initialFieldPrim[i*_nVar+j]=_VV(i,j);
 130         double *initialFieldCons = new double[_nVar*_Nmailles];
 131         for(int i=0; i<_Nmailles; i++)
 132                 primToCons(initialFieldPrim, i, initialFieldCons, i);
 133         for(int i =0; i<_Nmailles;i++)
 134                 for(int j =0; j<_nVar;j++)
 135                         _UU(i,j)=initialFieldCons[i*_nVar+j];
 136
 137         /**********Petsc structures:  ****************/
 138
 139         //creation de la matrice
 140         if(_timeScheme == Implicit)
 141                 MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
 142
 143         //creation des vecteurs
 144         VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uext);
 145         VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Vext);
 146         VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uextdiff);
 147         //        VecCreateSeq(PETSC_COMM_SELF, _nVar*_Nmailles, &_conservativeVars);
 148         VecCreate(PETSC_COMM_SELF, &_conservativeVars);
 149         VecSetSizes(_conservativeVars,PETSC_DECIDE,_nVar*_Nmailles);
 150         VecSetBlockSize(_conservativeVars,_nVar);
 151         VecSetFromOptions(_conservativeVars);
 152         VecDuplicate(_conservativeVars, &_old);
 153         VecDuplicate(_conservativeVars, &_newtonVariation);
 154         VecDuplicate(_conservativeVars, &_b);
 155         VecDuplicate(_conservativeVars, &_primitiveVars);
 156
 157         if(_isScaling)
 158         {
 159                 VecDuplicate(_conservativeVars, &_vecScaling);
 160                 VecDuplicate(_conservativeVars, &_invVecScaling);
 161                 VecDuplicate(_conservativeVars, &_bScaling);
 162
 163                 _blockDiag = new PetscScalar[_nVar];
 164                 _invBlockDiag = new PetscScalar[_nVar];
 165         }
 166
 167
 168         int *indices = new int[_Nmailles];
 169         for(int i=0; i<_Nmailles; i++)
 170                 indices[i] = i;
 171         VecSetValuesBlocked(_conservativeVars, _Nmailles, indices, initialFieldCons, INSERT_VALUES);
 172         VecAssemblyBegin(_conservativeVars);
 173         VecAssemblyEnd(_conservativeVars);
 174         VecCopy(_conservativeVars, _old);
 175         VecAssemblyBegin(_old);
 176         VecAssemblyEnd(_old);
 177         VecSetValuesBlocked(_primitiveVars, _Nmailles, indices, initialFieldPrim, INSERT_VALUES);
 178         VecAssemblyBegin(_primitiveVars);
 179         VecAssemblyEnd(_primitiveVars);
 180         if(_system)
 181         {
 182                 cout << "Variables primitives initiales : " << endl;
 183                 VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_WORLD);
 184                 cout << endl;
 185                 cout<<"Variables conservatives initiales : "<<endl;
 186                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_SELF);
 187         }
 188
 189         delete[] initialFieldPrim;
 190         delete[] initialFieldCons;
 191         delete[] indices;
 192
 193         //Linear solver
 194         KSPCreate(PETSC_COMM_SELF, &_ksp);
 195         KSPSetType(_ksp, _ksptype);
 196         // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
 197         KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
 198         KSPGetPC(_ksp, &_pc);
 199         PCSetType(_pc, _pctype);
 200
 201         _initializedMemory=true;
 202         save();//save initial data
 203 }
 204
 205 bool ProblemFluid::initTimeStep(double dt){
 206         _dt = dt;
 207         return _dt>0;//No need to call MatShift as the linear system matrix is filled at each Newton iteration (unlike linear problem)
 208 }
 209
 210 bool ProblemFluid::iterateTimeStep(bool &converged)
 211 {
 212         bool stop=false;
 213
 214         if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 215                 _maxvp=0.;
 216                 computeTimeStep(stop);//This compute timestep is just to update the linear system. The time step was imposed befor starting the Newton iterations
 217         }
 218         if(stop){//Le compute time step ne s'est pas bien passé
 219                 cout<<"ComputeTimeStep failed"<<endl;
 220                 converged=false;
 221                 return false;
 222         }
 223
 224         computeNewtonVariation();
 225
 226         //converged=convergence des iterations
 227         if(_timeScheme == Explicit)
 228                 converged=true;
 229         else{//Implicit scheme
 230
 231                 KSPGetIterationNumber(_ksp, &_PetscIts);
 232                 if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
 233                         _MaxIterLinearSolver = _PetscIts;
 234
 235                 KSPConvergedReason reason;
 236                 KSPGetConvergedReason(_ksp,&reason);
 237
 238                 if(reason<0)//solving the linear system failed
 239                 {
 240                         cout<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
 241                         cout<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
 242                         *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
 243                         *_runLogFile<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
 244                         converged=false;
 245                         return false;
 246                 }
 247                 else{//solving the linear system succeeded
 248                         //Calcul de la variation relative Uk+1-Uk
 249                         _erreur_rel = 0.;
 250                         double x, dx;
 251                         int I;
 252                         for(int j=1; j<=_Nmailles; j++)
 253                         {
 254                                 for(int k=0; k<_nVar; k++)
 255                                 {
 256                                         I = (j-1)*_nVar + k;
 257                                         VecGetValues(_newtonVariation, 1, &I, &dx);
 258                                         VecGetValues(_conservativeVars, 1, &I, &x);
 259                                         if (fabs(x)*fabs(x)< _precision)
 260                                         {
 261                                                 if(_erreur_rel < fabs(dx))
 262                                                         _erreur_rel = fabs(dx);
 263                                         }
 264                                         else if(_erreur_rel < fabs(dx/x))
 265                                                 _erreur_rel = fabs(dx/x);
 266                                 }
 267                         }
 268                 }
 269                 converged = _erreur_rel <= _precision_Newton;
 270         }
 271
 272         double relaxation=1;//Uk+1=Uk+relaxation*daltaU
 273
 274         VecAXPY(_conservativeVars,  relaxation, _newtonVariation);
 275
 276         //mise a jour du champ primitif
 277         updatePrimitives();
 278
 279         if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
 280         {
 281                 cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 282                 *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 283                 converged=false;
 284                 return false;
 285         }
 286         if(_system)
 287         {
 288                 cout<<"Vecteur Ukp1-Uk "<<endl;
 289                 VecView(_newtonVariation,  PETSC_VIEWER_STDOUT_SELF);
 290                 cout << "Nouvel etat courant Uk de l'iteration Newton: " << endl;
 291                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_SELF);
 292         }
 293
 294         if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
 295                 if(_minm1<-_precision || _minm2<-_precision)
 296                 {
 297                         cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 298                         *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 299                 }
 300
 301                 if (_nbVpCplx>0){
 302                         cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 303                         *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 304                 }
 305         }
 306         _minm1=1e30;
 307         _minm2=1e30;
 308         _nbMaillesNeg=0;
 309         _nbVpCplx =0;
 310         _part_imag_max=0;
 311
 312         return true;
 313 }
 314
 315 double ProblemFluid::computeTimeStep(bool & stop){
 316         VecZeroEntries(_b);
 317
 318         if(_restartWithNewTimeScheme)//This is a change of time scheme during a simulation
 319         {
 320                 if(_timeScheme == Implicit)
 321                         MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
 322                 else
 323                         MatDestroy(&_A);
 324                 _restartWithNewTimeScheme=false;
 325         }
 326         if(_timeScheme == Implicit)
 327                 MatZeroEntries(_A);
 328
 329         VecAssemblyBegin(_b);
 330         VecAssemblyBegin(_conservativeVars);
 331         std::vector< int > idCells;
 332         PetscInt idm, idn, size = 1;
 333
 334         long nbFaces = _mesh.getNumberOfFaces();
 335         Face Fj;
 336         Cell Ctemp1,Ctemp2;
 337         string nameOfGroup;
 338
 339         for (int j=0; j<nbFaces;j++){
 340                 Fj = _mesh.getFace(j);
 341                 _isBoundary=Fj.isBorder();
 342                 idCells = Fj.getCellsId();
 343
 344                 // If Fj is on the boundary
 345                 if (_isBoundary)
 346                 {
 347                         for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
 348                         {
 349                                 // compute the normal vector corresponding to face j : from Ctemp1 outward
 350                                 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
 351                                 if (_Ndim >1){
 352                                         for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
 353                                         {//we look for l the index of the face Fj for the cell Ctemp1
 354                                                 if (j == Ctemp1.getFacesId()[l])
 355                                                 {
 356                                                         for (int idim = 0; idim < _Ndim; ++idim)
 357                                                                 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
 358                                                         break;
 359                                                 }
 360                                         }
 361                                 }else{ // _Ndim = 1, build normal vector (bug cdmath)
 362                                         if(!_sectionFieldSet)
 363                                         {
 364                                                 if (Fj.x()<Ctemp1.x())
 365                                                         _vec_normal[0] = -1;
 366                                                 else
 367                                                         _vec_normal[0] = 1;
 368                                         }
 369                                         else
 370                                         {
 371                                                 if(idCells[0]==0)
 372                                                         _vec_normal[0] = -1;
 373                                                 else//idCells[0]==31
 374                                                         _vec_normal[0] = 1;
 375                                         }
 376                                 }
 377                                 if(_verbose && _nbTimeStep%_freqSave ==0)
 378                                 {
 379                                         cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
 380                                         for(int p=0; p<_Ndim; p++)
 381                                                 cout << _vec_normal[p] << ",";
 382                                         cout << "). "<<endl;
 383                                 }
 384                                 nameOfGroup = Fj.getGroupName();
 385                                 _porosityi=_porosityField(idCells[k]);
 386                                 _porosityj=_porosityi;
 387                                 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
 388                                 convectionState(idCells[k],0,true);
 389                                 convectionMatrices();
 390                                 diffusionStateAndMatrices(idCells[k], 0, true);
 391                                 // compute 1/dxi
 392                                 if (_Ndim > 1)
 393                                         _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
 394                                 else
 395                                         _inv_dxi = 1/Ctemp1.getMeasure();
 396
 397                                 addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
 398                                 addDiffusionToSecondMember(idCells[k],-1,true);
 399                                 addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
 400
 401                                 if(_timeScheme == Implicit){
 402                                         for(int l=0; l<_nVar*_nVar;l++){
 403                                                 _AroeMinusImplicit[l] *= _inv_dxi;
 404                                                 _Diffusion[l] *=_inv_dxi*_inv_dxi;
 405                                         }
 406
 407                                         jacobian(idCells[k],nameOfGroup,_vec_normal);
 408                                         jacobianDiff(idCells[k],nameOfGroup);
 409                                         if(_verbose && _nbTimeStep%_freqSave ==0){
 410                                                 cout << "Matrice Jacobienne CL convection:" << endl;
 411                                                 for(int p=0; p<_nVar; p++){
 412                                                         for(int q=0; q<_nVar; q++)
 413                                                                 cout << _Jcb[p*_nVar+q] << "\t";
 414                                                         cout << endl;
 415                                                 }
 416                                                 cout << endl;
 417                                                 cout << "Matrice Jacobienne CL diffusion:" << endl;
 418                                                 for(int p=0; p<_nVar; p++){
 419                                                         for(int q=0; q<_nVar; q++)
 420                                                                 cout << _JcbDiff[p*_nVar+q] << "\t";
 421                                                         cout << endl;
 422                                                 }
 423                                                 cout << endl;
 424                                         }
 425                                         idm = idCells[k];
 426                                         Polynoms Poly;
 427                                         //calcul et insertion de A^-*Jcb
 428                                         Poly.matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
 429                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
 430
 431                                         if(_verbose)
 432                                                 displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
 433
 434                                         //insertion de -A^-
 435                                         for(int k=0; k<_nVar*_nVar;k++){
 436                                                 _AroeMinusImplicit[k] *= -1;
 437                                         }
 438                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 439                                         if(_verbose)
 440                                                 displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
 441
 442                                         //calcul et insertion de D*JcbDiff
 443                                         Poly.matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
 444                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
 445                                         for(int k=0; k<_nVar*_nVar;k++)
 446                                                 _Diffusion[k] *= -1;
 447                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 448                                 }
 449                         }
 450                 } else  if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
 451                         // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
 452                         Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
 453                         Ctemp2 = _mesh.getCell(idCells[1]);
 454                         if (_Ndim >1){
 455                                 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
 456                                         if (j == Ctemp1.getFacesId()[l]){
 457                                                 for (int idim = 0; idim < _Ndim; ++idim)
 458                                                         _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
 459                                                 break;
 460                                         }
 461                                 }
 462                         }else{ // _Ndim = 1, build normal vector (bug cdmath)
 463                                 if(!_sectionFieldSet)
 464                                 {
 465                                         if (Fj.x()<Ctemp1.x())
 466                                                 _vec_normal[0] = -1;
 467                                         else
 468                                                 _vec_normal[0] = 1;
 469                                 }
 470                                 else
 471                                 {
 472                                         if(idCells[0]>idCells[1])
 473                                                 _vec_normal[0] = -1;
 474                                         else
 475                                                 _vec_normal[0] = 1;
 476                                 }
 477                         }
 478                         if(_verbose && _nbTimeStep%_freqSave ==0)
 479                         {
 480                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 481                                 cout<<" Normal vector= ";
 482                                 for (int idim = 0; idim < _Ndim; ++idim)
 483                                         cout<<_vec_normal[idim]<<", ";
 484                                 cout<<endl;
 485                         }
 486                         _porosityi=_porosityField(idCells[0]);
 487                         _porosityj=_porosityField(idCells[1]);
 488                         convectionState(idCells[0],idCells[1],false);
 489                         convectionMatrices();
 490                         diffusionStateAndMatrices(idCells[0], idCells[1], false);
 491
 492                         // compute 1/dxi and 1/dxj
 493                         if (_Ndim > 1)
 494                         {
 495                                 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
 496                                 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
 497                         }
 498                         else
 499                         {
 500                                 _inv_dxi = 1/Ctemp1.getMeasure();
 501                                 _inv_dxj = 1/Ctemp2.getMeasure();
 502                         }
 503
 504                         addConvectionToSecondMember( idCells[0],idCells[1], false);
 505                         addDiffusionToSecondMember( idCells[0],idCells[1], false);
 506                         addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
 507
 508                         if(_timeScheme == Implicit){
 509                                 for(int k=0; k<_nVar*_nVar;k++)
 510                                 {
 511                                         _AroeMinusImplicit[k] *= _inv_dxi;
 512                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
 513                                 }
 514                                 idm = idCells[0];
 515                                 idn = idCells[1];
 516                                 //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNb<<endl;
 517                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
 518                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
 519
 520                                 if(_verbose){
 521                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
 522                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
 523                                 }
 524                                 for(int k=0;k<_nVar*_nVar;k++){
 525                                         _AroeMinusImplicit[k] *= -1;
 526                                         _Diffusion[k] *= -1;
 527                                 }
 528                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 529                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 530                                 if(_verbose){
 531                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
 532                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
 533                                 }
 534                                 for(int k=0; k<_nVar*_nVar;k++)
 535                                 {
 536                                         _AroePlusImplicit[k]  *= _inv_dxj;
 537                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;
 538                                 }
 539                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
 540                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
 541                                 if(_verbose)
 542                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
 543
 544                                 for(int k=0;k<_nVar*_nVar;k++){
 545                                         _AroePlusImplicit[k] *= -1;
 546                                         _Diffusion[k] *= -1;
 547                                 }
 548                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
 549                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
 550
 551                                 if(_verbose)
 552                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
 553                         }
 554                 }
 555                 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
 556                         if(_verbose && _nbTimeStep%_freqSave ==0)
 557                                 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
 558                         *_runLogFile<<"Warning: treatment of a junction node"<<endl;
 559
 560                         if(!_sectionFieldSet)
 561                         {
 562                                 _runLogFile->close();
 563                                 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
 564                         }
 565                         int largestSectionCellIndex=0;
 566                         for(int i=1;i<Fj.getNumberOfCells();i++){
 567                                 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
 568                                         largestSectionCellIndex=i;
 569                         }
 570                         idm = idCells[largestSectionCellIndex];
 571                         Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
 572                         _porosityi=_porosityField(idm);
 573
 574                         if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
 575                                 _vec_normal[0] = 1;
 576                         else//j==16
 577                                 _vec_normal[0] = -1;
 578                         if(_verbose && _nbTimeStep%_freqSave ==0)
 579                         {
 580                                 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
 581                                 cout << " ; vecteur normal=(";
 582                                 for(int p=0; p<_Ndim; p++)
 583                                         cout << _vec_normal[p] << ",";
 584                                 cout << "). "<<endl;
 585                         }
 586                         for(int i=0;i<Fj.getNumberOfCells();i++){
 587                                 if(i != largestSectionCellIndex){
 588                                         idn = idCells[i];
 589                                         Ctemp2 = _mesh.getCell(idn);
 590                                         _porosityj=_porosityField(idn);
 591                                         convectionState(idm,idn,false);
 592                                         convectionMatrices();
 593                                         diffusionStateAndMatrices(idm, idn,false);
 594
 595                                         if(_verbose && _nbTimeStep%_freqSave ==0)
 596                                                 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
 597
 598                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
 599                                         _inv_dxj = 1/Ctemp2.getMeasure();
 600
 601                                         addConvectionToSecondMember(idm,idn, false);
 602                                         _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
 603                                         addDiffusionToSecondMember(idm,idn, false);
 604                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
 605                                         addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
 606
 607                                         if(_timeScheme == Implicit){
 608                                                 for(int k=0; k<_nVar*_nVar;k++)
 609                                                 {
 610                                                         _AroeMinusImplicit[k] *= _inv_dxi;
 611                                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
 612                                                 }
 613                                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
 614                                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
 615
 616                                                 if(_verbose){
 617                                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
 618                                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
 619                                                 }
 620                                                 for(int k=0;k<_nVar*_nVar;k++){
 621                                                         _AroeMinusImplicit[k] *= -1;
 622                                                         _Diffusion[k] *= -1;
 623                                                 }
 624                                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 625                                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 626                                                 if(_verbose){
 627                                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
 628                                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
 629                                                 }
 630                                                 for(int k=0; k<_nVar*_nVar;k++)
 631                                                 {
 632                                                         _AroePlusImplicit[k] *= _inv_dxj;
 633                                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
 634                                                 }
 635                                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
 636                                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
 637                                                 if(_verbose)
 638                                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
 639
 640                                                 for(int k=0;k<_nVar*_nVar;k++){
 641                                                         _AroePlusImplicit[k] *= -1;
 642                                                         _Diffusion[k] *= -1;
 643                                                 }
 644                                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
 645                                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
 646
 647                                                 if(_verbose)
 648                                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
 649                                         }
 650                                 }
 651                         }
 652                 }
 653                 else
 654                 {
 655                         cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
 656                         _runLogFile->close();
 657                         throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
 658                 }
 659
 660         }
 661         VecAssemblyEnd(_conservativeVars);
 662         VecAssemblyEnd(_b);
 663
 664         if(_timeScheme == Implicit){
 665                 for(int imaille = 0; imaille<_Nmailles; imaille++)
 666                         MatSetValuesBlocked(_A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
 667
 668                 if(_verbose && _nbTimeStep%_freqSave ==0)
 669                         displayMatrix(_GravityImplicitationMatrix,_nVar,"Gravity matrix:");
 670
 671                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 672                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 673                 if(_verbose && _nbTimeStep%_freqSave ==0){
 674                         cout << "Fin ProblemFluid.cxx : matrice implicite"<<endl;
 675                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 676                         cout << endl;
 677                 }
 678         }
 679
 680         stop=false;
 681
 682
 683         /*
 684         if(_nbTimeStep+1<_cfl)
 685                 return (_nbTimeStep+1)*_minl/_maxvp;
 686         else
 687          */
 688         return _cfl*_minl/_maxvp;
 689 }
 690
 691 void ProblemFluid::computeNewtonVariation()
 692 {
 693         if(_verbose)
 694         {
 695                 cout<<"Vecteur courant Uk "<<endl;
 696                 VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
 697                 cout << endl;
 698         }
 699         if(_timeScheme == Explicit)
 700         {
 701                 VecCopy(_b,_newtonVariation);
 702                 VecScale(_newtonVariation, _dt);
 703                 if(_verbose && _nbTimeStep%_freqSave ==0)
 704                 {
 705                         cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
 706                         VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
 707                         cout << endl;
 708                 }
 709         }
 710         else
 711         {
 712                 if(_verbose)
 713                 {
 714                         cout << "Matrice du système linéaire avant contribution delta t" << endl;
 715                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 716                         cout << endl;
 717                         cout << "Second membre du système linéaire avant contribution delta t" << endl;
 718                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 719                         cout << endl;
 720                 }
 721                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 722                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 723
 724                 VecAXPY(_b, 1/_dt, _old);
 725                 VecAXPY(_b, -1/_dt, _conservativeVars);
 726                 MatShift(_A, 1/_dt);
 727
 728 #if PETSC_VERSION_GREATER_3_5
 729                 KSPSetOperators(_ksp, _A, _A);
 730 #else
 731                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 732 #endif
 733
 734                 if(_verbose)
 735                 {
 736                         cout << "Matrice du système linéaire après contribution delta t" << endl;
 737                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 738                         cout << endl;
 739                         cout << "Second membre du système linéaire après contribution delta t" << endl;
 740                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 741                         cout << endl;
 742                 }
 743
 744                 if(_conditionNumber)
 745                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 746                 if(!_isScaling)
 747                 {
 748                         KSPSolve(_ksp, _b, _newtonVariation);
 749                 }
 750                 else
 751                 {
 752                         computeScaling(_maxvp);
 753                         int indice;
 754                         VecAssemblyBegin(_vecScaling);
 755                         VecAssemblyBegin(_invVecScaling);
 756                         for(int imaille = 0; imaille<_Nmailles; imaille++)
 757                         {
 758                                 indice = imaille;
 759                                 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
 760                                 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
 761                         }
 762                         VecAssemblyEnd(_vecScaling);
 763                         VecAssemblyEnd(_invVecScaling);
 764
 765                         if(_system)
 766                         {
 767                                 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
 768                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 769                                 cout << endl;
 770                                 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
 771                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 772                                 cout << endl;
 773                         }
 774                         MatDiagonalScale(_A,_vecScaling,_invVecScaling);
 775                         if(_system)
 776                         {
 777                                 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
 778                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 779                                 cout << endl;
 780                         }
 781                         VecCopy(_b,_bScaling);
 782                         VecPointwiseMult(_b,_vecScaling,_bScaling);
 783                         if(_system)
 784                         {
 785                                 cout << "Produit du second membre par le preconditionneur bloc diagonal  a gauche" << endl;
 786                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 787                                 cout << endl;
 788                         }
 789
 790                         KSPSolve(_ksp,_b, _bScaling);
 791                         VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
 792                 }
 793                 if(_verbose)
 794                 {
 795                         cout << "solution du systeme lineaire local:" << endl;
 796                         VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
 797                         cout << endl;
 798                 }
 799         }
 800 }
 801
 802 void ProblemFluid::validateTimeStep()
 803 {
 804         if(_system)
 805         {
 806                 cout<<" Vecteur Un"<<endl;
 807                 VecView(_old,  PETSC_VIEWER_STDOUT_WORLD);
 808                 cout<<" Vecteur Un+1"<<endl;
 809                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
 810         }
 811         VecAXPY(_old,  -1, _conservativeVars);//old contient old-courant
 812
 813         //Calcul de la variation Un+1-Un
 814         _erreur_rel= 0;
 815         double x, dx;
 816         int I;
 817
 818         for(int j=1; j<=_Nmailles; j++)
 819         {
 820                 for(int k=0; k<_nVar; k++)
 821                 {
 822                         I = (j-1)*_nVar + k;
 823                         VecGetValues(_old, 1, &I, &dx);
 824                         VecGetValues(_conservativeVars, 1, &I, &x);
 825                         if (fabs(x)< _precision)
 826                         {
 827                                 if(_erreur_rel < fabs(dx))
 828                                         _erreur_rel = fabs(dx);
 829                         }
 830                         else if(_erreur_rel < fabs(dx/x))
 831                                 _erreur_rel = fabs(dx/x);
 832                 }
 833         }
 834
 835         _isStationary =_erreur_rel <_precision;
 836
 837         VecCopy(_conservativeVars, _old);
 838
 839         if(_verbose && _nbTimeStep%_freqSave ==0){
 840                 if(!_usePrimitiveVarsInNewton)
 841                         testConservation();
 842                 cout <<"Valeur propre locale max: " << _maxvp << endl;
 843         }
 844
 845         if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
 846                 //Find minimum and maximum void fractions
 847                 double alphamin=1e30;
 848                 double alphamax=-1e30;
 849                 double T, Tmax=-1e30;
 850                 int J=_nVar-1;
 851                 I = 0;
 852                 for(int j=0; j<_Nmailles; j++)
 853                 {
 854                         VecGetValues(_primitiveVars, 1, &I, &x);//extract void fraction
 855                         if(x>alphamax)
 856                                 alphamax=x;
 857                         if(x<alphamin)
 858                                 alphamin=x;
 859                         I += _nVar;
 860                         VecGetValues(_primitiveVars, 1, &J, &T);//extract void fraction
 861                         if(T>Tmax)
 862                                 Tmax=T;
 863                         J += _nVar;
 864                 }
 865                 cout<<"Alpha min = " << alphamin << ", Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
 866                 cout<<endl;
 867                 *_runLogFile<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
 868         }
 869
 870         _time+=_dt;
 871         _nbTimeStep++;
 872         if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 873                 save();
 874 }
 875
 876 void ProblemFluid::abortTimeStep(){
 877         _dt = 0;
 878 }
 879
 880 void ProblemFluid::addConvectionToSecondMember
 881 (               const int &i,
 882                 const int &j, bool isBord, string groupname
 883 )
 884 {
 885         if(_verbose && _nbTimeStep%_freqSave ==0)
 886                 cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
 887
 888         //extraction des valeurs
 889         for(int k=0; k<_nVar; k++)
 890                 _idm[k] = _nVar*i + k;
 891         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 892         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
 893
 894         if(!isBord){
 895                 for(int k=0; k<_nVar; k++)
 896                         _idn[k] = _nVar*j + k;
 897                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
 898                 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
 899         }
 900         else{
 901                 for(int k=0; k<_nVar; k++)
 902                         _idn[k] = k;
 903                 VecGetValues(_Uext, _nVar, _idn, _Uj);
 904                 consToPrim(_Uj, _Vj,_porosityj);
 905         }
 906         _idm[0] = i;
 907         _idn[0] = j;
 908
 909         if(_verbose && _nbTimeStep%_freqSave ==0)
 910         {
 911                 cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
 912                 for(int q=0; q<_nVar; q++)
 913                         cout << _Ui[q] << endl;
 914                 cout << endl;
 915                 cout << "addConvectionToSecondMember : état j= " << j << ", _Uj=" << endl;
 916                 for(int q=0; q<_nVar; q++)
 917                         cout << _Uj[q] << endl;
 918                 cout << endl;
 919         }
 920         if(_nonLinearFormulation!=reducedRoe){//VFRoe, Roe or VFFC
 921                 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar),Fij(_nVar);
 922                 for(int i1=0;i1<_nVar;i1++)
 923                 {
 924                         Ui(i1)=_Ui[i1];
 925                         Uj(i1)=_Uj[i1];
 926                         Vi(i1)=_Vi[i1];
 927                         Vj(i1)=_Vj[i1];
 928                 }
 929                 Vector normale(_Ndim);
 930                 for(int i1=0;i1<_Ndim;i1++)
 931                         normale(i1)=_vec_normal[i1];
 932
 933                 Matrix signAroe(_nVar);
 934                 for(int i1=0;i1<_nVar;i1++)
 935                         for(int i2=0;i2<_nVar;i2++)
 936                                 signAroe(i1,i2)=_signAroe[i1*_nVar+i2];
 937
 938                 Matrix absAroe(_nVar);
 939                 for(int i1=0;i1<_nVar;i1++)
 940                         for(int i2=0;i2<_nVar;i2++)
 941                                 absAroe(i1,i2)=_absAroe[i1*_nVar+i2];
 942
 943                 if(_nonLinearFormulation==VFRoe)//VFRoe
 944                 {
 945                         Vector Uij(_nVar),Vij(_nVar);
 946                         double porosityij=sqrt(_porosityi*_porosityj);
 947
 948                         Uij=(Ui+Uj)/2+signAroe*(Ui-Uj)/2;
 949
 950                         for(int i1=0;i1<_nVar;i1++)
 951                                 _Uij[i1]=Uij(i1);
 952
 953                         consToPrim(_Uij, _Vij,porosityij);
 954
 955                         applyVFRoeLowMachCorrections(isBord, groupname);
 956
 957                         for(int i1=0;i1<_nVar;i1++)
 958                         {
 959                                 Uij(i1)=_Uij[i1];
 960                                 Vij(i1)=_Vij[i1];
 961                         }
 962
 963                         Fij=convectionFlux(Uij,Vij,normale,porosityij);
 964
 965                         if(_verbose && _nbTimeStep%_freqSave ==0)
 966                         {
 967                                 cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
 968                                 cout<<Uij<<endl;
 969                                 cout<<"Etat interfacial primitif "<<i<<", "<<j<< endl;
 970                                 cout<<Vij<<endl;
 971                         }
 972                 }
 973                 else //Roe or VFFC
 974                 {
 975                         if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
 976                                 Fij=staggeredVFFCFlux();
 977                         else
 978                         {
 979                                 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
 980                                 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
 981                                 if(_nonLinearFormulation==VFFC)//VFFC
 982                                         Fij=(Fi+Fj)/2+signAroe*(Fi-Fj)/2;
 983                                 else if(_nonLinearFormulation==Roe)//Roe
 984                                         Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
 985
 986                                 if(_verbose && _nbTimeStep%_freqSave ==0)
 987                                 {
 988                                         cout<<"Flux cellule "<<i<<" = "<< endl;
 989                                         cout<<Fi<<endl;
 990                                         cout<<"Flux cellule "<<j<<" = "<< endl;
 991                                         cout<<Fj<<endl;
 992                                 }
 993                         }
 994                 }
 995                 for(int i1=0;i1<_nVar;i1++)
 996                         _phi[i1]=-Fij(i1)*_inv_dxi;
 997                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
 998                 if(_verbose && _nbTimeStep%_freqSave ==0)
 999                 {
1000                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1001                         cout<<"Flux interfacial "<<i<<", "<<j<< endl;
1002                         cout<<Fij<<endl;
1003                         cout << "Contribution convection à " << i << ", -Fij(i1)*_inv_dxi= "<<endl;
1004                         for(int q=0; q<_nVar; q++)
1005                                 cout << _phi[q] << endl;
1006                         cout << endl;
1007                 }
1008                 if(!isBord){
1009                         for(int i1=0;i1<_nVar;i1++)
1010                                 _phi[i1]*=-_inv_dxj/_inv_dxi;
1011                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1012                         if(_verbose && _nbTimeStep%_freqSave ==0)
1013                         {
1014                                 cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
1015                                 for(int q=0; q<_nVar; q++)
1016                                         cout << _phi[q] << endl;
1017                                 cout << endl;
1018                         }
1019                 }
1020         }else //_nonLinearFormulation==reducedRoe)
1021         {
1022                 for(int k=0; k<_nVar; k++)
1023                         _temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
1024                 Polynoms Poly;
1025                 Poly.matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
1026                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1027
1028                 if(_verbose && _nbTimeStep%_freqSave ==0)
1029                 {
1030                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1031                         cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
1032                         for(int q=0; q<_nVar; q++)
1033                                 cout << _temp[q] << endl;
1034                         cout << endl;
1035                         cout << "Contribution convection à " << i << ", A^-*(Ui - Uj)*_inv_dxi= "<<endl;
1036                         for(int q=0; q<_nVar; q++)
1037                                 cout << _phi[q] << endl;
1038                         cout << endl;
1039                 }
1040
1041                 if(!isBord)
1042                 {
1043                         for(int k=0; k<_nVar; k++)
1044                                 _temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
1045                         Poly.matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
1046                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1047
1048                         if(_verbose && _nbTimeStep%_freqSave ==0)
1049                         {
1050                                 cout << "Contribution convection à  " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
1051                                 for(int q=0; q<_nVar; q++)
1052                                         cout << _phi[q] << endl;
1053                                 cout << endl;
1054                         }
1055                 }
1056         }
1057         if(_verbose && _nbTimeStep%_freqSave ==0)
1058         {
1059                 cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
1060                 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
1061                 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
1062                 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
1063                 displayMatrix(_signAroe, _nVar,"Signe de la matrice de Roe");
1064         }
1065 }
1066
1067 void ProblemFluid::addSourceTermToSecondMember
1068 (       const int i, int nbVoisinsi,
1069                 const int j, int nbVoisinsj,
1070                 bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
1071 {
1072         if(_verbose && _nbTimeStep%_freqSave ==0)
1073                 cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
1074
1075         _idm[0] = i*_nVar;
1076         for(int k=1; k<_nVar;k++)
1077                 _idm[k] = _idm[k-1] + 1;
1078         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1079         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1080         sourceVector(_Si,_Ui,_Vi,i);
1081
1082         if (_verbose && _nbTimeStep%_freqSave ==0)
1083         {
1084                 cout << "Terme source S(Ui), i= " << i<<endl;
1085                 for(int q=0; q<_nVar; q++)
1086                         cout << _Si[q] << endl;
1087                 cout << endl;
1088         }
1089         if(!isBord){
1090                 for(int k=0; k<_nVar; k++)
1091                         _idn[k] = _nVar*j + k;
1092                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1093                 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1094                 sourceVector(_Sj,_Uj,_Vj,j);
1095         }else
1096         {
1097                 for(int k=0; k<_nVar; k++)
1098                         _idn[k] = k;
1099                 VecGetValues(_Uext, _nVar, _idn, _Uj);
1100                 consToPrim(_Uj, _Vj,_porosityj);
1101                 sourceVector(_Sj,_Uj,_Vj,i);
1102         }
1103         if (_verbose && _nbTimeStep%_freqSave ==0)
1104         {
1105                 if(!isBord)
1106                         cout << "Terme source S(Uj), j= " << j<<endl;
1107                 else
1108                         cout << "Terme source S(Uj), cellule fantôme "<<endl;
1109                 for(int q=0; q<_nVar; q++)
1110                         cout << _Sj[q] << endl;
1111                 cout << endl;
1112         }
1113
1114         //Compute pressure loss vector
1115         double K;
1116         if(_pressureLossFieldSet){
1117                 K=_pressureLossField(ij);
1118                 pressureLossVector(_pressureLossVector, K,_Ui, _Vi, _Uj, _Vj);
1119         }
1120         else{
1121                 K=0;
1122                 for(int k=0; k<_nVar;k++)
1123                         _pressureLossVector[k]=0;
1124         }
1125         if (_verbose && _nbTimeStep%_freqSave ==0)
1126         {
1127                 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
1128                 for(int k=0; k<_nVar;k++)
1129                         cout<< _pressureLossVector[k]<<", ";
1130                 cout<<endl;
1131         }
1132         //Contribution of the porosityField gradient:
1133         if(!isBord)
1134                 porosityGradientSourceVector();
1135         else{
1136                 for(int k=0; k<_nVar;k++)
1137                         _porosityGradientSourceVector[k]=0;
1138         }
1139
1140         if (_verbose && _nbTimeStep%_freqSave ==0)
1141         {
1142                 if(!isBord)
1143                         cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
1144                 else
1145                         cout<<"interface frontière i= "<<i<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<endl;
1146                 cout<<"Gradient de porosite à l'interface"<<endl;
1147                 for(int k=0; k<_nVar;k++)
1148                         cout<< _porosityGradientSourceVector[k]<<", ";
1149                 cout<<endl;
1150         }
1151         Polynoms Poly;
1152         if(!isBord){
1153                 if(_wellBalancedCorrection){
1154                         for(int k=0; k<_nVar;k++)
1155                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1156                         Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1157                         for(int k=0; k<_nVar;k++){
1158                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1159                                 _Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
1160                         }
1161                         if (_verbose && _nbTimeStep%_freqSave ==0)
1162                         {
1163                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,j), j="<<j<<endl;
1164                                 for(int q=0; q<_nVar; q++)
1165                                         cout << _Si[q] << endl;
1166                                 cout << "Contribution au terme source Sj de la cellule j= " << j<<" venant  (après décentrement) de la face (i,j), i="<<i<<endl;
1167                                 for(int q=0; q<_nVar; q++)
1168                                         cout << _Sj[q] << endl;
1169                                 cout << endl;
1170                                 cout<<"ratio surface sur volume i = "<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1171                                 cout<<"ratio surface sur volume j = "<<mesureFace/_perimeters(j)<<" perimeter = "<< _perimeters(j)<<endl;
1172                                 cout << endl;
1173                         }
1174                 }
1175                 else{
1176                         for(int k=0; k<_nVar;k++){
1177                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
1178                                 _Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
1179                         }
1180                         if (_verbose && _nbTimeStep%_freqSave ==0)
1181                         {
1182                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant  de la face (i,j), j="<<j<<endl;
1183                                 for(int q=0; q<_nVar; q++)
1184                                         cout << _Si[q] << endl;
1185                                 cout << "Contribution au terme source Sj de la cellule j = " << j<<" venant  de la face (i,j), i="<<i <<endl;
1186                                 for(int q=0; q<_nVar; q++)
1187                                         cout << _Sj[q] << endl;
1188                                 cout << endl;
1189                         }
1190                 }
1191                 _idn[0] = j;
1192                 VecSetValuesBlocked(_b, 1, _idn, _Sj, ADD_VALUES);
1193         }else{
1194                 if(_wellBalancedCorrection){
1195                         for(int k=0; k<_nVar;k++)
1196                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1197                         Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1198                         for(int k=0; k<_nVar;k++)
1199                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1200                         if (_verbose && _nbTimeStep%_freqSave ==0)
1201                         {
1202                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,bord)"<<endl;
1203                                 for(int q=0; q<_nVar; q++)
1204                                         cout << _Si[q] << endl;
1205                                 cout<<"ratio surface sur volume i ="<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1206                                 cout << endl;
1207                         }
1208                 }
1209                 else
1210                 {
1211                         for(int k=0; k<_nVar;k++)
1212                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
1213                         if (_verbose && _nbTimeStep%_freqSave ==0)
1214                         {
1215                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
1216                                 for(int q=0; q<_nVar; q++)
1217                                         cout << _Si[q] << endl;
1218                                 cout << endl;
1219                         }
1220                 }
1221         }
1222         _idm[0] = i;
1223         VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
1224
1225         if(_verbose && _nbTimeStep%_freqSave ==0 && _wellBalancedCorrection)
1226                 displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
1227 }
1228
1229 void ProblemFluid::updatePrimitives()
1230 {
1231         VecAssemblyBegin(_primitiveVars);
1232         for(int i=1; i<=_Nmailles; i++)
1233         {
1234                 _idm[0] = _nVar*( (i-1));
1235                 for(int k=1; k<_nVar; k++)
1236                         _idm[k] = _idm[k-1] + 1;
1237
1238                 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1239                 if(_verbose && _nbTimeStep%_freqSave ==0)
1240                 {
1241                         cout << "ProblemFluid::updatePrimitives() cell " << i-1 << endl;
1242                         cout << "Ui = ";
1243                         for(int q=0; q<_nVar; q++)
1244                                 cout << _Ui[q]  << "\t";
1245                         cout << endl;
1246                         cout << "Vi = ";
1247                         for(int q=0; q<_nVar; q++)
1248                                 cout << _Vi[q]  << "\t";
1249                         cout << endl;
1250                 }
1251
1252                 consToPrim(_Ui,_Vi,_porosityField(i-1));
1253                 if(_nbPhases==2 && _Psat>-1e30){//Cas simulation flashing
1254                         double pressure;
1255                         VecGetValues(_primitiveVars, 1, _idm+1, &pressure);
1256                         _dp_over_dt(i-1)=(_Vi[1]-pressure)/_dt;//pn+1-pn
1257                 }
1258                 _idm[0] = i-1;
1259                 VecSetValuesBlocked(_primitiveVars, 1, _idm, _Vi, INSERT_VALUES);
1260         }
1261         VecAssemblyEnd(_primitiveVars);
1262
1263         if(_system)
1264         {
1265                 cout << "Nouvelles variables primitives : " << endl;
1266                 VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_WORLD);
1267                 cout << endl;
1268         }
1269 }
1270
1271 void ProblemFluid::updateConservatives()
1272 {
1273         VecAssemblyBegin(_conservativeVars);
1274         for(int i=1; i<=_Nmailles; i++)
1275         {
1276                 _idm[0] = _nVar*( (i-1));
1277                 for(int k=1; k<_nVar; k++)
1278                         _idm[k] = _idm[k-1] + 1;
1279
1280                 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1281
1282                 primToCons(_Vi,0,_Ui,0);
1283                 _idm[0] = i-1;
1284                 VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
1285
1286                 if(_verbose && _nbTimeStep%_freqSave ==0)
1287                 {
1288                         cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
1289                         cout << "Vi = ";
1290                         for(int q=0; q<_nVar; q++)
1291                                 cout << _Vi[q]  << "\t";
1292                         cout << endl;
1293                         cout << "Ui = ";
1294                         for(int q=0; q<_nVar; q++)
1295                                 cout << _Ui[q]  << "\t";
1296                         cout << endl;
1297                 }
1298         }
1299         VecAssemblyEnd(_conservativeVars);
1300
1301         if(_system)
1302         {
1303                 cout << "Nouvelles variables conservatives : " << endl;
1304                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
1305                 cout << endl;
1306         }
1307 }
1308
1309 vector< complex<double> > ProblemFluid::getRacines(vector< double > pol_car){
1310         int degre_polynom = pol_car.size()-1;
1311         double tmp;
1312         vector< complex< double > > valeurs_propres (degre_polynom);
1313         vector< double > zeror(degre_polynom);
1314         vector< double > zeroi(degre_polynom);
1315         for(int j=0; j<(degre_polynom+1)/2; j++){//coefficients in order of decreasing powers for rpoly
1316                 tmp=pol_car[j];
1317                 pol_car[j] =pol_car[degre_polynom-j];
1318                 pol_car[degre_polynom-j]=tmp;
1319         }
1320
1321         //Calcul des racines du polynome
1322         roots_polynoms roots;
1323         int size_vp = roots.rpoly(&pol_car[0],degre_polynom,&zeror[0],&zeroi[0]);
1324
1325         //On ordonne les valeurs propres
1326         if(zeror[1]<zeror[0])
1327         {
1328                 tmp=zeror[0];
1329                 zeror[0]=zeror[1];
1330                 zeror[1]=tmp;
1331                 tmp=zeroi[0];
1332                 zeroi[0]=zeroi[1];
1333                 zeroi[1]=tmp;
1334         }
1335
1336         if(size_vp == degre_polynom) {
1337                 for(int ct =0; ct<degre_polynom; ct++)
1338                         valeurs_propres[ct] = complex<double> (zeror[ct],zeroi[ct]);  //vp non triviales
1339         }
1340         else {
1341                 cout << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1342                 *_runLogFile << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1343                 cout<<"Coefficients polynome caracteristique: "<<endl;
1344                 for(int ct =0; ct<degre_polynom+1; ct++)
1345                         cout<<pol_car[ct]<< " , " <<endl;
1346
1347                 *_runLogFile<<"getRacines computation failed"<<endl;
1348                 _runLogFile->close();
1349                 throw CdmathException("getRacines computation failed");
1350         }
1351
1352         return valeurs_propres;
1353 }
1354
1355 void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1356 {
1357         Polynoms Poly;
1358         int nbVp_dist=valeurs_propres_dist.size();
1359         vector< complex< double > > y (nbVp_dist,0);
1360         for( int i=0 ; i<nbVp_dist ; i++)
1361                 y[i] = Poly.abs_generalise(valeurs_propres_dist[i]);
1362         Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
1363         if(_verbose && _nbTimeStep%_freqSave ==0)
1364         {
1365                 cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1366                 for(int ct =0; ct<nbVp_dist; ct++)
1367                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1368                 cout<< endl;
1369                 cout<<"ProblemFluid::AbsMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1370                 for(int ct =0; ct<nbVp_dist; ct++)
1371                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1372                 cout<< endl;
1373                 displayMatrix(_absAroe,_nVar,"Valeur absolue de la matrice de Roe");
1374         }
1375 }
1376
1377 void ProblemFluid::SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1378 {
1379         int nbVp_dist=valeurs_propres_dist.size();
1380         vector< complex< double > > y (nbVp_dist,0);
1381         for( int i=0 ; i<nbVp_dist ; i++)
1382         {
1383                 if(valeurs_propres_dist[i].real()>0)
1384                         y[i] = 1;
1385                 else if(valeurs_propres_dist[i].real()<0)
1386                         y[i] = -1;
1387                 else
1388                         y[i] = 0;
1389         }
1390         Polynoms Poly;
1391         Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
1392         if(_verbose && _nbTimeStep%_freqSave ==0)
1393         {
1394                 cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1395                 for(int ct =0; ct<nbVp_dist; ct++)
1396                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1397                 cout<< endl;
1398                 cout<<"ProblemFluid::SigneMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1399                 for(int ct =0; ct<nbVp_dist; ct++)
1400                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1401                 cout<< endl;
1402                 displayMatrix(_signAroe,_nVar,"Signe matrice de Roe");
1403         }
1404 }
1405 void ProblemFluid::InvMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1406 {
1407         int nbVp_dist=valeurs_propres_dist.size();
1408         vector< complex< double > > y (nbVp_dist,0);
1409         Polynoms Poly;
1410         for( int i=0 ; i<nbVp_dist ; i++)
1411         {
1412                 if(Poly.module(valeurs_propres_dist[i])>_precision)
1413                         y[i] = 1./valeurs_propres_dist[i];
1414                 else
1415                         y[i] = 1./_precision;
1416         }
1417         Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
1418         if(_verbose && _nbTimeStep%_freqSave ==0)
1419         {
1420                 cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
1421                 for(int ct =0; ct<nbVp_dist; ct++)
1422                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1423                 cout<< endl;
1424                 cout<<"ProblemFluid::InvMatriceRoe : Valeurs à atteindre pour l'interpolation"<<endl;
1425                 for(int ct =0; ct<nbVp_dist; ct++)
1426                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1427                 cout<< endl;
1428                 displayMatrix(_invAroe,_nVar,"Inverse matrice de Roe");
1429         }
1430 }
1431
1432 Field ProblemFluid::getConservativeField() const
1433 {
1434         if(!_initializedMemory)
1435         {
1436                 _runLogFile->close();
1437                 throw CdmathException("ProblemFluid::getConservativeField call initialize() first");
1438         }
1439         return _UU;
1440 }
1441
1442 Field ProblemFluid::getPrimitiveField() const
1443 {
1444         if(!_initializedMemory)
1445         {
1446                 _runLogFile->close();
1447                 throw CdmathException("ProblemFluid::getPrimitiveField call initialize() first");
1448         }
1449         return _VV;
1450 }
1451
1452 void ProblemFluid::terminate(){
1453
1454         delete[] _AroePlus;
1455         delete[] _Diffusion;
1456         delete[] _GravityImplicitationMatrix;
1457         delete[] _AroeMinus;
1458         delete[] _Aroe;
1459         delete[] _absAroe;
1460         delete[] _signAroe;
1461         delete[] _invAroe;
1462         delete[] _AroeImplicit;
1463         delete[] _AroeMinusImplicit;
1464         delete[] _AroePlusImplicit;
1465         delete[] _absAroeImplicit;
1466         delete[] _phi;
1467         delete[] _Jcb;
1468         delete[] _JcbDiff;
1469         delete[] _a;
1470         delete[] _primToConsJacoMat;
1471
1472         delete[] _l;
1473         delete[] _r;
1474         delete[] _Uroe;
1475         delete[] _Udiff;
1476         delete[] _temp;
1477         delete[] _externalStates;
1478         delete[] _idm;
1479         delete[] _idn;
1480         delete[] _vec_normal;
1481         delete[] _Ui;
1482         delete[] _Uj;
1483         delete[] _Vi;
1484         delete[] _Vj;
1485         if(_nonLinearFormulation==VFRoe){
1486                 delete[] _Uij;
1487                 delete[] _Vij;
1488         }
1489         delete[] _Si;
1490         delete[] _Sj;
1491         delete[] _pressureLossVector;
1492         delete[] _porosityGradientSourceVector;
1493         if(_isScaling)
1494         {
1495                 delete[] _blockDiag;
1496                 delete[] _invBlockDiag;
1497
1498                 VecDestroy(&_vecScaling);
1499                 VecDestroy(&_invVecScaling);
1500                 VecDestroy(&_bScaling);
1501         }
1502
1503         VecDestroy(&_conservativeVars);
1504         VecDestroy(&_newtonVariation);
1505         VecDestroy(&_b);
1506         VecDestroy(&_primitiveVars);
1507         VecDestroy(&_Uext);
1508         VecDestroy(&_Vext);
1509         VecDestroy(&_Uextdiff);
1510
1511         //      PCDestroy(_pc);
1512         KSPDestroy(&_ksp);
1513         for(int i=0;i<_nbPhases;i++)
1514                 delete _fluides[i];
1515 }