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