CoreFlows/Models/src/SinglePhase.cxx

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