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