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