CoreFlows/Models/src/IsothermalTwoFluid.cxx

   1 /*
   2  * IsothermalTwoFluid.cxx
   3  *
   4  *  Created on: Sep 16, 2014
   5  *      Author: tn236279
   6  */
   7
   8 #include "IsothermalTwoFluid.hxx"
   9
  10 using namespace std;
  11
  12 IsothermalTwoFluid::IsothermalTwoFluid(pressureEstimate pEstimate, int dim){
  13         _Ndim=dim;
  14         _nVar=2*(_Ndim+1);
  15         _nbPhases = 2;
  16         _dragCoeffs=vector<double>(2,0);
  17         _fluides.resize(2);
  18         if (pEstimate==around1bar300K)//EOS at 1 bar and 300K
  19         {
  20                 cout<<"Fluid is air-water mixture around 1 bar and 300 K (27°C)"<<endl;
  21                 *_runLogFile<<"Fluid is air-water mixture around 1 bar and 300 K (27°C)"<<endl;
  22                 _Temperature=300;//Constant temperature of the model
  23                 _internalEnergy1=2.22e5;//nitrogen internal energy at 1bar, 300K
  24                 _internalEnergy2=1.12e5;//water internal energy at 1 bar, 300K
  25                 _fluides[0] = new StiffenedGas(1.4,743,_Temperature,_internalEnergy1);  //ideal gas law for nitrogen at pressure 1 bar and temperature 27°C, c_v=743
  26                 _fluides[1] = new StiffenedGas(996,1e5,_Temperature,_internalEnergy2,1501,4130);  //stiffened gas law for water at pressure 1 bar and temperature 27°C
  27         }
  28         else//EOS at 155 bars and 618K
  29         {
  30                 cout<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
  31                 *_runLogFile<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
  32                 _Temperature=618;//Constant temperature of the model
  33                 _internalEnergy1=2.44e6;//Gas internal energy at saturation at 155 bar
  34                 _internalEnergy2=1.6e6;//water internal energy at saturation at 155 bar
  35                 _fluides[0] = new StiffenedGas(102,1.55e7,_Temperature,_internalEnergy1, 433,3633);  //stiffened gas law for Gas at pressure 155 bar and temperature 345°C:
  36                 _fluides[1] = new StiffenedGas(594,1.55e7,_Temperature,_internalEnergy2, 621,3100);  //stiffened gas law for water at pressure 155 bar and temperature 345°C:
  37         }
  38         _intPressCoeff=1.5;
  39 }
  40
  41 void IsothermalTwoFluid::initialize(){
  42         cout<<"Initialising the isothermal two-fluid model"<<endl;
  43         *_runLogFile<<"Initialising the isothermal two-fluid model"<<endl;
  44
  45         _Uroe = new double[_nVar+1];
  46
  47         _guessalpha = _VV(0,0);
  48
  49         _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
  50         for(int i=0; i<_Ndim; i++)
  51         {
  52                 _gravite[i+1]=_GravityField3d[i];
  53                 _gravite[i+1 +_Ndim+1]=_GravityField3d[i];
  54         }
  55         _GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
  56
  57         if(_saveVelocity){
  58                 _Vitesse1=Field("Gas velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
  59                 _Vitesse2=Field("Liquid velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
  60         }
  61         if(_entropicCorrection)
  62                 _entropicShift=vector<double>(3,0);
  63
  64         ProblemFluid::initialize();
  65 }
  66
  67 void IsothermalTwoFluid::convectionState( const long &i, const long &j, const bool &IsBord){
  68         //sortie: WRoe en (alpha, p, u1, u2, dm1,dm2,dalpha1,dp)z
  69         //entree: _conservativeVars en (alpha1 rho1, alpha1 rho1 u1, alpha2 rho2, alpha2 rho2 u2)
  70
  71         // _l always inside the domain (index i)
  72         // _r is maybe the boundary cell (negative index)
  73         // _l and _r are primative vectors (alp, P, u1, u2)
  74
  75         _idm[0] = _nVar*i;
  76         for(int k=1; k<_nVar; k++)
  77                 _idm[k] = _idm[k-1] + 1;
  78         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
  79
  80         _idm[0] = _nVar*j;
  81         for(int k=1; k<_nVar; k++)
  82                 _idm[k] = _idm[k-1] + 1;
  83         if(IsBord)
  84                 VecGetValues(_Uext, _nVar, _idm, _Uj);
  85         else
  86                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
  87         if(_verbose && _nbTimeStep%_freqSave ==0)
  88         {
  89                 cout<<"Convection Left state cell " << i<< ": "<<endl;
  90                 for(int k =0; k<_nVar; k++)
  91                         cout<< _Ui[k]<<endl;
  92                 cout<<"Convection Right state cell " << j<< ": "<<endl;
  93                 for(int k =0; k<_nVar; k++)
  94                         cout<< _Uj[k]<<endl;
  95         }
  96         //       if(_Ui[0]<-(_precision) || _Uj[0]<-(_precision) || _Ui[_Ndim+1]<-(_precision) || _Uj[_Ndim+1]<-(_precision))
  97         //      {
  98         //        cout<<"!!!!!!!!!!!!!!!!!!!!!!!! Masse partielle negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
  99         //        cout<< "valeurs a gauche: "<<_Ui[0]<<", "<<_Ui[_Ndim+1]<<"valeurs a droite: "<<_Uj[0]<<", "<<_Uj[_Ndim+1]<<endl;
 100         //        throw CdmathException(" Masse partielle negative, arret de calcul");
 101         //      }
 102         //       else
 103         {
 104                 _Ui[0]=max(0.,_Ui[0]);// mass1 a gauche
 105                 _Uj[0]=max(0.,_Uj[0]);// mass1 a droite
 106                 _Ui[_Ndim+1]=max(0.,_Ui[_Ndim+1]);// mass2 a gauche
 107                 _Uj[_Ndim+1]=max(0.,_Uj[_Ndim+1]);// mass2 a droite
 108         }
 109
 110         PetscScalar ri1, ri2, rj1, rj2, xi, xj;
 111         _idm[0] = _nVar*i;
 112         for(int k=1; k<_nVar; k++)
 113                 _idm[k] = _idm[k-1] + 1;
 114         VecGetValues(_primitiveVars, _nVar, _idm, _l);
 115
 116         if(_verbose && _nbTimeStep%_freqSave ==0)
 117         {
 118                 cout<<"Etat de Roe, etat primitif gauche: "<<endl;
 119                 for(int i =0; i<_nVar; i++)
 120                         cout<< _l[i]<<endl;
 121         }
 122
 123         // boundary : compute the _r
 124         if(IsBord)
 125         {
 126                 _guessalpha=_l[0];
 127                 consToPrim(_Uj, _r);
 128         }
 129         // inside the domain : extract from the primative vector
 130         else
 131         {
 132                 _idm[0] = _nVar*j;
 133                 for(int k=1; k<_nVar; k++)
 134                         _idm[k] = _idm[k-1] + 1;
 135                 VecGetValues(_primitiveVars, _nVar, _idm, _r);
 136         }
 137
 138         if(_verbose && _nbTimeStep%_freqSave ==0)
 139         {
 140                 cout<<"Etat de Roe, etat primitif droite: "<<endl;
 141                 for(int i =0; i<_nVar; i++)
 142                         cout<< _r[i]<<endl;
 143         }
 144
 145         // Using Toumi linearisation (read in the article of Toumi) : alpha^{Roe}_l
 146         if(2-_l[0]-_r[0] > _precision)
 147                 _Uroe[0] = 1- 2*(1-_l[0])*(1-_r[0])/(2-_l[0]-_r[0]);
 148         // Using an average : (alp_l + alp_r)/2 : suggestion of Michael (no theory)
 149         else
 150                 _Uroe[0] = (_l[0]+_r[0])/2;
 151         // Pressure is computed as function of alp and P_l, P_r (Toumi article)
 152         if(_l[0]+_r[0] > _precision)
 153                 _Uroe[1] = (_l[1]*_l[0]+_r[1]*_r[0])/(_l[0]+_r[0]);
 154         else
 155                 _Uroe[1] = (_l[1]*(1-_l[0])+_r[1]*(1-_r[0]))/(2-_l[0]-_r[0]);
 156         // i :left, j : right (U is normally conservative variable)
 157         ri1 = sqrt(_Ui[0]); ri2 = sqrt(_Ui[_Ndim+1]);
 158         rj1 = sqrt(_Uj[0]); rj2 = sqrt(_Uj[_Ndim+1]);
 159         // Roe average formula of the velocities
 160         for(int k=0;k<_Ndim;k++)
 161         {
 162                 xi = _Ui[k+1];
 163                 xj = _Uj[k+1];
 164                 // avoid dividing by zero, if mass is zero, do not consider the distribution of such a phase
 165                 if(ri1>_precision && rj1>_precision)
 166                         _Uroe[2+k] = (xi/ri1 + xj/rj1)/(ri1 + rj1);
 167                 else if(ri1<_precision && rj1>_precision)
 168                         _Uroe[2+k] =  xj/_Uj[0];
 169                 else if(ri1>_precision && rj1<_precision)
 170                         _Uroe[2+k] =  xi/_Ui[0];
 171                 else
 172                         _Uroe[2+k] =(_Ui[k+1+_Ndim+1]/ri2 + _Uj[k+1+_Ndim+1]/rj2)/(ri2 + rj2);
 173
 174                 xi = _Ui[k+1+_Ndim+1];
 175                 xj = _Uj[k+1+_Ndim+1];
 176                 if(ri2>_precision && rj2>_precision)
 177                         _Uroe[2+k+_Ndim] = (xi/ri2 + xj/rj2)/(ri2 + rj2);
 178                 else if(ri2<_precision && rj2>_precision)
 179                         _Uroe[2+k+_Ndim] = xj/_Uj[_Ndim+1];
 180                 else if(ri2>_precision && rj2<_precision)
 181                         _Uroe[2+k+_Ndim] = xi/_Ui[_Ndim+1];
 182                 else
 183                         _Uroe[2+k+_Ndim] = (xi/ri1 + xj/rj1)/(ri1 + rj1);
 184         }
 185         if(_verbose && _nbTimeStep%_freqSave ==0)
 186         {
 187                 cout<<"Etat de Roe calcule: "<<endl;
 188                 for(int k=0;k<_nVar; k++)//At this point _Uroe[_nVar] is not yet set
 189                         cout<< _Uroe[k]<<endl;
 190         }
 191 }
 192
 193 void IsothermalTwoFluid::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
 194         //sortie: matrices et etat de diffusion (m_v, q_v, m_l, q_l)
 195         _idm[0] = _nVar*i;
 196         for(int k=1; k<_nVar; k++)
 197                 _idm[k] = _idm[k-1] + 1;
 198
 199         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 200         _idm[0] = _nVar*j;
 201         for(int k=1; k<_nVar; k++)
 202                 _idm[k] = _idm[k-1] + 1;
 203
 204         if(IsBord)
 205                 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
 206         else
 207                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 208
 209         for(int k=0; k<_nVar; k++)
 210                 _Udiff[k] = (_Ui[k]+_Uj[k])/2;
 211         double q_2=0;
 212         for (int i = 0; i<_Ndim;i++)
 213                 q_2+=_Udiff[i+1]*_Udiff[i+1];
 214
 215         if(_timeScheme==Implicit)
 216         {
 217                 for(int i=0; i<_nVar*_nVar;i++)
 218                         _Diffusion[i] = 0;
 219                 double mu1 = _fluides[0]->getViscosity(_Temperature);
 220                 double mu2 = _fluides[1]->getViscosity(_Temperature);
 221                 for(int idim=1;idim<_Ndim+1;idim++)
 222                 {
 223                         _Diffusion[idim*_nVar] =  mu1*_Udiff[idim]/(_Udiff[0]*_Udiff[0]);
 224                         _Diffusion[idim*_nVar+idim] = -mu1/_Udiff[0];
 225                         _Diffusion[(idim+_Ndim+1)*_nVar] =  mu2*_Udiff[idim+_Ndim+1]/(_Udiff[_Ndim+1]*_Udiff[_Ndim+1]);
 226                         _Diffusion[(idim+_Ndim+1)*_nVar+idim+_Ndim+1] = -mu2/_Udiff[_Ndim+1];
 227                 }
 228         }
 229 }
 230
 231 void IsothermalTwoFluid::convectionMatrices()
 232 {
 233         //entree: URoe = alpha, p, u1, u2 + ajout dpi pour calcul flux ultérieur
 234         //sortie: matrices Roe+  et Roe- +Roe si scheme centre
 235
 236         if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)
 237                 throw CdmathException("Implicitation with primitive variables not yet available for IsothermalTwoFluid model");
 238
 239         /*Definitions */
 240         complex< double > tmp;
 241         double u1_n, u1_2, u2_n, u2_2, u_r2;
 242         u1_2 = 0; u2_2=0;
 243         u1_n = 0; u2_n=0;
 244         // relative velocity
 245         u_r2=0;
 246         for(int i=0;i<_Ndim;i++)
 247         {
 248                 u1_2 += _Uroe[2+i]*_Uroe[2+i];
 249                 u1_n += _Uroe[2+i]*_vec_normal[i];
 250                 u2_2 += _Uroe[2+i+_Ndim]*_Uroe[2+i+_Ndim];
 251                 u2_n += _Uroe[2+i+_Ndim]*_vec_normal[i];
 252                 u_r2 += (_Uroe[2+i]-_Uroe[2+i+_Ndim])*(_Uroe[2+i]-_Uroe[2+i+_Ndim]);
 253         }
 254         //Ancienne construction Mat Roe (Dm1,Dm2,Dalp,Dp)
 255         double alpha = _Uroe[0];
 256         double p = _Uroe[1];
 257         double rho1 = _fluides[0]->getDensity(p, _Temperature);
 258         double rho2 = _fluides[1]->getDensity(p, _Temperature);
 259         double dpi1 = intPressDef(alpha, u_r2, rho1, rho2);
 260         double dpi2 = dpi1;
 261         double invcsq1 = 1/_fluides[0]->vitesseSonTemperature(_Temperature,rho1);
 262         invcsq1*=invcsq1;
 263         double invcsq2 = 1/_fluides[1]->vitesseSonTemperature(_Temperature,rho2);
 264         invcsq2*=invcsq2;
 265         double g2 = 1/(alpha*rho2*invcsq1+(1-alpha)*rho1*invcsq2);
 266         double g2press=g2, g2alpha=g2;
 267         //saving dpi value for flux calculation later
 268         _Uroe[_nVar]=dpi1 ;
 269
 270         /***********Calcul des valeurs propres ********/
 271         vector< double > pol_car= Polynoms::polynome_caracteristique(alpha,  1-alpha, u1_n, u2_n, rho1, rho2,invcsq1,invcsq2, dpi1, dpi2, g2press, g2alpha, g2,_precision);
 272         for(int ct=0;ct<4;ct++){
 273                 if (abs(pol_car[ct])<_precision*_precision)
 274                         pol_car[ct]=0;
 275         }
 276         if(_verbose && _nbTimeStep%_freqSave ==0)
 277         {
 278                 cout<<"pol caract= "<<endl;
 279                 for(int i =0; i<5; i++)
 280                         cout<<pol_car[i]<<"  ";
 281                 cout<< endl;
 282                 cout<<"alpha= "<<alpha<<", p= " << p << ", rho1= " << rho1<< ", rho2= " << rho2<< ", c1= " <<sqrt(1/invcsq1)<<
 283                                 ", c2= " <<sqrt(1/invcsq2)<<endl;
 284                 cout<< "u1_n= "<<u1_n<<", u2_n= "<<u2_n<< ", dpi1= "<<dpi1<< ", dpi2= "<<dpi2<< endl;
 285         }
 286         vector< complex<double> > valeurs_propres = getRacines(pol_car);
 287
 288         //On ajoute les valeurs propres triviales
 289         if(_Ndim>1)
 290         {
 291                 if( !Polynoms::belongTo(u1_n,valeurs_propres, _precision) )
 292                         valeurs_propres.push_back(u1_n);//vp vapor energy
 293                 if( !Polynoms::belongTo(u2_n,valeurs_propres, _precision) )
 294                         valeurs_propres.push_back(u2_n);//vp liquid energy
 295         }
 296         bool doubleeigenval = norm(valeurs_propres[0] - valeurs_propres[1])<_precision;//norm= suqare of the magnitude
 297         if(doubleeigenval)
 298         {
 299                 valeurs_propres[0] = valeurs_propres[valeurs_propres.size()-1];
 300                 valeurs_propres.pop_back();
 301         }
 302
 303         int taille_vp = Polynoms::new_tri_selectif(valeurs_propres,valeurs_propres.size(),_precision);//valeurs_propres.size();//
 304
 305         _maxvploc=0;
 306         for(int i =0; i<taille_vp; i++)
 307                 if(fabs(valeurs_propres[i].real())>_maxvploc)
 308                         _maxvploc=fabs(valeurs_propres[i].real());
 309         if(_maxvploc>_maxvp)
 310                 _maxvp=_maxvploc;
 311
 312         int existVpCplx = 0,pos_conj;
 313         double vp_imag_iter;
 314         for (int ct=0; ct<taille_vp; ct++) {
 315                 vp_imag_iter = valeurs_propres[ct].imag();
 316                 if ( fabs(vp_imag_iter) > _precision ) {
 317                         existVpCplx +=1;
 318                         if ( _part_imag_max < fabs(vp_imag_iter))
 319                                 _part_imag_max = fabs(vp_imag_iter);
 320                         //On cherhe le conjugue
 321                         pos_conj = ct+1;
 322                         while(pos_conj<taille_vp && fabs(valeurs_propres[pos_conj].imag()+vp_imag_iter)>_precision)
 323                                 pos_conj++;
 324                         if(pos_conj!=ct+1 && pos_conj<taille_vp )
 325                         {
 326                                 tmp=valeurs_propres[ct+1];
 327                                 valeurs_propres[ct+1]=valeurs_propres[pos_conj];
 328                                 valeurs_propres[pos_conj] = tmp;
 329                                 ct++;
 330                         }
 331                 }
 332         }
 333         if (existVpCplx >0)
 334                 _nbVpCplx +=1;
 335
 336         //on ordonne les deux premieres valeurs
 337         /*
 338         if(valeurs_propres[1].real()<valeurs_propres[0].real())
 339         {
 340                 tmp=valeurs_propres[0];
 341                 valeurs_propres[0]=valeurs_propres[1];
 342                 valeurs_propres[1]=tmp;
 343         }
 344          */
 345         if(_verbose && _nbTimeStep%_freqSave ==0)
 346         {
 347                 cout<<" Vp apres tri " << valeurs_propres.size()<<endl;
 348                 for(int ct =0; ct<taille_vp; ct++)
 349                         cout<< "("<<valeurs_propres[ct].real()<< ", " <<valeurs_propres[ct].imag() <<")  ";
 350                 cout<< endl;
 351         }
 352
 353         /******** Construction de la matrice de Roe *********/
 354         //lignes de masse
 355         for(int i=0; i<_nVar*_nVar;i++)
 356                 _Aroe[i]=0;
 357
 358         for(int idim=0; idim<_Ndim;idim++)
 359         {
 360                 _Aroe[1+idim]=_vec_normal[idim];
 361                 _Aroe[1+idim+_Ndim+1]=0;
 362                 _Aroe[(_Ndim+1)*_nVar+1+idim]=0;
 363                 _Aroe[(_Ndim+1)*_nVar+1+idim+_Ndim+1]=_vec_normal[idim];
 364         }
 365         //lignes de qdm
 366         for(int idim=0; idim<_Ndim;idim++)
 367         {
 368                 //premiere colonne (masse gaz)
 369                 _Aroe[                 (1+idim)*_nVar]=   (alpha *rho2*g2press+dpi1*(1-alpha)*invcsq2*g2alpha)*_vec_normal[idim] - u1_n*_Uroe[2+idim];
 370                 _Aroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar]=((1-alpha)*rho2*g2press-dpi2*(1-alpha)*invcsq2*g2alpha)*_vec_normal[idim];
 371                 //colonnes intermediaires
 372                 for(int jdim=0; jdim<_Ndim;jdim++)
 373                 {
 374                         _Aroe[                 (1+idim)*_nVar + jdim + 1]         = _Uroe[      2+idim]*_vec_normal[jdim];
 375                         _Aroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar + jdim + 1+_Ndim+1] = _Uroe[_Ndim+2+idim]*_vec_normal[jdim];
 376                 }
 377                 //matrice identite
 378                 _Aroe[                 (1+idim)*_nVar +          idim + 1] += u1_n;
 379                 _Aroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar + _Ndim+1+ idim + 1] += u2_n;
 380                 //troisieme colonne (masse liquide)
 381                 _Aroe[                (1+idim)*_nVar + _Ndim+1]=   (alpha *rho1*g2press-dpi1*alpha*invcsq1*g2alpha)*_vec_normal[idim];
 382                 _Aroe[(_Ndim+1)*_nVar+(1+idim)*_nVar + _Ndim+1]=((1-alpha)*rho1*g2press+dpi2*alpha*invcsq1*g2alpha)*_vec_normal[idim] - u2_n*_Uroe[1+idim+_Ndim+1];
 383         }
 384
 385         /******* Construction des matrices de decentrement *****/
 386         if(_spaceScheme == centered){
 387                 if(_entropicCorrection)
 388                         throw CdmathException("IsothermalTwoFluid::convectionMatrices: entropic scheme not available for centered scheme");
 389                 for(int i=0; i<_nVar*_nVar;i++)
 390                         _absAroe[i]=0;
 391         }
 392         if( _spaceScheme ==staggered){
 393                 if(_entropicCorrection)//To do: study entropic correction for staggered
 394                         throw CdmathException("IsothermalTwoFluid::convectionMatrices: entropic scheme not yet available for staggered scheme");
 395                 /******** Construction du decentrement du type decale *********/
 396                 //lignes de masse
 397                 for(int i=0; i<_nVar*_nVar;i++)
 398                         _absAroe[i]=0;
 399
 400                 for(int idim=0; idim<_Ndim;idim++)
 401                 {
 402                         _absAroe[1+idim]=_vec_normal[idim];
 403                         _absAroe[1+idim+_Ndim+1]=0;
 404                         _absAroe[(_Ndim+1)*_nVar+1+idim]=0;
 405                         _absAroe[(_Ndim+1)*_nVar+1+idim+_Ndim+1]=_vec_normal[idim];
 406                 }
 407                 //lignes de qdm
 408                 for(int idim=0; idim<_Ndim;idim++)
 409                 {
 410                         //premiere colonne (masse gaz)
 411                         _absAroe[                 (1+idim)*_nVar]=   (-alpha *rho2*g2press+dpi1*(1-alpha)*invcsq2*g2alpha)*_vec_normal[idim] - u1_n*_Uroe[2+idim];
 412                         _absAroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar]=(-(1-alpha)*rho2*g2press-dpi2*(1-alpha)*invcsq2*g2alpha)*_vec_normal[idim];
 413                         //colonnes intermediaires
 414                         for(int jdim=0; jdim<_Ndim;jdim++)
 415                         {
 416                                 _absAroe[                 (1+idim)*_nVar + jdim + 1]         = _Uroe[      2+idim]*_vec_normal[jdim];
 417                                 _absAroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar + jdim + 1+_Ndim+1] = _Uroe[_Ndim+2+idim]*_vec_normal[jdim];
 418                         }
 419                         //matrice identite
 420                         _absAroe[                 (1+idim)*_nVar +          idim + 1] += u1_n;
 421                         _absAroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar + _Ndim+1+ idim + 1] += u2_n;
 422                         //troisieme colonne (masse liquide)
 423                         _absAroe[                (1+idim)*_nVar + _Ndim+1]=   (-alpha *rho1*g2press-dpi1*alpha*invcsq1*g2alpha)*_vec_normal[idim];
 424                         _absAroe[(_Ndim+1)*_nVar+(1+idim)*_nVar + _Ndim+1]=(-(1-alpha)*rho1*g2press+dpi2*alpha*invcsq1*g2alpha)*_vec_normal[idim] - u2_n*_Uroe[1+idim+_Ndim+1];
 425                 }
 426
 427                 double signu1,signu2;
 428                 if(u1_n>0)
 429                         signu1=1;
 430                 else if (u1_n<0)
 431                         signu1=-1;
 432                 else
 433                         signu1=0;
 434                 if(u2_n>0)
 435                         signu2=1;
 436                 else if (u2_n<0)
 437                         signu2=-1;
 438                 else
 439                         signu2=0;
 440
 441                 for(int i=0; i<(1+_Ndim)*_nVar;i++){
 442                         _absAroe[i] *= signu1;
 443                         _absAroe[i+(1+_Ndim)*_nVar] *= signu2;
 444                 }
 445         }
 446         if(_spaceScheme==upwind || _spaceScheme ==lowMach)
 447         {
 448                 vector< complex< double > > y (taille_vp,0);
 449                 for( int i=0 ; i<taille_vp ; i++)
 450                         y[i] = Polynoms::abs_generalise(valeurs_propres[i]);
 451
 452                 if(_entropicCorrection)
 453                 {
 454                         entropicShift(_vec_normal);
 455                         y[0] +=_entropicShift[0];
 456                         y[taille_vp-1] +=_entropicShift[2];
 457                         for( int i=1 ; i<taille_vp-1 ; i++)
 458                                 y[i] +=_entropicShift[1];
 459                 }
 460
 461                 Polynoms::abs_par_interp_directe(taille_vp,valeurs_propres, _Aroe, _nVar,_precision, _absAroe,y);
 462                 if( _spaceScheme ==pressureCorrection){
 463                         for( int i=0 ; i<_Ndim ; i++)
 464                                 for( int j=0 ; j<_Ndim ; j++){
 465                                         _absAroe[(1+i)*_nVar+1+j]-=alpha*(valeurs_propres[1].real()-valeurs_propres[0].real())/2*_vec_normal[i]*_vec_normal[j];
 466                                         _absAroe[(2+_Ndim+i)*_nVar+2+_Ndim+j]-=(1-alpha)*(valeurs_propres[1].real()-valeurs_propres[0].real())/2*_vec_normal[i]*_vec_normal[j];
 467                                 }
 468                 }
 469                 else if( _spaceScheme ==lowMach){
 470                         double M=max(fabs(u1_2),fabs(u2_2))/_maxvploc;
 471                         for( int i=0 ; i<_Ndim ; i++)
 472                                 for( int j=0 ; j<_Ndim ; j++){
 473                                         _absAroe[(1+i)*_nVar+1+j]-=(1-M)*alpha*(valeurs_propres[1].real()-valeurs_propres[0].real())/2*_vec_normal[i]*_vec_normal[j];
 474                                         _absAroe[(2+_Ndim+i)*_nVar+2+_Ndim+j]-=(1-M)*(1-alpha)*(valeurs_propres[1].real()-valeurs_propres[0].real())/2*_vec_normal[i]*_vec_normal[j];
 475                                 }
 476                 }
 477
 478         }
 479         //Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source
 480         vector< complex< double > > valeurs_propres_dist(taille_vp,0);
 481         for( int i=0 ; i<taille_vp ; i++)
 482                 valeurs_propres_dist[i] = valeurs_propres[i];
 483
 484         if(_entropicCorrection || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
 485                 InvMatriceRoe( valeurs_propres_dist);
 486                 Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
 487         }
 488         else if (_spaceScheme==upwind)//upwind sans entropic
 489                 SigneMatriceRoe( valeurs_propres_dist);
 490         else if (_spaceScheme==centered)//centre  sans entropic
 491         {
 492                 for(int i=0; i<_nVar*_nVar;i++)
 493                         _signAroe[i] = 0;
 494         }
 495         else  if(_spaceScheme ==staggered )
 496         {
 497                 double signu1,signu2;
 498                 if(u1_n>0)
 499                         signu1=1;
 500                 else if (u1_n<0)
 501                         signu1=-1;
 502                 else
 503                         signu1=0;
 504                 if(u2_n>0)
 505                         signu2=1;
 506                 else if (u2_n<0)
 507                         signu2=-1;
 508                 else
 509                         signu2=0;
 510                 for(int i=0; i<_nVar*_nVar;i++)
 511                         _signAroe[i] = 0;
 512                 _signAroe[0] = signu1;
 513                 _signAroe[(1+_Ndim)*_nVar +1+_Ndim] = signu2;
 514                 for(int i=2; i<_nVar-1;i++){
 515                         _signAroe[i*_nVar+i] = -signu1;
 516                         _signAroe[(i+1+_Ndim)*_nVar+i+1+_Ndim] = -signu2;
 517                 }
 518         }
 519         else
 520                 throw CdmathException("IsothermalTwoFluid::convectionMatrices: well balanced option not treated");
 521
 522         for(int i=0; i<_nVar*_nVar;i++)
 523         {
 524                 _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
 525                 _AroePlus[i]  = (_Aroe[i]+_absAroe[i])/2;
 526         }
 527         if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//Implicitation using primitive variables
 528                 for(int i=0; i<_nVar*_nVar;i++)
 529                         _AroeMinusImplicit[i] = (_AroeImplicit[i]-_absAroeImplicit[i])/2;
 530         else
 531                 for(int i=0; i<_nVar*_nVar;i++)
 532                         _AroeMinusImplicit[i] = _AroeMinus[i];
 533
 534         if(_verbose && _nbTimeStep%_freqSave ==0)
 535         {
 536                 cout<<endl<<"Matrice de Roe"<<endl;
 537                 for(int i=0; i<_nVar;i++)
 538                 {
 539                         for(int j=0; j<_nVar;j++)
 540                                 cout << _Aroe[i*_nVar+j]<< " , ";
 541                         cout<<endl;
 542                 }
 543                 cout<<"Valeur absolue matrice de Roe"<<endl;
 544                 for(int i=0; i<_nVar;i++){
 545                         for(int j=0; j<_nVar;j++)
 546                                 cout<<_absAroe[i*_nVar+j]<<" , ";
 547                         cout<<endl;
 548                 }
 549                 cout<<endl<<"Matrice _AroeMinus"<<endl;
 550                 for(int i=0; i<_nVar;i++)
 551                 {
 552                         for(int j=0; j<_nVar;j++)
 553                                 cout << _AroeMinus[i*_nVar+j]<< " , ";
 554                         cout<<endl;
 555                 }
 556                 cout<<endl<<"Matrice _AroePlus"<<endl;
 557                 for(int i=0; i<_nVar;i++)
 558                 {
 559                         for(int j=0; j<_nVar;j++)
 560                                 cout << _AroePlus[i*_nVar+j]<< " , ";
 561                         cout<<endl;
 562                 }
 563         }
 564 }
 565
 566 void IsothermalTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *normale){
 567         //To do controle signe des vitesses pour CL entree/sortie
 568         int k;
 569         _idm[0] = _nVar*j;
 570         for(k=1; k<_nVar; k++)
 571                 _idm[k] = _idm[k-1] + 1;
 572
 573         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
 574         double q1_n=0, q2_n=0;//quantité de mouvement normale à la face limite
 575         for(k=0; k<_Ndim; k++){
 576                 q1_n+=_externalStates[(k+1)]*normale[k];
 577                 q2_n+=_externalStates[(k+1+1+_Ndim)]*normale[k];
 578         }
 579
 580         if(_verbose && _nbTimeStep%_freqSave ==0)
 581         {
 582                 cout << "Boundary conditions for group "<< nameOfGroup<< ", inner cell j= "<<j << " face unit normal vector "<<endl;
 583                 for(k=0; k<_Ndim; k++){
 584                         cout<<normale[k]<<", ";
 585                 }
 586                 cout<<endl;
 587         }
 588         if (_limitField[nameOfGroup].bcType==Wall){
 589
 590                 //Pour la convection, inversion du sens des vitesses
 591                 for(k=0; k<_Ndim; k++){
 592                         _externalStates[(k+1)]-= 2*q1_n*normale[k];
 593                         _externalStates[(k+1+1+_Ndim)]-= 2*q2_n*normale[k];
 594                 }
 595
 596                 _idm[0] = 0;
 597                 for(k=1; k<_nVar; k++)
 598                         _idm[k] = _idm[k-1] + 1;
 599
 600                 VecAssemblyBegin(_Uext);
 601                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 602                 VecAssemblyEnd(_Uext);
 603
 604                 //Pour la diffusion, paroi à vitesses imposees
 605                 _externalStates[1]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
 606                 _externalStates[2+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_x[1];
 607                 if(_Ndim>1)
 608                 {
 609                         _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
 610                         _externalStates[3+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_y[1];
 611                         if(_Ndim==3)
 612                         {
 613                                 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
 614                                 _externalStates[4+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_z[1];
 615                         }
 616                 }
 617                 VecAssemblyBegin(_Uextdiff);
 618                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 619                 VecAssemblyEnd(_Uextdiff);
 620         }
 621         else if (_limitField[nameOfGroup].bcType==Neumann){
 622                 _idm[0] = 0;
 623                 for(k=1; k<_nVar; k++)
 624                         _idm[k] = _idm[k-1] + 1;
 625
 626                 VecAssemblyBegin(_Uext);
 627                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 628                 VecAssemblyEnd(_Uext);
 629
 630                 VecAssemblyBegin(_Uextdiff);
 631                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 632                 VecAssemblyEnd(_Uextdiff);
 633         }
 634         else if (_limitField[nameOfGroup].bcType==Inlet){
 635                 _idm[0] = _nVar*j;
 636                 for(k=1; k<_nVar; k++)
 637                         _idm[k] = _idm[k-1] + 1;
 638
 639                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 640                 double alpha=_limitField[nameOfGroup].alpha;
 641                 double pression=_Vj[1];
 642                 double T=_Temperature;
 643                 double rho_v=_fluides[0]->getDensity(pression,T);
 644                 double rho_l=_fluides[1]->getDensity(pression,T);
 645                 _externalStates[0]=alpha*rho_v;
 646                 _externalStates[1]=alpha*rho_v*(_limitField[nameOfGroup].v_x[0]);
 647                 _externalStates[1+_Ndim]=(1-alpha)*rho_l;
 648                 _externalStates[2+_Ndim]=(1-alpha)*rho_l*(_limitField[nameOfGroup].v_x[1]);
 649                 if(_Ndim>1)
 650                 {
 651                         _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
 652                         _externalStates[3+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_y[1];
 653                         if(_Ndim==3)
 654                         {
 655                                 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
 656                                 _externalStates[4+_Ndim]=_externalStates[1+_Ndim]*_limitField[nameOfGroup].v_z[1];
 657                         }
 658                 }
 659                 _idm[0] = 0;
 660                 for(k=1; k<_nVar; k++)
 661                         _idm[k] = _idm[k-1] + 1;
 662                 VecAssemblyBegin(_Uext);
 663                 VecAssemblyBegin(_Uextdiff);
 664                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 665                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 666                 VecAssemblyEnd(_Uext);
 667                 VecAssemblyEnd(_Uextdiff);
 668         }
 669         else if (_limitField[nameOfGroup].bcType==InletPressure){
 670                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 671                 Cell Cj=_mesh.getCell(j);
 672                 double hydroPress=Cj.x()*_GravityField3d[0];
 673                 if(_Ndim>1){
 674                         hydroPress+=Cj.y()*_GravityField3d[1];
 675                         if(_Ndim>2)
 676                                 hydroPress+=Cj.z()*_GravityField3d[2];
 677                 }
 678                 hydroPress*=_externalStates[0]+_externalStates[_Ndim];//multiplication by rho the total density
 679
 680                 //Building the external state
 681                 _idm[0] = _nVar*j;
 682                 for(k=1; k<_nVar; k++)
 683                         _idm[k] = _idm[k-1] + 1;
 684
 685                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 686                 double pression=_limitField[nameOfGroup].p+hydroPress;
 687                 double alpha=_limitField[nameOfGroup].alpha;
 688                 double T=_Temperature;
 689                 double rho_v=_fluides[0]->getDensity(pression,T);
 690                 double rho_l=_fluides[1]->getDensity(pression,T);
 691                 _externalStates[0]=alpha*rho_v;
 692                 _externalStates[1+_Ndim]=(1-alpha)*rho_l;
 693
 694                 for(k=1;k<1+_Ndim;k++){
 695                         _externalStates[k]=_externalStates[0]*_Vj[1+k];
 696                         _externalStates[k+1+_Ndim]=_externalStates[1+_Ndim]*_Vj[k+1+_Ndim];
 697                 }
 698                 _idm[0] = 0;
 699                 for(k=1; k<_nVar; k++)
 700                         _idm[k] = _idm[k-1] + 1;
 701                 VecAssemblyBegin(_Uext);
 702                 VecAssemblyBegin(_Uextdiff);
 703                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 704                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 705                 VecAssemblyEnd(_Uext);
 706                 VecAssemblyEnd(_Uextdiff);
 707
 708                 if(_verbose && _nbTimeStep%_freqSave ==0)
 709                 {
 710                         cout<<"Etat fantôme InletPressure"<<endl;
 711                         for(int k=0;k<_nVar;k++)
 712                                 cout<<_externalStates[k]<<endl;
 713                 }
 714         }
 715         else if (_limitField[nameOfGroup].bcType==Outlet){
 716                 _idm[0] = _nVar*j;
 717                 for(k=1; k<_nVar; k++)
 718                         _idm[k] = _idm[k-1] + 1;
 719
 720                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 721                 Cell Cj=_mesh.getCell(j);
 722                 double hydroPress=Cj.x()*_GravityField3d[0];
 723                 if(_Ndim>1){
 724                         hydroPress+=Cj.y()*_GravityField3d[1];
 725                         if(_Ndim>2)
 726                                 hydroPress+=Cj.z()*_GravityField3d[2];
 727                 }
 728                 hydroPress*=_externalStates[0]+_externalStates[_Ndim];//multiplication by rho the total density
 729
 730                 //Building the external state
 731                 VecGetValues(_primitiveVars, _nVar, _idm,_Vj);
 732                 double pression_int=_Vj[1];
 733                 double pression_ext=_limitField[nameOfGroup].p+hydroPress;
 734                 double T=_Vj[_nVar-1];
 735                 double rho_v_int=_fluides[0]->getDensity(pression_int,T);
 736                 double rho_l_int=_fluides[1]->getDensity(pression_int,T);
 737                 double rho_v_ext=_fluides[0]->getDensity(pression_ext,T);
 738                 double rho_l_ext=_fluides[1]->getDensity(pression_ext,T);
 739
 740                 for(k=0;k<1+_Ndim;k++){
 741                         _externalStates[k]*=rho_v_ext/rho_v_int;
 742                         _externalStates[k+1+_Ndim]*=rho_l_ext/rho_l_int;
 743                 }
 744                 _idm[0] = 0;
 745                 for(k=1; k<_nVar; k++)
 746                         _idm[k] = _idm[k-1] + 1;
 747                 VecAssemblyBegin(_Uext);
 748                 VecAssemblyBegin(_Uextdiff);
 749                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 750                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 751                 VecAssemblyEnd(_Uext);
 752                 VecAssemblyEnd(_Uextdiff);
 753         }
 754         else{
 755                 cout<<"!!!!!!!!!!!!!!!!! Error IsothermalTwoFluid::setBoundaryState !!!!!!!!!!"<<endl;
 756                 cout<<"!!!!!!!!!!!!! Boundary condition not set for boundary named"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 757                 cout<<"Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
 758                 throw CdmathException("Unknown boundary condition");
 759         }
 760 }
 761
 762 void IsothermalTwoFluid::addDiffusionToSecondMember
 763 (               const int &i,
 764                 const int &j,
 765                 bool isBord)
 766
 767 {
 768         double mu1 = _fluides[0]->getViscosity(_Temperature);
 769         double mu2 = _fluides[1]->getViscosity(_Temperature);
 770
 771         if(mu1==0 && mu2==0)
 772                 return;
 773
 774         //extraction des valeurs
 775         _idm[0] = _nVar*i;
 776         for(int k=1; k<_nVar; k++)
 777                 _idm[k] = _idm[k-1] + 1;
 778
 779         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
 780         if (_verbose && _nbTimeStep%_freqSave ==0)
 781         {
 782                 cout << "Contribution diffusion: variables primitives maille " << i<<endl;
 783                 for(int q=0; q<_nVar; q++)
 784                 {
 785                         cout << _Vi[q] << endl;
 786                 }
 787                 cout << endl;
 788         }
 789
 790         if(!isBord ){
 791                 for(int k=0; k<_nVar; k++)
 792                         _idm[k] = _nVar*j + k;
 793                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 794         }
 795         else
 796         {
 797                 for(int k=0; k<_nVar; k++)
 798                         _idm[k] = k;
 799
 800                 VecGetValues(_Uextdiff, _nVar, _idm, _phi);
 801                 consToPrim(_phi,_Vj);
 802         }
 803
 804         if (_verbose && _nbTimeStep%_freqSave ==0)
 805         {
 806                 cout << "Contribution diffusion: variables primitives maille " <<j <<endl;
 807                 for(int q=0; q<_nVar; q++)
 808                 {
 809                         cout << _Vj[q] << endl;
 810                 }
 811                 cout << endl;
 812         }
 813
 814
 815         double alpha=(_Vj[0]+_Vi[0])/2;
 816         //on n'a pas de contribution sur la masse
 817         _phi[0]=0;
 818         _phi[_Ndim+1]=0;
 819         //contribution visqueuse sur la quantite de mouvement
 820         for(int k=1; k<_Ndim+1; k++)
 821         {
 822                 _phi[k]         = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*alpha*mu1*(_Vj[k] - _Vi[k]);//attention car primitif=alpha p u1 u2
 823                 _phi[k+_Ndim+1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*(1-alpha)*mu2*(_Vj[1+k+_Ndim] - _Vi[1+k+_Ndim]);
 824         }
 825
 826         _idm[0] = i;
 827         VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
 828
 829         if(_verbose && _nbTimeStep%_freqSave ==0)
 830         {
 831                 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
 832                 for(int q=0; q<_nVar; q++)
 833                 {
 834                         cout << _phi[q] << endl;
 835                 }
 836                 cout << endl;
 837         }
 838
 839         if(!isBord)
 840         {
 841                 //On change de signe pour l'autre contribution
 842                 for(int k=0; k<_nVar; k++)
 843                         _phi[k] *= -_inv_dxj/_inv_dxi;
 844                 _idm[0] = j;
 845
 846                 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
 847         }
 848
 849         if(_verbose && _nbTimeStep%_freqSave ==0)
 850         {
 851                 cout << "Contribution diffusion au 2nd membre pour la maille  " << j << ": "<<endl;
 852                 for(int q=0; q<_nVar; q++)
 853                 {
 854                         cout << _phi[q] << endl;
 855                 }
 856                 cout << endl;
 857
 858                 if(_timeScheme==Implicit)
 859                 {
 860                         cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
 861                         for(int i=0; i<_nVar; i++)
 862                         {
 863                                 for(int j=0; j<_nVar; j++)
 864                                         cout << _Diffusion[i*_nVar+j]<<", ";
 865                                 cout << endl;
 866                         }
 867                         cout << endl;
 868                 }
 869         }
 870 }
 871
 872 void IsothermalTwoFluid::sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i)
 873 {
 874         double m1=Ui[0],m2=Ui[1+_Ndim],alpha=Vi[0], P=Vi[1];
 875         double norm_ur=0, Gamma;
 876         for(int k=0; k<_Ndim; k++)
 877                 norm_ur+=(Vi[2+k]-Vi[2+k+_Ndim])*(Vi[2+k]-Vi[2+k+_Ndim]);
 878         norm_ur=sqrt(norm_ur);
 879
 880         if(i>=0 &&_Temperature>_Tsat && alpha<1-_precision)
 881                 Gamma=_heatPowerField(i)/_latentHeat;
 882         else
 883                 Gamma=0;
 884         for(int k=1; k<_Ndim+1; k++)
 885         {
 886                 //cout<<"Vi[1+"<<k+_Ndim<<"]="<<Vi[1+k+_Ndim]<<endl;
 887                 Si[k] =_gravite[k]*m1 -_dragCoeffs[0]*norm_ur*(Vi[1+k]-Vi[1+k+_Ndim])+ Gamma*Vi[1+k+_Ndim];//interfacial velocity= ul
 888                 Si[k+_Ndim+1] =_gravite[k+_Ndim+1]*m2 + _dragCoeffs[0]*norm_ur*(Vi[1+k]-Vi[1+k+_Ndim])- Gamma*Vi[1+k+_Ndim];
 889         }
 890         if(true){//heated boiling
 891                 Si[0]=Gamma;
 892                 Si[1+_Ndim]=-Gamma;
 893         }
 894         else if (P<_Psat && alpha<1-_precision){//flash boiling
 895                 Si[0]=-_dHsatl_over_dp*_dp_over_dt(i)/_latentHeat;
 896                 Si[1+_Ndim]=_dHsatl_over_dp*_dp_over_dt(i)/_latentHeat;
 897         }
 898         else{
 899                 Si[0]=0;
 900                 Si[1+_Ndim]=0;
 901         }
 902
 903         if(_timeScheme==Implicit)
 904         {
 905                 for(int i=0; i<_nVar*_nVar;i++)
 906                         _GravityImplicitationMatrix[i] = 0;
 907                 for(int i=0; i<_nVar/2;i++)
 908                         _GravityImplicitationMatrix[i*_nVar]=-_gravite[i];
 909                 for(int i=_nVar/2; i<_nVar;i++)
 910                         _GravityImplicitationMatrix[i*_nVar+_nVar/2]=-_gravite[i];
 911         }
 912
 913         if(_verbose && _nbTimeStep%_freqSave ==0)
 914         {
 915                 cout<<"IsothermalTwoFluid::sourceVector"<<endl;
 916                 cout<<"Ui="<<endl;
 917                 for(int k=0;k<_nVar;k++)
 918                         cout<<Ui[k]<<", ";
 919                 cout<<endl;
 920                 cout<<"Vi="<<endl;
 921                 for(int k=0;k<_nVar;k++)
 922                         cout<<Vi[k]<<", ";
 923                 cout<<endl;
 924                 cout<<"Si="<<endl;
 925                 for(int k=0;k<_nVar;k++)
 926                         cout<<Si[k]<<", ";
 927                 cout<<endl;
 928                 if(_timeScheme==Implicit)
 929                         displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
 930         }
 931 }
 932 void IsothermalTwoFluid::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
 933 {
 934         double norm_u1=0, u1_n=0, norm_u2=0, u2_n=0, m1, m2;
 935         for(int i=0;i<_Ndim;i++){
 936                 u1_n += _Uroe[2+i]      *_vec_normal[i];
 937                 u2_n += _Uroe[2+i+_Ndim]*_vec_normal[i];
 938         }
 939         pressureLoss[0]=0;
 940         pressureLoss[1+_Ndim]=0;
 941         if(u1_n>0){
 942                 for(int i=0;i<_Ndim;i++)
 943                         norm_u1 += Vi[1+i]*Vi[1+i];
 944                 norm_u1=sqrt(norm_u1);
 945                 m1=Ui[0];
 946                 for(int i=0;i<_Ndim;i++)
 947                         pressureLoss[1+i]=-K*m1*norm_u1*Vi[1+i];
 948         }
 949         else{
 950                 for(int i=0;i<_Ndim;i++)
 951                         norm_u1 += Vj[1+i]*Vj[1+i];
 952                 norm_u1=sqrt(norm_u1);
 953                 m1=Uj[0];
 954                 for(int i=0;i<_Ndim;i++)
 955                         pressureLoss[1+i]=-K*m1*norm_u1*Vj[1+i];
 956         }
 957         if(u2_n>0){
 958                 for(int i=0;i<_Ndim;i++)
 959                         norm_u2 += Vi[2+i+_Ndim]*Vi[2+i+_Ndim];
 960                 norm_u2=sqrt(norm_u2);
 961                 m2=Ui[1+_Ndim];
 962                 for(int i=0;i<_Ndim;i++)
 963                         pressureLoss[2+i+_Ndim]=-K*m2*norm_u2*Vi[2+i+_Ndim];
 964         }
 965         else{
 966                 for(int i=0;i<_Ndim;i++)
 967                         norm_u2 += Vj[2+i+_Ndim]*Vj[2+i+_Ndim];
 968                 norm_u2=sqrt(norm_u2);
 969                 m2=Uj[1+_Ndim];
 970                 for(int i=0;i<_Ndim;i++)
 971                         pressureLoss[2+i+_Ndim]=-K*m2*norm_u2*Vj[2+i+_Ndim];
 972         }
 973         if(_verbose && _nbTimeStep%_freqSave ==0)
 974         {
 975                 cout<<"IsothermalTwoFluid::pressureLossVector K= "<<K<<endl;
 976                 cout<<"Ui="<<endl;
 977                 for(int k=0;k<_nVar;k++)
 978                         cout<<Ui[k]<<", ";
 979                 cout<<endl;
 980                 cout<<"Vi="<<endl;
 981                 for(int k=0;k<_nVar;k++)
 982                         cout<<Vi[k]<<", ";
 983                 cout<<endl;
 984                 cout<<"Uj="<<endl;
 985                 for(int k=0;k<_nVar;k++)
 986                         cout<<Uj[k]<<", ";
 987                 cout<<endl;
 988                 cout<<"Vj="<<endl;
 989                 for(int k=0;k<_nVar;k++)
 990                         cout<<Vj[k]<<", ";
 991                 cout<<endl;
 992                 cout<<"pressureLoss="<<endl;
 993                 for(int k=0;k<_nVar;k++)
 994                         cout<<pressureLoss[k]<<", ";
 995                 cout<<endl;
 996         }
 997 }
 998
 999 void IsothermalTwoFluid::porosityGradientSourceVector()
1000 {
1001         double u1_ni=0, u1_nj=0, u2_ni=0, u2_nj=0, rho1i, rho2i, rho1j, rho2j, pi=_Vi[1], pj=_Vj[1], pij1, pij2, alphaij=_Uroe[0];
1002         for(int i=0;i<_Ndim;i++) {
1003                 u1_ni += _Vi[2+i]*_vec_normal[i];
1004                 u2_ni += _Vi[2+_Ndim+i]*_vec_normal[i];
1005                 u1_nj += _Vj[2+i]*_vec_normal[i];
1006                 u2_nj += _Vj[2+_Ndim+i]*_vec_normal[i];
1007         }
1008         _porosityGradientSourceVector[0]=0;
1009         _porosityGradientSourceVector[1+_Ndim]=0;
1010         rho1i = _fluides[0]->getDensity(pi, _Temperature);
1011         rho2i = _fluides[1]->getDensity(pi, _Temperature);
1012         rho1j = _fluides[0]->getDensity(pj, _Temperature);
1013         rho2j = _fluides[1]->getDensity(pj, _Temperature);
1014         pij1=(pi+pj)/2+rho1i*rho1j/2/(rho1i+rho1j)*(u1_ni-u1_nj)*(u1_ni-u1_nj);
1015         pij2=(pi+pj)/2+rho2i*rho2j/2/(rho2i+rho2j)*(u2_ni-u2_nj)*(u2_ni-u2_nj);
1016         for(int i=0;i<_Ndim;i++){
1017                 _porosityGradientSourceVector[1+i]      =alphaij*pij1*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1018                 _porosityGradientSourceVector[2+_Ndim+i]=alphaij*pij2*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1019         }
1020 }
1021
1022 /* Funtion of equations of states */
1023 double IsothermalTwoFluid::ecartPression(double m1,double m2, double alpha, double e1, double e2){
1024         if(alpha>_precision*_precision&& alpha<1-_precision*_precision)
1025                 return _fluides[0]->getPressure((m1/alpha)*e1,m1/alpha) - _fluides[1]->getPressure((m2/(1-alpha))*e2,m2/(1-alpha));
1026         else if(alpha<=_precision*_precision)
1027         {
1028                 //cout<<"Warning ecartPression, alpha close to 0"<<endl;
1029                 //       if(m1<_precision*_precision*_precision)
1030                 //      return 0;
1031                 //       else
1032                 return alpha/_precision*_precision;
1033         }
1034         else
1035         {
1036                 //cout<<"Warning ecartPression, alpha close to 1"<<endl;
1037                 //       if(m2<_precision*_precision*_precision)
1038                 //      return 0;
1039                 //       else
1040                 return -(1-alpha)/_precision*_precision;
1041         }
1042 }
1043
1044 double IsothermalTwoFluid::ecartPressionDerivee(double m1,double m2, double alpha, double e1, double e2){
1045         if(alpha>_precision*_precision && alpha<1-_precision*_precision )
1046                 return -(m1/alpha)/alpha*e1*_fluides[0]->getPressureDerivativeRhoE() - (m2/(1-alpha))/(1-alpha)*e2*_fluides[1]->getPressureDerivativeRhoE();
1047         else if (alpha<_precision*_precision)
1048                 return -1/_precision*_precision;
1049         else
1050                 return 1/_precision*_precision;
1051 }
1052
1053 double IsothermalTwoFluid::intPressDef(double alpha, double u_r2, double rho1, double rho2)
1054 {
1055         return  _intPressCoeff*alpha*(1-alpha)*rho1*rho2*u_r2/( alpha*rho2+(1-alpha)*rho1)
1056                         +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
1057                         *(alpha*alpha*rho2-(1-alpha)*(1-alpha)*rho1)
1058                         *(alpha*alpha*rho2*rho2/(_fluides[0]->vitesseSonTemperature(_Temperature,rho1)*_fluides[0]->vitesseSonTemperature(_Temperature,rho1))
1059                                         -(1-alpha)*(1-alpha)*rho1*rho1/(_fluides[1]->vitesseSonTemperature(_Temperature,rho2)*_fluides[1]->vitesseSonTemperature(_Temperature,rho2)));
1060 }
1061
1062 void IsothermalTwoFluid::entropicShift(double* n)
1063 {
1064         vector< double > pol_car;
1065
1066         /*Left values */
1067         double u1_n = 0, u1_2 = 0, u2_n = 0, u2_2 = 0, u_r2 = 0;
1068         for(int i=0;i<_Ndim;i++)
1069         {
1070                 u1_2 += _l[2+i]*_l[2+i];
1071                 u1_n += _l[2+i]*n[i];
1072                 u2_2 += _l[1+i+1+_Ndim]*_l[1+i+1+_Ndim];
1073                 u2_n += _l[1+i+1+_Ndim]*n[i];
1074                 u_r2 += (_l[2+i]-_l[1+i+1+_Ndim])*(_l[2+i]-_l[1+i+1+_Ndim]);
1075         }
1076         double alpha = _l[0];
1077         double p = _l[1];
1078         double rho1 = _fluides[0]->getDensity(p, _Temperature);
1079         double rho2 = _fluides[1]->getDensity(p, _Temperature);
1080         double dpi1 = intPressDef(alpha, u_r2, rho1, rho2);
1081         double dpi2 = dpi1;
1082         double invcsq1 = 1/_fluides[0]->vitesseSonTemperature(_Temperature,rho1);
1083         invcsq1*=invcsq1;
1084         double invcsq2 = 1/_fluides[1]->vitesseSonTemperature(_Temperature,rho2);
1085         invcsq2*=invcsq2;
1086
1087         pol_car= Polynoms::polynome_caracteristique(alpha,  1-alpha, u1_n, u2_n, rho1, rho2, invcsq1, invcsq2, dpi1, dpi2);
1088         for(int ct=0;ct<4;ct++){
1089                 if (abs(pol_car[ct])<_precision*_precision)
1090                         pol_car[ct]=0;
1091         }
1092         vector< complex<double> > vp_left = getRacines(pol_car);
1093         int taille_vp_left = Polynoms::new_tri_selectif(vp_left,vp_left.size(),_precision);
1094
1095         if(_verbose && _nbTimeStep%_freqSave ==0)
1096         {
1097                 cout<<"Entropic shift left eigenvalues: "<<endl;
1098                 for(unsigned int ct =0; ct<vp_left.size(); ct++)
1099                         cout<<"("<< vp_left[ct].real()<< ", " <<vp_left[ct].imag() << ")";
1100                 cout<<endl;
1101         }
1102
1103         /*right values */
1104         u1_2 = 0; u2_2=0;
1105         u1_n = 0; u2_n=0;
1106         u_r2=0;
1107         for(int i=0;i<_Ndim;i++)
1108         {
1109                 u1_2 += _r[2+i]*_r[2+i];
1110                 u1_n += _r[2+i]*n[i];
1111                 u2_2 += _r[1+i+1+_Ndim]*_r[1+i+1+_Ndim];
1112                 u2_n += _r[1+i+1+_Ndim]*n[i];
1113                 u_r2 += (_r[2+i]-_r[1+i+1+_Ndim])*(_r[2+i]-_r[1+i+1+_Ndim]);
1114         }
1115         alpha = _r[0];
1116         p = _r[1];
1117         rho1 = _fluides[0]->getDensity(p, _Temperature);
1118         rho2 = _fluides[1]->getDensity(p, _Temperature);
1119         dpi1 = intPressDef(alpha, u_r2, rho1, rho2);
1120         dpi2 = dpi1;
1121         invcsq1 = 1/_fluides[0]->vitesseSonTemperature(_Temperature,rho1);
1122         invcsq1*=invcsq1;
1123         invcsq2 = 1/_fluides[1]->vitesseSonTemperature(_Temperature,rho2);
1124         invcsq2*=invcsq2;
1125
1126         pol_car= Polynoms::polynome_caracteristique(alpha,  1-alpha, u1_n, u2_n, rho1, rho2, invcsq1, invcsq2, dpi1, dpi2);
1127         for(int ct=0;ct<4;ct++){
1128                 if (abs(pol_car[ct])<_precision*_precision)
1129                         pol_car[ct]=0;
1130         }
1131         vector< complex<double> > vp_right = getRacines(pol_car);
1132         int taille_vp_right = Polynoms::new_tri_selectif(vp_right,vp_right.size(),_precision);
1133
1134         if(_verbose && _nbTimeStep%_freqSave ==0)
1135         {
1136                 cout<<"Entropic shift right eigenvalues: "<<endl;
1137                 for(unsigned int ct =0; ct<vp_right.size(); ct++)
1138                         cout<<"("<<vp_right[ct].real()<< ", " <<vp_right[ct].imag() <<")";
1139                 cout<< endl;
1140         }
1141         _entropicShift[0] = abs(vp_left[0]-vp_right[0]);
1142         _entropicShift[2] = abs(vp_left[taille_vp_left-1]-vp_right[taille_vp_right-1]);
1143         _entropicShift[1]=0;
1144         for(int i=1;i<min(taille_vp_right-1,taille_vp_left-1);i++)
1145                 _entropicShift[1] = max(_entropicShift[1],abs(vp_left[i]-vp_right[i]));
1146         if(_verbose && _nbTimeStep%_freqSave ==0)
1147         {
1148                 cout<<"eigenvalue jumps "<<endl;
1149                 cout<< _entropicShift[0] << ", " << _entropicShift[1] << ", "<< _entropicShift[2] <<endl;
1150         }
1151 }
1152
1153 void IsothermalTwoFluid::computeScaling(double maxvp)
1154 {
1155         //      double alphaScaling;
1156         //      if(_guessalpha>_precision && _guessalpha<1-_precision)
1157         //              alphaScaling=_guessalpha;
1158         //      else
1159         //              alphaScaling=0.5;
1160
1161         _blockDiag[0]=1;//alphaScaling;
1162         _invBlockDiag[0]=1;//_blockDiag[0];
1163         _blockDiag[1+_Ndim]=1;//-alphaScaling;
1164         _invBlockDiag[1+_Ndim]=1.0;//_blockDiag[1+_Ndim];
1165         for(int q=1; q<_Ndim+1; q++)
1166         {
1167                 _blockDiag[q]=1/(maxvp*maxvp);
1168                 _invBlockDiag[q]=1;//_blockDiag[q];
1169                 _blockDiag[q+1+_Ndim]=1/(maxvp*maxvp);
1170                 _invBlockDiag[q+1+_Ndim]=1;//_blockDiag[q+1+_Ndim];
1171         }
1172 }
1173
1174 void IsothermalTwoFluid::jacobian(const int &j, string nameOfGroup,double * normale)
1175 {
1176         int k;
1177         for(k=0; k<_nVar*_nVar;k++)
1178                 _Jcb[k] = 0;//No implicitation at this stage
1179
1180         // loop on boundaries
1181         if (_limitField[nameOfGroup].bcType==Wall)
1182         {
1183                 for(k=0; k<_nVar;k++)
1184                         _Jcb[k*_nVar + k] = 1;
1185                 for(k=1; k<1+_Ndim;k++)
1186                         for(int l=1; l<1+_Ndim;l++){
1187                                 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
1188                                 _Jcb[(k+1+_Ndim)*_nVar + l+1+_Ndim] -= 2*normale[k-1]*normale[l-1];
1189                         }
1190         }
1191         //not wall
1192         else if (_limitField[nameOfGroup].bcType==Inlet)
1193         {
1194                 _Jcb[0] = 1;
1195                 _Jcb[(1+_Ndim)*_nVar +1+_Ndim] = 1;
1196                 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];
1197                 _Jcb[(2+_Ndim)*_nVar +1+_Ndim] =_limitField[nameOfGroup].v_x[1];
1198                 if(_Ndim>1)
1199                 {
1200                         _Jcb[2*_nVar]=_limitField[nameOfGroup].v_y[0];
1201                         _Jcb[(3+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_y[1];
1202                         if(_Ndim==3)
1203                         {
1204                                 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1205                                 _Jcb[(4+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_z[1];
1206
1207                         }
1208                 }
1209         }
1210         // not wall, not inlet
1211         else if (_limitField[nameOfGroup].bcType==Outlet){
1212                 _idm[0] = j*_nVar;
1213                 for(k=1; k<_nVar;k++)
1214                 {_idm[k] = _idm[k-1] + 1;}
1215                 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1216                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1217         }
1218         else  if (_limitField[nameOfGroup].bcType!=Neumann && _limitField[nameOfGroup].bcType!=InletPressure)// not wall, not inlet, not outlet
1219         {
1220                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1221                 throw CdmathException("IsothermalTwoFluid::jacobianDiff: This boundary condition is not treated");
1222         }
1223 }
1224
1225 void IsothermalTwoFluid::jacobianDiff(const int &j, string nameOfGroup)
1226 {
1227
1228         int k;
1229         for(k=0; k<_nVar*_nVar;k++)
1230                 _JcbDiff[k] = 0;
1231         if (_limitField[nameOfGroup].bcType==Wall){
1232                 _JcbDiff[0] = 1;
1233                 _JcbDiff[(1+_Ndim)*_nVar +1+_Ndim] = 1;
1234                 _JcbDiff[_nVar]=_limitField[nameOfGroup].v_x[0];
1235                 _JcbDiff[(2+_Ndim)*_nVar +1+_Ndim] =_limitField[nameOfGroup].v_x[1];
1236                 if(_Ndim>1)
1237                 {
1238                         _JcbDiff[2*_nVar]=_limitField[nameOfGroup].v_y[0];
1239                         _JcbDiff[(3+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_y[1];
1240
1241                         if(_Ndim==3)
1242                         {
1243                                 _JcbDiff[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1244                                 _JcbDiff[(4+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_z[1];
1245                         }
1246                 }
1247         } else if (_limitField[nameOfGroup].bcType==Inlet){
1248                 /*
1249                 _JcbDiff[0] = 1;
1250                 _JcbDiff[(1+_Ndim)*_nVar +1+_Ndim] = 1;
1251                 _JcbDiff[_nVar]=_limitField[nameOfGroup].v_x[0];
1252                 _JcbDiff[(2+_Ndim)*_nVar +1+_Ndim] =_limitField[nameOfGroup].v_x[1];
1253                 if(_Ndim>1)
1254                 {
1255                         _JcbDiff[2*_nVar]=_limitField[nameOfGroup].v_y[0];
1256                         _JcbDiff[(3+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_y[1];
1257
1258                         if(_Ndim==3)
1259                         {
1260                                 _JcbDiff[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1261                                 _JcbDiff[(4+_Ndim)*_nVar +1+_Ndim]= _limitField[nameOfGroup].v_z[1];
1262                         }
1263                 }
1264                  */
1265         } else if (_limitField[nameOfGroup].bcType==Outlet){
1266                 /*
1267                 //extraction de l etat courant et primitives
1268                 _idm[0] = j*_nVar;
1269                 for(k=1; k<_nVar;k++)
1270                 {_idm[k] = _idm[k-1] + 1;}
1271                 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1272                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1273                  */
1274         }
1275         else if (_limitField[nameOfGroup].bcType!=Neumann && _limitField[nameOfGroup].bcType!=InletPressure){
1276                 cout<<"Condition  limite non traitee pour le bord "<<nameOfGroup<< endl;
1277                 throw CdmathException("IsothermalTwoFluid::jacobianDiff: Condition  limite non traitee");
1278         }
1279 }
1280
1281 void IsothermalTwoFluid::primToCons(const double *P, const int &i, double *W, const int &j){
1282         //P=alpha,p,u1,u2
1283         //W=m1,q1,m2,q2
1284         W[j*_nVar] = P[i*_nVar]*_fluides[0]->getDensity(P[i*_nVar+1],_Temperature);
1285         W[j*_nVar+1+_Ndim] = (1-P[i*_nVar])*_fluides[1]->getDensity(P[i*_nVar+1],_Temperature);
1286         for(int k=0; k<_Ndim; k++)
1287         {
1288                 W[j*_nVar+(k+1)] = W[j*_nVar]*P[i*_nVar+(k+2)];
1289                 W[j*_nVar+(k+1)+1+_Ndim] = W[j*_nVar+1+_Ndim]*P[i*_nVar+(k+2)+_Ndim];
1290         }
1291
1292 }
1293
1294 void IsothermalTwoFluid::consToPrim(const double *Wcons, double* Wprim,double porosity)//To do: treat porosity
1295 {
1296         //P=alpha,p,u1,u2
1297         //W=m1,q1,m2,q2
1298         double m1=Wcons[0];
1299         double m2=Wcons[_Ndim+1];
1300         double e1 = _internalEnergy1;
1301         double e2 = _internalEnergy2;
1302
1303         _minm1=min(m1,_minm1);
1304         _minm2=min(m2,_minm2);
1305         if(m1<-_precision || m2<-_precision)
1306                 _nbMaillesNeg+=1;
1307         if(fabs(m1)<_precision*_precision)
1308         {
1309                 Wprim[0]=0;
1310                 Wprim[1]=_fluides[1]->getPressure(m2*e2,m2);
1311
1312                 for(int idim=0; idim<_Ndim; idim++)
1313                 {
1314                         Wprim[2+_Ndim+idim] = Wcons[2+_Ndim+idim]/Wcons[1+_Ndim];
1315                         Wprim[2+idim] = Wprim[2+_Ndim+idim];
1316                 }
1317         }
1318         else if(fabs(m2)<_precision*_precision)
1319         {
1320                 Wprim[0]=1;
1321                 Wprim[1]=_fluides[0]->getPressure(m1*e1,m1);
1322                 for(int idim=0; idim<_Ndim; idim++)
1323                 {
1324                         Wprim[2+idim] = Wcons[1+idim]/Wcons[0];
1325                         Wprim[2+_Ndim+idim] = Wprim[2+idim];
1326                 }
1327         }
1328         else
1329         {
1330                 for(int idim=0; idim<_Ndim; idim++)
1331                 {
1332                         Wprim[2+idim] = Wcons[1+idim]/Wcons[0];
1333                         Wprim[2+_Ndim+idim] = Wcons[2+_Ndim+idim]/Wcons[1+_Ndim];
1334                 }
1335                 double alphanewton, alphainf=0, alphasup=1;
1336                 int iterMax=50, iter=0;
1337
1338                 double   dp=ecartPression( m1, m2, _guessalpha, e1, e2);
1339                 double   dpprim=ecartPressionDerivee( m1, m2, _guessalpha, e1, e2);
1340
1341                 while(fabs(dp)>1e3*_precision && iter<iterMax)
1342                 {
1343
1344                         //  if (dp>0)
1345                         //              alphainf = _guessalpha;
1346                         //            else
1347                         //              alphansup = _guessalpha;
1348
1349                         //            _guessalpha = (alphainf+alphansup)/2;
1350
1351                         if (dp>0)
1352                                 alphainf = _guessalpha;
1353                         else
1354                                 alphasup = _guessalpha;
1355
1356                         alphanewton = _guessalpha-dp/dpprim;
1357
1358                         if(alphanewton<=alphainf)
1359                         {
1360                                 _guessalpha = (9*alphainf+alphasup)/10;
1361                                 //cout<< "dichotomie"<<endl;
1362                         }
1363                         else if(alphanewton>=alphasup)
1364                         {
1365                                 _guessalpha = (alphainf+9*alphasup)/10;
1366                                 //cout<< "dichotomie"<<endl;
1367                         }
1368                         else
1369                                 _guessalpha=alphanewton;
1370
1371
1372                         if(_verbose && _nbTimeStep%_freqSave ==0)
1373                                 cout<<"consToPrim diphasique iter= " <<iter<<" dp= " << dp<<" dpprim= " << dpprim<< " _guessalpha= " << _guessalpha<<endl;
1374                         dp=ecartPression( m1, m2, _guessalpha, e1, e2);
1375                         dpprim=ecartPressionDerivee( m1, m2, _guessalpha, e1, e2);
1376
1377                         iter++;
1378                 }
1379                 if(_verbose && _nbTimeStep%_freqSave ==0)
1380                         cout<<endl;
1381
1382                 if(iter>=iterMax)
1383                 {
1384                         _err_press_max = max(_err_press_max,fabs(dp/1.e5));
1385                         //            if(_guessalpha>0.5)
1386                         //              _guessalpha=1;
1387                         //            else
1388                         //              _guessalpha=0;
1389                 }
1390                 if(_guessalpha>0.5)
1391                         Wprim[1]=_fluides[0]->getPressure((m1/_guessalpha)*e1,m1/_guessalpha);
1392                 else
1393                         Wprim[1]=_fluides[1]->getPressure((m2/(1-_guessalpha))*e2,m2/(1-_guessalpha));
1394                 Wprim[0]=_guessalpha;
1395         }
1396         if (Wprim[1]<0){
1397                 cout << "pressure = "<< Wprim[1] << " < 0 " << endl;
1398                 cout << "Conservative state = ";
1399                 for(int k=0; k<_nVar; k++){
1400                         cout<<Wcons[k]<<", ";
1401                 }
1402                 cout<<endl;
1403                 *_runLogFile<< "IsothermalTwoFluid::consToPrim: negative pressure = "<< Wprim[1] << " < 0 " << endl;
1404                 throw CdmathException("IsothermalTwoFluid::consToPrim: negative pressure");
1405         }
1406 }
1407
1408 Vector IsothermalTwoFluid::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
1409         if(_verbose){
1410                 cout<<"Ucons"<<endl;
1411                 cout<<U<<endl;
1412                 cout<<"Vprim"<<endl;
1413                 cout<<V<<endl;
1414         }
1415
1416         double phim1=U(0);//phi alpha1 rho1
1417         double phim2=U(1+_Ndim);//phi alpha2 rho2
1418         Vector phiq1(_Ndim),phiq2(_Ndim);//phi alpha1 rho1 u1, phi alpha2 rho2 u2
1419         for(int i=0;i<_Ndim;i++){
1420                 phiq1(i)=U(1+i);
1421                 phiq2(i)=U(2+_Ndim+i);
1422         }
1423         double alpha=V(0);
1424         double pression=V(1);
1425         Vector vitesse1(_Ndim),vitesse2(_Ndim);
1426         for(int i=0;i<_Ndim;i++){
1427                 vitesse1(i)=V(2+i);
1428                 vitesse2(i)=V(2+_Ndim+i);
1429         }
1430
1431         double vitesse1n=vitesse1*normale;
1432         double vitesse2n=vitesse2*normale;
1433
1434         double alpha_roe = _Uroe[0];//Toumi formula
1435         // interfacial pressure term (hyperbolic correction)
1436         double dpi = _Uroe[_nVar];
1437
1438         Vector F(_nVar);
1439         F(0)=phim1*vitesse1n;
1440         F(1+_Ndim)=phim2*vitesse2n;
1441         for(int i=0;i<_Ndim;i++){
1442                 F(1+i)=phim1*vitesse1n*vitesse1(i)+(alpha_roe*pression+dpi*alpha)*porosity*normale(i);
1443                 F(2+_Ndim+i)=phim2*vitesse2n*vitesse2(i)+((1-alpha_roe)*pression-dpi*alpha)*normale(i)*porosity;
1444         }
1445
1446         if(_verbose){
1447                 cout<<"Flux F(U,V)"<<endl;
1448                 cout<<F<<endl;
1449         }
1450
1451         return F;
1452 }
1453
1454 Vector IsothermalTwoFluid::staggeredVFFCFlux()
1455 {
1456         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1457                 throw CdmathException("IsothermalTwoFluid::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
1458         else//_spaceScheme==staggered
1459         {
1460                 Vector Fij(_nVar);
1461                 double alpha_roe = _Uroe[0];//Toumi formula
1462                 // interfacial pressure term (hyperbolic correction)
1463                 double dpi = _Uroe[_nVar];
1464
1465                 double u1ijn=0, u2ijn=0, phialphaq1n=0, phialphaq2n=0;
1466                 for(int idim=0;idim<_Ndim;idim++){//URoe = alpha, p, u1, u2, dpi
1467                         u1ijn+=_vec_normal[idim]*_Uroe[2+idim];
1468                         u2ijn+=_vec_normal[idim]*_Uroe[2+_Ndim+idim];
1469                 }
1470                 if(u1ijn>=0)
1471                 {
1472                         for(int idim=0;idim<_Ndim;idim++)
1473                                 phialphaq1n+=_vec_normal[idim]*_Ui[1+idim];//phi alpha rho u n
1474                         Fij(0)=phialphaq1n;
1475                         for(int idim=0;idim<_Ndim;idim++)
1476                                 Fij(1+idim)=phialphaq1n*_Vi[2+idim]+(alpha_roe*_Vj[1]*_porosityj+dpi*_Vi[0]*_porosityi)*_vec_normal[idim];
1477                 }
1478                 else
1479                 {
1480                         for(int idim=0;idim<_Ndim;idim++)
1481                                 phialphaq2n+=_vec_normal[idim]*_Uj[1+idim];//phi alpha rho u n
1482                         Fij(0)=phialphaq2n;
1483                         for(int idim=0;idim<_Ndim;idim++)
1484                                 Fij(1+idim)=phialphaq2n*_Vj[2+idim]+(alpha_roe*_Vi[1]*_porosityi+dpi*_Vj[0]*_porosityj)*_vec_normal[idim];
1485                 }
1486
1487                 if(u2ijn>=0)
1488                 {
1489                         for(int idim=0;idim<_Ndim;idim++)
1490                                 phialphaq2n+=_vec_normal[idim]*_Ui[2+_Ndim+idim];//phi alpha rho u n
1491                         Fij(1+_Ndim)=phialphaq2n;
1492                         for(int idim=0;idim<_Ndim;idim++)
1493                                 Fij(2+_Ndim+idim)=phialphaq2n*_Vi[2+_Ndim+idim]+((1-alpha_roe)*_Vj[1]*_porosityj-dpi*_Vi[0]*_porosityi)*_vec_normal[idim];
1494                 }
1495                 else
1496                 {
1497                         for(int idim=0;idim<_Ndim;idim++)
1498                                 phialphaq2n+=_vec_normal[idim]*_Uj[2+_Ndim+idim];//phi alpha rho u n
1499                         Fij(1+_Ndim)=phialphaq2n;
1500                         for(int idim=0;idim<_Ndim;idim++)
1501                                 Fij(2+_Ndim+idim)=phialphaq2n*_Vj[2+_Ndim+idim]+((1-alpha_roe)*_Vi[1]*_porosityi-dpi*_Vj[0]*_porosityj)*_vec_normal[idim];
1502                 }
1503                 return Fij;
1504         }
1505 }
1506
1507 void IsothermalTwoFluid::applyVFRoeLowMachCorrections(bool isBord, string groupname)
1508 {
1509         if(_nonLinearFormulation!=VFRoe)
1510                 throw CdmathException("IsothermalTwoFluid::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
1511         else//_nonLinearFormulation==VFRoe
1512         {
1513                 if(_spaceScheme==lowMach){
1514                         double u1_2=0, u2_2=0;
1515                         for(int i=0;i<_Ndim;i++){
1516                                 u1_2 += _Uroe[2+i]*_Uroe[2+i];
1517                                 u2_2 += _Uroe[2+i+_Ndim]*_Uroe[2+i+_Ndim];
1518                         }
1519
1520                         double  c = _maxvploc;//mixture sound speed
1521                         double M=max(sqrt(u1_2),sqrt(u2_2))/c;//Mach number
1522                         _Vij[1]=M*_Vij[1]+(1-M)*(_Vi[1]+_Vj[1])/2;
1523                         primToCons(_Vij,0,_Uij,0);
1524                 }
1525                 else if(_spaceScheme==pressureCorrection)
1526                 {
1527                         if(_pressureCorrectionOrder>2)
1528                                 throw CdmathException("IsothermalTwoFluid::applyVFRoeLowMachCorrections pressure correction order can be only 1 or 2 for Isothermal two-fluid model");
1529
1530                         double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;//mean velocities
1531                         double rho1 =  _fluides[0]->getDensity(_Uroe[1],_Temperature);
1532                         double rho2 =  _fluides[1]->getDensity(_Uroe[1],_Temperature);
1533                         double m1=_Uroe[0]*rho1,  m2=(1-_Uroe[0])*rho2;
1534                         double rhom=m1+m2;
1535                         for(int i=0;i<_Ndim;i++)
1536                         {
1537                                 norm_uij += (m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim])*(m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim]);
1538                                 uij_n += (m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim])*_vec_normal[i];
1539                                 ui_n += _Vi[2+i]*_vec_normal[i];
1540                                 uj_n += _Vj[2+i]*_vec_normal[i];
1541                         }
1542                         norm_uij=sqrt(norm_uij)/rhom;
1543                         uij_n/=rhom;
1544                         if(norm_uij>_precision)//avoid division by zero
1545                                 _Vij[1]=(_Vi[1]+_Vj[1])/2 + uij_n/norm_uij*(_Vj[1]-_Vi[1])/4 - rhom*norm_uij*(uj_n-ui_n)/4;
1546                         else
1547                                 _Vij[1]=(_Vi[1]+_Vj[1])/2                                    - rhom*norm_uij*(uj_n-ui_n)/4;
1548
1549                 }
1550                 else if(_spaceScheme==staggered)
1551                 {
1552                         double qij_n=0;
1553                         double rho1 =  _fluides[0]->getDensity(_Uroe[1],_Uroe[_nVar-1]);
1554                         double rho2 =  _fluides[1]->getDensity(_Uroe[1],_Uroe[_nVar-1]);
1555                         double m1=_Uroe[0]*rho1,  m2=(1-_Uroe[0])*rho2;
1556                         double rhom=m1+m2;
1557                         for(int i=0;i<_Ndim;i++)
1558                                 qij_n += (m1*_Uroe[2+i]+m2*_Uroe[2+i+_Ndim])*_vec_normal[i];
1559
1560                         if(qij_n>=0){
1561                                 _Vij[0]=_Vi[0];
1562                                 _Vij[1]=_Vj[1];
1563                                 for(int i=0;i<_Ndim;i++)
1564                                 {
1565                                         _Vij[2+i]               =_Vi[2+i];
1566                                         _Vij[2+i+_Ndim] =_Vi[2+i+_Ndim];
1567                                 }
1568                         }
1569                         else{
1570                                 _Vij[0]=_Vj[0];
1571                                 _Vij[1]=_Vi[1];
1572                                 for(int i=0;i<_Ndim;i++)
1573                                 {
1574                                         _Vij[2+i]               =_Vj[2+i];
1575                                         _Vij[2+i+_Ndim] =_Vj[2+i+_Ndim];
1576                                 }
1577                         }
1578                         primToCons(_Vij,0,_Uij,0);
1579                 }
1580         }
1581 }
1582
1583 void IsothermalTwoFluid::testConservation()
1584 {
1585         double SUM, DELTA, x;
1586         int I;
1587         for(int i=0; i<_nVar; i++)
1588         {
1589                 {
1590                         if(i == 0)
1591                                 cout << "Masse totale phase " << 0 <<" (kg): ";
1592                         else if( i == 1+_Ndim)
1593                                 cout << "Masse totale phase " << 1 <<" (kg): ";
1594                         else
1595                         {
1596                                 if(i < 1+_Ndim)
1597                                         cout << "Quantite de mouvement totale phase 0 (kg.m.s^-1): ";
1598                                 else
1599                                         cout << "Quantite de mouvement totale phase 1 (kg.m.s^-1): ";
1600                         }
1601                 }
1602                 SUM = 0;
1603                 DELTA = 0;
1604                 I = i;
1605                 for(int j=0; j<_Nmailles; j++)
1606                 {
1607                         VecGetValues(_old, 1, &I, &x);//on recupere la valeur du champ
1608                         SUM += x*_mesh.getCell(j).getMeasure();
1609                         VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
1610                         DELTA += x*_mesh.getCell(j).getMeasure();
1611                         I += _nVar;
1612                 }
1613                 if(fabs(SUM)>_precision)
1614                         cout << SUM<< ", variation relative: " << fabs(DELTA /SUM)  << endl;
1615                 else
1616                         cout << " a une somme nulle,  variation absolue: " << fabs(DELTA) << endl;
1617         }
1618 }
1619
1620 void IsothermalTwoFluid::save(){
1621         string prim(_path+"/IsothermalTwoFluidPrim_");
1622         string cons(_path+"/IsothermalTwoFluidCons_");
1623         prim+=_fileName;
1624         cons+=_fileName;
1625
1626         PetscInt Ii;
1627         for (long i = 0; i < _Nmailles; i++){
1628                 /* j = 0 : void fraction
1629                            j = 1 : pressure
1630                            j = 2, 3, 4: velocity phase 1
1631                            j = 5, 6, 7: velocity phase 2 */
1632                 for (int j = 0; j < _nVar; j++){
1633                         Ii = i*_nVar +j;
1634                         VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
1635                 }
1636         }
1637         if(_saveConservativeField){
1638                 for (long i = 0; i < _Nmailles; i++){
1639                         for (int j = 0; j < _nVar; j++){
1640                                 Ii = i*_nVar +j;
1641                                 //                              cout<<"i= "<<i<< " j= "<<j<<" UU(i,j)= "<<_UU(i,j)<<endl;
1642                                 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
1643                         }
1644                 }
1645                 _UU.setTime(_time,_nbTimeStep);
1646         }
1647         _VV.setTime(_time,_nbTimeStep);
1648         if (_nbTimeStep ==0 || _restartWithNewFileName){
1649                 string prim_suppress ="rm -rf "+prim+"_*";
1650                 string cons_suppress ="rm -rf "+cons+"_*";
1651                 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
1652                 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
1653                 _VV.setInfoOnComponent(0,"Void_fraction");
1654                 _VV.setInfoOnComponent(1,"Pressure_(Pa)");
1655                 _VV.setInfoOnComponent(2,"Velocity1_x_m/s");
1656                 if (_Ndim>1)
1657                         _VV.setInfoOnComponent(3,"Velocity1_y_m/s");
1658                 if (_Ndim>2)
1659                         _VV.setInfoOnComponent(4,"Velocity1_z_m/s");
1660
1661                 _VV.setInfoOnComponent(2+_Ndim,"Velocity2_x_m/s");
1662                 if (_Ndim>1)
1663                         _VV.setInfoOnComponent(3+_Ndim,"Velocity2_y_m/s");
1664                 if (_Ndim>2)
1665                         _VV.setInfoOnComponent(4+_Ndim,"Velocity2_z_m/s");
1666
1667                 switch(_saveFormat)
1668                 {
1669                 case VTK :
1670                         _VV.writeVTK(prim);
1671                         break;
1672                 case MED :
1673                         _VV.writeMED(prim);
1674                         break;
1675                 case CSV :
1676                         _VV.writeCSV(prim);
1677                         break;
1678                 }
1679
1680                 if(_saveConservativeField){
1681                         _UU.setInfoOnComponent(0,"Partial_density1");// (kg/m^3)
1682                         _UU.setInfoOnComponent(1,"Momentum1_x");// phase1  (kg/m^2/s)
1683                         if (_Ndim>1)
1684                                 _UU.setInfoOnComponent(2,"Momentum1_y");// phase1 (kg/m^2/s)
1685                         if (_Ndim>2)
1686                                 _UU.setInfoOnComponent(3,"Momentum1_z");// phase1 (kg/m^2/s)
1687
1688                         _UU.setInfoOnComponent(1+_Ndim,"Partial_density2");// phase2 (kg/m^3)
1689                         _UU.setInfoOnComponent(2+_Ndim,"Momentum2_x");// phase2 (kg/m^2/s)
1690                         if (_Ndim>1)
1691                                 _UU.setInfoOnComponent(3+_Ndim,"Momentum2_y");// phase2 (kg/m^2/s)
1692                         if (_Ndim>2)
1693                                 _UU.setInfoOnComponent(4+_Ndim,"Momentum2_z");// phase2 (kg/m^2/s)
1694
1695                         switch(_saveFormat)
1696                         {
1697                         case VTK :
1698                                 _UU.writeVTK(cons);
1699                                 break;
1700                         case MED :
1701                                 _UU.writeMED(cons);
1702                                 break;
1703                         case CSV :
1704                                 _UU.writeCSV(cons);
1705                                 break;
1706                         }
1707                 }
1708         }
1709         else{
1710                 switch(_saveFormat)
1711                 {
1712                 case VTK :
1713                         _VV.writeVTK(prim,false);
1714                         break;
1715                 case MED :
1716                         _VV.writeMED(prim,false);
1717                         break;
1718                 case CSV :
1719                         _VV.writeCSV(prim);
1720                         break;
1721                 }
1722
1723                 if(_saveConservativeField){
1724                         switch(_saveFormat)
1725                         {
1726                         case VTK :
1727                                 _UU.writeVTK(cons,false);
1728                                 break;
1729                         case MED :
1730                                 _UU.writeMED(cons,false);
1731                                 break;
1732                         case CSV :
1733                                 _UU.writeCSV(cons);
1734                                 break;
1735                         }
1736                 }
1737         }
1738         if(_saveVelocity){
1739                 for (long i = 0; i < _Nmailles; i++){
1740                         // j = 0 : concentration, j=1 : pressure; j = _nVar - 1: temperature; j = 2,..,_nVar-2: velocity
1741                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
1742                         {
1743                                 Ii = i*_nVar +2+j;
1744                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse1(i,j));
1745                                 Ii=i*_nVar +2+j+_Ndim;
1746                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse2(i,j));
1747                         }
1748                         for (int j = _Ndim; j < 3; j++){//On met à zero les composantes de vitesse si la dimension est <3
1749                                 _Vitesse1(i,j)=0;
1750                                 _Vitesse2(i,j)=0;
1751                         }
1752                 }
1753                 _Vitesse1.setTime(_time,_nbTimeStep);
1754                 _Vitesse2.setTime(_time,_nbTimeStep);
1755                 if (_nbTimeStep ==0 || _restartWithNewFileName){
1756                         _Vitesse1.setInfoOnComponent(0,"Velocity_x_(m/s)");
1757                         _Vitesse1.setInfoOnComponent(1,"Velocity_y_(m/s)");
1758                         _Vitesse1.setInfoOnComponent(2,"Velocity_z_(m/s)");
1759
1760                         _Vitesse2.setInfoOnComponent(0,"Velocity_x_(m/s)");
1761                         _Vitesse2.setInfoOnComponent(1,"Velocity_y_(m/s)");
1762                         _Vitesse2.setInfoOnComponent(2,"Velocity_z_(m/s)");
1763
1764                         switch(_saveFormat)
1765                         {
1766                         case VTK :
1767                                 _Vitesse1.writeVTK(prim+"_GasVelocity");
1768                                 _Vitesse2.writeVTK(prim+"_LiquidVelocity");
1769                                 break;
1770                         case MED :
1771                                 _Vitesse1.writeMED(prim+"_GasVelocity");
1772                                 _Vitesse2.writeMED(prim+"_LiquidVelocity");
1773                                 break;
1774                         case CSV :
1775                                 _Vitesse1.writeCSV(prim+"_GasVelocity");
1776                                 _Vitesse2.writeCSV(prim+"_LiquidVelocity");
1777                                 break;
1778                         }
1779                 }
1780                 else{
1781                         switch(_saveFormat)
1782                         {
1783                         case VTK :
1784                                 _Vitesse1.writeVTK(prim+"_GasVelocity",false);
1785                                 _Vitesse2.writeVTK(prim+"_LiquidVelocity",false);
1786                                 break;
1787                         case MED :
1788                                 _Vitesse1.writeMED(prim+"_GasVelocity",false);
1789                                 _Vitesse2.writeMED(prim+"_LiquidVelocity",false);
1790                                 break;
1791                         case CSV :
1792                                 _Vitesse1.writeCSV(prim+"_GasVelocity");
1793                                 _Vitesse2.writeCSV(prim+"_LiquidVelocity");
1794                                 break;
1795                         }
1796                 }
1797         }
1798
1799         if (_restartWithNewFileName)
1800                 _restartWithNewFileName=false;
1801 }
1802