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