CoreFlows/Models/src/FiveEqsTwoFluid.cxx

   1 /*
   2  * FiveEqsTwoFluid.cxx
   3  *
   4  *  Created on: Jan 28, 2015
   5  *      Author: M. Ndjinga
   6  */
   7 #include "FiveEqsTwoFluid.hxx"
   8 #include <cstdlib>
   9 #include "StiffenedGas.hxx"
  10
  11 using namespace std;
  12
  13 extern "C" int dgeev_(char *jobvl, char *jobvr, int *n, double *
  14                 a, int *lda, double *wr, double *wi, double *vl,
  15                 int *ldvl, double *vr, int *ldvr, double *work,
  16                 int *lwork, int *info);
  17
  18
  19 FiveEqsTwoFluid::FiveEqsTwoFluid(pressureEstimate pEstimate, int dim){
  20         _Ndim=dim;
  21         _nVar=2*(_Ndim+1)+1;
  22         _nbPhases = 2;
  23         _dragCoeffs=vector<double>(2,0);
  24         _fluides.resize(2);
  25         //Ne pas utiliser la loi de Stephane Dellacherie mais la stiffened gas standard
  26         if (pEstimate==around1bar300K)//EOS at 1 bar and 373K
  27         {
  28                 cout<<"Fluid is water-Gas mixture around saturation point 1 bar and 373 K (100°C)"<<endl;
  29                 *_runLogFile<<"Fluid is water-Gas mixture around saturation point 1 bar and 373 K (100°C)"<<endl;
  30                 _fluides[0] = new StiffenedGas(1.34,1555,373,2.5e6);  //ideal gas law for Gas at pressure 1 bar and temperature 100°C: eref1=2.5e6
  31                 _fluides[1] = new StiffenedGas(958,1e5,373,4.2e5,1543,3769);  //stiffened gas law for water at pressure 1 bar and temperature 100°C: eref2=5e5
  32                 _Tsat=373;//saturation temperature at 1 bar
  33                 _hsatl=4.2e5;//water enthalpy at saturation at 1 bar
  34                 _hsatv=2.5e6;//Gas enthalpy at saturation at 1 bar
  35         }
  36         else//EOS at 155 bars and 618K
  37         {
  38                 cout<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
  39                 *_runLogFile<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
  40                 _fluides[0] = new StiffenedGas(102,1.55e7,618,2.44e6, 433,3633);  //stiffened gas law for Gas at pressure 155 bar and temperature 345°C: eref1=2.4e6
  41                 _fluides[1] = new StiffenedGas(594,1.55e7,618,1.6e6, 621,3100);  //stiffened gas law for water at pressure 155 bar and temperature 345°C: eref2=1.6e6
  42                 _Tsat=618;//saturation temperature at 155 bars
  43                 _hsatl=1.63e6;//water enthalpy at saturation at 155 bars
  44                 _hsatv=2.6e6;//Gas enthalpy at saturation at 155 bars
  45         }
  46         _latentHeat=_hsatv-_hsatl;
  47         _intPressCoeff=1.5;
  48
  49         _fileName = "SolverlabFiveEquationTwoFluid";
  50     PetscPrintf(PETSC_COMM_WORLD,"\n Five equation two-fluid problem for two phase flow\n");
  51 }
  52
  53 void FiveEqsTwoFluid::initialize()
  54 {
  55         cout<<"\n Initialising the five equation two fluid model"<<endl;
  56         *_runLogFile<<"\n Initialising the five equation two fluid model"<<endl;
  57
  58         if(static_cast<StiffenedGas*>(_fluides[0])==NULL || static_cast<StiffenedGas*>(_fluides[1])==NULL)
  59                 throw CdmathException("!!!!!!!!FiveEqsTwoFluid::initialize: both phase must have stiffened gas EOS");
  60
  61         _Uroe = new double[_nVar+1];
  62
  63         _lCon = new PetscScalar[_nVar];//should be deleted in ::terminate
  64         _rCon = new PetscScalar[_nVar];//should be deleted in ::terminate
  65         _JacoMat = new PetscScalar[_nVar*_nVar];//should be deleted in ::terminate
  66
  67         _gravite = vector<double>(_nVar,0);//Not to be confused with _GravityField3d (size _Ndim). _gravite (size _Nvar) is usefull for dealing with source term and implicitation of gravity vector
  68         for(int i=0; i<_Ndim; i++)
  69         {
  70                 _gravite[i+1]=_GravityField3d[i];
  71                 _gravite[i+1 +_Ndim+1]=_GravityField3d[i];
  72         }
  73         _GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
  74
  75         if(_saveVelocity){
  76                 _Vitesse1=Field("Gas velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
  77                 _Vitesse2=Field("Liquid velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
  78         }
  79
  80         if(_entropicCorrection)
  81                 _entropicShift=vector<double>(_nVar);
  82
  83         ProblemFluid::initialize();
  84 }
  85
  86 void FiveEqsTwoFluid::convectionState( const long &i, const long &j, const bool &IsBord){
  87         //sortie: WRoe en (alpha, p, u1, u2, T, dm1,dm2,dalpha1,dp)
  88         //entree: _conservativeVars en (rho1, rho1 u1, rho2, rho2 u2)
  89
  90         _idm[0] = _nVar*i;
  91         for(int k=1; k<_nVar; k++)
  92                 _idm[k] = _idm[k-1] + 1;
  93         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
  94
  95         _idm[0] = _nVar*j;
  96         for(int k=1; k<_nVar; k++)
  97                 _idm[k] = _idm[k-1] + 1;
  98         if(IsBord)
  99                 VecGetValues(_Uext, _nVar, _idm, _Uj);
 100         else
 101                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 102
 103         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 104         {
 105                 cout<<"Convection Left state cell " << i<< ": "<<endl;
 106                 for(int k =0; k<_nVar; k++)
 107                         cout<< _Ui[k]<<endl;
 108                 cout<<"Convection Right state cell " << j<< ": "<<endl;
 109                 for(int k =0; k<_nVar; k++)
 110                         cout<< _Uj[k]<<endl;
 111         }
 112         /*
 113         if(_Ui[0]<-(_precision) || _Uj[0]<-(_precision) || _Ui[_Ndim+1]<-(_precision) || _Uj[_Ndim+1]<-(_precision))
 114         {
 115                 cout<<"Warning: masse partielle negative!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 116                 cout<< "valeurs a gauche: "<<_Ui[0]<<", "<<_Ui[_Ndim+1]<<", valeurs a droite: "<<_Uj[0]<<", "<<_Uj[_Ndim+1]<<endl;
 117                 //        throw CdmathException(" Masse partielle negative, arret de calcul");
 118         }
 119         else
 120                 _Ui[0]=max(0.,_Ui[0]);
 121                 _Uj[0]=max(0.,_Uj[0]);
 122                 _Ui[_Ndim+1]=max(0.,_Ui[_Ndim+1]);
 123                 _Uj[_Ndim+1]=max(0.,_Uj[_Ndim+1]);
 124          */
 125
 126         PetscScalar ri1, ri2, rj1, rj2, xi, xj;
 127         //get _l and _r the primitive  states left and right of the interface
 128         _idm[0] = _nVar*i;
 129         for(int k=1; k<_nVar; k++)
 130                 _idm[k] = _idm[k-1] + 1;
 131         VecGetValues(_primitiveVars, _nVar, _idm, _l);
 132
 133         if(IsBord)
 134         {
 135                 //cout<<"_r is border"<<endl;
 136                 //consToPrim(_Uj, _r);
 137                 _idm[0] = 0;
 138                 for(int k=1; k<_nVar; k++)
 139                         _idm[k] = _idm[k-1] + 1;
 140                 VecGetValues(_Vext, _nVar, _idm, _r);
 141         }
 142         else
 143         {
 144                 _idm[0] = _nVar*j;
 145                 for(int k=1; k<_nVar; k++)
 146                         _idm[k] = _idm[k-1] + 1;
 147                 VecGetValues(_primitiveVars, _nVar, _idm, _r);
 148         }
 149         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 150         {
 151                 cout<<"_l: "<<endl;
 152                 for(int k=0;k<_nVar; k++)
 153                         cout<< _l[k]<<endl;
 154                 cout<<"_r: "<<endl;
 155                 for(int k=0;k<_nVar; k++)
 156                         cout<< _r[k]<<endl;
 157         }
 158         // _Uroe[0] = \tilde{\alpha_v} = 1 - \tilde{\alpha_l} (formula of Toumi)
 159         if(2-_l[0]-_r[0] > _precision)
 160                 _Uroe[0] = 1- 2*(1-_l[0])*(1-_r[0])/(2-_l[0]-_r[0]);
 161         else
 162                 _Uroe[0] = (_l[0]+_r[0])/2;
 163
 164         if(_l[0]+_r[0] > _precision)
 165                 _Uroe[1] = (_l[1]*_l[0]+_r[1]*_r[0])/(_l[0]+_r[0]);
 166         else
 167                 _Uroe[1] = (_l[1]*(1-_l[0])+_r[1]*(1-_r[0]))/(2-_l[0]-_r[0]);
 168
 169         ri1 = sqrt(_Ui[0]); ri2 = sqrt(_Ui[_Ndim+1]);
 170         rj1 = sqrt(_Uj[0]); rj2 = sqrt(_Uj[_Ndim+1]);
 171         for(int k=0;k<_Ndim;k++)
 172         {
 173                 xi = _Ui[k+1];
 174                 xj = _Uj[k+1];
 175                 if(ri1>_precision && rj1>_precision)
 176                         _Uroe[2+k] = (xi/ri1 + xj/rj1)/(ri1 + rj1);
 177                 else if(ri1<_precision && rj1>_precision)
 178                         _Uroe[2+k] =  xj/_Uj[0];
 179                 else if(ri1>_precision && rj1<_precision)
 180                         _Uroe[2+k] =  xi/_Ui[0];
 181                 else
 182                         _Uroe[2+k] =(_Ui[k+1+_Ndim+1]/ri2 + _Uj[k+1+_Ndim+1]/rj2)/(ri2 + rj2);
 183
 184                 xi = _Ui[k+1+_Ndim+1];
 185                 xj = _Uj[k+1+_Ndim+1];
 186                 if(ri2>_precision && rj2>_precision)
 187                         _Uroe[1+k+_Ndim+1] = (xi/ri2 + xj/rj2)/(ri2 + rj2);
 188                 else if(ri2<_precision && rj2>_precision)
 189                         _Uroe[1+k+_Ndim+1] = xj/_Uj[_Ndim+1];
 190                 else if(ri2>_precision && rj2<_precision)
 191                         _Uroe[1+k+_Ndim+1] = xi/_Ui[_Ndim+1];
 192                 else
 193                         _Uroe[1+k+_Ndim+1] = (xi/ri1 + xj/rj1)/(ri1 + rj1);
 194         }
 195         _Uroe[_nVar-1]=.5*(_l[_nVar-1]+_r[_nVar-1]);
 196
 197         //Fin du remplissage dans la fonction convectionMatrices
 198
 199         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 200         {
 201                 cout<<"Etat de Roe calcule: "<<endl;
 202                 for(int k=0;k<_nVar; k++)
 203                         cout<< _Uroe[k]<<endl;
 204         }
 205 }
 206
 207 void FiveEqsTwoFluid::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
 208         //sortie: matrices et etat Diffusion (alpha1 rho1, q1, alpha2 rho2, q2,T)
 209         _idm[0] = _nVar*i;
 210         for(int k=1; k<_nVar; k++)
 211                 _idm[k] = _idm[k-1] + 1;
 212
 213         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 214         _idm[0] = _nVar*j;
 215         for(int k=1; k<_nVar; k++)
 216                 _idm[k] = _idm[k-1] + 1;
 217
 218         if(IsBord)
 219                 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
 220         else
 221                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 222
 223         for(int k=0; k<_nVar; k++)
 224                 _Udiff[k] = (_Ui[k]+_Uj[k])/2;
 225         double q1_2,q2_2=0;
 226         for (int i = 0; i<_Ndim;i++){
 227                 q1_2+=_Udiff[        i+1]*_Udiff[        i+1];
 228                 q2_2+=_Udiff[1+_Ndim+i+1]*_Udiff[1+_Ndim+i+1];
 229         }
 230         consToPrim(_Udiff,_phi);
 231         _Udiff[_nVar-1]=_phi[_nVar-1];
 232         double alpha=_phi[0];
 233         double Tm=_phi[_nVar-1];
 234         double mu1 = _fluides[0]->getViscosity(Tm);
 235         double mu2 = _fluides[1]->getViscosity(Tm);
 236         double lambda = alpha*_fluides[0]->getConductivity(Tm)+(1-alpha)*_fluides[1]->getConductivity(Tm);
 237         double Cv1= _fluides[0]->constante("Cv");
 238         double Cv2= _fluides[1]->constante("Cv");
 239
 240         if(_timeScheme==Implicit)
 241         {
 242                 for(int i=0; i<_nVar*_nVar;i++)
 243                         _Diffusion[i] = 0;
 244                 for(int idim=1;idim<_Ndim+1;idim++)
 245                 {
 246                         if(alpha>_precision){
 247                                 _Diffusion[idim*_nVar] = alpha* mu1*_Udiff[idim]/(_Udiff[0]*_Udiff[0]);
 248                                 _Diffusion[idim*_nVar+idim] = -alpha* mu1/_Udiff[0];
 249                         }
 250                         if(1-alpha>_precision){
 251                                 _Diffusion[(idim+_Ndim+1)*_nVar] =  (1-alpha)* mu2*_Udiff[idim+_Ndim+1]/(_Udiff[_Ndim+1]*_Udiff[_Ndim+1]);
 252                                 _Diffusion[(idim+_Ndim+1)*_nVar+idim+_Ndim+1] = -(1-alpha)* mu2/_Udiff[_Ndim+1];
 253                         }
 254                 }
 255                 /*//Should correct the formula before using
 256                 int i = (_nVar-1)*_nVar;
 257                 _Diffusion[i]=lambda*(Tm/_Udiff[0]-q1_2/(2*Cv1*_Udiff[0]*_Udiff[0]*_Udiff[0]));
 258                 _Diffusion[i+1+_Ndim]=lambda*(Tm/_Udiff[1+_Ndim]-q2_2/(2*Cv2*_Udiff[1+_Ndim]*_Udiff[1+_Ndim]*_Udiff[1+_Ndim]));
 259                 for(int k=1;k<1+_Ndim;k++)
 260                 {
 261                         _Diffusion[i+k]= lambda*_Udiff[k]/(_Udiff[0]*_Udiff[0]*Cv1);
 262                         _Diffusion[i+k+1+_Ndim]= lambda*_Udiff[k+1+_Ndim]/(_Udiff[1+_Ndim]*_Udiff[+1+_Ndim]*Cv2);
 263                 }
 264                 _Diffusion[_nVar*_nVar-1]=-lambda/(_Udiff[0]*Cv1+_Udiff[1+_Ndim]*Cv2);
 265                  */
 266         }
 267 }
 268
 269 void FiveEqsTwoFluid::sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i)
 270 {
 271         double m1=Ui[0],m2=Ui[1+_Ndim],rho=m1+m2, rhoE=Ui[_nVar-1], T=Vi[_nVar-1],alpha=Vi[0], P=Vi[1];
 272         double norm_ur=0,norm_u1sq=0,norm_u2sq=0, Gamma;
 273
 274         for(int k=0; k<_Ndim; k++){
 275                 norm_ur+=(Vi[2+k]-Vi[2+k+_Ndim])*(Vi[2+k]-Vi[2+k+_Ndim]);
 276                 norm_u1sq+=Vi[2+k]*Vi[2+k];
 277                 norm_u2sq+=Vi[2+k+_Ndim]*Vi[2+k+_Ndim];
 278         }
 279         norm_ur=sqrt(norm_ur);
 280         double h=(rhoE-0.5*m1*norm_u1sq-0.5*m2*norm_u2sq+P)/rho;
 281         double u_int[_Ndim];
 282         for(int k=0; k<_Ndim; k++)
 283                 u_int[k] = 0.5*(Vi[2+k]+Vi[2+k+_Ndim]);
 284         //              u_int[k] = Vi[0]*Vi[2+k+_Ndim] + (1-Vi[0])*Vi[2+k];
 285         if(i>=0 && T>_Tsat && alpha<1-_precision)//if(i>=0 && _hsatv>h  && h>_hsatl && alpha<1-_precision)
 286                 Gamma=_heatPowerField(i)/_latentHeat;
 287         else//boundary cell, no phase change
 288                 Gamma=0;
 289
 290         for(int k=1; k<_Ndim+1; k++)
 291         {
 292                 Si[k] =_gravite[k]*m1-_dragCoeffs[0]*norm_ur*(Vi[1+k]-Vi[1+k+_Ndim]) + Gamma*u_int[k-1];//interfacial velocity= ul
 293                 Si[k+_Ndim+1] =_gravite[k+_Ndim+1]*m2+ _dragCoeffs[0]*norm_ur*(Vi[1+k]-Vi[1+k+_Ndim]) - Gamma*u_int[k-1];
 294         }
 295         if(true){//heated boiling
 296                 Si[0]=Gamma;
 297                 Si[1+_Ndim]=-Gamma;
 298         }
 299         else if (P<_Psat && alpha<1-_precision){//flash boiling
 300                 Si[0]=-_dHsatl_over_dp*_dp_over_dt(i)/_latentHeat;
 301                 Si[1+_Ndim]=_dHsatl_over_dp*_dp_over_dt(i)/_latentHeat;
 302         }
 303         else{
 304                 Si[0]=0;
 305                 Si[1+_Ndim]=0;
 306         }
 307         if(i>=0)
 308                 Si[_nVar-1]=_heatPowerField(i);
 309         else//boundary cell, no heating
 310                 Si[_nVar-1]=0;
 311         for(int k=0; k<_Ndim; k++)
 312                 Si[_nVar-1] +=_GravityField3d[k]*(Ui[1+k]+Ui[2+k+_Ndim])-_dragCoeffs[0]*norm_ur*(Vi[2+k]-Vi[2+k+_Ndim])*(Vi[2+k]-Vi[2+k+_Ndim]);
 313
 314         if(_timeScheme==Implicit)
 315         {
 316                 for(int i=0; i<_nVar*_nVar;i++)
 317                         _GravityImplicitationMatrix[i] = 0;
 318                 for(int i=0; i<_nVar/2;i++)
 319                         _GravityImplicitationMatrix[i*_nVar]=-_gravite[i];
 320                 for(int i=_nVar/2; i<_nVar;i++)
 321                         _GravityImplicitationMatrix[i*_nVar+_nVar/2]=-_gravite[i];
 322         }
 323
 324         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 325         {
 326                 cout<<"FiveEqsTwoFluid::sourceVector"<<endl;
 327                 cout<<"Ui="<<endl;
 328                 for(int k=0;k<_nVar;k++)
 329                         cout<<Ui[k]<<", ";
 330                 cout<<endl;
 331                 cout<<"Vi="<<endl;
 332                 for(int k=0;k<_nVar;k++)
 333                         cout<<Vi[k]<<", ";
 334                 cout<<endl;
 335                 cout<<"Si="<<endl;
 336                 for(int k=0;k<_nVar;k++)
 337                         cout<<Si[k]<<", ";
 338                 cout<<endl;
 339                 if(_timeScheme==Implicit)
 340                         displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
 341         }
 342 }
 343
 344 void FiveEqsTwoFluid::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
 345 {
 346         double norm_u1=0, u1_n=0, norm_u2=0, u2_n=0, m1, m2;
 347         for(int i=0;i<_Ndim;i++){
 348                 u1_n += _Uroe[1+i]      *_vec_normal[i];
 349                 u2_n += _Uroe[1+i+_Ndim]*_vec_normal[i];
 350         }
 351         pressureLoss[0]=0;
 352         pressureLoss[1+_Ndim]=0;
 353         if(u1_n>0){
 354                 for(int i=0;i<_Ndim;i++)
 355                         norm_u1 += Vi[1+i]*Vi[1+i];
 356                 norm_u1=sqrt(norm_u1);
 357                 m1=Ui[0];
 358                 for(int i=0;i<_Ndim;i++)
 359                         pressureLoss[1+i]=-K*m1*norm_u1*Vi[1+i];
 360         }
 361         else{
 362                 for(int i=0;i<_Ndim;i++)
 363                         norm_u1 += Vj[1+i]*Vj[1+i];
 364                 norm_u1=sqrt(norm_u1);
 365                 m1=Uj[0];
 366                 for(int i=0;i<_Ndim;i++)
 367                         pressureLoss[1+i]=-K*m1*norm_u1*Vj[1+i];
 368         }
 369         if(u2_n>0){
 370                 for(int i=0;i<_Ndim;i++)
 371                         norm_u2 += Vi[2+i+_Ndim]*Vi[2+i+_Ndim];
 372                 norm_u2=sqrt(norm_u2);
 373                 m2=Ui[1+_Ndim];
 374                 for(int i=0;i<_Ndim;i++)
 375                         pressureLoss[2+i+_Ndim]=-K*m2*norm_u2*Vi[2+i+_Ndim];
 376         }
 377         else{
 378                 for(int i=0;i<_Ndim;i++)
 379                         norm_u2 += Vj[2+i+_Ndim]*Vj[2+i+_Ndim];
 380                 norm_u2=sqrt(norm_u2);
 381                 m2=Uj[1+_Ndim];
 382                 for(int i=0;i<_Ndim;i++)
 383                         pressureLoss[2+i+_Ndim]=-K*m2*norm_u2*Vj[2+i+_Ndim];
 384         }
 385         pressureLoss[_nVar-1]=-K*(m1*norm_u1*norm_u1*norm_u1+m2*norm_u2*norm_u2*norm_u2);
 386
 387         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 388         {
 389                 cout<<"FiveEqsTwoFluid::pressureLossVector K= "<<K<<endl;
 390                 cout<<"Ui="<<endl;
 391                 for(int k=0;k<_nVar;k++)
 392                         cout<<Ui[k]<<", ";
 393                 cout<<endl;
 394                 cout<<"Vi="<<endl;
 395                 for(int k=0;k<_nVar;k++)
 396                         cout<<Vi[k]<<", ";
 397                 cout<<endl;
 398                 cout<<"Uj="<<endl;
 399                 for(int k=0;k<_nVar;k++)
 400                         cout<<Uj[k]<<", ";
 401                 cout<<endl;
 402                 cout<<"Vj="<<endl;
 403                 for(int k=0;k<_nVar;k++)
 404                         cout<<Vj[k]<<", ";
 405                 cout<<endl;
 406                 cout<<"pressureLoss="<<endl;
 407                 for(int k=0;k<_nVar;k++)
 408                         cout<<pressureLoss[k]<<", ";
 409                 cout<<endl;
 410         }
 411 }
 412
 413 void FiveEqsTwoFluid::porosityGradientSourceVector()
 414 {
 415         double u1_ni=0, u1_nj=0, u2_ni=0, u2_nj=0, rho1i, rho2i, rho1j, rho2j, pi=_Vi[1], pj=_Vj[1], Ti=_Vi[_nVar-1], Tj=_Vj[_nVar-1];
 416         double pij1, pij2, alphaij=_Uroe[0];
 417         for(int i=0;i<_Ndim;i++) {
 418                 u1_ni += _Vi[2+i]*_vec_normal[i];
 419                 u2_ni += _Vi[2+_Ndim+i]*_vec_normal[i];
 420                 u1_nj += _Vj[2+i]*_vec_normal[i];
 421                 u2_nj += _Vj[2+_Ndim+i]*_vec_normal[i];
 422         }
 423         _porosityGradientSourceVector[0]=0;
 424         _porosityGradientSourceVector[1+_Ndim]=0;
 425         rho1i = _fluides[0]->getDensity(pi, Ti);
 426         rho2i = _fluides[1]->getDensity(pi, Ti);
 427         rho1j = _fluides[0]->getDensity(pj, Tj);
 428         rho2j = _fluides[1]->getDensity(pj, Tj);
 429         pij1=(pi+pj)/2+rho1i*rho1j/2/(rho1i+rho1j)*(u1_ni-u1_nj)*(u1_ni-u1_nj);
 430         pij2=(pi+pj)/2+rho2i*rho2j/2/(rho2i+rho2j)*(u2_ni-u2_nj)*(u2_ni-u2_nj);
 431         for(int i=0;i<_Ndim;i++){
 432                 _porosityGradientSourceVector[1+i]      =alphaij*pij1*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
 433                 _porosityGradientSourceVector[2+_Ndim+i]=alphaij*pij2*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
 434         }
 435         _porosityGradientSourceVector[_nVar-1]=0;
 436 }
 437
 438 double FiveEqsTwoFluid::intPressDef(double alpha, double u_r2, double rho1, double rho2, double Temperature)
 439 {
 440         return  _intPressCoeff*alpha*(1-alpha)*rho1*rho2*u_r2/( alpha*rho2+(1-alpha)*rho1);
 441         +alpha*(1-alpha)*rho1*rho2*u_r2/((alpha*rho2+(1-alpha)*rho1)*(alpha*rho2+(1-alpha)*rho1)*(alpha*rho2+(1-alpha)*rho1)*(alpha*rho2+(1-alpha)*rho1))*u_r2
 442                         *(alpha*alpha*rho2-(1-alpha)*(1-alpha)*rho1)
 443                         *(alpha*alpha*rho2*rho2/(_fluides[0]->vitesseSonTemperature(Temperature,rho1)*_fluides[0]->vitesseSonTemperature(Temperature,rho1))
 444                                         -(1-alpha)*(1-alpha)*rho1*rho1/(_fluides[1]->vitesseSonTemperature(Temperature,rho2)*_fluides[1]->vitesseSonTemperature(Temperature,rho2)));
 445 }
 446
 447 Vector FiveEqsTwoFluid::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
 448         if(_verbose){
 449                 cout<<"Ucons"<<endl;
 450                 cout<<U<<endl;
 451                 cout<<"Vprim"<<endl;
 452                 cout<<V<<endl;
 453         }
 454
 455         double phim1=U(0);//phi alpha1 rho1
 456         double phim2=U(1+_Ndim);//phi alpha2 rho2
 457         Vector phiq1(_Ndim),phiq2(_Ndim);//phi alpha1 rho1 u1, phi alpha2 rho2 u2
 458         for(int i=0;i<_Ndim;i++){
 459                 phiq1(i)=U(1+i);
 460                 phiq2(i)=U(2+_Ndim+i);
 461         }
 462         double alpha=V(0);
 463         double pression=V(1);
 464         Vector vitesse1(_Ndim),vitesse2(_Ndim);
 465         for(int i=0;i<_Ndim;i++){
 466                 vitesse1(i)=V(2+i);
 467                 vitesse2(i)=V(2+_Ndim+i);
 468         }
 469         double Temperature= V(_nVar-1);
 470
 471         double vitesse1n=vitesse1*normale;
 472         double vitesse2n=vitesse2*normale;
 473         double rho1=_fluides[0]->getDensity(pression,Temperature);
 474         double rho2=_fluides[1]->getDensity(pression,Temperature);
 475         double e1_int=_fluides[0]->getInternalEnergy(Temperature,rho1);
 476         double e2_int=_fluides[1]->getInternalEnergy(Temperature,rho2);
 477
 478         double alpha_roe = _Uroe[0];//Toumi formula
 479         // interfacial pressure term (hyperbolic correction)
 480         double dpi = _Uroe[_nVar];
 481
 482         Vector F(_nVar);
 483         F(0)=phim1*vitesse1n;
 484         F(1+_Ndim)=phim2*vitesse2n;
 485         for(int i=0;i<_Ndim;i++){
 486                 F(1+i)=phim1*vitesse1n*vitesse1(i)+(alpha_roe*pression+dpi*alpha)*porosity*normale(i);
 487                 F(2+_Ndim+i)=phim2*vitesse2n*vitesse2(i)+((1-alpha_roe)*pression+dpi*(1-alpha))*normale(i)*porosity;
 488         }
 489         F(_nVar-1)=phim1*(e1_int+0.5*vitesse1*vitesse1+pression/rho1)*vitesse1n+phim2*(e2_int+0.5*vitesse2*vitesse2+pression/rho2)*vitesse2n;
 490
 491         if(_verbose){
 492                 cout<<"Flux F(U,V)"<<endl;
 493                 cout<<F<<endl;
 494         }
 495
 496         return F;
 497 }
 498
 499 void FiveEqsTwoFluid::convectionJacobianMatrix(double *V, double *n)
 500 {
 501         complex< double > tmp;
 502         // enter : V(nVar) : primitive variables
 503         //            n    : normal vector
 504         double alp = V[0];
 505         double p = V[1];
 506         double Tm = V[_nVar-1];
 507         double rho1 = _fluides[0]->getDensity(p, Tm);
 508         double rho2 = _fluides[1]->getDensity(p, Tm);
 509         double ur_2 = 0;
 510         for (int idim=0; idim<_Ndim; idim++){
 511                 ur_2 += (V[2+idim]-V[2+idim+_Ndim])*(V[2+idim]-V[2+idim+_Ndim]);
 512         }
 513         // interfacial pressure term (hyperbolic correction)
 514         double dpi1 = intPressDef(alp,ur_2, rho1, rho2,Tm);
 515         double dpi2 = dpi1;
 516
 517         /********Prepare the parameters to compute the Jacobian Matrix********/
 518         /**** coefficients a, b, c ****/
 519         double inv_a1_2,inv_a2_2,b1,c1,a2,b2,c2;
 520         double e1,e2;
 521         e1 = _fluides[0]->getInternalEnergy(V[_nVar-1],rho1);// primitive variable _l[_nVar-1]=Tm
 522         e2 = _fluides[1]->getInternalEnergy(V[_nVar-1],rho2);
 523         inv_a1_2 = static_cast<StiffenedGas*>(_fluides[0])->getDiffDensPress(e1);
 524         inv_a2_2 = static_cast<StiffenedGas*>(_fluides[1])->getDiffDensPress(e2);
 525         //double getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r);
 526         b1 = static_cast<StiffenedGas*>(_fluides[0])->getDiffDensInternalEnergy(p,e1);
 527         b2 = static_cast<StiffenedGas*>(_fluides[1])->getDiffDensInternalEnergy(p,e2);
 528         //double getJumpInternalEnergyTemperature();
 529         c1 = static_cast<StiffenedGas*>(_fluides[0])->getDiffInternalEnergyTemperature();
 530         c2 = static_cast<StiffenedGas*>(_fluides[1])->getDiffInternalEnergyTemperature();
 531         /**** coefficients eta,  varrho_2 ****/
 532         double eta[_Ndim], varrho_2;
 533         // prefix m is arithmetic mean
 534         double m1,m2, eta_n;
 535         double u1[_Ndim],u2[_Ndim], alp_u1[_Ndim],alp_u2[_Ndim];
 536         m1 = alp*rho1;
 537         m2 = (1-alp)*rho2;
 538         varrho_2 =1/((alp*rho2)*inv_a1_2+((1-alp)*rho1)*inv_a2_2);
 539         eta_n = 0.;
 540         for (int idim=0; idim<_Ndim; idim++){
 541                 u1[idim] = V[idim+2];
 542                 u2[idim] = V[_Ndim+idim+2];
 543                 alp_u1[idim] = alp*u1[idim];
 544                 alp_u2[idim] = (1-alp)*u2[idim];
 545                 eta_n += (alp_u1[idim]*(1-p/rho1*inv_a1_2)+alp_u2[idim]*(1-p/rho2*inv_a2_2))*n[idim];
 546         }
 547         double eta_varrho_2n = eta_n*varrho_2;
 548         /**** compute jump of Delta T, Delta e1, Delta e2 ****/
 549         double inv_cm = 1/(c1*m1+c2*m2);
 550         double DeltaT [_nVar], Delta_e1[_nVar], Delta_e2[_nVar];
 551         // initialize DeltaT
 552         DeltaT[0] =-e1;
 553         DeltaT[1+_Ndim] =-e2;
 554         DeltaT[_nVar-1] = 1 ;
 555         for (int idim=0; idim<_Ndim; idim++){
 556                 DeltaT[idim+1] = 0.;
 557                 DeltaT[1+_Ndim+idim+1] = 0;
 558         }
 559         for (int idim=0; idim<_Ndim; idim++){
 560                 // wrt mass gas
 561                 DeltaT[0] += 0.5*u1[idim]*u1[idim];
 562                 // wrt mass liquid
 563                 DeltaT[_Ndim+1] += 0.5*u2[idim]*u2[idim];
 564                 // wrt momentum gass
 565                 DeltaT[idim+1] += - u1[idim];//*n[idim]
 566                 // wrt momentum liquid
 567                 DeltaT[_Ndim+idim+2] += - u2[idim];//*n[idim]
 568         }
 569         // finalize  DeltaT, Delta_e1 and Delta_e2
 570         for (int i =0; i< _nVar; ++i){
 571                 DeltaT[i] = inv_cm*DeltaT[i];
 572                 Delta_e1[i] = c1*DeltaT[i];
 573                 Delta_e2[i] = c2*DeltaT[i];
 574         }
 575         /**** compute differential flux (energy equation) A5 ****/
 576
 577         double dF5[_nVar];
 578         // initialize A5
 579         for (int i=0; i<_nVar; i++){
 580                 dF5[i]=0;
 581         }
 582         dF5[0] = eta_varrho_2n*rho2; // mass gas
 583         dF5[_Ndim+1] = eta_varrho_2n*rho1; // mass liquid
 584         for (int idim=0; idim<_Ndim; idim++){
 585                 // momentum gas
 586                 dF5[idim+1]= (e1+p/rho1)*n[idim];
 587                 // momentum liquid
 588                 dF5[_Ndim+idim+2]=(e2+p/rho2)*n[idim];
 589         }
 590         // assign the value of A5 (last row of the Roe matrix)
 591         for (int idim=0; idim<_Ndim; idim++){
 592                 for (int jdim=0; jdim<_Ndim;jdim++){
 593                         dF5[0] -= u1[idim]*u1[jdim]*u1[jdim]*n[idim];// -uin * ujn^2
 594                         dF5[_Ndim+1] -= u2[idim]*u2[jdim]*u2[jdim]*n[idim];
 595                         //momentum gas
 596                         dF5[idim+1] += u1[idim]*u1[jdim]*n[jdim]+0.5*(u1[jdim]*u1[jdim])*n[idim];
 597                         //momentum liquid
 598                         dF5[_Ndim+idim+2] += u2[idim]*u2[jdim]*n[jdim]+0.5*(u2[jdim]*u2[jdim])*n[idim];
 599                 }
 600         }
 601         // final dF5
 602         double coef_e1, coef_e2;
 603         coef_e1 = - eta_varrho_2n*alp*rho2*b1;
 604         coef_e2 = - eta_varrho_2n*(1-alp)*rho1*b2;
 605         for (int idim=0; idim<_Ndim; idim++){
 606                 coef_e1 += (alp*rho1 - alp*p*b1/rho1)*u1[idim]*n[idim];
 607                 coef_e2 += ((1-alp)*rho2 - (1-alp)*p*b2/rho2)*u2[idim]*n[idim];
 608         }
 609         for (int i =0; i< _nVar; i++){
 610                 dF5[i] += coef_e1*Delta_e1[i] + coef_e2*Delta_e2[i];
 611         }
 612         /******** Construction de la matrice J *********/
 613         //lignes de masse
 614         for(int i=0; i<_nVar*_nVar;i++)
 615                 _JacoMat[i]=0.;
 616
 617         for(int idim=0; idim<_Ndim;idim++)
 618         {
 619                 _JacoMat[1+idim]=n[idim];
 620                 _JacoMat[1+idim+_Ndim+1]=0.;
 621                 _JacoMat[(_Ndim+1)*_nVar+1+idim]=0.;
 622                 _JacoMat[(_Ndim+1)*_nVar+1+idim+_Ndim+1]=n[idim];
 623         }
 624         // Roe Matrix new version
 625         for(int idim=0; idim<_Ndim;idim++)
 626                 for (int jdim=0; jdim<_Ndim;jdim++){
 627                         // momentum gas (neglect delta alpha and delta P)
 628                         _JacoMat[ (1+idim)*_nVar] += -u1[idim]*u1[jdim]*n[jdim];
 629                         _JacoMat[(1+idim)*_nVar+jdim+1] += u1[idim]*n[jdim];
 630                         _JacoMat[(1+idim)*_nVar+idim+1] += u1[jdim]*n[jdim];
 631                         // momentum liquid (neglect delta alpha and delta P)
 632                         _JacoMat[(2+_Ndim+idim)*_nVar+_Ndim+1] += -u2[idim]*u2[jdim]*n[jdim];
 633                         _JacoMat[(2+_Ndim+idim)*_nVar+_Ndim+1+jdim+1] += u2[idim]*n[jdim];
 634                         _JacoMat[(2+_Ndim+idim)*_nVar+_Ndim+1+idim+1] += u2[jdim]*n[jdim];
 635                 }
 636         // update \Delta alpha
 637         /*
 638          * (alpha *rho2*varrho_2+dpi1*(1-alpha)*inv_a2_2*varrho_2)*n[idim]
 639          */
 640         for (int idim=0; idim<_Ndim; idim++){
 641                 _JacoMat[ (1+idim)*_nVar] += dpi1*varrho_2*(1-alp)*inv_a2_2*n[idim];
 642                 _JacoMat[ (1+idim)*_nVar+_Ndim+1] += -dpi1*varrho_2*alp*inv_a1_2*n[idim];
 643                 _JacoMat[(2+_Ndim+idim)*_nVar] += - dpi2*varrho_2*(1-alp)*inv_a2_2*n[idim];
 644                 _JacoMat[(2+_Ndim+idim)*_nVar+_Ndim+1] += dpi2*varrho_2*alp*inv_a1_2*n[idim];
 645                 for  (int i=0; i<_nVar; i++){
 646                         _JacoMat[ (1+idim)*_nVar+i]+=dpi1*varrho_2*(-alp*(1-alp)*inv_a2_2*b1*Delta_e1[i]+alp*(1-alp)*inv_a1_2*b2*Delta_e2[i])*n[idim];
 647                         _JacoMat[(_Ndim+1)*_nVar+ (1+idim)*_nVar+i]+=-dpi2*varrho_2*(-alp*(1-alp)*inv_a2_2*b1*Delta_e1[i]+alp*(1-alp)*inv_a1_2*b2*Delta_e2[i])*n[idim];
 648                 }
 649         }
 650         // update \Delta P
 651         for (int idim=0; idim<_Ndim; idim++){
 652                 _JacoMat[ (1+idim)*_nVar] += alp*varrho_2*rho2*n[idim];
 653                 _JacoMat[ (1+idim)*_nVar+_Ndim+1] +=alp* varrho_2*rho1*n[idim];
 654                 _JacoMat[(2+_Ndim+idim)*_nVar] +=  (1-alp)*varrho_2*rho2*n[idim];
 655                 _JacoMat[(2+_Ndim+idim)*_nVar+_Ndim+1] += (1-alp)* varrho_2*rho1*n[idim];
 656                 for  (int i=0; i<_nVar; i++){
 657                         _JacoMat[ (1+idim)*_nVar+i]+=alp*varrho_2*(-alp*rho2*b1*Delta_e1[i] -(1-alp)*rho1*b2*Delta_e2[i])*n[idim];
 658                         _JacoMat[(_Ndim+1)*_nVar+ (1+idim)*_nVar+i]+=(1-alp)*varrho_2*(-alp*rho2*b1*Delta_e1[i] -(1-alp)*rho1*b2*Delta_e2[i])*n[idim];
 659                 }
 660         }
 661         // last row (total energy)
 662         for (int i=0; i<_nVar; i++){
 663                 _JacoMat[(2*_Ndim+2)*_nVar +i] = dF5[i];
 664         }
 665 }
 666
 667 void FiveEqsTwoFluid::convectionMatrices()
 668 {
 669         //entree: URoe = alpha, p, u1, u2, T + ajout dpi
 670         //sortie: matrices Roe+  et Roe- +Roe si schéma centre
 671
 672         if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)
 673                 throw CdmathException("Implicitation with primitive variables not yet available for FiveEqsTwoFluid model");
 674
 675         /*Definitions */
 676         complex< double > tmp;
 677         double  u_r2 = 0, u1_n=0, u2_n=0;
 678         // u1_n = u1.n;u2_n = u2.n (scalar product)
 679         for(int i=0;i<_Ndim;i++)
 680         {
 681                 //u_r2 += (_Uroe[2+i]-_Uroe[1+i+1+_Ndim])*(_Uroe[2+i]-_Uroe[1+i+1+_Ndim]);
 682                 u_r2 += 0.5*((_l[2+i]-_l[2+i+_Ndim])*(_l[2+i]-_l[2+i+_Ndim])+(_r[2+i]-_r[2+i+_Ndim])*(_r[2+i]-_r[2+i+_Ndim])); //Kieu
 683                 u1_n += _Uroe[2+i]*_vec_normal[i];
 684                 u2_n += _Uroe[1+i+1+_Ndim]*_vec_normal[i];
 685         }
 686
 687         //double alpha = (_l[0]+_r[0])*0.5;//Kieu formula
 688         double alpha = _Uroe[0];//Toumi formula
 689         //double p = _Uroe[1];//Toumi pressure
 690
 691         /***********Calcul des valeurs propres ********/
 692
 693         // ********Prepare the parameters to compute the Roe Matrix******** //
 694         // **** coefficients eta,  varrho_2 **** //
 695         double eta[_Ndim], varrho_2;
 696         double rho1_l =  _fluides[0]->getDensity(_l[1], _l[_Ndim*2+2]);//(p,T)_Ndim*2+2
 697         double rho2_l =  _fluides[1]->getDensity(_l[1], _l[_Ndim*2+2]);
 698         double rho1_r =  _fluides[0]->getDensity(_r[1], _r[_Ndim*2+2]);
 699         double rho2_r =  _fluides[1]->getDensity(_r[1], _r[_Ndim*2+2]);
 700         // **** coefficients a, b, c **** //
 701         double inv_a1_2,inv_a2_2,b1,c1,a2,b2,c2;
 702         double e1_l,e1_r,e2_l,e2_r;
 703         e1_l = _fluides[0]->getInternalEnergy(_l[_nVar-1],rho1_l);// primitive variable _l[_nVar-1]=Tm
 704         e2_l = _fluides[1]->getInternalEnergy(_l[_nVar-1],rho2_l);
 705         e1_r = _fluides[0]->getInternalEnergy(_r[_nVar-1],rho1_r);
 706         e2_r = _fluides[1]->getInternalEnergy(_r[_nVar-1],rho2_r);
 707         inv_a1_2 = static_cast<StiffenedGas*>(_fluides[0])->getJumpDensPress(e1_l,e1_r);
 708         inv_a2_2 = static_cast<StiffenedGas*>(_fluides[1])->getJumpDensPress(e2_l,e2_r);
 709         //double getJumpDensInternalEnergy(const double p_l,const double p_r,const double e_l,const double e_r);
 710         b1 = static_cast<StiffenedGas*>(_fluides[0])->getJumpDensInternalEnergy(_l[1],_r[1],e1_l,e1_r);
 711         b2 = static_cast<StiffenedGas*>(_fluides[1])->getJumpDensInternalEnergy(_l[1],_r[1],e2_l,e2_r);
 712         //double getJumpInternalEnergyTemperature();
 713         c1 = static_cast<StiffenedGas*>(_fluides[0])->getJumpInternalEnergyTemperature();
 714         c2 = static_cast<StiffenedGas*>(_fluides[1])->getJumpInternalEnergyTemperature();
 715
 716         // prefix m is arithmetic mean
 717         double m_alp1,m_rho1,m_rho2,m_P,m_e1,m_e2,m_m1,m_m2, eta_n;
 718         double m_u1[_Ndim],m_u2[_Ndim], m_alp_u1[_Ndim],m_alp_u2[_Ndim];
 719         double u1_l[_Ndim], u2_l[_Ndim],u1_r[_Ndim],u2_r[_Ndim];
 720         double u1l_2 = 0, u1r_2 = 0, u2l_2 = 0, u2r_2 = 0;
 721         m_alp1 =  0.5*(_l[0]+_r[0]);
 722         m_rho1 = 0.5*(rho1_l+rho1_r);
 723         m_rho2 = 0.5*(rho2_l+rho2_r);
 724         m_P = 0.5*(_l[1]+_r[1]);
 725         m_e1 = 0.5*(e1_l+e1_r);
 726         m_e2 = 0.5*(e2_l+e2_r);
 727         m_m1 = 0.5*(_l[0]*rho1_l+_r[0]*rho1_r);
 728         m_m2 = 0.5*((1-_l[0])*rho2_l+(1-_r[0])*rho2_r);
 729         varrho_2 =1/((m_alp1*m_rho2)*inv_a1_2+((1-m_alp1)*m_rho1)*inv_a2_2);
 730         eta_n = 0.;
 731
 732         for (int idim=0; idim<_Ndim; idim++){
 733                 u1_l[idim] = _l[idim+2];
 734                 u1_r[idim] = _r[idim+2];
 735                 u2_l[idim] = _l[_Ndim+idim+2];
 736                 u2_r[idim] = _r[_Ndim+idim+2];
 737                 m_u1[idim] = 0.5*(u1_l[idim] + u1_r[idim]);
 738                 m_u2[idim] = 0.5*(u2_l[idim] + u2_r[idim]);
 739                 m_alp_u1[idim] = 0.5*(_l[0]*u1_l[idim]+_r[0]*u1_r[idim]);
 740                 m_alp_u2[idim] = 0.5*((1-_l[0])*u2_l[idim]+(1-_r[0])*u2_r[idim]);
 741                 eta_n += (m_alp_u1[idim]*(1-m_P/m_rho1*inv_a1_2)+m_alp_u2[idim]*(1-m_P/m_rho2*inv_a2_2))*_vec_normal[idim];
 742         }
 743
 744         double eta_varrho_2n = eta_n*varrho_2;
 745         //  **** compute jump of Delta T, Delta e1, Delta e2 **** //
 746         for  (int idim=0; idim<_Ndim; idim++){
 747                 u1l_2 += u1_l[idim]*u1_l[idim];
 748                 u1r_2 += u1_r[idim]*u1_r[idim];
 749                 u2l_2 += u2_l[idim]*u2_l[idim];
 750                 u2r_2 += u2_r[idim]*u2_r[idim];
 751         }
 752         double inv_m_cm = 1/(c1*m_m1+c2*m_m2);
 753         double DeltaT [_nVar], Delta_e1[_nVar], Delta_e2[_nVar];
 754         // initialize DeltaT
 755         for (int i=0; i<_nVar; i++){
 756                 DeltaT[i] = 0.;
 757         }
 758         DeltaT[0] += -m_e1;
 759         DeltaT[1+_Ndim] += -m_e2;
 760         DeltaT[_nVar-1] += 1 ;
 761         for (int idim=0; idim<_Ndim; idim++){
 762                 // wrt mass gas
 763                 DeltaT[0] += 0.5*_l[idim+2]     *_r[idim+2];//0.5*\tilde{u_g}^2
 764                 // wrt mass liquid
 765                 DeltaT[_Ndim+1] += 0.5*_l[_Ndim+idim+2] *_r[_Ndim+idim+2];
 766                 // wrt momentum gas
 767                 DeltaT[idim+1] += - m_u1[idim];
 768                 // wrt momentum liquid
 769                 DeltaT[_Ndim+idim+2] += - m_u2[idim];
 770         }
 771
 772         // finalize  DeltaT, Delta_e1 and Delta_e2
 773         for (int i =0; i< _nVar; ++i){
 774                 DeltaT[i] = inv_m_cm*DeltaT[i];
 775                 Delta_e1[i] = c1*DeltaT[i];
 776                 Delta_e2[i] = c2*DeltaT[i];
 777         }
 778
 779         // *** compute jump flux (energy equation) A5 *** //
 780
 781         double A5[_nVar];
 782         // initialize A5
 783         for (int i=0; i<_nVar; i++){
 784                 A5[i]=0;
 785         }
 786         A5[0] = eta_varrho_2n*m_rho2; // mass gas
 787         A5[_Ndim+1] = eta_varrho_2n*m_rho1; // mass liquid
 788         for (int idim=0; idim<_Ndim; idim++){
 789                 // momentum gas
 790                 A5[idim+1] = (m_e1+m_P/m_rho1)*_vec_normal[idim];
 791                 // momentum liquid
 792                 A5[_Ndim+idim+2] = (m_e2+m_P/m_rho2)*_vec_normal[idim];
 793         }
 794         // assign the value of A5 (last row of the Roe matrix)
 795         for (int idim=0; idim<_Ndim; idim++){
 796                 for (int jdim=0; jdim<_Ndim;jdim++){
 797                         A5[0] += 0.5*(0.5*(_l[idim+2]*_l[jdim+2]*_l[jdim+2]+_r[idim+2]*_r[jdim+2]*_r[jdim+2])*_vec_normal[idim]);// m_(uin*uj^2)
 798                         A5[0] -= m_u1[idim]*m_u1[jdim]*m_u1[jdim]*_vec_normal[idim];// -m_uin * m_ujn^2
 799                         A5[0] -= 0.5*(m_u1[idim]*_vec_normal[idim]*0.5*(_l[jdim+2]*_l[jdim+2]+_r[jdim+2]*_r[jdim+2]));//-0.5*m_uin*m_uj^2
 800                         A5[_Ndim+1]+= 0.5*(0.5*(_l[_Ndim+idim+2]*_l[_Ndim+jdim+2]*_l[_Ndim+jdim+2]+_r[_Ndim+idim+2]*_r[_Ndim+jdim+2]*_r[_Ndim+jdim+2])*_vec_normal[idim]);// m_(uin*uj^2)
 801                         A5[_Ndim+1] -= m_u2[idim]*m_u2[jdim]*m_u2[jdim]*_vec_normal[idim];// -m_uin * m_ujn^2
 802                         A5[_Ndim+1] -= 0.5*(m_u2[idim]*0.5*(_l[_Ndim+jdim+2]*_l[_Ndim+jdim+2]+_r[_Ndim+jdim+2]*_r[_Ndim+jdim+2]))*_vec_normal[idim];//-0.5*m_uin*m_uj^2
 803                         //momentum gas
 804                         A5[idim+1] += m_u1[jdim]*_vec_normal[jdim]*m_u1[idim]+0.5*_vec_normal[idim]*0.5*(_l[jdim+2]*_l[jdim+2]+_r[jdim+2]*_r[jdim+2]);
 805                         //momentum liquid
 806                         A5[_Ndim+idim+2] += m_u2[jdim]*_vec_normal[jdim]*m_u2[idim]+0.5*_vec_normal[idim]*0.5*(_l[_Ndim+jdim+2]*_l[_Ndim+jdim+2]+_r[_Ndim+jdim+2]*_r[_Ndim+jdim+2]);
 807                 }
 808         }
 809
 810         // final A5
 811         double coef_e1, coef_e2;
 812         coef_e1 = - eta_varrho_2n*m_alp1*m_rho2*b1;
 813         coef_e2 = - eta_varrho_2n*(1-m_alp1)*m_rho1*b2;
 814         for (int idim=0; idim<_Ndim; idim++){
 815                 coef_e1 += (0.5*(_l[0]*rho1_l*_l[idim+2]+_r[0]*rho1_r*_r[idim+2]) - m_alp_u1[idim]*m_P*b1/m_rho1)*_vec_normal[idim];
 816                 coef_e2 += (0.5*((1-_l[0])*rho2_l*_l[_Ndim+idim+2]+(1-_r[0])*rho2_r*_r[_Ndim+idim+2])-m_alp_u2[idim]*m_P*b2/m_rho2)*_vec_normal[idim];
 817         }
 818         for (int i =0; i< _nVar; i++){
 819                 A5[i] += coef_e1*Delta_e1[i] + coef_e2*Delta_e2[i];
 820         }
 821         // ******* Construction de la matrice de Roe ******** //
 822         // interfacial pressure correction
 823         double T=_Uroe[_nVar-1];
 824         double dpi1 = intPressDef(alpha, u_r2, m_rho1,m_rho2,T);
 825         double dpi2 = dpi1;
 826         //saving dpi value for flux calculation later
 827         _Uroe[_nVar]=dpi1 ;
 828
 829         //lignes de masse
 830         for(int i=0; i<_nVar*_nVar;i++)
 831                 _Aroe[i]=0;
 832         //      alpha = 0.; dpi1 = 0.; dpi2 = 0.;
 833         for(int idim=0; idim<_Ndim;idim++)
 834         {
 835                 _Aroe[1+idim]=_vec_normal[idim];
 836                 _Aroe[1+idim+_Ndim+1]=0;
 837                 _Aroe[(_Ndim+1)*_nVar+1+idim]=0;
 838                 _Aroe[(_Ndim+1)*_nVar+1+idim+_Ndim+1]=_vec_normal[idim];
 839         }
 840
 841         // Take into account the convection term in the momentum eqts
 842         for(int idim=0; idim<_Ndim;idim++)
 843                 for (int jdim=0; jdim<_Ndim;jdim++){
 844                         // momentum gas (neglect delta alpha and delta P)
 845                         _Aroe[ (1+idim)*_nVar] += (0.5*(_l[2+idim]*_l[2+jdim]+_r[2+idim]*_r[2+jdim])-2*m_u1[idim]*m_u1[jdim])*_vec_normal[jdim];
 846                         _Aroe[(1+idim)*_nVar+jdim+1] += m_u1[idim]*_vec_normal[jdim];
 847                         _Aroe[(1+idim)*_nVar+idim+1] += m_u1[jdim]*_vec_normal[jdim];
 848                         // momentum liquid (neglect delta alpha and delta P)
 849                         _Aroe[(2+_Ndim+idim)*_nVar+_Ndim+1] += (0.5*(_l[_Ndim+2+idim]*_l[_Ndim+2+jdim]+_r[_Ndim+2+idim]*_r[_Ndim+2+jdim])-2*m_u2[idim]*m_u2[jdim])*_vec_normal[jdim];
 850                         _Aroe[(2+_Ndim+idim)*_nVar+_Ndim+1+jdim+1] += m_u2[idim]*_vec_normal[jdim];
 851                         _Aroe[(2+_Ndim+idim)*_nVar+_Ndim+1+idim+1] += m_u2[jdim]*_vec_normal[jdim];
 852                 }
 853
 854         // update \Delta alpha
 855         for (int idim=0; idim<_Ndim; idim++){
 856                 _Aroe[ (1+idim)*_nVar] += dpi1*varrho_2*(1-m_alp1)*inv_a2_2*_vec_normal[idim];
 857                 _Aroe[ (1+idim)*_nVar+_Ndim+1] += -dpi1*varrho_2*(m_alp1)*inv_a1_2*_vec_normal[idim];
 858                 _Aroe[(2+_Ndim+idim)*_nVar] += - dpi2*varrho_2*(1-m_alp1)*inv_a2_2*_vec_normal[idim];
 859                 _Aroe[(2+_Ndim+idim)*_nVar+_Ndim+1] += dpi2*varrho_2*(m_alp1)*inv_a1_2*_vec_normal[idim];
 860                 for  (int i=0; i<_nVar; i++){
 861                         _Aroe[ (1+idim)*_nVar+i]+=dpi1*varrho_2*(-m_alp1*(1-m_alp1)*inv_a2_2*b1*Delta_e1[i]+m_alp1*(1-m_alp1)*inv_a1_2*b2*Delta_e2[i])*_vec_normal[idim];
 862                         _Aroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar+i]+=-dpi2*varrho_2*(-m_alp1*(1-m_alp1)*inv_a2_2*b1*Delta_e1[i]+m_alp1*(1-m_alp1)*inv_a1_2*b2*Delta_e2[i])*_vec_normal[idim];
 863                 }
 864         }
 865         // update \Delta P
 866         for (int idim=0; idim<_Ndim; idim++){
 867                 _Aroe[ (1+idim)*_nVar] += alpha*varrho_2*m_rho2*_vec_normal[idim];
 868                 _Aroe[ (1+idim)*_nVar+_Ndim+1] += alpha* varrho_2*m_rho1*_vec_normal[idim];
 869                 _Aroe[(2+_Ndim+idim)*_nVar] += (1-alpha)*varrho_2*m_rho2*_vec_normal[idim];
 870                 _Aroe[(2+_Ndim+idim)*_nVar+_Ndim+1] += (1-alpha)* varrho_2*m_rho1*_vec_normal[idim];
 871                 for  (int i=0; i<_nVar; i++){
 872                         _Aroe[ (1+idim)*_nVar+i] += alpha*varrho_2*(-m_alp1*m_rho2*b1*Delta_e1[i] -(1-m_alp1)*m_rho1*b2*Delta_e2[i])*_vec_normal[idim];
 873                         _Aroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar+i] += (1-alpha)*varrho_2*(-m_alp1*m_rho2*b1*Delta_e1[i] -(1-m_alp1)*m_rho1*b2*Delta_e2[i])*_vec_normal[idim];
 874                 }
 875         }
 876
 877         // last row (total energy)
 878         for (int i=0; i<_nVar; i++){
 879                 _Aroe[(2*_Ndim+2)*_nVar +i] += A5[i];
 880         }
 881
 882         int LDA = _nVar;
 883         int LDVL = _nVar;
 884         int LDVR = _nVar;
 885         int LWORK = 50*_nVar;
 886         char jobvl[]="N", jobvr[]="N";
 887         double WORK[LWORK], Aroe[_nVar*_nVar],egvaReal[_nVar],egvaImag[_nVar],
 888         egVectorL[_nVar*_nVar],egVectorR[_nVar*_nVar];
 889         int info=0;
 890         double maxvploc=0;
 891         std::vector< std::complex<double> > valeurs_propres_dist;
 892         int taille_vp ;
 893
 894         for (int i=0; i<_nVar*_nVar; i++)
 895                 Aroe[i] = _Aroe[i];
 896
 897         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
 898         {
 899                 cout<<endl<<"Matrice de Roe"<<endl;
 900                 for(int i=0; i<_nVar;i++)
 901                 {
 902                         for(int j=0; j<_nVar;j++)
 903                                 cout << _Aroe[i*_nVar+j]<< " , ";
 904                         cout<<endl;
 905                 }
 906         }
 907         /******** Compute the eigenvalues and eigenvectors of Roe Matrix (using lapack)*********/
 908         dgeev_(jobvl, jobvr,  &_nVar,
 909                         Aroe,&LDA,egvaReal,egvaImag, egVectorL,
 910                         &LDVL,egVectorR,
 911                         &LDVR, WORK,&LWORK,
 912                         &info);
 913         if(info <0){
 914                 cout<<"FiveEqsTwoFluid::convectionMatrices: error dgeev_ : argument "<<-info<<" invalid"<<endl;
 915                 *_runLogFile<<"FiveEqsTwoFluid::convectionMatrices: error dgeev_ : argument "<<-info<<" invalid"<<endl;
 916                 throw CdmathException("FiveEqsTwoFluid::convectionMatrices: dgeev_ unsuccessful computation of the eigenvalues ");
 917         }
 918         else if(info >0)
 919         {
 920                 cout<<"Warning FiveEqsTwoFluid::convectionMatrices: dgeev_ did not compute all the eigenvalues, trying Rusanov scheme "<<endl;
 921                 cout<<"Converged eigenvalues are ";
 922                 for(int i=info; i<_nVar; i++)
 923                         cout<<"("<< egvaReal[i]<<","<< egvaImag[i]<<"), ";
 924                 cout<<endl;
 925
 926                 _maxvploc=0;
 927                 for(int i =info; i<_nVar; i++)
 928                         if (fabs(egvaReal[i])>_maxvploc)
 929                                 _maxvploc=fabs(egvaReal[i]);
 930                 if(_maxvploc>_maxvp)
 931                         _maxvp=_maxvploc;
 932
 933                 valeurs_propres_dist=std::vector< std::complex<double> > (1,maxvploc);
 934                 taille_vp =1;
 935         }
 936         else{
 937                 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
 938                 {
 939                         for(int i=0; i<_nVar; i++)
 940                                 cout<<" Vp real part " << egvaReal[i]<<", Imaginary part " << egvaImag[i]<<endl;
 941                 }
 942
 943                 std::vector< std::complex<double> > valeurs_propres(_nVar);
 944
 945                 bool complexEigenvalues = false;
 946                 for(int j=0; j<_nVar; j++){
 947                         if (max(_l[0],_r[0])<_precision && abs(egvaImag[j])<_precision )// Kieu test Monophase
 948                                 egvaImag[j] = 0.;
 949                         else
 950                                 if (egvaImag[j] >_precision){// Kieu
 951                                         complexEigenvalues = true;
 952                                 }
 953                         if (abs(_l[0]-_r[0])<_precision*_precision && fabs(egvaImag[j])<_precision)// Kieu interfacial pressure
 954                                 egvaImag[j] = 0.;
 955                         valeurs_propres[j] = complex<double>(egvaReal[j],egvaImag[j]);
 956                 }
 957
 958                 taille_vp =Polynoms::new_tri_selectif(valeurs_propres,valeurs_propres.size(),_precision);
 959
 960                 valeurs_propres_dist=vector< complex< double > >(taille_vp);
 961                 for( int i=0 ; i<taille_vp ; i++)
 962                         valeurs_propres_dist[i] = valeurs_propres[i];
 963                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 964                 {
 965                         cout<<" Vp apres tri " << valeurs_propres_dist.size()<<endl;
 966                         for(int ct =0; ct<taille_vp; ct++)
 967                                 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<")  ";
 968                         cout<< endl;
 969                 }
 970
 971                 for(int i =0; i<taille_vp; i++){
 972                         if (fabs(valeurs_propres_dist[i].real())>maxvploc)
 973                                 maxvploc=fabs(valeurs_propres_dist[i].real());
 974                         if(maxvploc>_maxvp)
 975                                 _maxvp=maxvploc;
 976                 }
 977
 978                 int existVpCplx = 0,pos_conj;
 979                 double vp_imag_iter;
 980                 for (int ct=0; ct<taille_vp; ct++) {
 981                         vp_imag_iter = valeurs_propres_dist[ct].imag();
 982                         if ( fabs(vp_imag_iter) > 100*_precision ) {
 983                                 existVpCplx +=1;
 984                                 if ( _part_imag_max < fabs(vp_imag_iter))
 985                                         _part_imag_max = fabs(vp_imag_iter);
 986                                 //On cherhe le conjugue
 987                                 pos_conj = ct+1;
 988                                 while(pos_conj<taille_vp && fabs(valeurs_propres_dist[pos_conj].imag()+vp_imag_iter)>_precision)
 989                                         pos_conj++;
 990                                 if(pos_conj!=ct+1 && pos_conj<taille_vp )
 991                                 {
 992                                         tmp=valeurs_propres_dist[ct+1];
 993                                         valeurs_propres_dist[ct+1]=valeurs_propres_dist[pos_conj];
 994                                         valeurs_propres_dist[pos_conj] = tmp;
 995                                         ct++;
 996                                 }
 997                         }
 998                 }
 999                 if (existVpCplx >0)
1000                         _nbVpCplx +=1;
1001         }
1002         /******* Construction des matrices de decentrement *******/
1003         if(_spaceScheme == centered )
1004         {
1005                 if(_entropicCorrection)
1006                         throw CdmathException("FiveEqsTwoFluid::convectionMatrices: entropic scheme not available for centered scheme");
1007                 for(int i=0; i<_nVar*_nVar;i++)
1008                         _absAroe[i]=0;
1009                 //  if(alpha<_precision || alpha>1-_precision)//rusanov
1010                 //          for(int i=0; i<_nVar;i++)
1011                 //            _absAroe[i*_nVar+i]=maxvploc;
1012         }
1013         if( _spaceScheme ==staggered){//To do: study entropic correction for staggered
1014                 if(_entropicCorrection)//To do: study entropic correction for staggered
1015                         throw CdmathException("FiveEqsTwoFluid::convectionMatrices: entropic scheme not yet available for staggered scheme");
1016                 /******** Construction de la matrice de decentrement staggered *********/
1017                 /***** Compute eta_n **************/
1018                 eta_n = 0.;
1019                 for (int idim=0; idim<_Ndim; idim++){
1020                         eta_n += (m_alp_u1[idim]*(-1-m_P/m_rho1*inv_a1_2)+m_alp_u2[idim]*(-1-m_P/m_rho2*inv_a2_2))*_vec_normal[idim];
1021                 }
1022                 double eta_varrho_2n = eta_n*varrho_2;
1023                 /**** compute jump flux (energy equation) A5 ****/
1024
1025                 double A5[_nVar];
1026                 // initialize A5
1027                 for (int i=0; i<_nVar; i++){
1028                         A5[i]=0;
1029                 }
1030                 A5[0] = eta_varrho_2n*m_rho2; // mass gas
1031                 A5[_Ndim+1] = eta_varrho_2n*m_rho1; // mass liquid
1032                 // assign the value of A5 (last row of the Roe matrix)
1033                 for (int idim=0; idim<_Ndim; idim++){
1034                         for (int jdim=0; jdim<_Ndim;jdim++){
1035                                 A5[0] += 0.5*(0.5*(_l[idim+2]*_l[jdim+2]*_l[jdim+2]+_r[idim+2]*_r[jdim+2]*_r[jdim+2])*_vec_normal[idim]);// m_(uin*uj^2)
1036                                 A5[0] -= m_u1[idim]*m_u1[jdim]*m_u1[jdim]*_vec_normal[idim];// -m_uin * m_ujn^2
1037                                 A5[0] -= 0.5*(m_u1[idim]*_vec_normal[idim]*0.5*(_l[jdim+2]*_l[jdim+2]+_r[jdim+2]*_r[jdim+2]));//-0.5*m_uin*m_uj^2
1038                                 A5[_Ndim+1]+= 0.5*(0.5*(_l[_Ndim+idim+2]*_l[_Ndim+jdim+2]*_l[_Ndim+jdim+2]+_r[_Ndim+idim+2]*_r[_Ndim+jdim+2]*_r[_Ndim+jdim+2])*_vec_normal[idim]);// m_(uin*uj^2)
1039                                 A5[_Ndim+1] -= m_u2[idim]*m_u2[jdim]*m_u2[jdim]*_vec_normal[idim];// -m_uin * m_ujn^2
1040                                 A5[_Ndim+1] -= 0.5*(m_u2[idim]*_vec_normal[idim]*0.5*(_l[_Ndim+jdim+2]*_l[_Ndim+jdim+2]+_r[_Ndim+jdim+2]*_r[_Ndim+jdim+2]));//-0.5*m_uin*m_uj^2
1041                         }
1042                 }
1043                 // final A5
1044                 double coef_e1, coef_e2;
1045                 coef_e1 = - eta_varrho_2n*m_alp1*m_rho2*b1;
1046                 coef_e2 = - eta_varrho_2n*(1-m_alp1)*m_rho1*b2;
1047                 for (int idim=0; idim<_Ndim; idim++){
1048                         coef_e1 += (0.5*(_l[0]*rho1_l*_l[idim+2]+_r[0]*rho1_r*_r[idim+2]) - m_alp_u1[idim]*m_P*b1/m_rho1)*_vec_normal[idim];
1049                         coef_e2 += (0.5*((1-_l[0])*rho2_l*_l[_Ndim+idim+2]+(1-_r[0])*rho2_r*_r[_Ndim+idim+2])-m_alp_u2[idim]*m_P*b2/m_rho2)*_vec_normal[idim];
1050                 }
1051                 for (int i =0; i< _nVar; i++){
1052                         A5[i] += coef_e1*Delta_e1[i] + coef_e2*Delta_e2[i];
1053                 }
1054                 /** Début remplissage matrice décentrement staggered **/
1055                 //lignes de masse
1056                 for(int i=0; i<_nVar*_nVar;i++)
1057                         _absAroe[i]=0;
1058
1059                 //Mass lines
1060                 for(int idim=0; idim<_Ndim;idim++)
1061                 {
1062                         _absAroe[1+idim]=_vec_normal[idim];
1063                         _absAroe[1+idim+_Ndim+1]=0;
1064                         _absAroe[(_Ndim+1)*_nVar+1+idim]=0;
1065                         _absAroe[(_Ndim+1)*_nVar+1+idim+_Ndim+1]=_vec_normal[idim];
1066                 }
1067                 //Contribution of convection (rho u\times u) in the momentum equations
1068                 for(int idim=0; idim<_Ndim;idim++)
1069                         for (int jdim=0; jdim<_Ndim;jdim++){
1070                                 // momentum gas (neglect delta alpha and delta P)
1071                                 _absAroe[ (1+idim)*_nVar] += (0.5*(_l[2+idim]*_l[2+jdim]+_r[2+idim]*_r[2+jdim])-2*m_u1[idim]*m_u1[jdim])*_vec_normal[jdim];
1072                                 _absAroe[(1+idim)*_nVar+jdim+1] += m_u1[idim]*_vec_normal[jdim];
1073                                 _absAroe[(1+idim)*_nVar+idim+1] += m_u1[jdim]*_vec_normal[jdim];
1074                                 // momentum liquid (neglect delta alpha and delta P)
1075                                 _absAroe[(2+_Ndim+idim)*_nVar+_Ndim+1] += (0.5*(_l[_Ndim+2+idim]*_l[_Ndim+2+jdim]+_r[_Ndim+2+idim]*_r[_Ndim+2+jdim])-2*m_u2[idim]*m_u2[jdim])*_vec_normal[jdim];
1076                                 _absAroe[(2+_Ndim+idim)*_nVar+_Ndim+1+jdim+1] += m_u2[idim]*_vec_normal[jdim];
1077                                 _absAroe[(2+_Ndim+idim)*_nVar+_Ndim+1+idim+1] += m_u2[jdim]*_vec_normal[jdim];
1078                         }
1079                 // contribution of interfacial pressure in momentum equation
1080                 /*
1081                  * (alpha *rho2*varrho_2+dpi1*(1-alpha)*inv_a2_2*varrho_2)*_vec_normal[idim]
1082                  */
1083                 for (int idim=0; idim<_Ndim; idim++){
1084                         _absAroe[ (1+idim)*_nVar] += dpi1*varrho_2*(1-m_alp1)*inv_a2_2*_vec_normal[idim];
1085                         _absAroe[ (1+idim)*_nVar+_Ndim+1] += -dpi1*varrho_2*(m_alp1)*inv_a1_2*_vec_normal[idim];
1086                         _absAroe[(2+_Ndim+idim)*_nVar] += - dpi2*varrho_2*(1-m_alp1)*inv_a2_2*_vec_normal[idim];
1087                         _absAroe[(2+_Ndim+idim)*_nVar+_Ndim+1] += dpi2*varrho_2*(m_alp1)*inv_a1_2*_vec_normal[idim];
1088                         for  (int i=0; i<_nVar; i++){
1089                                 _absAroe[ (1+idim)*_nVar+i]+=dpi1*varrho_2*(-m_alp1*inv_a2_2*b1*Delta_e1[i]+(1-m_alp1)*inv_a1_2*b2*Delta_e2[i])*_vec_normal[idim];
1090                                 _absAroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar+i]+=-dpi2*varrho_2*(-m_alp1*inv_a2_2*b1*Delta_e1[i]+(1-m_alp1)*inv_a1_2*b2*Delta_e2[i])*_vec_normal[idim];
1091                         }
1092                 }
1093                 // contribution of pressure gradient in momentum equation
1094                 for (int idim=0; idim<_Ndim; idim++){
1095                         _absAroe[ (1+idim)*_nVar] -= alpha*varrho_2*m_rho2*_vec_normal[idim];
1096                         _absAroe[ (1+idim)*_nVar+_Ndim+1] -=alpha* varrho_2*m_rho1*_vec_normal[idim];
1097                         _absAroe[(2+_Ndim+idim)*_nVar] -=  (1-alpha)*varrho_2*m_rho2*_vec_normal[idim];
1098                         _absAroe[(2+_Ndim+idim)*_nVar+_Ndim+1] -= (1-alpha)* varrho_2*m_rho1*_vec_normal[idim];
1099                         for  (int i=0; i<_nVar; i++){
1100                                 _absAroe[ (1+idim)*_nVar+i]-=alpha*varrho_2*(-m_alp1*m_rho2*b1*Delta_e1[i] -(1-m_alp1)*m_rho1*b2*Delta_e2[i])*_vec_normal[idim];
1101                                 _absAroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar+i]-=(1-alpha)*varrho_2*(-m_alp1*m_rho2*b1*Delta_e1[i] -(1-m_alp1)*m_rho1*b2*Delta_e2[i])*_vec_normal[idim];
1102                         }
1103                 }
1104                 // last row (total energy) To be changed soon
1105                 for (int i=0; i<_nVar; i++){
1106                         _Aroe[(2*_Ndim+2)*_nVar +i] = A5[i];
1107                 }
1108                 double signu1,signu2;
1109                 if(u1_n>0)
1110                         signu1=1;
1111                 else if (u1_n<0)
1112                         signu1=-1;
1113                 else
1114                         signu1=0;
1115                 if(u2_n>0)
1116                         signu2=1;
1117                 else if (u2_n<0)
1118                         signu2=-1;
1119                 else
1120                         signu2=0;
1121
1122                 for(int i=0; i<(1+_Ndim)*_nVar;i++){
1123                         _absAroe[i] *= signu1;
1124                         _absAroe[i+(1+_Ndim)*_nVar] *= signu2;
1125                 }
1126
1127         }
1128         if(_spaceScheme ==upwind || _spaceScheme==pressureCorrection || _spaceScheme ==lowMach)//calcul de la valeur absolue
1129         {               //on ordonne les deux premieres valeurs
1130                 if(valeurs_propres_dist[1].real()<valeurs_propres_dist[0].real())
1131                 {
1132                         tmp=valeurs_propres_dist[0];
1133                         valeurs_propres_dist[0]=valeurs_propres_dist[1];
1134                         valeurs_propres_dist[1]=tmp;
1135                 }
1136                 vector< complex< double > > y (taille_vp,0);
1137                 for( int i=0 ; i<taille_vp ; i++)
1138                         y[i] = Polynoms::abs_generalise(valeurs_propres_dist[i]);
1139
1140                 if(_entropicCorrection)
1141                 {
1142                         double entShift0 = 0;
1143                         double entShift1 = 0;
1144                         entropicShift(_vec_normal,entShift0,entShift1);
1145                         //cout<<"entShift0= "<<entShift0<<endl;
1146                         for( int i=0 ; i<taille_vp ; i++)
1147                         {
1148                                 //cout<<"y["<<i<<"]="<<y[i].real()<<endl;
1149                                 y[i] += max(entShift0,entShift1);
1150                         }
1151                 }
1152
1153                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1154                 {
1155                         cout<<"valeurs propres"<<endl;
1156                         for( int i=0 ; i<taille_vp ; i++)
1157                                 cout<<valeurs_propres_dist[i] <<", "<<endl;
1158                         cout<<"valeurs à atteindre"<<endl;
1159                         for( int i=0 ; i<taille_vp ; i++)
1160                                 cout<<y[i] <<", "<<endl;
1161                 }
1162                 Polynoms::abs_par_interp_directe(taille_vp,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
1163
1164                 if( _spaceScheme ==pressureCorrection){
1165                         for( int i=0 ; i<_Ndim ; i++)
1166                                 for( int j=0 ; j<_Ndim ; j++){
1167                                         _absAroe[(1+i)*_nVar+1+j]-=alpha*(valeurs_propres_dist[1].real()-valeurs_propres_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
1168                                         _absAroe[(2+_Ndim+i)*_nVar+2+_Ndim+j]-=(1-alpha)*(valeurs_propres_dist[1].real()-valeurs_propres_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
1169                                 }
1170                 }
1171                 else if( _spaceScheme ==lowMach){
1172                         double M=max(fabs(u1_n),fabs(u2_n))/maxvploc;
1173                         for( int i=0 ; i<_Ndim ; i++)
1174                                 for( int j=0 ; j<_Ndim ; j++){
1175                                         _absAroe[(1+i)*_nVar+1+j]-=(1-M)*alpha*(valeurs_propres_dist[1].real()-valeurs_propres_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
1176                                         _absAroe[(2+_Ndim+i)*_nVar+2+_Ndim+j]-=(1-M)*(1-alpha)*(valeurs_propres_dist[1].real()-valeurs_propres_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
1177                                 }
1178                 }
1179         }
1180
1181         //Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source
1182
1183         if(_entropicCorrection || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
1184                 InvMatriceRoe( valeurs_propres_dist);
1185                 Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
1186         }
1187         else if (_spaceScheme==upwind)//upwind sans entropic
1188                 SigneMatriceRoe( valeurs_propres_dist);
1189         else if (_spaceScheme==centered)//centre  sans entropic
1190         {
1191                 for(int i=0; i<_nVar*_nVar;i++)
1192                         _signAroe[i] = 0;
1193         }
1194         else  if(_spaceScheme ==staggered )
1195         {
1196                 double signu1,signu2;
1197                 if(u1_n>0)
1198                         signu1=1;
1199                 else if (u1_n<0)
1200                         signu1=-1;
1201                 else
1202                         signu1=0;
1203                 if(u2_n>0)
1204                         signu2=1;
1205                 else if (u2_n<0)
1206                         signu2=-1;
1207                 else
1208                         signu2=0;
1209                 for(int i=0; i<_nVar*_nVar;i++)
1210                         _signAroe[i] = 0;
1211                 _signAroe[0] = signu1;
1212                 _signAroe[(1+_Ndim)*_nVar +1+_Ndim] = signu2;
1213                 for(int i=2; i<_nVar-1;i++){
1214                         _signAroe[i*_nVar+i] = -signu1;
1215                         _signAroe[(i+1+_Ndim)*_nVar+i+1+_Ndim] = -signu2;
1216                 }
1217                 //_signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
1218         }
1219         else
1220                 throw CdmathException("FiveEqsTwoFluid::convectionMatrices: well balanced option not treated");
1221
1222         for(int i=0; i<_nVar*_nVar;i++)
1223         {
1224                 _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
1225                 _AroePlus[i]  = (_Aroe[i]+_absAroe[i])/2;
1226         }
1227         if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//Implicitation using primitive variables
1228                 for(int i=0; i<_nVar*_nVar;i++)
1229                         _AroeMinusImplicit[i] = (_AroeImplicit[i]-_absAroeImplicit[i])/2;
1230         else
1231                 for(int i=0; i<_nVar*_nVar;i++)
1232                         _AroeMinusImplicit[i] = _AroeMinus[i];
1233
1234         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1235         {
1236                 cout<<"Matrice de Roe"<<endl;
1237                 for(int i=0; i<_nVar;i++){
1238                         for(int j=0; j<_nVar;j++)
1239                                 cout<<_Aroe[i*_nVar+j]<<" , ";
1240                         cout<<endl;
1241                 }
1242                 cout<<"Valeur absolue matrice de Roe"<<endl;
1243                 for(int i=0; i<_nVar;i++){
1244                         for(int j=0; j<_nVar;j++)
1245                                 cout<<_absAroe[i*_nVar+j]<<" , ";
1246                         cout<<endl;
1247                 }
1248                 cout<<"Signe matrice de Roe"<<endl;
1249                 for(int i=0; i<_nVar;i++){
1250                         for(int j=0; j<_nVar;j++)
1251                                 cout<<_signAroe[i*_nVar+j]<<" , ";
1252                         cout<<endl;
1253                 }
1254                 cout<<endl<<"Matrice _AroeMinus"<<endl;
1255                 for(int i=0; i<_nVar;i++)
1256                 {
1257                         for(int j=0; j<_nVar;j++)
1258                                 cout << _AroeMinus[i*_nVar+j]<< " , ";
1259                         cout<<endl;
1260                 }
1261                 cout<<endl<<"Matrice _AroePlus"<<endl;
1262                 for(int i=0; i<_nVar;i++)
1263                 {
1264                         for(int j=0; j<_nVar;j++)
1265                                 cout << _AroePlus[i*_nVar+j]<< " , ";
1266                         cout<<endl;
1267                 }
1268         }
1269 }
1270
1271 void FiveEqsTwoFluid::jacobianDiff(const int &j, string nameOfGroup)
1272 {
1273         int k;
1274         for(k=0; k<_nVar*_nVar;k++)
1275                 _JcbDiff[k] = 0;
1276         if (_limitField[nameOfGroup].bcType==Wall){
1277                 _idm[0] = _nVar*j;
1278                 for(k=1; k<_nVar; k++)
1279                         _idm[k] = _idm[k-1] + 1;
1280                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1281                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1282
1283                 double pression=_Vj[1];//pressure inside
1284                 double T=_Vj[_nVar-1];//temperature outside
1285                 double rho_v=_fluides[0]->getDensity(pression,T);
1286                 double rho_l=_fluides[1]->getDensity(pression,T);
1287                 _JcbDiff[0] = 1;
1288                 _JcbDiff[(1+_Ndim)*_nVar +1+_Ndim] = 1;
1289                 _JcbDiff[_nVar]=_limitField[nameOfGroup].v_x[0];
1290                 _JcbDiff[(2+_Ndim)*_nVar +1+_Ndim] =_limitField[nameOfGroup].v_x[1];
1291                 double v2_v=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
1292                 double v2_l=_limitField[nameOfGroup].v_x[1]*_limitField[nameOfGroup].v_x[1];
1293                 if(_Ndim>1)
1294                 {
1295                         _JcbDiff[2*_nVar]=_limitField[nameOfGroup].v_y[0];
1296                         _JcbDiff[(3+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_y[1];
1297                         v2_v+=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
1298                         v2_l+=_limitField[nameOfGroup].v_y[1]*_limitField[nameOfGroup].v_y[1];
1299                         if(_Ndim==3)
1300                         {
1301                                 _JcbDiff[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1302                                 _JcbDiff[(4+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_z[1];
1303                                 v2_v+=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
1304                                 v2_l+=_limitField[nameOfGroup].v_z[1]*_limitField[nameOfGroup].v_z[1];
1305                         }
1306                 }
1307                 _JcbDiff[(_nVar-1)*_nVar]=         _fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho_v)+0.5*v2_v;
1308                 _JcbDiff[(_nVar-1)*_nVar +1+_Ndim]=_fluides[1]->getInternalEnergy(_limitField[nameOfGroup].T,rho_l)+0.5*v2_l;
1309         }
1310         else if (_limitField[nameOfGroup].bcType==Inlet){
1311                 /*
1312                 _JcbDiff[0] = 1;
1313                 _JcbDiff[(1+_Ndim)*_nVar +1+_Ndim] = 1;
1314                 _JcbDiff[_nVar]=_limitField[nameOfGroup].v_x[0];
1315                 _JcbDiff[(2+_Ndim)*_nVar +1+_Ndim] =_limitField[nameOfGroup].v_x[1];
1316                 double v2_v=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
1317                 double v2_l=_limitField[nameOfGroup].v_x[1]*_limitField[nameOfGroup].v_x[1];
1318                 if(_Ndim>1)
1319                 {
1320                         _JcbDiff[2*_nVar]=_limitField[nameOfGroup].v_y[0];
1321                         _JcbDiff[(3+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_y[1];
1322                         v2_v+=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
1323                         v2_l+=_limitField[nameOfGroup].v_y[1]*_limitField[nameOfGroup].v_y[1];
1324                         if(_Ndim==3)
1325                         {
1326                                 _JcbDiff[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1327                                 _JcbDiff[(4+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_z[1];
1328                                 v2_v+=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
1329                                 v2_l+=_limitField[nameOfGroup].v_z[1]*_limitField[nameOfGroup].v_z[1];
1330                         }
1331                 }
1332                 _JcbDiff[(_nVar-1)*_nVar]=         _fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho_v)+0.5*v2_v;
1333                 _JcbDiff[(_nVar-1)*_nVar +1+_Ndim]=_fluides[1]->getInternalEnergy(_limitField[nameOfGroup].T,rho_l)+0.5*v2_l;
1334                  */
1335         } else if (_limitField[nameOfGroup].bcType==Outlet){
1336                 //extraction de l etat courant et primitives
1337                 /*
1338                 _idm[0] = j*_nVar;
1339                 for(k=1; k<_nVar;k++)
1340                 {_idm[k] = _idm[k-1] + 1;}
1341                 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1342                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1343                  */
1344         }
1345         else if (_limitField[nameOfGroup].bcType!=Neumann && _limitField[nameOfGroup].bcType!=InletPressure){
1346                 cout<<"Condition  limite non traitee pour le bord "<<nameOfGroup<< endl;
1347                 throw CdmathException("FiveEqsTwoFluid::jacobianDiff: Condition  limite non traitee");
1348         }
1349 }
1350
1351
1352 void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *normale){
1353         //To do controle signe des vitesses pour CL entree/sortie
1354         int k;
1355         double v1_2=0, v2_2=0, q1_n=0, q2_n=0, u1_n=0, u2_n=0;//quantités de mouvement et vitesses normales à la face limite;
1356         double v1[_Ndim], v2[_Ndim];
1357
1358         _idm[0] = _nVar*j;
1359         for(k=1; k<_nVar; k++)
1360                 _idm[k] = _idm[k-1] + 1;
1361
1362         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
1363         VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1364
1365         for(k=0; k<_Ndim; k++){
1366                 q1_n+=_externalStates[(k+1)]*normale[k];
1367                 q2_n+=_externalStates[(k+1+1+_Ndim)]*normale[k];
1368                 u1_n+=_Vj[(k+2)]*normale[k];
1369                 u2_n+=_Vj[(k+2+_Ndim)]*normale[k];
1370         }
1371
1372         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1373         {
1374                 cout << "Boundary conditions for group "<< nameOfGroup<< ", inner cell j= "<<j << " face unit normal vector "<<endl;
1375                 for(k=0; k<_Ndim; k++){
1376                         cout<<normale[k]<<", ";
1377                 }
1378                 cout<<endl;
1379         }
1380
1381         if (_limitField[nameOfGroup].bcType==Wall){
1382
1383                 //Pour la convection, inversion du sens des vitesses
1384                 for(k=0; k<_Ndim; k++){
1385                         _externalStates[(k+1)]-= 2*q1_n*normale[k];
1386                         _externalStates[(k+1+1+_Ndim)]-= 2*q2_n*normale[k];
1387                         _Vj[(k+2)]-= 2*u1_n*normale[k];
1388                         _Vj[(k+2+_Ndim)]-= 2*u2_n*normale[k];
1389                 }
1390                 _idm[0] = 0;
1391                 for(k=1; k<_nVar; k++)
1392                         _idm[k] = _idm[k-1] + 1;
1393
1394                 VecAssemblyBegin(_Uext);
1395                 VecAssemblyBegin(_Vext);
1396                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
1397                 VecSetValues(_Vext, _nVar, _idm, _Vj, INSERT_VALUES);
1398                 VecAssemblyEnd(_Uext);
1399                 VecAssemblyEnd(_Vext);
1400
1401                 //Pour la diffusion, paroi à vitesses et temperature imposees
1402                 double pression=_Vj[1];//pressure inside
1403                 double T=_Vj[_nVar-1];//temperature outside
1404                 double rho_v=_fluides[0]->getDensity(pression,T);
1405                 double rho_l=_fluides[1]->getDensity(pression,T);
1406                 _externalStates[1]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
1407                 _externalStates[2+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_x[1];
1408                 if(_Ndim>1)
1409                 {
1410                         _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
1411                         _externalStates[3+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_y[1];
1412                         if(_Ndim==3)
1413                         {
1414                                 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
1415                                 _externalStates[4+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_z[1];
1416                         }
1417                 }
1418                 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho_v) + v1_2/2)
1419                                                                                                                                                                                                                                                                                                                                                                         +_externalStates[1+_Ndim]*(_fluides[1]->getInternalEnergy(_limitField[nameOfGroup].T,rho_l) + v2_2/2);
1420                 VecAssemblyBegin(_Uextdiff);
1421                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
1422                 VecAssemblyEnd(_Uextdiff);
1423         }
1424         else if (_limitField[nameOfGroup].bcType==Neumann){
1425                 _idm[0] = _nVar*j;
1426                 for(k=1; k<_nVar; k++)
1427                         _idm[k] = _idm[k-1] + 1;
1428
1429                 _idm[0] = 0;
1430                 for(k=1; k<_nVar; k++)
1431                         _idm[k] = _idm[k-1] + 1;
1432
1433                 VecAssemblyBegin(_Uext);
1434                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
1435                 VecAssemblyEnd(_Uext);
1436
1437                 VecAssemblyBegin(_Vext);
1438                 VecSetValues(_Vext, _nVar, _idm, _Vj, INSERT_VALUES);
1439                 VecAssemblyEnd(_Vext);
1440
1441                 VecAssemblyBegin(_Uextdiff);
1442                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
1443                 VecAssemblyEnd(_Uextdiff);
1444         }
1445         else if (_limitField[nameOfGroup].bcType==Inlet){
1446                 _idm[0] = _nVar*j;
1447                 for(k=1; k<_nVar; k++)
1448                         _idm[k] = _idm[k-1] + 1;
1449
1450                 double alpha=_limitField[nameOfGroup].alpha;//void fraction outside
1451                 double pression=_Vj[1];//pressure inside
1452                 double T=_limitField[nameOfGroup].T;//temperature outside
1453                 double rho_v=_fluides[0]->getDensity(pression,T);
1454                 double rho_l=_fluides[1]->getDensity(pression,T);
1455                 //cout<<"Inlet alpha= "<<alpha<<" pression= "<<pression<<" temperature= "<<T<<" velocity gas "<<_limitField[nameOfGroup].v_x[0]<<" velocity liq "<<_limitField[nameOfGroup].v_x[1]<<endl;
1456
1457                 _externalStates[0]=alpha*rho_v;
1458                 _externalStates[1 + _Ndim] = (1-alpha)*rho_l;
1459                 v1[0] = _limitField[nameOfGroup].v_x[0];
1460                 v2[0] = _limitField[nameOfGroup].v_x[1];
1461                 if (_Ndim >1){
1462                         v1[1] = _limitField[nameOfGroup].v_y[0];
1463                         v2[1] = _limitField[nameOfGroup].v_y[1];
1464                 }
1465                 if (_Ndim >2){
1466                         v1[2] = _limitField[nameOfGroup].v_z[0];
1467                         v2[2] = _limitField[nameOfGroup].v_z[1];
1468                 }
1469                 for (int idim=0;idim<_Ndim;idim++){
1470                         _externalStates[1 + idim] = v1[idim]* _externalStates[0]; // phase 1
1471                         _externalStates[2 + _Ndim + idim] = v2[idim]* _externalStates[1+_Ndim]; // phase 2
1472                         v1_2 += v1[idim]*v1[idim];
1473                         v2_2 += v2[idim]*v2[idim];
1474                         _Vj[2+idim] = v1[idim];
1475                         _Vj[2+_Ndim+idim] = v2[idim];
1476                 }
1477                 _externalStates[_nVar-1] = _externalStates[0]      *(_fluides[0]->getInternalEnergy(T,rho_v) + v1_2/2)
1478                                                                                                                                                                                                                                                                                                                                         +_externalStates[1+_Ndim]*(_fluides[1]->getInternalEnergy(T,rho_l) + v2_2/2);
1479                 // _Vj external primitives
1480                 _Vj[0] = alpha;
1481                 _Vj[_nVar-1] = T;
1482
1483                 _idm[0] = 0;
1484                 for(k=1; k<_nVar; k++)
1485                         _idm[k] = _idm[k-1] + 1;
1486
1487                 VecAssemblyBegin(_Uext);
1488                 VecAssemblyBegin(_Vext);
1489                 VecAssemblyBegin(_Uextdiff);
1490                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
1491                 VecSetValues(_Vext, _nVar, _idm, _Vj, INSERT_VALUES);
1492                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
1493                 VecAssemblyEnd(_Uext);
1494                 VecAssemblyEnd(_Vext);
1495                 VecAssemblyEnd(_Uextdiff);
1496         }
1497         else if (_limitField[nameOfGroup].bcType==InletPressure){
1498                 _idm[0] = _nVar*j;
1499                 for(k=1; k<_nVar; k++)
1500                         _idm[k] = _idm[k-1] + 1;
1501
1502                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
1503                 Cell Cj=_mesh.getCell(j);
1504                 double hydroPress=Cj.x()*_GravityField3d[0];
1505                 if(_Ndim>1){
1506                         hydroPress+=Cj.y()*_GravityField3d[1];
1507                         if(_Ndim>2)
1508                                 hydroPress+=Cj.z()*_GravityField3d[2];
1509                 }
1510                 hydroPress*=_externalStates[0]+_externalStates[_Ndim];//multiplication by rho the total density
1511
1512                 //Building the external state
1513                 double alpha=_limitField[nameOfGroup].alpha;
1514                 double pression=_limitField[nameOfGroup].p+hydroPress;
1515                 double T=_limitField[nameOfGroup].T;
1516                 double rho_v=_fluides[0]->getDensity(pression,T);
1517                 double rho_l=_fluides[1]->getDensity(pression,T);
1518                 _externalStates[0]=alpha*rho_v;
1519                 _externalStates[1+_Ndim]=(1-alpha)*rho_l;
1520
1521                 for(int idim=0;idim<_Ndim;idim++){
1522                         _externalStates[idim+1]=_externalStates[0]*_Vj[idim+2];
1523                         _externalStates[idim+2+_Ndim]=_externalStates[1+_Ndim]*_Vj[idim+2+_Ndim];
1524                         v1_2+=_Vj[2+idim]*_Vj[2+idim];
1525                         v2_2+=_Vj[2+_Ndim+idim]*_Vj[2+_Ndim+idim];
1526                 }
1527                 _externalStates[_nVar-1]=    alpha *rho_v*(_fluides[0]->getInternalEnergy(T,rho_v)+v1_2/2)
1528                                                                                                                                                                                                                                                                                                                                                 +(1-alpha)*rho_l*(_fluides[1]->getInternalEnergy(T,rho_l)+v2_2/2);
1529                 // _Vj external primitives
1530                 _Vj[0] = alpha;
1531                 _Vj[1] = pression;
1532                 _Vj[_nVar-1] = T;
1533
1534                 _idm[0] = 0;
1535                 for(k=1; k<_nVar; k++)
1536                         _idm[k] = _idm[k-1] + 1;
1537                 VecAssemblyBegin(_Uext);
1538                 VecAssemblyBegin(_Vext);
1539                 VecAssemblyBegin(_Uextdiff);
1540                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
1541                 VecSetValues(_Vext, _nVar, _idm, _Vj, INSERT_VALUES);
1542                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
1543                 VecAssemblyEnd(_Uext);
1544                 VecAssemblyEnd(_Vext);
1545                 VecAssemblyEnd(_Uextdiff);
1546         }
1547         else if (_limitField[nameOfGroup].bcType==Outlet){
1548                 _idm[0] = _nVar*j;
1549                 for(k=1; k<_nVar; k++)
1550                         _idm[k] = _idm[k-1] + 1;
1551
1552                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
1553                 Cell Cj=_mesh.getCell(j);
1554                 double hydroPress=Cj.x()*_GravityField3d[0];
1555                 if(_Ndim>1){
1556                         hydroPress+=Cj.y()*_GravityField3d[1];
1557                         if(_Ndim>2)
1558                                 hydroPress+=Cj.z()*_GravityField3d[2];
1559                 }
1560                 hydroPress*=_externalStates[0]+_externalStates[_Ndim];//multiplication by rho the total density
1561
1562                 //Building the external state
1563                 double pression_int=_Vj[1];
1564                 double pression_ext=_limitField[nameOfGroup].p+hydroPress;
1565                 double T=_Vj[_nVar-1];
1566                 double rho_v_int=_fluides[0]->getDensity(pression_int,T);
1567                 double rho_l_int=_fluides[1]->getDensity(pression_int,T);
1568                 double rho_v_ext=_fluides[0]->getDensity(pression_ext,T);
1569                 double rho_l_ext=_fluides[1]->getDensity(pression_ext,T);
1570
1571                 for(k=0;k<1+_Ndim;k++){
1572                         _externalStates[k]*=rho_v_ext/rho_v_int;
1573                         _externalStates[k+1+_Ndim]*=rho_l_ext/rho_l_int;
1574                 }
1575                 double alpha=_Vj[0];
1576                 //cout<<"Outlet alpha= "<<alpha<<" pression int= "<<pression_int<<" pression ext= "<<pression_ext<<" temperature= "<<T<<" velocity gas "<<_Uj[2]<<" velocity liq "<<_Uj[2+_Ndim]<<endl;
1577                 for(int idim=0;idim<_Ndim;idim++){
1578                         v1_2+=_Vj[2+idim]*_Vj[2+idim];
1579                         v2_2+=_Vj[2+_Ndim+idim]*_Vj[2+_Ndim+idim];
1580                 }
1581                 _externalStates[_nVar-1]=alpha*rho_v_ext*(_fluides[0]->getInternalEnergy(T,rho_v_int)+v1_2/2)+(1-alpha)*rho_l_ext*(_fluides[1]->getInternalEnergy(T,rho_l_int)+v2_2/2);
1582
1583                 // _Vj external primitives
1584                 _Vj[1] = pression_ext;
1585
1586                 _idm[0] = 0;
1587                 for(k=1; k<_nVar; k++)
1588                         _idm[k] = _idm[k-1] + 1;
1589                 VecAssemblyBegin(_Uext);
1590                 VecAssemblyBegin(_Vext);
1591                 VecAssemblyBegin(_Uextdiff);
1592                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
1593                 VecSetValues(_Vext, _nVar, _idm, _Vj, INSERT_VALUES);
1594                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
1595                 VecAssemblyEnd(_Uext);
1596                 VecAssemblyEnd(_Vext);
1597                 VecAssemblyEnd(_Uextdiff);
1598         }
1599         else{
1600                 cout<<"Boundary condition not set for boundary named "<<nameOfGroup<<endl;
1601                 cout<<"Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
1602                 throw CdmathException("Unknown boundary condition");
1603         }
1604 }
1605
1606 void FiveEqsTwoFluid::addDiffusionToSecondMember
1607 (               const int &i,
1608                 const int &j,
1609                 bool isBord)
1610 {
1611         double Tm=_Udiff[_nVar-1];
1612         double lambdal=_fluides[1]->getConductivity(Tm);
1613         double lambdav=_fluides[0]->getConductivity(Tm);
1614         double mu1 =_fluides[0]->getViscosity(Tm);
1615         double mu2 = _fluides[1]->getViscosity(Tm);
1616
1617         if(mu1==0 && mu2 ==0 && lambdav==0 && lambdal==0 && _heatTransfertCoeff==0)
1618                 return;
1619
1620         //extraction des valeurs
1621         _idm[0] = _nVar*i;
1622         for(int k=1; k<_nVar; k++)
1623                 _idm[k] = _idm[k-1] + 1;
1624
1625         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1626         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1627         {
1628                 cout << "Contribution diffusion: variables primitives maille " << i<<endl;
1629                 for(int q=0; q<_nVar; q++)
1630                 {
1631                         cout << _Vi[q] << endl;
1632                 }
1633                 cout << endl;
1634         }
1635
1636         if(!isBord ){
1637                 for(int k=0; k<_nVar; k++)
1638                         _idm[k] = _nVar*j + k;
1639
1640                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1641         }
1642         else
1643         {
1644                 lambdal=max(lambdal,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
1645
1646                 for(int k=0; k<_nVar; k++)
1647                         _idm[k] = k;
1648                 VecGetValues(_Uextdiff, _nVar, _idm, _phi);
1649                 consToPrim(_phi,_Vj);
1650         }
1651
1652         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1653         {
1654                 cout << "Contribution diffusion: variables primitives maille " <<j <<endl;
1655                 for(int q=0; q<_nVar; q++)
1656                 {
1657                         cout << _Vj[q] << endl;
1658                 }
1659                 cout << endl;
1660         }
1661         double alpha=(_Vj[0]+_Vi[0])/2;
1662         double lambda = alpha*lambdav+(1-alpha)*lambdal;
1663         //on n'a pas de contribution sur la masse
1664         _phi[0]=0;
1665         _phi[_Ndim+1]=0;
1666         //contribution visqueuse sur la quantite de mouvement
1667         for(int k=1; k<_Ndim+1; k++)
1668         {
1669                 _phi[k]         = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu1*alpha    *(_Vj[2+k] - _Vi[2+k]);//attention car primitif=alpha p u1 u2 T
1670                 _phi[k+_Ndim+1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu2*(1-alpha)*(_Vj[2+k+_Ndim] - _Vi[2+k+_Ndim]);
1671         }
1672         _phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda  *(_Vj[_nVar-1] - _Vi[_nVar-1]);
1673         _idm[0] = i;
1674         VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1675
1676         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1677         {
1678                 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
1679                 for(int q=0; q<_nVar; q++)
1680                 {
1681                         cout << _phi[q] << endl;
1682                 }
1683                 cout << endl;
1684         }
1685
1686         if(!isBord)
1687         {
1688                 //On change de signe pour l'autre contribution
1689                 for(int k=0; k<_nVar; k++)
1690                         _phi[k] *= -_inv_dxj/_inv_dxi;
1691                 _idm[0] = j;
1692
1693                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1694         }
1695
1696         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1697         {
1698                 cout << "Contribution diffusion au 2nd membre pour la maille  " << j << ": "<<endl;
1699                 for(int q=0; q<_nVar; q++)
1700                 {
1701                         cout << _phi[q] << endl;
1702                 }
1703                 cout << endl;
1704
1705                 if(_timeScheme==Implicit)
1706                 {
1707                         cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
1708                         for(int i=0; i<_nVar; i++)
1709                         {
1710                                 for(int j=0; j<_nVar; j++)
1711                                         cout << _Diffusion[i*_nVar+j]<<", ";
1712                                 cout << endl;
1713                         }
1714                         cout << endl;
1715                 }
1716         }
1717 }
1718
1719 void FiveEqsTwoFluid::jacobian(const int &j, string nameOfGroup, double * normale)
1720 {
1721         int k;
1722         for(k=0; k<_nVar*_nVar;k++)
1723                 _Jcb[k] = 0;//No implicitation at this stage
1724
1725         // loop of boundary types
1726         if (_limitField[nameOfGroup].bcType==Wall)
1727         {
1728                 for(k=0; k<_nVar;k++)
1729                         _Jcb[k*_nVar + k] = 1;
1730                 for(k=1; k<1+_Ndim;k++)
1731                         for(int l=1; l<1+_Ndim;l++){
1732                                 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
1733                                 _Jcb[(k+1+_Ndim)*_nVar + l+1+_Ndim] -= 2*normale[k-1]*normale[l-1];
1734                         }
1735         }
1736         //not wall
1737         else if (_limitField[nameOfGroup].bcType==Inlet)
1738         {
1739                 _Jcb[0] = 1;
1740                 _Jcb[(1+_Ndim)*_nVar +1+_Ndim] = 1;
1741                 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];
1742                 _Jcb[(2+_Ndim)*_nVar +1+_Ndim] =_limitField[nameOfGroup].v_x[1];
1743                 if(_Ndim>1)
1744                 {
1745                         _Jcb[2*_nVar]=_limitField[nameOfGroup].v_y[0];
1746                         _Jcb[(3+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_y[1];
1747                         if(_Ndim==3)
1748                         {
1749                                 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1750                                 _Jcb[(4+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_z[1];
1751
1752                         }
1753                 }
1754         }
1755         // not wall, not inlet
1756         else if (_limitField[nameOfGroup].bcType==Outlet){
1757                 _idm[0] = j*_nVar;
1758                 for(k=1; k<_nVar;k++)
1759                 {_idm[k] = _idm[k-1] + 1;}
1760                 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1761                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1762         }
1763         else  if (_limitField[nameOfGroup].bcType!=Neumann && _limitField[nameOfGroup].bcType!=InletPressure)// not wall, not inlet, not outlet
1764         {
1765                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1766                 throw CdmathException("FiveEqs::jacobianDiff: This boundary condition is not treated");
1767         }
1768 }
1769
1770
1771 void FiveEqsTwoFluid::primToCons(const double *P, const int &i, double *W, const int &j){
1772         //P= alpha, p, u1, u2, T
1773         //W=m1,q1,m2,q2,rhoE =alpha1*rho1*(e1+u1^2/2)+alpha2*rho2*(e2+u2^2/2)
1774         double alpha=P[i*_nVar];
1775         double pression=P[i*_nVar+1];
1776         double temperature=P[i*_nVar+_nVar-1];
1777         double rho_v=_fluides[0]->getDensity(pression,temperature);
1778         double rho_l=_fluides[1]->getDensity(pression,temperature);
1779         double u1_sq=0, u2_sq=0;
1780
1781         W[j*_nVar] = alpha*rho_v;
1782         W[j*_nVar+1+_Ndim] = (1-alpha)*rho_l;
1783         // alpha*rho*u
1784         for(int k=0; k<_Ndim; k++)
1785         {
1786                 W[j*_nVar+(k+1)] = W[j*_nVar]*P[i*_nVar+(k+2)];//alpha1*rho1*u1
1787                 W[j*_nVar+(k+1)+1+_Ndim] = W[j*_nVar+1+_Ndim]*P[i*_nVar+(k+2)+_Ndim];//alpha2*rho2*u2
1788         }
1789         // total energy
1790         W[j*_nVar+_nVar-1] = W[j*(_nVar)]* _fluides[0]->getInternalEnergy(temperature,rho_v)+W[j*(_nVar)+1+_Ndim]* _fluides[1]->getInternalEnergy(temperature,rho_l);
1791         for(int k=0; k<_Ndim; k++){
1792                 u1_sq+=P[i*_nVar+(k+2)]*P[i*_nVar+(k+2)];
1793                 u2_sq+=P[i*_nVar+(k+2)+_Ndim]*P[i*_nVar+(k+2)+_Ndim];
1794         }
1795         W[j*_nVar+_nVar-1] += (W[j*_nVar]*u1_sq+W[j*_nVar+1+_Ndim]*u2_sq)*0.5;
1796 }
1797
1798 void FiveEqsTwoFluid::consToPrim(const double *Wcons, double* Wprim,double porosity)//To do: treat porosity
1799 {
1800         //Wprim= alpha, p, u1, u2, T
1801         //Wcons=m1,q1,m2,q2,rhoE
1802         double m_v=Wcons[0];
1803         double m_l=Wcons[1+_Ndim];
1804         double q1_sq = 0,q2_sq = 0;
1805         _minm1=min(m_v,_minm1);
1806         _minm2=min(m_l,_minm2);
1807
1808         if(m_v<-_precision || m_l<-_precision){
1809                 _nbMaillesNeg+=1;
1810         }
1811         for(int k=0;k<_Ndim;k++){
1812                 q1_sq += Wcons[k+1]*Wcons[k+1];
1813                 q2_sq += Wcons[k+1+1+_Ndim]*Wcons[k+1+1+_Ndim];
1814         }
1815         if(Wcons[0]>0)//_precision*_precision*_precision)
1816                 q1_sq /= Wcons[0];      //alpha1 rho1 u1²
1817         else
1818                 q1_sq = 0;
1819         if(Wcons[1+_Ndim]>0)//_precision*_precision*_precision)
1820                 q2_sq /= Wcons[1+_Ndim];        //alpha2 rho2 u1²
1821         else
1822                 q2_sq = 0;
1823         double rho_m_e_m=Wcons[_nVar-1] -0.5*(q1_sq+q2_sq);
1824         //calcul de la temperature et de la pression pour une loi stiffened gas
1825         double temperature= (rho_m_e_m-m_v*static_cast<StiffenedGas*>(_fluides[0])->getInternalEnergy(0)-m_l*static_cast<StiffenedGas*>(_fluides[1])->getInternalEnergy(0))/(m_v*_fluides[0]->constante("cv")+m_l*_fluides[1]->constante("cv"));
1826         double e_v=static_cast<StiffenedGas*>(_fluides[0])->getInternalEnergy(temperature);
1827         double e_l=static_cast<StiffenedGas*>(_fluides[1])->getInternalEnergy(temperature);
1828         double gamma_v=_fluides[0]->constante("gamma");
1829         double gamma_l=_fluides[1]->constante("gamma");
1830         double Pinf_v=- gamma_v*_fluides[0]->constante("p0");
1831         double Pinf_l=- gamma_l*_fluides[1]->constante("p0");
1832         double a=1;
1833         double b=-(Pinf_v+m_v*(gamma_v-1)*e_v+Pinf_l+m_l*(gamma_l-1)*e_l);
1834         double c=Pinf_v*Pinf_l+Pinf_v*m_l*(gamma_l-1)*e_l+ Pinf_l*m_v*(gamma_v-1)*e_v;
1835         double delta= b*b-4*a*c;
1836         if(delta<0){
1837                 cout<<"delta= "<<delta<<" <0"<<endl;
1838                 *_runLogFile<<"delta= "<<delta<<" <0"<<endl;
1839                 throw CdmathException("FiveEqsTwoFluid::consToPrim: Failed to compute pressure");
1840         }
1841         double pression=(-b+sqrt(delta))/(2*a);
1842         if (pression < 1){
1843                 cout << "pressure = "<< pression << " < 1 Pa " << endl;
1844                 cout << "Conservative state = ";
1845                 for(int k=0; k<_nVar; k++){
1846                         cout<<Wcons[k]<<", ";
1847                 }
1848                 cout<<endl;
1849                 *_runLogFile << "FiveEqsTwoFluid::consToPrim: Failed to compute pressure = "<< pression << " < 1 Pa " << endl;
1850                 throw CdmathException("FiveEqsTwoFluid::consToPrim: Failed to compute pressure");
1851         }
1852
1853         double rho_v=_fluides[0]->getDensity(pression,temperature);
1854         double alpha=m_v/rho_v;
1855         Wprim[0]= alpha;
1856         Wprim[1] =  pression;
1857         for(int k=0;k<_Ndim;k++){//vitesses
1858                 if(Wcons[0]>0)//_precision*_precision*_precision)
1859                         Wprim[k+2] = Wcons[k+1]/Wcons[0];
1860                 else
1861                         Wprim[k+2] = Wcons[k+2+_Ndim]/Wcons[1+_Ndim];
1862                 if(Wcons[1+_Ndim]>0)//_precision*_precision*_precision)
1863                         Wprim[k+2+_Ndim] = Wcons[k+2+_Ndim]/Wcons[1+_Ndim];
1864                 else
1865                         Wprim[k+2+_Ndim] = Wcons[k+1]/Wcons[0];
1866         }
1867         Wprim[_nVar-1] = temperature;
1868 }
1869
1870 void FiveEqsTwoFluid::entropicShift(double* n, double& vpcorr0, double& vpcorr1)
1871 {
1872
1873         // parameters of function dgeev_ (compute the eigenvalues)
1874         int LDA, LDVL,LWORK, SDIM,LDVR;
1875         LDA = _nVar;
1876         LDVL = _nVar;
1877         LDVR = _nVar;
1878         LWORK = 20*_nVar;
1879         char jobvl[]="N", jobvr[]="N";
1880         double WORK[LWORK], JacoMat[_nVar*_nVar],egvaReal[_nVar],egvaImag[_nVar],egVectorL[_nVar*_nVar],egVectorR[_nVar*_nVar];
1881         int info_l = 0, info_r = 0;
1882
1883         /******** Left: Compute the eigenvalues and eigenvectors of Jacobian Matrix (using lapack)********/
1884         convectionJacobianMatrix(_l, n);
1885         for (int i=0; i<_nVar*_nVar; i++){
1886                 JacoMat[i] = _JacoMat[i];
1887         }
1888         dgeev_(jobvl, jobvl,  &_nVar,
1889                         JacoMat,&LDA,egvaReal,egvaImag, egVectorL,
1890                         &LDVL,egVectorR,
1891                         &LDVR, WORK,&LWORK,
1892                         &info_l);
1893
1894         //      /******** Right: Compute the eigenvalues and eigenvectors of Jacobian Matrix (using lapack)********/
1895         convectionJacobianMatrix(_r, n);
1896         for (int i=0; i<_nVar*_nVar; i++){
1897                 JacoMat[i] = _JacoMat[i];
1898         }
1899         dgeev_(jobvl, jobvl,  &_nVar,
1900                         JacoMat,&LDA,egvaReal,egvaImag, egVectorL,
1901                         &LDVL,egVectorR,
1902                         &LDVR, WORK,&LWORK,
1903                         &info_r);
1904
1905         if (info_l < 0 || info_r < 0)
1906         {
1907                 cout<<"Warning FiveEqsTwoFluid::entropicShift: dgeev_ did not compute all the eigenvalues, trying heuristic entropy correction "<<endl;
1908                 double u1l_n=0, u2l_n=0, u1r_n=0, u2r_n=0;
1909                 for (int idim=0; idim<_Ndim; idim++){
1910                         u1l_n= _l[2+idim]      *n[idim];
1911                         u2l_n= _l[2+idim+_Ndim]*n[idim];
1912                         u1r_n= _r[2+idim]      *n[idim];
1913                         u2r_n= _r[2+idim+_Ndim]*n[idim];
1914                 }
1915
1916                 vpcorr0 =max(fabs(u1l_n-u1r_n),fabs(u2l_n-u2r_n));
1917                 vpcorr1 = vpcorr0;
1918         }
1919         else
1920         {
1921                 std::vector< std::complex<double> > eigValuesLeft(_nVar);
1922                 std::vector< std::complex<double> > eigValuesRight(_nVar);
1923                 for(int j=0; j<_nVar; j++){
1924                         eigValuesLeft[j] = complex<double>(egvaReal[j],egvaImag[j]);
1925                         eigValuesRight[j] = complex<double>(egvaReal[j],egvaImag[j]);
1926                 }
1927                 int sizeLeft =  Polynoms::new_tri_selectif(eigValuesLeft, eigValuesLeft.size(), _precision);
1928                 int sizeRight =  Polynoms::new_tri_selectif(eigValuesRight, eigValuesRight.size(), _precision);
1929                 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1930                 {
1931                         cout<<" Eigenvalue of JacoMat Left: " << endl;
1932                         for(int i=0; i<sizeLeft; i++)
1933                                 cout<<eigValuesLeft[i] << ", "<<endl;
1934                 }
1935                 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1936                 {
1937                         cout<<" Eigenvalue of JacoMat Right: " << endl;
1938                         for(int i=0; i<sizeRight; i++)
1939                                 cout<<eigValuesRight[i] << ", "<<endl;
1940                 }
1941                 vpcorr0 = 0;
1942                 for (int i=1; i<min(sizeLeft,sizeRight)-1; i++)
1943                         vpcorr0 = max(vpcorr0, abs(eigValuesRight[i]-eigValuesLeft[i]));// Kieu
1944                 vpcorr1 = vpcorr0;
1945         }
1946 }
1947
1948 void FiveEqsTwoFluid::entropicShift(double* n)
1949 {
1950
1951         // parameters of function dgeev_ (compute the eigenvalues)
1952         int LDA, LDVL,LWORK, SDIM,LDVR;
1953         LDA = _nVar;
1954         LDVL = _nVar;
1955         LDVR = _nVar;
1956         LWORK = 20*_nVar;
1957         char jobvl[]="N", jobvr[]="N";
1958         double WORK[LWORK], JacoMat[_nVar*_nVar],egvaReal[_nVar],egvaImag[_nVar],egVectorL[_nVar*_nVar],egVectorR[_nVar*_nVar];
1959         int info = 0;
1960         /******** Left: Compute the eigenvalues and eigenvectors of Jacobian Matrix (using lapack)********/
1961         convectionJacobianMatrix(_l, n);
1962         for (int i=0; i<_nVar*_nVar; i++){
1963                 JacoMat[i] = _JacoMat[i];
1964         }
1965         dgeev_(jobvl, jobvl,  &_nVar,
1966                         JacoMat,&LDA,egvaReal,egvaImag, egVectorL,
1967                         &LDVL,egVectorR,
1968                         &LDVR, WORK,&LWORK,
1969                         &info);
1970         if (info != 0){
1971                 *_runLogFile<<"FiveEqsTwoFluid::JacoMat: dgeev_ unsuccessful computation of the eigenvalues of Jacobian Matrix (left)"<<endl;
1972                 throw CdmathException(
1973                                 "FiveEqsTwoFluid::JacoMat: dgeev_ unsuccessful computation of the eigenvalues of Jacobian Matrix (left)");
1974         }
1975
1976         std::vector< std::complex<double> > eigValuesLeft(_nVar,0.);
1977         for(int j=0; j<_nVar; j++){
1978                 eigValuesLeft[j] = complex<double>(egvaReal[j],egvaImag[j]);
1979         }
1980         //      /******** Right: Compute the eigenvalues and eigenvectors of Jacobian Matrix (using lapack)********/
1981         convectionJacobianMatrix(_r, n);
1982         for (int i=0; i<_nVar*_nVar; i++){
1983                 JacoMat[i] = _JacoMat[i];
1984         }
1985         dgeev_(jobvl, jobvl,  &_nVar,
1986                         JacoMat,&LDA,egvaReal,egvaImag, egVectorL,
1987                         &LDVL,egVectorR,
1988                         &LDVR, WORK,&LWORK,
1989                         &info);
1990         if (info != 0)
1991                 {
1992                 *_runLogFile<<"FiveEqsTwoFluid::entropicShift: dgeev_ unsuccessful computation of the eigenvalues of Jacobian Matrix (right)"<<endl;
1993                 throw CdmathException(
1994                                 "FiveEqsTwoFluid::entropicShift: dgeev_ unsuccessful computation of the eigenvalues of Jacobian Matrix (right)");
1995                 }
1996
1997         std::vector< std::complex<double> > eigValuesRight(_nVar,0.);
1998         for(int j=0; j<_nVar; j++){
1999                 eigValuesRight[j] = complex<double>(egvaReal[j],egvaImag[j]);
2000         }
2001         int sizeLeft =  Polynoms::new_tri_selectif(eigValuesLeft, eigValuesLeft.size(), _precision);
2002         int sizeRight =  Polynoms::new_tri_selectif(eigValuesRight, eigValuesRight.size(), _precision);
2003         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
2004         {
2005                 cout<<" Eigenvalue of JacoMat Left: " << endl;
2006                 for(int i=0; i<sizeLeft; i++)
2007                         cout<<eigValuesLeft[i] << ", "<<endl;
2008                 cout<<" Eigenvalue of JacoMat Right: " << endl;
2009                 for(int i=0; i<sizeRight; i++)
2010                         cout<<eigValuesRight[i] << ", "<<endl;
2011         }
2012         for (int i=0; i<min(sizeLeft,sizeRight)-1; i++)
2013                 _entropicShift[i]= abs(eigValuesRight[i]-eigValuesLeft[i]);
2014 }
2015
2016 Vector FiveEqsTwoFluid::staggeredVFFCFlux()
2017 {
2018         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2019                 throw CdmathException("IsothermalTwoFluid::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
2020         else//_spaceScheme==staggered
2021         {
2022                 Vector Fij(_nVar);
2023                 Fij(_nVar-1)=0;
2024                 double alpha_roe = _Uroe[0];//Toumi formula
2025                 // interfacial pressure term (hyperbolic correction)
2026                 double dpi = _Uroe[_nVar];
2027
2028                 double u1ijn=0, u2ijn=0, phialphaq1n=0, phialphaq2n=0,vitesse1n=0,vitesse2n=0;
2029                 for(int idim=0;idim<_Ndim;idim++){//URoe = alpha, p, u1, u2, T, dpi
2030                         u1ijn+=_vec_normal[idim]*_Uroe[2+idim];
2031                         u2ijn+=_vec_normal[idim]*_Uroe[2+_Ndim+idim];
2032                 }
2033                 if(u1ijn>=0)
2034                 {
2035                         for(int idim=0;idim<_Ndim;idim++){
2036                                 phialphaq1n+=_vec_normal[idim]*_Ui[1+idim];//phi alpha rho u n
2037                                 vitesse1n  +=_vec_normal[idim]*_Vi[2+idim];
2038                         }
2039                         Fij(0)=phialphaq1n;
2040                         for(int idim=0;idim<_Ndim;idim++)
2041                                 Fij(1+idim)=phialphaq1n*_Vi[2+idim]+(alpha_roe*_Vj[1]*_porosityj+dpi*_Vi[0]*_porosityi)*_vec_normal[idim];
2042
2043                         double pressioni=_Vi[1];
2044                         double Temperaturei= _Vi[_nVar-1];
2045                         double rho1=_fluides[0]->getDensity(pressioni,Temperaturei);
2046                         double e1_int=_fluides[0]->getInternalEnergy(Temperaturei,rho1);
2047                         Fij(_nVar-1)+=_Ui[0]*(e1_int+0.5*vitesse1n*vitesse1n+_Vj[1]/rho1)*vitesse1n;
2048                 }
2049                 else
2050                 {
2051                         for(int idim=0;idim<_Ndim;idim++){
2052                                 phialphaq2n+=_vec_normal[idim]*_Uj[1+idim];//phi alpha rho u n
2053                                 vitesse1n  +=_vec_normal[idim]*_Vj[2+idim];
2054                         }
2055                         Fij(0)=phialphaq2n;
2056                         for(int idim=0;idim<_Ndim;idim++)
2057                                 Fij(1+idim)=phialphaq2n*_Vj[2+idim]+(alpha_roe*_Vi[1]*_porosityi+dpi*_Vj[0]*_porosityj)*_vec_normal[idim];
2058
2059                         double pressionj=_Vj[1];
2060                         double Temperaturej= _Vj[_nVar-1];
2061                         double rho1=_fluides[0]->getDensity(pressionj,Temperaturej);
2062                         double e1_int=_fluides[0]->getInternalEnergy(Temperaturej,rho1);
2063                         Fij(_nVar-1)+=_Uj[0]*(e1_int+0.5*vitesse1n*vitesse1n+_Vi[1]/rho1)*vitesse1n;
2064                 }
2065
2066                 if(u2ijn>=0)
2067                 {
2068                         for(int idim=0;idim<_Ndim;idim++){
2069                                 phialphaq2n+=_vec_normal[idim]*_Ui[2+_Ndim+idim];//phi alpha rho u n
2070                                 vitesse2n  +=_vec_normal[idim]*_Vi[2+idim+_Ndim];
2071                         }
2072                         Fij(1+_Ndim)=phialphaq2n;
2073                         for(int idim=0;idim<_Ndim;idim++)
2074                                 Fij(2+_Ndim+idim)=phialphaq2n*_Vi[2+_Ndim+idim]+((1-alpha_roe)*_Vj[1]*_porosityj-dpi*_Vi[0]*_porosityi)*_vec_normal[idim];
2075
2076                         double pressioni=_Vi[1];
2077                         double Temperaturei= _Vi[_nVar-1];
2078                         double rho2=_fluides[1]->getDensity(pressioni,Temperaturei);
2079                         double e2_int=_fluides[1]->getInternalEnergy(Temperaturei,rho2);
2080                         Fij(_nVar-1)+=_Ui[1+_Ndim]*(e2_int+0.5*vitesse2n*vitesse2n+_Vj[1]/rho2)*vitesse2n;
2081                 }
2082                 else
2083                 {
2084                         for(int idim=0;idim<_Ndim;idim++){
2085                                 phialphaq2n+=_vec_normal[idim]*_Uj[2+_Ndim+idim];//phi alpha rho u n
2086                                 vitesse2n  +=_vec_normal[idim]*_Vj[2+idim+_Ndim];
2087                         }
2088                         Fij(1+_Ndim)=phialphaq2n;
2089                         for(int idim=0;idim<_Ndim;idim++)
2090                                 Fij(2+_Ndim+idim)=phialphaq2n*_Vj[2+_Ndim+idim]+((1-alpha_roe)*_Vi[1]*_porosityi-dpi*_Vj[0]*_porosityj)*_vec_normal[idim];
2091
2092                         double pressionj=_Vj[1];
2093                         double Temperaturej= _Vj[_nVar-1];
2094                         double rho2=_fluides[1]->getDensity(pressionj,Temperaturej);
2095                         double e2_int=_fluides[1]->getInternalEnergy(Temperaturej,rho2);
2096                         Fij(_nVar-1)+=_Uj[1+_Ndim]*(e2_int+0.5*vitesse2n*vitesse2n+_Vi[1]/rho2)*vitesse2n;
2097                 }
2098                 return Fij;
2099         }
2100 }
2101
2102 void FiveEqsTwoFluid::applyVFRoeLowMachCorrections(bool isBord, string groupname)
2103 {
2104         if(_nonLinearFormulation!=VFRoe)
2105                 throw CdmathException("FiveEqsTwoFluid::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
2106         else//_nonLinearFormulation==VFRoe
2107         {
2108                 if(_spaceScheme==lowMach){
2109                         double u1_2=0, u2_2=0;
2110                         for(int i=0;i<_Ndim;i++){
2111                                 u1_2 += _Uroe[2+i]*_Uroe[2+i];
2112                                 u2_2 += _Uroe[2+i+_Ndim]*_Uroe[2+i+_Ndim];
2113                         }
2114
2115                         double  c = _maxvploc;//mixture sound speed
2116                         double M=max(sqrt(u1_2),sqrt(u2_2))/c;//Mach number
2117                         _Vij[1]=M*_Vij[1]+(1-M)*(_Vi[1]+_Vj[1])/2;
2118                         primToCons(_Vij,0,_Uij,0);
2119                 }
2120                 else if(_spaceScheme==pressureCorrection)
2121                 {
2122                         if(_pressureCorrectionOrder>2)
2123                                 throw CdmathException("FiveEqsTwoFluid::applyVFRoeLowMachCorrections pressure correction order can be only 1 or 2 for five equation two-fluid model");
2124
2125                         double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;//mean velocities
2126                         double rho1 =  _fluides[0]->getDensity(_Uroe[1],_Uroe[_nVar-1]);
2127                         double rho2 =  _fluides[1]->getDensity(_Uroe[1],_Uroe[_nVar-1]);
2128                         double m1=_Uroe[0]*rho1,  m2=(1-_Uroe[0])*rho2;
2129                         double rhom=m1+m2;
2130                         for(int i=0;i<_Ndim;i++)
2131                         {
2132                                 norm_uij += (m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim])*(m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim]);
2133                                 uij_n += (m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim])*_vec_normal[i];
2134                                 ui_n += _Vi[2+i]*_vec_normal[i];
2135                                 uj_n += _Vj[2+i]*_vec_normal[i];
2136                         }
2137                         norm_uij=sqrt(norm_uij)/rhom;
2138                         uij_n/=rhom;
2139                         if(norm_uij>_precision)//avoid division by zero
2140                                 _Vij[1]=(_Vi[1]+_Vj[1])/2 + uij_n/norm_uij*(_Vj[1]-_Vi[1])/4 - rhom*norm_uij*(uj_n-ui_n)/4;
2141                         else
2142                                 _Vij[1]=(_Vi[1]+_Vj[1])/2                                    - rhom*norm_uij*(uj_n-ui_n)/4;
2143                 }
2144                 else if(_spaceScheme==staggered)
2145                 {
2146                         double qij_n=0;
2147                         double rho1 =  _fluides[0]->getDensity(_Uroe[1],_Uroe[_nVar-1]);
2148                         double rho2 =  _fluides[1]->getDensity(_Uroe[1],_Uroe[_nVar-1]);
2149                         double m1=_Uroe[0]*rho1,  m2=(1-_Uroe[0])*rho2;
2150                         double rhom=m1+m2;
2151                         for(int i=0;i<_Ndim;i++)
2152                                 qij_n += (m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim])*_vec_normal[i];
2153
2154                         if(qij_n>=0){
2155                                 _Vij[0]=_Vi[0];
2156                                 _Vij[1]=_Vj[1];
2157                                 for(int i=0;i<_Ndim;i++)
2158                                 {
2159                                         _Vij[2+i]               =_Vi[2+i];
2160                                         _Vij[2+i+_Ndim] =_Vi[2+i+_Ndim];
2161                                 }
2162                                 _Vij[_nVar-1]=_Vi[_nVar-1];
2163                         }
2164                         else{
2165                                 _Vij[0]=_Vj[0];
2166                                 _Vij[1]=_Vi[1];
2167                                 for(int i=0;i<_Ndim;i++)
2168                                 {
2169                                         _Vij[2+i]               =_Vj[2+i];
2170                                         _Vij[2+i+_Ndim] =_Vj[2+i+_Ndim];
2171                                 }
2172                                 _Vij[_nVar-1]=_Vj[_nVar-1];
2173                         }
2174                         primToCons(_Vij,0,_Uij,0);
2175                 }
2176         }
2177 }
2178
2179 void FiveEqsTwoFluid::computeScaling(double maxvp)
2180 {
2181         _blockDiag[0]=1;//alphaScaling;
2182         _invBlockDiag[0]=1;//_blockDiag[0];
2183         _blockDiag[1+_Ndim]=1;//-alphaScaling;
2184         _invBlockDiag[1+_Ndim]=1.0;//_blockDiag[1+_Ndim];
2185         for(int q=1; q<_Ndim+1; q++)
2186         {
2187                 _blockDiag[q]=1/maxvp;
2188                 _invBlockDiag[q]=1/_blockDiag[q];
2189                 _blockDiag[q+1+_Ndim]=1/maxvp;
2190                 _invBlockDiag[q+1+_Ndim]=1/_blockDiag[q+1+_Ndim];
2191         }
2192         _blockDiag[_nVar - 1]=1/(maxvp*maxvp);//1
2193         _invBlockDiag[_nVar - 1]=  1./_blockDiag[_nVar - 1] ;// 1.;//
2194 }
2195
2196 void FiveEqsTwoFluid::testConservation()
2197 {
2198         double SUM, DELTA, x;
2199         int I;
2200         for(int i=0; i<_nVar; i++)
2201         {
2202                 {
2203                         if(i == 0)
2204                                 cout << "Masse totale phase " << 0 <<" (kg): ";
2205                         else if( i == 1+_Ndim)
2206                                 cout << "Masse totale phase " << 1 <<" (kg): ";
2207                         else if( i == _nVar-1)
2208                                 cout << "Energie totale " <<" (J.m^-3): ";
2209                         else
2210                         {
2211                                 if(i < 1+_Ndim)
2212                                         cout << "Quantite de mouvement totale phase 0 (kg.m.s^-1): ";
2213                                 else
2214                                         cout << "Quantite de mouvement totale phase 1 (kg.m.s^-1): ";
2215                         }
2216                 }
2217                 SUM = 0;
2218                 DELTA = 0;
2219                 I =  i;
2220                 for(int j=0; j<_Nmailles; j++)
2221                 {
2222                         VecGetValues(_old, 1, &I, &x);//on recupere la valeur du champ
2223                         SUM += x*_mesh.getCell(j).getMeasure();
2224                         VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
2225                         DELTA += x*_mesh.getCell(j).getMeasure();
2226                         I += _nVar;
2227                 }
2228                 if(fabs(SUM)>_precision)
2229                         cout << SUM << ", variation relative: " << fabs(DELTA /SUM)  << endl;
2230                 else
2231                         cout << " a une somme nulle,  variation absolue: " << fabs(DELTA) << endl;
2232         }
2233 }
2234
2235 void FiveEqsTwoFluid::save(){
2236     PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results at time step number %d \n\n", _nbTimeStep);
2237     *_runLogFile<< "Saving numerical results at time step number "<< _nbTimeStep << endl<<endl;
2238
2239         string prim(_path+"/FiveEqsTwoFluidPrim_");
2240         string cons(_path+"/FiveEqsTwoFluidCons_");
2241         prim+=_fileName;
2242         cons+=_fileName;
2243
2244         PetscInt Ii, lj;
2245         for (PetscInt i = 0; i < _Nmailles; i++){
2246                 /* j = 0 : void fraction
2247                            j = 1 : pressure
2248                            j = 2, 3, 4: velocity phase 1
2249                            j = 5, 6, 7: velocity phase 2
2250                            j = 8 : temperature */
2251                 for (int j = 0; j < _nVar; j++){
2252                         Ii = i*_nVar +j;
2253                         VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
2254                 }
2255         }
2256         if(_saveConservativeField){
2257                 for (long i = 0; i < _Nmailles; i++){
2258                         for (int j = 0; j < _nVar; j++){
2259                                 Ii = i*_nVar +j;
2260                                 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
2261                         }
2262                 }
2263                 _UU.setTime(_time,_nbTimeStep+1);
2264         }
2265         _VV.setTime(_time,_nbTimeStep+1);
2266
2267         if (_nbTimeStep ==0 || _restartWithNewFileName){
2268                 string prim_suppress ="rm -rf "+prim+"_*";
2269                 string cons_suppress ="rm -rf "+cons+"_*";
2270                 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
2271                 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
2272                 _VV.setInfoOnComponent(0,"Void_fraction");
2273                 _VV.setInfoOnComponent(1,"Pressure_(Pa)");
2274                 _VV.setInfoOnComponent(2,"Velocity1_x_m/s");
2275
2276                 if (_Ndim>1)
2277                         _VV.setInfoOnComponent(3,"Velocity1_y_m/s");
2278                 if (_Ndim>2)
2279                         _VV.setInfoOnComponent(4,"Velocity1_z_m/s");
2280                 _VV.setInfoOnComponent(2+_Ndim,"Velocity2_x_m/s");
2281                 if (_Ndim>1)
2282                         _VV.setInfoOnComponent(3+_Ndim,"Velocity2_y_m/s");
2283                 if (_Ndim>2)
2284                         _VV.setInfoOnComponent(4+_Ndim,"Velocity2_z_m/s");
2285                 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
2286
2287                 switch(_saveFormat)
2288                 {
2289                 case VTK :
2290                         _VV.writeVTK(prim);
2291                         break;
2292                 case MED :
2293                         _VV.writeMED(prim);
2294                         break;
2295                 case CSV :
2296                         _VV.writeCSV(prim);
2297                         break;
2298                 }
2299                 if(_saveConservativeField){
2300                         _UU.setInfoOnComponent(0,"Partial_density1");// (kg/m^3)
2301                         _UU.setInfoOnComponent(1,"Momentum1_x");// phase1 (kg/m^2/s)
2302                         if (_Ndim>1)
2303                                 _UU.setInfoOnComponent(2,"Momentum1_y");// phase1 (kg/m^2/s)
2304                         if (_Ndim>2)
2305                                 _UU.setInfoOnComponent(3,"Momentum1_z");// phase1  (kg/m^2/s)
2306                         _UU.setInfoOnComponent(1+_Ndim,"Partial_density2");// phase2 (kg/m^3)
2307                         _UU.setInfoOnComponent(2+_Ndim,"Momentum2_x");// phase2 (kg/m^2/s)
2308
2309                         if (_Ndim>1)
2310                                 _UU.setInfoOnComponent(3+_Ndim,"Momentum2_y");// phase2 (kg/m^2/s)
2311                         if (_Ndim>2)
2312                                 _UU.setInfoOnComponent(4+_Ndim,"Momentum2_z");// phase2 (kg/m^2/s)
2313                         _UU.setInfoOnComponent(_nVar-1,"Total_energy");
2314
2315                         switch(_saveFormat)
2316                         {
2317                         case VTK :
2318                                 _UU.writeVTK(cons);
2319                                 break;
2320                         case MED :
2321                                 _UU.writeMED(cons);
2322                                 break;
2323                         case CSV :
2324                                 _UU.writeCSV(cons);
2325                                 break;
2326                         }
2327                 }
2328         }
2329         else{
2330                 switch(_saveFormat)
2331                 {
2332                 case VTK :
2333                         _VV.writeVTK(prim,false);
2334                         break;
2335                 case MED :
2336                         _VV.writeMED(prim,false);
2337                         break;
2338                 case CSV :
2339                         _VV.writeCSV(prim);
2340                         break;
2341                 }
2342
2343                 if(_saveConservativeField){
2344                         switch(_saveFormat)
2345                         {
2346                         case VTK :
2347                                 _UU.writeVTK(cons,false);
2348                                 break;
2349                         case MED :
2350                                 _UU.writeMED(cons,false);
2351                                 break;
2352                         case CSV :
2353                                 _UU.writeCSV(cons);
2354                                 break;
2355                         }
2356                 }
2357         }
2358         if(_saveVelocity){
2359                 for (long i = 0; i < _Nmailles; i++){
2360                         // j = 0 : concentration, j=1 : pressure; j = _nVar - 1: temperature; j = 2,..,_nVar-2: velocity
2361                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
2362                         {
2363                                 Ii = i*_nVar +2+j;
2364                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse1(i,j));
2365                                 Ii=i*_nVar +2+j+_Ndim;
2366                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse2(i,j));
2367                         }
2368                         for (int j = _Ndim; j < 3; j++){//On met à zero les composantes de vitesse si la dimension est <3
2369                                 _Vitesse1(i,j)=0;
2370                                 _Vitesse2(i,j)=0;
2371                         }
2372                 }
2373                 _Vitesse1.setTime(_time,_nbTimeStep);
2374                 _Vitesse2.setTime(_time,_nbTimeStep);
2375                 if (_nbTimeStep ==0 || _restartWithNewFileName){
2376                         _Vitesse1.setInfoOnComponent(0,"Velocity_x_(m/s)");
2377                         _Vitesse1.setInfoOnComponent(1,"Velocity_y_(m/s)");
2378                         _Vitesse1.setInfoOnComponent(2,"Velocity_z_(m/s)");
2379
2380                         _Vitesse2.setInfoOnComponent(0,"Velocity_x_(m/s)");
2381                         _Vitesse2.setInfoOnComponent(1,"Velocity_y_(m/s)");
2382                         _Vitesse2.setInfoOnComponent(2,"Velocity_z_(m/s)");
2383
2384                         switch(_saveFormat)
2385                         {
2386                         case VTK :
2387                                 _Vitesse1.writeVTK(prim+"_GasVelocity");
2388                                 _Vitesse2.writeVTK(prim+"_LiquidVelocity");
2389                                 break;
2390                         case MED :
2391                                 _Vitesse1.writeMED(prim+"_GasVelocity");
2392                                 _Vitesse2.writeMED(prim+"_LiquidVelocity");
2393                                 break;
2394                         case CSV :
2395                                 _Vitesse1.writeCSV(prim+"_GasVelocity");
2396                                 _Vitesse2.writeCSV(prim+"_LiquidVelocity");
2397                                 break;
2398                         }
2399                 }
2400                 else{
2401                         switch(_saveFormat)
2402                         {
2403                         case VTK :
2404                                 _Vitesse1.writeVTK(prim+"_GasVelocity",false);
2405                                 _Vitesse2.writeVTK(prim+"_LiquidVelocity",false);
2406                                 break;
2407                         case MED :
2408                                 _Vitesse1.writeMED(prim+"_GasVelocity",false);
2409                                 _Vitesse2.writeMED(prim+"_LiquidVelocity",false);
2410                                 break;
2411                         case CSV :
2412                                 _Vitesse1.writeCSV(prim+"_GasVelocity");
2413                                 _Vitesse2.writeCSV(prim+"_LiquidVelocity");
2414                                 break;
2415                         }
2416                 }
2417         }
2418
2419         if (_restartWithNewFileName)
2420                 _restartWithNewFileName=false;
2421 }