CoreFlows/Models/src/SinglePhase.cxx

   1 /*
   2  * SinglePhase.cxx
   3  *
   4  *  Created on: Sep 15, 2014
   5  *      Author: tn236279
   6  */
   7
   8 #include "SinglePhase.hxx"
   9
  10 using namespace std;
  11
  12 SinglePhase::SinglePhase(phaseType fluid, pressureEstimate pEstimate, int dim, bool useDellacherieEOS){
  13         _Ndim=dim;
  14         _nVar=_Ndim+2;
  15         _nbPhases  = 1;
  16         _dragCoeffs=vector<double>(1,0);
  17         _fluides.resize(1);
  18         _useDellacherieEOS=useDellacherieEOS;
  19         _saveAllFields=false;
  20
  21         if(pEstimate==around1bar300K){//EOS at 1 bar and 300K
  22                 _Tref=300;
  23                 _Pref=1e5;
  24                 if(fluid==Gas){
  25                         cout<<"Fluid is air around 1 bar and 300 K (27°C)"<<endl;
  26                         *_runLogFile<<"Fluid is air around 1 bar and 300 K (27°C)"<<endl;
  27                         _fluides[0] = new StiffenedGas(1.4,743,_Tref,2.22e5);  //ideal gas law for nitrogen at pressure 1 bar and temperature 27°C, e=2.22e5, c_v=743
  28                 }
  29                 else{
  30                         cout<<"Fluid is water around 1 bar and 300 K (27°C)"<<endl;
  31                         *_runLogFile<<"Fluid is water around 1 bar and 300 K (27°C)"<<endl;
  32                         _fluides[0] = new StiffenedGas(996,_Pref,_Tref,1.12e5,1501,4130);  //stiffened gas law for water at pressure 1 bar and temperature 27°C, e=1.12e5, c_v=4130
  33                 }
  34         }
  35         else{//EOS at 155 bars and 618K
  36                 _Tref=618;
  37                 _Pref=155e5;
  38                 if(fluid==Gas){
  39                         cout<<"Fluid is Gas around saturation point 155 bars and 618 K (345°C)"<<endl;
  40                         *_runLogFile<<"Fluid is Gas around saturation point 155 bars and 618 K (345°C)"<<endl;
  41                         _fluides[0] = new StiffenedGas(102,_Pref,_Tref,2.44e6, 433,3633);  //stiffened gas law for Gas at pressure 155 bar and temperature 345°C
  42                 }
  43                 else{//To do : change to normal regime: 155 bars and 573K
  44                         cout<<"Fluid is water around saturation point 155 bars and 618 K (345°C)"<<endl;
  45                         *_runLogFile<<"Fluid is water around saturation point 155 bars and 618 K (345°C)"<<endl;
  46                         if(_useDellacherieEOS)
  47                                 _fluides[0]= new StiffenedGasDellacherie(2.35,1e9,-1.167e6,1816); //stiffened gas law for water from S. Dellacherie
  48                         else
  49                                 _fluides[0]= new StiffenedGas(594.,_Pref,_Tref,1.6e6, 621.,3100.);  //stiffened gas law for water at pressure 155 bar, and temperature 345°C
  50                 }
  51         }
  52 }
  53 void SinglePhase::initialize(){
  54         cout<<"Initialising the Navier-Stokes model"<<endl;
  55         *_runLogFile<<"Initialising the Navier-Stokes model"<<endl;
  56
  57         _Uroe = new double[_nVar];
  58         _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
  59         for(int i=0; i<_Ndim; i++)
  60                 _gravite[i+1]=_GravityField3d[i];
  61
  62         _GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
  63         if(_saveVelocity || _saveAllFields)
  64                 _Vitesse=Field("Velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
  65
  66         if(_saveAllFields)
  67         {
  68                 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
  69                 _Pressure=Field("Pressure",CELLS,_mesh,1);
  70                 _Density=Field("Density",CELLS,_mesh,1);
  71                 _Temperature=Field("Temperature",CELLS,_mesh,1);
  72                 _MachNumber=Field("MachNumber",CELLS,_mesh,1);
  73                 _VitesseX=Field("Velocity x",CELLS,_mesh,1);
  74                 if(_Ndim>1)
  75                 {
  76                         _VitesseY=Field("Velocity y",CELLS,_mesh,1);
  77                         if(_Ndim>2)
  78                                 _VitesseZ=Field("Velocity z",CELLS,_mesh,1);
  79                 }
  80         }
  81
  82         if(_entropicCorrection)
  83                 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
  84
  85         ProblemFluid::initialize();
  86 }
  87
  88 bool SinglePhase::iterateTimeStep(bool &converged)
  89 {
  90         if(_timeScheme == Explicit || !_usePrimitiveVarsInNewton)
  91                 return ProblemFluid::iterateTimeStep(converged);
  92         else
  93         {
  94                 bool stop=false;
  95
  96                 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
  97                         _maxvp=0;
  98                         computeTimeStep(stop);
  99                 }
 100                 if(stop){//Le compute time step ne s'est pas bien passé
 101                         cout<<"ComputeTimeStep failed"<<endl;
 102                         converged=false;
 103                         return false;
 104                 }
 105
 106                 computeNewtonVariation();
 107
 108                 //converged=convergence des iterations
 109                 if(_timeScheme == Explicit)
 110                         converged=true;
 111                 else{//Implicit scheme
 112
 113                         KSPGetIterationNumber(_ksp, &_PetscIts);
 114                         if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
 115                                 _MaxIterLinearSolver = _PetscIts;
 116                         if(_PetscIts>=_maxPetscIts)//solving the linear system failed
 117                         {
 118                                 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 119                                 *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 120                                 converged=false;
 121                                 return false;
 122                         }
 123                         else{//solving the linear system succeeded
 124                                 //Calcul de la variation relative Uk+1-Uk
 125                                 _erreur_rel = 0.;
 126                                 double x, dx;
 127                                 int I;
 128                                 for(int j=1; j<=_Nmailles; j++)
 129                                 {
 130                                         for(int k=0; k<_nVar; k++)
 131                                         {
 132                                                 I = (j-1)*_nVar + k;
 133                                                 VecGetValues(_newtonVariation, 1, &I, &dx);
 134                                                 VecGetValues(_primitiveVars, 1, &I, &x);
 135                                                 if (fabs(x)*fabs(x)< _precision)
 136                                                 {
 137                                                         if(_erreur_rel < fabs(dx))
 138                                                                 _erreur_rel = fabs(dx);
 139                                                 }
 140                                                 else if(_erreur_rel < fabs(dx/x))
 141                                                         _erreur_rel = fabs(dx/x);
 142                                         }
 143                                 }
 144                         }
 145                         converged = _erreur_rel <= _precision_Newton;
 146                 }
 147
 148                 double relaxation=1;//Vk+1=Vk+relaxation*deltaV
 149
 150                 VecAXPY(_primitiveVars,  relaxation, _newtonVariation);
 151
 152                 //mise a jour du champ primitif
 153                 updateConservatives();
 154
 155                 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
 156                 {
 157                         cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 158                         *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 159                         converged=false;
 160                         return false;
 161                 }
 162                 if(_system)
 163                 {
 164                         cout<<"Vecteur Vkp1-Vk "<<endl;
 165                         VecView(_newtonVariation,  PETSC_VIEWER_STDOUT_SELF);
 166                         cout << "Nouvel etat courant Vk de l'iteration Newton: " << endl;
 167                         VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_SELF);
 168                 }
 169
 170                 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
 171                         if(_minm1<-_precision || _minm2<-_precision)
 172                         {
 173                                 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 174                                 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 175                         }
 176
 177                         if (_nbVpCplx>0){
 178                                 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 179                                 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 180                         }
 181                 }
 182                 _minm1=1e30;
 183                 _minm2=1e30;
 184                 _nbMaillesNeg=0;
 185                 _nbVpCplx =0;
 186                 _part_imag_max=0;
 187
 188                 return true;
 189         }
 190 }
 191 void SinglePhase::computeNewtonVariation()
 192 {
 193         if(!_usePrimitiveVarsInNewton)
 194                 ProblemFluid::computeNewtonVariation();
 195         else
 196         {
 197                 if(_verbose)
 198                 {
 199                         cout<<"Vecteur courant Vk "<<endl;
 200                         VecView(_primitiveVars,PETSC_VIEWER_STDOUT_SELF);
 201                         cout << endl;
 202                         cout << "Matrice du système linéaire avant contribution delta t" << endl;
 203                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 204                         cout << endl;
 205                         cout << "Second membre du système linéaire avant contribution delta t" << endl;
 206                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 207                         cout << endl;
 208                 }
 209                 if(_timeScheme == Explicit)
 210                 {
 211                         VecCopy(_b,_newtonVariation);
 212                         VecScale(_newtonVariation, _dt);
 213                         if(_verbose && _nbTimeStep%_freqSave ==0)
 214                         {
 215                                 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
 216                                 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
 217                                 cout << endl;
 218                         }
 219                 }
 220                 else
 221                 {
 222                         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 223
 224                         VecAXPY(_b, 1/_dt, _old);
 225                         VecAXPY(_b, -1/_dt, _conservativeVars);
 226
 227                         for(int imaille = 0; imaille<_Nmailles; imaille++){
 228                                 _idm[0] = _nVar*imaille;
 229                                 for(int k=1; k<_nVar; k++)
 230                                         _idm[k] = _idm[k-1] + 1;
 231                                 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
 232                                 primToConsJacobianMatrix(_Vi);
 233                                 for(int k=0; k<_nVar*_nVar; k++)
 234                                         _primToConsJacoMat[k]*=1/_dt;
 235                                 MatSetValuesBlocked(_A, 1, &imaille, 1, &imaille, _primToConsJacoMat, ADD_VALUES);
 236                         }
 237                         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 238
 239 #if PETSC_VERSION_GREATER_3_5
 240                         KSPSetOperators(_ksp, _A, _A);
 241 #else
 242                         KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 243 #endif
 244
 245                         if(_verbose)
 246                         {
 247                                 cout << "Matrice du système linéaire" << endl;
 248                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 249                                 cout << endl;
 250                                 cout << "Second membre du système linéaire" << endl;
 251                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 252                                 cout << endl;
 253                         }
 254
 255                         if(_conditionNumber)
 256                                 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 257                         if(!_isScaling)
 258                         {
 259                                 KSPSolve(_ksp, _b, _newtonVariation);
 260                         }
 261                         else
 262                         {
 263                                 computeScaling(_maxvp);
 264                                 int indice;
 265                                 VecAssemblyBegin(_vecScaling);
 266                                 VecAssemblyBegin(_invVecScaling);
 267                                 for(int imaille = 0; imaille<_Nmailles; imaille++)
 268                                 {
 269                                         indice = imaille;
 270                                         VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
 271                                         VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
 272                                 }
 273                                 VecAssemblyEnd(_vecScaling);
 274                                 VecAssemblyEnd(_invVecScaling);
 275
 276                                 if(_system)
 277                                 {
 278                                         cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
 279                                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 280                                         cout << endl;
 281                                         cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
 282                                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 283                                         cout << endl;
 284                                 }
 285                                 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
 286                                 if(_system)
 287                                 {
 288                                         cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
 289                                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 290                                         cout << endl;
 291                                 }
 292                                 VecCopy(_b,_bScaling);
 293                                 VecPointwiseMult(_b,_vecScaling,_bScaling);
 294                                 if(_system)
 295                                 {
 296                                         cout << "Produit du second membre par le preconditionneur bloc diagonal  a gauche" << endl;
 297                                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 298                                         cout << endl;
 299                                 }
 300
 301                                 KSPSolve(_ksp,_b, _bScaling);
 302                                 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
 303                         }
 304                         if(_system)
 305                         {
 306                                 cout << "solution du systeme lineaire local:" << endl;
 307                                 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
 308                                 cout << endl;
 309                         }
 310                 }
 311         }
 312 }
 313 void SinglePhase::convectionState( const long &i, const long &j, const bool &IsBord){
 314         //First conservative state then further down we will compute interface (Roe) state and then compute primitive state
 315         _idm[0] = _nVar*i;
 316         for(int k=1; k<_nVar; k++)
 317                 _idm[k] = _idm[k-1] + 1;
 318         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 319
 320         _idm[0] = _nVar*j;
 321         for(int k=1; k<_nVar; k++)
 322                 _idm[k] = _idm[k-1] + 1;
 323         if(IsBord)
 324                 VecGetValues(_Uext, _nVar, _idm, _Uj);
 325         else
 326                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 327         if(_verbose && _nbTimeStep%_freqSave ==0)
 328         {
 329                 cout<<"Convection Left state cell " << i<< ": "<<endl;
 330                 for(int k =0; k<_nVar; k++)
 331                         cout<< _Ui[k]<<endl;
 332                 cout<<"Convection Right state cell " << j<< ": "<<endl;
 333                 for(int k =0; k<_nVar; k++)
 334                         cout<< _Uj[k]<<endl;
 335         }
 336         if(_Ui[0]<0||_Uj[0]<0)
 337         {
 338                 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 339                 *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 340                 _runLogFile->close();
 341                 throw CdmathException("densite negative, arret de calcul");
 342         }
 343
 344         //Computation of Roe state
 345         PetscScalar ri, rj, xi, xj, pi, pj;
 346         PetscInt Ii;
 347         ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
 348         rj = sqrt(_Uj[0]);
 349         _Uroe[0] = ri*rj;       //moyenne geometrique des densites
 350         if(_verbose && _nbTimeStep%_freqSave ==0)
 351                 cout << "Densité moyenne Roe  gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
 352         for(int k=0;k<_Ndim;k++){
 353                 xi = _Ui[k+1];
 354                 xj = _Uj[k+1];
 355                 _Uroe[1+k] = (xi/ri + xj/rj)/(ri + rj);
 356                 //"moyenne" des vitesses
 357                 if(_verbose && _nbTimeStep%_freqSave ==0)
 358                         cout << "Vitesse de Roe composante "<< k<<"  gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+1] << endl;
 359         }
 360         // H = (rho E + p)/rho
 361         xi = _Ui[_nVar-1];//phi rho E
 362         xj = _Uj[_nVar-1];
 363         Ii = i*_nVar; // correct Kieu
 364         VecGetValues(_primitiveVars, 1, &Ii, &pi);// _primitiveVars pour p
 365         if(IsBord)
 366         {
 367                 double q_2 = 0;
 368                 for(int k=1;k<=_Ndim;k++)
 369                         q_2 += _Uj[k]*_Uj[k];
 370                 q_2 /= _Uj[0];  //phi rho u²
 371                 pj =  _fluides[0]->getPressure((_Uj[(_Ndim+2)-1]-q_2/2)/_porosityj,_Uj[0]/_porosityj);
 372         }
 373         else
 374         {
 375                 Ii = j*_nVar; // correct Kieu
 376                 VecGetValues(_primitiveVars, 1, &Ii, &pj);
 377         }
 378         xi = (xi + pi)/(ri*ri);
 379         xj = (xj + pj)/(rj*rj);
 380         _Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
 381         //on se donne l enthalpie ici
 382         if(_verbose && _nbTimeStep%_freqSave ==0)
 383                 cout << "Enthalpie totale de Roe H  gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
 384
 385         if(_verbose && _nbTimeStep%_freqSave ==0)
 386         {
 387                 cout<<"Convection interfacial state"<<endl;
 388                 for(int k=0;k<_nVar;k++)
 389                         cout<< _Uroe[k]<<" , "<<endl;
 390         }
 391
 392         //Extraction of primitive states
 393         _idm[0] = _nVar*i; // Kieu
 394         for(int k=1; k<_nVar; k++)
 395                 _idm[k] = _idm[k-1] + 1;
 396
 397         VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
 398         if (_verbose && _nbTimeStep%_freqSave ==0)
 399         {
 400                 cout << "Convection state: variables primitives maille " << i<<endl;
 401                 for(int q=0; q<_nVar; q++)
 402                         cout << _Vi[q] << endl;
 403                 cout << endl;
 404         }
 405
 406         if(!IsBord ){
 407                 for(int k=0; k<_nVar; k++)
 408                         _idn[k] = _nVar*j + k;
 409
 410                 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
 411         }
 412         else
 413         {
 414                 for(int k=0; k<_nVar; k++)
 415                         _idn[k] = k;
 416
 417                 VecGetValues(_Uextdiff, _nVar, _idn, _phi);
 418                 consToPrim(_phi,_Vj,1);
 419         }
 420
 421         if (_verbose && _nbTimeStep%_freqSave ==0)
 422         {
 423                 cout << "Convection state: variables primitives maille " <<j <<endl;
 424                 for(int q=0; q<_nVar; q++)
 425                         cout << _Vj[q] << endl;
 426                 cout << endl;
 427         }
 428 }
 429
 430 void SinglePhase::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
 431         //sortie: matrices et etat de diffusion (rho, q, T)
 432         _idm[0] = _nVar*i;
 433         for(int k=1; k<_nVar; k++)
 434                 _idm[k] = _idm[k-1] + 1;
 435
 436         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 437         for(int k=0; k<_nVar; k++)
 438                 _idn[k] = k;
 439
 440         if(IsBord)
 441                 VecGetValues(_Uextdiff, _nVar, _idn, _Uj);
 442         else
 443                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 444
 445         if(_verbose && _nbTimeStep%_freqSave ==0)
 446         {
 447                 cout << "SinglePhase::diffusionStateAndMatrices cellule gauche" << i << endl;
 448                 cout << "Ui = ";
 449                 for(int q=0; q<_nVar; q++)
 450                         cout << _Ui[q]  << "\t";
 451                 cout << endl;
 452                 cout << "SinglePhase::diffusionStateAndMatrices cellule droite" << j << endl;
 453                 cout << "Uj = ";
 454                 for(int q=0; q<_nVar; q++)
 455                         cout << _Uj[q]  << "\t";
 456                 cout << endl;
 457         }
 458
 459         for(int k=0; k<_nVar; k++)
 460                 _Udiff[k] = (_Ui[k]/_porosityi+_Uj[k]/_porosityj)/2;
 461
 462         if(_verbose && _nbTimeStep%_freqSave ==0)
 463         {
 464                 cout << "SinglePhase::diffusionStateAndMatrices conservative diffusion state" << endl;
 465                 cout << "_Udiff = ";
 466                 for(int q=0; q<_nVar; q++)
 467                         cout << _Udiff[q]  << "\t";
 468                 cout << endl;
 469                 cout << "porosite gauche= "<<_porosityi<< ", porosite droite= "<<_porosityj<<endl;
 470         }
 471         consToPrim(_Udiff,_phi,1);
 472         _Udiff[_nVar-1]=_phi[_nVar-1];
 473
 474         if(_timeScheme==Implicit)
 475         {
 476                 double q_2=0;
 477                 for (int i = 0; i<_Ndim;i++)
 478                         q_2+=_Udiff[i+1]*_Udiff[i+1];
 479                 double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
 480                 double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
 481                 double Cv= _fluides[0]->constante("Cv");
 482                 for(int i=0; i<_nVar*_nVar;i++)
 483                         _Diffusion[i] = 0;
 484                 for(int i=1;i<(_nVar-1);i++)
 485                 {
 486                         _Diffusion[i*_nVar] =  mu*_Udiff[i]/(_Udiff[0]*_Udiff[0]);
 487                         _Diffusion[i*_nVar+i] = -mu/_Udiff[0];
 488                 }
 489                 int i = (_nVar-1)*_nVar;
 490                 _Diffusion[i]=lambda*(_Udiff[_nVar-1]/_Udiff[0]-q_2/(2*Cv*_Udiff[0]*_Udiff[0]*_Udiff[0]));
 491                 for(int k=1;k<(_nVar-1);k++)
 492                 {
 493                         _Diffusion[i+k]= lambda*_Udiff[k]/(_Udiff[0]*_Udiff[0]*Cv);
 494                 }
 495                 _Diffusion[_nVar*_nVar-1]=-lambda/(_Udiff[0]*Cv);
 496         }
 497 }
 498 void SinglePhase::setBoundaryState(string nameOfGroup, const int &j,double *normale){
 499         _idm[0] = _nVar*j;
 500         for(int k=1; k<_nVar; k++)
 501                 _idm[k] = _idm[k-1] + 1;
 502
 503         double porosityj=_porosityField(j);
 504
 505         if(_verbose && _nbTimeStep%_freqSave ==0)
 506         {
 507                 cout << "setBoundaryState for group "<< nameOfGroup << ", inner cell j= "<<j<< " face unit normal vector "<<endl;
 508                 for(int k=0; k<_Ndim; k++){
 509                         cout<<normale[k]<<", ";
 510                 }
 511                 cout<<endl;
 512         }
 513
 514         if (_limitField[nameOfGroup].bcType==Wall){
 515                 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
 516                 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
 517                 for(int k=0; k<_Ndim; k++)
 518                         q_n+=_externalStates[(k+1)]*normale[k];
 519
 520                 //Pour la convection, inversion du sens de la vitesse normale
 521                 for(int k=0; k<_Ndim; k++)
 522                         _externalStates[(k+1)]-= 2*q_n*normale[k];
 523
 524                 _idm[0] = 0;
 525                 for(int k=1; k<_nVar; k++)
 526                         _idm[k] = _idm[k-1] + 1;
 527
 528                 VecAssemblyBegin(_Uext);
 529                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 530                 VecAssemblyEnd(_Uext);
 531
 532                 //Pour la diffusion, paroi à vitesse et temperature imposees
 533                 _idm[0] = _nVar*j;
 534                 for(int k=1; k<_nVar; k++)
 535                         _idm[k] = _idm[k-1] + 1;
 536                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//L'état fantome contient à présent les variables primitives internes
 537                 double pression=_externalStates[0];
 538                 double T=_limitField[nameOfGroup].T;
 539                 double rho=_fluides[0]->getDensity(pression,T);
 540
 541                 _externalStates[0]=porosityj*rho;
 542                 _externalStates[1]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
 543                 double v2=0;
 544                 v2 +=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
 545                 if(_Ndim>1)
 546                 {
 547                         v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 548                         _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
 549                         if(_Ndim==3)
 550                         {
 551                                 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
 552                                 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 553                         }
 554                 }
 555                 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
 556                 _idm[0] = 0;
 557                 for(int k=1; k<_nVar; k++)
 558                         _idm[k] = _idm[k-1] + 1;
 559                 VecAssemblyBegin(_Uextdiff);
 560                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 561                 VecAssemblyEnd(_Uextdiff);
 562         }
 563         else if (_limitField[nameOfGroup].bcType==Neumann){
 564                 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On prend l'état fantôme égal à l'état interne (conditions limites de Neumann)
 565
 566                 _idm[0] = 0;
 567                 for(int k=1; k<_nVar; k++)
 568                         _idm[k] = _idm[k-1] + 1;
 569
 570                 VecAssemblyBegin(_Uext);
 571                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 572                 VecAssemblyEnd(_Uext);
 573
 574                 VecAssemblyBegin(_Uextdiff);
 575                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 576                 VecAssemblyEnd(_Uextdiff);
 577         }
 578         else if (_limitField[nameOfGroup].bcType==Inlet){
 579                 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne (conditions limites de Neumann)
 580                 double q_int_n=0;//q_int_n=composante normale de la quantité de mouvement à la face frontière;
 581                 for(int k=0; k<_Ndim; k++)
 582                         q_int_n+=_externalStates[(k+1)]*normale[k];//On calcule la vitesse normale sortante
 583
 584                 double q_ext_n=_limitField[nameOfGroup].v_x[0]*normale[0];
 585                 if(_Ndim>1)
 586                         {
 587                                 q_ext_n+=_limitField[nameOfGroup].v_y[0]*normale[1];
 588                                 if(_Ndim>2)
 589                                                 q_ext_n+=_limitField[nameOfGroup].v_z[0]*normale[2];
 590                         }
 591
 592                 if(q_int_n+q_ext_n<=0){//Interfacial velocity goes inward
 593                         VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//On met à jour l'état fantome avec les variables primitives internes
 594                         double pression=_externalStates[0];
 595                         double T=_limitField[nameOfGroup].T;
 596                         double rho=_fluides[0]->getDensity(pression,T);
 597
 598                         _externalStates[0]=porosityj*rho;//Composante fantome de masse
 599                         _externalStates[1]=_externalStates[0]*(_limitField[nameOfGroup].v_x[0]);//Composante fantome de qdm x
 600                         double v2=0;
 601                         v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
 602                         if(_Ndim>1)
 603                         {
 604                                 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 605                                 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];//Composante fantome de qdm y
 606                                 if(_Ndim==3)
 607                                 {
 608                                         _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];//Composante fantome de qdm z
 609                                         v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 610                                 }
 611                         }
 612                         _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);//Composante fantome de l'nrj
 613                 }
 614                 else if(_nbTimeStep%_freqSave ==0)
 615                         cout<< "Warning : fluid possibly going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
 616
 617                 _idm[0] = 0;
 618                 for(int k=1; k<_nVar; k++)
 619                         _idm[k] = _idm[k-1] + 1;
 620                 VecAssemblyBegin(_Uext);
 621                 VecAssemblyBegin(_Uextdiff);
 622                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 623                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 624                 VecAssemblyEnd(_Uext);
 625                 VecAssemblyEnd(_Uextdiff);
 626         }
 627         else if (_limitField[nameOfGroup].bcType==InletRotationVelocity){
 628                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
 629                 double u_int_n=0;//u_int_n=composante normale de la vitesse intérieure à la face frontière;
 630                 for(int k=0; k<_Ndim; k++)
 631                         u_int_n+=_externalStates[(k+1)]*normale[k];
 632
 633                 double u_ext_n=0;
 634         //ghost velocity
 635         if(_Ndim>1)
 636         {
 637             Vector omega(3);
 638             omega[0]=_limitField[nameOfGroup].v_x[0];
 639             omega[1]=_limitField[nameOfGroup].v_y[0];
 640
 641             Vector Normale(3);
 642             Normale[0]=normale[0];
 643             Normale[1]=normale[1];
 644
 645             if(_Ndim==3)
 646             {
 647                 omega[2]=_limitField[nameOfGroup].v_z[0];
 648                 Normale[2]=normale[2];
 649             }
 650
 651             Vector tangent_vel=omega%Normale;
 652                         u_ext_n=-0.01*tangent_vel.norm();
 653                         //Changing external state velocity
 654             for(int k=0; k<_Ndim; k++)
 655                 _externalStates[(k+1)]=u_ext_n*normale[k] + tangent_vel[k];
 656         }
 657
 658                 if(u_ext_n + u_int_n <= 0){
 659                         double pression=_externalStates[0];
 660                         double T=_limitField[nameOfGroup].T;
 661                         double rho=_fluides[0]->getDensity(pression,T);
 662
 663                         double v2=0;
 664                         v2 +=_externalStates[1]*_externalStates[1];//v_x*v_x
 665                         _externalStates[0]=porosityj*rho;
 666                         _externalStates[1]*=_externalStates[0];
 667                         if(_Ndim>1)
 668                         {
 669                                 v2 +=_externalStates[2]*_externalStates[2];//+v_y*v_y
 670                                 _externalStates[2]*=_externalStates[0];
 671                                 if(_Ndim==3)
 672                                 {
 673                                         v2 +=_externalStates[3]*_externalStates[3];//+v_z*v_z
 674                                         _externalStates[3]*=_externalStates[0];
 675                                 }
 676                         }
 677                         _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
 678                 }
 679                 else if(_nbTimeStep%_freqSave ==0)
 680                 {
 681                         /*
 682                          * cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
 683                          */
 684                         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On définit l'état fantôme avec l'état interne
 685                         if(_nbTimeStep%_freqSave ==0)
 686                                 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
 687
 688                         //Changing external state momentum
 689             for(int k=0; k<_Ndim; k++)
 690                 _externalStates[(k+1)]-=2*_externalStates[0]*u_int_n*normale[k];
 691                 }
 692
 693                 _idm[0] = 0;
 694                 for(int k=1; k<_nVar; k++)
 695                         _idm[k] = _idm[k-1] + 1;
 696                 VecAssemblyBegin(_Uext);
 697                 VecAssemblyBegin(_Uextdiff);
 698                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 699                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 700                 VecAssemblyEnd(_Uext);
 701                 VecAssemblyEnd(_Uextdiff);
 702         }
 703         else if (_limitField[nameOfGroup].bcType==InletPressure){
 704                 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
 705
 706                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 707                 Cell Cj=_mesh.getCell(j);
 708                 double hydroPress=Cj.x()*_GravityField3d[0];
 709                 if(_Ndim>1){
 710                         hydroPress+=Cj.y()*_GravityField3d[1];
 711                         if(_Ndim>2)
 712                                 hydroPress+=Cj.z()*_GravityField3d[2];
 713                 }
 714                 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
 715
 716                 //Building the primitive external state
 717                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
 718                 double u_n=0;//u_n=vitesse normale à la face frontière;
 719                 for(int k=0; k<_Ndim; k++)
 720                         u_n+=_externalStates[(k+1)]*normale[k];
 721
 722                 if(u_n<=0)
 723                 {
 724                         _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress,_limitField[nameOfGroup].T);
 725
 726                 //Contribution from the tangential velocity
 727                 if(_Ndim>1)
 728                 {
 729                     Vector omega(3);
 730                     omega[0]=_limitField[nameOfGroup].v_x[0];
 731                     omega[1]=_limitField[nameOfGroup].v_y[0];
 732
 733                     Vector Normale(3);
 734                     Normale[0]=normale[0];
 735                     Normale[1]=normale[1];
 736
 737                     if(_Ndim==3)
 738                     {
 739                         omega[2]=_limitField[nameOfGroup].v_z[0];
 740                         Normale[2]=normale[2];
 741                     }
 742
 743                     Vector tangent_vel=omega%Normale;
 744
 745                                 //Changing external state velocity
 746                     for(int k=0; k<_Ndim; k++)
 747                         _externalStates[(k+1)]=u_n*normale[k] + abs(u_n)*tangent_vel[k];
 748                 }
 749                 }
 750                 else{
 751                         /*
 752                         if(_nbTimeStep%_freqSave ==0)
 753                                 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Neumann boundary condition for velocity and temperature (only pressure value is imposed as in outlet BC)."<<endl;
 754                         _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
 755                         */
 756                         if(_nbTimeStep%_freqSave ==0)
 757                                 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
 758                         _externalStates[0]=porosityj*_fluides[0]->getDensity(_externalStates[0]+hydroPress, _externalStates[_nVar-1]);
 759                         //Changing external state velocity
 760             for(int k=0; k<_Ndim; k++)
 761                 _externalStates[(k+1)]-=2*u_n*normale[k];
 762                 }
 763
 764                 double v2=0;
 765                 for(int k=0; k<_Ndim; k++)
 766                 {
 767                         v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
 768                         _externalStates[(k+1)]*=_externalStates[0] ;//qdm component
 769                 }
 770                 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);//nrj component
 771
 772
 773                 _idm[0] = 0;
 774                 for(int k=1; k<_nVar; k++)
 775                         _idm[k] = _idm[k-1] + 1;
 776                 VecAssemblyBegin(_Uext);
 777                 VecAssemblyBegin(_Uextdiff);
 778                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 779                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 780                 VecAssemblyEnd(_Uext);
 781                 VecAssemblyEnd(_Uextdiff);
 782         }
 783         else if (_limitField[nameOfGroup].bcType==Outlet){
 784                 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
 785                 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
 786                 for(int k=0; k<_Ndim; k++)
 787                         q_n+=_externalStates[(k+1)]*normale[k];
 788
 789                 if(q_n < -_precision &&  _nbTimeStep%_freqSave ==0)
 790         {
 791                         cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<" with flow rate "<< q_n<<endl;
 792             cout<< "Applying Neumann boundary condition for velocity and temperature"<<endl;
 793         }
 794                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 795                 Cell Cj=_mesh.getCell(j);
 796                 double hydroPress=Cj.x()*_GravityField3d[0];
 797                 if(_Ndim>1){
 798                         hydroPress+=Cj.y()*_GravityField3d[1];
 799                         if(_Ndim>2)
 800                                 hydroPress+=Cj.z()*_GravityField3d[2];
 801                 }
 802                 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
 803
 804                 if(_verbose && _nbTimeStep%_freqSave ==0)
 805                 {
 806                         cout<<"Cond lim outlet densite= "<<_externalStates[0]<<" gravite= "<<_GravityField3d[0]<<" Cj.x()= "<<Cj.x()<<endl;
 807                         cout<<"Cond lim outlet pression ref= "<<_limitField[nameOfGroup].p<<" pression hydro= "<<hydroPress<<" total= "<<_limitField[nameOfGroup].p+hydroPress<<endl;
 808                 }
 809                 //Building the external state
 810                 _idm[0] = _nVar*j;// Kieu
 811                 for(int k=1; k<_nVar; k++)
 812                         _idm[k] = _idm[k-1] + 1;
 813                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
 814
 815                 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
 816                 double v2=0;
 817                 for(int k=0; k<_Ndim; k++)
 818                 {
 819                         v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
 820                         _externalStates[(k+1)]*=_externalStates[0] ;
 821                 }
 822                 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);
 823                 _idm[0] = 0;
 824                 for(int k=1; k<_nVar; k++)
 825                         _idm[k] = _idm[k-1] + 1;
 826                 VecAssemblyBegin(_Uext);
 827                 VecAssemblyBegin(_Uextdiff);
 828                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 829                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 830                 VecAssemblyEnd(_Uext);
 831                 VecAssemblyEnd(_Uextdiff);
 832         }else {
 833                 cout<<"Boundary condition not set for boundary named "<<nameOfGroup<< " _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
 834                 cout<<"Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
 835                 *_runLogFile<<"Boundary condition not set for boundary named. Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
 836                 _runLogFile->close();
 837                 throw CdmathException("Unknown boundary condition");
 838         }
 839 }
 840
 841 void SinglePhase::convectionMatrices()
 842 {
 843         //entree: URoe = rho, u, H
 844         //sortie: matrices Roe+  et Roe-
 845
 846         if(_verbose && _nbTimeStep%_freqSave ==0)
 847                 cout<<"SinglePhase::convectionMatrices()"<<endl;
 848
 849         double u_n=0, u_2=0;//vitesse normale et carré du module
 850
 851         for(int i=0;i<_Ndim;i++)
 852         {
 853                 u_2 += _Uroe[1+i]*_Uroe[1+i];
 854                 u_n += _Uroe[1+i]*_vec_normal[i];
 855         }
 856
 857         vector<complex<double> > vp_dist(3,0);
 858
 859         if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
 860         {
 861                 staggeredVFFCMatricesConservativeVariables(u_n);//Computation of classical upwinding matrices
 862                 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//For use in implicit matrix
 863                         staggeredVFFCMatricesPrimitiveVariables(u_n);
 864         }
 865         else
 866         {
 867                 Vector vitesse(_Ndim);
 868                 for(int idim=0;idim<_Ndim;idim++)
 869                         vitesse[idim]=_Uroe[1+idim];
 870
 871                 double  c, H, K, k;
 872                 /***********Calcul des valeurs propres ********/
 873                 H = _Uroe[_nVar-1];
 874                 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
 875                 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
 876                 K = u_2*k/2; //g-1/2 *|u|²
 877
 878                 vp_dist[0]=u_n-c;vp_dist[1]=u_n;vp_dist[2]=u_n+c;
 879
 880                 _maxvploc=fabs(u_n)+c;
 881                 if(_maxvploc>_maxvp)
 882                         _maxvp=_maxvploc;
 883
 884                 if(_verbose && _nbTimeStep%_freqSave ==0)
 885                         cout<<"SinglePhase::convectionMatrices Eigenvalues "<<u_n-c<<" , "<<u_n<<" , "<<u_n+c<<endl;
 886
 887                 RoeMatrixConservativeVariables( u_n, H,vitesse,k,K);
 888
 889                 /******** Construction des matrices de decentrement ********/
 890                 if( _spaceScheme ==centered){
 891                         if(_entropicCorrection)
 892                         {
 893                                 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for centered scheme"<<endl;
 894                                 _runLogFile->close();
 895                                 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for centered scheme");
 896                         }
 897
 898                         for(int i=0; i<_nVar*_nVar;i++)
 899                                 _absAroe[i] = 0;
 900                 }
 901                 else if(_spaceScheme == upwind || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
 902                         if(_entropicCorrection)
 903                                 entropicShift(_vec_normal);
 904                         else
 905                                 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
 906
 907                         vector< complex< double > > y (3,0);
 908                         for( int i=0 ; i<3 ; i++)
 909                                 y[i] = Polynoms::abs_generalise(vp_dist[i])+_entropicShift[i];
 910                         Polynoms::abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
 911
 912                         if( _spaceScheme ==pressureCorrection)
 913                                 for( int i=0 ; i<_Ndim ; i++)
 914                                         for( int j=0 ; j<_Ndim ; j++)
 915                                                 _absAroe[(1+i)*_nVar+1+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 916                         else if( _spaceScheme ==lowMach){
 917                                 double M=sqrt(u_2)/c;
 918                                 for( int i=0 ; i<_Ndim ; i++)
 919                                         for( int j=0 ; j<_Ndim ; j++)
 920                                                 _absAroe[(1+i)*_nVar+1+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 921                         }
 922                 }
 923                 else if( _spaceScheme ==staggered ){
 924                         if(_entropicCorrection)//To do: study entropic correction for staggered
 925                         {
 926                                 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme"<<endl;
 927                                 _runLogFile->close();
 928                                 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme");
 929                         }
 930
 931                         staggeredRoeUpwindingMatrixConservativeVariables( u_n, H, vitesse, k, K);
 932                 }
 933                 else
 934                 {
 935                         *_runLogFile<<"SinglePhase::convectionMatrices: scheme not treated"<<endl;
 936                         _runLogFile->close();
 937                         throw CdmathException("SinglePhase::convectionMatrices: scheme not treated");
 938                 }
 939
 940                 for(int i=0; i<_nVar*_nVar;i++)
 941                 {
 942                         _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
 943                         _AroePlus[i]  = (_Aroe[i]+_absAroe[i])/2;
 944                 }
 945                 if(_timeScheme==Implicit)
 946                 {
 947                         if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
 948                         {
 949                                 _Vij[0]=_fluides[0]->getPressureFromEnthalpy(_Uroe[_nVar-1]-u_2/2, _Uroe[0]);//pressure
 950                                 _Vij[_nVar-1]=_fluides[0]->getTemperatureFromPressure( _Vij[0], _Uroe[0]);//Temperature
 951                                 for(int idim=0;idim<_Ndim; idim++)
 952                                         _Vij[1+idim]=_Uroe[1+idim];
 953                                 primToConsJacobianMatrix(_Vij);
 954                                 Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
 955                                 Polynoms::matrixProduct(_AroePlus,  _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
 956                         }
 957                         else
 958                                 for(int i=0; i<_nVar*_nVar;i++)
 959                                 {
 960                                         _AroeMinusImplicit[i] = _AroeMinus[i];
 961                                         _AroePlusImplicit[i]  = _AroePlus[i];
 962                                 }
 963                 }
 964                 if(_verbose && _nbTimeStep%_freqSave ==0)
 965                 {
 966                         displayMatrix(_Aroe, _nVar,"Matrice de Roe");
 967                         displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
 968                         displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
 969                         displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
 970                 }
 971         }
 972
 973         if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
 974         {
 975                 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
 976                 displayMatrix(_AroePlusImplicit,  _nVar,"Matrice _AroePlusImplicit");
 977         }
 978
 979         /*********Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source*****/
 980         if(_entropicCorrection)
 981         {
 982                 InvMatriceRoe( vp_dist);
 983                 Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
 984         }
 985         else if (_spaceScheme==upwind || (_spaceScheme ==pressureCorrection ) || (_spaceScheme ==lowMach ))//upwind sans entropic
 986                 SigneMatriceRoe( vp_dist);
 987         else if (_spaceScheme==centered)//centre  sans entropic
 988                 for(int i=0; i<_nVar*_nVar;i++)
 989                         _signAroe[i] = 0;
 990         else if( _spaceScheme ==staggered )//à tester
 991         {
 992                 double signu;
 993                 if(u_n>0)
 994                         signu=1;
 995                 else if (u_n<0)
 996                         signu=-1;
 997                 else
 998                         signu=0;
 999                 for(int i=0; i<_nVar*_nVar;i++)
1000                         _signAroe[i] = 0;
1001                 _signAroe[0] = signu;
1002                 for(int i=1; i<_nVar-1;i++)
1003                         _signAroe[i*_nVar+i] = -signu;
1004                 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
1005         }
1006         else
1007         {
1008                 *_runLogFile<<"SinglePhase::convectionMatrices: well balanced option not treated"<<endl;
1009                 _runLogFile->close();
1010                 throw CdmathException("SinglePhase::convectionMatrices: well balanced option not treated");
1011         }
1012 }
1013 void SinglePhase::computeScaling(double maxvp)
1014 {
1015         _blockDiag[0]=1;
1016         _invBlockDiag[0]=1;
1017         for(int q=1; q<_nVar-1; q++)
1018         {
1019                 _blockDiag[q]=1./maxvp;//
1020                 _invBlockDiag[q]= maxvp;//1.;//
1021         }
1022         _blockDiag[_nVar - 1]=(_fluides[0]->constante("gamma")-1)/(maxvp*maxvp);//1
1023         _invBlockDiag[_nVar - 1]=  1./_blockDiag[_nVar - 1] ;// 1.;//
1024 }
1025
1026 void SinglePhase::addDiffusionToSecondMember
1027 (               const int &i,
1028                 const int &j,
1029                 bool isBord)
1030
1031 {
1032         double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
1033         double mu     = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
1034
1035         if(isBord )
1036                 lambda=max(lambda,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
1037
1038         if(lambda==0 && mu ==0 && _heatTransfertCoeff==0)
1039                 return;
1040
1041         _idm[0] = 0;
1042         for(int k=1; k<_nVar; k++)
1043                 _idm[k] = _idm[k-1] + 1;
1044         if(isBord)
1045         {
1046                 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
1047                 _inv_dxj=_inv_dxi;
1048         }
1049
1050         //on n'a pas de contribution sur la masse
1051         _phi[0]=0;
1052         //contribution visqueuse sur la quantite de mouvement
1053         for(int k=1; k<_nVar-1; k++)
1054                 _phi[k] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu*(_porosityj*_Vj[k] - _porosityi*_Vi[k]);
1055         //contribution visqueuse sur l'energie
1056         _phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda*(_porosityj*_Vj[_nVar-1] - _porosityi*_Vi[_nVar-1]);
1057
1058         _idm[0] = i;
1059         VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1060
1061         if(_verbose && _nbTimeStep%_freqSave ==0)
1062         {
1063                 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
1064                 for(int q=0; q<_nVar; q++)
1065                         cout << _phi[q] << endl;
1066                 cout << endl;
1067         }
1068
1069         if(!isBord)
1070         {
1071                 //On change de signe pour l'autre contribution
1072                 for(int k=0; k<_nVar; k++)
1073                         _phi[k] *= -_inv_dxj/_inv_dxi;
1074
1075                 _idn[0] = j;
1076                 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1077                 if(_verbose && _nbTimeStep%_freqSave ==0)
1078                 {
1079                         cout << "Contribution diffusion au 2nd membre pour la maille  " << j << ": "<<endl;
1080                         for(int q=0; q<_nVar; q++)
1081                                 cout << _phi[q] << endl;
1082                         cout << endl;
1083                 }
1084         }
1085
1086         if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
1087         {
1088                 cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
1089                 for(int i=0; i<_nVar; i++)
1090                 {
1091                         for(int j=0; j<_nVar; j++)
1092                                 cout << _Diffusion[i*_nVar+j]<<", ";
1093                         cout << endl;
1094                 }
1095                 cout << endl;
1096         }
1097 }
1098
1099 void SinglePhase::sourceVector(PetscScalar * Si, PetscScalar * Ui, PetscScalar * Vi, int i)
1100 {
1101         double phirho=Ui[0], T=Vi[_nVar-1];
1102         double norm_u=0;
1103         for(int k=0; k<_Ndim; k++)
1104                 norm_u+=Vi[1+k]*Vi[1+k];
1105         norm_u=sqrt(norm_u);
1106         if(T>_Tsat)
1107                 Si[0]=_heatPowerField(i)/_latentHeat;
1108         else
1109                 Si[0]=0;
1110         for(int k=1; k<_nVar-1; k++)
1111                 Si[k]  =(_gravite[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*phirho;
1112
1113         Si[_nVar-1]=_heatPowerField(i);
1114
1115         for(int k=0; k<_Ndim; k++)
1116                 Si[_nVar-1] +=(_GravityField3d[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*Vi[1+k]*phirho;
1117
1118         if(_timeScheme==Implicit)
1119         {
1120                 for(int k=0; k<_nVar*_nVar;k++)
1121                         _GravityImplicitationMatrix[k] = 0;
1122                 if(!_usePrimitiveVarsInNewton)
1123                         for(int k=0; k<_nVar;k++)
1124                                 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
1125                 else
1126                 {
1127                         double pression=Vi[0];
1128                         getDensityDerivatives( pression, T, norm_u*norm_u);
1129                         for(int k=0; k<_nVar;k++)
1130                         {
1131                                 _GravityImplicitationMatrix[k*_nVar+0]      =-_gravite[k]*_drho_sur_dp;
1132                                 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
1133                         }
1134                 }
1135         }
1136         if(_verbose && _nbTimeStep%_freqSave ==0)
1137         {
1138                 cout<<"SinglePhase::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 SinglePhase::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
1157 {
1158         double norm_u=0, u_n=0, rho;
1159         for(int i=0;i<_Ndim;i++)
1160                 u_n += _Uroe[1+i]*_vec_normal[i];
1161
1162         pressureLoss[0]=0;
1163         if(u_n>0){
1164                 for(int i=0;i<_Ndim;i++)
1165                         norm_u += Vi[1+i]*Vi[1+i];
1166                 norm_u=sqrt(norm_u);
1167                 rho=Ui[0];
1168                 for(int i=0;i<_Ndim;i++)
1169                         pressureLoss[1+i]=-1/2*K*rho*norm_u*Vi[1+i];
1170         }
1171         else{
1172                 for(int i=0;i<_Ndim;i++)
1173                         norm_u += Vj[1+i]*Vj[1+i];
1174                 norm_u=sqrt(norm_u);
1175                 rho=Uj[0];
1176                 for(int i=0;i<_Ndim;i++)
1177                         pressureLoss[1+i]=-1/2*K*rho*norm_u*Vj[1+i];
1178         }
1179         pressureLoss[_nVar-1]=-1/2*K*rho*norm_u*norm_u*norm_u;
1180
1181         if(_verbose && _nbTimeStep%_freqSave ==0)
1182         {
1183                 cout<<"SinglePhase::pressureLossVector K= "<<K<<endl;
1184                 cout<<"Ui="<<endl;
1185                 for(int k=0;k<_nVar;k++)
1186                         cout<<Ui[k]<<", ";
1187                 cout<<endl;
1188                 cout<<"Vi="<<endl;
1189                 for(int k=0;k<_nVar;k++)
1190                         cout<<Vi[k]<<", ";
1191                 cout<<endl;
1192                 cout<<"Uj="<<endl;
1193                 for(int k=0;k<_nVar;k++)
1194                         cout<<Uj[k]<<", ";
1195                 cout<<endl;
1196                 cout<<"Vj="<<endl;
1197                 for(int k=0;k<_nVar;k++)
1198                         cout<<Vj[k]<<", ";
1199                 cout<<endl;
1200                 cout<<"pressureLoss="<<endl;
1201                 for(int k=0;k<_nVar;k++)
1202                         cout<<pressureLoss[k]<<", ";
1203                 cout<<endl;
1204         }
1205 }
1206
1207 void SinglePhase::porosityGradientSourceVector()
1208 {
1209         double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[0], pj=_Vj[0],pij;
1210         for(int i=0;i<_Ndim;i++) {
1211                 u_ni += _Vi[1+i]*_vec_normal[i];
1212                 u_nj += _Vj[1+i]*_vec_normal[i];
1213         }
1214         _porosityGradientSourceVector[0]=0;
1215         rhoj=_Uj[0]/_porosityj;
1216         rhoi=_Ui[0]/_porosityi;
1217         pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
1218         for(int i=0;i<_Ndim;i++)
1219                 _porosityGradientSourceVector[1+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1220         _porosityGradientSourceVector[_nVar-1]=0;
1221 }
1222
1223 void SinglePhase::jacobian(const int &j, string nameOfGroup,double * normale)
1224 {
1225         if(_verbose && _nbTimeStep%_freqSave ==0)
1226                 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
1227
1228         int k;
1229         for(k=0; k<_nVar*_nVar;k++)
1230                 _Jcb[k] = 0;//No implicitation at this stage
1231
1232         _idm[0] = _nVar*j;
1233         for(k=1; k<_nVar; k++)
1234                 _idm[k] = _idm[k-1] + 1;
1235         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
1236         double q_n=0;//quantité de mouvement normale à la paroi
1237         for(k=0; k<_Ndim; k++)
1238                 q_n+=_externalStates[(k+1)]*normale[k];
1239
1240         // loop of boundary types
1241         if (_limitField[nameOfGroup].bcType==Wall)
1242         {
1243                 for(k=0; k<_nVar;k++)
1244                         _Jcb[k*_nVar + k] = 1;
1245                 for(k=1; k<_nVar-1;k++)
1246                         for(int l=1; l<_nVar-1;l++)
1247                                 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
1248         }
1249         else if (_limitField[nameOfGroup].bcType==Inlet)
1250         {
1251                 return;
1252                 if(q_n<0){
1253                         double v[_Ndim], ve[_Ndim], v2, ve2;
1254
1255                         _idm[0] = _nVar*j;
1256                         for(k=1; k<_nVar; k++)
1257                                 _idm[k] = _idm[k-1] + 1;
1258                         VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1259                         VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1260
1261                         ve[0] = _limitField[nameOfGroup].v_x[0];
1262                         v[0]=_Vj[1];
1263                         ve2 = ve[0]*ve[0];
1264                         v2 = v[0]*v[0];
1265                         if (_Ndim >1){
1266                                 ve[1] = _limitField[nameOfGroup].v_y[0];
1267                                 v[1]=_Vj[2];
1268                                 ve2 += ve[1]*ve[1];
1269                                 v2 = v[1]*v[1];
1270                         }
1271                         if (_Ndim >2){
1272                                 ve[2] = _limitField[nameOfGroup].v_z[0];
1273                                 v[2]=_Vj[3];
1274                                 ve2 += ve[2]*ve[2];
1275                                 v2 = v[2]*v[2];
1276                         }
1277                         double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,_Uj[0]);
1278                         double total_energy=internal_energy+ve2/2;
1279
1280                         //Mass line
1281                         _Jcb[0]=v2/(2*internal_energy);
1282                         for(k=0; k<_Ndim;k++)
1283                                 _Jcb[1+k]=-v[k]/internal_energy;
1284                         _Jcb[_nVar-1]=1/internal_energy;
1285                         //Momentum lines
1286                         for(int l =1;l<1+_Ndim;l++){
1287                                 _Jcb[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1288                                 for(k=0; k<_Ndim;k++)
1289                                         _Jcb[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1290                                 _Jcb[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1291                         }
1292                         //Energy line
1293                         _Jcb[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1294                         for(k=0; k<_Ndim;k++)
1295                                 _Jcb[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1296                         _Jcb[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1297                 }
1298                 else
1299                         for(k=0;k<_nVar;k++)
1300                                 _Jcb[k*_nVar+k]=1;
1301                 //Old jacobian
1302                 /*
1303                  _Jcb[0] = 1;
1304                 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];//Kieu
1305                 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
1306                 if(_Ndim>1)
1307                 {
1308                         _Jcb[2*_nVar]= _limitField[nameOfGroup].v_y[0];
1309                         v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
1310                         if(_Ndim==3){
1311                                 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1312                                 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
1313                         }
1314                 }
1315                 _Jcb[(_nVar-1)*_nVar]=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2;
1316                  */
1317         }
1318         else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0){
1319                 return;
1320                 double v[_Ndim], v2=0;
1321                 _idm[0] = _nVar*j;
1322                 for(k=1; k<_nVar; k++)
1323                         _idm[k] = _idm[k-1] + 1;
1324                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1325
1326                 for(k=0; k<_Ndim;k++){
1327                         v[k]=_Vj[1+k];
1328                         v2+=v[k]*v[k];
1329                 }
1330
1331                 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _limitField[nameOfGroup].T);
1332                 double rho_int = _externalStates[0];
1333                 double density_ratio=rho_ext/rho_int;
1334                 //Momentum lines
1335                 for(int l =1;l<1+_Ndim;l++){
1336                         _Jcb[l*_nVar]=-density_ratio*v[l-1];
1337                         _Jcb[l*_nVar+l]=density_ratio;
1338                 }
1339                 //Energy lines
1340                 _Jcb[(_nVar-1)*_nVar]=-v2*density_ratio;
1341                 for(k=0; k<_Ndim;k++)
1342                         _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k];
1343         }
1344         // not wall, not inlet, not inletPressure
1345         else if(_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0))
1346         {
1347                 return;
1348                 double v[_Ndim], v2=0;
1349                 _idm[0] = _nVar*j;
1350                 for(k=1; k<_nVar; k++)
1351                         _idm[k] = _idm[k-1] + 1;
1352                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1353
1354                 for(k=0; k<_Ndim;k++){
1355                         v[k]=_Vj[1+k];
1356                         v2+=v[k]*v[k];
1357                 }
1358
1359                 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _externalStates[_nVar-1]);
1360                 double rho_int = _externalStates[0];
1361                 double density_ratio=rho_ext/rho_int;
1362                 double internal_energy=_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho_int);
1363                 double total_energy=internal_energy+v2/2;
1364
1365                 //Mass line
1366                 _Jcb[0]=density_ratio*(1-v2/(2*internal_energy));
1367                 for(k=0; k<_Ndim;k++)
1368                         _Jcb[1+k]=density_ratio*v[k]/internal_energy;
1369                 _Jcb[_nVar-1]=-density_ratio/internal_energy;
1370                 //Momentum lines
1371                 for(int l =1;l<1+_Ndim;l++){
1372                         _Jcb[l*_nVar]=density_ratio*v2*v[l-1]/(2*internal_energy);
1373                         for(k=0; k<_Ndim;k++)
1374                                 _Jcb[l*_nVar+1+k]=density_ratio*v[k]*v[l-1]/internal_energy;
1375                         _Jcb[l*_nVar+1+k]-=density_ratio;
1376                         _Jcb[l*_nVar+_nVar-1]=-density_ratio*v[l-1]/internal_energy;
1377                 }
1378                 //Energy line
1379                 _Jcb[(_nVar-1)*_nVar]=density_ratio*v2*total_energy/(2*internal_energy);
1380                 for(k=0; k<_Ndim;k++)
1381                         _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k]*total_energy/internal_energy;
1382                 _Jcb[(_nVar-1)*_nVar+_nVar-1]=density_ratio*(1-total_energy/internal_energy);
1383                 //Old jacobian
1384                 /*
1385                 int idim,jdim;
1386                 double cd = 1,cn=0,p0, gamma;
1387                 _idm[0] = j*_nVar;// Kieu
1388                 for(k=1; k<_nVar;k++)
1389                         _idm[k] = _idm[k-1] + 1;
1390                 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1391                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1392
1393                 // compute the common numerator and common denominator
1394                 p0=_fluides[0]->constante("p0");
1395                 gamma =_fluides[0]->constante("gamma");
1396                 cn =_limitField[nameOfGroup].p +p0;
1397                 cd = _phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0;
1398                 cd*=cd;
1399                 cd*=(gamma-1);
1400                 //compute the v2
1401                 for(k=1; k<_nVar-1;k++)
1402                         v2+=_externalStates[k]*_externalStates[k];
1403                 // drho_ext/dU
1404                 _JcbDiff[0] = cn*(_phi[_nVar-1] -v2 -p0)/cd;
1405                 for(k=1; k<_nVar-1;k++)
1406                         _JcbDiff[k]=cn*_phi[k]/cd;
1407                 _JcbDiff[_nVar-1]= -cn*_phi[0]/cd;
1408                 //dq_ext/dU
1409                 for(idim=0; idim<_Ndim;idim++)
1410                 {
1411                         //premiere colonne
1412                         _JcbDiff[(1+idim)*_nVar]=-(v2*cn*_phi[idim+1])/(2*cd);
1413                         //colonnes intermediaire
1414                         for(jdim=0; jdim<_Ndim;jdim++)
1415                         {
1416                                 _JcbDiff[(1+idim)*_nVar + jdim + 1] =_externalStates[idim+1]*_phi[jdim+1];
1417                                 _JcbDiff[(1+idim)*_nVar + jdim + 1]*=cn/cd;
1418                         }
1419                         //matrice identite*cn*(rhoe- p0)
1420                         _JcbDiff[(1+idim)*_nVar + idim + 1] +=( cn*(_phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0))/cd;
1421
1422                         //derniere colonne
1423                         _JcbDiff[(1+idim)*_nVar + _nVar-1]=-_phi[idim+1]*cn/cd;
1424                 }
1425                 //drhoE/dU
1426                 _JcbDiff[_nVar*(_nVar-1)] = -(v2*_phi[_nVar -1]*cn)/(2*cd);
1427                 for(int idim=0; idim<_Ndim;idim++)
1428                         _JcbDiff[_nVar*(_nVar-1)+idim+1]=_externalStates[idim+1]*_phi[_nVar -1]*cn/cd;
1429                 _JcbDiff[_nVar*_nVar -1] = -(v2/2+p0)*cn/cd;
1430                  */
1431         }
1432         else  if (_limitField[nameOfGroup].bcType!=Neumann)// not wall, not inlet, not outlet
1433         {
1434                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1435                 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1436                 _runLogFile->close();
1437                 throw CdmathException("SinglePhase::jacobian: This boundary condition is not treated");
1438         }
1439 }
1440
1441 //calcule la jacobienne pour une CL de diffusion
1442 void  SinglePhase::jacobianDiff(const int &j, string nameOfGroup)
1443 {
1444         if(_verbose && _nbTimeStep%_freqSave ==0)
1445                 cout<<"Jacobienne condition limite diffusion bord "<< nameOfGroup<<endl;
1446
1447         int k;
1448         for(k=0; k<_nVar*_nVar;k++)
1449                 _JcbDiff[k] = 0;
1450
1451         if (_limitField[nameOfGroup].bcType==Wall){
1452                 double v[_Ndim], ve[_Ndim], v2, ve2;
1453
1454                 _idm[0] = _nVar*j;
1455                 for(k=1; k<_nVar; k++)
1456                         _idm[k] = _idm[k-1] + 1;
1457                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1458                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1459
1460                 ve[0] = _limitField[nameOfGroup].v_x[0];
1461                 v[0]=_Vj[1];
1462                 ve2 = ve[0]*ve[0];
1463                 v2 = v[0]*v[0];
1464                 if (_Ndim >1){
1465                         ve[1] = _limitField[nameOfGroup].v_y[0];
1466                         v[1]=_Vj[2];
1467                         ve2 += ve[1]*ve[1];
1468                         v2 = v[1]*v[1];
1469                 }
1470                 if (_Ndim >2){
1471                         ve[2] = _limitField[nameOfGroup].v_z[0];
1472                         v[2]=_Vj[3];
1473                         ve2 += ve[2]*ve[2];
1474                         v2 = v[2]*v[2];
1475                 }
1476                 double rho=_Uj[0];
1477                 double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho);
1478                 double total_energy=internal_energy+ve2/2;
1479
1480                 //Mass line
1481                 _JcbDiff[0]=v2/(2*internal_energy);
1482                 for(k=0; k<_Ndim;k++)
1483                         _JcbDiff[1+k]=-v[k]/internal_energy;
1484                 _JcbDiff[_nVar-1]=1/internal_energy;
1485                 //Momentum lines
1486                 for(int l =1;l<1+_Ndim;l++){
1487                         _JcbDiff[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1488                         for(k=0; k<_Ndim;k++)
1489                                 _JcbDiff[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1490                         _JcbDiff[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1491                 }
1492                 //Energy line
1493                 _JcbDiff[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1494                 for(k=0; k<_Ndim;k++)
1495                         _JcbDiff[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1496                 _JcbDiff[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1497         }
1498         else if (_limitField[nameOfGroup].bcType==Outlet || _limitField[nameOfGroup].bcType==Neumann
1499                         ||_limitField[nameOfGroup].bcType==Inlet || _limitField[nameOfGroup].bcType==InletPressure)
1500         {
1501                 for(k=0;k<_nVar;k++)
1502                         _JcbDiff[k*_nVar+k]=1;
1503         }
1504         else{
1505                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1506                 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1507                 _runLogFile->close();
1508                 throw CdmathException("SinglePhase::jacobianDiff: This boundary condition is not recognised");
1509         }
1510 }
1511
1512 void SinglePhase::primToCons(const double *P, const int &i, double *W, const int &j){
1513         //cout<<"SinglePhase::primToCons i="<<i<<", j="<<j<<", P[i*(_Ndim+2)]="<<P[i*(_Ndim+2)]<<", P[i*(_Ndim+2)+(_Ndim+1)]="<<P[i*(_Ndim+2)+(_Ndim+1)]<<endl;
1514
1515         double rho=_fluides[0]->getDensity(P[i*(_Ndim+2)], P[i*(_Ndim+2)+(_Ndim+1)]);
1516         W[j*(_Ndim+2)] =  _porosityField(j)*rho;//phi*rho
1517         for(int k=0; k<_Ndim; k++)
1518                 W[j*(_Ndim+2)+(k+1)] = W[j*(_Ndim+2)]*P[i*(_Ndim+2)+(k+1)];//phi*rho*u
1519
1520         W[j*(_Ndim+2)+(_Ndim+1)] = W[j*(_Ndim+2)]*_fluides[0]->getInternalEnergy(P[i*(_Ndim+2)+ (_Ndim+1)],rho);//rho*e
1521         for(int k=0; k<_Ndim; k++)
1522                 W[j*(_Ndim+2)+(_Ndim+1)] += W[j*(_Ndim+2)]*P[i*(_Ndim+2)+(k+1)]*P[i*(_Ndim+2)+(k+1)]*0.5;//phi*rho*e+0.5*phi*rho*u^2
1523 }
1524
1525 void SinglePhase::primToConsJacobianMatrix(double *V)
1526 {
1527         double pression=V[0];
1528         double temperature=V[_nVar-1];
1529         double vitesse[_Ndim];
1530         for(int idim=0;idim<_Ndim;idim++)
1531                 vitesse[idim]=V[1+idim];
1532         double v2=0;
1533         for(int idim=0;idim<_Ndim;idim++)
1534                 v2+=vitesse[idim]*vitesse[idim];
1535
1536         double rho=_fluides[0]->getDensity(pression,temperature);
1537         double gamma=_fluides[0]->constante("gamma");
1538         double Pinf=_fluides[0]->constante("p0");
1539         double q=_fluides[0]->constante("q");
1540         double cv=_fluides[0]->constante("cv");
1541
1542         for(int k=0;k<_nVar*_nVar; k++)
1543                 _primToConsJacoMat[k]=0;
1544
1545         if(             !_useDellacherieEOS)
1546         {
1547                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1548                 double e=fluide0->getInternalEnergy(temperature);
1549                 double E=e+0.5*v2;
1550
1551                 _primToConsJacoMat[0]=1/((gamma-1)*(e-q));
1552                 _primToConsJacoMat[_nVar-1]=-rho*cv/(e-q);
1553
1554                 for(int idim=0;idim<_Ndim;idim++)
1555                 {
1556                         _primToConsJacoMat[_nVar+idim*_nVar]=vitesse[idim]/((gamma-1)*(e-q));
1557                         _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1558                         _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cv/(e-q);
1559                 }
1560                 _primToConsJacoMat[(_nVar-1)*_nVar]=E/((gamma-1)*(e-q));
1561                 for(int idim=0;idim<_Ndim;idim++)
1562                         _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1563                 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cv*(1-E/(e-q));
1564         }
1565         else if(        _useDellacherieEOS)
1566         {
1567                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1568                 double h=fluide0->getEnthalpy(temperature);
1569                 double H=h+0.5*v2;
1570                 double cp=_fluides[0]->constante("cp");
1571
1572                 _primToConsJacoMat[0]=gamma/((gamma-1)*(h-q));
1573                 _primToConsJacoMat[_nVar-1]=-rho*cp/(h-q);
1574
1575                 for(int idim=0;idim<_Ndim;idim++)
1576                 {
1577                         _primToConsJacoMat[_nVar+idim*_nVar]=gamma*vitesse[idim]/((gamma-1)*(h-q));
1578                         _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1579                         _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cp/(h-q);
1580                 }
1581                 _primToConsJacoMat[(_nVar-1)*_nVar]=gamma*H/((gamma-1)*(h-q))-1;
1582                 for(int idim=0;idim<_Ndim;idim++)
1583                         _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1584                 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cp*(1-H/(h-q));
1585         }
1586         else
1587         {
1588                 *_runLogFile<<"SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie"<<endl;
1589                 _runLogFile->close();
1590                 throw CdmathException("SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
1591         }
1592
1593         if(_verbose && _nbTimeStep%_freqSave ==0)
1594         {
1595                 cout<<" SinglePhase::primToConsJacobianMatrix" << endl;
1596                 displayVector(_Vi,_nVar," _Vi " );
1597                 cout<<" Jacobienne primToCons: " << endl;
1598                 displayMatrix(_primToConsJacoMat,_nVar," Jacobienne primToCons: ");
1599         }
1600 }
1601
1602 void SinglePhase::consToPrim(const double *Wcons, double* Wprim,double porosity)
1603 {
1604         double q_2 = 0;
1605         for(int k=1;k<=_Ndim;k++)
1606                 q_2 += Wcons[k]*Wcons[k];
1607         q_2 /= Wcons[0];        //phi rho u²
1608         double rhoe=(Wcons[(_Ndim+2)-1]-q_2/2)/porosity;
1609         double rho=Wcons[0]/porosity;
1610         Wprim[0] =  _fluides[0]->getPressure(rhoe,rho);//pressure p
1611         if (Wprim[0]<0){
1612                 cout << "pressure = "<< Wprim[0] << " < 0 " << endl;
1613                 *_runLogFile<< "pressure = "<< Wprim[0] << " < 0 " << endl;
1614                 _runLogFile->close();
1615                 throw CdmathException("SinglePhase::consToPrim: negative pressure");
1616         }
1617         for(int k=1;k<=_Ndim;k++)
1618                 Wprim[k] = Wcons[k]/Wcons[0];//velocity u
1619         Wprim[(_Ndim+2)-1] =  _fluides[0]->getTemperatureFromPressure(Wprim[0],Wcons[0]/porosity);
1620
1621         if(_verbose && _nbTimeStep%_freqSave ==0)
1622         {
1623                 cout<<"ConsToPrim Vecteur conservatif"<<endl;
1624                 for(int k=0;k<_nVar;k++)
1625                         cout<<Wcons[k]<<endl;
1626                 cout<<"ConsToPrim Vecteur primitif"<<endl;
1627                 for(int k=0;k<_nVar;k++)
1628                         cout<<Wprim[k]<<endl;
1629         }
1630 }
1631
1632 void SinglePhase::entropicShift(double* n)//TO do: make sure _Vi and _Vj are well set
1633 {
1634         /*Left and right values */
1635         double ul_n = 0, ul_2=0, ur_n=0,        ur_2 = 0; //valeurs de vitesse normale et |u|² a droite et a gauche
1636         for(int i=0;i<_Ndim;i++)
1637         {
1638                 ul_n += _Vi[1+i]*n[i];
1639                 ul_2 += _Vi[1+i]*_Vi[1+i];
1640                 ur_n += _Vj[1+i]*n[i];
1641                 ur_2 += _Vj[1+i]*_Vj[1+i];
1642         }
1643
1644         double cl = _fluides[0]->vitesseSonEnthalpie(_Vi[_Ndim+1]-ul_2/2);//vitesse du son a l'interface
1645         double cr = _fluides[0]->vitesseSonEnthalpie(_Vj[_Ndim+1]-ur_2/2);//vitesse du son a l'interface
1646
1647         _entropicShift[0]=fabs(ul_n-cl - (ur_n-cr));
1648         _entropicShift[1]=fabs(ul_n     - ur_n);
1649         _entropicShift[2]=fabs(ul_n+cl - (ur_n+cr));
1650
1651         if(_verbose && _nbTimeStep%_freqSave ==0)
1652         {
1653                 cout<<"Entropic shift left eigenvalues: "<<endl;
1654                 cout<<"("<< ul_n-cl<< ", " <<ul_n<<", "<<ul_n+cl << ")";
1655                 cout<<endl;
1656                 cout<<"Entropic shift right eigenvalues: "<<endl;
1657                 cout<<"("<< ur_n-cr<< ", " <<ur_n<<", "<<ur_n+cr << ")";
1658                 cout<< endl;
1659                 cout<<"eigenvalue jumps "<<endl;
1660                 cout<< _entropicShift[0] << ", " << _entropicShift[1] << ", "<< _entropicShift[2] <<endl;
1661         }
1662 }
1663
1664 Vector SinglePhase::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
1665         if(_verbose && _nbTimeStep%_freqSave ==0)
1666         {
1667                 cout<<"SinglePhase::convectionFlux start"<<endl;
1668                 cout<<"Ucons"<<endl;
1669                 cout<<U<<endl;
1670                 cout<<"Vprim"<<endl;
1671                 cout<<V<<endl;
1672         }
1673
1674         double phirho=U(0);//phi rho
1675         Vector phiq(_Ndim);//phi rho u
1676         for(int i=0;i<_Ndim;i++)
1677                 phiq(i)=U(1+i);
1678
1679         double pression=V(0);
1680         Vector vitesse(_Ndim);
1681         for(int i=0;i<_Ndim;i++)
1682                 vitesse(i)=V(1+i);
1683         double Temperature= V(1+_Ndim);
1684
1685         double vitessen=vitesse*normale;
1686         double rho=phirho/porosity;
1687         double e_int=_fluides[0]->getInternalEnergy(Temperature,rho);
1688
1689         Vector F(_nVar);
1690         F(0)=phirho*vitessen;
1691         for(int i=0;i<_Ndim;i++)
1692                 F(1+i)=phirho*vitessen*vitesse(i)+pression*normale(i)*porosity;
1693         F(1+_Ndim)=phirho*(e_int+0.5*vitesse*vitesse+pression/rho)*vitessen;
1694
1695         if(_verbose && _nbTimeStep%_freqSave ==0)
1696         {
1697                 cout<<"SinglePhase::convectionFlux end"<<endl;
1698                 cout<<"Flux F(U,V)"<<endl;
1699                 cout<<F<<endl;
1700         }
1701
1702         return F;
1703 }
1704
1705 void SinglePhase::RoeMatrixConservativeVariables(double u_n, double H,Vector velocity, double k, double K)
1706 {
1707         /******** Construction de la matrice de Roe *********/
1708         //premiere ligne (masse)
1709         _Aroe[0]=0;
1710         for(int idim=0; idim<_Ndim;idim++)
1711                 _Aroe[1+idim]=_vec_normal[idim];
1712         _Aroe[_nVar-1]=0;
1713
1714         //lignes intermadiaires(qdm)
1715         for(int idim=0; idim<_Ndim;idim++)
1716         {
1717                 //premiere colonne
1718                 _Aroe[(1+idim)*_nVar]=K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1719                 //colonnes intermediaires
1720                 for(int jdim=0; jdim<_Ndim;jdim++)
1721                         _Aroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]-k*_vec_normal[idim]*_Uroe[1+jdim];
1722                 //matrice identite
1723                 _Aroe[(1+idim)*_nVar + idim + 1] += u_n;
1724                 //derniere colonne
1725                 _Aroe[(1+idim)*_nVar + _nVar-1]=k*_vec_normal[idim];
1726         }
1727
1728         //derniere ligne (energie)
1729         _Aroe[_nVar*(_nVar-1)] = (K - H)*u_n;
1730         for(int idim=0; idim<_Ndim;idim++)
1731                 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - k*u_n*_Uroe[idim+1];
1732         _Aroe[_nVar*_nVar -1] = (1 + k)*u_n;
1733 }
1734 void SinglePhase::convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector vitesse)
1735 {
1736         //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1737         //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
1738         //EOS is more involved with primitive variables
1739         // call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
1740         _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1741         for(int i=0;i<_Ndim;i++)
1742                 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1743         _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1744         for(int i=0;i<_Ndim;i++)
1745         {
1746                 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]+_vec_normal[i];
1747                 for(int j=0;j<_Ndim;j++)
1748                         _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1749                 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1750                 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1751         }
1752         _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1753         for(int i=0;i<_Ndim;i++)
1754                 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1755         _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1756 }
1757 void SinglePhase::staggeredRoeUpwindingMatrixConservativeVariables( double u_n, double H,Vector velocity, double k, double K)
1758 {
1759         //Calcul de décentrement de type décalé pour formulation de Roe
1760         if(fabs(u_n)>_precision)
1761         {
1762                 //premiere ligne (masse)
1763                 _absAroe[0]=0;
1764                 for(int idim=0; idim<_Ndim;idim++)
1765                         _absAroe[1+idim]=_vec_normal[idim];
1766                 _absAroe[_nVar-1]=0;
1767
1768                 //lignes intermadiaires(qdm)
1769                 for(int idim=0; idim<_Ndim;idim++)
1770                 {
1771                         //premiere colonne
1772                         _absAroe[(1+idim)*_nVar]=-K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1773                         //colonnes intermediaires
1774                         for(int jdim=0; jdim<_Ndim;jdim++)
1775                                 _absAroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]+k*_vec_normal[idim]*_Uroe[1+jdim];
1776                         //matrice identite
1777                         _absAroe[(1+idim)*_nVar + idim + 1] += u_n;
1778                         //derniere colonne
1779                         _absAroe[(1+idim)*_nVar + _nVar-1]=-k*_vec_normal[idim];
1780                 }
1781
1782                 //derniere ligne (energie)
1783                 _absAroe[_nVar*(_nVar-1)] = (-K - H)*u_n;
1784                 for(int idim=0; idim<_Ndim;idim++)
1785                         _absAroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] + k*u_n*_Uroe[idim+1];
1786                 _absAroe[_nVar*_nVar -1] = (1 - k)*u_n;
1787
1788                 double signu;
1789                 if(u_n>0)
1790                         signu=1;
1791                 else if (u_n<0)
1792                         signu=-1;
1793
1794                 for(int i=0; i<_nVar*_nVar;i++)
1795                         _absAroe[i] *= signu;
1796         }
1797         else//umn=0 ->centered scheme
1798         {
1799                 for(int i=0; i<_nVar*_nVar;i++)
1800                         _absAroe[i] =0;
1801         }
1802 }
1803 void SinglePhase::staggeredRoeUpwindingMatrixPrimitiveVariables(double rho, double u_n,double H, Vector vitesse)
1804 {
1805         //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1806         //Calcul de décentrement de type décalé pour formulation Roe
1807         _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1808         for(int i=0;i<_Ndim;i++)
1809                 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1810         _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1811         for(int i=0;i<_Ndim;i++)
1812         {
1813                 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]-_vec_normal[i];
1814                 for(int j=0;j<_Ndim;j++)
1815                         _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1816                 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1817                 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1818         }
1819         _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1820         for(int i=0;i<_Ndim;i++)
1821                 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1822         _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1823 }
1824
1825 Vector SinglePhase::staggeredVFFCFlux()
1826 {
1827         if(_verbose && _nbTimeStep%_freqSave ==0)
1828                 cout<<"SinglePhase::staggeredVFFCFlux() start"<<endl;
1829
1830         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1831         {
1832                 *_runLogFile<< "SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding, pressure = "<<  endl;
1833                 _runLogFile->close();
1834                 throw CdmathException("SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
1835         }
1836         else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1837         {
1838                 Vector Fij(_nVar);
1839
1840                 double uijn=0, phiqn=0, uin=0, ujn=0;
1841                 for(int idim=0;idim<_Ndim;idim++)
1842                 {
1843                         uijn+=_vec_normal[idim]*_Uroe[1+idim];//URoe = rho, u, H
1844                         uin +=_vec_normal[idim]*_Ui[1+idim];
1845                         ujn +=_vec_normal[idim]*_Uj[1+idim];
1846                 }
1847
1848                 if( (uin>0 && ujn >0) || (uin>=0 && ujn <=0 && uijn>0) ) // formerly (uijn>_precision)
1849                 {
1850                         for(int idim=0;idim<_Ndim;idim++)
1851                                 phiqn+=_vec_normal[idim]*_Ui[1+idim];//phi rho u n
1852                         Fij(0)=phiqn;
1853                         for(int idim=0;idim<_Ndim;idim++)
1854                                 Fij(1+idim)=phiqn*_Vi[1+idim]+_Vj[0]*_vec_normal[idim]*_porosityj;
1855                         Fij(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[0]*sqrt(_porosityj/_porosityi));
1856                 }
1857                 else if( (uin<0 && ujn <0) || (uin>=0 && ujn <=0 && uijn<0) ) // formerly (uijn<-_precision)
1858                 {
1859                         for(int idim=0;idim<_Ndim;idim++)
1860                                 phiqn+=_vec_normal[idim]*_Uj[1+idim];//phi rho u n
1861                         Fij(0)=phiqn;
1862                         for(int idim=0;idim<_Ndim;idim++)
1863                                 Fij(1+idim)=phiqn*_Vj[1+idim]+_Vi[0]*_vec_normal[idim]*_porosityi;
1864                         Fij(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[0]*sqrt(_porosityi/_porosityj));
1865                 }
1866                 else//case (uin<=0 && ujn >=0) or (uin>=0 && ujn <=0 && uijn==0), apply centered scheme
1867                 {
1868                         Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar);
1869                         Vector normale(_Ndim);
1870                         for(int i1=0;i1<_Ndim;i1++)
1871                                 normale(i1)=_vec_normal[i1];
1872                         for(int i1=0;i1<_nVar;i1++)
1873                         {
1874                                 Ui(i1)=_Ui[i1];
1875                                 Uj(i1)=_Uj[i1];
1876                                 Vi(i1)=_Vi[i1];
1877                                 Vj(i1)=_Vj[i1];
1878                         }
1879                         Fi=convectionFlux(Ui,Vi,normale,_porosityi);
1880                         Fj=convectionFlux(Uj,Vj,normale,_porosityj);
1881                         Fij=(Fi+Fj)/2;//+_maxvploc*(Ui-Uj)/2;
1882                 }
1883                 if(_verbose && _nbTimeStep%_freqSave ==0)
1884                 {
1885                         cout<<"SinglePhase::staggeredVFFCFlux() endf uijn="<<uijn<<endl;
1886                         cout<<Fij<<endl;
1887                 }
1888                 return Fij;
1889         }
1890 }
1891
1892 void SinglePhase::staggeredVFFCMatricesConservativeVariables(double un)//vitesse normale de Roe en entrée
1893 {
1894         if(_verbose && _nbTimeStep%_freqSave ==0)
1895                 cout<<"SinglePhase::staggeredVFFCMatrices()"<<endl;
1896
1897         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1898         {
1899                 *_runLogFile<< "SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding"<< endl;
1900                 _runLogFile->close();
1901                 throw CdmathException("SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding");
1902         }
1903         else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1904         {
1905                 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
1906                 double H;//enthalpie totale (expression particulière)
1907                 consToPrim(_Ui,_Vi,_porosityi);
1908                 consToPrim(_Uj,_Vj,_porosityj);
1909                 for(int i=0;i<_Ndim;i++)
1910                 {
1911                         ui_2 += _Vi[1+i]*_Vi[1+i];
1912                         ui_n += _Vi[1+i]*_vec_normal[i];
1913                         uj_2 += _Vj[1+i]*_Vj[1+i];
1914                         uj_n += _Vj[1+i]*_vec_normal[i];
1915                 }
1916
1917                 double rhoi,pj, Ei, rhoj;
1918                 double cj, Kj, kj;//dérivées de la pression
1919                 rhoi=_Ui[0]/_porosityi;
1920                 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
1921                 pj=_Vj[0];
1922                 rhoj=_Uj[0]/_porosityj;
1923                 H =Ei+pj/rhoi;
1924                 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
1925                 kj = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1926                 Kj = uj_2*kj/2; //g-1/2 *|u|²
1927
1928                 double pi, Ej;
1929                 double ci, Ki, ki;//dérivées de la pression
1930                 Ej= _Uj[_Ndim+1]/rhoj;
1931                 pi=_Vi[0];
1932                 H =Ej+pi/rhoj;
1933                 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
1934                 ki = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1935                 Ki = ui_2*ki/2; //g-1/2 *|u|²
1936
1937                 if(un>_precision)
1938                 {
1939                         /***********Calcul des valeurs propres ********/
1940                         vector<complex<double> > vp_dist(3,0);
1941                         vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
1942
1943                         _maxvploc=fabs(ui_n)+cj;
1944                         if(_maxvploc>_maxvp)
1945                                 _maxvp=_maxvploc;
1946
1947                         if(_verbose && _nbTimeStep%_freqSave ==0)
1948                                 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
1949
1950                         /******** Construction de la matrice A^+ *********/
1951                         //premiere ligne (masse)
1952                         _AroePlus[0]=0;
1953                         for(int idim=0; idim<_Ndim;idim++)
1954                                 _AroePlus[1+idim]=_vec_normal[idim];
1955                         _AroePlus[_nVar-1]=0;
1956
1957                         //lignes intermadiaires(qdm)
1958                         for(int idim=0; idim<_Ndim;idim++)
1959                         {
1960                                 //premiere colonne
1961                                 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
1962                                 //colonnes intermediaires
1963                                 for(int jdim=0; jdim<_Ndim;jdim++)
1964                                         _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim];
1965                                 //matrice identite
1966                                 _AroePlus[(1+idim)*_nVar + idim + 1] += ui_n;
1967                                 //derniere colonne
1968                                 _AroePlus[(1+idim)*_nVar + _nVar-1]=0;
1969                         }
1970
1971                         //derniere ligne (energie)
1972                         _AroePlus[_nVar*(_nVar-1)] = - H*ui_n;
1973                         for(int idim=0; idim<_Ndim;idim++)
1974                                 _AroePlus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
1975                         _AroePlus[_nVar*_nVar -1] = ui_n;
1976
1977                         /******** Construction de la matrice A^- *********/
1978                         //premiere ligne (masse)
1979                         _AroeMinus[0]=0;
1980                         for(int idim=0; idim<_Ndim;idim++)
1981                                 _AroeMinus[1+idim]=0;
1982                         _AroeMinus[_nVar-1]=0;
1983
1984                         //lignes intermadiaires(qdm)
1985                         for(int idim=0; idim<_Ndim;idim++)
1986                         {
1987                                 //premiere colonne
1988                                 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] ;
1989                                 //colonnes intermediaires
1990                                 for(int jdim=0; jdim<_Ndim;jdim++)
1991                                         _AroeMinus[(1+idim)*_nVar + jdim + 1] = -kj*_vec_normal[idim]*_Vj[1+jdim];
1992                                 //matrice identite
1993                                 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0;
1994                                 //derniere colonne
1995                                 _AroeMinus[(1+idim)*_nVar + _nVar-1]=kj*_vec_normal[idim];
1996                         }
1997
1998                         //derniere ligne (energie)
1999                         _AroeMinus[_nVar*(_nVar-1)] = Kj *ui_n;
2000                         for(int idim=0; idim<_Ndim;idim++)
2001                                 _AroeMinus[_nVar*(_nVar-1)+idim+1]= - kj*uj_n*_Vi[idim+1];
2002                         _AroeMinus[_nVar*_nVar -1] = kj*ui_n;
2003                 }
2004                 else if(un<-_precision)
2005                 {
2006                         /***********Calcul des valeurs propres ********/
2007                         vector<complex<double> > vp_dist(3,0);
2008                         vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2009
2010                         _maxvploc=fabs(uj_n)+ci;
2011                         if(_maxvploc>_maxvp)
2012                                 _maxvp=_maxvploc;
2013
2014                         if(_verbose && _nbTimeStep%_freqSave ==0)
2015                                 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2016
2017                         /******** Construction de la matrice A^+ *********/
2018                         //premiere ligne (masse)
2019                         _AroePlus[0]=0;
2020                         for(int idim=0; idim<_Ndim;idim++)
2021                                 _AroePlus[1+idim]=0;
2022                         _AroePlus[_nVar-1]=0;
2023
2024                         //lignes intermadiaires(qdm)
2025                         for(int idim=0; idim<_Ndim;idim++)
2026                         {
2027                                 //premiere colonne
2028                                 _AroePlus[(1+idim)*_nVar]=Ki*_vec_normal[idim] ;
2029                                 //colonnes intermediaires
2030                                 for(int jdim=0; jdim<_Ndim;jdim++)
2031                                         _AroePlus[(1+idim)*_nVar + jdim + 1] = -ki*_vec_normal[idim]*_Vi[1+jdim];
2032                                 //matrice identite
2033                                 _AroePlus[(1+idim)*_nVar + idim + 1] += 0;
2034                                 //derniere colonne
2035                                 _AroePlus[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2036                         }
2037
2038                         //derniere ligne (energie)
2039                         _AroePlus[_nVar*(_nVar-1)] = Ki *uj_n;
2040                         for(int idim=0; idim<_Ndim;idim++)
2041                                 _AroePlus[_nVar*(_nVar-1)+idim+1]= - ki*ui_n*_Vj[idim+1];
2042                         _AroePlus[_nVar*_nVar -1] =  ki*uj_n;
2043
2044                         /******** Construction de la matrice A^- *********/
2045                         //premiere ligne (masse)
2046                         _AroeMinus[0]=0;
2047                         for(int idim=0; idim<_Ndim;idim++)
2048                                 _AroeMinus[1+idim]=_vec_normal[idim];
2049                         _AroeMinus[_nVar-1]=0;
2050
2051                         //lignes intermadiaires(qdm)
2052                         for(int idim=0; idim<_Ndim;idim++)
2053                         {
2054                                 //premiere colonne
2055                                 _AroeMinus[(1+idim)*_nVar]= - uj_n*_Vj[1+idim];
2056                                 //colonnes intermediaires
2057                                 for(int jdim=0; jdim<_Ndim;jdim++)
2058                                         _AroeMinus[(1+idim)*_nVar + jdim + 1] = _Vj[1+idim]*_vec_normal[jdim];
2059                                 //matrice identite
2060                                 _AroeMinus[(1+idim)*_nVar + idim + 1] += uj_n;
2061                                 //derniere colonne
2062                                 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0;
2063                         }
2064
2065                         //derniere ligne (energie)
2066                         _AroeMinus[_nVar*(_nVar-1)] = - H*uj_n;
2067                         for(int idim=0; idim<_Ndim;idim++)
2068                                 _AroeMinus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
2069                         _AroeMinus[_nVar*_nVar -1] = uj_n;
2070                 }
2071                 else//case nil velocity on the interface, apply centered scheme
2072                 {
2073                         double u_n=0, u_2=0;//vitesse normale et carré du module
2074                         for(int i=0;i<_Ndim;i++)
2075                         {
2076                                 u_2 += _Uroe[1+i]*_Uroe[1+i];
2077                                 u_n += _Uroe[1+i]*_vec_normal[i];
2078                         }
2079                         Vector vitesse(_Ndim);
2080                         for(int idim=0;idim<_Ndim;idim++)
2081                                 vitesse[idim]=_Uroe[1+idim];
2082
2083                         double  c, H, K, k;
2084                         /***********Calcul des valeurs propres ********/
2085                         H = _Uroe[_nVar-1];
2086                         c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
2087                         k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
2088                         K = u_2*k/2; //g-1/2 *|u|²
2089
2090                         _maxvploc=fabs(u_n)+c;
2091                         if(_maxvploc>_maxvp)
2092                                 _maxvp=_maxvploc;
2093
2094                         /******** Construction de la matrice A^+ *********/
2095                         //premiere ligne (masse)
2096                         _AroePlus[0]=0;
2097                         for(int idim=0; idim<_Ndim;idim++)
2098                                 _AroePlus[1+idim]=0;
2099                         _AroePlus[_nVar-1]=0;
2100
2101                         //lignes intermadiaires(qdm)
2102                         for(int idim=0; idim<_Ndim;idim++)
2103                         {
2104                                 //premiere colonne
2105                                 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
2106                                 //colonnes intermediaires
2107                                 for(int jdim=0; jdim<_Ndim;jdim++)
2108                                         _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim]-0.5*ki*_vec_normal[idim]*_Vi[1+jdim];
2109                                 //matrice identite
2110                                 _AroePlus[(1+idim)*_nVar + idim + 1] += 0.5*ui_n;
2111                                 //derniere colonne
2112                                 _AroePlus[(1+idim)*_nVar + _nVar-1]=0.5*ki*_vec_normal[idim];
2113                         }
2114
2115                         //derniere ligne (energie)
2116                         _AroePlus[_nVar*(_nVar-1)] = 0;
2117                         for(int idim=0; idim<_Ndim;idim++)
2118                                 _AroePlus[_nVar*(_nVar-1)+idim+1]=0 ;
2119                         _AroePlus[_nVar*_nVar -1] = 0;
2120
2121                         /******** Construction de la matrice A^- *********/
2122                         //premiere ligne (masse)
2123                         _AroeMinus[0]=0;
2124                         for(int idim=0; idim<_Ndim;idim++)
2125                                 _AroeMinus[1+idim]=0;
2126                         _AroeMinus[_nVar-1]=0;
2127
2128                         //lignes intermadiaires(qdm)
2129                         for(int idim=0; idim<_Ndim;idim++)
2130                         {
2131                                 //premiere colonne
2132                                 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] - uj_n*_Vj[1+idim];
2133                                 //colonnes intermediaires
2134                                 for(int jdim=0; jdim<_Ndim;jdim++)
2135                                         _AroeMinus[(1+idim)*_nVar + jdim + 1] = -0.5*kj*_vec_normal[idim]*_Vj[1+jdim];
2136                                 //matrice identite
2137                                 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0.5*uj_n;
2138                                 //derniere colonne
2139                                 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0.5*kj*_vec_normal[idim];
2140                         }
2141
2142                         //derniere ligne (energie)
2143                         _AroeMinus[_nVar*(_nVar-1)] = 0;
2144                         for(int idim=0; idim<_Ndim;idim++)
2145                                 _AroeMinus[_nVar*(_nVar-1)+idim+1]= 0;
2146                         _AroeMinus[_nVar*_nVar -1] = 0;
2147                 }
2148         }
2149         if(_timeScheme==Implicit)
2150                 for(int i=0; i<_nVar*_nVar;i++)
2151                 {
2152                         _AroeMinusImplicit[i] = _AroeMinus[i];
2153                         _AroePlusImplicit[i]  = _AroePlus[i];
2154                 }
2155
2156         /******** Construction de la matrices Aroe *********/
2157         /*
2158         //premiere ligne (masse)
2159         _Aroe[0]=0;
2160         for(int idim=0; idim<_Ndim;idim++)
2161                 _Aroe[1+idim]=_vec_normal[idim];
2162         _Aroe[_nVar-1]=0;
2163
2164         //lignes intermadiaires(qdm)
2165         for(int idim=0; idim<_Ndim;idim++)
2166         {
2167                 //premiere colonne
2168                 _Aroe[(1+idim)*_nVar]=Ki*_vec_normal[idim] - uj_n*_Uj[1+idim];
2169                 //colonnes intermediaires
2170                 for(int jdim=0; jdim<_Ndim;jdim++)
2171                         _Aroe[(1+idim)*_nVar + jdim + 1] = _Uj[1+idim]*_vec_normal[jdim]-ki*_vec_normal[idim]*_Ui[1+jdim];
2172                 //matrice identite
2173                 _Aroe[(1+idim)*_nVar + idim + 1] += uj_n;
2174                 //derniere colonne
2175                 _Aroe[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2176         }
2177
2178         //derniere ligne (energie)
2179         _Aroe[_nVar*(_nVar-1)] = (Ki - H)*uj_n;
2180         for(int idim=0; idim<_Ndim;idim++)
2181                 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - ki*ui_n*_Uj[idim+1];
2182         _Aroe[_nVar*_nVar -1] = (1 + ki)*uj_n;
2183          */
2184 }
2185
2186 void SinglePhase::staggeredVFFCMatricesPrimitiveVariables(double un)//vitesse normale de Roe en entrée
2187 {
2188         if(_verbose && _nbTimeStep%_freqSave ==0)
2189                 cout<<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables()"<<endl;
2190
2191         if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2192         {
2193                 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding" << endl;
2194                 _runLogFile->close();
2195                 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding");
2196         }
2197         else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2198         {
2199                 double ui_n=0., ui_2=0., uj_n=0., uj_2=0.;//vitesse normale et carré du module
2200                 double H;//enthalpie totale (expression particulière)
2201                 consToPrim(_Ui,_Vi,_porosityi);
2202                 consToPrim(_Uj,_Vj,_porosityj);
2203
2204                 for(int i=0;i<_Ndim;i++)
2205                 {
2206                         ui_2 += _Vi[1+i]*_Vi[1+i];
2207                         ui_n += _Vi[1+i]*_vec_normal[i];
2208                         uj_2 += _Vj[1+i]*_Vj[1+i];
2209                         uj_n += _Vj[1+i]*_vec_normal[i];
2210                 }
2211
2212                 if(_verbose && _nbTimeStep%_freqSave ==0){
2213                         cout <<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables " << endl;
2214                         cout<<"Vecteur primitif _Vi" << endl;
2215                         for(int i=0;i<_nVar;i++)
2216                                 cout<<_Vi[i]<<", ";
2217                         cout<<endl;
2218                         cout<<"Vecteur primitif _Vj" << endl;
2219                         for(int i=0;i<_nVar;i++)
2220                                 cout<<_Vj[i]<<", ";
2221                         cout<<endl;
2222                 }
2223
2224                 double gamma=_fluides[0]->constante("gamma");
2225                 double q=_fluides[0]->constante("q");
2226
2227                 if(fabs(un)>_precision)//non zero velocity on the interface
2228                 {
2229                         if(     !_useDellacherieEOS)
2230                         {
2231                                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2232                                 double cv=fluide0->constante("cv");
2233
2234                                 if(un>_precision)
2235                                 {
2236                                         double rhoi,rhoj,pj, Ei, ei;
2237                                         double cj;//vitesse du son pour calcul valeurs propres
2238                                         rhoi=_Ui[0]/_porosityi;
2239                                         Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2240                                         ei=Ei-0.5*ui_2;
2241                                         pj=_Vj[0];
2242                                         rhoj=_Uj[0]/_porosityj;
2243                                         cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2244
2245                                         /***********Calcul des valeurs propres ********/
2246                                         vector<complex<double> > vp_dist(3,0);
2247                                         vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2248
2249                                         _maxvploc=fabs(ui_n)+cj;
2250                                         if(_maxvploc>_maxvp)
2251                                                 _maxvp=_maxvploc;
2252
2253                                         if(_verbose && _nbTimeStep%_freqSave ==0)
2254                                                 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2255
2256                                         /******** Construction de la matrice A^+ *********/
2257                                         //premiere ligne (masse)
2258                                         _AroePlusImplicit[0]=ui_n/((gamma-1)*(ei-q));
2259                                         for(int idim=0; idim<_Ndim;idim++)
2260                                                 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2261                                         _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cv/(ei-q);
2262
2263                                         //lignes intermadiaires(qdm)
2264                                         for(int idim=0; idim<_Ndim;idim++)
2265                                         {
2266                                                 //premiere colonne
2267                                                 _AroePlusImplicit[(1+idim)*_nVar]=ui_n/((gamma-1)*(ei-q))*_Vi[1+idim];
2268                                                 //colonnes intermediaires
2269                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2270                                                         _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2271                                                 //matrice identite
2272                                                 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2273                                                 //derniere colonne
2274                                                 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cv/(ei-q)*_Vi[1+idim];
2275                                         }
2276
2277                                         //derniere ligne (energie)
2278                                         _AroePlusImplicit[_nVar*(_nVar-1)] = Ei*ui_n/((gamma-1)*(ei-q));
2279                                         for(int idim=0; idim<_Ndim;idim++)
2280                                                 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2281                                         _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Ei/(ei-q))*cv;
2282
2283                                         /******** Construction de la matrice A^- *********/
2284                                         //premiere ligne (masse)
2285                                         _AroeMinusImplicit[0]=0;
2286                                         for(int idim=0; idim<_Ndim;idim++)
2287                                                 _AroeMinusImplicit[1+idim]=0;
2288                                         _AroeMinusImplicit[_nVar-1]=0;
2289
2290                                         //lignes intermadiaires(qdm)
2291                                         for(int idim=0; idim<_Ndim;idim++)
2292                                         {
2293                                                 //premiere colonne
2294                                                 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2295                                                 //colonnes intermediaires
2296                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2297                                                         _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2298                                                 //matrice identite
2299                                                 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2300                                                 //derniere colonne
2301                                                 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2302                                         }
2303
2304                                         //derniere ligne (energie)
2305                                         _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2306                                         for(int idim=0; idim<_Ndim;idim++)
2307                                                 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2308                                         _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2309                                 }
2310                                 else if(un<-_precision)
2311                                 {
2312                                         double pi, Ej, rhoi, rhoj, ej;
2313                                         double ci;//vitesse du son pour calcul valeurs propres
2314                                         rhoj=_Uj[0]/_porosityj;
2315                                         Ej= _Uj[_Ndim+1]/rhoj;
2316                                         ej=Ej-0.5*uj_2;
2317                                         pi=_Vi[0];
2318                                         rhoi=_Ui[0]/_porosityi;
2319                                         ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2320
2321                                         /***********Calcul des valeurs propres ********/
2322                                         vector<complex<double> > vp_dist(3,0);
2323                                         vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2324
2325                                         _maxvploc=fabs(uj_n)+ci;
2326                                         if(_maxvploc>_maxvp)
2327                                                 _maxvp=_maxvploc;
2328
2329                                         if(_verbose && _nbTimeStep%_freqSave ==0)
2330                                                 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2331
2332                                         /******** Construction de la matrice A^+ *********/
2333                                         //premiere ligne (masse)
2334                                         _AroePlusImplicit[0]=0;
2335                                         for(int idim=0; idim<_Ndim;idim++)
2336                                                 _AroePlusImplicit[1+idim]=0;
2337                                         _AroePlusImplicit[_nVar-1]=0;
2338
2339                                         //lignes intermadiaires(qdm)
2340                                         for(int idim=0; idim<_Ndim;idim++)
2341                                         {
2342                                                 //premiere colonne
2343                                                 _AroePlusImplicit[(1+idim)*_nVar]=0;
2344                                                 //colonnes intermediaires
2345                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2346                                                         _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2347                                                 //matrice identite
2348                                                 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2349                                                 //derniere colonne
2350                                                 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2351                                         }
2352
2353                                         //derniere ligne (energie)
2354                                         _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2355                                         for(int idim=0; idim<_Ndim;idim++)
2356                                                 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2357                                         _AroePlusImplicit[_nVar*_nVar -1] =  0;
2358
2359                                         /******** Construction de la matrice A^- *********/
2360                                         //premiere ligne (masse)
2361                                         _AroeMinusImplicit[0]=uj_n/((gamma-1)*(ej-q));
2362                                         for(int idim=0; idim<_Ndim;idim++)
2363                                                 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2364                                         _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cv/(ej-q);
2365
2366                                         //lignes intermadiaires(qdm)
2367                                         for(int idim=0; idim<_Ndim;idim++)
2368                                         {
2369                                                 //premiere colonne
2370                                                 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n/((gamma-1)*(ej-q))*_Vj[1+idim];
2371                                                 //colonnes intermediaires
2372                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2373                                                         _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2374                                                 //matrice identite
2375                                                 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2376                                                 //derniere colonne
2377                                                 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cv/(ej-q)*_Vj[1+idim];
2378                                         }
2379
2380                                         //derniere ligne (energie)
2381                                         _AroeMinusImplicit[_nVar*(_nVar-1)] = Ej*uj_n/((gamma-1)*(ej-q));
2382                                         for(int idim=0; idim<_Ndim;idim++)
2383                                                 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2384                                         _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Ej/(ej-q))*cv;
2385                                 }
2386                                 else
2387                                 {
2388                                         *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2389                                         _runLogFile->close();
2390                                         throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2391                                 }
2392                         }
2393                         else if(_useDellacherieEOS )
2394                         {
2395                                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2396                                 double cp=fluide0->constante("cp");
2397
2398                                 if(un>_precision)
2399                                 {
2400                                         double rhoi,rhoj,pj, Ei, hi, Hi;
2401                                         double cj;//vitesse du son pour calcul valeurs propres
2402                                         rhoi=_Ui[0]/_porosityi;
2403                                         Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2404                                         Hi=Ei+_Vi[0]/rhoi;
2405                                         hi=Ei-0.5*ui_2;
2406                                         pj=_Vj[0];
2407                                         rhoj=_Uj[0]/_porosityj;
2408                                         cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2409
2410                                         /***********Calcul des valeurs propres ********/
2411                                         vector<complex<double> > vp_dist(3,0);
2412                                         vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2413
2414                                         _maxvploc=fabs(ui_n)+cj;
2415                                         if(_maxvploc>_maxvp)
2416                                                 _maxvp=_maxvploc;
2417
2418                                         if(_verbose && _nbTimeStep%_freqSave ==0)
2419                                                 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2420
2421                                         /******** Construction de la matrice A^+ *********/
2422                                         //premiere ligne (masse)
2423                                         _AroePlusImplicit[0]=ui_n*gamma/((gamma-1)*(hi-q));
2424                                         for(int idim=0; idim<_Ndim;idim++)
2425                                                 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2426                                         _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cp/(hi-q);
2427
2428                                         //lignes intermadiaires(qdm)
2429                                         for(int idim=0; idim<_Ndim;idim++)
2430                                         {
2431                                                 //premiere colonne
2432                                                 _AroePlusImplicit[(1+idim)*_nVar]=ui_n*gamma/((gamma-1)*(hi-q))*_Vi[1+idim];
2433                                                 //colonnes intermediaires
2434                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2435                                                         _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2436                                                 //matrice identite
2437                                                 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2438                                                 //derniere colonne
2439                                                 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cp/(hi-q)*_Vi[1+idim];
2440                                         }
2441
2442                                         //derniere ligne (energie)
2443                                         _AroePlusImplicit[_nVar*(_nVar-1)] = ui_n*(Hi*gamma/((gamma-1)*(hi-q))-1);
2444                                         for(int idim=0; idim<_Ndim;idim++)
2445                                                 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2446                                         _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Hi/(hi-q))*cp;
2447
2448                                         /******** Construction de la matrice A^- *********/
2449                                         //premiere ligne (masse)
2450                                         _AroeMinusImplicit[0]=0;
2451                                         for(int idim=0; idim<_Ndim;idim++)
2452                                                 _AroeMinusImplicit[1+idim]=0;
2453                                         _AroeMinusImplicit[_nVar-1]=0;
2454
2455                                         //lignes intermadiaires(qdm)
2456                                         for(int idim=0; idim<_Ndim;idim++)
2457                                         {
2458                                                 //premiere colonne
2459                                                 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2460                                                 //colonnes intermediaires
2461                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2462                                                         _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2463                                                 //matrice identite
2464                                                 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2465                                                 //derniere colonne
2466                                                 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2467                                         }
2468
2469                                         //derniere ligne (energie)
2470                                         _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2471                                         for(int idim=0; idim<_Ndim;idim++)
2472                                                 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2473                                         _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2474                                 }
2475                                 else if(un<-_precision)
2476                                 {
2477                                         double pi, Ej, rhoi,rhoj, Hj, hj;
2478                                         double ci;//vitesse du son pour calcul valeurs propres
2479                                         rhoj=_Uj[0]/_porosityj;
2480                                         Ej= _Uj[_Ndim+1]/rhoj;
2481                                         Hj=Ej+_Vj[0]/rhoj;
2482                                         hj=Ej-0.5*uj_2;
2483                                         pi=_Vi[0];
2484                                         rhoi=_Ui[0]/_porosityi;
2485                                         ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2486
2487                                         /***********Calcul des valeurs propres ********/
2488                                         vector<complex<double> > vp_dist(3,0);
2489                                         vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2490
2491                                         _maxvploc=fabs(uj_n)+ci;
2492                                         if(_maxvploc>_maxvp)
2493                                                 _maxvp=_maxvploc;
2494
2495                                         if(_verbose && _nbTimeStep%_freqSave ==0)
2496                                                 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2497
2498                                         /******** Construction de la matrice A^+ *********/
2499                                         //premiere ligne (masse)
2500                                         _AroePlusImplicit[0]=0;
2501                                         for(int idim=0; idim<_Ndim;idim++)
2502                                                 _AroePlusImplicit[1+idim]=0;
2503                                         _AroePlusImplicit[_nVar-1]=0;
2504
2505                                         //lignes intermadiaires(qdm)
2506                                         for(int idim=0; idim<_Ndim;idim++)
2507                                         {
2508                                                 //premiere colonne
2509                                                 _AroePlusImplicit[(1+idim)*_nVar]=0;
2510                                                 //colonnes intermediaires
2511                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2512                                                         _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2513                                                 //matrice identite
2514                                                 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2515                                                 //derniere colonne
2516                                                 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2517                                         }
2518
2519                                         //derniere ligne (energie)
2520                                         _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2521                                         for(int idim=0; idim<_Ndim;idim++)
2522                                                 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2523                                         _AroePlusImplicit[_nVar*_nVar -1] =  0;
2524
2525                                         /******** Construction de la matrice A^- *********/
2526                                         //premiere ligne (masse)
2527                                         _AroeMinusImplicit[0]=uj_n*gamma/((gamma-1)*(hj-q));
2528                                         for(int idim=0; idim<_Ndim;idim++)
2529                                                 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2530                                         _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cp/(hj-q);
2531
2532                                         //lignes intermadiaires(qdm)
2533                                         for(int idim=0; idim<_Ndim;idim++)
2534                                         {
2535                                                 //premiere colonne
2536                                                 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n*gamma/((gamma-1)*(hj-q))*_Vj[1+idim];
2537                                                 //colonnes intermediaires
2538                                                 for(int jdim=0; jdim<_Ndim;jdim++)
2539                                                         _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2540                                                 //matrice identite
2541                                                 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2542                                                 //derniere colonne
2543                                                 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cp/(hj-q)*_Vj[1+idim];
2544                                         }
2545
2546                                         //derniere ligne (energie)
2547                                         _AroeMinusImplicit[_nVar*(_nVar-1)] = uj_n*(Hj*gamma/((gamma-1)*(hj-q))-1);
2548                                         for(int idim=0; idim<_Ndim;idim++)
2549                                                 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2550                                         _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Hj/(hj-q))*cp;
2551                                 }
2552                                 else
2553                                 {
2554                                         *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2555                                         _runLogFile->close();
2556                                         throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2557                                 }
2558                         }
2559                         else
2560                         {
2561                                 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2562                                 _runLogFile->close();
2563                                 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2564                         }
2565                 }
2566                 else//case nil velocity on the interface, apply centered scheme
2567                 {
2568                         primToConsJacobianMatrix(_Vj);
2569                         Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
2570                         primToConsJacobianMatrix(_Vi);
2571                         Polynoms::matrixProduct(_AroePlus,  _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
2572                 }
2573         }
2574 }
2575 void SinglePhase::applyVFRoeLowMachCorrections(bool isBord, string groupname)
2576 {
2577         if(_nonLinearFormulation!=VFRoe)
2578         {
2579                 *_runLogFile<< "SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation" << endl;
2580                 _runLogFile->close();
2581                 throw CdmathException("SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
2582         }
2583         else//_nonLinearFormulation==VFRoe
2584         {
2585                 if(_spaceScheme==lowMach){
2586                         double u_2=0;
2587                         for(int i=0;i<_Ndim;i++)
2588                                 u_2 += _Uroe[1+i]*_Uroe[1+i];
2589                         double  c = _maxvploc;//vitesse du son a l'interface
2590                         double M=sqrt(u_2)/c;//Mach number
2591                         _Vij[0]=M*_Vij[0]+(1-M)*(_Vi[0]+_Vj[0])/2;
2592                 }
2593                 else if(_spaceScheme==pressureCorrection)
2594                 {//order 1 : no correction, oarder 2 : correction everywhere, order 3 : correction only inside, orders 4 and 5 : special correction at boundaries
2595                         if(_pressureCorrectionOrder==2 || (!isBord && _pressureCorrectionOrder==3) || (!isBord && _pressureCorrectionOrder==4) || (!isBord && _pressureCorrectionOrder==5) )
2596                         {
2597                                 double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;
2598                                 for(int i=0;i<_Ndim;i++)
2599                                 {
2600                                         norm_uij += _Uroe[1+i]*_Uroe[1+i];
2601                                         uij_n += _Uroe[1+i]*_vec_normal[i];
2602                                         ui_n += _Vi[1+i]*_vec_normal[i];
2603                                         uj_n += _Vj[1+i]*_vec_normal[i];
2604                                 }
2605                                 norm_uij=sqrt(norm_uij);
2606                                 if(norm_uij>_precision)//avoid division by zero
2607                                         _Vij[0]=(_Vi[0]+_Vj[0])/2 + uij_n/norm_uij*(_Vj[0]-_Vi[0])/4 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
2608                                 else
2609                                         _Vij[0]=(_Vi[0]+_Vj[0])/2                                    - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
2610                         }
2611                         else if(_pressureCorrectionOrder==4 && isBord)
2612                                 _Vij[0]=_Vi[0];
2613                         else if(_pressureCorrectionOrder==5 && isBord)
2614                         {
2615                                 double g_n=0;//scalar product of gravity and normal vector
2616                                 for(int i=0;i<_Ndim;i++)
2617                                         g_n += _GravityField3d[i]*_vec_normal[i];
2618                                 _Vij[0]=_Vi[0]- _Ui[0]*g_n/_inv_dxi/2;
2619                         }
2620                 }
2621                 else if(_spaceScheme==staggered)
2622                 {
2623                         double uij_n=0;
2624                         for(int i=0;i<_Ndim;i++)
2625                                 uij_n += _Uroe[1+i]*_vec_normal[i];
2626
2627                         if(uij_n>_precision){
2628                                 _Vij[0]=_Vj[0];
2629                                 for(int i=0;i<_Ndim;i++)
2630                                         _Vij[1+i]=_Vi[1+i];
2631                                 _Vij[_nVar-1]=_Vi[_nVar-1];
2632                         }
2633                         else if(uij_n<-_precision){
2634                                 _Vij[0]=_Vi[0];
2635                                 for(int i=0;i<_Ndim;i++)
2636                                         _Vij[1+i]=_Vj[1+i];
2637                                 _Vij[_nVar-1]=_Vj[_nVar-1];
2638                         }
2639                         else{
2640                                 _Vij[0]=(_Vi[0]+_Vi[0])/2;
2641                                 for(int i=0;i<_Ndim;i++)
2642                                         _Vij[1+i]=(_Vj[1+i]+_Vj[1+i])/2;
2643                                 _Vij[_nVar-1]=(_Vj[_nVar-1]+_Vj[_nVar-1])/2;
2644                         }
2645                 }
2646                 primToCons(_Vij,0,_Uij,0);
2647         }
2648 }
2649
2650 void SinglePhase::testConservation()
2651 {
2652         double SUM, DELTA, x;
2653         int I;
2654         for(int i=0; i<_nVar; i++)
2655         {
2656                 {
2657                         if(i == 0)
2658                                 cout << "Masse totale (kg): ";
2659                         else
2660                         {
2661                                 if(i == _nVar-1)
2662                                         cout << "Energie totale (J): ";
2663                                 else
2664                                         cout << "Quantite de mouvement totale (kg.m.s^-1): ";
2665                         }
2666                 }
2667                 SUM = 0;
2668                 I =  i;
2669                 DELTA = 0;
2670                 for(int j=0; j<_Nmailles; j++)
2671                 {
2672                         if(!_usePrimitiveVarsInNewton)
2673                                 VecGetValues(_conservativeVars, 1, &I, &x);//on recupere la valeur du champ
2674                         else
2675                                 VecGetValues(_primitiveVars, 1, &I, &x);//on recupere la valeur du champ
2676                         SUM += x*_mesh.getCell(j).getMeasure();
2677                         VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
2678                         DELTA += x*_mesh.getCell(j).getMeasure();
2679                         I += _nVar;
2680                 }
2681                 if(fabs(SUM)>_precision)
2682                         cout << SUM << ", variation relative: " << fabs(DELTA /SUM)  << endl;
2683                 else
2684                         cout << " a une somme quasi nulle,  variation absolue: " << fabs(DELTA) << endl;
2685         }
2686 }
2687
2688 void SinglePhase::getDensityDerivatives( double pressure, double temperature, double v2)
2689 {
2690         double rho=_fluides[0]->getDensity(pressure,temperature);
2691         double gamma=_fluides[0]->constante("gamma");
2692         double q=_fluides[0]->constante("q");
2693
2694         if(     !_useDellacherieEOS)
2695         {
2696                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2697                 double e = fluide0->getInternalEnergy(temperature);
2698                 double cv=fluide0->constante("cv");
2699                 double E=e+0.5*v2;
2700
2701                 _drho_sur_dp=1/((gamma-1)*(e-q));
2702                 _drho_sur_dT=-rho*cv/(e-q);
2703                 _drhoE_sur_dp=E/((gamma-1)*(e-q));
2704                 _drhoE_sur_dT=rho*cv*(1-E/(e-q));
2705         }
2706         else if(_useDellacherieEOS )
2707         {
2708                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2709                 double h=fluide0->getEnthalpy(temperature);
2710                 double H=h+0.5*v2;
2711                 double cp=fluide0->constante("cp");
2712
2713                 _drho_sur_dp=gamma/((gamma-1)*(h-q));
2714                 _drho_sur_dT=-rho*cp/(h-q);
2715                 _drhoE_sur_dp=gamma*H/((gamma-1)*(h-q))-1;
2716                 _drhoE_sur_dT=rho*cp*(1-H/(h-q));
2717         }
2718         else
2719         {
2720                 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2721                 _runLogFile->close();
2722                 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2723         }
2724
2725         if(_verbose && _nbTimeStep%_freqSave ==0)
2726         {
2727                 cout<<"_drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
2728                 cout<<"_drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
2729         }
2730 }
2731 void SinglePhase::save(){
2732         string prim(_path+"/SinglePhasePrim_");///Results
2733         string cons(_path+"/SinglePhaseCons_");
2734         string allFields(_path+"/");
2735         prim+=_fileName;
2736         cons+=_fileName;
2737         allFields+=_fileName;
2738
2739         PetscInt Ii;
2740         for (long i = 0; i < _Nmailles; i++){
2741                 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2742                 for (int j = 0; j < _nVar; j++){
2743                         Ii = i*_nVar +j;
2744                         VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
2745                 }
2746         }
2747         if(_saveConservativeField){
2748                 for (long i = 0; i < _Nmailles; i++){
2749                         // j = 0 : density; j = _nVar - 1 : energy j = 1,..,_nVar-2: momentum
2750                         for (int j = 0; j < _nVar; j++){
2751                                 Ii = i*_nVar +j;
2752                                 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
2753                         }
2754                 }
2755                 _UU.setTime(_time,_nbTimeStep);
2756         }
2757         _VV.setTime(_time,_nbTimeStep);
2758
2759         // create mesh and component info
2760         if (_nbTimeStep ==0 || _restartWithNewFileName){
2761                 string prim_suppress ="rm -rf "+prim+"_*";
2762                 string cons_suppress ="rm -rf "+cons+"_*";
2763
2764                 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
2765                 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
2766
2767                 if(_saveConservativeField){
2768                         _UU.setInfoOnComponent(0,"Density_(kg/m^3)");
2769                         _UU.setInfoOnComponent(1,"Momentum_x");// (kg/m^2/s)
2770                         if (_Ndim>1)
2771                                 _UU.setInfoOnComponent(2,"Momentum_y");// (kg/m^2/s)
2772                         if (_Ndim>2)
2773                                 _UU.setInfoOnComponent(3,"Momentum_z");// (kg/m^2/s)
2774
2775                         _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
2776
2777                         switch(_saveFormat)
2778                         {
2779                         case VTK :
2780                                 _UU.writeVTK(cons);
2781                                 break;
2782                         case MED :
2783                                 _UU.writeMED(cons);
2784                                 break;
2785                         case CSV :
2786                                 _UU.writeCSV(cons);
2787                                 break;
2788                         }
2789                 }
2790                 _VV.setInfoOnComponent(0,"Pressure_(Pa)");
2791                 _VV.setInfoOnComponent(1,"Velocity_x_(m/s)");
2792                 if (_Ndim>1)
2793                         _VV.setInfoOnComponent(2,"Velocity_y_(m/s)");
2794                 if (_Ndim>2)
2795                         _VV.setInfoOnComponent(3,"Velocity_z_(m/s)");
2796                 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
2797
2798                 switch(_saveFormat)
2799                 {
2800                 case VTK :
2801                         _VV.writeVTK(prim);
2802                         break;
2803                 case MED :
2804                         _VV.writeMED(prim);
2805                         break;
2806                 case CSV :
2807                         _VV.writeCSV(prim);
2808                         break;
2809                 }
2810
2811         }
2812         // do not create mesh
2813         else{
2814                 switch(_saveFormat)
2815                 {
2816                 case VTK :
2817                         _VV.writeVTK(prim,false);
2818                         break;
2819                 case MED :
2820                         _VV.writeMED(prim,false);
2821                         break;
2822                 case CSV :
2823                         _VV.writeCSV(prim);
2824                         break;
2825                 }
2826                 if(_saveConservativeField){
2827                         switch(_saveFormat)
2828                         {
2829                         case VTK :
2830                                 _UU.writeVTK(cons,false);
2831                                 break;
2832                         case MED :
2833                                 _UU.writeMED(cons,false);
2834                                 break;
2835                         case CSV :
2836                                 _UU.writeCSV(cons);
2837                                 break;
2838                         }
2839                 }
2840         }
2841         if(_saveVelocity || _saveAllFields){
2842                 for (long i = 0; i < _Nmailles; i++){
2843                         // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2844                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
2845                         {
2846                                 int Ii = i*_nVar +1+j;
2847                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
2848                         }
2849                         for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
2850                                 _Vitesse(i,j)=0;
2851                 }
2852                 _Vitesse.setTime(_time,_nbTimeStep);
2853                 if (_nbTimeStep ==0 || _restartWithNewFileName){
2854                         _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
2855                         _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
2856                         _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
2857
2858                         switch(_saveFormat)
2859                         {
2860                         case VTK :
2861                                 _Vitesse.writeVTK(prim+"_Velocity");
2862                                 break;
2863                         case MED :
2864                                 _Vitesse.writeMED(prim+"_Velocity");
2865                                 break;
2866                         case CSV :
2867                                 _Vitesse.writeCSV(prim+"_Velocity");
2868                                 break;
2869                         }
2870                 }
2871                 else{
2872                         switch(_saveFormat)
2873                         {
2874                         case VTK :
2875                                 _Vitesse.writeVTK(prim+"_Velocity",false);
2876                                 break;
2877                         case MED :
2878                                 _Vitesse.writeMED(prim+"_Velocity",false);
2879                                 break;
2880                         case CSV :
2881                                 _Vitesse.writeCSV(prim+"_Velocity");
2882                                 break;
2883                         }
2884                 }
2885         }
2886
2887         if(_saveAllFields)
2888         {
2889                 double p,T,rho, h, vx,vy,vz,v2;
2890                 int Ii;
2891                 for (long i = 0; i < _Nmailles; i++){
2892                         Ii = i*_nVar;
2893                         VecGetValues(_conservativeVars,1,&Ii,&rho);
2894                         Ii = i*_nVar;
2895                         VecGetValues(_primitiveVars,1,&Ii,&p);
2896                         Ii = i*_nVar +_nVar-1;
2897                         VecGetValues(_primitiveVars,1,&Ii,&T);
2898                         Ii = i*_nVar + 1;
2899                         VecGetValues(_primitiveVars,1,&Ii,&vx);
2900                         if(_Ndim>1)
2901                         {
2902                                 Ii = i*_nVar + 2;
2903                                 VecGetValues(_primitiveVars,1,&Ii,&vy);
2904                                 if(_Ndim>2){
2905                                         Ii = i*_nVar + 3;
2906                                         VecGetValues(_primitiveVars,1,&Ii,&vz);
2907                                 }
2908                         }
2909
2910                         h   = _fluides[0]->getEnthalpy(T,rho);
2911
2912                         _Enthalpy(i)=h;
2913                         _Density(i)=rho;
2914                         _Pressure(i)=p;
2915                         _Temperature(i)=T;
2916                         _VitesseX(i)=vx;
2917                         v2=vx*vx;
2918                         if(_Ndim>1)
2919                         {
2920                                 _VitesseY(i)=vy;
2921                                 v2+=vy*vy;
2922                                 if(_Ndim>2)
2923                                 {
2924                                         _VitesseZ(i)=vz;
2925                                         v2+=vz*vz;
2926                                 }
2927                         }
2928                         _MachNumber(i)=sqrt(v2)/_fluides[0]->vitesseSonEnthalpie(h);
2929                 }
2930                 _Enthalpy.setTime(_time,_nbTimeStep);
2931                 _Density.setTime(_time,_nbTimeStep);
2932                 _Pressure.setTime(_time,_nbTimeStep);
2933                 _Temperature.setTime(_time,_nbTimeStep);
2934                 _MachNumber.setTime(_time,_nbTimeStep);
2935                 _VitesseX.setTime(_time,_nbTimeStep);
2936                 if(_Ndim>1)
2937                 {
2938                         _VitesseY.setTime(_time,_nbTimeStep);
2939                         if(_Ndim>2)
2940                                 _VitesseZ.setTime(_time,_nbTimeStep);
2941                 }
2942                 if (_nbTimeStep ==0 || _restartWithNewFileName){
2943                         switch(_saveFormat)
2944                         {
2945                         case VTK :
2946                                 _Enthalpy.writeVTK(allFields+"_Enthalpy");
2947                                 _Density.writeVTK(allFields+"_Density");
2948                                 _Pressure.writeVTK(allFields+"_Pressure");
2949                                 _Temperature.writeVTK(allFields+"_Temperature");
2950                                 _MachNumber.writeVTK(allFields+"_MachNumber");
2951                                 _VitesseX.writeVTK(allFields+"_VelocityX");
2952                                 if(_Ndim>1)
2953                                 {
2954                                         _VitesseY.writeVTK(allFields+"_VelocityY");
2955                                         if(_Ndim>2)
2956                                                 _VitesseZ.writeVTK(allFields+"_VelocityZ");
2957                                 }
2958                                 break;
2959                         case MED :
2960                                 _Enthalpy.writeMED(allFields+"_Enthalpy");
2961                                 _Density.writeMED(allFields+"_Density");
2962                                 _Pressure.writeMED(allFields+"_Pressure");
2963                                 _Temperature.writeMED(allFields+"_Temperature");
2964                                 _MachNumber.writeMED(allFields+"_MachNumber");
2965                                 _VitesseX.writeMED(allFields+"_VelocityX");
2966                                 if(_Ndim>1)
2967                                 {
2968                                         _VitesseY.writeMED(allFields+"_VelocityY");
2969                                         if(_Ndim>2)
2970                                                 _VitesseZ.writeMED(allFields+"_VelocityZ");
2971                                 }
2972                                 break;
2973                         case CSV :
2974                                 _Enthalpy.writeCSV(allFields+"_Enthalpy");
2975                                 _Density.writeCSV(allFields+"_Density");
2976                                 _Pressure.writeCSV(allFields+"_Pressure");
2977                                 _Temperature.writeCSV(allFields+"_Temperature");
2978                                 _MachNumber.writeCSV(allFields+"_MachNumber");
2979                                 _VitesseX.writeCSV(allFields+"_VelocityX");
2980                                 if(_Ndim>1)
2981                                 {
2982                                         _VitesseY.writeCSV(allFields+"_VelocityY");
2983                                         if(_Ndim>2)
2984                                                 _VitesseZ.writeCSV(allFields+"_VelocityZ");
2985                                 }
2986                                 break;
2987                         }
2988                 }
2989                 else{
2990                         switch(_saveFormat)
2991                         {
2992                         case VTK :
2993                                 _Enthalpy.writeVTK(allFields+"_Enthalpy",false);
2994                                 _Density.writeVTK(allFields+"_Density",false);
2995                                 _Pressure.writeVTK(allFields+"_Pressure",false);
2996                                 _Temperature.writeVTK(allFields+"_Temperature",false);
2997                                 _MachNumber.writeVTK(allFields+"_MachNumber",false);
2998                                 _VitesseX.writeVTK(allFields+"_VelocityX",false);
2999                                 if(_Ndim>1)
3000                                 {
3001                                         _VitesseY.writeVTK(allFields+"_VelocityY",false);
3002                                         if(_Ndim>2)
3003                                                 _VitesseZ.writeVTK(allFields+"_VelocityZ",false);
3004                                 }
3005                                 break;
3006                         case MED :
3007                                 _Enthalpy.writeMED(allFields+"_Enthalpy",false);
3008                                 _Density.writeMED(allFields+"_Density",false);
3009                                 _Pressure.writeMED(allFields+"_Pressure",false);
3010                                 _Temperature.writeMED(allFields+"_Temperature",false);
3011                                 _MachNumber.writeMED(allFields+"_MachNumber",false);
3012                                 _VitesseX.writeMED(allFields+"_VelocityX",false);
3013                                 if(_Ndim>1)
3014                                 {
3015                                         _VitesseY.writeMED(allFields+"_VelocityY",false);
3016                                         if(_Ndim>2)
3017                                                 _VitesseZ.writeMED(allFields+"_VelocityZ",false);
3018                                 }
3019                                 break;
3020                         case CSV :
3021                                 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3022                                 _Density.writeCSV(allFields+"_Density");
3023                                 _Pressure.writeCSV(allFields+"_Pressure");
3024                                 _Temperature.writeCSV(allFields+"_Temperature");
3025                                 _MachNumber.writeCSV(allFields+"_MachNumber");
3026                                 _VitesseX.writeCSV(allFields+"_VelocityX");
3027                                 if(_Ndim>1)
3028                                 {
3029                                         _VitesseY.writeCSV(allFields+"_VelocityY");
3030                                         if(_Ndim>2)
3031                                                 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3032                                 }
3033                                 break;
3034                         }
3035                 }
3036         }
3037
3038         if(_isStationary)
3039         {
3040                 prim+="_Stat";
3041                 cons+="_Stat";
3042
3043                 switch(_saveFormat)
3044                 {
3045                 case VTK :
3046                         _VV.writeVTK(prim);
3047                         break;
3048                 case MED :
3049                         _VV.writeMED(prim);
3050                         break;
3051                 case CSV :
3052                         _VV.writeCSV(prim);
3053                         break;
3054                 }
3055
3056                 if(_saveConservativeField){
3057                         switch(_saveFormat)
3058                         {
3059                         case VTK :
3060                                 _UU.writeVTK(cons);
3061                                 break;
3062                         case MED :
3063                                 _UU.writeMED(cons);
3064                                 break;
3065                         case CSV :
3066                                 _UU.writeCSV(cons);
3067                                 break;
3068                         }
3069                 }
3070
3071                 if(_saveVelocity || _saveAllFields){
3072                         switch(_saveFormat)
3073                         {
3074                         case VTK :
3075                                 _Vitesse.writeVTK(prim+"_Velocity");
3076                                 break;
3077                         case MED :
3078                                 _Vitesse.writeMED(prim+"_Velocity");
3079                                 break;
3080                         case CSV :
3081                                 _Vitesse.writeCSV(prim+"_Velocity");
3082                                 break;
3083                         }
3084                 }
3085         }
3086
3087         if (_restartWithNewFileName)
3088                 _restartWithNewFileName=false;
3089 }
3090
3091 Field& SinglePhase::getPressureField()
3092 {
3093         if(!_initializedMemory)
3094                 throw CdmathException("SinglePhase::getPressureField, Call initialize first");
3095
3096         if(!_saveAllFields)
3097         {
3098                 _Pressure=Field("Pressure",CELLS,_mesh,1);
3099                 int Ii;
3100                 for (long i = 0; i < _Nmailles; i++){
3101                         Ii = i*_nVar;
3102                         VecGetValues(_primitiveVars,1,&Ii,&_Pressure(i));
3103                 }
3104                 _Pressure.setTime(_time,_nbTimeStep);
3105         }
3106         return _Pressure;
3107 }
3108
3109 Field& SinglePhase::getTemperatureField()
3110 {
3111         if(!_initializedMemory)
3112                 throw CdmathException("SinglePhase::getTemperatureField, Call initialize first");
3113
3114         if(!_saveAllFields)
3115         {
3116                 _Temperature=Field("Temperature",CELLS,_mesh,1);
3117                 int Ii;
3118                 for (long i = 0; i < _Nmailles; i++){
3119                         Ii = i*_nVar +_nVar-1;
3120                         VecGetValues(_primitiveVars,1,&Ii,&_Temperature(i));
3121                 }
3122                 _Temperature.setTime(_time,_nbTimeStep);
3123         }
3124         return _Temperature;
3125 }
3126
3127 Field& SinglePhase::getVelocityField()
3128 {
3129         if(!_initializedMemory)
3130                 throw CdmathException("SinglePhase::getVelocityField, Call initialize first");
3131
3132         if(!_saveAllFields )
3133         {
3134                 _Vitesse=Field("Vitesse",CELLS,_mesh,3);
3135                 int Ii;
3136                 for (long i = 0; i < _Nmailles; i++)
3137                 {
3138                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
3139                         {
3140                                 int Ii = i*_nVar +1+j;
3141                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
3142                         }
3143                         for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
3144                                 _Vitesse(i,j)=0;
3145                 }
3146                 _Vitesse.setTime(_time,_nbTimeStep);
3147                 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
3148                 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
3149                 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
3150         }
3151
3152         return _Vitesse;
3153 }
3154
3155 Field& SinglePhase::getMachNumberField()
3156 {
3157         if(!_initializedMemory)
3158                 throw CdmathException("SinglePhase::getMachNumberField, Call initialize first");
3159
3160         if(!_saveAllFields )
3161         {
3162                 _MachNumber=Field("Mach number",CELLS,_mesh,1);
3163                 int Ii;
3164                 double p,T,rho,h, temp, u2=0;
3165                 for (long i = 0; i < _Nmailles; i++){
3166                         Ii = i*_nVar;
3167                         VecGetValues(_primitiveVars,1,&Ii,&p);
3168                         Ii = i*_nVar +_nVar-1;
3169                         VecGetValues(_primitiveVars,1,&Ii,&T);
3170
3171                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
3172                         {
3173                                 int Ii = i*_nVar +1+j;
3174                                 VecGetValues(_primitiveVars,1,&Ii,&temp);
3175                                 u2+=temp*temp;
3176                         }
3177
3178                         rho=_fluides[0]->getDensity(p,T);
3179                         h  =_fluides[0]->getEnthalpy(T,rho);
3180                         _MachNumber[i]  =sqrt(u2)/_fluides[0]->vitesseSonEnthalpie(h);
3181                         //cout<<"u="<<sqrt(u2)<<", c= "<<_fluides[0]->vitesseSonEnthalpie(h)<<", MachNumberField[i] = "<<MachNumberField[i] <<endl;
3182                 }
3183                 _MachNumber.setTime(_time,_nbTimeStep);
3184         }
3185         //cout<<", MachNumberField = "<<MachNumberField <<endl;
3186
3187         return _MachNumber;
3188 }
3189
3190 Field& SinglePhase::getVelocityXField()
3191 {
3192         if(!_initializedMemory)
3193                 throw CdmathException("SinglePhase::getVelocityXField, Call initialize first");
3194
3195         if(!_saveAllFields )
3196         {
3197                 _VitesseX=Field("Velocity X",CELLS,_mesh,1);
3198                 int Ii;
3199                 for (long i = 0; i < _Nmailles; i++)
3200                 {
3201                         int Ii = i*_nVar +1;
3202                         VecGetValues(_primitiveVars,1,&Ii,&_VitesseX(i));
3203                 }
3204                 _VitesseX.setTime(_time,_nbTimeStep);
3205                 _VitesseX.setInfoOnComponent(0,"Velocity_x_(m/s)");
3206         }
3207
3208         return _VitesseX;
3209 }
3210
3211 Field& SinglePhase::getVelocityYField()
3212 {
3213         if(!_initializedMemory)
3214                 throw CdmathException("SinglePhase::getVelocityYField, Call initialize first");
3215
3216         if(_Ndim<2)
3217         throw CdmathException("SinglePhase::getVelocityYField() error : dimension should be at least 2");
3218         else
3219                 if(!_saveAllFields )
3220                 {
3221                         _VitesseY=Field("Velocity Y",CELLS,_mesh,1);
3222                         int Ii;
3223                         for (long i = 0; i < _Nmailles; i++)
3224                         {
3225                                 int Ii = i*_nVar +2;
3226                                 VecGetValues(_primitiveVars,1,&Ii,&_VitesseY(i));
3227                         }
3228                         _VitesseY.setTime(_time,_nbTimeStep);
3229                         _VitesseY.setInfoOnComponent(0,"Velocity_y_(m/s)");
3230                 }
3231
3232                 return _VitesseY;
3233 }
3234
3235 Field& SinglePhase::getVelocityZField()
3236 {
3237         if(!_initializedMemory)
3238                 throw CdmathException("SinglePhase::getVelocityZField, Call initialize first");
3239
3240         if(_Ndim<3)
3241         throw CdmathException("SinglePhase::getvelocityZField() error : dimension should be 3");
3242         else
3243                 if(!_saveAllFields )
3244                 {
3245                         _VitesseZ=Field("Velocity Z",CELLS,_mesh,1);
3246                         int Ii;
3247                         for (long i = 0; i < _Nmailles; i++)
3248                         {
3249                                 int Ii = i*_nVar +3;
3250                                 VecGetValues(_primitiveVars,1,&Ii,&_VitesseZ(i));
3251                         }
3252                         _VitesseZ.setTime(_time,_nbTimeStep);
3253                         _VitesseZ.setInfoOnComponent(0,"Velocity_z_(m/s)");
3254                 }
3255
3256                 return _VitesseZ;
3257 }
3258
3259 Field& SinglePhase::getDensityField()
3260 {
3261         if(!_initializedMemory)
3262                 throw CdmathException("SinglePhase::getDensityField, Call initialize first");
3263
3264         if(!_saveAllFields )
3265         {
3266                 _Density=Field("Density",CELLS,_mesh,1);
3267                 int Ii;
3268                 for (long i = 0; i < _Nmailles; i++){
3269                         Ii = i*_nVar;
3270                         VecGetValues(_conservativeVars,1,&Ii,&_Density(i));
3271                 }
3272                 _Density.setTime(_time,_nbTimeStep);
3273         }
3274         return _Density;
3275 }
3276
3277 Field& SinglePhase::getMomentumField()//not yet managed by parameter _saveAllFields
3278 {
3279         if(!_initializedMemory)
3280                 throw CdmathException("SinglePhase::getMomentumField, Call initialize first");
3281
3282         _Momentum=Field("Momentum",CELLS,_mesh,_Ndim);
3283         int Ii;
3284         for (long i = 0; i < _Nmailles; i++)
3285                 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de qdm
3286                 {
3287                         int Ii = i*_nVar +1+j;
3288                         VecGetValues(_conservativeVars,1,&Ii,&_Momentum(i,j));
3289                 }
3290         _Momentum.setTime(_time,_nbTimeStep);
3291
3292         return _Momentum;
3293 }
3294
3295 Field& SinglePhase::getTotalEnergyField()//not yet managed by parameter _saveAllFields
3296 {
3297         if(!_initializedMemory)
3298                 throw CdmathException("SinglePhase::getTotalEnergyField, Call initialize first");
3299
3300         _TotalEnergy=Field("TotalEnergy",CELLS,_mesh,1);
3301         int Ii;
3302         for (long i = 0; i < _Nmailles; i++){
3303                 Ii = i*_nVar +_nVar-1;
3304                 VecGetValues(_conservativeVars,1,&Ii,&_TotalEnergy(i));
3305         }
3306         _TotalEnergy.setTime(_time,_nbTimeStep);
3307
3308         return _TotalEnergy;
3309 }
3310
3311 Field& SinglePhase::getEnthalpyField()
3312 {
3313         if(!_initializedMemory)
3314                 throw CdmathException("SinglePhase::getEnthalpyField, Call initialize first");
3315
3316         if(!_saveAllFields )
3317         {
3318                 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
3319                 int Ii;
3320                 double p,T,rho;
3321                 for (long i = 0; i < _Nmailles; i++){
3322                         Ii = i*_nVar;
3323                         VecGetValues(_primitiveVars,1,&Ii,&p);
3324                         Ii = i*_nVar +_nVar-1;
3325                         VecGetValues(_primitiveVars,1,&Ii,&T);
3326
3327                         rho=_fluides[0]->getDensity(p,T);
3328                         _Enthalpy(i)=_fluides[0]->getEnthalpy(T,rho);
3329                 }
3330                 _Enthalpy.setTime(_time,_nbTimeStep);
3331         }
3332
3333         return _Enthalpy;
3334 }
3335
3336 vector<string> SinglePhase::getOutputFieldsNames()
3337 {
3338         vector<string> result(8);
3339
3340         result[0]="Pressure";
3341         result[1]="Velocity";
3342         result[2]="Temperature";
3343         result[3]="Density";
3344         result[4]="Momentum";
3345         result[5]="TotalEnergy";
3346         result[6]="Enthalpy";
3347         result[7]="VelocityX";
3348
3349         return result;
3350 }
3351
3352 Field& SinglePhase::getOutputField(const string& nameField )
3353 {
3354         if(nameField=="pressure" || nameField=="Pressure" || nameField=="PRESSURE" || nameField=="PRESSION" || nameField=="Pression"  || nameField=="pression" )
3355                 return getPressureField();
3356         else if(nameField=="velocity" || nameField=="Velocity" || nameField=="VELOCITY" || nameField=="Vitesse" || nameField=="VITESSE" || nameField=="vitesse" )
3357                 return getVelocityField();
3358         else if(nameField=="velocityX" || nameField=="VelocityX" || nameField=="VELOCITYX" || nameField=="VitesseX" || nameField=="VITESSEX" || nameField=="vitesseX" )
3359                 return getVelocityXField();
3360         else if(nameField=="velocityY" || nameField=="VelocityY" || nameField=="VELOCITYY" || nameField=="VitesseY" || nameField=="VITESSEY" || nameField=="vitesseY" )
3361                 return getVelocityYField();
3362         else if(nameField=="velocityZ" || nameField=="VelocityZ" || nameField=="VELOCITYZ" || nameField=="VitesseZ" || nameField=="VITESSEZ" || nameField=="vitesseZ" )
3363                 return getVelocityZField();
3364         else if(nameField=="temperature" || nameField=="Temperature" || nameField=="TEMPERATURE" || nameField=="temperature" )
3365                 return getTemperatureField();
3366         else if(nameField=="density" || nameField=="Density" || nameField=="DENSITY" || nameField=="Densite" || nameField=="DENSITE" || nameField=="densite" )
3367                 return getDensityField();
3368         else if(nameField=="momentum" || nameField=="Momentum" || nameField=="MOMENTUM" || nameField=="Qdm" || nameField=="QDM" || nameField=="qdm" )
3369                 return getMomentumField();
3370         else if(nameField=="enthalpy" || nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" || nameField=="enthalpie" )
3371                 return getEnthalpyField();
3372         else if(nameField=="totalenergy" || nameField=="TotalEnergy" || nameField=="TOTALENERGY" || nameField=="ENERGIETOTALE" || nameField=="EnergieTotale" || nameField=="energietotale" )
3373                 return getTotalEnergyField();
3374     else
3375     {
3376         cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first to check" << endl;
3377         throw CdmathException("SinglePhase::getOutputField error : Unknown Field name");
3378     }
3379 }