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*deltaU
 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                                         //calcul et insertion de A^-*Jcb
 427                                         Polynoms::matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
 428                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
 429
 430                                         if(_verbose)
 431                                                 displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
 432
 433                                         //insertion de -A^-
 434                                         for(int k=0; k<_nVar*_nVar;k++){
 435                                                 _AroeMinusImplicit[k] *= -1;
 436                                         }
 437                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 438                                         if(_verbose)
 439                                                 displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
 440
 441                                         //calcul et insertion de D*JcbDiff
 442                                         Polynoms::matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
 443                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
 444                                         for(int k=0; k<_nVar*_nVar;k++)
 445                                                 _Diffusion[k] *= -1;
 446                                         MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 447                                 }
 448                         }
 449                 } else  if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
 450                         // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
 451                         Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
 452                         Ctemp2 = _mesh.getCell(idCells[1]);
 453                         if (_Ndim >1){
 454                                 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
 455                                         if (j == Ctemp1.getFacesId()[l]){
 456                                                 for (int idim = 0; idim < _Ndim; ++idim)
 457                                                         _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
 458                                                 break;
 459                                         }
 460                                 }
 461                         }else{ // _Ndim = 1, build normal vector (bug cdmath)
 462                                 if(!_sectionFieldSet)
 463                                 {
 464                                         if (Fj.x()<Ctemp1.x())
 465                                                 _vec_normal[0] = -1;
 466                                         else
 467                                                 _vec_normal[0] = 1;
 468                                 }
 469                                 else
 470                                 {
 471                                         if(idCells[0]>idCells[1])
 472                                                 _vec_normal[0] = -1;
 473                                         else
 474                                                 _vec_normal[0] = 1;
 475                                 }
 476                         }
 477                         if(_verbose && _nbTimeStep%_freqSave ==0)
 478                         {
 479                                 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
 480                                 cout<<" Normal vector= ";
 481                                 for (int idim = 0; idim < _Ndim; ++idim)
 482                                         cout<<_vec_normal[idim]<<", ";
 483                                 cout<<endl;
 484                         }
 485                         _porosityi=_porosityField(idCells[0]);
 486                         _porosityj=_porosityField(idCells[1]);
 487                         convectionState(idCells[0],idCells[1],false);
 488                         convectionMatrices();
 489                         diffusionStateAndMatrices(idCells[0], idCells[1], false);
 490
 491                         // compute 1/dxi and 1/dxj
 492                         if (_Ndim > 1)
 493                         {
 494                                 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
 495                                 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
 496                         }
 497                         else
 498                         {
 499                                 _inv_dxi = 1/Ctemp1.getMeasure();
 500                                 _inv_dxj = 1/Ctemp2.getMeasure();
 501                         }
 502
 503                         addConvectionToSecondMember( idCells[0],idCells[1], false);
 504                         addDiffusionToSecondMember( idCells[0],idCells[1], false);
 505                         addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
 506
 507                         if(_timeScheme == Implicit){
 508                                 for(int k=0; k<_nVar*_nVar;k++)
 509                                 {
 510                                         _AroeMinusImplicit[k] *= _inv_dxi;
 511                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
 512                                 }
 513                                 idm = idCells[0];
 514                                 idn = idCells[1];
 515                                 //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNb<<endl;
 516                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
 517                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
 518
 519                                 if(_verbose){
 520                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
 521                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
 522                                 }
 523                                 for(int k=0;k<_nVar*_nVar;k++){
 524                                         _AroeMinusImplicit[k] *= -1;
 525                                         _Diffusion[k] *= -1;
 526                                 }
 527                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 528                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 529                                 if(_verbose){
 530                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
 531                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
 532                                 }
 533                                 for(int k=0; k<_nVar*_nVar;k++)
 534                                 {
 535                                         _AroePlusImplicit[k]  *= _inv_dxj;
 536                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;
 537                                 }
 538                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
 539                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
 540                                 if(_verbose)
 541                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
 542
 543                                 for(int k=0;k<_nVar*_nVar;k++){
 544                                         _AroePlusImplicit[k] *= -1;
 545                                         _Diffusion[k] *= -1;
 546                                 }
 547                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
 548                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
 549
 550                                 if(_verbose)
 551                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
 552                         }
 553                 }
 554                 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
 555                         if(_verbose && _nbTimeStep%_freqSave ==0)
 556                                 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
 557                         *_runLogFile<<"Warning: treatment of a junction node"<<endl;
 558
 559                         if(!_sectionFieldSet)
 560                         {
 561                                 _runLogFile->close();
 562                                 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
 563                         }
 564                         int largestSectionCellIndex=0;
 565                         for(int i=1;i<Fj.getNumberOfCells();i++){
 566                                 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
 567                                         largestSectionCellIndex=i;
 568                         }
 569                         idm = idCells[largestSectionCellIndex];
 570                         Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
 571                         _porosityi=_porosityField(idm);
 572
 573                         if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
 574                                 _vec_normal[0] = 1;
 575                         else//j==16
 576                                 _vec_normal[0] = -1;
 577                         if(_verbose && _nbTimeStep%_freqSave ==0)
 578                         {
 579                                 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
 580                                 cout << " ; vecteur normal=(";
 581                                 for(int p=0; p<_Ndim; p++)
 582                                         cout << _vec_normal[p] << ",";
 583                                 cout << "). "<<endl;
 584                         }
 585                         for(int i=0;i<Fj.getNumberOfCells();i++){
 586                                 if(i != largestSectionCellIndex){
 587                                         idn = idCells[i];
 588                                         Ctemp2 = _mesh.getCell(idn);
 589                                         _porosityj=_porosityField(idn);
 590                                         convectionState(idm,idn,false);
 591                                         convectionMatrices();
 592                                         diffusionStateAndMatrices(idm, idn,false);
 593
 594                                         if(_verbose && _nbTimeStep%_freqSave ==0)
 595                                                 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
 596
 597                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
 598                                         _inv_dxj = 1/Ctemp2.getMeasure();
 599
 600                                         addConvectionToSecondMember(idm,idn, false);
 601                                         _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
 602                                         addDiffusionToSecondMember(idm,idn, false);
 603                                         _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
 604                                         addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
 605
 606                                         if(_timeScheme == Implicit){
 607                                                 for(int k=0; k<_nVar*_nVar;k++)
 608                                                 {
 609                                                         _AroeMinusImplicit[k] *= _inv_dxi;
 610                                                         _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
 611                                                 }
 612                                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
 613                                                 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
 614
 615                                                 if(_verbose){
 616                                                         displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
 617                                                         displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
 618                                                 }
 619                                                 for(int k=0;k<_nVar*_nVar;k++){
 620                                                         _AroeMinusImplicit[k] *= -1;
 621                                                         _Diffusion[k] *= -1;
 622                                                 }
 623                                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
 624                                                 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
 625                                                 if(_verbose){
 626                                                         displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
 627                                                         displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
 628                                                 }
 629                                                 for(int k=0; k<_nVar*_nVar;k++)
 630                                                 {
 631                                                         _AroePlusImplicit[k] *= _inv_dxj;
 632                                                         _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
 633                                                 }
 634                                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
 635                                                 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
 636                                                 if(_verbose)
 637                                                         displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
 638
 639                                                 for(int k=0;k<_nVar*_nVar;k++){
 640                                                         _AroePlusImplicit[k] *= -1;
 641                                                         _Diffusion[k] *= -1;
 642                                                 }
 643                                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
 644                                                 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
 645
 646                                                 if(_verbose)
 647                                                         displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
 648                                         }
 649                                 }
 650                         }
 651                 }
 652                 else
 653                 {
 654                         cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
 655                         _runLogFile->close();
 656                         throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
 657                 }
 658
 659         }
 660         VecAssemblyEnd(_conservativeVars);
 661         VecAssemblyEnd(_b);
 662
 663         if(_timeScheme == Implicit){
 664                 for(int imaille = 0; imaille<_Nmailles; imaille++)
 665                         MatSetValuesBlocked(_A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
 666
 667                 if(_verbose && _nbTimeStep%_freqSave ==0)
 668                         displayMatrix(_GravityImplicitationMatrix,_nVar,"Gravity matrix:");
 669
 670                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 671                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 672                 if(_verbose && _nbTimeStep%_freqSave ==0){
 673                         cout << "Fin ProblemFluid.cxx : matrice implicite"<<endl;
 674                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 675                         cout << endl;
 676                 }
 677         }
 678
 679         stop=false;
 680
 681         /*
 682         if(_nbTimeStep+1<_cfl)
 683                 return (_nbTimeStep+1)*_minl/_maxvp;
 684         else
 685          */
 686         return _cfl*_minl/_maxvp;
 687 }
 688
 689 void ProblemFluid::computeNewtonVariation()
 690 {
 691         if(_verbose)
 692         {
 693                 cout<<"Vecteur courant Uk "<<endl;
 694                 VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
 695                 cout << endl;
 696         }
 697         if(_timeScheme == Explicit)
 698         {
 699                 VecCopy(_b,_newtonVariation);
 700                 VecScale(_newtonVariation, _dt);
 701                 if(_verbose && _nbTimeStep%_freqSave ==0)
 702                 {
 703                         cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
 704                         VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
 705                         cout << endl;
 706                 }
 707         }
 708         else
 709         {
 710                 if(_verbose)
 711                 {
 712                         cout << "Matrice du système linéaire avant contribution delta t" << endl;
 713                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 714                         cout << endl;
 715                         cout << "Second membre du système linéaire avant contribution delta t" << endl;
 716                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 717                         cout << endl;
 718                 }
 719                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 720                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 721
 722                 VecAXPY(_b, 1/_dt, _old);
 723                 VecAXPY(_b, -1/_dt, _conservativeVars);
 724                 MatShift(_A, 1/_dt);
 725
 726 #if PETSC_VERSION_GREATER_3_5
 727                 KSPSetOperators(_ksp, _A, _A);
 728 #else
 729                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 730 #endif
 731
 732                 if(_verbose)
 733                 {
 734                         cout << "Matrice du système linéaire après contribution delta t" << endl;
 735                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 736                         cout << endl;
 737                         cout << "Second membre du système linéaire après contribution delta t" << endl;
 738                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 739                         cout << endl;
 740                 }
 741
 742                 if(_conditionNumber)
 743                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 744                 if(!_isScaling)
 745                 {
 746                         KSPSolve(_ksp, _b, _newtonVariation);
 747                 }
 748                 else
 749                 {
 750                         computeScaling(_maxvp);
 751                         int indice;
 752                         VecAssemblyBegin(_vecScaling);
 753                         VecAssemblyBegin(_invVecScaling);
 754                         for(int imaille = 0; imaille<_Nmailles; imaille++)
 755                         {
 756                                 indice = imaille;
 757                                 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
 758                                 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
 759                         }
 760                         VecAssemblyEnd(_vecScaling);
 761                         VecAssemblyEnd(_invVecScaling);
 762
 763                         if(_system)
 764                         {
 765                                 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
 766                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 767                                 cout << endl;
 768                                 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
 769                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 770                                 cout << endl;
 771                         }
 772                         MatDiagonalScale(_A,_vecScaling,_invVecScaling);
 773                         if(_system)
 774                         {
 775                                 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
 776                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 777                                 cout << endl;
 778                         }
 779                         VecCopy(_b,_bScaling);
 780                         VecPointwiseMult(_b,_vecScaling,_bScaling);
 781                         if(_system)
 782                         {
 783                                 cout << "Produit du second membre par le preconditionneur bloc diagonal  a gauche" << endl;
 784                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 785                                 cout << endl;
 786                         }
 787
 788                         KSPSolve(_ksp,_b, _bScaling);
 789                         VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
 790                 }
 791                 if(_verbose)
 792                 {
 793                         cout << "solution du systeme lineaire local:" << endl;
 794                         VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
 795                         cout << endl;
 796                 }
 797         }
 798 }
 799
 800 void ProblemFluid::validateTimeStep()
 801 {
 802         if(_system)
 803         {
 804                 cout<<" Vecteur Un"<<endl;
 805                 VecView(_old,  PETSC_VIEWER_STDOUT_WORLD);
 806                 cout<<" Vecteur Un+1"<<endl;
 807                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
 808         }
 809         VecAXPY(_old,  -1, _conservativeVars);//old contient old-courant
 810
 811         //Calcul de la variation Un+1-Un
 812         _erreur_rel= 0;
 813         double x, dx;
 814         int I;
 815
 816         for(int j=1; j<=_Nmailles; j++)
 817         {
 818                 for(int k=0; k<_nVar; k++)
 819                 {
 820                         I = (j-1)*_nVar + k;
 821                         VecGetValues(_old, 1, &I, &dx);
 822                         VecGetValues(_conservativeVars, 1, &I, &x);
 823                         if (fabs(x)< _precision)
 824                         {
 825                                 if(_erreur_rel < fabs(dx))
 826                                         _erreur_rel = fabs(dx);
 827                         }
 828                         else if(_erreur_rel < fabs(dx/x))
 829                                 _erreur_rel = fabs(dx/x);
 830                 }
 831         }
 832
 833         _isStationary =_erreur_rel <_precision;
 834
 835         VecCopy(_conservativeVars, _old);
 836
 837         if(_verbose && _nbTimeStep%_freqSave ==0){
 838                 if(!_usePrimitiveVarsInNewton)
 839                         testConservation();
 840                 cout <<"Valeur propre locale max: " << _maxvp << endl;
 841         }
 842
 843         if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
 844                 //Find minimum and maximum void fractions
 845                 double alphamin=1e30;
 846                 double alphamax=-1e30;
 847                 double T, Tmax=-1e30;
 848                 int J=_nVar-1;
 849                 I = 0;
 850                 for(int j=0; j<_Nmailles; j++)
 851                 {
 852                         VecGetValues(_primitiveVars, 1, &I, &x);//extract void fraction
 853                         if(x>alphamax)
 854                                 alphamax=x;
 855                         if(x<alphamin)
 856                                 alphamin=x;
 857                         I += _nVar;
 858                         VecGetValues(_primitiveVars, 1, &J, &T);//extract void fraction
 859                         if(T>Tmax)
 860                                 Tmax=T;
 861                         J += _nVar;
 862                 }
 863                 cout<<"Alpha min = " << alphamin << ", Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
 864                 cout<<endl;
 865                 *_runLogFile<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
 866         }
 867
 868         _time+=_dt;
 869         _nbTimeStep++;
 870         if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 871                 save();
 872 }
 873
 874 void ProblemFluid::abortTimeStep(){
 875         _dt = 0;
 876 }
 877
 878 void ProblemFluid::addConvectionToSecondMember
 879 (               const int &i,
 880                 const int &j, bool isBord, string groupname
 881 )
 882 {
 883         if(_verbose && _nbTimeStep%_freqSave ==0)
 884                 cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
 885
 886         //extraction des valeurs
 887         for(int k=0; k<_nVar; k++)
 888                 _idm[k] = _nVar*i + k;
 889         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 890         VecGetValues(_primitiveVars,    _nVar, _idm, _Vi);
 891
 892         if(!isBord){
 893                 for(int k=0; k<_nVar; k++)
 894                         _idn[k] = _nVar*j + k;
 895                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
 896                 VecGetValues(_primitiveVars,    _nVar, _idn, _Vj);
 897         }
 898         else{
 899                 for(int k=0; k<_nVar; k++)
 900                         _idn[k] = k;
 901                 VecGetValues(_Uext, _nVar, _idn, _Uj);
 902                 consToPrim(_Uj, _Vj,_porosityj);
 903         }
 904         _idm[0] = i;
 905         _idn[0] = j;
 906
 907         if(_verbose && _nbTimeStep%_freqSave ==0)
 908         {
 909                 cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
 910                 for(int q=0; q<_nVar; q++)
 911                         cout << _Ui[q] << endl;
 912                 cout << endl;
 913                 cout << "addConvectionToSecondMember : état j= " << j << ", _Uj=" << endl;
 914                 for(int q=0; q<_nVar; q++)
 915                         cout << _Uj[q] << endl;
 916                 cout << endl;
 917         }
 918         if(_nonLinearFormulation!=reducedRoe){//VFRoe, Roe or VFFC
 919                 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar),Fij(_nVar);
 920                 for(int i1=0;i1<_nVar;i1++)
 921                 {
 922                         Ui(i1)=_Ui[i1];
 923                         Uj(i1)=_Uj[i1];
 924                         Vi(i1)=_Vi[i1];
 925                         Vj(i1)=_Vj[i1];
 926                 }
 927                 Vector normale(_Ndim);
 928                 for(int i1=0;i1<_Ndim;i1++)
 929                         normale(i1)=_vec_normal[i1];
 930
 931                 Matrix signAroe(_nVar);
 932                 for(int i1=0;i1<_nVar;i1++)
 933                         for(int i2=0;i2<_nVar;i2++)
 934                                 signAroe(i1,i2)=_signAroe[i1*_nVar+i2];
 935
 936                 Matrix absAroe(_nVar);
 937                 for(int i1=0;i1<_nVar;i1++)
 938                         for(int i2=0;i2<_nVar;i2++)
 939                                 absAroe(i1,i2)=_absAroe[i1*_nVar+i2];
 940
 941                 if(_nonLinearFormulation==VFRoe)//VFRoe
 942                 {
 943                         Vector Uij(_nVar),Vij(_nVar);
 944                         double porosityij=sqrt(_porosityi*_porosityj);
 945
 946                         Uij=(Ui+Uj)/2+signAroe*(Ui-Uj)/2;
 947
 948                         for(int i1=0;i1<_nVar;i1++)
 949                                 _Uij[i1]=Uij(i1);
 950
 951                         consToPrim(_Uij, _Vij,porosityij);
 952
 953                         applyVFRoeLowMachCorrections(isBord, groupname);
 954
 955                         for(int i1=0;i1<_nVar;i1++)
 956                         {
 957                                 Uij(i1)=_Uij[i1];
 958                                 Vij(i1)=_Vij[i1];
 959                         }
 960
 961                         Fij=convectionFlux(Uij,Vij,normale,porosityij);
 962
 963                         if(_verbose && _nbTimeStep%_freqSave ==0)
 964                         {
 965                                 cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
 966                                 cout<<Uij<<endl;
 967                                 cout<<"Etat interfacial primitif "<<i<<", "<<j<< endl;
 968                                 cout<<Vij<<endl;
 969                         }
 970                 }
 971                 else //Roe or VFFC
 972                 {
 973                         if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
 974                                 Fij=staggeredVFFCFlux();
 975                         else
 976                         {
 977                                 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
 978                                 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
 979                                 if(_nonLinearFormulation==VFFC)//VFFC
 980                                         Fij=(Fi+Fj)/2+signAroe*(Fi-Fj)/2;
 981                                 else if(_nonLinearFormulation==Roe)//Roe
 982                                         Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
 983
 984                                 if(_verbose && _nbTimeStep%_freqSave ==0)
 985                                 {
 986                                         cout<<"Flux cellule "<<i<<" = "<< endl;
 987                                         cout<<Fi<<endl;
 988                                         cout<<"Flux cellule "<<j<<" = "<< endl;
 989                                         cout<<Fj<<endl;
 990                                 }
 991                         }
 992                 }
 993                 for(int i1=0;i1<_nVar;i1++)
 994                         _phi[i1]=-Fij(i1)*_inv_dxi;
 995                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
 996                 if(_verbose && _nbTimeStep%_freqSave ==0)
 997                 {
 998                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
 999                         cout<<"Flux interfacial "<<i<<", "<<j<< endl;
1000                         cout<<Fij<<endl;
1001                         cout << "Contribution convection à " << i << ", -Fij(i1)*_inv_dxi= "<<endl;
1002                         for(int q=0; q<_nVar; q++)
1003                                 cout << _phi[q] << endl;
1004                         cout << endl;
1005                 }
1006                 if(!isBord){
1007                         for(int i1=0;i1<_nVar;i1++)
1008                                 _phi[i1]*=-_inv_dxj/_inv_dxi;
1009                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1010                         if(_verbose && _nbTimeStep%_freqSave ==0)
1011                         {
1012                                 cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
1013                                 for(int q=0; q<_nVar; q++)
1014                                         cout << _phi[q] << endl;
1015                                 cout << endl;
1016                         }
1017                 }
1018         }else //_nonLinearFormulation==reducedRoe)
1019         {
1020                 for(int k=0; k<_nVar; k++)
1021                         _temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
1022                 Polynoms::matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
1023                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1024
1025                 if(_verbose && _nbTimeStep%_freqSave ==0)
1026                 {
1027                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1028                         cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
1029                         for(int q=0; q<_nVar; q++)
1030                                 cout << _temp[q] << endl;
1031                         cout << endl;
1032                         cout << "Contribution convection à " << i << ", A^-*(Ui - Uj)*_inv_dxi= "<<endl;
1033                         for(int q=0; q<_nVar; q++)
1034                                 cout << _phi[q] << endl;
1035                         cout << endl;
1036                 }
1037
1038                 if(!isBord)
1039                 {
1040                         for(int k=0; k<_nVar; k++)
1041                                 _temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
1042                         Polynoms::matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
1043                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1044
1045                         if(_verbose && _nbTimeStep%_freqSave ==0)
1046                         {
1047                                 cout << "Contribution convection à  " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
1048                                 for(int q=0; q<_nVar; q++)
1049                                         cout << _phi[q] << endl;
1050                                 cout << endl;
1051                         }
1052                 }
1053         }
1054         if(_verbose && _nbTimeStep%_freqSave ==0)
1055         {
1056                 cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
1057                 displayMatrix(_absAroe,   _nVar,"Valeur absolue matrice de Roe");
1058                 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
1059                 displayMatrix(_AroePlus,  _nVar,"Matrice _AroePlus");
1060                 displayMatrix(_signAroe,  _nVar,"Signe de la matrice de Roe");
1061         }
1062 }
1063
1064 void ProblemFluid::addSourceTermToSecondMember
1065 (       const int i, int nbVoisinsi,
1066                 const int j, int nbVoisinsj,
1067                 bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
1068 {
1069         if(_verbose && _nbTimeStep%_freqSave ==0)
1070                 cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
1071
1072         _idm[0] = i*_nVar;
1073         for(int k=1; k<_nVar;k++)
1074                 _idm[k] = _idm[k-1] + 1;
1075         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1076         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1077         sourceVector(_Si,_Ui,_Vi,i);
1078
1079         if (_verbose && _nbTimeStep%_freqSave ==0)
1080         {
1081                 cout << "Terme source S(Ui), i= " << i<<endl;
1082                 for(int q=0; q<_nVar; q++)
1083                         cout << _Si[q] << endl;
1084                 cout << endl;
1085         }
1086         if(!isBord){
1087                 for(int k=0; k<_nVar; k++)
1088                         _idn[k] = _nVar*j + k;
1089                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1090                 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1091                 sourceVector(_Sj,_Uj,_Vj,j);
1092         }else
1093         {
1094                 for(int k=0; k<_nVar; k++)
1095                         _idn[k] = k;
1096                 VecGetValues(_Uext, _nVar, _idn, _Uj);
1097                 consToPrim(_Uj, _Vj,_porosityj);
1098                 sourceVector(_Sj,_Uj,_Vj,i);
1099         }
1100         if (_verbose && _nbTimeStep%_freqSave ==0)
1101         {
1102                 if(!isBord)
1103                         cout << "Terme source S(Uj), j= " << j<<endl;
1104                 else
1105                         cout << "Terme source S(Uj), cellule fantôme "<<endl;
1106                 for(int q=0; q<_nVar; q++)
1107                         cout << _Sj[q] << endl;
1108                 cout << endl;
1109         }
1110
1111         //Compute pressure loss vector
1112         double K;
1113         if(_pressureLossFieldSet){
1114                 K=_pressureLossField(ij);
1115                 pressureLossVector(_pressureLossVector, K,_Ui, _Vi, _Uj, _Vj);
1116         }
1117         else{
1118                 K=0;
1119                 for(int k=0; k<_nVar;k++)
1120                         _pressureLossVector[k]=0;
1121         }
1122         if (_verbose && _nbTimeStep%_freqSave ==0)
1123         {
1124                 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
1125                 for(int k=0; k<_nVar;k++)
1126                         cout<< _pressureLossVector[k]<<", ";
1127                 cout<<endl;
1128         }
1129         //Contribution of the porosityField gradient:
1130         if(!isBord)
1131                 porosityGradientSourceVector();
1132         else{
1133                 for(int k=0; k<_nVar;k++)
1134                         _porosityGradientSourceVector[k]=0;
1135         }
1136
1137         if (_verbose && _nbTimeStep%_freqSave ==0)
1138         {
1139                 if(!isBord)
1140                         cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
1141                 else
1142                         cout<<"interface frontière i= "<<i<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<endl;
1143                 cout<<"Gradient de porosite à l'interface"<<endl;
1144                 for(int k=0; k<_nVar;k++)
1145                         cout<< _porosityGradientSourceVector[k]<<", ";
1146                 cout<<endl;
1147         }
1148
1149         if(!isBord){
1150                 if(_wellBalancedCorrection){
1151                         for(int k=0; k<_nVar;k++)
1152                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1153                         Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1154                         for(int k=0; k<_nVar;k++){
1155                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1156                                 _Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
1157                         }
1158                         if (_verbose && _nbTimeStep%_freqSave ==0)
1159                         {
1160                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,j), j="<<j<<endl;
1161                                 for(int q=0; q<_nVar; q++)
1162                                         cout << _Si[q] << endl;
1163                                 cout << "Contribution au terme source Sj de la cellule j= " << j<<" venant  (après décentrement) de la face (i,j), i="<<i<<endl;
1164                                 for(int q=0; q<_nVar; q++)
1165                                         cout << _Sj[q] << endl;
1166                                 cout << endl;
1167                                 cout<<"ratio surface sur volume i = "<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1168                                 cout<<"ratio surface sur volume j = "<<mesureFace/_perimeters(j)<<" perimeter = "<< _perimeters(j)<<endl;
1169                                 cout << endl;
1170                         }
1171                 }
1172                 else{
1173                         for(int k=0; k<_nVar;k++){
1174                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
1175                                 _Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
1176                         }
1177                         if (_verbose && _nbTimeStep%_freqSave ==0)
1178                         {
1179                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant  de la face (i,j), j="<<j<<endl;
1180                                 for(int q=0; q<_nVar; q++)
1181                                         cout << _Si[q] << endl;
1182                                 cout << "Contribution au terme source Sj de la cellule j = " << j<<" venant  de la face (i,j), i="<<i <<endl;
1183                                 for(int q=0; q<_nVar; q++)
1184                                         cout << _Sj[q] << endl;
1185                                 cout << endl;
1186                         }
1187                 }
1188                 _idn[0] = j;
1189                 VecSetValuesBlocked(_b, 1, _idn, _Sj, ADD_VALUES);
1190         }else{
1191                 if(_wellBalancedCorrection){
1192                         for(int k=0; k<_nVar;k++)
1193                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1194                         Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1195                         for(int k=0; k<_nVar;k++)
1196                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1197                         if (_verbose && _nbTimeStep%_freqSave ==0)
1198                         {
1199                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,bord)"<<endl;
1200                                 for(int q=0; q<_nVar; q++)
1201                                         cout << _Si[q] << endl;
1202                                 cout<<"ratio surface sur volume i ="<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1203                                 cout << endl;
1204                         }
1205                 }
1206                 else
1207                 {
1208                         for(int k=0; k<_nVar;k++)
1209                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
1210                         if (_verbose && _nbTimeStep%_freqSave ==0)
1211                         {
1212                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
1213                                 for(int q=0; q<_nVar; q++)
1214                                         cout << _Si[q] << endl;
1215                                 cout << endl;
1216                         }
1217                 }
1218         }
1219         _idm[0] = i;
1220         VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
1221
1222         if(_verbose && _nbTimeStep%_freqSave ==0 && _wellBalancedCorrection)
1223                 displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
1224 }
1225
1226 void ProblemFluid::updatePrimitives()
1227 {
1228         VecAssemblyBegin(_primitiveVars);
1229         for(int i=1; i<=_Nmailles; i++)
1230         {
1231                 _idm[0] = _nVar*( (i-1));
1232                 for(int k=1; k<_nVar; k++)
1233                         _idm[k] = _idm[k-1] + 1;
1234
1235                 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1236                 consToPrim(_Ui,_Vi,_porosityField(i-1));
1237                 if(_verbose && _nbTimeStep%_freqSave ==0)
1238                 {
1239                         cout << "ProblemFluid::updatePrimitives() cell " << i-1 << endl;
1240                         cout << "Ui = ";
1241                         for(int q=0; q<_nVar; q++)
1242                                 cout << _Ui[q]  << "\t";
1243                         cout << endl;
1244                         cout << "Vi = ";
1245                         for(int q=0; q<_nVar; q++)
1246                                 cout << _Vi[q]  << "\t";
1247                         cout << endl;
1248                 }
1249
1250                 if(_nbPhases==2 && _Psat>-1e30){//Cas simulation flashing
1251                         double pressure;
1252                         VecGetValues(_primitiveVars, 1, _idm+1, &pressure);
1253                         _dp_over_dt(i-1)=(_Vi[1]-pressure)/_dt;//pn+1-pn
1254                 }
1255                 _idm[0] = i-1;
1256                 VecSetValuesBlocked(_primitiveVars, 1, _idm, _Vi, INSERT_VALUES);
1257         }
1258         VecAssemblyEnd(_primitiveVars);
1259
1260         if(_system)
1261         {
1262                 cout << "Nouvelles variables primitives : " << endl;
1263                 VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_WORLD);
1264                 cout << endl;
1265         }
1266 }
1267
1268 void ProblemFluid::updateConservatives()
1269 {
1270         VecAssemblyBegin(_conservativeVars);
1271         for(int i=1; i<=_Nmailles; i++)
1272         {
1273                 _idm[0] = _nVar*( (i-1));
1274                 for(int k=1; k<_nVar; k++)
1275                         _idm[k] = _idm[k-1] + 1;
1276
1277                 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1278
1279                 primToCons(_Vi,0,_Ui,0);
1280                 _idm[0] = i-1;
1281                 VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
1282
1283                 if(_verbose && _nbTimeStep%_freqSave ==0)
1284                 {
1285                         cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
1286                         cout << "Vi = ";
1287                         for(int q=0; q<_nVar; q++)
1288                                 cout << _Vi[q]  << "\t";
1289                         cout << endl;
1290                         cout << "Ui = ";
1291                         for(int q=0; q<_nVar; q++)
1292                                 cout << _Ui[q]  << "\t";
1293                         cout << endl;
1294                 }
1295         }
1296         VecAssemblyEnd(_conservativeVars);
1297
1298         if(_system)
1299         {
1300                 cout << "Nouvelles variables conservatives : " << endl;
1301                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
1302                 cout << endl;
1303         }
1304 }
1305
1306 vector< complex<double> > ProblemFluid::getRacines(vector< double > pol_car){
1307         int degre_polynom = pol_car.size()-1;
1308         double tmp;
1309         vector< complex< double > > valeurs_propres (degre_polynom);
1310         vector< double > zeror(degre_polynom);
1311         vector< double > zeroi(degre_polynom);
1312         for(int j=0; j<(degre_polynom+1)/2; j++){//coefficients in order of decreasing powers for rpoly
1313                 tmp=pol_car[j];
1314                 pol_car[j] =pol_car[degre_polynom-j];
1315                 pol_car[degre_polynom-j]=tmp;
1316         }
1317
1318         //Calcul des racines du polynome
1319         roots_polynoms roots;
1320         int size_vp = roots.rpoly(&pol_car[0],degre_polynom,&zeror[0],&zeroi[0]);
1321
1322         //On ordonne les valeurs propres
1323         if(zeror[1]<zeror[0])
1324         {
1325                 tmp=zeror[0];
1326                 zeror[0]=zeror[1];
1327                 zeror[1]=tmp;
1328                 tmp=zeroi[0];
1329                 zeroi[0]=zeroi[1];
1330                 zeroi[1]=tmp;
1331         }
1332
1333         if(size_vp == degre_polynom) {
1334                 for(int ct =0; ct<degre_polynom; ct++)
1335                         valeurs_propres[ct] = complex<double> (zeror[ct],zeroi[ct]);  //vp non triviales
1336         }
1337         else {
1338                 cout << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1339                 *_runLogFile << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1340                 cout<<"Coefficients polynome caracteristique: "<<endl;
1341                 for(int ct =0; ct<degre_polynom+1; ct++)
1342                         cout<<pol_car[ct]<< " , " <<endl;
1343
1344                 *_runLogFile<<"getRacines computation failed"<<endl;
1345                 _runLogFile->close();
1346                 throw CdmathException("getRacines computation failed");
1347         }
1348
1349         return valeurs_propres;
1350 }
1351
1352 void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1353 {
1354         int nbVp_dist=valeurs_propres_dist.size();
1355         vector< complex< double > > y (nbVp_dist,0);
1356         for( int i=0 ; i<nbVp_dist ; i++)
1357                 y[i] = Polynoms::abs_generalise(valeurs_propres_dist[i]);
1358         Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
1359         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1360         {
1361                 cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1362                 for(int ct =0; ct<nbVp_dist; ct++)
1363                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1364                 cout<< endl;
1365                 cout<<"ProblemFluid::AbsMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1366                 for(int ct =0; ct<nbVp_dist; ct++)
1367                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1368                 cout<< endl;
1369                 displayMatrix(_absAroe,_nVar,"Valeur absolue de la matrice de Roe");
1370         }
1371 }
1372
1373 void ProblemFluid::SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1374 {
1375         int nbVp_dist=valeurs_propres_dist.size();
1376         vector< complex< double > > y (nbVp_dist,0);
1377         for( int i=0 ; i<nbVp_dist ; i++)
1378         {
1379                 if(valeurs_propres_dist[i].real()>0)
1380                         y[i] = 1;
1381                 else if(valeurs_propres_dist[i].real()<0)
1382                         y[i] = -1;
1383                 else
1384                         y[i] = 0;
1385         }
1386
1387         Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
1388         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1389         {
1390                 cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1391                 for(int ct =0; ct<nbVp_dist; ct++)
1392                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1393                 cout<< endl;
1394                 cout<<"ProblemFluid::SigneMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1395                 for(int ct =0; ct<nbVp_dist; ct++)
1396                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1397                 cout<< endl;
1398                 displayMatrix(_signAroe,_nVar,"Signe matrice de Roe");
1399         }
1400 }
1401 void ProblemFluid::InvMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1402 {
1403         int nbVp_dist=valeurs_propres_dist.size();
1404         vector< complex< double > > y (nbVp_dist,0);
1405
1406         for( int i=0 ; i<nbVp_dist ; i++)
1407         {
1408                 if(Polynoms::module(valeurs_propres_dist[i])>_precision)
1409                         y[i] = 1./valeurs_propres_dist[i];
1410                 else
1411                         y[i] = 1./_precision;
1412         }
1413         Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
1414         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1415         {
1416                 cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
1417                 for(int ct =0; ct<nbVp_dist; ct++)
1418                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1419                 cout<< endl;
1420                 cout<<"ProblemFluid::InvMatriceRoe : Valeurs à atteindre pour l'interpolation"<<endl;
1421                 for(int ct =0; ct<nbVp_dist; ct++)
1422                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1423                 cout<< endl;
1424                 displayMatrix(_invAroe,_nVar,"Inverse matrice de Roe");
1425         }
1426 }
1427
1428 Field ProblemFluid::getConservativeField() const
1429 {
1430         if(!_initializedMemory)
1431         {
1432                 _runLogFile->close();
1433                 throw CdmathException("ProblemFluid::getConservativeField call initialize() first");
1434         }
1435         return _UU;
1436 }
1437
1438 Field ProblemFluid::getPrimitiveField() const
1439 {
1440         if(!_initializedMemory)
1441         {
1442                 _runLogFile->close();
1443                 throw CdmathException("ProblemFluid::getPrimitiveField call initialize() first");
1444         }
1445         return _VV;
1446 }
1447
1448 void ProblemFluid::terminate(){
1449
1450         delete[] _AroePlus;
1451         delete[] _Diffusion;
1452         delete[] _GravityImplicitationMatrix;
1453         delete[] _AroeMinus;
1454         delete[] _Aroe;
1455         delete[] _absAroe;
1456         delete[] _signAroe;
1457         delete[] _invAroe;
1458         delete[] _AroeImplicit;
1459         delete[] _AroeMinusImplicit;
1460         delete[] _AroePlusImplicit;
1461         delete[] _absAroeImplicit;
1462         delete[] _phi;
1463         delete[] _Jcb;
1464         delete[] _JcbDiff;
1465         delete[] _a;
1466         delete[] _primToConsJacoMat;
1467
1468         delete[] _l;
1469         delete[] _r;
1470         delete[] _Uroe;
1471         delete[] _Udiff;
1472         delete[] _temp;
1473         delete[] _externalStates;
1474         delete[] _idm;
1475         delete[] _idn;
1476         delete[] _vec_normal;
1477         delete[] _Ui;
1478         delete[] _Uj;
1479         delete[] _Vi;
1480         delete[] _Vj;
1481         if(_nonLinearFormulation==VFRoe){
1482                 delete[] _Uij;
1483                 delete[] _Vij;
1484         }
1485         delete[] _Si;
1486         delete[] _Sj;
1487         delete[] _pressureLossVector;
1488         delete[] _porosityGradientSourceVector;
1489         if(_isScaling)
1490         {
1491                 delete[] _blockDiag;
1492                 delete[] _invBlockDiag;
1493
1494                 VecDestroy(&_vecScaling);
1495                 VecDestroy(&_invVecScaling);
1496                 VecDestroy(&_bScaling);
1497         }
1498
1499         VecDestroy(&_conservativeVars);
1500         VecDestroy(&_newtonVariation);
1501         VecDestroy(&_b);
1502         VecDestroy(&_primitiveVars);
1503         VecDestroy(&_Uext);
1504         VecDestroy(&_Vext);
1505         VecDestroy(&_Uextdiff);
1506
1507         //      PCDestroy(_pc);
1508         KSPDestroy(&_ksp);
1509         for(int i=0;i<_nbPhases;i++)
1510                 delete _fluides[i];
1511 }