CoreFlows/Models/src/DriftModel.cxx

   1 /*
   2  * DriftModel.cxx
   3  *
   4  *  Created on: 1 janv. 2015
   5  *      Author: Michael Ndjinga
   6  */
   7 #include "DriftModel.hxx"
   8 #include "StiffenedGas.hxx"
   9
  10 using namespace std;
  11
  12 DriftModel::DriftModel(pressureEstimate pEstimate, int dim, bool useDellacherieEOS){
  13         _Ndim=dim;
  14         _nVar=_Ndim+3;
  15         _nbPhases  = 2;
  16         _dragCoeffs=vector<double>(2,0);
  17         _fluides.resize(2);
  18         _saveAllFields=false;
  19
  20         if( pEstimate==around1bar300K){//EOS at 1 bar and 373K
  21                 cout<<"Fluid is water-Gas mixture around saturation point 1 bar and 373 K (100°C)"<<endl;
  22                 *_runLogFile<<"Fluid is water-Gas mixture around saturation point 1 bar and 373 K (100°C)"<<endl;
  23                 _Tsat=373;//saturation temperature at 1 bar
  24                 double esatv=2.505e6;//Gas internal energy at saturation at 1 bar
  25                 double esatl=4.174e5;//water internal energy at saturation at 1 bar
  26                 double cv_v=1555;//Gas specific heat capacity at saturation at 1 bar
  27                 double cv_l=3770;//water specific heat capacity at saturation at 1 bar
  28                 double gamma_v=1.337;//Gas heat capacity ratio at saturation at 1 bar
  29                 double rho_sat_l=958;//water density at saturation at 1 bar
  30                 double sound_speed_l=1543;//water sound speed at saturation at 1 bar
  31                 _fluides[0] = new StiffenedGas(gamma_v,cv_v,_Tsat,esatv);  //ideal gas law for Gas at pressure 1 bar and temperature 100°C, gamma=1.34
  32                 _fluides[1] = new StiffenedGas(rho_sat_l,1e5,_Tsat,esatl,sound_speed_l,cv_l);  //stiffened gas law for water at pressure 1 bar and temperature 100°C
  33                 _hsatl=4.175e5;//water enthalpy at saturation at 1 bar
  34                 _hsatv=2.675e6;//Gas enthalpy at saturation at 1 bar
  35
  36                 _useDellacherieEOS=false;
  37         }
  38         else{//EOS at 155 bars and 618K
  39                 cout<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
  40                 *_runLogFile<<"Fluid is water-Gas mixture around saturation point 155 bars and 618 K (345°C)"<<endl;
  41                 _useDellacherieEOS=useDellacherieEOS;
  42                 if(useDellacherieEOS)
  43                 {
  44                         _Tsat=656;//saturation temperature used in Dellacherie EOS
  45                         _hsatl=1.633e6;//water enthalpy at saturation at 155 bars
  46                         _hsatv=3.006e6;//Gas enthalpy at saturation at 155 bars
  47                         _fluides[0] = new StiffenedGasDellacherie(1.43,0  ,2.030255e6  ,1040.14); //stiffened gas law for Gas from S. Dellacherie
  48                         _fluides[1] = new StiffenedGasDellacherie(2.35,1e9,-1.167056e6,1816.2); //stiffened gas law for water from S. Dellacherie
  49                 }
  50                 else
  51                 {
  52                         double esatv=2.444e6;//Gas internal energy at saturation at 155 bar
  53                         double esatl=1.604e6;//water internal energy at saturation at 155 bar
  54                         double sound_speed_v=433;//Gas sound speed at saturation at 155 bar
  55                         double sound_speed_l=621;//water sound speed at saturation at 155 bar
  56                         double cv_v=3633;//Gas specific heat capacity at saturation at 155 bar
  57                         double cv_l=3100;//water specific heat capacity at saturation at 155 bar
  58                         double rho_sat_v=102;//Gas density at saturation at 155 bar
  59                         double rho_sat_l=594;//water density at saturation at 155 bar
  60                         _Tsat=618;//saturation temperature at 155 bars
  61                         _hsatl=1.63e6;//water enthalpy at saturation at 155 bars
  62                         _hsatv=2.6e6;//Gas enthalpy at saturation at 155 bars
  63                         _fluides[0] = new StiffenedGas(rho_sat_v,1.55e7,_Tsat,esatv, sound_speed_v,cv_v); //stiffened gas law for Gas at pressure 155 bar and temperature 345°C
  64                         _fluides[1] = new StiffenedGas(rho_sat_l,1.55e7,_Tsat,esatl, sound_speed_l,cv_l); //stiffened gas law for water at pressure 155 bar
  65                 }
  66         }
  67         _latentHeat=_hsatv-_hsatl;
  68         cout<<"Liquid saturation enthalpy "<< _hsatl<<" J/Kg"<<endl;
  69         *_runLogFile<<"Liquid saturation enthalpy "<< _hsatl<<" J/Kg"<<endl;
  70         cout<<"Vapour saturation enthalpy "<< _hsatv<<" J/Kg"<<endl;
  71         *_runLogFile<<"Vapour saturation enthalpy "<< _hsatv<<" J/Kg"<<endl;
  72         cout<<"Latent heat "<< _latentHeat<<endl;
  73         *_runLogFile<<"Latent heat "<< _latentHeat<<endl;
  74
  75         _fileName = "SolverlabDriftModel";
  76     PetscPrintf(PETSC_COMM_WORLD,"\n Drift model problem for two phase flow\n");
  77 }
  78
  79 void DriftModel::initialize(){
  80         cout<<"\n Initialising the drift model"<<endl;
  81         *_runLogFile<<"\n Initialising the drift model"<<endl;
  82
  83         _Uroe = new double[_nVar];
  84         _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
  85         for(int i=0; i<_Ndim; i++)
  86                 _gravite[i+2]=_GravityField3d[i];
  87
  88         _GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
  89
  90         if(_saveVelocity)
  91                 _Vitesse=Field("Velocity",CELLS,_mesh,3);//Forcement en dimension 3 (3 composantes) pour le posttraitement des lignes de courant
  92
  93         if(_saveAllFields)
  94         {
  95                 _VoidFraction=Field("Void fraction",CELLS,_mesh,1);
  96                 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
  97                 _Concentration=Field("Concentration",CELLS,_mesh,1);
  98                 _Pressure=Field("Pressure",CELLS,_mesh,1);
  99                 _mixtureDensity=Field("Mixt density",CELLS,_mesh,1);
 100                 _Temperature=Field("Temperature",CELLS,_mesh,1);
 101                 _DensiteLiquide=Field("Liquid density",CELLS,_mesh,1);
 102                 _DensiteVapeur=Field("Steam density",CELLS,_mesh,1);
 103                 _EnthalpieLiquide=Field("Liquid enthalpy",CELLS,_mesh,1);
 104                 _EnthalpieVapeur=Field("Steam enthalpy",CELLS,_mesh,1);
 105                 _VitesseX=Field("Velocity x",CELLS,_mesh,1);
 106                 if(_Ndim>1)
 107                 {
 108                         _VitesseY=Field("Velocity y",CELLS,_mesh,1);
 109                         if(_Ndim>2)
 110                                 _VitesseZ=Field("Velocity z",CELLS,_mesh,1);
 111                 }
 112         }
 113
 114         if(_entropicCorrection)
 115                 _entropicShift=vector<double>(3);//at most 3 distinct eigenvalues
 116
 117         ProblemFluid::initialize();
 118 }
 119
 120 bool DriftModel::iterateTimeStep(bool &converged)
 121 {
 122         if(_timeScheme == Explicit || !_usePrimitiveVarsInNewton)
 123                 return ProblemFluid::iterateTimeStep(converged);
 124         else
 125         {
 126                 bool stop=false;
 127
 128                 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
 129                         _maxvp=0;
 130                         computeTimeStep(stop);
 131                 }
 132                 if(stop){//Le compute time step ne s'est pas bien passé
 133                         cout<<"ComputeTimeStep failed"<<endl;
 134                         converged=false;
 135                         return false;
 136                 }
 137
 138                 computeNewtonVariation();
 139
 140                 //converged=convergence des iterations
 141                 if(_timeScheme == Explicit)
 142                         converged=true;
 143                 else{//Implicit scheme
 144
 145                         KSPGetIterationNumber(_ksp, &_PetscIts);
 146                         if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
 147                                 _MaxIterLinearSolver = _PetscIts;
 148                         if(_PetscIts>=_maxPetscIts)//solving the linear system failed
 149                         {
 150                                 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 151                                 //*_runLogFileogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 152                                 converged=false;
 153                                 return false;
 154                         }
 155                         else{//solving the linear system succeeded
 156                                 //Calcul de la variation relative Uk+1-Uk
 157                                 _erreur_rel = 0.;
 158                                 double x, dx;
 159                                 int I;
 160                                 for(int j=1; j<=_Nmailles; j++)
 161                                 {
 162                                         for(int k=0; k<_nVar; k++)
 163                                         {
 164                                                 I = (j-1)*_nVar + k;
 165                                                 VecGetValues(_newtonVariation, 1, &I, &dx);
 166                                                 VecGetValues(_primitiveVars, 1, &I, &x);
 167                                                 if (fabs(x)*fabs(x)< _precision)
 168                                                 {
 169                                                         if(_erreur_rel < fabs(dx))
 170                                                                 _erreur_rel = fabs(dx);
 171                                                 }
 172                                                 else if(_erreur_rel < fabs(dx/x))
 173                                                         _erreur_rel = fabs(dx/x);
 174                                         }
 175                                 }
 176                         }
 177                         converged = _erreur_rel <= _precision_Newton;
 178                 }
 179
 180                 double relaxation=1;//Vk+1=Vk+relaxation*deltaV
 181
 182                 VecAXPY(_primitiveVars,  relaxation, _newtonVariation);
 183
 184                 //mise a jour du champ primitif
 185                 updateConservatives();
 186
 187                 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
 188                 {
 189                         cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 190                         *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 191                         converged=false;
 192                         return false;
 193                 }
 194                 if(_system)
 195                 {
 196                         cout<<"Vecteur Vkp1-Vk "<<endl;
 197                         VecView(_newtonVariation,  PETSC_VIEWER_STDOUT_SELF);
 198                         cout << "Nouvel etat courant Vk de l'iteration Newton: " << endl;
 199                         VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_SELF);
 200                 }
 201
 202                 if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
 203                         if(_minm1<-_precision || _minm2<-_precision)
 204                         {
 205                                 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 206                                 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 207                         }
 208
 209                         if (_nbVpCplx>0){
 210                                 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 211                                 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 212                         }
 213                 }
 214                 _minm1=1e30;
 215                 _minm2=1e30;
 216                 _nbMaillesNeg=0;
 217                 _nbVpCplx =0;
 218                 _part_imag_max=0;
 219
 220                 return true;
 221         }
 222 }
 223 void DriftModel::computeNewtonVariation()
 224 {
 225         if(_timeScheme==Explicit || !_usePrimitiveVarsInNewton)
 226                 ProblemFluid::computeNewtonVariation();
 227         else
 228         {
 229                 if(_verbose)
 230                 {
 231                         cout<<"Vecteur courant Vk "<<endl;
 232                         VecView(_primitiveVars,PETSC_VIEWER_STDOUT_SELF);
 233                         cout << endl;
 234                 }
 235                 if(_timeScheme == Explicit)
 236                 {
 237                         VecCopy(_b,_newtonVariation);
 238                         VecScale(_newtonVariation, _dt);
 239                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 240                         {
 241                                 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
 242                                 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
 243                                 cout << endl;
 244                         }
 245                 }
 246                 else
 247                 {
 248                         if(_verbose)
 249                         {
 250                                 cout << "Matrice du système linéaire avant contribution delta t" << endl;
 251                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 252                                 cout << endl;
 253                                 cout << "Second membre du système linéaire avant contribution delta t" << endl;
 254                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 255                                 cout << endl;
 256                         }
 257                         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 258
 259                         VecAXPY(_b, 1/_dt, _old);
 260                         VecAXPY(_b, -1/_dt, _conservativeVars);
 261
 262                         for(int imaille = 0; imaille<_Nmailles; imaille++){
 263                                 _idm[0] = _nVar*imaille;
 264                                 for(int k=1; k<_nVar; k++)
 265                                         _idm[k] = _idm[k-1] + 1;
 266                                 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
 267                                 primToConsJacobianMatrix(_Vi);
 268                                 for(int k=0; k<_nVar*_nVar; k++)
 269                                         _primToConsJacoMat[k]*=1/_dt;
 270                                 MatSetValuesBlocked(_A, 1, &imaille, 1, &imaille, _primToConsJacoMat, ADD_VALUES);
 271                         }
 272                         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 273
 274 #if PETSC_VERSION_GREATER_3_5
 275                         KSPSetOperators(_ksp, _A, _A);
 276 #else
 277                         KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 278 #endif
 279
 280                         if(_verbose)
 281                         {
 282                                 cout << "Matrice du système linéaire après contribution delta t" << endl;
 283                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 284                                 cout << endl;
 285                                 cout << "Second membre du système linéaire après contribution delta t" << endl;
 286                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 287                                 cout << endl;
 288                         }
 289
 290                         if(_conditionNumber)
 291                                 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 292                         if(!_isScaling)
 293                         {
 294                                 KSPSolve(_ksp, _b, _newtonVariation);
 295                         }
 296                         else
 297                         {
 298                                 computeScaling(_maxvp);
 299                                 int indice;
 300                                 VecAssemblyBegin(_vecScaling);
 301                                 VecAssemblyBegin(_invVecScaling);
 302                                 for(int imaille = 0; imaille<_Nmailles; imaille++)
 303                                 {
 304                                         indice = imaille;
 305                                         VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
 306                                         VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
 307                                 }
 308                                 VecAssemblyEnd(_vecScaling);
 309                                 VecAssemblyEnd(_invVecScaling);
 310
 311                                 if(_system)
 312                                 {
 313                                         cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
 314                                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 315                                         cout << endl;
 316                                         cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
 317                                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 318                                         cout << endl;
 319                                 }
 320                                 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
 321                                 if(_system)
 322                                 {
 323                                         cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
 324                                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 325                                         cout << endl;
 326                                 }
 327                                 VecCopy(_b,_bScaling);
 328                                 VecPointwiseMult(_b,_vecScaling,_bScaling);
 329                                 if(_system)
 330                                 {
 331                                         cout << "Produit du second membre par le preconditionneur bloc diagonal  a gauche" << endl;
 332                                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 333                                         cout << endl;
 334                                 }
 335
 336                                 KSPSolve(_ksp,_b, _bScaling);
 337                                 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
 338                         }
 339                         if(_verbose)
 340                         {
 341                                 cout << "solution du systeme lineaire local:" << endl;
 342                                 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
 343                                 cout << endl;
 344                         }
 345                 }
 346         }
 347 }
 348
 349 void DriftModel::convectionState( const long &i, const long &j, const bool &IsBord){
 350         //      sortie: etat de Roe rho, cm, v, H,T
 351
 352         _idm[0] = _nVar*i;
 353         for(int k=1; k<_nVar; k++)
 354                 _idm[k] = _idm[k-1] + 1;
 355         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 356
 357         _idm[0] = _nVar*j;
 358         for(int k=1; k<_nVar; k++)
 359                 _idm[k] = _idm[k-1] + 1;
 360         if(IsBord)
 361                 VecGetValues(_Uext, _nVar, _idm, _Uj);
 362         else
 363                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 364         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 365         {
 366                 cout<<"DriftModel::convectionState Left state cell " << i<< ": "<<endl;
 367                 for(int k =0; k<_nVar; k++)
 368                         cout<< _Ui[k]<<endl;
 369                 cout<<"DriftModel::convectionState Right state cell " << j<< ": "<<endl;
 370                 for(int k =0; k<_nVar; k++)
 371                         cout<< _Uj[k]<<endl;
 372         }
 373         if(_Ui[0]<0||_Uj[0]<0)
 374         {
 375                 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!densite de melange negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 376                 *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!!densite de melange negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 377                 _runLogFile->close();
 378                 throw CdmathException("densite negative, arret de calcul");
 379         }
 380         PetscScalar ri, rj, xi, xj, pi, pj;
 381         PetscInt Ii;
 382         ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
 383         rj = sqrt(_Uj[0]);
 384
 385         _Uroe[0] = ri*rj/sqrt(_porosityi*_porosityj);   //moyenne geometrique des densites
 386         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 387                 cout << "Densité moyenne Roe gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
 388         xi = _Ui[1]/_Ui[0];//cm gauche
 389         xj = _Uj[1]/_Uj[0];//cm droite
 390
 391         _Uroe[1] = (xi*ri + xj*rj)/(ri + rj);//moyenne de Roe des concentrations
 392         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 393                 cout << "Concentration de Roe  gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[1] << endl;
 394         for(int k=0;k<_Ndim;k++){
 395                 xi = _Ui[k+2];//phi rho u gauche
 396                 xj = _Uj[k+2];//phi rho u droite
 397                 _Uroe[2+k] = (xi/ri + xj/rj)/(ri + rj);
 398                 //"moyenne" des vitesses
 399                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 400                         cout << "Vitesse de Roe composante "<< k<<"  gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+2] << endl;
 401         }
 402         // H = (rho E + p)/rho
 403         xi = _Ui[_nVar-1];//phi rho E
 404         xj = _Uj[_nVar-1];
 405         Ii = i*_nVar+1; // correct Kieu
 406         VecGetValues(_primitiveVars, 1, &Ii, &pi);// _primitiveVars pour p
 407         if(IsBord)
 408         {
 409                 consToPrim(_Uj,_Vj,_porosityj);
 410                 pj =  _Vj[1];
 411         }
 412         else
 413         {
 414                 Ii = j*_nVar+1; // correct Kieu
 415                 VecGetValues(_primitiveVars, 1, &Ii, &pj);
 416         }
 417         xi = (xi + pi)/_Ui[0];
 418         xj = (xj + pj)/_Uj[0];
 419         _Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
 420         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 421                 cout << "Enthalpie totale de Roe H  gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
 422         // Moyenne de Roe de Tg et Td
 423         Ii = i*_nVar+_nVar-1;
 424         VecGetValues(_primitiveVars, 1, &Ii, &xi);// _primitiveVars pour T
 425         if(IsBord)
 426         {
 427                 //consToPrim(_Uj,_Vj,_porosityj);//Fonction déjà appelée
 428                 xj =  _Vj[_nVar-1];
 429         }
 430         else
 431         {
 432                 Ii = j*_nVar+_nVar-1;
 433                 VecGetValues(_primitiveVars, 1, &Ii, &xj);
 434         }
 435         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 436         {
 437                 cout<<"Convection interfacial state"<<endl;
 438                 for(int k=0;k<_nVar;k++)
 439                         cout<< _Uroe[k]<<" , "<<endl;
 440         }
 441 }
 442
 443 void DriftModel::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
 444         //sortie: matrices et etat de diffusion (rho, rho cm, q, T)
 445
 446         _idm[0] = _nVar*i;
 447         for(int k=1; k<_nVar; k++)
 448                 _idm[k] = _idm[k-1] + 1;
 449
 450         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 451         _idm[0] = _nVar*j;
 452         for(int k=1; k<_nVar; k++)
 453                 _idm[k] = _idm[k-1] + 1;
 454
 455         if(IsBord)
 456                 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
 457         else
 458                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 459
 460         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 461         {
 462                 cout << "DriftModel::diffusionStateAndMatrices cellule gauche" << i << endl;
 463                 cout << "Ui = ";
 464                 for(int q=0; q<_nVar; q++)
 465                         cout << _Ui[q]  << "\t";
 466                 cout << endl;
 467                 cout << "DriftModel::diffusionStateAndMatrices cellule droite" << j << endl;
 468                 cout << "Uj = ";
 469                 for(int q=0; q<_nVar; q++)
 470                         cout << _Uj[q]  << "\t";
 471                 cout << endl;
 472         }
 473
 474         for(int k=0; k<_nVar; k++)
 475                 _Udiff[k] = (_Ui[k]/_porosityi+_Uj[k]/_porosityj)/2;
 476         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 477         {
 478                 cout << "DriftModel::diffusionStateAndMatrices conservative diffusion state" << endl;
 479                 cout << "_Udiff = ";
 480                 for(int q=0; q<_nVar; q++)
 481                         cout << _Udiff[q]  << "\t";
 482                 cout << endl;
 483                 cout << "porosite gauche= "<<_porosityi<< ", porosite droite= "<<_porosityj<<endl;
 484         }
 485
 486         consToPrim(_Udiff,_phi,1);
 487         _Udiff[_nVar-1]=_phi[_nVar-1];
 488
 489         double Tm=_phi[_nVar-1];
 490         double RhomEm=_Udiff[_nVar-1];
 491         _Udiff[_nVar-1]=Tm;
 492
 493         if(_timeScheme==Implicit)
 494         {
 495                 double q_2=0;
 496                 for (int l = 0; l<_Ndim;l++)
 497                         q_2+=_Udiff[l+2]*_Udiff[l+2];
 498                 double pression=_phi[1];
 499                 double m_v=_Udiff[1];
 500                 double m_l=_Udiff[0]-_Udiff[1];
 501                 double rho_v=_fluides[0]->getDensity(pression,Tm);
 502                 double rho_l=_fluides[1]->getDensity(pression,Tm);
 503                 double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
 504
 505                 for(int l=0; l<_nVar*_nVar;l++)
 506                         _Diffusion[l] = 0;
 507                 double mu = alpha_v*_fluides[0]->getViscosity(Tm)+alpha_l*_fluides[1]->getViscosity(Tm);
 508                 for(int l=2;l<(_nVar-1);l++)
 509                 {
 510                         _Diffusion[l*_nVar] =  mu*_Udiff[l]/(_Udiff[0]*_Udiff[0]);
 511                         _Diffusion[l*_nVar+l] = -mu/_Udiff[0];
 512                 }
 513                 double lambda = alpha_v*_fluides[0]->getConductivity(Tm)+alpha_l*_fluides[1]->getConductivity(Tm);
 514                 double C_v=  alpha_v*_fluides[0]->constante("cv");
 515                 double C_l=      alpha_l*_fluides[1]->constante("cv");
 516                 double ev0=_fluides[0]->getInternalEnergy(0,rho_v);//Corriger
 517                 double el0=_fluides[1]->getInternalEnergy(0,rho_l);//Corriger
 518                 double Rhomem=RhomEm-0.5*q_2/(_Udiff[0]*_Udiff[0]);
 519                 int q = (_nVar-1)*_nVar;
 520                 //Formules a verifier
 521                 _Diffusion[q]=lambda*((Rhomem-m_v*ev0-m_l*el0)*C_l/((m_v*C_v+m_l*C_l)*(m_v*C_v+m_l*C_l))+el0/(m_v*C_v+m_l*C_l)-q_2/(2*_Udiff[0]*_Udiff[0]*(m_v*C_v+m_l*C_l)));
 522                 _Diffusion[q+1]=lambda*((Rhomem-m_v*ev0-m_l*el0)*(C_v-C_l)/((m_v*C_v+m_l*C_l)*(m_v*C_v+m_l*C_l))+(ev0+el0)/(m_v*C_v+m_l*C_l));
 523                 for(int k=2;k<(_nVar-1);k++)
 524                 {
 525                         _Diffusion[q+k]= lambda*_Udiff[k]/(_Udiff[0]*(m_v*C_v+m_l*C_l));
 526                 }
 527                 _Diffusion[_nVar*_nVar-1]=-lambda/(m_v*C_v+m_l*C_l);
 528                 /*Affichages */
 529                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 530                 {
 531                         cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
 532                         displayMatrix( _Diffusion,_nVar," Matrice de diffusion ");
 533                 }
 534         }
 535 }
 536
 537 void DriftModel::setBoundaryState(string nameOfGroup, const int &j,double *normale){
 538         int k;
 539         double v2=0, q_n=0;//q_n=quantité de mouvement normale à la face interne
 540         _idm[0] = _nVar*j;
 541         for(k=1; k<_nVar; k++)
 542                 _idm[k] = _idm[k-1] + 1;
 543         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
 544         for(k=0; k<_Ndim; k++)
 545                 q_n+=_externalStates[(k+2)]*normale[k];
 546
 547         double porosityj=_porosityField(j);
 548
 549         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 550         {
 551                 cout << "setBoundaryState for group "<< nameOfGroup<< ", inner cell j= "<<j << " face unit normal vector "<<endl;
 552                 for(k=0; k<_Ndim; k++){
 553                         cout<<normale[k]<<", ";
 554                 }
 555                 cout<<endl;
 556         }
 557
 558         if (_limitField[nameOfGroup].bcType==Wall || _limitField[nameOfGroup].bcType==InnerWall){
 559                 //Pour la convection, inversion du sens de la vitesse normale
 560                 for(k=0; k<_Ndim; k++)
 561                         _externalStates[(k+2)]-= 2*q_n*normale[k];
 562
 563                 _idm[0] = 0;
 564                 for(k=1; k<_nVar; k++)
 565                         _idm[k] = _idm[k-1] + 1;
 566
 567                 VecAssemblyBegin(_Uext);
 568                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 569                 VecAssemblyEnd(_Uext);
 570
 571                 //Pour la diffusion, paroi à vitesse et temperature imposees
 572                 _idm[0] = _nVar*j;
 573                 for(k=1; k<_nVar; k++)
 574                         _idm[k] = _idm[k-1] + 1;
 575                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 576                 double concentration=_Vj[0];
 577                 double pression=_Vj[1];
 578                 double T=_limitField[nameOfGroup].T;
 579                 double rho_v=_fluides[0]->getDensity(pression,T);
 580                 double rho_l=_fluides[1]->getDensity(pression,T);
 581                 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
 582                 {
 583                         cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
 584                         cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 585                         *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 586                         _runLogFile->close();
 587                         throw CdmathException("DriftModel::setBoundaryState: Inlet, impossible to compute mixture density, division by zero");
 588                 }
 589
 590                 _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
 591                 _externalStates[1]=concentration*_externalStates[0];
 592
 593                 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
 594                 v2 +=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
 595                 if(_Ndim>1)
 596                 {
 597                         v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 598                         _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
 599                         if(_Ndim==3)
 600                         {
 601                                 _externalStates[4]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
 602                                 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 603                         }
 604                 }
 605                 _externalStates[_nVar-1] = _externalStates[1]*_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho_v)
 606                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          +(_externalStates[0]-_externalStates[1])*_fluides[1]->getInternalEnergy(_limitField[nameOfGroup].T,rho_l) + _externalStates[0]*v2/2;
 607                 _idm[0] = 0;
 608                 for(k=1; k<_nVar; k++)
 609                         _idm[k] = _idm[k-1] + 1;
 610                 VecAssemblyBegin(_Uextdiff);
 611                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 612                 VecAssemblyEnd(_Uextdiff);
 613         }
 614         else if (_limitField[nameOfGroup].bcType==Neumann){
 615                 _idm[0] = 0;
 616                 for(k=1; k<_nVar; k++)
 617                         _idm[k] = _idm[k-1] + 1;
 618
 619                 VecAssemblyBegin(_Uext);
 620                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 621                 VecAssemblyEnd(_Uext);
 622
 623                 VecAssemblyBegin(_Uextdiff);
 624                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 625                 VecAssemblyEnd(_Uextdiff);
 626         }
 627         else if (_limitField[nameOfGroup].bcType==Inlet){
 628
 629                 if(q_n<=0){
 630                         VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 631                         double concentration=_limitField[nameOfGroup].conc;
 632                         double pression=_Vj[1];
 633                         double T=_limitField[nameOfGroup].T;
 634                         double rho_v=_fluides[0]->getDensity(pression,T);
 635                         double rho_l=_fluides[1]->getDensity(pression,T);
 636                         if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
 637                         {
 638                                 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
 639                                 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 640                                 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 641                                 _runLogFile->close();
 642                                 throw CdmathException("DriftModel::setBoundaryState: Inlet, impossible to compute mixture density, division by zero");
 643                         }
 644
 645                         _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
 646                         _externalStates[1]=concentration*_externalStates[0];
 647                         _externalStates[2]=_externalStates[0]*(_limitField[nameOfGroup].v_x[0]);
 648                         v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
 649                         if(_Ndim>1)
 650                         {
 651                                 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 652                                 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
 653                                 if(_Ndim==3)
 654                                 {
 655                                         _externalStates[4]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
 656                                         v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 657                                 }
 658                         }
 659                         _externalStates[_nVar-1] = _externalStates[1]*_fluides[0]->getInternalEnergy(T,rho_v)+(_externalStates[0]-_externalStates[1])*_fluides[1]->getInternalEnergy(T,rho_l) + _externalStates[0]*v2/2;
 660                 }
 661                 else if((_nbTimeStep-1)%_freqSave ==0)
 662                         cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
 663
 664                 _idm[0] = 0;
 665                 for(k=1; k<_nVar; k++)
 666                         _idm[k] = _idm[k-1] + 1;
 667
 668                 VecAssemblyBegin(_Uext);
 669                 VecAssemblyBegin(_Uextdiff);
 670                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 671                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 672                 VecAssemblyEnd(_Uext);
 673                 VecAssemblyEnd(_Uextdiff);
 674         }
 675         else if (_limitField[nameOfGroup].bcType==InletPressure){
 676
 677                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 678                 Cell Cj=_mesh.getCell(j);
 679                 double hydroPress=Cj.x()*_GravityField3d[0];
 680                 if(_Ndim>1){
 681                         hydroPress+=Cj.y()*_GravityField3d[1];
 682                         if(_Ndim>2)
 683                                 hydroPress+=Cj.z()*_GravityField3d[2];
 684                 }
 685                 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho
 686
 687                 //Building the external state
 688                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 689                 double concentration, Tm;
 690                 if(q_n<=0){
 691                         concentration=_limitField[nameOfGroup].conc;
 692                         Tm=_limitField[nameOfGroup].T;
 693                 }
 694                 else{
 695                         if((_nbTimeStep-1)%_freqSave ==0)
 696                                 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Neumann boundary condition for concentration, velocity and temperature"<<endl;
 697                         concentration=_Vj[0];
 698                         Tm=_Vj[_nVar-1];
 699                 }
 700
 701                 double pression=_limitField[nameOfGroup].p + hydroPress;
 702                 double rho_v=_fluides[0]->getDensity(pression, Tm);
 703                 double rho_l=_fluides[1]->getDensity(pression, Tm);
 704                 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
 705                 {
 706                         cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
 707                         cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 708                         *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 709                         _runLogFile->close();
 710                         throw CdmathException("DriftModel::jacobian: Inlet, impossible to compute mixture density, division by zero");
 711                 }
 712
 713                 _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
 714                 _externalStates[1]=concentration*_externalStates[0];
 715                 double mv=_externalStates[1], ml=_externalStates[0]-_externalStates[1];
 716                 for(k=0; k<_Ndim; k++)
 717                 {
 718                         v2+=_Vj[k+2]*_Vj[k+2];
 719                         _externalStates[k+2]=_externalStates[0]*_Vj[(k+2)] ;
 720                 }
 721                 _externalStates[_nVar-1] = mv*_fluides[0]->getInternalEnergy(Tm,rho_v)+ml*_fluides[1]->getInternalEnergy(Tm,rho_l) +_externalStates[0]* v2/2;
 722
 723                 _idm[0] = 0;
 724                 for(k=1; k<_nVar; k++)
 725                         _idm[k] = _idm[k-1] + 1;
 726                 VecAssemblyBegin(_Uext);
 727                 VecAssemblyBegin(_Uextdiff);
 728                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 729                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 730                 VecAssemblyEnd(_Uext);
 731                 VecAssemblyEnd(_Uextdiff);
 732         }
 733         else if (_limitField[nameOfGroup].bcType==Outlet){
 734                 if(q_n<=0  &&  (_nbTimeStep-1)%_freqSave ==0)
 735                         cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition for concentration, velocity and temperature"<<endl;
 736
 737                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 738                 Cell Cj=_mesh.getCell(j);
 739                 double hydroPress=(Cj.x()-_gravityReferencePoint[0])*_GravityField3d[0];
 740                 if(_Ndim>1){
 741                         hydroPress+=(Cj.y()-_gravityReferencePoint[1])*_GravityField3d[1];
 742                         if(_Ndim>2)
 743                                 hydroPress+=(Cj.z()-_gravityReferencePoint[2])*_GravityField3d[2];
 744                 }
 745                 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho
 746
 747                 //Building the external state
 748                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 749
 750                 double concentration=_Vj[0];
 751                 double pression=_limitField[nameOfGroup].p+hydroPress;
 752                 double Tm=_Vj[_nVar-1];
 753                 double rho_v=_fluides[0]->getDensity(pression, Tm);
 754                 double rho_l=_fluides[1]->getDensity(pression, Tm);
 755                 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
 756                 {
 757                         cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
 758                         cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 759                         *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
 760                         _runLogFile->close();
 761                         throw CdmathException("DriftModel::jacobian: Inlet, impossible to compute mixture density, division by zero");
 762                 }
 763
 764                 _externalStates[0]=porosityj*rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
 765                 _externalStates[1]=concentration*_externalStates[0];
 766                 double mv=_externalStates[1], ml=_externalStates[0]-_externalStates[1];
 767                 for(k=0; k<_Ndim; k++)
 768                 {
 769                         v2+=_Vj[k+2]*_Vj[k+2];
 770                         _externalStates[k+2]=_externalStates[0]*_Vj[(k+2)] ;
 771                 }
 772                 _externalStates[_nVar-1] = mv*_fluides[0]->getInternalEnergy(Tm,rho_v)+ml*_fluides[1]->getInternalEnergy(Tm,rho_l) +_externalStates[0]* v2/2;
 773
 774                 _idm[0] = 0;
 775                 for(k=1; k<_nVar; k++)
 776                         _idm[k] = _idm[k-1] + 1;
 777                 VecAssemblyBegin(_Uext);
 778                 VecAssemblyBegin(_Uextdiff);
 779                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 780                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 781                 VecAssemblyEnd(_Uext);
 782                 VecAssemblyEnd(_Uextdiff);
 783         }
 784         else {
 785                 cout<<"!!!!!!!!!!!!!!! Error DriftModel::setBoundaryState !!!!!!!!!!"<<endl;
 786                 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 787                 cout<<"Accepted boundary condition are Neumann, Wall, InnerWall, Inlet, and Outlet"<<endl;
 788                 *_runLogFile<<"Boundary condition not set for boundary named "<<nameOfGroup<<"Accepted boundary condition are Neumann,Wall, InnerWall, Inlet, and Outlet"<<endl;
 789                 _runLogFile->close();
 790                 throw CdmathException("Unknown boundary condition");
 791         }
 792 }
 793
 794 void DriftModel::convectionMatrices()
 795 {
 796         //entree: URoe = rho, cm, u, H
 797         //sortie: matrices Roe+  et Roe- +Roe si scheme centre
 798
 799         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 800                 cout<<"DriftModel::convectionMatrices()"<<endl;
 801
 802         double umn=0, u_2=0; //valeur de u.normale et |u|²
 803         for(int i=0;i<_Ndim;i++)
 804         {
 805                 u_2 += _Uroe[2+i]*_Uroe[2+i];
 806                 umn += _Uroe[2+i]*_vec_normal[i];//vitesse normale
 807         }
 808
 809         vector<complex<double> > vp_dist(3);
 810
 811         if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case where we need no Roe matrix
 812         {
 813                 staggeredVFFCMatricesConservativeVariables(umn);//Computation of classical upwinding matrices
 814                 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//For use in implicit matrix
 815                         staggeredVFFCMatricesPrimitiveVariables(umn);
 816         }
 817         else//In this case we build a Roe matrix
 818         {
 819                 double rhom=_Uroe[0];
 820                 double cm=_Uroe[1];
 821                 double Hm=_Uroe[_nVar-1];
 822                 double hm=Hm-0.5*u_2;
 823                 double m_v=cm*rhom, m_l=rhom-m_v;
 824                 double Tm;
 825                 Vector vitesse(_Ndim);
 826                 for(int idim=0;idim<_Ndim;idim++)
 827                         vitesse[idim]=_Uroe[2+idim];
 828
 829                 if(cm<_precision)//pure liquid
 830                         Tm=_fluides[1]->getTemperatureFromEnthalpy(hm,rhom);
 831                 else if(cm>1-_precision)
 832                         Tm=_fluides[0]->getTemperatureFromEnthalpy(hm,rhom);
 833                 else//Hypothèse de saturation
 834                         Tm=_Tsat;
 835
 836                 double pression= getMixturePressure(cm, rhom, Tm);
 837
 838                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 839                         cout<<"Roe state rhom="<<rhom<<" cm= "<<cm<<" Hm= "<<Hm<<" Tm= "<<Tm<<" pression "<<pression<<endl;
 840
 841                 getMixturePressureDerivatives( m_v, m_l, pression, Tm);//EOS is involved to express pressure jump and sound speed
 842                 if(_kappa*hm+_khi+cm*_ksi<0){
 843                         *_runLogFile<<"Calcul matrice de Roe: vitesse du son complexe"<<endl;
 844                         _runLogFile->close();
 845                         throw CdmathException("Calcul matrice de Roe: vitesse du son complexe");
 846                 }
 847                 double am=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
 848                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 849                         cout<<"DriftModel::convectionMatrices, sound speed am= "<<am<<endl;
 850
 851                 //On remplit la matrice de Roe
 852                 double ecin=0.5*u_2;
 853                 RoeMatrixConservativeVariables( cm, umn, ecin, Hm,vitesse);
 854
 855                 //On remplit les valeurs propres
 856                 vp_dist[0]=umn+am;
 857                 vp_dist[1]=umn-am;
 858                 vp_dist[2]=umn;
 859
 860                 _maxvploc=fabs(umn)+am;
 861                 if(_maxvploc>_maxvp)
 862                         _maxvp=_maxvploc;
 863
 864                 /* Construction des matrices de decentrement */
 865                 if(_spaceScheme== centered){
 866                         if(_entropicCorrection)
 867                         {
 868                                 *_runLogFile<<"DriftModel::convectionMatrices: entropy schemes not yet available for drift model"<<endl;
 869                                 _runLogFile->close();
 870                                 throw CdmathException("DriftModel::convectionMatrices: entropy schemes not yet available for drift model");
 871                         }
 872
 873                         for(int i=0; i<_nVar*_nVar;i++)
 874                                 _absAroe[i] = 0;
 875                         if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)
 876                                 for(int i=0; i<_nVar*_nVar;i++)
 877                                         _absAroeImplicit[i] = 0;
 878                 }
 879                 else if( _spaceScheme ==staggered){
 880                         //Calcul de décentrement de type décalé pour formulation Roe
 881                         staggeredRoeUpwindingMatrixConservativeVariables( cm, umn, ecin, Hm, vitesse);
 882                         //staggeredRoeUpwindingMatrixPrimitiveVariables( cm, umn, ecin, Hm, vitesse);
 883                 }
 884                 else if(_spaceScheme == upwind || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach ){
 885                         if(_entropicCorrection)
 886                                 entropicShift(_vec_normal);
 887                         else
 888                                 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
 889
 890                         vector< complex< double > > y (3,0);
 891                         for( int i=0 ; i<3 ; i++)
 892                                 y[i] = Polynoms::abs_generalise(vp_dist[i])+1*_entropicShift[i];
 893                         Polynoms::abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
 894
 895                         if( _spaceScheme ==pressureCorrection)
 896                         {
 897                                 for( int i=0 ; i<_Ndim ; i++)
 898                                         for( int j=0 ; j<_Ndim ; j++)
 899                                                 _absAroe[(2+i)*_nVar+2+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 900                         }
 901                         else if( _spaceScheme ==lowMach){
 902                                 double M=sqrt(u_2)/am;
 903                                 for( int i=0 ; i<_Ndim ; i++)
 904                                         for( int j=0 ; j<_Ndim ; j++)
 905                                                 _absAroe[(2+i)*_nVar+2+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 906                         }
 907                 }
 908                 else
 909                 {
 910                         *_runLogFile<<"DriftModel::convectionMatrices: scheme not treated"<<endl;
 911                         _runLogFile->close();
 912                         throw CdmathException("DriftModel::convectionMatrices: scheme not treated");
 913                 }
 914
 915                 for(int i=0; i<_nVar*_nVar;i++)
 916                 {
 917                         _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
 918                         _AroePlus[i]  = (_Aroe[i]+_absAroe[i])/2;
 919                 }
 920                 if(_timeScheme==Implicit)
 921                 {
 922                         if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
 923                         {
 924                                 _Vij[0]=cm;
 925                                 _Vij[1]=pression;//pressure
 926                                 _Vij[_nVar-1]=Tm;//Temperature
 927                                 for(int idim=0;idim<_Ndim; idim++)
 928                                         _Vij[2+idim]=_Uroe[2+idim];
 929                                 primToConsJacobianMatrix(_Vij);
 930                                 Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
 931                                 Polynoms::matrixProduct(_AroePlus,  _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
 932                         }
 933                         else
 934                                 for(int i=0; i<_nVar*_nVar;i++)
 935                                 {
 936                                         _AroeMinusImplicit[i] = _AroeMinus[i];
 937                                         _AroePlusImplicit[i]  = _AroePlus[i];
 938                                 }
 939                 }
 940                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
 941                 {
 942                         displayMatrix(_Aroe, _nVar,"Matrice de Roe");
 943                         displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
 944                         displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
 945                         displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
 946                 }
 947         }
 948
 949         if(_verbose && (_nbTimeStep-1)%_freqSave ==0 && _timeScheme==Implicit)
 950         {
 951                 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
 952                 displayMatrix(_AroePlusImplicit, _nVar,"Matrice _AroePlusImplicit");
 953         }
 954
 955         /*******Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source***/
 956         if(_entropicCorrection)
 957         {
 958                 InvMatriceRoe( vp_dist);
 959                 Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
 960         }
 961         else if (_spaceScheme==upwind  || (_spaceScheme ==lowMach) || (_spaceScheme ==pressureCorrection))//upwind sans entropic
 962                 SigneMatriceRoe( vp_dist);
 963         else if(_spaceScheme== centered)//centre  sans entropic
 964                 for(int i=0; i<_nVar*_nVar;i++)
 965                         _signAroe[i] = 0;
 966         else if(_spaceScheme ==staggered )
 967         {
 968                 double signu;
 969                 if(umn>0)
 970                         signu=1;
 971                 else if (umn<0)
 972                         signu=-1;
 973                 else
 974                         signu=0;
 975                 for(int i=0; i<_nVar*_nVar;i++)
 976                         _signAroe[i] = 0;
 977                 _signAroe[0] = signu;
 978                 _signAroe[1+_nVar] = signu;
 979                 for(int i=2; i<_nVar-1;i++)
 980                         _signAroe[i*_nVar+i] = -signu;
 981                 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
 982         }
 983         else
 984         {
 985                 *_runLogFile<<"DriftModel::convectionMatrices: well balanced option not treated"<<endl;
 986                 _runLogFile->close();
 987                 throw CdmathException("DriftModel::convectionMatrices: well balanced option not treated");
 988         }
 989
 990         if(_verbose && (_nbTimeStep-1)%_freqSave ==0 && _timeScheme==Implicit)
 991                 displayMatrix(_signAroe, _nVar,"Signe de la matrice de Roe");
 992 }
 993
 994 void DriftModel::addDiffusionToSecondMember
 995 (               const int &i,
 996                 const int &j,
 997                 bool isBord)
 998 {
 999         double Tm=_Udiff[_nVar-1];
1000         double lambdal=_fluides[1]->getConductivity(Tm);
1001         double lambdav = _fluides[0]->getConductivity(Tm);
1002         double mu_l = _fluides[1]->getViscosity(Tm);
1003         double mu_v = _fluides[0]->getViscosity(Tm);
1004
1005         if(mu_v==0 && mu_l ==0 && lambdav==0 && lambdal==0 && _heatTransfertCoeff==0)
1006                 return;
1007
1008         //extraction des valeurs
1009         for(int k=0; k<_nVar; k++)
1010                 _idm[k] = _nVar*i + k;
1011
1012         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1013         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1014         {
1015                 cout << "Contribution diffusion: variables primitives maille " << i<<endl;
1016                 for(int q=0; q<_nVar; q++)
1017                         cout << _Vi[q] << endl;
1018                 cout << endl;
1019         }
1020
1021         if(!isBord ){
1022                 for(int k=0; k<_nVar; k++)
1023                         _idn[k] = _nVar*j + k;
1024
1025                 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1026         }
1027         else
1028         {
1029                 lambdal=max(lambdal,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
1030                 for(int k=0; k<_nVar; k++)
1031                         _idn[k] = k;
1032
1033                 VecGetValues(_Uextdiff, _nVar, _idn, _Uj);
1034                 consToPrim(_Uj,_Vj,1);
1035         }
1036         if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1037         {
1038                 cout << "Contribution diffusion: variables primitives maille " <<j <<endl;
1039                 for(int q=0; q<_nVar; q++)
1040                         cout << _Vj[q] << endl;
1041                 cout << endl;
1042         }
1043         double pression=(_Vi[1]+_Vj[1])/2;//ameliorer car traitement different pour pression et temperature
1044         double m_v=_Udiff[1];
1045         double rho_v=_fluides[0]->getDensity(pression,Tm);
1046         double rho_l=_fluides[1]->getDensity(pression,Tm);
1047         double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
1048         double mu = alpha_v*mu_v+alpha_l*mu_l;
1049         double lambda = alpha_v*lambdav+alpha_l*lambdal;
1050
1051         //pas de diffusion sur la masse totale
1052         _phi[0]=0;
1053         //on n'a pas encore mis la contribution sur la masse
1054         _phi[1]=0;
1055         //contribution visqueuse sur la quantite de mouvement
1056         for(int k=2; k<_nVar-1; k++)
1057                 _phi[k] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu*(_porosityj*_Vj[k] - _porosityi*_Vi[k]);
1058
1059         //contribution visqueuse sur l'energie
1060         _phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda*(_porosityj*_Vj[_nVar-1] - _porosityi*_Vi[_nVar-1]);
1061
1062         _idm[0] = i;
1063         VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1064
1065         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1066         {
1067                 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
1068                 for(int q=0; q<_nVar; q++)
1069                         cout << _phi[q] << endl;
1070                 cout << endl;
1071         }
1072
1073         if(!isBord)
1074         {
1075                 //On change de signe pour l'autre contribution
1076                 for(int k=0; k<_nVar; k++)
1077                         _phi[k] *= -_inv_dxj/_inv_dxi;
1078                 _idn[0] = j;
1079
1080                 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1081
1082                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1083                 {
1084                         cout << "Contribution diffusion au 2nd membre pour la maille  " << j << ": "<<endl;
1085                         for(int q=0; q<_nVar; q++)
1086                                 cout << _phi[q] << endl;
1087                         cout << endl;
1088                 }
1089         }
1090 }
1091
1092
1093 void DriftModel::sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i)
1094 {
1095         double phirho=Ui[0],phim1=Ui[1],phim2=phirho-phim1,phirhoE=Ui[_nVar-1], cv=Vi[0],  P=Vi[1], T=Vi[_nVar-1];
1096         double norm_u=0;
1097         for(int k=0; k<_Ndim; k++)
1098                 norm_u+=Vi[2+k]*Vi[2+k];
1099         norm_u=sqrt(norm_u);
1100         double h=(phirhoE-0.5*phirho*norm_u*norm_u+P*_porosityField(i))/phirho;//e+p/rho
1101
1102         Si[0]=0;
1103         //if(T>_Tsat && cv<1-_precision)
1104         if(_hsatv>h  && h>_hsatl && cv<1-_precision)//heated boiling//
1105                 Si[1]=_heatPowerField(i)/_latentHeat*_porosityField(i);//phi*gamma
1106         else if (P<_Psat && cv<1-_precision)//flash boiling
1107                 Si[1]=-_dHsatl_over_dp*_dp_over_dt(i)/_latentHeat;
1108         else
1109                 Si[1]=0;
1110         for(int k=2; k<_nVar-1; k++)
1111                 Si[k]  =_gravite[k]*phirho-(phim1*_dragCoeffs[0]+phim2*_dragCoeffs[1])*norm_u*Vi[k];
1112
1113         Si[_nVar-1]=_heatPowerField(i);
1114
1115         for(int k=0; k<_Ndim; k++)
1116                 Si[_nVar-1] +=(_GravityField3d[k]*phirho-(phim1*_dragCoeffs[0]+phim2*_dragCoeffs[1])*norm_u*Vi[2+k])*Vi[2+k];
1117
1118         if(_timeScheme==Implicit)
1119         {
1120                 for(int k=0; k<_nVar*_nVar;k++)
1121                         _GravityImplicitationMatrix[k] = 0;
1122                 if(!_usePrimitiveVarsInNewton)
1123                         for(int k=0; k<_nVar;k++)
1124                                 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
1125                 else
1126                 {
1127                         getDensityDerivatives( cv, P, T ,norm_u*norm_u);
1128                         for(int k=0; k<_nVar;k++)
1129                         {
1130                                 _GravityImplicitationMatrix[k*_nVar+0]      =-_gravite[k]*_drho_sur_dcv;
1131                                 _GravityImplicitationMatrix[k*_nVar+1]      =-_gravite[k]*_drho_sur_dp;
1132                                 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
1133                         }
1134                 }
1135         }
1136
1137         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1138         {
1139                 cout<<"DriftModel::sourceVector"<<endl;
1140                 cout<<"Ui="<<endl;
1141                 for(int k=0;k<_nVar;k++)
1142                         cout<<Ui[k]<<", ";
1143                 cout<<endl;
1144                 cout<<"Vi="<<endl;
1145                 for(int k=0;k<_nVar;k++)
1146                         cout<<Vi[k]<<", ";
1147                 cout<<endl;
1148                 cout<<"Si="<<endl;
1149                 for(int k=0;k<_nVar;k++)
1150                         cout<<Si[k]<<", ";
1151                 cout<<endl;
1152                 if(_timeScheme==Implicit)
1153                         displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
1154         }
1155 }
1156
1157 void DriftModel::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
1158 {
1159         double norm_u=0, u_n=0, phirho;
1160         for(int i=0;i<_Ndim;i++)
1161                 u_n += _Uroe[2+i]*_vec_normal[i];
1162
1163         pressureLoss[0]=0;
1164         pressureLoss[1]=0;
1165         if(u_n>0){
1166                 for(int i=0;i<_Ndim;i++)
1167                         norm_u += Vi[2+i]*Vi[2+i];
1168                 norm_u=sqrt(norm_u);
1169                 phirho=Ui[0];
1170                 for(int i=0;i<_Ndim;i++)
1171                         pressureLoss[2+i]=-K*phirho*norm_u*Vi[2+i];
1172         }
1173         else{
1174                 for(int i=0;i<_Ndim;i++)
1175                         norm_u += Vj[2+i]*Vj[2+i];
1176                 norm_u=sqrt(norm_u);
1177                 phirho=Uj[0];
1178                 for(int i=0;i<_Ndim;i++)
1179                         pressureLoss[2+i]=-K*phirho*norm_u*Vj[2+i];
1180         }
1181         pressureLoss[_nVar-1]=-K*phirho*norm_u*norm_u*norm_u;
1182
1183
1184         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1185         {
1186                 cout<<"DriftModel::pressureLossVector K= "<<K<<endl;
1187                 cout<<"pressure loss vector phirho= "<< phirho<<" norm_u= "<<norm_u<<endl;
1188                 cout<<"Ui="<<endl;
1189                 for(int k=0;k<_nVar;k++)
1190                         cout<<Ui[k]<<", ";
1191                 cout<<endl;
1192                 cout<<"Vi="<<endl;
1193                 for(int k=0;k<_nVar;k++)
1194                         cout<<Vi[k]<<", ";
1195                 cout<<endl;
1196                 cout<<"Uj="<<endl;
1197                 for(int k=0;k<_nVar;k++)
1198                         cout<<Uj[k]<<", ";
1199                 cout<<endl;
1200                 cout<<"Vj="<<endl;
1201                 for(int k=0;k<_nVar;k++)
1202                         cout<<Vj[k]<<", ";
1203                 cout<<endl;
1204                 cout<<"pressureLoss="<<endl;
1205                 for(int k=0;k<_nVar;k++)
1206                         cout<<pressureLoss[k]<<", ";
1207                 cout<<endl;
1208         }
1209 }
1210 void DriftModel::porosityGradientSourceVector()
1211 {
1212         double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[1], pj=_Vj[1],pij;
1213         for(int i=0;i<_Ndim;i++) {
1214                 u_ni += _Vi[2+i]*_vec_normal[i];
1215                 u_nj += _Vj[2+i]*_vec_normal[i];
1216         }
1217         _porosityGradientSourceVector[0]=0;
1218         _porosityGradientSourceVector[1]=0;
1219         rhoj=_Uj[0]/_porosityj;
1220         rhoi=_Ui[0]/_porosityi;
1221         pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
1222         for(int i=0;i<_Ndim;i++)
1223                 _porosityGradientSourceVector[2+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1224         _porosityGradientSourceVector[_nVar-1]=0;
1225 }
1226
1227 void DriftModel::jacobian(const int &j, string nameOfGroup,double * normale)
1228 {
1229         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1230                 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
1231
1232         int k;
1233         for(k=0; k<_nVar*_nVar;k++)
1234                 _Jcb[k] = 0;//No implicitation at this stage
1235
1236         _idm[0] = _nVar*j;
1237         for(k=1; k<_nVar; k++)
1238                 _idm[k] = _idm[k-1] + 1;
1239         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
1240         double q_n=0;//quantité de mouvement normale à la paroi
1241         for(k=0; k<_Ndim; k++)
1242                 q_n+=_externalStates[(k+1)]*normale[k];
1243
1244         // loop of boundary types
1245         if (_limitField[nameOfGroup].bcType==Wall || _limitField[nameOfGroup].bcType==InnerWall)
1246         {
1247                 for(k=0; k<_nVar;k++)
1248                         _Jcb[k*_nVar + k] = 1;
1249                 for(k=2; k<_nVar-1;k++)
1250                         for(int l=2; l<_nVar-1;l++)
1251                                 _Jcb[k*_nVar + l] -= 2*normale[k-2]*normale[l-2];
1252         }
1253         else if (_limitField[nameOfGroup].bcType==Inlet)
1254         {
1255                 return;//For the moment no implicitation, should debug inlet
1256                 if(q_n<0){
1257                         //Boundary state quantities
1258                         double v[_Ndim], ve[_Ndim], v2, ve2;
1259                         VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1260                         double concentration=_limitField[nameOfGroup].conc;
1261                         double pression=_Vj[1];
1262                         double T=_limitField[nameOfGroup].T;
1263                         double rho_v=_fluides[0]->getDensity(pression,T);
1264                         double rho_l=_fluides[1]->getDensity(pression,T);
1265                         if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1266                         {
1267                                 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1268                                 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1269                                 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1270                                 _runLogFile->close();
1271                                 throw CdmathException("DriftModel::jacobian: Inlet, impossible to compute mixture density, division by zero");
1272                         }
1273
1274                         double rho_e=rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
1275                         double e_v=_fluides[0]->getInternalEnergy(T,rho_v);
1276                         double e_l=_fluides[1]->getInternalEnergy(T,rho_l);
1277                         ve[0] = _limitField[nameOfGroup].v_x[0];
1278                         v[0]=_Vj[1];
1279                         ve2 = ve[0]*ve[0];
1280                         v2 = v[0]*v[0];
1281                         if (_Ndim >1){
1282                                 ve[1] = _limitField[nameOfGroup].v_y[0];
1283                                 v[1]=_Vj[2];
1284                                 ve2 += ve[1]*ve[1];
1285                                 v2 = v[1]*v[1];
1286                         }
1287                         if (_Ndim >2){
1288                                 ve[2] = _limitField[nameOfGroup].v_z[0];
1289                                 v[2]=_Vj[3];
1290                                 ve2 += ve[2]*ve[2];
1291                                 v2 = v[2]*v[2];
1292                         }
1293                         double E_e = concentration*e_v+(1-concentration)*e_l + ve2/2;
1294
1295                         //Pressure differential
1296                         double gamma_v =_fluides[0]->constante("gamma");
1297                         double gamma_l =_fluides[1]->constante("gamma");
1298                         double omega=concentration/((gamma_v-1)*rho_v*rho_v*e_v)+(1-concentration)/((gamma_l-1)*rho_l*rho_l*e_l);
1299                         double rhom=_externalStates[0];
1300                         double m_v=concentration*rhom, m_l=rhom-m_v;
1301                         getMixturePressureDerivatives( m_v, m_l, pression, T);
1302
1303                         double CoeffCol1 = rho_e*rho_e*omega*(_khi+_kappa*v2/2);
1304                         double CoeffCol2 = rho_e*rho_e*omega*_ksi;
1305                         double CoeffCol3et4 = rho_e*rho_e*omega*_kappa;
1306
1307                         //Mass line
1308                         _Jcb[0]=CoeffCol1;
1309                         _Jcb[1]=CoeffCol2;
1310                         for(k=0;k<_Ndim;k++)
1311                                 _Jcb[2+k]=-CoeffCol3et4*v[k];
1312                         _Jcb[_nVar-1]=CoeffCol3et4;
1313                         //vapour mass line
1314                         for(k=0;k<_nVar;k++)
1315                                 _Jcb[_nVar+k]=_Jcb[k]*concentration;
1316                         //Momentum lines
1317                         for(int l=0;l<_Ndim;l++)
1318                                 for(k=0;k<_nVar;k++)
1319                                         _Jcb[(2+l)*_nVar+k]=_Jcb[k]*ve[l];
1320                         //Energy line
1321                         for(k=0;k<_nVar;k++)
1322                                 _Jcb[(_nVar-1)*_nVar+k]=_Jcb[k]*E_e;
1323                 }
1324                 else
1325                         for(k=0;k<_nVar;k++)
1326                                 _Jcb[k*_nVar+k]=1;
1327         }
1328         else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0)
1329         {
1330                 return;//For the moment no implicitation, should debug inletPressure
1331                 //Boundary state quantities
1332                 double v[_Ndim], v2=0;
1333                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1334                 double concentration=_limitField[nameOfGroup].conc;
1335                 double pression=_limitField[nameOfGroup].p;
1336                 double T=_limitField[nameOfGroup].T;
1337                 double rho_v=_fluides[0]->getDensity(pression,T);
1338                 double rho_l=_fluides[1]->getDensity(pression,T);
1339                 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1340                 {
1341                         cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1342                         cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1343                         *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1344                         _runLogFile->close();
1345                         throw CdmathException("DriftModel::jacobian: InletPressure, impossible to compute mixture density, division by zero");
1346                 }
1347
1348                 double rho_ext=rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
1349                 double rho_int=_externalStates[0];
1350                 double ratio_densite=rho_ext/rho_int;
1351                 for(k=0;k<_Ndim;k++){
1352                         v[k]=_Vj[2+k];
1353                         v2+=v[k]*v[k];
1354                 }
1355                 //Momentum lines
1356                 for(int l=0;l<_Ndim;l++){
1357                         _Jcb[(2+l)*_nVar]=-ratio_densite*v[l];
1358                         _Jcb[(2+l)*_nVar+2+l]=ratio_densite;
1359                 }
1360                 //Energy line
1361                 _Jcb[(_nVar-1)*_nVar]=-ratio_densite*v2;
1362                 for(int l=0;l<_Ndim;l++)
1363                         _Jcb[(_nVar-1)*_nVar+2+l]=ratio_densite*v[l];
1364
1365         }
1366         else if (_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0)){
1367                 return;//For the moment no implicitation, should debug inletPressure
1368                 //Boundary state quantities
1369                 double v[_Ndim], v2;
1370                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1371                 double concentration=_Vj[0];
1372                 double pression=_limitField[nameOfGroup].p;
1373                 double T=_Vj[_nVar-1];
1374                 double rho_v=_fluides[0]->getDensity(pression,T);
1375                 double rho_l=_fluides[1]->getDensity(pression,T);
1376                 if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1377                 {
1378                         cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1379                         cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1380                         *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1381                         _runLogFile->close();
1382                         throw CdmathException("DriftModel::jacobian: Outlet, impossible to compute mixture density, division by zero");
1383                 }
1384
1385                 double rho_ext=rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v);
1386                 double rho_int=_externalStates[0];
1387                 double e_v=_fluides[0]->getInternalEnergy(T,rho_v);
1388                 double e_l=_fluides[1]->getInternalEnergy(T,rho_l);
1389                 double ratio_densite=rho_ext/rho_int;
1390                 for(k=0;k<_Ndim;k++){
1391                         v[k]=_Vj[2+k];
1392                         v2+=v[k]*v[k];
1393                 }
1394                 double E_m = concentration*e_v+(1-concentration)*e_l + v2/2;//total energy
1395
1396                 //Temperature differential
1397                 double C_v=  _fluides[0]->constante("cv");
1398                 double C_l=      _fluides[1]->constante("cv");
1399                 double ev0=_fluides[0]->getInternalEnergy(0,rho_v);//Corriger
1400                 double el0=_fluides[1]->getInternalEnergy(0,rho_l);//Corriger
1401
1402                 double omega=concentration*C_v/(rho_v*e_v)+(1-concentration)*C_l/(rho_l*e_l);
1403                 double rhomem=_externalStates[0]*(concentration*e_v+(1-concentration)*e_l);
1404                 double m_v=concentration*rho_ext, m_l=rho_ext-m_v;
1405                 double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
1406                 double denom=m_v *C_v+m_l* C_l;
1407
1408                 double khi=(m_v*(ev0*C_l-el0*C_v)-rhomem*C_l)/(denom*denom);
1409                 double ksi=(rho_ext*(el0*C_v-ev0*C_l)+rhomem*(C_l-C_v))/(denom*denom);
1410                 double kappa=1/denom;
1411
1412                 double CoeffCol1 = rho_int*rho_int*omega*(khi+kappa*v2/2)+ratio_densite;//The +ratio_densite helps filling the lines other than total mass
1413                 double CoeffCol2 = rho_int*rho_int*omega*ksi;
1414                 double CoeffCol3et4 = rho_int*rho_int*omega*kappa;
1415
1416                 //Mass line
1417                 _Jcb[0]=-CoeffCol1;
1418                 _Jcb[1]=-CoeffCol2;
1419                 for(k=0;k<_Ndim;k++)
1420                         _Jcb[2+k]=CoeffCol3et4*v[k];
1421                 _Jcb[_nVar-1]=-CoeffCol3et4;
1422                 //vapour mass line
1423                 for(k=0;k<_nVar;k++)
1424                         _Jcb[_nVar+k]=_Jcb[k]*concentration;
1425                 //Momentum lines
1426                 for(int l=0;l<_Ndim;l++)
1427                         for(k=0;k<_nVar;k++)
1428                                 _Jcb[(2+l)*_nVar+k]=_Jcb[k]*v[l];
1429                 //Energy line
1430                 for(k=0;k<_nVar;k++)
1431                         _Jcb[(_nVar-1)*_nVar+k]=_Jcb[k]*E_m;
1432
1433                 //adding the remaining diagonal term
1434                 for(k=0;k<_nVar;k++)
1435                         _Jcb[k*_nVar+k]+=ratio_densite;
1436
1437         }
1438         else  if (_limitField[nameOfGroup].bcType!=Neumann)
1439         {
1440                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1441                 _runLogFile->close();
1442                 throw CdmathException("DriftModel::jacobian: The boundary condition is not recognised: neither inlet, inltPressure, outlet, wall, InnerWall, nor neumann");
1443         }
1444 }
1445
1446 //calcule la jacobienne pour une CL de diffusion
1447 void  DriftModel::jacobianDiff(const int &j, string nameOfGroup)
1448 {
1449         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1450                 cout<<"Jacobienne condition limite diffusion bord "<< nameOfGroup<<endl;
1451
1452
1453         double v2=0,cd = 1,cn=0,p0, gamma;
1454         int idim,jdim,k;
1455
1456         for(k=0; k<_nVar*_nVar;k++)
1457                 _JcbDiff[k] = 0;
1458
1459         if (_limitField[nameOfGroup].bcType==Wall || _limitField[nameOfGroup].bcType==InnerWall){
1460         }
1461         else if (_limitField[nameOfGroup].bcType==Outlet || _limitField[nameOfGroup].bcType==Neumann
1462                         ||_limitField[nameOfGroup].bcType==Inlet || _limitField[nameOfGroup].bcType==InletPressure)
1463         {
1464                 for(k=0;k<_nVar;k++)
1465                         _JcbDiff[k*_nVar+k]=1;
1466         }
1467         else{
1468                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1469                 _runLogFile->close();
1470                 throw CdmathException("DriftModel::jacobianDiff: This boundary condition is not recognised");
1471         }
1472 }
1473
1474
1475 void DriftModel::computeScaling(double maxvp)
1476 {
1477         //      if(_spaceScheme!=staggered)
1478         {
1479                 _blockDiag[0]=1;
1480                 _invBlockDiag[0]=1;
1481                 _blockDiag[1]=1;
1482                 _invBlockDiag[1]=1;
1483                 for(int q=2; q<_nVar-1; q++)
1484                 {
1485                         _blockDiag[q]=1./maxvp;//
1486                         _invBlockDiag[q]= maxvp;//1.;//
1487                 }
1488                 _blockDiag[_nVar - 1]=1/(maxvp*maxvp);//1
1489                 _invBlockDiag[_nVar - 1]=  1./_blockDiag[_nVar - 1] ;// 1.;//
1490         }
1491         /*
1492         else{//_spaceScheme==staggered
1493                 _blockDiag[0]=maxvp*maxvp;
1494                 _invBlockDiag[0]=1/_blockDiag[0];
1495                 _blockDiag[1]=maxvp*maxvp;
1496                 _invBlockDiag[1]=1/_blockDiag[1];
1497                 for(int q=2; q<_nVar-1; q++)
1498                 {
1499                         _blockDiag[q]=1;//
1500                         _invBlockDiag[q]= 1;//1.;//
1501                 }
1502                 _blockDiag[_nVar - 1]=1;//1
1503                 _invBlockDiag[_nVar - 1]=  1./_blockDiag[_nVar - 1] ;// 1.;//
1504         }
1505          */
1506 }
1507 Vector DriftModel::computeExtendedPrimState(double *V)
1508 {
1509         Vector Vext(7+2*_Ndim);
1510
1511         double C=V[0], P=V[1], T=V[_nVar-1];
1512         double rho_v=_fluides[0]->getDensity(P,T);
1513         double rho_l=_fluides[1]->getDensity(P,T);
1514         double e_v=_fluides[0]->getInternalEnergy(T,rho_v);
1515         double e_l=_fluides[1]->getInternalEnergy(T,rho_l);
1516         if(fabs(rho_l*C+rho_v*(1-C))<_precision)
1517         {
1518                 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<", concentration= "<<C<<endl;
1519                 *_runLogFile<<"DriftModel::computeExtendedPrimState: impossible to compute void fraction, division by zero"<<endl;
1520                 _runLogFile->close();
1521                 throw CdmathException("DriftModel::computeExtendedPrimState: impossible to compute void fraction, division by zero");
1522         }
1523
1524         _externalStates[0]=rho_v*rho_l/(C*rho_l+(1-C)*rho_v);
1525         double alpha_v=rho_l*C/(rho_l*C+rho_v*(1-C)), alpha_l=1-alpha_v;
1526         double rho_m=alpha_v*rho_v+alpha_l*rho_l;
1527         double h_m=(alpha_v*rho_v*e_v+alpha_l*rho_l*e_l+P)/rho_m;
1528         Vector Vit(_Ndim);
1529         for(int i=0;i<_Ndim;i++)
1530                 Vit(i)=V[2+i];
1531         Vector u_r=relative_velocity(C, Vit,rho_m);
1532
1533         Vext(0)=C;
1534         Vext(1)=P;
1535         for(int i=0;i<_Ndim;i++)
1536                 Vext(2+i)=Vit(i);
1537         Vext(2+_Ndim)=T;
1538         Vext(3+_Ndim)=alpha_v;
1539         for(int i=0;i<_Ndim;i++)
1540                 Vext(4+_Ndim+i)=u_r(i);
1541         Vext((4+2*_Ndim))=rho_v;
1542         Vext((5+2*_Ndim))=rho_l;
1543         Vext(6+2*_Ndim)=h_m;
1544
1545         return Vext;
1546 }
1547
1548
1549 void DriftModel::primToCons(const double *P, const int &i, double *W, const int &j){
1550         double concentration=P[i*_nVar];
1551         double pression=P[i*_nVar+1];
1552         double temperature=P[i*_nVar+_nVar-1];
1553         double rho_v=_fluides[0]->getDensity(pression,temperature);
1554         double rho_l=_fluides[1]->getDensity(pression,temperature);
1555         double e_v = _fluides[0]->getInternalEnergy(temperature,rho_v);
1556         double e_l = _fluides[1]->getInternalEnergy(temperature,rho_l);
1557         if(fabs(concentration*rho_l+(1-concentration)*rho_v)<_precision)
1558         {
1559                 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<endl;
1560                 cout<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1561                 *_runLogFile<<"concentration*rho_l+(1-concentration)*rho_v= "<<concentration*rho_l+(1-concentration)*rho_v<<endl;
1562                 *_runLogFile<<"DriftModel::primToCons: impossible to compute mixture density, division by zero"<<endl;
1563                 _runLogFile->close();
1564                 throw CdmathException("DriftModel::primToCons: impossible to compute mixture density, division by zero");
1565         }
1566         W[j*(_nVar)]  =(rho_v*rho_l/(concentration*rho_l+(1-concentration)*rho_v))*_porosityField(j);//phi*rho
1567         W[j*(_nVar)+1]  =W[j*(_nVar)] *concentration;//phi *rho*c_v
1568         for(int k=0; k<_Ndim; k++)
1569                 W[j*_nVar+(k+2)] = W[j*(_nVar)] *P[i*_nVar+(k+2)];//phi*rho*u
1570
1571         W[j*_nVar+_nVar-1] = W[j*(_nVar)+1]* e_v + (W[j*(_nVar)]-W[j*(_nVar)+1])* e_l;//phi rhom e_m
1572         for(int k=0; k<_Ndim; k++)
1573                 W[j*_nVar+_nVar-1] += W[j*_nVar]*P[i*_nVar+(k+2)]*P[i*_nVar+(k+2)]*0.5;//phi rhom e_m + phi rho u*u
1574
1575         if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
1576                 cout<<"DriftModel::primToCons: T="<<temperature<<", pression= "<<pression<<endl;
1577                 cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<", rhom= "<<W[j*(_nVar)]<<" e_v="<<e_v<<" e_l="<<e_l<<endl;
1578         }
1579 }
1580 void DriftModel::primToConsJacobianMatrix(double *V)
1581 {
1582         double concentration=V[0];
1583         double pression=V[1];
1584         double temperature=V[_nVar-1];
1585         double vitesse[_Ndim];
1586         for(int idim=0;idim<_Ndim;idim++)
1587                 vitesse[idim]=V[2+idim];
1588         double v2=0;
1589         for(int idim=0;idim<_Ndim;idim++)
1590                 v2+=vitesse[idim]*vitesse[idim];
1591
1592         double rho_v=_fluides[0]->getDensity(pression,temperature);
1593         double rho_l=_fluides[1]->getDensity(pression,temperature);
1594         double gamma_v=_fluides[0]->constante("gamma");
1595         double gamma_l=_fluides[1]->constante("gamma");
1596         double q_v=_fluides[0]->constante("q");
1597         double q_l=_fluides[1]->constante("q");
1598
1599         double rho=concentration*rho_v+(1-concentration)*rho_l;;
1600
1601         for(int k=0;k<_nVar*_nVar; k++)
1602                 _primToConsJacoMat[k]=0;
1603
1604         if(             !_useDellacherieEOS)
1605         {
1606                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1607                 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1608                 double e_v = fluide0->getInternalEnergy(temperature);
1609                 double e_l = fluide0->getInternalEnergy(temperature);
1610                 double cv_v=fluide0->constante("cv");
1611                 double cv_l=fluide1->constante("cv");
1612                 double e=concentration*e_v+(1-concentration)*e_l;
1613                 double E=e+0.5*v2;
1614
1615                 /******* Masse totale **********/
1616                 _primToConsJacoMat[0]=-rho*rho*(1/rho_v-1/rho_l);
1617                 _primToConsJacoMat[1]=rho*rho*(
1618                                 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1619                                 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1620                 );
1621                 _primToConsJacoMat[_nVar-1]=-rho*rho*(
1622                                 cv_v*   concentration /(rho_v*(e_v-q_v))
1623                                 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
1624                 );
1625
1626                 /******* Masse partielle **********/
1627                 _primToConsJacoMat[_nVar]=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
1628                 _primToConsJacoMat[_nVar+1]=concentration*rho*rho*(
1629                                 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1630                                 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1631                 );
1632                 _primToConsJacoMat[_nVar+_nVar-1]=-concentration*rho*rho*(
1633                                 cv_v*   concentration /(rho_v*(e_v-q_v))
1634                                 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
1635                 );
1636                 /******* Quantité de mouvement **********/
1637                 for(int idim=0;idim<_Ndim;idim++)
1638                 {
1639                         _primToConsJacoMat[2*_nVar+idim*_nVar]=-vitesse[idim]*rho*rho*(1/rho_v-1/rho_l);
1640                         _primToConsJacoMat[2*_nVar+idim*_nVar+1]=vitesse[idim]*rho*rho*(
1641                                         concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1642                                         +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1643                         );
1644                         _primToConsJacoMat[2*_nVar+idim*_nVar+2+idim]=rho;
1645                         _primToConsJacoMat[2*_nVar+idim*_nVar+_nVar-1]=-vitesse[idim]*rho*rho*(
1646                                         cv_v*   concentration /(rho_v*(e_v-q_v))
1647                                         +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
1648                         );
1649                 }
1650                 /******* Energie totale **********/
1651                 _primToConsJacoMat[(_nVar-1)*_nVar]=rho*(e_v-e_l)-E*rho*rho*(1/rho_v-1/rho_l);
1652                 _primToConsJacoMat[(_nVar-1)*_nVar+1]=E*rho*rho*(
1653                                 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
1654                                 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
1655                 );
1656                 for(int idim=0;idim<_Ndim;idim++)
1657                         _primToConsJacoMat[(_nVar-1)*_nVar+2+idim]=rho*vitesse[idim];
1658                 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*(cv_v*concentration + cv_l*(1-concentration))
1659                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         -rho*rho*E*( cv_v*   concentration /(rho_v*(e_v-q_v))
1660                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         +cv_l*(1-concentration)/(rho_l*(e_l-q_l)));
1661         }
1662         else if(_useDellacherieEOS)
1663         {
1664                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1665                 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1666                 double h_v=fluide0->getEnthalpy(temperature);
1667                 double h_l=fluide1->getEnthalpy(temperature);
1668                 double h=concentration*h_v+(1-concentration)*h_l;
1669                 double H=h+0.5*v2;
1670                 double cp_v=fluide0->constante("cp");
1671                 double cp_l=fluide1->constante("cp");
1672
1673                 /******* Masse totale **********/
1674                 _primToConsJacoMat[0]=-rho*rho*(1/rho_v-1/rho_l);
1675                 _primToConsJacoMat[1]=rho*rho*(
1676                                 gamma_v*   concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1677                                 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1678                 );
1679                 _primToConsJacoMat[_nVar-1]=-rho*rho*(
1680                                 cp_v*   concentration /(rho_v*(h_v-q_v))
1681                                 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
1682                 );
1683
1684                 /******* Masse partielle **********/
1685                 _primToConsJacoMat[_nVar]=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
1686                 _primToConsJacoMat[_nVar+1]=concentration*rho*rho*(
1687                                 gamma_v*   concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1688                                 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1689                 );
1690                 _primToConsJacoMat[_nVar+_nVar-1]=-concentration*rho*rho*(
1691                                 cp_v*   concentration /(rho_v*(h_v-q_v))
1692                                 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
1693                 );
1694                 /******* Quantité de mouvement **********/
1695                 for(int idim=0;idim<_Ndim;idim++)
1696                 {
1697                         _primToConsJacoMat[2*_nVar+idim*_nVar]=-vitesse[idim]*rho*rho*(1/rho_v-1/rho_l);
1698                         _primToConsJacoMat[2*_nVar+idim*_nVar+1]=vitesse[idim]*rho*rho*(
1699                                         gamma_v*   concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1700                                         +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1701                         );
1702                         _primToConsJacoMat[2*_nVar+idim*_nVar+2+idim]=rho;
1703                         _primToConsJacoMat[2*_nVar+idim*_nVar+_nVar-1]=-vitesse[idim]*rho*rho*(
1704                                         cp_v*   concentration /(rho_v*(h_v-q_v))
1705                                         +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
1706                         );
1707                 }
1708                 /******* Energie totale **********/
1709                 _primToConsJacoMat[(_nVar-1)*_nVar]=rho*(h_v-h_l)-H*rho*rho*(1/rho_v-1/rho_l);
1710                 _primToConsJacoMat[(_nVar-1)*_nVar+1]=H*rho*rho*(
1711                                 gamma_v*   concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
1712                                 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
1713                 )-1;
1714                 for(int idim=0;idim<_Ndim;idim++)
1715                         _primToConsJacoMat[(_nVar-1)*_nVar+2+idim]=rho*vitesse[idim];
1716                 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*(cp_v*concentration + cp_l*(1-concentration))
1717                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         -rho*rho*H*(cp_v*   concentration /(rho_v*(h_v-q_v))
1718                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         +cp_l*(1-concentration)/(rho_l*(h_l-q_l)));
1719         }
1720         else
1721         {
1722                 *_runLogFile<<"DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie"<<endl;
1723                 _runLogFile->close();
1724                 throw CdmathException("DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
1725         }
1726
1727         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1728         {
1729                 cout<<" DriftModel::primToConsJacobianMatrix" << endl;
1730                 displayVector(_Vi,_nVar," _Vi " );
1731                 cout<<"rho_v= "<<rho_v<<", rho_l= "<<rho_l<<endl;
1732                 displayMatrix(_primToConsJacoMat,_nVar," Jacobienne primToCons: ");
1733         }
1734 }
1735
1736 void DriftModel::consToPrim(const double *Wcons, double* Wprim,double porosity)
1737 {
1738         double m_v=Wcons[1]/porosity;
1739         double m_l=(Wcons[0]-Wcons[1])/porosity;
1740         _minm1=min(m_v,_minm1);
1741         _minm2=min(m_l,_minm2);
1742         if(m_v<-_precision || m_l<-_precision)
1743                 _nbMaillesNeg+=1;
1744         else if(m_v< 0 && m_v>-_precision )
1745                 m_v=0;
1746         else if(m_l< 0 && m_l>-_precision )
1747                 m_l=0;
1748         double concentration=m_v/(m_v+m_l);
1749         double q_2 = 0;
1750         for(int k=0;k<_Ndim;k++)
1751                 q_2 += Wcons[k+2]*Wcons[k+2];
1752         q_2 /= Wcons[0];        //phi rho u²
1753
1754         double rho_m_e_m=(Wcons[_nVar-1] -0.5*q_2)/porosity;
1755         double pression,Temperature;
1756
1757         if(concentration<_precision)
1758         {
1759                 pression=_fluides[1]->getPressure(rho_m_e_m,m_l);
1760                 Temperature=_fluides[1]->getTemperatureFromPressure(pression,m_l);
1761         }
1762         else if(concentration>1-_precision)
1763         {
1764                 pression=_fluides[0]->getPressure(rho_m_e_m,m_v);
1765                 Temperature=_fluides[0]->getTemperatureFromPressure(pression,m_v);
1766         }
1767         else//Hypothèses de saturation
1768         {
1769                 //Première approche
1770                 getMixturePressureAndTemperature(concentration, m_v+m_l, rho_m_e_m, pression, Temperature);
1771                 //cout<<"Temperature= "<<Temperature<<", pression= "<<pression<<endl;
1772                 //throw CdmathException("DriftModel::consToPrim: Apparition de vapeur");
1773
1774                 //Seconde approche : on impose la température directement
1775                 //Temperature=_Tsat;
1776                 //pression=getMixturePressure(concentration, m_v+m_l, Temperature);
1777
1778                 //Troisieme approche : on impose les enthalpies de saturation
1779                 //double rho_m_h_m= m_v*_hsatv + m_l*_hsatl;
1780                 //pression = rho_m_h_m - rho_m_e_m;
1781                 //Temperature = getMixtureTemperature(concentration, m_v+m_l, pression);
1782         }
1783
1784         if (pression<0){
1785                 cout << "pressure = "<< pression << " < 0 " << endl;
1786                 *_runLogFile << "pressure = "<< pression << " < 0 " << endl;
1787                 cout<<"Vecteur conservatif"<<endl;
1788                 *_runLogFile<<"Vecteur conservatif"<<endl;
1789                 for(int k=0;k<_nVar;k++){
1790                         cout<<Wcons[k]<<endl;
1791                         *_runLogFile<<Wcons[k]<<endl;
1792                 }
1793                 _runLogFile->close();
1794                 throw CdmathException("DriftModel::consToPrim: negative pressure");
1795         }
1796
1797         Wprim[0]=concentration;
1798         Wprim[1] =  pression;
1799         for(int k=0;k<_Ndim;k++)
1800                 Wprim[k+2] = Wcons[k+2]/Wcons[0];
1801         Wprim[_nVar-1] =  Temperature;
1802         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1803         {
1804                 cout<<"ConsToPrim Vecteur conservatif"<<endl;
1805                 for(int k=0;k<_nVar;k++)
1806                         cout<<Wcons[k]<<endl;
1807                 cout<<"ConsToPrim Vecteur primitif"<<endl;
1808                 for(int k=0;k<_nVar;k++)
1809                         cout<<Wprim[k]<<endl;
1810         }
1811 }
1812
1813 double DriftModel::getMixturePressure(double c_v, double rhom, double temperature)
1814 {
1815         double gamma_v=_fluides[0]->constante("gamma");
1816         double gamma_l=_fluides[1]->constante("gamma");
1817         double Pinf_v=_fluides[0]->constante("p0");
1818         double Pinf_l=_fluides[1]->constante("p0");
1819         double q_v=_fluides[0]->constante("q");
1820         double q_l=_fluides[1]->constante("q");
1821         double c_l=1-c_v;
1822         double a=1., b, c;
1823
1824         if(     !_useDellacherieEOS)
1825         {
1826                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1827                 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1828                 double e_v=fluide0->getInternalEnergy(temperature);
1829                 double e_l=fluide1->getInternalEnergy(temperature);
1830                 b=gamma_v*Pinf_v+gamma_l*Pinf_l -rhom*(c_l*(gamma_l-1)*(e_l-q_l)+c_v*(gamma_v-1)*(e_v-q_v));
1831                 c=      gamma_v*Pinf_v*gamma_l*Pinf_l
1832                                 -rhom*(c_l*(gamma_l-1)*(e_l-q_l)*gamma_v*Pinf_v + c_v*(gamma_v-1)*(e_v-q_v)*gamma_l*Pinf_l);
1833         }
1834         else if(_useDellacherieEOS)
1835         {
1836                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1837                 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1838                 double h_v=fluide0->getEnthalpy(temperature);
1839                 double h_l=fluide1->getEnthalpy(temperature);
1840                 b=Pinf_v+Pinf_l -rhom*(c_l*(gamma_l-1)*(h_l-q_l)/gamma_l+c_v*(gamma_v-1)*(h_v-q_v)/gamma_v);
1841                 c=Pinf_v*Pinf_l-rhom*(c_l*(gamma_l-1)*(h_l-q_l)*Pinf_v + c_v*(gamma_v-1)*(h_v-q_v)*Pinf_l);
1842         }
1843         else
1844         {
1845                 *_runLogFile<<"DriftModel::getMixturePressure: phases must have the same eos"<<endl;
1846                 _runLogFile->close();
1847                 throw CdmathException("DriftModel::getMixturePressure: phases must have the same eos");
1848         }
1849
1850         double delta= b*b-4*a*c;
1851
1852         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1853                 cout<<"DriftModel::getMixturePressure: a= "<<a<<", b= "<<b<<", c= "<<c<<", delta= "<<delta<<endl;
1854
1855         if(delta<0){
1856                 *_runLogFile<<"DriftModel::getMixturePressure: cannot compute pressure, delta<0"<<endl;
1857                 _runLogFile->close();
1858                 throw CdmathException("DriftModel::getMixturePressure: cannot compute pressure, delta<0");
1859         }
1860         else
1861                 return (-b+sqrt(delta))/(2*a);
1862 }
1863
1864 void DriftModel::getMixturePressureAndTemperature(double c_v, double rhom, double rhom_em, double& pression, double& temperature)
1865 {
1866         double gamma_v=_fluides[0]->constante("gamma");
1867         double gamma_l=_fluides[1]->constante("gamma");
1868         double Pinf_v=_fluides[0]->constante("p0");
1869         double Pinf_l=_fluides[1]->constante("p0");
1870         double q_v=_fluides[0]->constante("q");
1871         double q_l=_fluides[1]->constante("q");
1872         double c_l=1-c_v, m_v=c_v*rhom, m_l=rhom-m_v;
1873         double a, b, c, delta;
1874
1875         if(     !_useDellacherieEOS)
1876         {
1877                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1878                 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1879
1880                 temperature= (rhom_em-m_v*fluide0->getInternalEnergy(0)-m_l*fluide1->getInternalEnergy(0))
1881                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         /(m_v*fluide0->constante("cv")+m_l*fluide1->constante("cv"));
1882
1883                 double e_v=fluide0->getInternalEnergy(temperature);
1884                 double e_l=fluide1->getInternalEnergy(temperature);
1885                 a=1.;
1886                 b=gamma_v*Pinf_v+gamma_l*Pinf_l -rhom*(c_l*(gamma_l-1)*(e_l-q_l)+c_v*(gamma_v-1)*(e_v-q_v));
1887                 c=      gamma_v*Pinf_v*gamma_l*Pinf_l
1888                                 -rhom*(c_l*(gamma_l-1)*(e_l-q_l)*gamma_v*Pinf_v + c_v*(gamma_v-1)*(e_v-q_v)*gamma_l*Pinf_l);
1889
1890                 delta= b*b-4*a*c;
1891
1892                 if(delta<0){
1893                         *_runLogFile<<"DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0"<<endl;
1894                         _runLogFile->close();
1895                         throw CdmathException("DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0");
1896                 }
1897                 else
1898                         pression= (-b+sqrt(delta))/(2*a);
1899
1900         }
1901         else if(_useDellacherieEOS)
1902         {
1903                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1904                 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1905
1906                 double h0_v=fluide0->getEnthalpy(0);
1907                 double h0_l=fluide1->getEnthalpy(0);
1908                 double cp_v = _fluides[0]->constante("cp");
1909                 double cp_l = _fluides[1]->constante("cp");
1910                 double denom=m_v*cp_v+m_l*cp_l;
1911                 double num_v=cp_v*rhom_em+m_l*(cp_l* h0_v-cp_v* h0_l);
1912                 double num_l=cp_l*rhom_em+m_v*(cp_v* h0_l-cp_l* h0_v);
1913
1914                 a=gamma_v*gamma_l-(m_v*(gamma_v-1)*gamma_l*cp_v+m_l*(gamma_l-1)*gamma_v*cp_l)/denom;
1915                 b=gamma_v*gamma_l*(Pinf_v+Pinf_l)-(m_v*(gamma_v-1)*gamma_l*cp_v*Pinf_l+m_l*(gamma_l-1)*gamma_v*cp_l*Pinf_v)/denom
1916                                 -m_v*(gamma_v-1)*gamma_l*(num_v/denom -q_v)-m_l*(gamma_l-1)*gamma_v*(num_l/denom -q_l);
1917                 c=gamma_v*gamma_l*Pinf_v*Pinf_l
1918                                 -m_v*(gamma_v-1)*gamma_l*(num_v/denom -q_v)*Pinf_l-m_l*(gamma_l-1)*gamma_v*(num_l/denom -q_l)*Pinf_v;
1919
1920                 delta= b*b-4*a*c;
1921
1922                 if(delta<0)
1923                 {
1924                         *_runLogFile<<"DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0"<<endl;
1925                         _runLogFile->close();
1926                         throw CdmathException("DriftModel::getMixturePressureAndTemperature: cannot compute pressure, delta<0");
1927                 }
1928                 else
1929                         pression= (-b+sqrt(delta))/(2*a);
1930
1931                 temperature=(rhom_em+pression-m_v*h0_v-m_l*h0_l)/denom;
1932         }
1933         else
1934         {
1935                 _runLogFile->close();
1936                 throw CdmathException("DriftModel::getMixturePressureAndTemperature: phases must have the same eos");
1937         }
1938
1939
1940         if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
1941                 cout<<"DriftModel::getMixturePressureAndTemperature: a= "<<a<<", b= "<<b<<", c= "<<c<<", delta= "<<delta<<endl;
1942                 cout<<"pressure= "<<pression<<", temperature= "<<temperature<<endl;
1943         }
1944
1945 }
1946 double DriftModel::getMixtureTemperature(double c_v, double rhom, double pression)
1947 {
1948         double gamma_v=_fluides[0]->constante("gamma");
1949         double gamma_l=_fluides[1]->constante("gamma");
1950         double Pinf_v = _fluides[0]->constante("p0");
1951         double Pinf_l = _fluides[1]->constante("p0");
1952         double q_v=_fluides[0]->constante("q");
1953         double q_l=_fluides[1]->constante("q");
1954         double c_l=1-c_v;
1955
1956         if(     !_useDellacherieEOS)
1957         {
1958                 double cv_v = _fluides[0]->constante("cv");
1959                 double cv_l = _fluides[1]->constante("cv");
1960                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1961                 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
1962                 double e0_v=fluide0->getInternalEnergy(0);
1963                 double e0_l=fluide1->getInternalEnergy(0);
1964
1965                 double numerator=(pression+gamma_v*Pinf_v)*(pression+gamma_l*Pinf_l)/rhom
1966                                 - c_l*(pression+gamma_v*Pinf_v)*(gamma_l-1)*(e0_l-q_v)
1967                                 - c_v*(pression+gamma_l*Pinf_l)*(gamma_v-1)*(e0_v-q_l);
1968                 double denominator=  c_l*(pression+gamma_v*Pinf_v)*(gamma_l-1)*cv_l
1969                                 +c_v*(pression+gamma_l*Pinf_l)*(gamma_v-1)*cv_v;
1970                 return numerator/denominator;
1971         }
1972         else if(_useDellacherieEOS)
1973         {
1974                 double cp_v = _fluides[0]->constante("cp");
1975                 double cp_l = _fluides[1]->constante("cp");
1976                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1977                 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
1978                 double h0_v=fluide0->getEnthalpy(0);
1979                 double h0_l=fluide1->getEnthalpy(0);
1980
1981                 double numerator= gamma_v*(pression+Pinf_v)* gamma_l*(pression+Pinf_l)/rhom
1982                                 - c_l*gamma_v*(pression+Pinf_v)*(gamma_l-1)*(h0_l-q_l)
1983                                 - c_v*gamma_l*(pression+Pinf_l)*(gamma_v-1)*(h0_v-q_v);
1984                 double denominator=  c_l*gamma_v*(pression+Pinf_v)*(gamma_l-1)*cp_l
1985                                 +c_v*gamma_l*(pression+Pinf_l)*(gamma_v-1)*cp_v;
1986                 return numerator/denominator;
1987         }
1988         else
1989         {
1990                 *_runLogFile<<"DriftModel::getMixtureTemperature: phases must have the same eos"<<endl;
1991                 _runLogFile->close();
1992                 throw CdmathException("DriftModel::getMixtureTemperature: phases must have the same eos");
1993         }
1994
1995 }
1996
1997 void DriftModel::getMixturePressureDerivatives(double m_v, double m_l, double pression, double Tm)
1998 {
1999         double rho_v=_fluides[0]->getDensity(pression,Tm);
2000         double rho_l=_fluides[1]->getDensity(pression,Tm);
2001         double alpha_v=m_v/rho_v,alpha_l=1-alpha_v;
2002         double gamma_v=_fluides[0]->constante("gamma");
2003         double gamma_l=_fluides[1]->constante("gamma");
2004         double q_v=_fluides[0]->constante("q");
2005         double q_l=_fluides[1]->constante("q");
2006         double temp1, temp2, denom;
2007
2008         if(        !_useDellacherieEOS)
2009         {//Classical stiffened gas with linear law e(T)
2010                 double cv_v = _fluides[0]->constante("cv");
2011                 double cv_l = _fluides[1]->constante("cv");
2012
2013                 double e_v=_fluides[0]->getInternalEnergy(Tm, rho_v);//Actually does not depends on rho_v
2014                 double e_l=_fluides[1]->getInternalEnergy(Tm, rho_l);//Actually does not depends on rho_l
2015
2016                 //On estime les dérivées discrètes de la pression (cf doc)
2017                 temp1= m_v* cv_v + m_l* cv_l;
2018                 denom=(alpha_v/((gamma_v-1)*rho_v*(e_v-q_v))+alpha_l/((gamma_l-1)*rho_l*(e_l-q_l)))*temp1;
2019                 temp2=alpha_v*cv_v/ (e_v-q_v)+alpha_l*cv_l/ (e_l-q_l);
2020
2021                 _khi=(temp1/rho_l-e_l*temp2)/denom;
2022                 _ksi=(temp1*(rho_l-rho_v)/(rho_v*rho_l)+(e_l-e_v)*temp2)/denom;
2023                 _kappa=temp2/denom;
2024         }
2025
2026         else if( _useDellacherieEOS)
2027         {//S. Dellacherie stiffened gas with linear law h(T)
2028                 double cp_v = _fluides[0]->constante("cp");
2029                 double cp_l = _fluides[1]->constante("cp");
2030
2031                 double h_v=_fluides[0]->getEnthalpy(Tm, rho_v);//Actually does not depends on rho_v
2032                 double h_l=_fluides[1]->getEnthalpy(Tm, rho_l);//Actually does not depends on rho_l
2033
2034                 //On estime les dérivées discrètes de la pression (cf doc)
2035                 temp1= m_v* cp_v + m_l* cp_l;
2036                 temp2= alpha_v*cp_v/(h_v-q_v)+alpha_l*cp_l/(h_l-q_l);
2037                 //denom=temp1*(alpha_v/(pression+Pinf_v)+alpha_l/(pression+Pinf_l))-temp2;
2038                 denom=temp1*(alpha_v*gamma_v/((gamma_v-1)*rho_v*(h_v-q_v))+alpha_l*gamma_l/((gamma_l-1)*rho_l*(h_l-q_l)))-temp2;
2039
2040                 //On estime les dérivées discrètes de la pression (cf doc)
2041                 _khi=(temp1/rho_l - h_l*temp2)/denom;
2042                 _ksi=(temp1*(1/rho_v-1/rho_l)+(h_l-h_v)*temp2)/denom;
2043                 _kappa=temp2/denom;
2044         }
2045
2046         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2047         {
2048                 cout<<"rho_l= "<<rho_l<<", temp1= "<<temp1<<", temp2= "<<temp2<<", denom= "<<denom<<endl;
2049                 cout<<"khi= "<<_khi<<", ksi= "<<_ksi<<", kappa= "<<_kappa<<endl;
2050         }
2051 }
2052
2053 void DriftModel::entropicShift(double* n)//TO do: make sure _Vi and _Vj are well set
2054 {
2055         //_Vi is (cm, p, u, T)
2056         //_Ui is (rhom, rhom cm, rhom um, rhom Em)
2057         consToPrim(_Ui,_Vi,_porosityi);
2058         double umnl=0, ul_2=0, umnr=0, ur_2=0; //valeur de u.normale et |u|²
2059         for(int i=0;i<_Ndim;i++)
2060         {
2061                 ul_2 += _Vi[2+i]*_Vi[2+i];
2062                 umnl += _Vi[2+i]*n[i];//vitesse normale
2063                 ur_2 += _Vj[2+i]*_Vj[2+i];
2064                 umnr += _Vj[2+i]*n[i];//vitesse normale
2065         }
2066
2067         //Left sound speed
2068         double rhom=_Ui[0];
2069         double cm=_Vi[0];
2070         double Hm=(_Ui[_nVar-1]+_Vi[1])/rhom;
2071         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2072                 cout<<"Entropic shift left state rhom="<<rhom<<" cm= "<<cm<<"Hm= "<<Hm<<endl;
2073         double Tm=_Vi[_nVar-1];
2074         double hm=Hm-0.5*ul_2;
2075         double m_v=cm*rhom, m_l=rhom-m_v;
2076         double pression=getMixturePressure( cm, rhom, Tm);
2077
2078         getMixturePressureDerivatives( m_v, m_l, pression, Tm);
2079
2080         if(_kappa*hm+_khi+cm*_ksi<0)
2081         {
2082                 *_runLogFile<<"DriftModel::entropicShift : vitesse du son cellule gauche complexe"<<endl;
2083                 _runLogFile->close();
2084                 throw CdmathException("Valeurs propres jacobienne: vitesse du son complexe");
2085         }
2086         double aml=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
2087         //cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", am= "<<am<<endl;
2088
2089         //Right sound speed
2090         consToPrim(_Uj,_Vj,_porosityj);
2091         rhom=_Uj[0];
2092         cm=_Vj[0];
2093         Hm=(_Uj[_nVar-1]+_Vj[1])/rhom;
2094         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2095                 cout<<"Entropic shift right state rhom="<<rhom<<" cm= "<<cm<<"Hm= "<<Hm<<endl;
2096         Tm=_Vj[_nVar-1];
2097         hm=Hm-0.5*ul_2;
2098         m_v=cm*rhom;
2099         m_l=rhom-m_v;
2100         pression=getMixturePressure( cm, rhom, Tm);
2101
2102         getMixturePressureDerivatives( m_v, m_l, pression, Tm);
2103
2104         if(_kappa*hm+_khi+cm*_ksi<0){
2105                 *_runLogFile<<"DriftModel::entropicShift: vitesse du son cellule droite complexe"<<endl;
2106                 _runLogFile->close();
2107                 throw CdmathException("Valeurs propres jacobienne: vitesse du son complexe");
2108         }
2109
2110         double amr=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
2111         //cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", am= "<<am<<endl;
2112
2113         _entropicShift[0]=abs(umnl-aml - (umnr-amr));
2114         _entropicShift[1]=abs(umnl - umnr);
2115         _entropicShift[2]=abs(umnl+aml - (umnr+amr));
2116
2117         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2118         {
2119                 cout<<"Entropic shift left eigenvalues: "<<endl;
2120                 cout<<"("<< umnl-aml<< ", " <<umnl<<", "<<umnl+aml << ")";
2121                 cout<<endl;
2122                 cout<<"Entropic shift right eigenvalues: "<<endl;
2123                 cout<<"("<< umnr-amr<< ", " <<umnr<<", "<<umnr+amr << ")";
2124                 cout<< endl;
2125                 cout<<"eigenvalue jumps "<<endl;
2126                 cout<< _entropicShift[0] << ", " << _entropicShift[1] << ", "<< _entropicShift[2] <<endl;
2127         }
2128 }
2129
2130 Vector DriftModel::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
2131         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2132         {
2133                 cout<<"DriftModel::convectionFlux start"<<endl;
2134                 cout<<"Ucons"<<endl;
2135                 cout<<U<<endl;
2136                 cout<<"Vprim"<<endl;
2137                 cout<<V<<endl;
2138         }
2139
2140         double rho_m=U(0);//phi rhom
2141         double m_v=U(1), m_l=rho_m-m_v;//phi rhom cv, phi rhom cl
2142         Vector q_m(_Ndim);
2143         for(int i=0;i<_Ndim;i++)
2144                 q_m(i)=U(2+i);
2145
2146         double concentration_v=V(0);
2147         double concentration_l= 1 - concentration_v;
2148         double pression=V(1);
2149         Vector vitesse_melange(_Ndim);
2150         for(int i=0;i<_Ndim;i++)
2151                 vitesse_melange(i)=V(2+i);
2152         double Temperature= V(2+_Ndim);//(rho_m_e_m-m_v*internal_energy(0, e0v,c0v,T0)-m_l*internal_energy(0, e0l,c0l,T0))/(m_v*c0v+m_l*c0l);
2153
2154         Vector vitesse_relative=relative_velocity(concentration_v,vitesse_melange,rho_m);
2155         Vector vitesse_v=vitesse_melange+concentration_l*vitesse_relative;
2156         Vector vitesse_l=vitesse_melange-concentration_v*vitesse_relative;
2157         double vitesse_vn=vitesse_v*normale;
2158         double vitesse_ln=vitesse_l*normale;
2159         //double rho_m_e_m=rho_m_E_m -0.5*(q_m*vitesse_melange + rho_m*concentration_v*concentration_l*vitesse_relative*vitesse_relative)
2160         double rho_v=_fluides[0]->getDensity(pression,Temperature);
2161         double rho_l=_fluides[1]->getDensity(pression,Temperature);
2162         double e_v=_fluides[0]->getInternalEnergy(Temperature,rho_v);
2163         double e_l=_fluides[1]->getInternalEnergy(Temperature,rho_l);
2164
2165         Vector F(_nVar);
2166         F(0)=m_v*vitesse_vn+m_l*vitesse_ln;
2167         F(1)=m_v*vitesse_vn;
2168         for(int i=0;i<_Ndim;i++)
2169                 F(2+i)=m_v*vitesse_vn*vitesse_v(i)+m_l*vitesse_ln*vitesse_l(i)+pression*normale(i)*porosity;
2170         F(2+_Ndim)=m_v*(e_v+0.5*vitesse_v*vitesse_v+pression/rho_v)*vitesse_vn+m_l*(e_l+0.5*vitesse_l*vitesse_l+pression/rho_l)*vitesse_ln;
2171
2172         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2173         {
2174                 cout<<"DriftModel::convectionFlux end"<<endl;
2175                 cout<<"Flux F(U,V)"<<endl;
2176                 cout<<F<<endl;
2177         }
2178
2179         return F;
2180 }
2181
2182 Vector DriftModel::staggeredVFFCFlux()
2183 {
2184         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2185                 cout<<"DriftModel::staggeredVFFCFlux()start"<<endl;
2186
2187         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2188         {
2189                 *_runLogFile<<"DriftModel::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding"<<endl;
2190                 _runLogFile->close();
2191                 throw CdmathException("DriftModel::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
2192         }
2193         else//_spaceScheme==staggered
2194         {
2195                 Vector Fij(_nVar);
2196
2197                 double uijn=0, phiqn=0, uin=0, ujn=0;
2198                 for(int idim=0;idim<_Ndim;idim++)
2199                 {
2200                         uijn+=_vec_normal[idim]*_Uroe[2+idim];//URoe = rho, cm, u, H
2201                         uin+=_vec_normal[idim]*_Ui[2+idim];
2202                         ujn+=_vec_normal[idim]*_Uj[2+idim];
2203                 }
2204                 if( (uin>0 && ujn >0) || (uin>=0 && ujn <=0 && uijn>0) ) // formerly (uijn>_precision)
2205                 {
2206                         for(int idim=0;idim<_Ndim;idim++)
2207                                 phiqn+=_vec_normal[idim]*_Ui[2+idim];//phi rho u n
2208                         Fij(0)=phiqn;
2209                         Fij(1)=_Vi[0]*phiqn;
2210                         for(int idim=0;idim<_Ndim;idim++)
2211                                 Fij(2+idim)=phiqn*_Vi[2+idim]+_Vj[1]*_vec_normal[idim]*_porosityj;
2212                         Fij(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[1]*sqrt(_porosityj/_porosityi));
2213                 }
2214                 else if( (uin<0 && ujn <0) || (uin>=0 && ujn <=0 && uijn<0) ) // formerly (uijn<-_precision)
2215                 {
2216                         for(int idim=0;idim<_Ndim;idim++)
2217                                 phiqn+=_vec_normal[idim]*_Uj[2+idim];//phi rho u n
2218                         Fij(0)=phiqn;
2219                         Fij(1)=_Vj[0]*phiqn;
2220                         for(int idim=0;idim<_Ndim;idim++)
2221                                 Fij(2+idim)=phiqn*_Vj[2+idim]+_Vi[1]*_vec_normal[idim]*_porosityi;
2222                         Fij(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[1]*sqrt(_porosityi/_porosityj));
2223                 }
2224                 else//case (uin<=0 && ujn >=0) or (uin>=0 && ujn <=0 && uijn==0), apply centered scheme
2225                 {
2226                         Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar);
2227                         Vector normale(_Ndim);
2228                         for(int i1=0;i1<_Ndim;i1++)
2229                                 normale(i1)=_vec_normal[i1];
2230                         for(int i1=0;i1<_nVar;i1++)
2231                         {
2232                                 Ui(i1)=_Ui[i1];
2233                                 Uj(i1)=_Uj[i1];
2234                                 Vi(i1)=_Vi[i1];
2235                                 Vj(i1)=_Vj[i1];
2236                         }
2237                         Fi=convectionFlux(Ui,Vi,normale,_porosityi);
2238                         Fj=convectionFlux(Uj,Vj,normale,_porosityj);
2239                         Fij=(Fi+Fj)/2+_maxvploc*(Ui-Uj)/2;
2240
2241                         //case nil velocity on the interface, apply smoothed scheme
2242                         /*
2243                         double theta=uijn/(.01);
2244                         Vector F1(_nVar),F2(_nVar);
2245                         for(int idim=0;idim<_Ndim;idim++)
2246                                 phiqn+=_vec_normal[idim]*_Ui[2+idim];//phi rho u n
2247                         F1(0)=phiqn;
2248                         F1(1)=_Vi[0]*phiqn;
2249                         for(int idim=0;idim<_Ndim;idim++)
2250                                 F1(2+idim)=phiqn*_Vi[2+idim]+_Vj[1]*_vec_normal[idim]*_porosityj;
2251                         F1(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[1]*sqrt(_porosityj/_porosityi));
2252
2253                         for(int idim=0;idim<_Ndim;idim++)
2254                                 phiqn+=_vec_normal[idim]*_Uj[2+idim];//phi rho u n
2255                         F2(0)=phiqn;
2256                         F2(1)=_Vj[0]*phiqn;
2257                         for(int idim=0;idim<_Ndim;idim++)
2258                                 F2(2+idim)=phiqn*_Vj[2+idim]+_Vi[1]*_vec_normal[idim]*_porosityi;
2259                         F2(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[1]*sqrt(_porosityi/_porosityj));
2260
2261                         Fij=(1+theta)/2*F1+(1-theta)/2*F2;
2262                          */
2263                 }
2264                 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2265                 {
2266                         cout<<"DriftModel::staggeredVFFCFlux() end uijn="<<uijn<<endl;
2267                         cout<<Fij<<endl;
2268                 }
2269                 return Fij;
2270         }
2271 }
2272
2273 void DriftModel::staggeredVFFCMatricesConservativeVariables(double u_mn)
2274 {
2275         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2276                 cout<<"DriftModel::staggeredVFFCMatricesConservativeVariables()"<<endl;
2277
2278         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2279         {
2280                 *_runLogFile<<"DriftModel::staggeredVFFCMatricesConservativeVariables: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding"<<endl;
2281                 _runLogFile->close();
2282                 throw CdmathException("DriftModel::staggeredVFFCMatricesConservativeVariables: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding");
2283         }
2284         else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2285         {
2286                 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
2287                 double H;//enthalpie totale (expression particulière)
2288                 consToPrim(_Ui,_Vi,_porosityi);
2289                 consToPrim(_Uj,_Vj,_porosityj);
2290
2291                 for(int i=0;i<_Ndim;i++)
2292                 {
2293                         ui_2 += _Vi[2+i]*_Vi[2+i];
2294                         ui_n += _Vi[2+i]*_vec_normal[i];
2295                         uj_2 += _Vj[2+i]*_Vj[2+i];
2296                         uj_n += _Vj[2+i]*_vec_normal[i];
2297                 }
2298
2299                 double rhomi=_Ui[0]/_porosityi;
2300                 double mi_v=_Ui[1]/_porosityi;
2301                 double mi_l=rhomi-mi_v;
2302                 double cmi=_Vi[0];
2303                 double pi=_Vi[1];
2304                 double Tmi=_Vi[_Ndim+2];
2305                 double Emi=_Ui[_Ndim+2]/(rhomi*_porosityi);
2306                 double ecini=0.5*ui_2;
2307                 double emi=Emi-0.5*ui_2;
2308                 double hmi=emi+pi/rhomi;
2309
2310                 double pj=_Vj[1];
2311                 double rhomj=_Uj[0]/_porosityj;
2312                 double mj_v =_Uj[1]/_porosityj;
2313                 double mj_l=rhomj-mj_v;
2314                 double cmj=_Vj[0];
2315                 double Tmj=_Vj[_Ndim+2];
2316                 double Emj=_Uj[_Ndim+2]/(rhomj*_porosityj);
2317                 double ecinj=0.5*uj_2;
2318                 double emj=Emj-0.5*uj_2;
2319                 double hmj=emj+pj/rhomj;
2320
2321                 double Hm;//value depends on the sign of interfacial velocity
2322
2323                 if(u_mn>_precision)
2324                 {
2325                         Hm=Emi+pj/rhomi;
2326                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2327                                 cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Hm= "<<Hm<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
2328
2329                         /***********Calcul des valeurs propres ********/
2330                         vector<complex<double> > vp_dist(3,0);
2331                         getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2332                         if(_kappa*hmj+_khi+cmj*_ksi<0){
2333                                 *_runLogFile<<"staggeredVFFCMatricesConservativeVariables: vitesse du son complexe"<<endl;
2334                                 _runLogFile->close();
2335                                 throw CdmathException("staggeredVFFCMatricesConservativeVariables: vitesse du son complexe");
2336                         }
2337                         double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
2338
2339                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2340                                 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
2341
2342                         //On remplit les valeurs propres
2343                         vp_dist[0]=ui_n+amj;
2344                         vp_dist[1]=ui_n-amj;
2345                         vp_dist[2]=ui_n;
2346
2347                         _maxvploc=fabs(ui_n)+amj;
2348                         if(_maxvploc>_maxvp)
2349                                 _maxvp=_maxvploc;
2350
2351                         /******** Construction de la matrice A^+ *********/
2352                         _AroePlus[0*_nVar+0]=0;
2353                         _AroePlus[0*_nVar+1]=0;
2354                         for(int i=0;i<_Ndim;i++)
2355                                 _AroePlus[0*_nVar+2+i]=_vec_normal[i];
2356                         _AroePlus[0*_nVar+2+_Ndim]=0;
2357                         _AroePlus[1*_nVar+0]=-ui_n*cmi;
2358                         _AroePlus[1*_nVar+1]=ui_n;
2359                         for(int i=0;i<_Ndim;i++)
2360                                 _AroePlus[1*_nVar+2+i]=cmi*_vec_normal[i];
2361                         _AroePlus[1*_nVar+2+_Ndim]=0;
2362                         for(int i=0;i<_Ndim;i++)
2363                         {
2364                                 _AroePlus[(2+i)*_nVar+0]=-ui_n*_Vi[2+i];
2365                                 _AroePlus[(2+i)*_nVar+1]=0;
2366                                 for(int j=0;j<_Ndim;j++)
2367                                         _AroePlus[(2+i)*_nVar+2+j]=_Vi[2+i]*_vec_normal[j];
2368                                 _AroePlus[(2+i)*_nVar+2+i]+=ui_n;
2369                                 _AroePlus[(2+i)*_nVar+2+_Ndim]=0;
2370                         }
2371                         _AroePlus[(2+_Ndim)*_nVar+0]=-Hm*ui_n;
2372                         _AroePlus[(2+_Ndim)*_nVar+1]=0;
2373                         for(int i=0;i<_Ndim;i++)
2374                                 _AroePlus[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2375                         _AroePlus[(2+_Ndim)*_nVar+2+_Ndim]=ui_n;
2376
2377                         /******** Construction de la matrice A^- *********/
2378                         _AroeMinus[0*_nVar+0]=0;
2379                         _AroeMinus[0*_nVar+1]=0;
2380                         for(int i=0;i<_Ndim;i++)
2381                                 _AroeMinus[0*_nVar+2+i]=0;
2382                         _AroeMinus[0*_nVar+2+_Ndim]=0;
2383                         _AroeMinus[1*_nVar+0]=0;
2384                         _AroeMinus[1*_nVar+1]=0;
2385                         for(int i=0;i<_Ndim;i++)
2386                                 _AroeMinus[1*_nVar+2+i]=0;
2387                         _AroeMinus[1*_nVar+2+_Ndim]=0;
2388                         for(int i=0;i<_Ndim;i++)
2389                         {
2390                                 _AroeMinus[(2+i)*_nVar+0]=(_khi+_kappa*ecinj)*_vec_normal[i];
2391                                 _AroeMinus[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2392                                 for(int j=0;j<_Ndim;j++)
2393                                         _AroeMinus[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vj[2+j];
2394                                 _AroeMinus[(2+i)*_nVar+2+i]+=0;
2395                                 _AroeMinus[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2396                         }
2397                         _AroeMinus[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecinj)*ui_n;
2398                         _AroeMinus[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2399                         for(int i=0;i<_Ndim;i++)
2400                                 _AroeMinus[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vi[2+i];
2401                         _AroeMinus[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2402                 }
2403                 else if(u_mn<-_precision)
2404                 {
2405                         Hm=Emj+pi/rhomj;
2406                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2407                                 cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Hm= "<<Hm<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
2408
2409                         /***********Calcul des valeurs propres ********/
2410                         vector<complex<double> > vp_dist(3,0);
2411                         getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
2412                         if(_kappa*hmi+_khi+cmi*_ksi<0)
2413                         {
2414                                 *_runLogFile<<"staggeredVFFCMatricesConservativeVariables: vitesse du son complexe"<<endl;
2415                                 _runLogFile->close();
2416                                 throw CdmathException("staggeredVFFCMatricesConservativeVariables: vitesse du son complexe");
2417                         }
2418                         double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
2419                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2420                                 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
2421
2422                         //On remplit les valeurs propres
2423                         vp_dist[0]=uj_n+ami;
2424                         vp_dist[1]=uj_n-ami;
2425                         vp_dist[2]=uj_n;
2426
2427                         _maxvploc=fabs(uj_n)+ami;
2428                         if(_maxvploc>_maxvp)
2429                                 _maxvp=_maxvploc;
2430
2431                         /******** Construction de la matrice A^+ *********/
2432                         _AroePlus[0*_nVar+0]=0;
2433                         _AroePlus[0*_nVar+1]=0;
2434                         for(int i=0;i<_Ndim;i++)
2435                                 _AroePlus[0*_nVar+2+i]=0;
2436                         _AroePlus[0*_nVar+2+_Ndim]=0;
2437                         _AroePlus[1*_nVar+0]=0;
2438                         _AroePlus[1*_nVar+1]=0;
2439                         for(int i=0;i<_Ndim;i++)
2440                                 _AroePlus[1*_nVar+2+i]=0;
2441                         _AroePlus[1*_nVar+2+_Ndim]=0;
2442                         for(int i=0;i<_Ndim;i++)
2443                         {
2444                                 _AroePlus[(2+i)*_nVar+0]=(_khi+_kappa*ecini)*_vec_normal[i];
2445                                 _AroePlus[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2446                                 for(int j=0;j<_Ndim;j++)
2447                                         _AroePlus[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vi[2+j];
2448                                 _AroePlus[(2+i)*_nVar+2+i]+=0;
2449                                 _AroePlus[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2450                         }
2451                         _AroePlus[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecini)*ui_n;
2452                         _AroePlus[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2453                         for(int i=0;i<_Ndim;i++)
2454                                 _AroePlus[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vj[2+i];
2455                         _AroePlus[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2456
2457                         /******** Construction de la matrice A^- *********/
2458                         _AroeMinus[0*_nVar+0]=0;
2459                         _AroeMinus[0*_nVar+1]=0;
2460                         for(int i=0;i<_Ndim;i++)
2461                                 _AroeMinus[0*_nVar+2+i]=_vec_normal[i];
2462                         _AroeMinus[0*_nVar+2+_Ndim]=0;
2463                         _AroeMinus[1*_nVar+0]=-uj_n*cmj;
2464                         _AroeMinus[1*_nVar+1]=uj_n;
2465                         for(int i=0;i<_Ndim;i++)
2466                                 _AroeMinus[1*_nVar+2+i]=cmj*_vec_normal[i];
2467                         _AroeMinus[1*_nVar+2+_Ndim]=0;
2468                         for(int i=0;i<_Ndim;i++)
2469                         {
2470                                 _AroeMinus[(2+i)*_nVar+0]=-uj_n*_Vj[2+i];
2471                                 _AroeMinus[(2+i)*_nVar+1]=0;
2472                                 for(int j=0;j<_Ndim;j++)
2473                                         _AroeMinus[(2+i)*_nVar+2+j]=_Vj[2+i]*_vec_normal[j];
2474                                 _AroeMinus[(2+i)*_nVar+2+i]+=uj_n;
2475                                 _AroeMinus[(2+i)*_nVar+2+_Ndim]=0;
2476                         }
2477                         _AroeMinus[(2+_Ndim)*_nVar+0]=-Hm*uj_n;
2478                         _AroeMinus[(2+_Ndim)*_nVar+1]=0;
2479                         for(int i=0;i<_Ndim;i++)
2480                                 _AroeMinus[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2481                         _AroeMinus[(2+_Ndim)*_nVar+2+_Ndim]=uj_n;
2482                 }
2483                 else//case nil velocity on the interface
2484                 {
2485                         double rhom=_Uroe[0];
2486                         double cm=_Uroe[1];
2487                         double Hm=_Uroe[_nVar-1];
2488                         double m_v=cm*rhom, m_l=rhom-m_v;
2489                         double Tm;
2490                         double umn=0,u_2=0;
2491                         Vector vitesse(_Ndim);
2492                         for(int idim=0;idim<_Ndim;idim++)
2493                         {
2494                                 vitesse[idim]=_Uroe[2+idim];
2495                                 umn += _Uroe[2+idim]*_vec_normal[idim];//vitesse normale
2496                                 u_2 += _Uroe[2+idim]*_Uroe[2+idim];
2497                         }
2498                         double hm=Hm-0.5*u_2;
2499
2500                         if(cm<_precision)//pure liquid
2501                                 Tm=_fluides[1]->getTemperatureFromEnthalpy(hm,rhom);
2502                         else if(cm>1-_precision)
2503                                 Tm=_fluides[0]->getTemperatureFromEnthalpy(hm,rhom);
2504                         else//Hypothèse de saturation
2505                                 Tm=_Tsat;
2506
2507                         double pression= getMixturePressure(cm, rhom, Tm);
2508
2509                         getMixturePressureDerivatives( m_v, m_l, pression, Tm);//EOS is involved to express pressure jump and sound speed
2510                         if(_kappa*hm+_khi+cm*_ksi<0){
2511                                 *_runLogFile<<"Calcul matrice de Roe: vitesse du son complexe"<<endl;
2512                                 _runLogFile->close();
2513                                 throw CdmathException("Calcul matrice de Roe: vitesse du son complexe");
2514                         }
2515                         double am=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
2516                         _maxvploc=fabs(umn)+am;
2517                         if(_maxvploc>_maxvp)
2518                                 _maxvp=_maxvploc;
2519
2520                         //Rusanov scheme
2521                         /******* Construction de la matrice A^+ ********/
2522                         //A^+=0.5*jacobienne Flux (U_i)+0.5 _maxvploc Id
2523                         Hm=Emi+pi/rhomi;
2524                         getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
2525                         _AroePlus[0*_nVar+0]=0;
2526                         _AroePlus[0*_nVar+1]=0;
2527                         for(int i=0;i<_Ndim;i++)
2528                                 _AroePlus[0*_nVar+2+i]=0.5*_vec_normal[i];
2529                         _AroePlus[0*_nVar+2+_Ndim]=0;
2530                         _AroePlus[1*_nVar+0]=-0.5*ui_n*cmi;
2531                         _AroePlus[1*_nVar+1]=0.5*ui_n;
2532                         for(int i=0;i<_Ndim;i++)
2533                                 _AroePlus[1*_nVar+2+i]=0.5*cmi*_vec_normal[i];
2534                         _AroePlus[1*_nVar+2+_Ndim]=0;
2535                         for(int i=0;i<_Ndim;i++)
2536                         {
2537                                 _AroePlus[(2+i)*_nVar+0]=0.5*(_khi+_kappa*ecini)*_vec_normal[i]-0.5*ui_n*_Vi[2+i];
2538                                 _AroePlus[(2+i)*_nVar+1]=0.5*_ksi*_vec_normal[i];
2539                                 for(int j=0;j<_Ndim;j++)
2540                                         _AroePlus[(2+i)*_nVar+2+j]=0.5*_Vi[2+i]*_vec_normal[j]-0.5*_kappa*_vec_normal[i]*_Vi[2+j];
2541                                 _AroePlus[(2+i)*_nVar+2+i]+=0.5*ui_n;
2542                                 _AroePlus[(2+i)*_nVar+2+_Ndim]=0.5*_kappa*_vec_normal[i];
2543                         }
2544                         _AroePlus[(2+_Ndim)*_nVar+0]=-0.5*Hm*ui_n;
2545                         _AroePlus[(2+_Ndim)*_nVar+1]=0;
2546                         for(int i=0;i<_Ndim;i++)
2547                                 _AroePlus[(2+_Ndim)*_nVar+2+i]=0.5*Hm*_vec_normal[i];
2548                         _AroePlus[(2+_Ndim)*_nVar+2+_Ndim]=0.5*ui_n;
2549
2550                         for(int i=0;i<_nVar;i++)
2551                                 _AroePlus[i*_nVar+i]+=_maxvploc/2;
2552                         /****** Construction de la matrice A^- *******/
2553                         //A^-=0.5*jacobienne Flux (U_j)-0.5 _maxvploc Id
2554                         Hm=Emj+pj/rhomj;
2555                         getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2556                         _AroeMinus[0*_nVar+0]=0;
2557                         _AroeMinus[0*_nVar+1]=0;
2558                         for(int i=0;i<_Ndim;i++)
2559                                 _AroeMinus[0*_nVar+2+i]=0.5*_vec_normal[i];
2560                         _AroeMinus[0*_nVar+2+_Ndim]=0;
2561                         _AroeMinus[1*_nVar+0]=-0.5*uj_n*cmj;
2562                         _AroeMinus[1*_nVar+1]=0.5*uj_n;
2563                         for(int i=0;i<_Ndim;i++)
2564                                 _AroeMinus[1*_nVar+2+i]=0.5*cmj*_vec_normal[i];
2565                         _AroeMinus[1*_nVar+2+_Ndim]=0;
2566                         for(int i=0;i<_Ndim;i++)
2567                         {
2568                                 _AroeMinus[(2+i)*_nVar+0]=0.5*(_khi+_kappa*ecinj)*_vec_normal[i]-0.5*uj_n*_Vj[2+i];
2569                                 _AroeMinus[(2+i)*_nVar+1]=0.5*_ksi*_vec_normal[i];
2570                                 for(int j=0;j<_Ndim;j++)
2571                                         _AroeMinus[(2+i)*_nVar+2+j]=0.5*_Vj[2+i]*_vec_normal[j]-0.5*_kappa*_vec_normal[i]*_Vj[2+j];
2572                                 _AroeMinus[(2+i)*_nVar+2+i]+=0.5*uj_n;
2573                                 _AroeMinus[(2+i)*_nVar+2+_Ndim]=0.5*_kappa*_vec_normal[i];
2574                         }
2575                         _AroeMinus[(2+_Ndim)*_nVar+0]=-0.5*Hm*uj_n;
2576                         _AroeMinus[(2+_Ndim)*_nVar+1]=0;
2577                         for(int i=0;i<_Ndim;i++)
2578                                 _AroeMinus[(2+_Ndim)*_nVar+2+i]=0.5*Hm*_vec_normal[i];
2579                         _AroeMinus[(2+_Ndim)*_nVar+2+_Ndim]=0.5*uj_n;
2580
2581                         for(int i=0;i<_nVar;i++)
2582                                 _AroeMinus[i*_nVar+i]-=_maxvploc/2;
2583
2584                         /*
2585                         double _AroePlusi[_nVar*_nVar], _AroePlusj[_nVar*_nVar];
2586                         double _AroeMinusi[_nVar*_nVar], _AroeMinusj[_nVar*_nVar];
2587
2588                         //Matrices côté gauche
2589                         Hm=Emi+pj/rhomi;
2590                          ****** Construction de la matrice A^+ i (gauche) *******
2591                         _AroePlusi[0*_nVar+0]=0;
2592                         _AroePlusi[0*_nVar+1]=0;
2593                         for(int i=0;i<_Ndim;i++)
2594                                 _AroePlusi[0*_nVar+2+i]=_vec_normal[i];
2595                         _AroePlusi[0*_nVar+2+_Ndim]=0;
2596                         _AroePlusi[1*_nVar+0]=-ui_n*cmi;
2597                         _AroePlusi[1*_nVar+1]=ui_n;
2598                         for(int i=0;i<_Ndim;i++)
2599                                 _AroePlusi[1*_nVar+2+i]=cmi*_vec_normal[i];
2600                         _AroePlusi[1*_nVar+2+_Ndim]=0;
2601                         for(int i=0;i<_Ndim;i++)
2602                         {
2603                                 _AroePlusi[(2+i)*_nVar+0]=-ui_n*_Vi[2+i];
2604                                 _AroePlusi[(2+i)*_nVar+1]=0;
2605                                 for(int j=0;j<_Ndim;j++)
2606                                         _AroePlusi[(2+i)*_nVar+2+j]=_Vi[2+i]*_vec_normal[j];
2607                                 _AroePlusi[(2+i)*_nVar+2+i]+=ui_n;
2608                                 _AroePlusi[(2+i)*_nVar+2+_Ndim]=0;
2609                         }
2610                         _AroePlusi[(2+_Ndim)*_nVar+0]=-Hm*ui_n;
2611                         _AroePlusi[(2+_Ndim)*_nVar+1]=0;
2612                         for(int i=0;i<_Ndim;i++)
2613                                 _AroePlusi[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2614                         _AroePlusi[(2+_Ndim)*_nVar+2+_Ndim]=ui_n;
2615
2616                          ****** Construction de la matrice A^- i (coté gauche)*******
2617                         _AroeMinusi[0*_nVar+0]=0;
2618                         _AroeMinusi[0*_nVar+1]=0;
2619                         for(int i=0;i<_Ndim;i++)
2620                                 _AroeMinusi[0*_nVar+2+i]=0;
2621                         _AroeMinusi[0*_nVar+2+_Ndim]=0;
2622                         _AroeMinusi[1*_nVar+0]=0;
2623                         _AroeMinusi[1*_nVar+1]=0;
2624                         for(int i=0;i<_Ndim;i++)
2625                                 _AroeMinusi[1*_nVar+2+i]=0;
2626                         _AroeMinusi[1*_nVar+2+_Ndim]=0;
2627                         for(int i=0;i<_Ndim;i++)
2628                         {
2629                                 _AroeMinusi[(2+i)*_nVar+0]=(_khi+_kappa*ecinj)*_vec_normal[i];
2630                                 _AroeMinusi[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2631                                 for(int j=0;j<_Ndim;j++)
2632                                         _AroeMinusi[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vj[2+j];
2633                                 _AroeMinusi[(2+i)*_nVar+2+i]+=0;
2634                                 _AroeMinusi[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2635                         }
2636                         _AroeMinusi[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecinj)*ui_n;
2637                         _AroeMinusi[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2638                         for(int i=0;i<_Ndim;i++)
2639                                 _AroeMinusi[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vi[2+i];
2640                         _AroeMinusi[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2641
2642                         //Matrices côté droit
2643                         Hm=Emj+pi/rhomj;
2644
2645                          ****** Construction de la matrice A^+ j (coté droit)*******
2646                         _AroePlusj[0*_nVar+0]=0;
2647                         _AroePlusj[0*_nVar+1]=0;
2648                         for(int i=0;i<_Ndim;i++)
2649                                 _AroePlusj[0*_nVar+2+i]=0;
2650                         _AroePlusj[0*_nVar+2+_Ndim]=0;
2651                         _AroePlusj[1*_nVar+0]=0;
2652                         _AroePlusj[1*_nVar+1]=0;
2653                         for(int i=0;i<_Ndim;i++)
2654                                 _AroePlusj[1*_nVar+2+i]=0;
2655                         _AroePlusj[1*_nVar+2+_Ndim]=0;
2656                         for(int i=0;i<_Ndim;i++)
2657                         {
2658                                 _AroePlusj[(2+i)*_nVar+0]=(_khi+_kappa*ecini)*_vec_normal[i];
2659                                 _AroePlusj[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
2660                                 for(int j=0;j<_Ndim;j++)
2661                                         _AroePlusj[(2+i)*_nVar+2+j]=-_kappa*_vec_normal[i]*_Vi[2+j];
2662                                 _AroePlusj[(2+i)*_nVar+2+i]+=0;
2663                                 _AroePlusj[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
2664                         }
2665                         _AroePlusj[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecini)*ui_n;
2666                         _AroePlusj[(2+_Ndim)*_nVar+1]=_ksi*ui_n;
2667                         for(int i=0;i<_Ndim;i++)
2668                                 _AroePlusj[(2+_Ndim)*_nVar+2+i]=-_kappa*uj_n*_Vj[2+i];
2669                         _AroePlusj[(2+_Ndim)*_nVar+2+_Ndim]=_kappa*ui_n;
2670
2671                          ****** Construction de la matrice A^- j (coté droit)*******
2672                         _AroeMinusj[0*_nVar+0]=0;
2673                         _AroeMinusj[0*_nVar+1]=0;
2674                         for(int i=0;i<_Ndim;i++)
2675                                 _AroeMinusj[0*_nVar+2+i]=_vec_normal[i];
2676                         _AroeMinusj[0*_nVar+2+_Ndim]=0;
2677                         _AroeMinusj[1*_nVar+0]=-uj_n*cmj;
2678                         _AroeMinusj[1*_nVar+1]=uj_n;
2679                         for(int i=0;i<_Ndim;i++)
2680                                 _AroeMinusj[1*_nVar+2+i]=cmj*_vec_normal[i];
2681                         _AroeMinusj[1*_nVar+2+_Ndim]=0;
2682                         for(int i=0;i<_Ndim;i++)
2683                         {
2684                                 _AroeMinusj[(2+i)*_nVar+0]=-uj_n*_Vj[2+i];
2685                                 _AroeMinusj[(2+i)*_nVar+1]=0;
2686                                 for(int j=0;j<_Ndim;j++)
2687                                         _AroeMinusj[(2+i)*_nVar+2+j]=_Vj[2+i]*_vec_normal[j];
2688                                 _AroeMinusj[(2+i)*_nVar+2+i]+=uj_n;
2689                                 _AroeMinusj[(2+i)*_nVar+2+_Ndim]=0;
2690                         }
2691                         _AroeMinusj[(2+_Ndim)*_nVar+0]=-Hm*uj_n;
2692                         _AroeMinusj[(2+_Ndim)*_nVar+1]=0;
2693                         for(int i=0;i<_Ndim;i++)
2694                                 _AroeMinusj[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i];
2695                         _AroeMinusj[(2+_Ndim)*_nVar+2+_Ndim]=uj_n;
2696
2697                         double theta=u_mn/(.01);
2698                         for(int i=0; i<_nVar*_nVar;i++)
2699                         {
2700                                 _AroePlus[i] =(1+theta)/2*_AroePlusi[i] +(1-theta)/2*_AroePlusj[i];
2701                                 _AroeMinus[i]=(1+theta)/2*_AroeMinusi[i]+(1-theta)/2*_AroeMinusj[i];
2702                         }
2703                          */
2704                 }
2705         }
2706         for(int i=0; i<_nVar*_nVar;i++)
2707         {
2708                 _Aroe[i]=_AroePlus[i]+_AroeMinus[i];
2709                 _absAroe[i]=_AroePlus[i]-_AroeMinus[i];
2710         }
2711
2712         if(_timeScheme==Implicit)
2713                 for(int i=0; i<_nVar*_nVar;i++)
2714                 {
2715                         _AroeMinusImplicit[i] = _AroeMinus[i];
2716                         _AroePlusImplicit[i]  = _AroePlus[i];
2717                 }
2718 }
2719
2720 void DriftModel::staggeredVFFCMatricesPrimitiveVariables(double u_mn)
2721 {
2722         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2723                 cout<<"DriftModel::staggeredVFFCMatricesPrimitiveVariables()"<<endl;
2724
2725         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2726         {
2727                 *_runLogFile<< "DriftModel::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding" << endl;
2728                 _runLogFile->close();
2729                 throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding");
2730         }
2731         else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2732         {
2733                 //Calls to getDensityDerivatives needed
2734                 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
2735                 double H;//enthalpie totale (expression particulière)
2736                 consToPrim(_Ui,_Vi,_porosityi);
2737                 consToPrim(_Uj,_Vj,_porosityj);
2738
2739                 for(int i=0;i<_Ndim;i++)
2740                 {
2741                         ui_2 += _Vi[2+i]*_Vi[2+i];
2742                         ui_n += _Vi[2+i]*_vec_normal[i];
2743                         uj_2 += _Vj[2+i]*_Vj[2+i];
2744                         uj_n += _Vj[2+i]*_vec_normal[i];
2745                 }
2746
2747                 double rhomi=_Ui[0]/_porosityi;
2748                 double mi_v=_Ui[1]/_porosityi;
2749                 double mi_l=rhomi-mi_v;
2750                 double cmi=_Vi[0];
2751                 double pi=_Vi[1];
2752                 double Tmi=_Vi[_Ndim+2];
2753                 double Emi=_Ui[_Ndim+2]/(rhomi*_porosityi);
2754                 double ecini=0.5*ui_2;
2755                 double emi=Emi-0.5*ui_2;
2756                 double hmi=emi+pi/rhomi;
2757                 double rho_vi=_fluides[0]->getDensity(pi,Tmi);
2758                 double rho_li=_fluides[1]->getDensity(pi,Tmi);
2759
2760                 double pj=_Vj[1];
2761                 double rhomj=_Uj[0]/_porosityj;
2762                 double mj_v =_Uj[1]/_porosityj;
2763                 double mj_l=rhomj-mj_v;
2764                 double cmj=_Vj[0];
2765                 double Tmj=_Vj[_Ndim+2];
2766                 double Emj=_Uj[_Ndim+2]/(rhomj*_porosityj);
2767                 double ecinj=0.5*uj_2;
2768                 double emj=Emj-0.5*uj_2;
2769                 double hmj=emj+pj/rhomj;
2770                 double rho_vj=_fluides[0]->getDensity(pj,Tmj);
2771                 double rho_lj=_fluides[1]->getDensity(pj,Tmj);
2772
2773                 double gamma_v=_fluides[0]->constante("gamma");
2774                 double gamma_l=_fluides[1]->constante("gamma");
2775                 double q_v=_fluides[0]->constante("q");
2776                 double q_l=_fluides[1]->constante("q");
2777
2778                 if(fabs(u_mn)>_precision)//non zero velocity on the interface
2779                 {
2780                         if(             !_useDellacherieEOS)
2781                         {
2782                                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2783                                 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
2784                                 double cv_v=fluide0->constante("cv");
2785                                 double cv_l=fluide1->constante("cv");
2786
2787                                 double e_vi=_fluides[0]->getInternalEnergy(Tmi,rho_vi);
2788                                 double e_li=_fluides[1]->getInternalEnergy(Tmi,rho_li);
2789                                 double e_vj=_fluides[0]->getInternalEnergy(Tmj,rho_vj);
2790                                 double e_lj=_fluides[1]->getInternalEnergy(Tmj,rho_lj);
2791
2792                                 if(u_mn>_precision)
2793                                 {
2794                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2795                                                 cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Emi= "<<Emi<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
2796
2797                                         /***********Calcul des valeurs propres ********/
2798                                         vector<complex<double> > vp_dist(3,0);
2799                                         getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2800                                         if(_kappa*hmj+_khi+cmj*_ksi<0)
2801                                         {
2802                                                 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
2803                                                 _runLogFile->close();
2804                                                 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
2805                                         }
2806                                         double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
2807
2808                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2809                                                 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
2810
2811                                         //On remplit les valeurs propres
2812                                         vp_dist[0]=ui_n+amj;
2813                                         vp_dist[1]=ui_n-amj;
2814                                         vp_dist[2]=ui_n;
2815
2816                                         _maxvploc=fabs(ui_n)+amj;
2817                                         if(_maxvploc>_maxvp)
2818                                                 _maxvp=_maxvploc;
2819
2820                                         /******** Construction de la matrice A^+ *********/
2821                                         _AroePlusImplicit[0*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
2822                                         _AroePlusImplicit[0*_nVar+1]= rhomi*rhomi*ui_n*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)));
2823                                         for(int i=0;i<_Ndim;i++)
2824                                                 _AroePlusImplicit[0*_nVar+2+i]=rhomi*_vec_normal[i];
2825                                         _AroePlusImplicit[0*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cv_v*cmi/(rho_vi*(e_vi-q_v))+cv_l*(1-cmi)/(rho_li*(e_li-q_l)));
2826                                         _AroePlusImplicit[1*_nVar+0]=rhomi*ui_n-cmi*rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
2827                                         _AroePlusImplicit[1*_nVar+1]=-cmi*rhomi*rhomi*ui_n*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)));
2828                                         for(int i=0;i<_Ndim;i++)
2829                                                 _AroePlusImplicit[1*_nVar+2+i]=cmi*rhomi*_vec_normal[i];
2830                                         _AroePlusImplicit[1*_nVar+2+_Ndim]=-cmi*rhomi*rhomi*ui_n*(cv_v*cmi/(rho_vi*(e_vi-q_v))+cv_l*(1-cmi)/(rho_li*(e_li-q_l)));
2831                                         for(int i=0;i<_Ndim;i++)
2832                                         {
2833                                                 _AroePlusImplicit[(2+i)*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li)*_Vi[2+i];
2834                                                 _AroePlusImplicit[(2+i)*_nVar+1]=rhomi*rhomi*ui_n*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)))*_Vi[2+i];
2835                                                 for(int j=0;j<_Ndim;j++)
2836                                                         _AroePlusImplicit[(2+i)*_nVar+2+j]=rhomi*_Vi[2+i]*_vec_normal[j];
2837                                                 _AroePlusImplicit[(2+i)*_nVar+2+i]+=rhomi*ui_n;
2838                                                 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cv_v*cmi/(rho_vi*(e_vi-q_v))+(1-cmi)/(rho_li*(e_li-q_l)))*_Vi[2+i];
2839                                         }
2840                                         _AroePlusImplicit[(2+_Ndim)*_nVar+0]=ui_n*(rhomi*(e_vi-e_li)-Emi*rhomi*rhomi*(1/rho_vi-1/rho_li));
2841                                         _AroePlusImplicit[(2+_Ndim)*_nVar+1]=rhomi*rhomi*Emi*(cmi/(rho_vi*rho_vi*(gamma_v-1)*(e_vi-q_v))+(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(e_li-q_l)));
2842                                         for(int i=0;i<_Ndim;i++)
2843                                                 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomi*Emi+pj)*_vec_normal[i] + ui_n*rhomi*_Vi[2+i];
2844                                         _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=ui_n*(rhomi*(cv_v*cmi+cv_l*(1-cmi))-Emi*rhomi*rhomi*(cv_v*cmi/(rho_vi*(e_vi-q_v))+cv_l*(1-cmi)/(rho_li*(e_li-q_l))));
2845
2846
2847                                         /******** Construction de la matrice A^- *********/
2848                                         _AroeMinusImplicit[0*_nVar+0]=0;
2849                                         _AroeMinusImplicit[0*_nVar+1]=0;
2850                                         for(int i=0;i<_Ndim;i++)
2851                                                 _AroeMinusImplicit[0*_nVar+2+i]=0;
2852                                         _AroeMinusImplicit[0*_nVar+2+_Ndim]=0;
2853                                         _AroeMinusImplicit[1*_nVar+0]=0;
2854                                         _AroeMinusImplicit[1*_nVar+1]=0;
2855                                         for(int i=0;i<_Ndim;i++)
2856                                                 _AroeMinusImplicit[1*_nVar+2+i]=0;
2857                                         _AroeMinusImplicit[1*_nVar+2+_Ndim]=0;
2858                                         for(int i=0;i<_Ndim;i++)
2859                                         {
2860                                                 _AroeMinusImplicit[(2+i)*_nVar+0]=0;
2861                                                 _AroeMinusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
2862                                                 for(int j=0;j<_Ndim;j++)
2863                                                         _AroeMinusImplicit[(2+i)*_nVar+2+j]=0;
2864                                                 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=0;
2865                                                 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=0;
2866                                         }
2867                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=0;
2868                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=ui_n;
2869                                         for(int i=0;i<_Ndim;i++)
2870                                                 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=0;
2871                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
2872                                 }
2873                                 else if(u_mn<-_precision)
2874                                 {
2875                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2876                                                 cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Emj= "<<Emj<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
2877
2878                                         /***********Calcul des valeurs propres ********/
2879                                         vector<complex<double> > vp_dist(3,0);
2880                                         getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
2881                                         if(_kappa*hmi+_khi+cmi*_ksi<0)
2882                                         {
2883                                                 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
2884                                                 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
2885                                         }
2886                                         double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
2887                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2888                                                 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
2889
2890                                         //On remplit les valeurs propres
2891                                         vp_dist[0]=uj_n+ami;
2892                                         vp_dist[1]=uj_n-ami;
2893                                         vp_dist[2]=uj_n;
2894
2895                                         _maxvploc=fabs(uj_n)+ami;
2896                                         if(_maxvploc>_maxvp)
2897                                                 _maxvp=_maxvploc;
2898
2899                                         /******** Construction de la matrice A^+ *********/
2900                                         _AroePlusImplicit[0*_nVar+0]=0;
2901                                         _AroePlusImplicit[0*_nVar+1]=0;
2902                                         for(int i=0;i<_Ndim;i++)
2903                                                 _AroePlusImplicit[0*_nVar+2+i]=0;
2904                                         _AroePlusImplicit[0*_nVar+2+_Ndim]=0;
2905                                         _AroePlusImplicit[1*_nVar+0]=0;
2906                                         _AroePlusImplicit[1*_nVar+1]=0;
2907                                         for(int i=0;i<_Ndim;i++)
2908                                                 _AroePlusImplicit[1*_nVar+2+i]=0;
2909                                         _AroePlusImplicit[1*_nVar+2+_Ndim]=0;
2910                                         for(int i=0;i<_Ndim;i++)
2911                                         {
2912                                                 _AroePlusImplicit[(2+i)*_nVar+0]=0;
2913                                                 _AroePlusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
2914                                                 for(int j=0;j<_Ndim;j++)
2915                                                         _AroePlusImplicit[(2+i)*_nVar+2+j]=0;
2916                                                 _AroePlusImplicit[(2+i)*_nVar+2+i]+=0;
2917                                                 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=0;
2918                                         }
2919                                         _AroePlusImplicit[(2+_Ndim)*_nVar+0]=0;
2920                                         _AroePlusImplicit[(2+_Ndim)*_nVar+1]=uj_n;
2921                                         for(int i=0;i<_Ndim;i++)
2922                                                 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=0;
2923                                         _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
2924
2925                                         /******** Construction de la matrice A^- *********/
2926                                         _AroeMinusImplicit[0*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
2927                                         _AroeMinusImplicit[0*_nVar+1]= rhomj*rhomj*uj_n*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)));
2928                                         for(int i=0;i<_Ndim;i++)
2929                                                 _AroeMinusImplicit[0*_nVar+2+i]=rhomj*_vec_normal[i];
2930                                         _AroeMinusImplicit[0*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l)));
2931                                         _AroeMinusImplicit[1*_nVar+0]=rhomj*uj_n-cmj*rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
2932                                         _AroeMinusImplicit[1*_nVar+1]=-cmj*rhomj*rhomj*uj_n*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)));
2933                                         for(int i=0;i<_Ndim;i++)
2934                                                 _AroeMinusImplicit[1*_nVar+2+i]=cmj*rhomj*_vec_normal[i];
2935                                         _AroeMinusImplicit[1*_nVar+2+_Ndim]=-cmj*rhomj*rhomj*uj_n*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l)));
2936                                         for(int i=0;i<_Ndim;i++)
2937                                         {
2938                                                 _AroeMinusImplicit[(2+i)*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj)*_Vj[2+i];
2939                                                 _AroeMinusImplicit[(2+i)*_nVar+1]=rhomj*rhomj*uj_n*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)))*_Vj[2+i];
2940                                                 for(int j=0;j<_Ndim;j++)
2941                                                         _AroeMinusImplicit[(2+i)*_nVar+2+j]=rhomj*_Vj[2+i]*_vec_normal[j];
2942                                                 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=rhomj*uj_n;
2943                                                 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l)))*_Vj[2+i];
2944                                         }
2945                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=uj_n*(rhomj*(e_vj-e_lj)-Emj*rhomj*rhomj*(1/rho_vj-1/rho_lj));
2946                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=rhomj*rhomj*Emj*(cmj/(rho_vj*rho_vj*(gamma_v-1)*(e_vj-q_v))+(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(e_lj-q_l)));
2947                                         for(int i=0;i<_Ndim;i++)
2948                                                 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomj*Emj+pi)*_vec_normal[i] + uj_n*rhomj*_Vj[2+i];
2949                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=uj_n*(rhomj*(cv_v*cmj+cv_l*(1-cmj))-Emj*rhomj*rhomj*(cv_v*cmj/(rho_vj*(e_vj-q_v))+cv_l*(1-cmj)/(rho_lj*(e_lj-q_l))));
2950                                 }
2951                                 else
2952                                 {
2953                                         *_runLogFile<< "DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero" << endl;
2954                                         _runLogFile->close();
2955                                         throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero");
2956                                 }
2957                         }
2958                         else if(_useDellacherieEOS)
2959                         {
2960                                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2961                                 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
2962                                 double cp_v=fluide0->constante("cp");
2963                                 double cp_l=fluide1->constante("cp");
2964
2965                                 double h_vi=_fluides[0]->getEnthalpy(Tmi,rho_vi);
2966                                 double h_li=_fluides[1]->getEnthalpy(Tmi,rho_li);
2967                                 double h_vj=_fluides[0]->getEnthalpy(Tmj,rho_vj);
2968                                 double h_lj=_fluides[1]->getEnthalpy(Tmj,rho_lj);
2969                                 double Hmi=Emi+pi/rho_vi;
2970                                 double Hmj=Emj+pj/rho_vj;
2971
2972                                 if(u_mn>_precision)
2973                                 {
2974                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2975                                                 cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Hmi= "<<Hmi<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
2976
2977                                         /***********Calcul des valeurs propres ********/
2978                                         vector<complex<double> > vp_dist(3,0);
2979                                         getMixturePressureDerivatives( mj_v, mj_l, pj, Tmj);
2980                                         if(_kappa*hmj+_khi+cmj*_ksi<0)
2981                                         {
2982                                                 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
2983                                                 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
2984                                         }
2985                                         double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
2986
2987                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2988                                                 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
2989
2990                                         //On remplit les valeurs propres
2991                                         vp_dist[0]=ui_n+amj;
2992                                         vp_dist[1]=ui_n-amj;
2993                                         vp_dist[2]=ui_n;
2994
2995                                         _maxvploc=fabs(ui_n)+amj;
2996                                         if(_maxvploc>_maxvp)
2997                                                 _maxvp=_maxvploc;
2998
2999                                         /******** Construction de la matrice A^+ *********/
3000                                         _AroePlusImplicit[0*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
3001                                         _AroePlusImplicit[0*_nVar+1]= rhomi*rhomi*ui_n*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)));
3002                                         for(int i=0;i<_Ndim;i++)
3003                                                 _AroePlusImplicit[0*_nVar+2+i]=rhomi*_vec_normal[i];
3004                                         _AroePlusImplicit[0*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l)));
3005                                         _AroePlusImplicit[1*_nVar+0]=rhomi*ui_n-cmi*rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li);
3006                                         _AroePlusImplicit[1*_nVar+1]=-cmi*rhomi*rhomi*ui_n*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)));
3007                                         for(int i=0;i<_Ndim;i++)
3008                                                 _AroePlusImplicit[1*_nVar+2+i]=cmi*rhomi*_vec_normal[i];
3009                                         _AroePlusImplicit[1*_nVar+2+_Ndim]=-cmi*rhomi*rhomi*ui_n*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l)));
3010                                         for(int i=0;i<_Ndim;i++)
3011                                         {
3012                                                 _AroePlusImplicit[(2+i)*_nVar+0]=-rhomi*rhomi*ui_n*(1/rho_vi-1/rho_li)*_Vi[2+i];
3013                                                 _AroePlusImplicit[(2+i)*_nVar+1]=rhomi*rhomi*ui_n*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)))*_Vi[2+i];
3014                                                 for(int j=0;j<_Ndim;j++)
3015                                                         _AroePlusImplicit[(2+i)*_nVar+2+j]=rhomi*_Vi[2+i]*_vec_normal[j];
3016                                                 _AroePlusImplicit[(2+i)*_nVar+2+i]+=rhomi*ui_n;
3017                                                 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomi*rhomi*ui_n*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l)))*_Vi[2+i];
3018                                         }
3019                                         _AroePlusImplicit[(2+_Ndim)*_nVar+0]=ui_n*(rhomi*(h_vi-h_li)-Hmi*rhomi*rhomi*(1/rho_vi-1/rho_li));
3020                                         _AroePlusImplicit[(2+_Ndim)*_nVar+1]=rhomi*rhomi*Hmi*(gamma_v*cmi/(rho_vi*rho_vi*(gamma_v-1)*(h_vi-q_v))+gamma_l*(1-cmi)/(rho_li*rho_li*(gamma_l-1)*(h_li-q_l)));
3021                                         for(int i=0;i<_Ndim;i++)
3022                                                 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomi*Emi+pj)*_vec_normal[i] + ui_n*rhomi*_Vi[2+i];
3023                                         _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=ui_n*(rhomi*(cp_v*cmi+cp_l*(1-cmi))-Hmi*rhomi*rhomi*(cp_v*cmi/(rho_vi*(h_vi-q_v))+cp_l*(1-cmi)/(rho_li*(h_li-q_l))));
3024
3025
3026                                         /******** Construction de la matrice A^- *********/
3027                                         _AroeMinusImplicit[0*_nVar+0]=0;
3028                                         _AroeMinusImplicit[0*_nVar+1]=0;
3029                                         for(int i=0;i<_Ndim;i++)
3030                                                 _AroeMinusImplicit[0*_nVar+2+i]=0;
3031                                         _AroeMinusImplicit[0*_nVar+2+_Ndim]=0;
3032                                         _AroeMinusImplicit[1*_nVar+0]=0;
3033                                         _AroeMinusImplicit[1*_nVar+1]=0;
3034                                         for(int i=0;i<_Ndim;i++)
3035                                                 _AroeMinusImplicit[1*_nVar+2+i]=0;
3036                                         _AroeMinusImplicit[1*_nVar+2+_Ndim]=0;
3037                                         for(int i=0;i<_Ndim;i++)
3038                                         {
3039                                                 _AroeMinusImplicit[(2+i)*_nVar+0]=0;
3040                                                 _AroeMinusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
3041                                                 for(int j=0;j<_Ndim;j++)
3042                                                         _AroeMinusImplicit[(2+i)*_nVar+2+j]=0;
3043                                                 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=0;
3044                                                 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=0;
3045                                         }
3046                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=0;
3047                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=ui_n;
3048                                         for(int i=0;i<_Ndim;i++)
3049                                                 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=0;
3050                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
3051                                 }
3052                                 else if(u_mn<-_precision)
3053                                 {
3054                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
3055                                                 cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Hmj= "<<Hmj<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
3056
3057                                         /***********Calcul des valeurs propres ********/
3058                                         vector<complex<double> > vp_dist(3,0);
3059                                         getMixturePressureDerivatives( mi_v, mi_l, pi, Tmi);
3060                                         if(_kappa*hmi+_khi+cmi*_ksi<0)
3061                                         {
3062                                                 *_runLogFile<<"staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe"<<endl;
3063                                                 throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
3064                                         }
3065                                         double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
3066                                         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
3067                                                 cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
3068
3069                                         //On remplit les valeurs propres
3070                                         vp_dist[0]=uj_n+ami;
3071                                         vp_dist[1]=uj_n-ami;
3072                                         vp_dist[2]=uj_n;
3073
3074                                         _maxvploc=fabs(uj_n)+ami;
3075                                         if(_maxvploc>_maxvp)
3076                                                 _maxvp=_maxvploc;
3077
3078                                         /******** Construction de la matrice A^+ *********/
3079                                         _AroePlusImplicit[0*_nVar+0]=0;
3080                                         _AroePlusImplicit[0*_nVar+1]=0;
3081                                         for(int i=0;i<_Ndim;i++)
3082                                                 _AroePlusImplicit[0*_nVar+2+i]=0;
3083                                         _AroePlusImplicit[0*_nVar+2+_Ndim]=0;
3084                                         _AroePlusImplicit[1*_nVar+0]=0;
3085                                         _AroePlusImplicit[1*_nVar+1]=0;
3086                                         for(int i=0;i<_Ndim;i++)
3087                                                 _AroePlusImplicit[1*_nVar+2+i]=0;
3088                                         _AroePlusImplicit[1*_nVar+2+_Ndim]=0;
3089                                         for(int i=0;i<_Ndim;i++)
3090                                         {
3091                                                 _AroePlusImplicit[(2+i)*_nVar+0]=0;
3092                                                 _AroePlusImplicit[(2+i)*_nVar+1]=_vec_normal[i];
3093                                                 for(int j=0;j<_Ndim;j++)
3094                                                         _AroePlusImplicit[(2+i)*_nVar+2+j]=0;
3095                                                 _AroePlusImplicit[(2+i)*_nVar+2+i]+=0;
3096                                                 _AroePlusImplicit[(2+i)*_nVar+2+_Ndim]=0;
3097                                         }
3098                                         _AroePlusImplicit[(2+_Ndim)*_nVar+0]=0;
3099                                         _AroePlusImplicit[(2+_Ndim)*_nVar+1]=uj_n;
3100                                         for(int i=0;i<_Ndim;i++)
3101                                                 _AroePlusImplicit[(2+_Ndim)*_nVar+2+i]=0;
3102                                         _AroePlusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=0;
3103
3104                                         /******** Construction de la matrice A^- *********/
3105                                         _AroeMinusImplicit[0*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
3106                                         _AroeMinusImplicit[0*_nVar+1]= rhomj*rhomj*uj_n*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)));
3107                                         for(int i=0;i<_Ndim;i++)
3108                                                 _AroeMinusImplicit[0*_nVar+2+i]=rhomj*_vec_normal[i];
3109                                         _AroeMinusImplicit[0*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l)));
3110                                         _AroeMinusImplicit[1*_nVar+0]=rhomj*uj_n-cmj*rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj);
3111                                         _AroeMinusImplicit[1*_nVar+1]=-cmj*rhomj*rhomj*uj_n*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)));
3112                                         for(int i=0;i<_Ndim;i++)
3113                                                 _AroeMinusImplicit[1*_nVar+2+i]=cmj*rhomj*_vec_normal[i];
3114                                         _AroeMinusImplicit[1*_nVar+2+_Ndim]=-cmj*rhomj*rhomj*uj_n*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l)));
3115                                         for(int i=0;i<_Ndim;i++)
3116                                         {
3117                                                 _AroeMinusImplicit[(2+i)*_nVar+0]=-rhomj*rhomj*uj_n*(1/rho_vj-1/rho_lj)*_Vj[2+i];
3118                                                 _AroeMinusImplicit[(2+i)*_nVar+1]=rhomj*rhomj*uj_n*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)))*_Vj[2+i];
3119                                                 for(int j=0;j<_Ndim;j++)
3120                                                         _AroeMinusImplicit[(2+i)*_nVar+2+j]=rhomj*_Vj[2+i]*_vec_normal[j];
3121                                                 _AroeMinusImplicit[(2+i)*_nVar+2+i]+=rhomj*uj_n;
3122                                                 _AroeMinusImplicit[(2+i)*_nVar+2+_Ndim]=-rhomj*rhomj*uj_n*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l)))*_Vj[2+i];
3123                                         }
3124                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+0]=uj_n*(rhomj*(h_vj-h_lj)-Hmj*rhomj*rhomj*(1/rho_vj-1/rho_lj));
3125                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+1]=rhomj*rhomj*Hmj*(gamma_v*cmj/(rho_vj*rho_vj*(gamma_v-1)*(h_vj-q_v))+gamma_l*(1-cmj)/(rho_lj*rho_lj*(gamma_l-1)*(h_lj-q_l)));
3126                                         for(int i=0;i<_Ndim;i++)
3127                                                 _AroeMinusImplicit[(2+_Ndim)*_nVar+2+i]=(rhomj*Emj+pi)*_vec_normal[i] + uj_n*rhomj*_Vj[2+i];
3128                                         _AroeMinusImplicit[(2+_Ndim)*_nVar+2+_Ndim]=uj_n*(rhomj*(cp_v*cmj+cp_l*(1-cmj))-Hmj*rhomj*rhomj*(cp_v*cmj/(rho_vj*(h_vj-q_v))+cp_l*(1-cmj)/(rho_lj*(h_lj-q_l))));
3129                                 }
3130                                 else
3131                                 {
3132                                         *_runLogFile<< "DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero" << endl;
3133                                         _runLogFile->close();
3134                                         throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: velocity umn should be non zero");
3135                                 }
3136                         }
3137                         else
3138                                 throw CdmathException("DriftModel::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
3139                 }
3140                 else//case nil velocity on the interface, multiply by jacobian matrix
3141                 {
3142                         primToConsJacobianMatrix(_Vj);
3143                         Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
3144                         primToConsJacobianMatrix(_Vi);
3145                         Polynoms::matrixProduct(_AroePlus,  _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
3146                 }
3147         }
3148 }
3149
3150 void DriftModel::applyVFRoeLowMachCorrections(bool isBord, string nameOfGroup)
3151 {
3152         if(_nonLinearFormulation!=VFRoe)
3153                 throw CdmathException("DriftModel::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
3154         else//_nonLinearFormulation==VFRoe
3155         {
3156                 if(_spaceScheme==lowMach){
3157                         double u_2=0;
3158                         for(int i=0;i<_Ndim;i++)
3159                                 u_2 += _Uroe[2+i]*_Uroe[2+i];
3160
3161                         double  c = _maxvploc;//mixture sound speed
3162                         double M=sqrt(u_2)/c;//Mach number
3163                         _Vij[1]=M*_Vij[1]+(1-M)*(_Vi[1]+_Vj[1])/2;
3164                 }
3165                 else if(_spaceScheme==pressureCorrection)
3166                 {
3167                         /*
3168                          * option 1 : no pressure correction except for inner walls where p^*=p_int
3169                          * option 2 : Clerc pressure correction everywhere, except for inner walls where p^*=p_int
3170                          * option 3 : Clerc pressure correction everywhere, even for inner walls
3171                          * option 4 : Clerc pressure correction only inside, upwind on wall and innerwall boundaries
3172                          * option 5 : Clerc pressure correction inside the domain and special pressure correction involving gravity at wall and inner wall boundaries
3173                          * */
3174                         bool isWall=false, isInnerWall=false;
3175                         if(isBord)
3176                         {
3177                                 isWall = (_limitField[nameOfGroup].bcType==Wall);
3178                                 isInnerWall = (_limitField[nameOfGroup].bcType==InnerWall);
3179                         }
3180
3181                         if(     (!isInnerWall && _pressureCorrectionOrder==2)  || _pressureCorrectionOrder==3 ||
3182                                 (!isBord &&     _pressureCorrectionOrder==4) || (!isWall && !isInnerWall && _pressureCorrectionOrder==5)   )//Clerc pressure correction
3183                         {
3184                                 double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;
3185                                 for(int i=0;i<_Ndim;i++)
3186                                 {
3187                                         norm_uij += _Uroe[2+i]*_Uroe[2+i];
3188                                         uij_n += _Uroe[2+i]*_vec_normal[i];
3189                                         ui_n += _Vi[2+i]*_vec_normal[i];
3190                                         uj_n += _Vj[2+i]*_vec_normal[i];
3191                                 }
3192                                 norm_uij=sqrt(norm_uij);
3193                                         if(norm_uij>_precision)//avoid division by zero
3194                                         _Vij[1]=(_Vi[1]+_Vj[1])/2 + uij_n/norm_uij*(_Vj[1]-_Vi[1])/4 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
3195                                 else
3196                                         _Vij[1]=(_Vi[1]+_Vj[1])/2                                    - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
3197                         }
3198                         else if((isWall || isInnerWall) && _pressureCorrectionOrder==5)//correction de pression gravitaire
3199                         {
3200                                 double g_n=0;//scalar product of gravity and normal vector
3201                                 for(int i=0;i<_Ndim;i++)
3202                                         g_n += _GravityField3d[i]*_vec_normal[i];
3203                                 _Vij[1]=_Vi[1]+ _Ui[0]*g_n/_inv_dxi/2;
3204                         }
3205                         else if(isInnerWall && (_pressureCorrectionOrder==1 || _pressureCorrectionOrder==2) )
3206                                         _Vij[1]=_Vi[1];
3207                 }
3208                 else if(_spaceScheme==staggered)
3209                 {
3210                         double uij_n=0;
3211                         for(int i=0;i<_Ndim;i++)
3212                                 uij_n += _Uroe[1+i]*_vec_normal[i];
3213
3214                         if(uij_n>_precision){
3215                                 _Vij[0]=_Vi[0];
3216                                 _Vij[1]=_Vj[1];
3217                                 for(int i=0;i<_Ndim;i++)
3218                                         _Vij[2+i]=_Vi[2+i];
3219                                 _Vij[_nVar-1]=_Vi[_nVar-1];
3220                         }
3221                         else if(uij_n<-_precision){
3222                                 _Vij[0]=_Vj[0];
3223                                 _Vij[1]=_Vi[1];
3224                                 for(int i=0;i<_Ndim;i++)
3225                                         _Vij[2+i]=_Vj[2+i];
3226                                 _Vij[_nVar-1]=_Vj[_nVar-1];
3227                         }
3228                         else{
3229                                 _Vij[0]=(_Vi[0]+_Vj[0])/2;
3230                                 _Vij[1]=(_Vi[1]+_Vj[1])/2;
3231                                 for(int i=0;i<_Ndim;i++)
3232                                         _Vij[2+i]=(_Vi[2+i]+_Vj[2+i])/2;
3233                                 _Vij[_nVar-1]=(_Vi[_nVar-1]+_Vj[_nVar-1])/2;
3234                         }
3235                 }
3236                 primToCons(_Vij,0,_Uij,0);
3237         }
3238 }
3239 void DriftModel::RoeMatrixConservativeVariables(double cm, double umn,double ecin, double Hm,Vector vitesse)
3240 {
3241         //On remplit la matrice de Roe en variables conservatives : F(u_L)-F(U_R)=Aroe (U_L-U_R)
3242         _Aroe[0*_nVar+0]=0;
3243         _Aroe[0*_nVar+1]=0;
3244         for(int i=0;i<_Ndim;i++)
3245                 _Aroe[0*_nVar+2+i]=_vec_normal[i];
3246         _Aroe[0*_nVar+2+_Ndim]=0;
3247         _Aroe[1*_nVar+0]=-umn*cm;
3248         _Aroe[1*_nVar+1]=umn;
3249         for(int i=0;i<_Ndim;i++)
3250                 _Aroe[1*_nVar+2+i]=cm*_vec_normal[i];
3251         _Aroe[1*_nVar+2+_Ndim]=0;
3252         for(int i=0;i<_Ndim;i++)
3253         {
3254                 _Aroe[(2+i)*_nVar+0]=(_khi+_kappa*ecin)*_vec_normal[i]-umn*vitesse[i];
3255                 _Aroe[(2+i)*_nVar+1]=_ksi*_vec_normal[i];
3256                 for(int j=0;j<_Ndim;j++)
3257                         _Aroe[(2+i)*_nVar+2+j]=vitesse[i]*_vec_normal[j]-_kappa*_vec_normal[i]*vitesse[j];
3258                 _Aroe[(2+i)*_nVar+2+i]+=umn;
3259                 _Aroe[(2+i)*_nVar+2+_Ndim]=_kappa*_vec_normal[i];
3260         }
3261         _Aroe[(2+_Ndim)*_nVar+0]=(_khi+_kappa*ecin-Hm)*umn;
3262         _Aroe[(2+_Ndim)*_nVar+1]=_ksi*umn;
3263         for(int i=0;i<_Ndim;i++)
3264                 _Aroe[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i]-_kappa*umn*vitesse[i];
3265         _Aroe[(2+_Ndim)*_nVar+2+_Ndim]=(_kappa+1)*umn;
3266 }
3267
3268 void DriftModel::staggeredRoeUpwindingMatrixConservativeVariables(double cm, double umn,double ecin, double Hm, Vector vitesse)
3269 {
3270         //Calcul de décentrement de type décalé pour formulation Roe
3271         if(abs(umn)>_precision)//non zero velocity at the interface
3272         {
3273                 _absAroe[0*_nVar+0]=0;
3274                 _absAroe[0*_nVar+1]=0;
3275                 for(int i=0;i<_Ndim;i++)
3276                         _absAroe[0*_nVar+2+i]=_vec_normal[i];
3277                 _absAroe[0*_nVar+2+_Ndim]=0;
3278                 _absAroe[1*_nVar+0]=-umn*cm;
3279                 _absAroe[1*_nVar+1]=umn;
3280                 for(int i=0;i<_Ndim;i++)
3281                         _absAroe[1*_nVar+2+i]=cm*_vec_normal[i];
3282                 _absAroe[1*_nVar+2+_Ndim]=0;
3283                 for(int i=0;i<_Ndim;i++)
3284                 {
3285                         _absAroe[(2+i)*_nVar+0]=-(_khi+_kappa*ecin)*_vec_normal[i]-umn*vitesse[i];
3286                         _absAroe[(2+i)*_nVar+1]=-_ksi*_vec_normal[i];
3287                         for(int j=0;j<_Ndim;j++)
3288                                 _absAroe[(2+i)*_nVar+2+j]=vitesse[i]*_vec_normal[j]+_kappa*_vec_normal[i]*vitesse[j];
3289                         _absAroe[(2+i)*_nVar+2+i]+=umn;
3290                         _absAroe[(2+i)*_nVar+2+_Ndim]=-_kappa*_vec_normal[i];
3291                 }
3292                 _absAroe[(2+_Ndim)*_nVar+0]=(-_khi-_kappa*ecin-Hm)*umn;
3293                 _absAroe[(2+_Ndim)*_nVar+1]=-_ksi*umn;
3294                 for(int i=0;i<_Ndim;i++)
3295                         _absAroe[(2+_Ndim)*_nVar+2+i]=Hm*_vec_normal[i]+_kappa*umn*vitesse[i];
3296                 _absAroe[(2+_Ndim)*_nVar+2+_Ndim]=(-_kappa+1)*umn;
3297
3298                 double signu;
3299                 if(umn>_precision)
3300                         signu=1;
3301                 else // umn<-_precision
3302                         signu=-1;
3303
3304                 for(int i=0; i<_nVar*_nVar;i++)
3305                         _absAroe[i] *= signu;
3306         }
3307         else//umn=0>centered scheme
3308         {
3309                 for(int i=0; i<_nVar*_nVar;i++)
3310                         _absAroe[i] =0;
3311                 //for(int i=0; i<_nVar;i++)
3312                 //      _absAroe[i+_nVar*i] =_maxvploc;
3313         }
3314 }
3315
3316 void DriftModel::staggeredRoeUpwindingMatrixPrimitiveVariables(double concentration, double rhom, double umn, double Hm, Vector vitesse)
3317 {
3318         //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
3319         //Calcul de décentrement de type décalé pour formulation Roe
3320         _absAroeImplicit[0*_nVar+0]= _drho_sur_dcv*umn;
3321         _absAroeImplicit[0*_nVar+1]= _drho_sur_dp *umn;
3322         for(int i=0;i<_Ndim;i++)
3323                 _absAroeImplicit[0*_nVar+2+i]=rhom*_vec_normal[i];
3324         _absAroeImplicit[0*_nVar+2+_Ndim]=_drho_sur_dT*umn;
3325         _absAroeImplicit[1*_nVar+0]= _drhocv_sur_dcv*umn;
3326         _absAroeImplicit[1*_nVar+1]= _drhocv_sur_dp *umn;
3327         for(int i=0;i<_Ndim;i++)
3328                 _absAroeImplicit[1*_nVar+2+i]=rhom*concentration*_vec_normal[i];
3329         _absAroeImplicit[1*_nVar+2+_Ndim]=_drhocv_sur_dT*umn;
3330         for(int i=0;i<_Ndim;i++)
3331         {
3332                 _absAroeImplicit[(2+i)*_nVar+0]= _drho_sur_dcv*umn*vitesse[i];
3333                 _absAroeImplicit[(2+i)*_nVar+1]= _drho_sur_dp *umn*vitesse[i]-_vec_normal[i];
3334                 for(int j=0;j<_Ndim;j++)
3335                         _absAroeImplicit[(2+i)*_nVar+2+j]=rhom*vitesse[i]*_vec_normal[j];
3336                 _absAroeImplicit[(2+i)*_nVar+2+i]+=rhom*umn;
3337                 _absAroeImplicit[(2+i)*_nVar+2+_Ndim]=_drho_sur_dT*umn*vitesse[i];
3338         }
3339         _absAroeImplicit[(2+_Ndim)*_nVar+0]=  _drhoE_sur_dcv  *umn;
3340         _absAroeImplicit[(2+_Ndim)*_nVar+1]=( _drhoE_sur_dp+1)*umn;
3341         for(int i=0;i<_Ndim;i++)
3342                 _absAroeImplicit[(2+_Ndim)*_nVar+2+i]=rhom*(Hm*_vec_normal[i]+umn*vitesse[i]);
3343         _absAroeImplicit[(2+_Ndim)*_nVar+2+_Ndim]=_drhoE_sur_dT*umn;
3344 }
3345
3346 void DriftModel::convectionMatrixPrimitiveVariables(double concentration, double rhom, double umn, double Hm, Vector vitesse)
3347 {
3348         //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
3349         //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
3350         //EOS is more involved with primitive variables
3351         //first call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
3352         _AroeImplicit[0*_nVar+0]=_drho_sur_dcv*umn;
3353         _AroeImplicit[0*_nVar+1]=_drho_sur_dp*umn;
3354         for(int i=0;i<_Ndim;i++)
3355                 _AroeImplicit[0*_nVar+2+i]=rhom*_vec_normal[i];
3356         _AroeImplicit[0*_nVar+2+_Ndim]=_drho_sur_dT*umn;
3357         _AroeImplicit[1*_nVar+0]=_drhocv_sur_dcv*umn;
3358         _AroeImplicit[1*_nVar+1]=_drhocv_sur_dp*umn;
3359         for(int i=0;i<_Ndim;i++)
3360                 _AroeImplicit[1*_nVar+2+i]=rhom*concentration*_vec_normal[i];
3361         _AroeImplicit[1*_nVar+2+_Ndim]=_drhocv_sur_dT*umn;
3362         for(int i=0;i<_Ndim;i++)
3363         {
3364                 _AroeImplicit[(2+i)*_nVar+0]=_drho_sur_dcv*umn*vitesse[i];
3365                 _AroeImplicit[(2+i)*_nVar+1]=_drho_sur_dp *umn*vitesse[i]+_vec_normal[i];
3366                 for(int j=0;j<_Ndim;j++)
3367                         _AroeImplicit[(2+i)*_nVar+2+j]=rhom*vitesse[i]*_vec_normal[j];
3368                 _AroeImplicit[(2+i)*_nVar+2+i]+=rhom*umn;
3369                 _AroeImplicit[(2+i)*_nVar+2+_Ndim]=_drho_sur_dT*umn*vitesse[i];
3370         }
3371         _AroeImplicit[(2+_Ndim)*_nVar+0]= _drhoE_sur_dcv  *umn;
3372         _AroeImplicit[(2+_Ndim)*_nVar+1]=(_drhoE_sur_dp+1)*umn;
3373         for(int i=0;i<_Ndim;i++)
3374                 _AroeImplicit[(2+_Ndim)*_nVar+2+i]=rhom*(Hm*_vec_normal[i]+umn*vitesse[i]);
3375         _AroeImplicit[(2+_Ndim)*_nVar+2+_Ndim]=_drhoE_sur_dT*umn;
3376 }
3377
3378 void DriftModel::getDensityDerivatives(double concentration, double pression, double temperature,double v2)
3379 {
3380         //EOS is more involved with primitive variables
3381
3382         double rho_v=_fluides[0]->getDensity(pression,temperature);
3383         double rho_l=_fluides[1]->getDensity(pression,temperature);
3384         double gamma_v=_fluides[0]->constante("gamma");
3385         double gamma_l=_fluides[1]->constante("gamma");
3386         double q_v=_fluides[0]->constante("q");
3387         double q_l=_fluides[1]->constante("q");
3388
3389         double rho=concentration*rho_v+(1-concentration)*rho_l;;
3390
3391         if(     !_useDellacherieEOS)
3392         {
3393                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
3394                 StiffenedGas* fluide1=dynamic_cast<StiffenedGas*>(_fluides[1]);
3395                 double e_v = fluide0->getInternalEnergy(temperature);
3396                 double e_l = fluide1->getInternalEnergy(temperature);
3397                 double cv_v=fluide0->constante("cv");
3398                 double cv_l=fluide1->constante("cv");
3399                 double e=concentration*e_v+(1-concentration)*e_l;
3400                 double E=e+0.5*v2;
3401
3402                 _drho_sur_dcv=-rho*rho*(1/rho_v-1/rho_l);
3403                 _drho_sur_dp=rho*rho*(
3404                                 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
3405                                 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
3406                 );
3407                 _drho_sur_dT=-rho*rho*(
3408                                 cv_v*   concentration /(rho_v*(e_v-q_v))
3409                                 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
3410                 );
3411
3412                 _drhocv_sur_dcv=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
3413                 _drhocv_sur_dp=concentration*rho*rho*(
3414                                 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
3415                                 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
3416                 );
3417                 _drhocv_sur_dT=-concentration*rho*rho*(
3418                                 cv_v*   concentration /(rho_v*(e_v-q_v))
3419                                 +cv_l*(1-concentration)/(rho_l*(e_l-q_l))
3420                 );
3421                 _drhoE_sur_dcv=rho*(e_v-e_l)-E*rho*rho*(1/rho_v-1/rho_l);
3422                 _drhoE_sur_dp=E*rho*rho*(
3423                                 concentration /(rho_v*rho_v*(gamma_v-1)*(e_v-q_v))
3424                                 +(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(e_l-q_l))
3425                 );
3426                 _drhoE_sur_dT=rho*(cv_v*concentration + cv_l*(1-concentration))
3427                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         -rho*rho*E*( cv_v*   concentration /(rho_v*(e_v-q_v))
3428                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         +cv_l*(1-concentration)/(rho_l*(e_l-q_l)));
3429         }
3430         else if(_useDellacherieEOS)
3431         {
3432                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
3433                 StiffenedGasDellacherie* fluide1=dynamic_cast<StiffenedGasDellacherie*>(_fluides[1]);
3434                 double h_v=fluide0->getEnthalpy(temperature);
3435                 double h_l=fluide1->getEnthalpy(temperature);
3436                 double h=concentration*h_v+(1-concentration)*h_l;
3437                 double H=h+0.5*v2;
3438                 double cp_v=fluide0->constante("cp");
3439                 double cp_l=fluide1->constante("cp");
3440
3441                 _drho_sur_dcv=-rho*rho*(1/rho_v-1/rho_l);
3442                 _drho_sur_dp =rho*rho*(
3443                                 gamma_v*   concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
3444                                 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
3445                 );
3446                 _drho_sur_dT=-rho*rho*(
3447                                 cp_v*   concentration /(rho_v*(h_v-q_v))
3448                                 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
3449                 );
3450
3451                 _drhocv_sur_dcv=rho-concentration*rho*rho*(1/rho_v-1/rho_l);
3452                 _drhocv_sur_dp=    concentration*rho*rho*(
3453                                 gamma_v*   concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
3454                                 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
3455                 );
3456                 _drhocv_sur_dT=-concentration*rho*rho*(
3457                                 cp_v*   concentration /(rho_v*(h_v-q_v))
3458                                 +cp_l*(1-concentration)/(rho_l*(h_l-q_l))
3459                 );
3460                 _drhoE_sur_dcv=rho*(h_v-h_l)-H*rho*rho*(1/rho_v-1/rho_l);
3461                 _drhoE_sur_dp=H*rho*rho*(
3462                                 gamma_v*   concentration /(rho_v*rho_v*(gamma_v-1)*(h_v-q_v))
3463                                 +gamma_l*(1-concentration)/(rho_l*rho_l*(gamma_l-1)*(h_l-q_l))
3464                 )-1;
3465                 _drhoE_sur_dT=rho*(cp_v*concentration + cp_l*(1-concentration))
3466                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    -rho*rho*H*( cp_v*   concentration /(rho_v*(h_v-q_v))
3467                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    +cp_l*(1-concentration)/(rho_l*(h_l-q_l)));
3468         }
3469         else
3470                 throw CdmathException("DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
3471
3472         if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
3473         {
3474                 cout<<"_drho_sur_dcv= "<<_drho_sur_dcv<<", _drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
3475                 cout<<"_drhocv_sur_dcv= "<<_drhocv_sur_dcv<<", _drhocv_sur_dp= "<<_drhocv_sur_dp<<", _drhocv_sur_dT= "<<_drhocv_sur_dT<<endl;
3476                 cout<<"_drhoE_sur_dcv= "<<_drhoE_sur_dcv<<", _drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
3477         }
3478 }
3479
3480 void DriftModel::save(){
3481     PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results at time step number %d \n\n", _nbTimeStep);
3482     *_runLogFile<< "Saving numerical results at time step number "<< _nbTimeStep << endl<<endl;
3483
3484         string prim(_path+"/DriftModelPrim_");
3485         string cons(_path+"/DriftModelCons_");
3486         string allFields(_path+"/");
3487         prim+=_fileName;
3488         cons+=_fileName;
3489         allFields+=_fileName;
3490
3491         PetscInt Ii;
3492         for (long i = 0; i < _Nmailles; i++){
3493                 // j = 0 : concentration, j=1 : pressure; j = _nVar - 1: temperature; j = 2,..,_nVar-2: velocity
3494                 for (int j = 0; j < _nVar; j++){
3495                         Ii = i*_nVar +j;
3496                         VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
3497                 }
3498         }
3499         if(_saveConservativeField){
3500                 for (long i = 0; i < _Nmailles; i++){
3501                         // j = 0 : total density; j = 1 : vapour density; j = _nVar - 1 : energy j = 2,..,_nVar-2: momentum
3502                         for (int j = 0; j < _nVar; j++){
3503                                 Ii = i*_nVar +j;
3504                                 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
3505                         }
3506                 }
3507                 _UU.setTime(_time,_nbTimeStep);
3508         }
3509         _VV.setTime(_time,_nbTimeStep);
3510         // create mesh and component info
3511         if (_nbTimeStep ==0 || _restartWithNewFileName){
3512                 if (_restartWithNewFileName)
3513                         _restartWithNewFileName=false;
3514                 string suppress_previous_runs ="rm -rf *"+_fileName+"_*";
3515                 system(suppress_previous_runs.c_str());//Nettoyage des précédents calculs identiques
3516
3517                 if(_saveConservativeField){
3518                         _UU.setInfoOnComponent(0,"Total_Density");// (kg/m^3)
3519
3520                         _UU.setInfoOnComponent(1,"Partial_Density");// (kg/m^3)
3521                         _UU.setInfoOnComponent(2,"Momentum_x");// (kg/m^2/s)
3522                         if (_Ndim>1)
3523                                 _UU.setInfoOnComponent(3,"Momentum_y");// (kg/m^2/s)
3524                         if (_Ndim>2)
3525                                 _UU.setInfoOnComponent(4,"Momentum_z");// (kg/m^2/s)
3526
3527                         _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
3528                         switch(_saveFormat)
3529                         {
3530                         case VTK :
3531                                 _UU.writeVTK(cons);
3532                                 break;
3533                         case MED :
3534                                 _UU.writeMED(cons);
3535                                 break;
3536                         case CSV :
3537                                 _UU.writeCSV(cons);
3538                                 break;
3539                         }
3540                 }
3541
3542                 _VV.setInfoOnComponent(0,"Concentration");
3543                 _VV.setInfoOnComponent(1,"Pressure_(Pa)");
3544                 _VV.setInfoOnComponent(2,"Velocity_x_(m/s)");
3545                 if (_Ndim>1)
3546                         _VV.setInfoOnComponent(3,"Velocity_y_(m/s)");
3547                 if (_Ndim>2)
3548                         _VV.setInfoOnComponent(4,"Velocity_z_(m/s)");
3549                 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
3550                 switch(_saveFormat)
3551                 {
3552                 case VTK :
3553                         _VV.writeVTK(prim);
3554                         break;
3555                 case MED :
3556                         _VV.writeMED(prim);
3557                         break;
3558                 case CSV :
3559                         _VV.writeCSV(prim);
3560                         break;
3561                 }
3562         }
3563         // do not create mesh
3564         else{
3565                 switch(_saveFormat)
3566                 {
3567                 case VTK :
3568                         _VV.writeVTK(prim,false);
3569                         break;
3570                 case MED :
3571                         _VV.writeMED(prim,false);
3572                         break;
3573                 case CSV :
3574                         _VV.writeCSV(prim);
3575                         break;
3576                 }
3577                 if(_saveConservativeField){
3578                         switch(_saveFormat)
3579                         {
3580                         case VTK :
3581                                 _UU.writeVTK(cons,false);
3582                                 break;
3583                         case MED :
3584                                 _UU.writeMED(cons,false);
3585                                 break;
3586                         case CSV :
3587                                 _UU.writeCSV(cons);
3588                                 break;
3589                         }
3590                 }
3591         }
3592         if(_saveVelocity){
3593                 for (long i = 0; i < _Nmailles; i++){
3594                         // j = 0 : concentration, j=1 : pressure; j = _nVar - 1: temperature; j = 2,..,_nVar-2: velocity
3595                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
3596                         {
3597                                 int Ii = i*_nVar +2+j;
3598                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
3599                         }
3600                         for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
3601                                 _Vitesse(i,j)=0;
3602                 }
3603                 _Vitesse.setTime(_time,_nbTimeStep);
3604                 if (_nbTimeStep ==0 || _restartWithNewFileName){
3605                         _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
3606                         _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
3607                         _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
3608                         switch(_saveFormat)
3609                         {
3610                         case VTK :
3611                                 _Vitesse.writeVTK(prim+"_Velocity");
3612                                 break;
3613                         case MED :
3614                                 _Vitesse.writeMED(prim+"_Velocity");
3615                                 break;
3616                         case CSV :
3617                                 _Vitesse.writeCSV(prim+"_Velocity");
3618                                 break;
3619                         }
3620                 }
3621                 else{
3622                         switch(_saveFormat)
3623                         {
3624                         case VTK :
3625                                 _Vitesse.writeVTK(prim+"_Velocity",false);
3626                                 break;
3627                         case MED :
3628                                 _Vitesse.writeMED(prim+"_Velocity",false);
3629                                 break;
3630                         case CSV :
3631                                 _Vitesse.writeCSV(prim+"_Velocity");
3632                                 break;
3633                         }
3634                 }
3635         }
3636         if(_saveAllFields)
3637         {
3638                 double p,Tm,cv,alpha_v,rhom,rho_v,rho_l, m_v, m_l, h_v, h_l, vx,vy,vz;
3639                 int Ii;
3640                 for (long i = 0; i < _Nmailles; i++){
3641                         Ii = i*_nVar;
3642                         VecGetValues(_conservativeVars,1,&Ii,&rhom);
3643                         Ii = i*_nVar;
3644                         VecGetValues(_primitiveVars,1,&Ii,&cv);
3645                         Ii = i*_nVar +1;
3646                         VecGetValues(_primitiveVars,1,&Ii,&p);
3647                         Ii = i*_nVar +_nVar-1;
3648                         VecGetValues(_primitiveVars,1,&Ii,&Tm);
3649                         Ii = i*_nVar +2;
3650                         VecGetValues(_primitiveVars,1,&Ii,&vx);
3651                         if(_Ndim>1)
3652                         {
3653                                 Ii = i*_nVar +3;
3654                                 VecGetValues(_primitiveVars,1,&Ii,&vy);
3655                                 if(_Ndim>2){
3656                                         Ii = i*_nVar +4;
3657                                         VecGetValues(_primitiveVars,1,&Ii,&vz);
3658                                 }
3659                         }
3660
3661                         rho_v=_fluides[0]->getDensity(p,Tm);
3662                         rho_l=_fluides[1]->getDensity(p,Tm);
3663                         alpha_v=cv*rhom/rho_v;
3664                         m_v=cv*rhom;
3665                         m_l=(1-cv)*rhom;
3666                         h_v=_fluides[0]->getEnthalpy(Tm,rho_v);
3667                         h_l=_fluides[1]->getEnthalpy(Tm,rho_l);
3668
3669                         _VoidFraction(i)=alpha_v;
3670                         _Enthalpy(i)=(m_v*h_v+m_l*h_l)/rhom;
3671                         _Concentration(i)=cv;
3672                         _mixtureDensity(i)=rhom;
3673                         _Pressure(i)=p;
3674                         _Temperature(i)=Tm;
3675                         _DensiteLiquide(i)=rho_l;
3676                         _DensiteVapeur(i)=rho_v;
3677                         _EnthalpieLiquide(i)=h_l;
3678                         _EnthalpieVapeur(i)=h_v;
3679                         _VitesseX(i)=vx;
3680                         if(_Ndim>1)
3681                         {
3682                                 _VitesseY(i)=vy;
3683                                 if(_Ndim>2)
3684                                         _VitesseZ(i)=vz;
3685                         }
3686                 }
3687                 _VoidFraction.setTime(_time,_nbTimeStep);
3688                 _Enthalpy.setTime(_time,_nbTimeStep);
3689                 _Concentration.setTime(_time,_nbTimeStep);
3690                 _mixtureDensity.setTime(_time,_nbTimeStep);
3691                 _Pressure.setTime(_time,_nbTimeStep);
3692                 _Temperature.setTime(_time,_nbTimeStep);
3693                 _DensiteLiquide.setTime(_time,_nbTimeStep);
3694                 _DensiteVapeur.setTime(_time,_nbTimeStep);
3695                 _EnthalpieLiquide.setTime(_time,_nbTimeStep);
3696                 _EnthalpieVapeur.setTime(_time,_nbTimeStep);
3697                 _VitesseX.setTime(_time,_nbTimeStep);
3698                 if(_Ndim>1)
3699                 {
3700                         _VitesseY.setTime(_time,_nbTimeStep);
3701                         if(_Ndim>2)
3702                                 _VitesseZ.setTime(_time,_nbTimeStep);
3703                 }
3704                 if (_nbTimeStep ==0 || _restartWithNewFileName){
3705                         switch(_saveFormat)
3706                         {
3707                         case VTK :
3708                                 _VoidFraction.writeVTK(allFields+"_VoidFraction");
3709                                 _Enthalpy.writeVTK(allFields+"_Enthalpy");
3710                                 _Concentration.writeVTK(allFields+"_Concentration");
3711                                 _mixtureDensity.writeVTK(allFields+"_Density");
3712                                 _Pressure.writeVTK(allFields+"_Pressure");
3713                                 _Temperature.writeVTK(allFields+"_Temperature");
3714                                 _DensiteLiquide.writeVTK(allFields+"_LiquidDensity");
3715                                 _DensiteVapeur.writeVTK(allFields+"_SteamDensity");
3716                                 _EnthalpieLiquide.writeVTK(allFields+"_LiquidEnthalpy");
3717                                 _EnthalpieVapeur.writeVTK(allFields+"_SteamEnthalpy");
3718                                 _VitesseX.writeVTK(allFields+"_VelocityX");
3719                                 if(_Ndim>1)
3720                                 {
3721                                         _VitesseY.writeVTK(allFields+"_VelocityY");
3722                                         if(_Ndim>2)
3723                                                 _VitesseZ.writeVTK(allFields+"_VelocityZ");
3724                                 }
3725                                 break;
3726                         case MED :
3727                                 _VoidFraction.writeMED(allFields+"_VoidFraction");
3728                                 _Enthalpy.writeMED(allFields+"_Enthalpy");
3729                                 _Concentration.writeMED(allFields+"_Concentration");
3730                                 _mixtureDensity.writeMED(allFields+"_Density");
3731                                 _Pressure.writeMED(allFields+"_Pressure");
3732                                 _Temperature.writeMED(allFields+"_Temperature");
3733                                 _DensiteLiquide.writeMED(allFields+"_LiquidDensity");
3734                                 _DensiteVapeur.writeMED(allFields+"_SteamDensity");
3735                                 _EnthalpieLiquide.writeMED(allFields+"_LiquidEnthalpy");
3736                                 _EnthalpieVapeur.writeMED(allFields+"_SteamEnthalpy");
3737                                 _VitesseX.writeMED(allFields+"_VelocityX");
3738                                 if(_Ndim>1)
3739                                 {
3740                                         _VitesseY.writeMED(allFields+"_VelocityY");
3741                                         if(_Ndim>2)
3742                                                 _VitesseZ.writeMED(allFields+"_VelocityZ");
3743                                 }
3744                                 break;
3745                         case CSV :
3746                                 _VoidFraction.writeCSV(allFields+"_VoidFraction");
3747                                 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3748                                 _Concentration.writeCSV(allFields+"_Concentration");
3749                                 _mixtureDensity.writeCSV(allFields+"_Density");
3750                                 _Pressure.writeCSV(allFields+"_Pressure");
3751                                 _Temperature.writeCSV(allFields+"_Temperature");
3752                                 _DensiteLiquide.writeCSV(allFields+"_LiquidDensity");
3753                                 _DensiteVapeur.writeCSV(allFields+"_SteamDensity");
3754                                 _EnthalpieLiquide.writeCSV(allFields+"_LiquidEnthalpy");
3755                                 _EnthalpieVapeur.writeCSV(allFields+"_SteamEnthalpy");
3756                                 _VitesseX.writeCSV(allFields+"_VelocityX");
3757                                 if(_Ndim>1)
3758                                 {
3759                                         _VitesseY.writeCSV(allFields+"_VelocityY");
3760                                         if(_Ndim>2)
3761                                                 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3762                                 }
3763                                 break;
3764                         }
3765                 }
3766                 else{
3767                         switch(_saveFormat)
3768                         {
3769                         case VTK :
3770                                 _VoidFraction.writeVTK(allFields+"_VoidFraction",false);
3771                                 _Enthalpy.writeVTK(allFields+"_Enthalpy",false);
3772                                 _Concentration.writeVTK(allFields+"_Concentration",false);
3773                                 _mixtureDensity.writeVTK(allFields+"_Density",false);
3774                                 _Pressure.writeVTK(allFields+"_Pressure",false);
3775                                 _Temperature.writeVTK(allFields+"_Temperature",false);
3776                                 _DensiteLiquide.writeVTK(allFields+"_LiquidDensity",false);
3777                                 _DensiteVapeur.writeVTK(allFields+"_SteamDensity",false);
3778                                 _EnthalpieLiquide.writeVTK(allFields+"_LiquidEnthalpy",false);
3779                                 _EnthalpieVapeur.writeVTK(allFields+"_SteamEnthalpy",false);
3780                                 _VitesseX.writeVTK(allFields+"_VelocityX",false);
3781                                 if(_Ndim>1)
3782                                 {
3783                                         _VitesseY.writeVTK(allFields+"_VelocityY",false);
3784                                         if(_Ndim>2)
3785                                                 _VitesseZ.writeVTK(allFields+"_VelocityZ",false);
3786                                 }
3787                                 break;
3788                         case MED :
3789                                 _VoidFraction.writeMED(allFields+"_VoidFraction",false);
3790                                 _Enthalpy.writeMED(allFields+"_Enthalpy",false);
3791                                 _Concentration.writeMED(allFields+"_Concentration",false);
3792                                 _mixtureDensity.writeMED(allFields+"_Density",false);
3793                                 _Pressure.writeMED(allFields+"_Pressure",false);
3794                                 _Temperature.writeMED(allFields+"_Temperature",false);
3795                                 _DensiteLiquide.writeMED(allFields+"_LiquidDensity",false);
3796                                 _DensiteVapeur.writeMED(allFields+"_SteamDensity",false);
3797                                 _EnthalpieLiquide.writeMED(allFields+"_LiquidEnthalpy",false);
3798                                 _EnthalpieVapeur.writeMED(allFields+"_SteamEnthalpy",false);
3799                                 _VitesseX.writeMED(allFields+"_VelocityX",false);
3800                                 if(_Ndim>1)
3801                                 {
3802                                         _VitesseY.writeMED(allFields+"_VelocityY",false);
3803                                         if(_Ndim>2)
3804                                                 _VitesseZ.writeMED(allFields+"_VelocityZ",false);
3805                                 }
3806                                 break;
3807                         case CSV :
3808                                 _VoidFraction.writeCSV(allFields+"_VoidFraction");
3809                                 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3810                                 _Concentration.writeCSV(allFields+"_Concentration");
3811                                 _mixtureDensity.writeCSV(allFields+"_Density");
3812                                 _Pressure.writeCSV(allFields+"_Pressure");
3813                                 _Temperature.writeCSV(allFields+"_Temperature");
3814                                 _DensiteLiquide.writeCSV(allFields+"_LiquidDensity");
3815                                 _DensiteVapeur.writeCSV(allFields+"_SteamDensity");
3816                                 _EnthalpieLiquide.writeCSV(allFields+"_LiquidEnthalpy");
3817                                 _EnthalpieVapeur.writeCSV(allFields+"_SteamEnthalpy");
3818                                 _VitesseX.writeCSV(allFields+"_VelocityX");
3819                                 if(_Ndim>1)
3820                                 {
3821                                         _VitesseY.writeCSV(allFields+"_VelocityY");
3822                                         if(_Ndim>2)
3823                                                 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3824                                 }
3825                                 break;
3826                         }
3827                 }
3828         }
3829         if(_isStationary)
3830         {
3831                 prim+="_Stat";
3832                 cons+="_Stat";
3833
3834                 switch(_saveFormat)
3835                 {
3836                 case VTK :
3837                         _VV.writeVTK(prim);
3838                         break;
3839                 case MED :
3840                         _VV.writeMED(prim);
3841                         break;
3842                 case CSV :
3843                         _VV.writeCSV(prim);
3844                         break;
3845                 }
3846
3847                 if(_saveConservativeField){
3848                         switch(_saveFormat)
3849                         {
3850                         case VTK :
3851                                 _UU.writeVTK(cons);
3852                                 break;
3853                         case MED :
3854                                 _UU.writeMED(cons);
3855                                 break;
3856                         case CSV :
3857                                 _UU.writeCSV(cons);
3858                                 break;
3859                         }
3860                 }
3861
3862                 if(_saveVelocity){
3863                         switch(_saveFormat)
3864                         {
3865                         case VTK :
3866                                 _Vitesse.writeVTK(prim+"_Velocity");
3867                                 break;
3868                         case MED :
3869                                 _Vitesse.writeMED(prim+"_Velocity");
3870                                 break;
3871                         case CSV :
3872                                 _Vitesse.writeCSV(prim+"_Velocity");
3873                                 break;
3874                         }
3875                 }
3876
3877                 if(_saveAllFields)
3878                 {
3879                         allFields+="_Stat";
3880                         switch(_saveFormat)
3881                         {
3882                         case VTK :
3883                                 _VoidFraction.writeVTK(allFields+"_VoidFraction");
3884                                 _Enthalpy.writeVTK(allFields+"_Enthalpy");
3885                                 _Concentration.writeVTK(allFields+"_Concentration");
3886                                 _mixtureDensity.writeVTK(allFields+"_Density");
3887                                 _Pressure.writeVTK(allFields+"_Pressure");
3888                                 _Temperature.writeVTK(allFields+"_Temperature");
3889                                 _DensiteLiquide.writeVTK(allFields+"_LiquidDensity");
3890                                 _DensiteVapeur.writeVTK(allFields+"_SteamDensity");
3891                                 _EnthalpieLiquide.writeVTK(allFields+"_LiquidEnthalpy");
3892                                 _EnthalpieVapeur.writeVTK(allFields+"_SteamEnthalpy");
3893                                 _VitesseX.writeVTK(allFields+"_VelocityX");
3894                                 if(_Ndim>1)
3895                                 {
3896                                         _VitesseY.writeVTK(allFields+"_VelocityY");
3897                                         if(_Ndim>2)
3898                                                 _VitesseZ.writeVTK(allFields+"_VelocityZ");
3899                                 }
3900                                 break;
3901                         case MED :
3902                                 _VoidFraction.writeMED(allFields+"_VoidFraction");
3903                                 _Enthalpy.writeMED(allFields+"_Enthalpy");
3904                                 _Concentration.writeMED(allFields+"_Concentration");
3905                                 _mixtureDensity.writeMED(allFields+"_Density");
3906                                 _Pressure.writeMED(allFields+"_Pressure");
3907                                 _Temperature.writeMED(allFields+"_Temperature");
3908                                 _DensiteLiquide.writeMED(allFields+"_LiquidDensity");
3909                                 _DensiteVapeur.writeMED(allFields+"_SteamDensity");
3910                                 _EnthalpieLiquide.writeMED(allFields+"_LiquidEnthalpy");
3911                                 _EnthalpieVapeur.writeMED(allFields+"_SteamEnthalpy");
3912                                 _VitesseX.writeMED(allFields+"_VelocityX");
3913                                 if(_Ndim>1)
3914                                 {
3915                                         _VitesseY.writeMED(allFields+"_VelocityY");
3916                                         if(_Ndim>2)
3917                                                 _VitesseZ.writeMED(allFields+"_VelocityZ");
3918                                 }
3919                                 break;
3920                         case CSV :
3921                                 _VoidFraction.writeCSV(allFields+"_VoidFraction");
3922                                 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3923                                 _Concentration.writeCSV(allFields+"_Concentration");
3924                                 _mixtureDensity.writeCSV(allFields+"_Density");
3925                                 _Pressure.writeCSV(allFields+"_Pressure");
3926                                 _Temperature.writeCSV(allFields+"_Temperature");
3927                                 _DensiteLiquide.writeCSV(allFields+"_LiquidDensity");
3928                                 _DensiteVapeur.writeCSV(allFields+"_SteamDensity");
3929                                 _EnthalpieLiquide.writeCSV(allFields+"_LiquidEnthalpy");
3930                                 _EnthalpieVapeur.writeCSV(allFields+"_SteamEnthalpy");
3931                                 _VitesseX.writeCSV(allFields+"_VelocityX");
3932                                 if(_Ndim>1)
3933                                 {
3934                                         _VitesseY.writeCSV(allFields+"_VelocityY");
3935                                         if(_Ndim>2)
3936                                                 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3937                                 }
3938                                 break;
3939                         }
3940                 }
3941         }
3942
3943         if (_restartWithNewFileName)
3944                 _restartWithNewFileName=false;
3945 }
3946
3947 void DriftModel::testConservation()
3948 {
3949         double SUM, DELTA, x;
3950         int I;
3951         for(int i=0; i<_nVar; i++)
3952         {
3953                 {
3954                         if(i == 0)
3955                                 cout << "Masse totale (kg): ";
3956                         else if(i == 1)
3957                                 cout << "Masse partielle 1 (kg): ";
3958                         else
3959                         {
3960                                 if(i == _nVar-1)
3961                                         cout << "Energie totale (J): ";
3962                                 else
3963                                         cout << "Quantite de mouvement direction "<<i-1<<" (kg.m.s^-1): ";
3964                         }
3965                 }
3966                 SUM = 0;
3967                 DELTA = 0;
3968                 I=i;
3969                 for(int j=0; j<_Nmailles; j++)
3970                 {
3971                         if(!_usePrimitiveVarsInNewton)
3972                                 VecGetValues(_old, 1, &I, &x);//on recupere la valeur du champ
3973                         else
3974                                 VecGetValues(_primitiveVars, 1, &I, &x);//on recupere la valeur du champ
3975                         SUM += x*_mesh.getCell(j).getMeasure();
3976                         VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
3977                         DELTA += x*_mesh.getCell(j).getMeasure();
3978                         I += _nVar;
3979                 }
3980                 if(fabs(SUM)>_precision)
3981                         cout << SUM << ", variation relative: " << fabs(DELTA /SUM)  << endl;
3982                 else
3983                         cout << " a une somme nulle,  variation absolue: " << fabs(DELTA) << endl;
3984         }
3985 }
3986