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