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                                         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         if(_nbTimeStep+1<_cfl)
 684                 return (_nbTimeStep+1)*_minl/_maxvp;
 685         else
 686          */
 687         return _cfl*_minl/_maxvp;
 688 }
 689
 690 void ProblemFluid::computeNewtonVariation()
 691 {
 692         if(_verbose)
 693         {
 694                 cout<<"Vecteur courant Uk "<<endl;
 695                 VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
 696                 cout << endl;
 697         }
 698         if(_timeScheme == Explicit)
 699         {
 700                 VecCopy(_b,_newtonVariation);
 701                 VecScale(_newtonVariation, _dt);
 702                 if(_verbose && _nbTimeStep%_freqSave ==0)
 703                 {
 704                         cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
 705                         VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
 706                         cout << endl;
 707                 }
 708         }
 709         else
 710         {
 711                 if(_verbose)
 712                 {
 713                         cout << "Matrice du système linéaire avant contribution delta t" << endl;
 714                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 715                         cout << endl;
 716                         cout << "Second membre du système linéaire avant contribution delta t" << endl;
 717                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 718                         cout << endl;
 719                 }
 720                 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 721                 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 722
 723                 VecAXPY(_b, 1/_dt, _old);
 724                 VecAXPY(_b, -1/_dt, _conservativeVars);
 725                 MatShift(_A, 1/_dt);
 726
 727 #if PETSC_VERSION_GREATER_3_5
 728                 KSPSetOperators(_ksp, _A, _A);
 729 #else
 730                 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 731 #endif
 732
 733                 if(_verbose)
 734                 {
 735                         cout << "Matrice du système linéaire après contribution delta t" << endl;
 736                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 737                         cout << endl;
 738                         cout << "Second membre du système linéaire après contribution delta t" << endl;
 739                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 740                         cout << endl;
 741                 }
 742
 743                 if(_conditionNumber)
 744                         KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 745                 if(!_isScaling)
 746                 {
 747                         KSPSolve(_ksp, _b, _newtonVariation);
 748                 }
 749                 else
 750                 {
 751                         computeScaling(_maxvp);
 752                         int indice;
 753                         VecAssemblyBegin(_vecScaling);
 754                         VecAssemblyBegin(_invVecScaling);
 755                         for(int imaille = 0; imaille<_Nmailles; imaille++)
 756                         {
 757                                 indice = imaille;
 758                                 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
 759                                 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
 760                         }
 761                         VecAssemblyEnd(_vecScaling);
 762                         VecAssemblyEnd(_invVecScaling);
 763
 764                         if(_system)
 765                         {
 766                                 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
 767                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 768                                 cout << endl;
 769                                 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
 770                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 771                                 cout << endl;
 772                         }
 773                         MatDiagonalScale(_A,_vecScaling,_invVecScaling);
 774                         if(_system)
 775                         {
 776                                 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
 777                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 778                                 cout << endl;
 779                         }
 780                         VecCopy(_b,_bScaling);
 781                         VecPointwiseMult(_b,_vecScaling,_bScaling);
 782                         if(_system)
 783                         {
 784                                 cout << "Produit du second membre par le preconditionneur bloc diagonal  a gauche" << endl;
 785                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 786                                 cout << endl;
 787                         }
 788
 789                         KSPSolve(_ksp,_b, _bScaling);
 790                         VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
 791                 }
 792                 if(_verbose)
 793                 {
 794                         cout << "solution du systeme lineaire local:" << endl;
 795                         VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
 796                         cout << endl;
 797                 }
 798         }
 799 }
 800
 801 void ProblemFluid::validateTimeStep()
 802 {
 803         if(_system)
 804         {
 805                 cout<<" Vecteur Un"<<endl;
 806                 VecView(_old,  PETSC_VIEWER_STDOUT_WORLD);
 807                 cout<<" Vecteur Un+1"<<endl;
 808                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
 809         }
 810         VecAXPY(_old,  -1, _conservativeVars);//old contient old-courant
 811
 812         //Calcul de la variation Un+1-Un
 813         _erreur_rel= 0;
 814         double x, dx;
 815         int I;
 816
 817         for(int j=1; j<=_Nmailles; j++)
 818         {
 819                 for(int k=0; k<_nVar; k++)
 820                 {
 821                         I = (j-1)*_nVar + k;
 822                         VecGetValues(_old, 1, &I, &dx);
 823                         VecGetValues(_conservativeVars, 1, &I, &x);
 824                         if (fabs(x)< _precision)
 825                         {
 826                                 if(_erreur_rel < fabs(dx))
 827                                         _erreur_rel = fabs(dx);
 828                         }
 829                         else if(_erreur_rel < fabs(dx/x))
 830                                 _erreur_rel = fabs(dx/x);
 831                 }
 832         }
 833
 834         _isStationary =_erreur_rel <_precision;
 835
 836         VecCopy(_conservativeVars, _old);
 837
 838         if(_verbose && _nbTimeStep%_freqSave ==0){
 839                 if(!_usePrimitiveVarsInNewton)
 840                         testConservation();
 841                 cout <<"Valeur propre locale max: " << _maxvp << endl;
 842         }
 843
 844         if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
 845                 //Find minimum and maximum void fractions
 846                 double alphamin=1e30;
 847                 double alphamax=-1e30;
 848                 double T, Tmax=-1e30;
 849                 int J=_nVar-1;
 850                 I = 0;
 851                 for(int j=0; j<_Nmailles; j++)
 852                 {
 853                         VecGetValues(_primitiveVars, 1, &I, &x);//extract void fraction
 854                         if(x>alphamax)
 855                                 alphamax=x;
 856                         if(x<alphamin)
 857                                 alphamin=x;
 858                         I += _nVar;
 859                         VecGetValues(_primitiveVars, 1, &J, &T);//extract void fraction
 860                         if(T>Tmax)
 861                                 Tmax=T;
 862                         J += _nVar;
 863                 }
 864                 cout<<"Alpha min = " << alphamin << ", Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
 865                 cout<<endl;
 866                 *_runLogFile<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
 867         }
 868
 869         _time+=_dt;
 870         _nbTimeStep++;
 871         if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
 872                 save();
 873 }
 874
 875 void ProblemFluid::abortTimeStep(){
 876         _dt = 0;
 877 }
 878
 879 void ProblemFluid::addConvectionToSecondMember
 880 (               const int &i,
 881                 const int &j, bool isBord, string groupname
 882 )
 883 {
 884         if(_verbose && _nbTimeStep%_freqSave ==0)
 885                 cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
 886
 887         //extraction des valeurs
 888         for(int k=0; k<_nVar; k++)
 889                 _idm[k] = _nVar*i + k;
 890         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 891         VecGetValues(_primitiveVars,    _nVar, _idm, _Vi);
 892
 893         if(!isBord){
 894                 for(int k=0; k<_nVar; k++)
 895                         _idn[k] = _nVar*j + k;
 896                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
 897                 VecGetValues(_primitiveVars,    _nVar, _idn, _Vj);
 898         }
 899         else{
 900                 for(int k=0; k<_nVar; k++)
 901                         _idn[k] = k;
 902                 VecGetValues(_Uext, _nVar, _idn, _Uj);
 903                 consToPrim(_Uj, _Vj,_porosityj);
 904         }
 905         _idm[0] = i;
 906         _idn[0] = j;
 907
 908         if(_verbose && _nbTimeStep%_freqSave ==0)
 909         {
 910                 cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
 911                 for(int q=0; q<_nVar; q++)
 912                         cout << _Ui[q] << endl;
 913                 cout << endl;
 914                 cout << "addConvectionToSecondMember : état j= " << j << ", _Uj=" << endl;
 915                 for(int q=0; q<_nVar; q++)
 916                         cout << _Uj[q] << endl;
 917                 cout << endl;
 918         }
 919         if(_nonLinearFormulation!=reducedRoe){//VFRoe, Roe or VFFC
 920                 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar),Fij(_nVar);
 921                 for(int i1=0;i1<_nVar;i1++)
 922                 {
 923                         Ui(i1)=_Ui[i1];
 924                         Uj(i1)=_Uj[i1];
 925                         Vi(i1)=_Vi[i1];
 926                         Vj(i1)=_Vj[i1];
 927                 }
 928                 Vector normale(_Ndim);
 929                 for(int i1=0;i1<_Ndim;i1++)
 930                         normale(i1)=_vec_normal[i1];
 931
 932                 Matrix signAroe(_nVar);
 933                 for(int i1=0;i1<_nVar;i1++)
 934                         for(int i2=0;i2<_nVar;i2++)
 935                                 signAroe(i1,i2)=_signAroe[i1*_nVar+i2];
 936
 937                 Matrix absAroe(_nVar);
 938                 for(int i1=0;i1<_nVar;i1++)
 939                         for(int i2=0;i2<_nVar;i2++)
 940                                 absAroe(i1,i2)=_absAroe[i1*_nVar+i2];
 941
 942                 if(_nonLinearFormulation==VFRoe)//VFRoe
 943                 {
 944                         Vector Uij(_nVar),Vij(_nVar);
 945                         double porosityij=sqrt(_porosityi*_porosityj);
 946
 947                         Uij=(Ui+Uj)/2+signAroe*(Ui-Uj)/2;
 948
 949                         for(int i1=0;i1<_nVar;i1++)
 950                                 _Uij[i1]=Uij(i1);
 951
 952                         consToPrim(_Uij, _Vij,porosityij);
 953
 954                         applyVFRoeLowMachCorrections(isBord, groupname);
 955
 956                         for(int i1=0;i1<_nVar;i1++)
 957                         {
 958                                 Uij(i1)=_Uij[i1];
 959                                 Vij(i1)=_Vij[i1];
 960                         }
 961
 962                         Fij=convectionFlux(Uij,Vij,normale,porosityij);
 963
 964                         if(_verbose && _nbTimeStep%_freqSave ==0)
 965                         {
 966                                 cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
 967                                 cout<<Uij<<endl;
 968                                 cout<<"Etat interfacial primitif "<<i<<", "<<j<< endl;
 969                                 cout<<Vij<<endl;
 970                         }
 971                 }
 972                 else //Roe or VFFC
 973                 {
 974                         if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
 975                                 Fij=staggeredVFFCFlux();
 976                         else
 977                         {
 978                                 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
 979                                 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
 980                                 if(_nonLinearFormulation==VFFC)//VFFC
 981                                         Fij=(Fi+Fj)/2+signAroe*(Fi-Fj)/2;
 982                                 else if(_nonLinearFormulation==Roe)//Roe
 983                                         Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
 984
 985                                 if(_verbose && _nbTimeStep%_freqSave ==0)
 986                                 {
 987                                         cout<<"Flux cellule "<<i<<" = "<< endl;
 988                                         cout<<Fi<<endl;
 989                                         cout<<"Flux cellule "<<j<<" = "<< endl;
 990                                         cout<<Fj<<endl;
 991                                 }
 992                         }
 993                 }
 994                 for(int i1=0;i1<_nVar;i1++)
 995                         _phi[i1]=-Fij(i1)*_inv_dxi;
 996                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
 997                 if(_verbose && _nbTimeStep%_freqSave ==0)
 998                 {
 999                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1000                         cout<<"Flux interfacial "<<i<<", "<<j<< endl;
1001                         cout<<Fij<<endl;
1002                         cout << "Contribution convection à " << i << ", -Fij(i1)*_inv_dxi= "<<endl;
1003                         for(int q=0; q<_nVar; q++)
1004                                 cout << _phi[q] << endl;
1005                         cout << endl;
1006                 }
1007                 if(!isBord){
1008                         for(int i1=0;i1<_nVar;i1++)
1009                                 _phi[i1]*=-_inv_dxj/_inv_dxi;
1010                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1011                         if(_verbose && _nbTimeStep%_freqSave ==0)
1012                         {
1013                                 cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
1014                                 for(int q=0; q<_nVar; q++)
1015                                         cout << _phi[q] << endl;
1016                                 cout << endl;
1017                         }
1018                 }
1019         }else //_nonLinearFormulation==reducedRoe)
1020         {
1021                 for(int k=0; k<_nVar; k++)
1022                         _temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
1023                 Polynoms Poly;
1024                 Poly.matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
1025                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1026
1027                 if(_verbose && _nbTimeStep%_freqSave ==0)
1028                 {
1029                         cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1030                         cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
1031                         for(int q=0; q<_nVar; q++)
1032                                 cout << _temp[q] << endl;
1033                         cout << endl;
1034                         cout << "Contribution convection à " << i << ", A^-*(Ui - Uj)*_inv_dxi= "<<endl;
1035                         for(int q=0; q<_nVar; q++)
1036                                 cout << _phi[q] << endl;
1037                         cout << endl;
1038                 }
1039
1040                 if(!isBord)
1041                 {
1042                         for(int k=0; k<_nVar; k++)
1043                                 _temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
1044                         Poly.matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
1045                         VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1046
1047                         if(_verbose && _nbTimeStep%_freqSave ==0)
1048                         {
1049                                 cout << "Contribution convection à  " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
1050                                 for(int q=0; q<_nVar; q++)
1051                                         cout << _phi[q] << endl;
1052                                 cout << endl;
1053                         }
1054                 }
1055         }
1056         if(_verbose && _nbTimeStep%_freqSave ==0)
1057         {
1058                 cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
1059                 displayMatrix(_absAroe,   _nVar,"Valeur absolue matrice de Roe");
1060                 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
1061                 displayMatrix(_AroePlus,  _nVar,"Matrice _AroePlus");
1062                 displayMatrix(_signAroe,  _nVar,"Signe de la matrice de Roe");
1063         }
1064 }
1065
1066 void ProblemFluid::addSourceTermToSecondMember
1067 (       const int i, int nbVoisinsi,
1068                 const int j, int nbVoisinsj,
1069                 bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
1070 {
1071         if(_verbose && _nbTimeStep%_freqSave ==0)
1072                 cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
1073
1074         _idm[0] = i*_nVar;
1075         for(int k=1; k<_nVar;k++)
1076                 _idm[k] = _idm[k-1] + 1;
1077         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1078         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1079         sourceVector(_Si,_Ui,_Vi,i);
1080
1081         if (_verbose && _nbTimeStep%_freqSave ==0)
1082         {
1083                 cout << "Terme source S(Ui), i= " << i<<endl;
1084                 for(int q=0; q<_nVar; q++)
1085                         cout << _Si[q] << endl;
1086                 cout << endl;
1087         }
1088         if(!isBord){
1089                 for(int k=0; k<_nVar; k++)
1090                         _idn[k] = _nVar*j + k;
1091                 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1092                 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1093                 sourceVector(_Sj,_Uj,_Vj,j);
1094         }else
1095         {
1096                 for(int k=0; k<_nVar; k++)
1097                         _idn[k] = k;
1098                 VecGetValues(_Uext, _nVar, _idn, _Uj);
1099                 consToPrim(_Uj, _Vj,_porosityj);
1100                 sourceVector(_Sj,_Uj,_Vj,i);
1101         }
1102         if (_verbose && _nbTimeStep%_freqSave ==0)
1103         {
1104                 if(!isBord)
1105                         cout << "Terme source S(Uj), j= " << j<<endl;
1106                 else
1107                         cout << "Terme source S(Uj), cellule fantôme "<<endl;
1108                 for(int q=0; q<_nVar; q++)
1109                         cout << _Sj[q] << endl;
1110                 cout << endl;
1111         }
1112
1113         //Compute pressure loss vector
1114         double K;
1115         if(_pressureLossFieldSet){
1116                 K=_pressureLossField(ij);
1117                 pressureLossVector(_pressureLossVector, K,_Ui, _Vi, _Uj, _Vj);
1118         }
1119         else{
1120                 K=0;
1121                 for(int k=0; k<_nVar;k++)
1122                         _pressureLossVector[k]=0;
1123         }
1124         if (_verbose && _nbTimeStep%_freqSave ==0)
1125         {
1126                 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
1127                 for(int k=0; k<_nVar;k++)
1128                         cout<< _pressureLossVector[k]<<", ";
1129                 cout<<endl;
1130         }
1131         //Contribution of the porosityField gradient:
1132         if(!isBord)
1133                 porosityGradientSourceVector();
1134         else{
1135                 for(int k=0; k<_nVar;k++)
1136                         _porosityGradientSourceVector[k]=0;
1137         }
1138
1139         if (_verbose && _nbTimeStep%_freqSave ==0)
1140         {
1141                 if(!isBord)
1142                         cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
1143                 else
1144                         cout<<"interface frontière i= "<<i<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<endl;
1145                 cout<<"Gradient de porosite à l'interface"<<endl;
1146                 for(int k=0; k<_nVar;k++)
1147                         cout<< _porosityGradientSourceVector[k]<<", ";
1148                 cout<<endl;
1149         }
1150         Polynoms Poly;
1151         if(!isBord){
1152                 if(_wellBalancedCorrection){
1153                         for(int k=0; k<_nVar;k++)
1154                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1155                         Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1156                         for(int k=0; k<_nVar;k++){
1157                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1158                                 _Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
1159                         }
1160                         if (_verbose && _nbTimeStep%_freqSave ==0)
1161                         {
1162                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,j), j="<<j<<endl;
1163                                 for(int q=0; q<_nVar; q++)
1164                                         cout << _Si[q] << endl;
1165                                 cout << "Contribution au terme source Sj de la cellule j= " << j<<" venant  (après décentrement) de la face (i,j), i="<<i<<endl;
1166                                 for(int q=0; q<_nVar; q++)
1167                                         cout << _Sj[q] << endl;
1168                                 cout << endl;
1169                                 cout<<"ratio surface sur volume i = "<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1170                                 cout<<"ratio surface sur volume j = "<<mesureFace/_perimeters(j)<<" perimeter = "<< _perimeters(j)<<endl;
1171                                 cout << endl;
1172                         }
1173                 }
1174                 else{
1175                         for(int k=0; k<_nVar;k++){
1176                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
1177                                 _Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
1178                         }
1179                         if (_verbose && _nbTimeStep%_freqSave ==0)
1180                         {
1181                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant  de la face (i,j), j="<<j<<endl;
1182                                 for(int q=0; q<_nVar; q++)
1183                                         cout << _Si[q] << endl;
1184                                 cout << "Contribution au terme source Sj de la cellule j = " << j<<" venant  de la face (i,j), i="<<i <<endl;
1185                                 for(int q=0; q<_nVar; q++)
1186                                         cout << _Sj[q] << endl;
1187                                 cout << endl;
1188                         }
1189                 }
1190                 _idn[0] = j;
1191                 VecSetValuesBlocked(_b, 1, _idn, _Sj, ADD_VALUES);
1192         }else{
1193                 if(_wellBalancedCorrection){
1194                         for(int k=0; k<_nVar;k++)
1195                                 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1196                         Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1197                         for(int k=0; k<_nVar;k++)
1198                                 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1199                         if (_verbose && _nbTimeStep%_freqSave ==0)
1200                         {
1201                                 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant  (après décentrement) de la face (i,bord)"<<endl;
1202                                 for(int q=0; q<_nVar; q++)
1203                                         cout << _Si[q] << endl;
1204                                 cout<<"ratio surface sur volume i ="<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1205                                 cout << endl;
1206                         }
1207                 }
1208                 else
1209                 {
1210                         for(int k=0; k<_nVar;k++)
1211                                 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
1212                         if (_verbose && _nbTimeStep%_freqSave ==0)
1213                         {
1214                                 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
1215                                 for(int q=0; q<_nVar; q++)
1216                                         cout << _Si[q] << endl;
1217                                 cout << endl;
1218                         }
1219                 }
1220         }
1221         _idm[0] = i;
1222         VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
1223
1224         if(_verbose && _nbTimeStep%_freqSave ==0 && _wellBalancedCorrection)
1225                 displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
1226 }
1227
1228 void ProblemFluid::updatePrimitives()
1229 {
1230         VecAssemblyBegin(_primitiveVars);
1231         for(int i=1; i<=_Nmailles; i++)
1232         {
1233                 _idm[0] = _nVar*( (i-1));
1234                 for(int k=1; k<_nVar; k++)
1235                         _idm[k] = _idm[k-1] + 1;
1236
1237                 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1238                 consToPrim(_Ui,_Vi,_porosityField(i-1));
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                 if(_nbPhases==2 && _Psat>-1e30){//Cas simulation flashing
1253                         double pressure;
1254                         VecGetValues(_primitiveVars, 1, _idm+1, &pressure);
1255                         _dp_over_dt(i-1)=(_Vi[1]-pressure)/_dt;//pn+1-pn
1256                 }
1257                 _idm[0] = i-1;
1258                 VecSetValuesBlocked(_primitiveVars, 1, _idm, _Vi, INSERT_VALUES);
1259         }
1260         VecAssemblyEnd(_primitiveVars);
1261
1262         if(_system)
1263         {
1264                 cout << "Nouvelles variables primitives : " << endl;
1265                 VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_WORLD);
1266                 cout << endl;
1267         }
1268 }
1269
1270 void ProblemFluid::updateConservatives()
1271 {
1272         VecAssemblyBegin(_conservativeVars);
1273         for(int i=1; i<=_Nmailles; i++)
1274         {
1275                 _idm[0] = _nVar*( (i-1));
1276                 for(int k=1; k<_nVar; k++)
1277                         _idm[k] = _idm[k-1] + 1;
1278
1279                 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1280
1281                 primToCons(_Vi,0,_Ui,0);
1282                 _idm[0] = i-1;
1283                 VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
1284
1285                 if(_verbose && _nbTimeStep%_freqSave ==0)
1286                 {
1287                         cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
1288                         cout << "Vi = ";
1289                         for(int q=0; q<_nVar; q++)
1290                                 cout << _Vi[q]  << "\t";
1291                         cout << endl;
1292                         cout << "Ui = ";
1293                         for(int q=0; q<_nVar; q++)
1294                                 cout << _Ui[q]  << "\t";
1295                         cout << endl;
1296                 }
1297         }
1298         VecAssemblyEnd(_conservativeVars);
1299
1300         if(_system)
1301         {
1302                 cout << "Nouvelles variables conservatives : " << endl;
1303                 VecView(_conservativeVars,  PETSC_VIEWER_STDOUT_WORLD);
1304                 cout << endl;
1305         }
1306 }
1307
1308 vector< complex<double> > ProblemFluid::getRacines(vector< double > pol_car){
1309         int degre_polynom = pol_car.size()-1;
1310         double tmp;
1311         vector< complex< double > > valeurs_propres (degre_polynom);
1312         vector< double > zeror(degre_polynom);
1313         vector< double > zeroi(degre_polynom);
1314         for(int j=0; j<(degre_polynom+1)/2; j++){//coefficients in order of decreasing powers for rpoly
1315                 tmp=pol_car[j];
1316                 pol_car[j] =pol_car[degre_polynom-j];
1317                 pol_car[degre_polynom-j]=tmp;
1318         }
1319
1320         //Calcul des racines du polynome
1321         roots_polynoms roots;
1322         int size_vp = roots.rpoly(&pol_car[0],degre_polynom,&zeror[0],&zeroi[0]);
1323
1324         //On ordonne les valeurs propres
1325         if(zeror[1]<zeror[0])
1326         {
1327                 tmp=zeror[0];
1328                 zeror[0]=zeror[1];
1329                 zeror[1]=tmp;
1330                 tmp=zeroi[0];
1331                 zeroi[0]=zeroi[1];
1332                 zeroi[1]=tmp;
1333         }
1334
1335         if(size_vp == degre_polynom) {
1336                 for(int ct =0; ct<degre_polynom; ct++)
1337                         valeurs_propres[ct] = complex<double> (zeror[ct],zeroi[ct]);  //vp non triviales
1338         }
1339         else {
1340                 cout << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1341                 *_runLogFile << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1342                 cout<<"Coefficients polynome caracteristique: "<<endl;
1343                 for(int ct =0; ct<degre_polynom+1; ct++)
1344                         cout<<pol_car[ct]<< " , " <<endl;
1345
1346                 *_runLogFile<<"getRacines computation failed"<<endl;
1347                 _runLogFile->close();
1348                 throw CdmathException("getRacines computation failed");
1349         }
1350
1351         return valeurs_propres;
1352 }
1353
1354 void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1355 {
1356         Polynoms Poly;
1357         int nbVp_dist=valeurs_propres_dist.size();
1358         vector< complex< double > > y (nbVp_dist,0);
1359         for( int i=0 ; i<nbVp_dist ; i++)
1360                 y[i] = Poly.abs_generalise(valeurs_propres_dist[i]);
1361         Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
1362         if(_verbose && _nbTimeStep%_freqSave ==0)
1363         {
1364                 cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1365                 for(int ct =0; ct<nbVp_dist; ct++)
1366                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1367                 cout<< endl;
1368                 cout<<"ProblemFluid::AbsMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1369                 for(int ct =0; ct<nbVp_dist; ct++)
1370                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1371                 cout<< endl;
1372                 displayMatrix(_absAroe,_nVar,"Valeur absolue de la matrice de Roe");
1373         }
1374 }
1375
1376 void ProblemFluid::SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1377 {
1378         int nbVp_dist=valeurs_propres_dist.size();
1379         vector< complex< double > > y (nbVp_dist,0);
1380         for( int i=0 ; i<nbVp_dist ; i++)
1381         {
1382                 if(valeurs_propres_dist[i].real()>0)
1383                         y[i] = 1;
1384                 else if(valeurs_propres_dist[i].real()<0)
1385                         y[i] = -1;
1386                 else
1387                         y[i] = 0;
1388         }
1389         Polynoms Poly;
1390         Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
1391         if(_verbose && _nbTimeStep%_freqSave ==0)
1392         {
1393                 cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1394                 for(int ct =0; ct<nbVp_dist; ct++)
1395                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1396                 cout<< endl;
1397                 cout<<"ProblemFluid::SigneMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1398                 for(int ct =0; ct<nbVp_dist; ct++)
1399                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1400                 cout<< endl;
1401                 displayMatrix(_signAroe,_nVar,"Signe matrice de Roe");
1402         }
1403 }
1404 void ProblemFluid::InvMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1405 {
1406         int nbVp_dist=valeurs_propres_dist.size();
1407         vector< complex< double > > y (nbVp_dist,0);
1408         Polynoms Poly;
1409         for( int i=0 ; i<nbVp_dist ; i++)
1410         {
1411                 if(Poly.module(valeurs_propres_dist[i])>_precision)
1412                         y[i] = 1./valeurs_propres_dist[i];
1413                 else
1414                         y[i] = 1./_precision;
1415         }
1416         Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
1417         if(_verbose && _nbTimeStep%_freqSave ==0)
1418         {
1419                 cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
1420                 for(int ct =0; ct<nbVp_dist; ct++)
1421                         cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
1422                 cout<< endl;
1423                 cout<<"ProblemFluid::InvMatriceRoe : Valeurs à atteindre pour l'interpolation"<<endl;
1424                 for(int ct =0; ct<nbVp_dist; ct++)
1425                         cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<")  ";
1426                 cout<< endl;
1427                 displayMatrix(_invAroe,_nVar,"Inverse matrice de Roe");
1428         }
1429 }
1430
1431 Field ProblemFluid::getConservativeField() const
1432 {
1433         if(!_initializedMemory)
1434         {
1435                 _runLogFile->close();
1436                 throw CdmathException("ProblemFluid::getConservativeField call initialize() first");
1437         }
1438         return _UU;
1439 }
1440
1441 Field ProblemFluid::getPrimitiveField() const
1442 {
1443         if(!_initializedMemory)
1444         {
1445                 _runLogFile->close();
1446                 throw CdmathException("ProblemFluid::getPrimitiveField call initialize() first");
1447         }
1448         return _VV;
1449 }
1450
1451 void ProblemFluid::terminate(){
1452
1453         delete[] _AroePlus;
1454         delete[] _Diffusion;
1455         delete[] _GravityImplicitationMatrix;
1456         delete[] _AroeMinus;
1457         delete[] _Aroe;
1458         delete[] _absAroe;
1459         delete[] _signAroe;
1460         delete[] _invAroe;
1461         delete[] _AroeImplicit;
1462         delete[] _AroeMinusImplicit;
1463         delete[] _AroePlusImplicit;
1464         delete[] _absAroeImplicit;
1465         delete[] _phi;
1466         delete[] _Jcb;
1467         delete[] _JcbDiff;
1468         delete[] _a;
1469         delete[] _primToConsJacoMat;
1470
1471         delete[] _l;
1472         delete[] _r;
1473         delete[] _Uroe;
1474         delete[] _Udiff;
1475         delete[] _temp;
1476         delete[] _externalStates;
1477         delete[] _idm;
1478         delete[] _idn;
1479         delete[] _vec_normal;
1480         delete[] _Ui;
1481         delete[] _Uj;
1482         delete[] _Vi;
1483         delete[] _Vj;
1484         if(_nonLinearFormulation==VFRoe){
1485                 delete[] _Uij;
1486                 delete[] _Vij;
1487         }
1488         delete[] _Si;
1489         delete[] _Sj;
1490         delete[] _pressureLossVector;
1491         delete[] _porosityGradientSourceVector;
1492         if(_isScaling)
1493         {
1494                 delete[] _blockDiag;
1495                 delete[] _invBlockDiag;
1496
1497                 VecDestroy(&_vecScaling);
1498                 VecDestroy(&_invVecScaling);
1499                 VecDestroy(&_bScaling);
1500         }
1501
1502         VecDestroy(&_conservativeVars);
1503         VecDestroy(&_newtonVariation);
1504         VecDestroy(&_b);
1505         VecDestroy(&_primitiveVars);
1506         VecDestroy(&_Uext);
1507         VecDestroy(&_Vext);
1508         VecDestroy(&_Uextdiff);
1509
1510         //      PCDestroy(_pc);
1511         KSPDestroy(&_ksp);
1512         for(int i=0;i<_nbPhases;i++)
1513                 delete _fluides[i];
1514 }