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