CoreFlows/Models/src/SinglePhaseStaggered.cxx

   1 /*
   2  * SinglePhaseStaggered.cxx
   3  */
   4
   5 #include "SinglePhaseStaggered.hxx"
   6
   7 using namespace std;
   8
   9 void
  10 computeVelocityMCells(const Field& velocity,
  11                       Field& velocityMCells)
  12 {
  13     Mesh myMesh=velocity.getMesh();
  14     int nbCells=myMesh.getNumberOfCells();
  15
  16     for(int i=0;i<nbCells;i++)
  17     {
  18         std::vector< int > facesId=myMesh.getCell(i).getFacesId();
  19         velocityMCells(i)=(velocity(facesId[0])+velocity(facesId[1]))/2.;
  20     }
  21 }
  22
  23 SinglePhaseStaggered::SinglePhaseStaggered(phaseType fluid, pressureEstimate pEstimate, int dim, bool useDellacherieEOS){
  24         _Ndim=dim;
  25         _nVar=_Ndim+2;
  26         _nbPhases  = 1;
  27         _dragCoeffs=vector<double>(1,0);
  28         _fluides.resize(1);
  29         _useDellacherieEOS=useDellacherieEOS;
  30
  31         if(pEstimate==around1bar300K){//EOS at 1 bar and 300K
  32                 if(fluid==Gas){
  33                         cout<<"Fluid is air around 1 bar and 300 K (27°C)"<<endl;
  34                         *_runLogFile<<"Fluid is air around 1 bar and 300 K (27°C)"<<endl;
  35                         _fluides[0] = new StiffenedGas(1.4,743,300,2.22e5);  //ideal gas law for nitrogen at pressure 1 bar and temperature 27°C, e=2.22e5, c_v=743
  36                 }
  37                 else{
  38                         cout<<"Fluid is water around 1 bar and 300 K (27°C)"<<endl;
  39                         *_runLogFile<<"Fluid is water around 1 bar and 300 K (27°C)"<<endl;
  40                         _fluides[0] = new StiffenedGas(996,1e5,300,1.12e5,1501,4130);  //stiffened gas law for water at pressure 1 bar and temperature 27°C, e=1.12e5, c_v=4130
  41                 }
  42         }
  43         else{//EOS at 155 bars and 618K
  44                 if(fluid==Gas){
  45                         cout<<"Fluid is Gas around saturation point 155 bars and 618 K (345°C)"<<endl;
  46                         *_runLogFile<<"Fluid is Gas around saturation point 155 bars and 618 K (345°C)"<<endl;
  47                         _fluides[0] = new StiffenedGas(102,1.55e7,618,2.44e6, 433,3633);  //stiffened gas law for Gas at pressure 155 bar and temperature 345°C
  48                 }
  49                 else{//To do : change to normal regime: 155 bars and 573K
  50                         cout<<"Fluid is water around saturation point 155 bars and 618 K (345°C)"<<endl;
  51                         *_runLogFile<<"Fluid is water around saturation point 155 bars and 618 K (345°C)"<<endl;
  52                         if(_useDellacherieEOS)
  53                                 _fluides[0]= new StiffenedGasDellacherie(2.35,1e9,-1.167e6,1816); //stiffened gas law for water from S. Dellacherie
  54                         else
  55                                 _fluides[0]= new StiffenedGas(594.,1.55e7,618.,1.6e6, 621.,3100.);  //stiffened gas law for water at pressure 155 bar, and temperature 345°C
  56                 }
  57         }
  58 }
  59 void SinglePhaseStaggered::initialize(){
  60         cout<<"Initialising the Navier-Stokes model"<<endl;
  61         *_runLogFile<<"Initialising the Navier-Stokes model"<<endl;
  62
  63         _Uroe = new double[_nVar];
  64         _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
  65         for(int i=0; i<_Ndim; i++)
  66                 _gravite[i+1]=_GravityField3d[i];
  67
  68         _GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
  69         if(_saveVelocity)
  70                 _Vitesse=Field("Velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
  71
  72         if(_entropicCorrection)
  73                 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
  74
  75         ProblemFluid::initialize();
  76 }
  77
  78 bool SinglePhaseStaggered::iterateTimeStep(bool &converged)
  79 {
  80         if(_timeScheme == Explicit || !_usePrimitiveVarsInNewton)
  81                 ProblemFluid::iterateTimeStep(converged);
  82         else
  83         {
  84                 bool stop=false;
  85
  86                 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
  87                         _maxvp=0;
  88                         computeTimeStep(stop);
  89                 }
  90                 if(stop){//Le compute time step ne s'est pas bien passé
  91                         cout<<"ComputeTimeStep failed"<<endl;
  92                         converged=false;
  93                         return false;
  94                 }
  95
  96                 computeNewtonVariation();
  97
  98                 //converged=convergence des iterations
  99                 if(_timeScheme == Explicit)
 100                         converged=true;
 101                 else{//Implicit scheme
 102
 103                         KSPGetIterationNumber(_ksp, &_PetscIts);
 104                         if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
 105                                 _MaxIterLinearSolver = _PetscIts;
 106                         if(_PetscIts>=_maxPetscIts)//solving the linear system failed
 107                         {
 108                                 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 109                                 *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
 110                                 converged=false;
 111                                 return false;
 112                         }
 113                         else{//solving the linear system succeeded
 114                                 //Calcul de la variation relative Uk+1-Uk
 115                                 _erreur_rel = 0.;
 116                                 double x, dx;
 117                                 int I;
 118                                 for(int j=1; j<=_Nmailles; j++)
 119                                 {
 120                                         for(int k=0; k<_nVar; k++)
 121                                         {
 122                                                 I = (j-1)*_nVar + k;
 123                                                 VecGetValues(_newtonVariation, 1, &I, &dx);
 124                                                 VecGetValues(_primitiveVars, 1, &I, &x);
 125                                                 if (fabs(x)*fabs(x)< _precision)
 126                                                 {
 127                                                         if(_erreur_rel < fabs(dx))
 128                                                                 _erreur_rel = fabs(dx);
 129                                                 }
 130                                                 else if(_erreur_rel < fabs(dx/x))
 131                                                         _erreur_rel = fabs(dx/x);
 132                                         }
 133                                 }
 134                         }
 135                         converged = _erreur_rel <= _precision_Newton;
 136                 }
 137
 138                 double relaxation=1;//Vk+1=Vk+relaxation*deltaV
 139
 140                 VecAXPY(_primitiveVars,  relaxation, _newtonVariation);
 141
 142                 //mise a jour du champ primitif
 143                 updateConservatives();
 144
 145                 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
 146                 {
 147                         cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 148                         *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
 149                         converged=false;
 150                         return false;
 151                 }
 152                 if(_system)
 153                 {
 154                         cout<<"Vecteur Vkp1-Vk "<<endl;
 155                         VecView(_newtonVariation,  PETSC_VIEWER_STDOUT_SELF);
 156                         cout << "Nouvel etat courant Vk de l'iteration Newton: " << endl;
 157                         VecView(_primitiveVars,  PETSC_VIEWER_STDOUT_SELF);
 158                 }
 159
 160                 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
 161                         if(_minm1<-_precision || _minm2<-_precision)
 162                         {
 163                                 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 164                                 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
 165                         }
 166
 167                         if (_nbVpCplx>0){
 168                                 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 169                                 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
 170                         }
 171                 }
 172                 _minm1=1e30;
 173                 _minm2=1e30;
 174                 _nbMaillesNeg=0;
 175                 _nbVpCplx =0;
 176                 _part_imag_max=0;
 177
 178                 return true;
 179         }
 180 }
 181 void SinglePhaseStaggered::computeNewtonVariation()
 182 {
 183         if(!_usePrimitiveVarsInNewton)
 184                 ProblemFluid::computeNewtonVariation();
 185         else
 186         {
 187                 if(_verbose)
 188                 {
 189                         cout<<"Vecteur courant Vk "<<endl;
 190                         VecView(_primitiveVars,PETSC_VIEWER_STDOUT_SELF);
 191                         cout << endl;
 192                         cout << "Matrice du système linéaire avant contribution delta t" << endl;
 193                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 194                         cout << endl;
 195                         cout << "Second membre du système linéaire avant contribution delta t" << endl;
 196                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 197                         cout << endl;
 198                 }
 199                 if(_timeScheme == Explicit)
 200                 {
 201                         VecCopy(_b,_newtonVariation);
 202                         VecScale(_newtonVariation, _dt);
 203                         if(_verbose && _nbTimeStep%_freqSave ==0)
 204                         {
 205                                 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
 206                                 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
 207                                 cout << endl;
 208                         }
 209                 }
 210                 else
 211                 {
 212                         MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
 213
 214                         VecAXPY(_b, 1/_dt, _old);
 215                         VecAXPY(_b, -1/_dt, _conservativeVars);
 216
 217                         for(int imaille = 0; imaille<_Nmailles; imaille++){
 218                                 _idm[0] = _nVar*imaille;
 219                                 for(int k=1; k<_nVar; k++)
 220                                         _idm[k] = _idm[k-1] + 1;
 221                                 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
 222                                 primToConsJacobianMatrix(_Vi);
 223                                 for(int k=0; k<_nVar*_nVar; k++)
 224                                         _primToConsJacoMat[k]*=1/_dt;
 225                                 MatSetValuesBlocked(_A, 1, &imaille, 1, &imaille, _primToConsJacoMat, ADD_VALUES);
 226                         }
 227                         MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
 228
 229 #if PETSC_VERSION_GREATER_3_5
 230                         KSPSetOperators(_ksp, _A, _A);
 231 #else
 232                         KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
 233 #endif
 234
 235                         if(_verbose)
 236                         {
 237                                 cout << "Matrice du système linéaire" << endl;
 238                                 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 239                                 cout << endl;
 240                                 cout << "Second membre du système linéaire" << endl;
 241                                 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 242                                 cout << endl;
 243                         }
 244
 245                         if(_conditionNumber)
 246                                 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
 247                         if(!_isScaling)
 248                         {
 249                                 KSPSolve(_ksp, _b, _newtonVariation);
 250                         }
 251                         else
 252                         {
 253                                 computeScaling(_maxvp);
 254                                 int indice;
 255                                 VecAssemblyBegin(_vecScaling);
 256                                 VecAssemblyBegin(_invVecScaling);
 257                                 for(int imaille = 0; imaille<_Nmailles; imaille++)
 258                                 {
 259                                         indice = imaille;
 260                                         VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
 261                                         VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
 262                                 }
 263                                 VecAssemblyEnd(_vecScaling);
 264                                 VecAssemblyEnd(_invVecScaling);
 265
 266                                 if(_system)
 267                                 {
 268                                         cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
 269                                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 270                                         cout << endl;
 271                                         cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
 272                                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 273                                         cout << endl;
 274                                 }
 275                                 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
 276                                 if(_system)
 277                                 {
 278                                         cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
 279                                         MatView(_A,PETSC_VIEWER_STDOUT_SELF);
 280                                         cout << endl;
 281                                 }
 282                                 VecCopy(_b,_bScaling);
 283                                 VecPointwiseMult(_b,_vecScaling,_bScaling);
 284                                 if(_system)
 285                                 {
 286                                         cout << "Produit du second membre par le preconditionneur bloc diagonal  a gauche" << endl;
 287                                         VecView(_b, PETSC_VIEWER_STDOUT_SELF);
 288                                         cout << endl;
 289                                 }
 290
 291                                 KSPSolve(_ksp,_b, _bScaling);
 292                                 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
 293                         }
 294                         if(_system)
 295                         {
 296                                 cout << "solution du systeme lineaire local:" << endl;
 297                                 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
 298                                 cout << endl;
 299                         }
 300                 }
 301         }
 302 }
 303 void SinglePhaseStaggered::convectionState( const long &i, const long &j, const bool &IsBord){
 304
 305         _idm[0] = _nVar*i;
 306         for(int k=1; k<_nVar; k++)
 307                 _idm[k] = _idm[k-1] + 1;
 308         VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
 309
 310         _idm[0] = _nVar*j;
 311         for(int k=1; k<_nVar; k++)
 312                 _idm[k] = _idm[k-1] + 1;
 313         if(IsBord)
 314                 VecGetValues(_Uext, _nVar, _idm, _Uj);
 315         else
 316                 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 317         if(_verbose && _nbTimeStep%_freqSave ==0)
 318         {
 319                 cout<<"Convection Left state cell " << i<< ": "<<endl;
 320                 for(int k =0; k<_nVar; k++)
 321                         cout<< _Ui[k]<<endl;
 322                 cout<<"Convection Right state cell " << j<< ": "<<endl;
 323                 for(int k =0; k<_nVar; k++)
 324                         cout<< _Uj[k]<<endl;
 325         }
 326         if(_Ui[0]<0||_Uj[0]<0)
 327         {
 328                 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 329                 *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
 330                 _runLogFile->close();
 331                 throw CdmathException("densite negative, arret de calcul");
 332         }
 333         PetscScalar ri, rj, xi, xj, pi, pj;
 334         PetscInt Ii;
 335         ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
 336         rj = sqrt(_Uj[0]);
 337         _Uroe[0] = ri*rj;       //moyenne geometrique des densites
 338         if(_verbose && _nbTimeStep%_freqSave ==0)
 339                 cout << "Densité moyenne Roe  gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
 340         for(int k=0;k<_Ndim;k++){
 341                 xi = _Ui[k+1];
 342                 xj = _Uj[k+1];
 343                 _Uroe[1+k] = (xi/ri + xj/rj)/(ri + rj);
 344                 //"moyenne" des vitesses
 345                 if(_verbose && _nbTimeStep%_freqSave ==0)
 346                         cout << "Vitesse de Roe composante "<< k<<"  gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+1] << endl;
 347         }
 348         // H = (rho E + p)/rho
 349         xi = _Ui[_nVar-1];//phi rho E
 350         xj = _Uj[_nVar-1];
 351         Ii = i*_nVar; // correct Kieu
 352         VecGetValues(_primitiveVars, 1, &Ii, &pi);// _primitiveVars pour p
 353         if(IsBord)
 354         {
 355                 double q_2 = 0;
 356                 for(int k=1;k<=_Ndim;k++)
 357                         q_2 += _Uj[k]*_Uj[k];
 358                 q_2 /= _Uj[0];  //phi rho u²
 359                 pj =  _fluides[0]->getPressure((_Uj[(_Ndim+2)-1]-q_2/2)/_porosityj,_Uj[0]/_porosityj);
 360         }
 361         else
 362         {
 363                 Ii = j*_nVar; // correct Kieu
 364                 VecGetValues(_primitiveVars, 1, &Ii, &pj);
 365         }
 366         xi = (xi + pi)/(ri*ri);
 367         xj = (xj + pj)/(rj*rj);
 368         _Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
 369         //on se donne l enthalpie ici
 370         if(_verbose && _nbTimeStep%_freqSave ==0)
 371                 cout << "Enthalpie totale de Roe H  gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
 372
 373         if(_verbose && _nbTimeStep%_freqSave ==0)
 374         {
 375                 cout<<"Convection interfacial state"<<endl;
 376                 for(int k=0;k<_nVar;k++)
 377                         cout<< _Uroe[k]<<" , "<<endl;
 378         }
 379 }
 380
 381 void SinglePhaseStaggered::setBoundaryState(string nameOfGroup, const int &j,double *normale){
 382         int k;
 383         double v2=0, q_n=0;//q_n=quantité de mouvement normale à la face frontière;
 384         _idm[0] = _nVar*j;
 385         for(k=1; k<_nVar; k++)
 386                 _idm[k] = _idm[k-1] + 1;
 387
 388         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
 389         for(k=0; k<_Ndim; k++)
 390                 q_n+=_externalStates[(k+1)]*normale[k];
 391
 392         double porosityj=_porosityField(j);
 393
 394         if(_verbose && _nbTimeStep%_freqSave ==0)
 395         {
 396                 cout << "setBoundaryState for group "<< nameOfGroup << ", inner cell j= "<<j<< " face unit normal vector "<<endl;
 397                 for(k=0; k<_Ndim; k++){
 398                         cout<<normale[k]<<", ";
 399                 }
 400                 cout<<endl;
 401         }
 402
 403         if (_limitField[nameOfGroup].bcType==Wall){
 404                 //Pour la convection, inversion du sens de la vitesse normale
 405                 for(k=0; k<_Ndim; k++)
 406                         _externalStates[(k+1)]-= 2*q_n*normale[k];
 407
 408                 _idm[0] = 0;
 409                 for(k=1; k<_nVar; k++)
 410                         _idm[k] = _idm[k-1] + 1;
 411
 412                 VecAssemblyBegin(_Uext);
 413                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 414                 VecAssemblyEnd(_Uext);
 415
 416                 //Pour la diffusion, paroi à vitesse et temperature imposees
 417                 _idm[0] = _nVar*j;
 418                 for(k=1; k<_nVar; k++)
 419                         _idm[k] = _idm[k-1] + 1;
 420                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
 421                 double pression=_externalStates[0];
 422                 double T=_limitField[nameOfGroup].T;
 423                 double rho=_fluides[0]->getDensity(pression,T);
 424
 425                 _externalStates[0]=porosityj*rho;
 426                 _externalStates[1]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
 427                 v2 +=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
 428                 if(_Ndim>1)
 429                 {
 430                         v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 431                         _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
 432                         if(_Ndim==3)
 433                         {
 434                                 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
 435                                 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 436                         }
 437                 }
 438                 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
 439                 _idm[0] = 0;
 440                 for(k=1; k<_nVar; k++)
 441                         _idm[k] = _idm[k-1] + 1;
 442                 VecAssemblyBegin(_Uextdiff);
 443                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 444                 VecAssemblyEnd(_Uextdiff);
 445         }
 446         else if (_limitField[nameOfGroup].bcType==Neumann){
 447                 _idm[0] = 0;
 448                 for(k=1; k<_nVar; k++)
 449                         _idm[k] = _idm[k-1] + 1;
 450
 451                 VecAssemblyBegin(_Uext);
 452                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 453                 VecAssemblyEnd(_Uext);
 454
 455                 VecAssemblyBegin(_Uextdiff);
 456                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 457                 VecAssemblyEnd(_Uextdiff);
 458         }
 459         else if (_limitField[nameOfGroup].bcType==Inlet){
 460
 461                 if(q_n<=0){
 462                         VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
 463                         double pression=_externalStates[0];
 464                         double T=_limitField[nameOfGroup].T;
 465                         double rho=_fluides[0]->getDensity(pression,T);
 466
 467                         _externalStates[0]=porosityj*rho;
 468                         _externalStates[1]=_externalStates[0]*(_limitField[nameOfGroup].v_x[0]);
 469                         v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
 470                         if(_Ndim>1)
 471                         {
 472                                 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 473                                 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
 474                                 if(_Ndim==3)
 475                                 {
 476                                         _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
 477                                         v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 478                                 }
 479                         }
 480                         _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
 481                 }
 482                 else if(_nbTimeStep%_freqSave ==0)
 483                         cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
 484
 485                 _idm[0] = 0;
 486                 for(k=1; k<_nVar; k++)
 487                         _idm[k] = _idm[k-1] + 1;
 488                 VecAssemblyBegin(_Uext);
 489                 VecAssemblyBegin(_Uextdiff);
 490                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 491                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 492                 VecAssemblyEnd(_Uext);
 493                 VecAssemblyEnd(_Uextdiff);
 494         }
 495         else if (_limitField[nameOfGroup].bcType==InletPressure){
 496
 497                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 498                 Cell Cj=_mesh.getCell(j);
 499                 double hydroPress=Cj.x()*_GravityField3d[0];
 500                 if(_Ndim>1){
 501                         hydroPress+=Cj.y()*_GravityField3d[1];
 502                         if(_Ndim>2)
 503                                 hydroPress+=Cj.z()*_GravityField3d[2];
 504                 }
 505                 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
 506
 507                 //Building the external state
 508                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
 509                 if(q_n<=0){
 510                         _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress,_limitField[nameOfGroup].T);
 511                 }
 512                 else{
 513                         if(_nbTimeStep%_freqSave ==0)
 514                                 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Neumann boundary condition for velocity and temperature"<<endl;
 515                         _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
 516                 }
 517
 518                 for(k=0; k<_Ndim; k++)
 519                 {
 520                         v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
 521                         _externalStates[(k+1)]*=_externalStates[0] ;
 522                 }
 523                 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);
 524
 525
 526                 _idm[0] = 0;
 527                 for(k=1; k<_nVar; k++)
 528                         _idm[k] = _idm[k-1] + 1;
 529                 VecAssemblyBegin(_Uext);
 530                 VecAssemblyBegin(_Uextdiff);
 531                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 532                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 533                 VecAssemblyEnd(_Uext);
 534                 VecAssemblyEnd(_Uextdiff);
 535         }
 536         else if (_limitField[nameOfGroup].bcType==Outlet){
 537                 if(q_n<=0 &&  _nbTimeStep%_freqSave ==0)
 538                         cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition for velocity and temperature"<<endl;
 539
 540                 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
 541                 Cell Cj=_mesh.getCell(j);
 542                 double hydroPress=Cj.x()*_GravityField3d[0];
 543                 if(_Ndim>1){
 544                         hydroPress+=Cj.y()*_GravityField3d[1];
 545                         if(_Ndim>2)
 546                                 hydroPress+=Cj.z()*_GravityField3d[2];
 547                 }
 548                 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
 549
 550                 if(_verbose && _nbTimeStep%_freqSave ==0)
 551                 {
 552                         cout<<"Cond lim outlet densite= "<<_externalStates[0]<<" gravite= "<<_GravityField3d[0]<<" Cj.x()= "<<Cj.x()<<endl;
 553                         cout<<"Cond lim outlet pression ref= "<<_limitField[nameOfGroup].p<<" pression hydro= "<<hydroPress<<" total= "<<_limitField[nameOfGroup].p+hydroPress<<endl;
 554                 }
 555                 //Building the external state
 556                 _idm[0] = _nVar*j;// Kieu
 557                 for(k=1; k<_nVar; k++)
 558                         _idm[k] = _idm[k-1] + 1;
 559                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
 560
 561                 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
 562                 for(k=0; k<_Ndim; k++)
 563                 {
 564                         v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
 565                         _externalStates[(k+1)]*=_externalStates[0] ;
 566                 }
 567                 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);
 568                 _idm[0] = 0;
 569                 for(k=1; k<_nVar; k++)
 570                         _idm[k] = _idm[k-1] + 1;
 571                 VecAssemblyBegin(_Uext);
 572                 VecAssemblyBegin(_Uextdiff);
 573                 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
 574                 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
 575                 VecAssemblyEnd(_Uext);
 576                 VecAssemblyEnd(_Uextdiff);
 577         }else {
 578                 cout<<"Boundary condition not set for boundary named "<<nameOfGroup<<endl;
 579                 cout<<"Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
 580                 *_runLogFile<<"Boundary condition not set for boundary named. Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
 581                 _runLogFile->close();
 582                 throw CdmathException("Unknown boundary condition");
 583         }
 584 }
 585
 586 void SinglePhaseStaggered::convectionMatrices()
 587 {
 588         //entree: URoe = rho, u, H
 589         //sortie: matrices Roe+  et Roe-
 590
 591         if(_verbose && _nbTimeStep%_freqSave ==0)
 592                 cout<<"SinglePhaseStaggered::convectionMatrices()"<<endl;
 593
 594         double u_n=0, u_2=0;//vitesse normale et carré du module
 595
 596         for(int i=0;i<_Ndim;i++)
 597         {
 598                 u_2 += _Uroe[1+i]*_Uroe[1+i];
 599                 u_n += _Uroe[1+i]*_vec_normal[i];
 600         }
 601
 602         vector<complex<double> > vp_dist(3,0);
 603
 604         if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
 605         {
 606                 staggeredVFFCMatricesConservativeVariables(u_n);//Computation of classical upwinding matrices
 607                 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//For use in implicit matrix
 608                         staggeredVFFCMatricesPrimitiveVariables(u_n);
 609         }
 610         else
 611         {
 612                 Vector vitesse(_Ndim);
 613                 for(int idim=0;idim<_Ndim;idim++)
 614                         vitesse[idim]=_Uroe[1+idim];
 615
 616                 double  c, H, K, k;
 617                 /***********Calcul des valeurs propres ********/
 618                 H = _Uroe[_nVar-1];
 619                 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
 620                 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
 621                 K = u_2*k/2; //g-1/2 *|u|²
 622
 623                 vp_dist[0]=u_n-c;vp_dist[1]=u_n;vp_dist[2]=u_n+c;
 624
 625                 _maxvploc=fabs(u_n)+c;
 626                 if(_maxvploc>_maxvp)
 627                         _maxvp=_maxvploc;
 628
 629                 if(_verbose && _nbTimeStep%_freqSave ==0)
 630                         cout<<"SinglePhaseStaggered::convectionMatrices Eigenvalues "<<u_n-c<<" , "<<u_n<<" , "<<u_n+c<<endl;
 631
 632                 RoeMatrixConservativeVariables( u_n, H,vitesse,k,K);
 633
 634                 /******** Construction des matrices de decentrement ********/
 635                 if( _spaceScheme ==centered){
 636                         if(_entropicCorrection)
 637                         {
 638                                 *_runLogFile<<"SinglePhaseStaggered::convectionMatrices: entropy scheme not available for centered scheme"<<endl;
 639                                 _runLogFile->close();
 640                                 throw CdmathException("SinglePhaseStaggered::convectionMatrices: entropy scheme not available for centered scheme");
 641                         }
 642
 643                         for(int i=0; i<_nVar*_nVar;i++)
 644                                 _absAroe[i] = 0;
 645                 }
 646                 else if(_spaceScheme == upwind || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
 647                         if(_entropicCorrection)
 648                                 entropicShift(_vec_normal);
 649                         else
 650                                 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
 651
 652                         vector< complex< double > > y (3,0);
 653                         for( int i=0 ; i<3 ; i++)
 654                                 y[i] = Polynoms::abs_generalise(vp_dist[i])+_entropicShift[i];
 655                         Polynoms::abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
 656
 657                         if( _spaceScheme ==pressureCorrection)
 658                                 for( int i=0 ; i<_Ndim ; i++)
 659                                         for( int j=0 ; j<_Ndim ; j++)
 660                                                 _absAroe[(1+i)*_nVar+1+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 661                         else if( _spaceScheme ==lowMach){
 662                                 double M=sqrt(u_2)/c;
 663                                 for( int i=0 ; i<_Ndim ; i++)
 664                                         for( int j=0 ; j<_Ndim ; j++)
 665                                                 _absAroe[(1+i)*_nVar+1+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 666                         }
 667                 }
 668                 else if( _spaceScheme ==staggered ){
 669                         if(_entropicCorrection)//To do: study entropic correction for staggered
 670                         {
 671                                 *_runLogFile<<"SinglePhaseStaggered::convectionMatrices: entropy scheme not available for staggered scheme"<<endl;
 672                                 _runLogFile->close();
 673                                 throw CdmathException("SinglePhaseStaggered::convectionMatrices: entropy scheme not available for staggered scheme");
 674                         }
 675
 676                         staggeredRoeUpwindingMatrixConservativeVariables( u_n, H, vitesse, k, K);
 677                 }
 678                 else
 679                 {
 680                         *_runLogFile<<"SinglePhaseStaggered::convectionMatrices: scheme not treated"<<endl;
 681                         _runLogFile->close();
 682                         throw CdmathException("SinglePhaseStaggered::convectionMatrices: scheme not treated");
 683                 }
 684
 685                 for(int i=0; i<_nVar*_nVar;i++)
 686                 {
 687                         _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
 688                         _AroePlus[i]  = (_Aroe[i]+_absAroe[i])/2;
 689                 }
 690                 if(_timeScheme==Implicit)
 691                 {
 692                         if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
 693                         {
 694                                 _Vij[0]=_fluides[0]->getPressureFromEnthalpy(_Uroe[_nVar-1]-u_2/2, _Uroe[0]);//pressure
 695                                 _Vij[_nVar-1]=_fluides[0]->getTemperatureFromPressure( _Vij[0], _Uroe[0]);//Temperature
 696                                 for(int idim=0;idim<_Ndim; idim++)
 697                                         _Vij[1+idim]=_Uroe[1+idim];
 698                                 primToConsJacobianMatrix(_Vij);
 699                                 Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
 700                                 Polynoms::matrixProduct(_AroePlus,  _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
 701                         }
 702                         else
 703                                 for(int i=0; i<_nVar*_nVar;i++)
 704                                 {
 705                                         _AroeMinusImplicit[i] = _AroeMinus[i];
 706                                         _AroePlusImplicit[i]  = _AroePlus[i];
 707                                 }
 708                 }
 709                 if(_verbose && _nbTimeStep%_freqSave ==0)
 710                 {
 711                         displayMatrix(_Aroe, _nVar,"Matrice de Roe");
 712                         displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
 713                         displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
 714                         displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
 715                 }
 716         }
 717
 718         if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
 719         {
 720                 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
 721                 displayMatrix(_AroePlusImplicit,  _nVar,"Matrice _AroePlusImplicit");
 722         }
 723
 724         /*********Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source*****/
 725         if(_entropicCorrection)
 726         {
 727                 InvMatriceRoe( vp_dist);
 728                 Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
 729         }
 730         else if (_spaceScheme==upwind || (_spaceScheme ==pressureCorrection ) || (_spaceScheme ==lowMach ))//upwind sans entropic
 731                 SigneMatriceRoe( vp_dist);
 732         else if (_spaceScheme==centered)//centre  sans entropic
 733                 for(int i=0; i<_nVar*_nVar;i++)
 734                         _signAroe[i] = 0;
 735         else if( _spaceScheme ==staggered )//à tester
 736         {
 737                 double signu;
 738                 if(u_n>0)
 739                         signu=1;
 740                 else if (u_n<0)
 741                         signu=-1;
 742                 else
 743                         signu=0;
 744                 for(int i=0; i<_nVar*_nVar;i++)
 745                         _signAroe[i] = 0;
 746                 _signAroe[0] = signu;
 747                 for(int i=1; i<_nVar-1;i++)
 748                         _signAroe[i*_nVar+i] = -signu;
 749                 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
 750         }
 751         else
 752         {
 753                 *_runLogFile<<"SinglePhaseStaggered::convectionMatrices: well balanced option not treated"<<endl;
 754                 _runLogFile->close();
 755                 throw CdmathException("SinglePhaseStaggered::convectionMatrices: well balanced option not treated");
 756         }
 757 }
 758 void SinglePhaseStaggered::computeScaling(double maxvp)
 759 {
 760         _blockDiag[0]=1;
 761         _invBlockDiag[0]=1;
 762         for(int q=1; q<_nVar-1; q++)
 763         {
 764                 _blockDiag[q]=1./maxvp;//
 765                 _invBlockDiag[q]= maxvp;//1.;//
 766         }
 767         _blockDiag[_nVar - 1]=(_fluides[0]->constante("gamma")-1)/(maxvp*maxvp);//1
 768         _invBlockDiag[_nVar - 1]=  1./_blockDiag[_nVar - 1] ;// 1.;//
 769 }
 770
 771
 772 void SinglePhaseStaggered::sourceVector(PetscScalar * Si, PetscScalar * Ui, PetscScalar * Vi, int i)
 773 {
 774         double phirho=Ui[0], T=Vi[_nVar-1];
 775         double norm_u=0;
 776         for(int k=0; k<_Ndim; k++)
 777                 norm_u+=Vi[1+k]*Vi[1+k];
 778         norm_u=sqrt(norm_u);
 779         if(T>_Tsat)
 780                 Si[0]=_heatPowerField(i)/_latentHeat;
 781         else
 782                 Si[0]=0;
 783         for(int k=1; k<_nVar-1; k++)
 784                 Si[k]  =(_gravite[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*phirho;
 785
 786         Si[_nVar-1]=_heatPowerField(i);
 787
 788         for(int k=0; k<_Ndim; k++)
 789                 Si[_nVar-1] +=(_GravityField3d[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*Vi[1+k]*phirho;
 790
 791         if(_timeScheme==Implicit)
 792         {
 793                 for(int k=0; k<_nVar*_nVar;k++)
 794                         _GravityImplicitationMatrix[k] = 0;
 795                 if(!_usePrimitiveVarsInNewton)
 796                         for(int k=0; k<_nVar;k++)
 797                                 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
 798                 else
 799                 {
 800                         double pression=Vi[0];
 801                         getDensityDerivatives( pression, T, norm_u*norm_u);
 802                         for(int k=0; k<_nVar;k++)
 803                         {
 804                                 _GravityImplicitationMatrix[k*_nVar+0]      =-_gravite[k]*_drho_sur_dp;
 805                                 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
 806                         }
 807                 }
 808         }
 809         if(_verbose && _nbTimeStep%_freqSave ==0)
 810         {
 811                 cout<<"SinglePhaseStaggered::sourceVector"<<endl;
 812                 cout<<"Ui="<<endl;
 813                 for(int k=0;k<_nVar;k++)
 814                         cout<<Ui[k]<<", ";
 815                 cout<<endl;
 816                 cout<<"Vi="<<endl;
 817                 for(int k=0;k<_nVar;k++)
 818                         cout<<Vi[k]<<", ";
 819                 cout<<endl;
 820                 cout<<"Si="<<endl;
 821                 for(int k=0;k<_nVar;k++)
 822                         cout<<Si[k]<<", ";
 823                 cout<<endl;
 824                 if(_timeScheme==Implicit)
 825                         displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
 826         }
 827 }
 828
 829 void SinglePhaseStaggered::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
 830 {
 831         double norm_u=0, u_n=0, rho;
 832         for(int i=0;i<_Ndim;i++)
 833                 u_n += _Uroe[1+i]*_vec_normal[i];
 834
 835         pressureLoss[0]=0;
 836         if(u_n>0){
 837                 for(int i=0;i<_Ndim;i++)
 838                         norm_u += Vi[1+i]*Vi[1+i];
 839                 norm_u=sqrt(norm_u);
 840                 rho=Ui[0];
 841                 for(int i=0;i<_Ndim;i++)
 842                         pressureLoss[1+i]=-1/2*K*rho*norm_u*Vi[1+i];
 843         }
 844         else{
 845                 for(int i=0;i<_Ndim;i++)
 846                         norm_u += Vj[1+i]*Vj[1+i];
 847                 norm_u=sqrt(norm_u);
 848                 rho=Uj[0];
 849                 for(int i=0;i<_Ndim;i++)
 850                         pressureLoss[1+i]=-1/2*K*rho*norm_u*Vj[1+i];
 851         }
 852         pressureLoss[_nVar-1]=-1/2*K*rho*norm_u*norm_u*norm_u;
 853
 854         if(_verbose && _nbTimeStep%_freqSave ==0)
 855         {
 856                 cout<<"SinglePhaseStaggered::pressureLossVector K= "<<K<<endl;
 857                 cout<<"Ui="<<endl;
 858                 for(int k=0;k<_nVar;k++)
 859                         cout<<Ui[k]<<", ";
 860                 cout<<endl;
 861                 cout<<"Vi="<<endl;
 862                 for(int k=0;k<_nVar;k++)
 863                         cout<<Vi[k]<<", ";
 864                 cout<<endl;
 865                 cout<<"Uj="<<endl;
 866                 for(int k=0;k<_nVar;k++)
 867                         cout<<Uj[k]<<", ";
 868                 cout<<endl;
 869                 cout<<"Vj="<<endl;
 870                 for(int k=0;k<_nVar;k++)
 871                         cout<<Vj[k]<<", ";
 872                 cout<<endl;
 873                 cout<<"pressureLoss="<<endl;
 874                 for(int k=0;k<_nVar;k++)
 875                         cout<<pressureLoss[k]<<", ";
 876                 cout<<endl;
 877         }
 878 }
 879
 880 void SinglePhaseStaggered::porosityGradientSourceVector()
 881 {
 882         double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[0], pj=_Vj[0],pij;
 883         for(int i=0;i<_Ndim;i++) {
 884                 u_ni += _Vi[1+i]*_vec_normal[i];
 885                 u_nj += _Vj[1+i]*_vec_normal[i];
 886         }
 887         _porosityGradientSourceVector[0]=0;
 888         rhoj=_Uj[0]/_porosityj;
 889         rhoi=_Ui[0]/_porosityi;
 890         pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
 891         for(int i=0;i<_Ndim;i++)
 892                 _porosityGradientSourceVector[1+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
 893         _porosityGradientSourceVector[_nVar-1]=0;
 894 }
 895
 896 void SinglePhaseStaggered::jacobian(const int &j, string nameOfGroup,double * normale)
 897 {
 898         if(_verbose && _nbTimeStep%_freqSave ==0)
 899                 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
 900
 901         int k;
 902         for(k=0; k<_nVar*_nVar;k++)
 903                 _Jcb[k] = 0;//No implicitation at this stage
 904
 905         _idm[0] = _nVar*j;
 906         for(k=1; k<_nVar; k++)
 907                 _idm[k] = _idm[k-1] + 1;
 908         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
 909         double q_n=0;//quantité de mouvement normale à la paroi
 910         for(k=0; k<_Ndim; k++)
 911                 q_n+=_externalStates[(k+1)]*normale[k];
 912
 913         // loop of boundary types
 914         if (_limitField[nameOfGroup].bcType==Wall)
 915         {
 916                 for(k=0; k<_nVar;k++)
 917                         _Jcb[k*_nVar + k] = 1;
 918                 for(k=1; k<_nVar-1;k++)
 919                         for(int l=1; l<_nVar-1;l++)
 920                                 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
 921         }
 922         else if (_limitField[nameOfGroup].bcType==Inlet)
 923         {
 924                 return;
 925                 if(q_n<0){
 926                         double v[_Ndim], ve[_Ndim], v2, ve2;
 927
 928                         _idm[0] = _nVar*j;
 929                         for(k=1; k<_nVar; k++)
 930                                 _idm[k] = _idm[k-1] + 1;
 931                         VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 932                         VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 933
 934                         ve[0] = _limitField[nameOfGroup].v_x[0];
 935                         v[0]=_Vj[1];
 936                         ve2 = ve[0]*ve[0];
 937                         v2 = v[0]*v[0];
 938                         if (_Ndim >1){
 939                                 ve[1] = _limitField[nameOfGroup].v_y[0];
 940                                 v[1]=_Vj[2];
 941                                 ve2 += ve[1]*ve[1];
 942                                 v2 = v[1]*v[1];
 943                         }
 944                         if (_Ndim >2){
 945                                 ve[2] = _limitField[nameOfGroup].v_z[0];
 946                                 v[2]=_Vj[3];
 947                                 ve2 += ve[2]*ve[2];
 948                                 v2 = v[2]*v[2];
 949                         }
 950                         double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,_Uj[0]);
 951                         double total_energy=internal_energy+ve2/2;
 952
 953                         //Mass line
 954                         _Jcb[0]=v2/(2*internal_energy);
 955                         for(k=0; k<_Ndim;k++)
 956                                 _Jcb[1+k]=-v[k]/internal_energy;
 957                         _Jcb[_nVar-1]=1/internal_energy;
 958                         //Momentum lines
 959                         for(int l =1;l<1+_Ndim;l++){
 960                                 _Jcb[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
 961                                 for(k=0; k<_Ndim;k++)
 962                                         _Jcb[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
 963                                 _Jcb[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
 964                         }
 965                         //Energy line
 966                         _Jcb[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
 967                         for(k=0; k<_Ndim;k++)
 968                                 _Jcb[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
 969                         _Jcb[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
 970                 }
 971                 else
 972                         for(k=0;k<_nVar;k++)
 973                                 _Jcb[k*_nVar+k]=1;
 974                 //Old jacobian
 975                 /*
 976                  _Jcb[0] = 1;
 977                 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];//Kieu
 978                 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
 979                 if(_Ndim>1)
 980                 {
 981                         _Jcb[2*_nVar]= _limitField[nameOfGroup].v_y[0];
 982                         v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 983                         if(_Ndim==3){
 984                                 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
 985                                 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 986                         }
 987                 }
 988                 _Jcb[(_nVar-1)*_nVar]=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2;
 989                  */
 990         }
 991         else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0){
 992                 return;
 993                 double v[_Ndim], v2=0;
 994                 _idm[0] = _nVar*j;
 995                 for(k=1; k<_nVar; k++)
 996                         _idm[k] = _idm[k-1] + 1;
 997                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 998
 999                 for(k=0; k<_Ndim;k++){
1000                         v[k]=_Vj[1+k];
1001                         v2+=v[k]*v[k];
1002                 }
1003
1004                 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _limitField[nameOfGroup].T);
1005                 double rho_int = _externalStates[0];
1006                 double density_ratio=rho_ext/rho_int;
1007                 //Momentum lines
1008                 for(int l =1;l<1+_Ndim;l++){
1009                         _Jcb[l*_nVar]=-density_ratio*v[l-1];
1010                         _Jcb[l*_nVar+l]=density_ratio;
1011                 }
1012                 //Energy lines
1013                 _Jcb[(_nVar-1)*_nVar]=-v2*density_ratio;
1014                 for(k=0; k<_Ndim;k++)
1015                         _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k];
1016         }
1017         // not wall, not inlet, not inletPressure
1018         else if(_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0))
1019         {
1020                 return;
1021                 double v[_Ndim], v2=0;
1022                 _idm[0] = _nVar*j;
1023                 for(k=1; k<_nVar; k++)
1024                         _idm[k] = _idm[k-1] + 1;
1025                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1026
1027                 for(k=0; k<_Ndim;k++){
1028                         v[k]=_Vj[1+k];
1029                         v2+=v[k]*v[k];
1030                 }
1031
1032                 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _externalStates[_nVar-1]);
1033                 double rho_int = _externalStates[0];
1034                 double density_ratio=rho_ext/rho_int;
1035                 double internal_energy=_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho_int);
1036                 double total_energy=internal_energy+v2/2;
1037
1038                 //Mass line
1039                 _Jcb[0]=density_ratio*(1-v2/(2*internal_energy));
1040                 for(k=0; k<_Ndim;k++)
1041                         _Jcb[1+k]=density_ratio*v[k]/internal_energy;
1042                 _Jcb[_nVar-1]=-density_ratio/internal_energy;
1043                 //Momentum lines
1044                 for(int l =1;l<1+_Ndim;l++){
1045                         _Jcb[l*_nVar]=density_ratio*v2*v[l-1]/(2*internal_energy);
1046                         for(k=0; k<_Ndim;k++)
1047                                 _Jcb[l*_nVar+1+k]=density_ratio*v[k]*v[l-1]/internal_energy;
1048                         _Jcb[l*_nVar+1+k]-=density_ratio;
1049                         _Jcb[l*_nVar+_nVar-1]=-density_ratio*v[l-1]/internal_energy;
1050                 }
1051                 //Energy line
1052                 _Jcb[(_nVar-1)*_nVar]=density_ratio*v2*total_energy/(2*internal_energy);
1053                 for(k=0; k<_Ndim;k++)
1054                         _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k]*total_energy/internal_energy;
1055                 _Jcb[(_nVar-1)*_nVar+_nVar-1]=density_ratio*(1-total_energy/internal_energy);
1056                 //Old jacobian
1057                 /*
1058                 int idim,jdim;
1059                 double cd = 1,cn=0,p0, gamma;
1060                 _idm[0] = j*_nVar;// Kieu
1061                 for(k=1; k<_nVar;k++)
1062                         _idm[k] = _idm[k-1] + 1;
1063                 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1064                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1065
1066                 // compute the common numerator and common denominator
1067                 p0=_fluides[0]->constante("p0");
1068                 gamma =_fluides[0]->constante("gamma");
1069                 cn =_limitField[nameOfGroup].p +p0;
1070                 cd = _phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0;
1071                 cd*=cd;
1072                 cd*=(gamma-1);
1073                 //compute the v2
1074                 for(k=1; k<_nVar-1;k++)
1075                         v2+=_externalStates[k]*_externalStates[k];
1076                 // drho_ext/dU
1077                 _JcbDiff[0] = cn*(_phi[_nVar-1] -v2 -p0)/cd;
1078                 for(k=1; k<_nVar-1;k++)
1079                         _JcbDiff[k]=cn*_phi[k]/cd;
1080                 _JcbDiff[_nVar-1]= -cn*_phi[0]/cd;
1081                 //dq_ext/dU
1082                 for(idim=0; idim<_Ndim;idim++)
1083                 {
1084                         //premiere colonne
1085                         _JcbDiff[(1+idim)*_nVar]=-(v2*cn*_phi[idim+1])/(2*cd);
1086                         //colonnes intermediaire
1087                         for(jdim=0; jdim<_Ndim;jdim++)
1088                         {
1089                                 _JcbDiff[(1+idim)*_nVar + jdim + 1] =_externalStates[idim+1]*_phi[jdim+1];
1090                                 _JcbDiff[(1+idim)*_nVar + jdim + 1]*=cn/cd;
1091                         }
1092                         //matrice identite*cn*(rhoe- p0)
1093                         _JcbDiff[(1+idim)*_nVar + idim + 1] +=( cn*(_phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0))/cd;
1094
1095                         //derniere colonne
1096                         _JcbDiff[(1+idim)*_nVar + _nVar-1]=-_phi[idim+1]*cn/cd;
1097                 }
1098                 //drhoE/dU
1099                 _JcbDiff[_nVar*(_nVar-1)] = -(v2*_phi[_nVar -1]*cn)/(2*cd);
1100                 for(int idim=0; idim<_Ndim;idim++)
1101                         _JcbDiff[_nVar*(_nVar-1)+idim+1]=_externalStates[idim+1]*_phi[_nVar -1]*cn/cd;
1102                 _JcbDiff[_nVar*_nVar -1] = -(v2/2+p0)*cn/cd;
1103                  */
1104         }
1105         else  if (_limitField[nameOfGroup].bcType!=Neumann)// not wall, not inlet, not outlet
1106         {
1107                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1108                 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1109                 _runLogFile->close();
1110                 throw CdmathException("SinglePhaseStaggered::jacobian: This boundary condition is not treated");
1111         }
1112 }
1113
1114
1115 Vector SinglePhaseStaggered::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
1116         if(_verbose && _nbTimeStep%_freqSave ==0)
1117         {
1118                 cout<<"SinglePhaseStaggered::convectionFlux start"<<endl;
1119                 cout<<"Ucons"<<endl;
1120                 cout<<U<<endl;
1121                 cout<<"Vprim"<<endl;
1122                 cout<<V<<endl;
1123         }
1124
1125         double phirho=U(0);//phi rho
1126         Vector phiq(_Ndim);//phi rho u
1127         for(int i=0;i<_Ndim;i++)
1128                 phiq(i)=U(1+i);
1129
1130         double pression=V(0);
1131         Vector vitesse(_Ndim);
1132         for(int i=0;i<_Ndim;i++)
1133                 vitesse(i)=V(1+i);
1134         double Temperature= V(1+_Ndim);
1135
1136         double vitessen=vitesse*normale;
1137         double rho=phirho/porosity;
1138         double e_int=_fluides[0]->getInternalEnergy(Temperature,rho);
1139
1140         Vector F(_nVar);
1141         F(0)=phirho*vitessen;
1142         for(int i=0;i<_Ndim;i++)
1143                 F(1+i)=phirho*vitessen*vitesse(i)+pression*normale(i)*porosity;
1144         F(1+_Ndim)=phirho*(e_int+0.5*vitesse*vitesse+pression/rho)*vitessen;
1145
1146         if(_verbose && _nbTimeStep%_freqSave ==0)
1147         {
1148                 cout<<"SinglePhaseStaggered::convectionFlux end"<<endl;
1149                 cout<<"Flux F(U,V)"<<endl;
1150                 cout<<F<<endl;
1151         }
1152
1153         return F;
1154 }
1155
1156 void SinglePhaseStaggered::convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector vitesse)
1157 {
1158         //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1159         //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
1160         //EOS is more involved with primitive variables
1161         // call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
1162         _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1163         for(int i=0;i<_Ndim;i++)
1164                 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1165         _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1166         for(int i=0;i<_Ndim;i++)
1167         {
1168                 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]+_vec_normal[i];
1169                 for(int j=0;j<_Ndim;j++)
1170                         _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1171                 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1172                 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1173         }
1174         _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1175         for(int i=0;i<_Ndim;i++)
1176                 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1177         _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1178 }
1179
1180 void SinglePhaseStaggered::getDensityDerivatives( double pressure, double temperature, double v2)
1181 {
1182         double rho=_fluides[0]->getDensity(pressure,temperature);
1183         double gamma=_fluides[0]->constante("gamma");
1184         double q=_fluides[0]->constante("q");
1185
1186         if(     !_useDellacherieEOS)
1187         {
1188                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1189                 double e = fluide0->getInternalEnergy(temperature);
1190                 double cv=fluide0->constante("cv");
1191                 double E=e+0.5*v2;
1192
1193                 _drho_sur_dp=1/((gamma-1)*(e-q));
1194                 _drho_sur_dT=-rho*cv/(e-q);
1195                 _drhoE_sur_dp=E/((gamma-1)*(e-q));
1196                 _drhoE_sur_dT=rho*cv*(1-E/(e-q));
1197         }
1198         else if(_useDellacherieEOS )
1199         {
1200                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1201                 double h=fluide0->getEnthalpy(temperature);
1202                 double H=h+0.5*v2;
1203                 double cp=fluide0->constante("cp");
1204
1205                 _drho_sur_dp=gamma/((gamma-1)*(h-q));
1206                 _drho_sur_dT=-rho*cp/(h-q);
1207                 _drhoE_sur_dp=gamma*H/((gamma-1)*(h-q))-1;
1208                 _drhoE_sur_dT=rho*cp*(1-H/(h-q));
1209         }
1210         else
1211         {
1212                 *_runLogFile<< "SinglePhaseStaggered::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
1213                 _runLogFile->close();
1214                 throw CdmathException("SinglePhaseStaggered::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
1215         }
1216
1217         if(_verbose && _nbTimeStep%_freqSave ==0)
1218         {
1219                 cout<<"_drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
1220                 cout<<"_drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
1221         }
1222 }
1223 void SinglePhaseStaggered::save(){
1224         string prim(_path+"/SinglePhaseStaggeredPrim_");///Results
1225         string cons(_path+"/SinglePhaseStaggeredCons_");
1226         prim+=_fileName;
1227         cons+=_fileName;
1228
1229         PetscInt Ii;
1230         for (long i = 0; i < _Nmailles; i++){
1231                 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
1232                 for (int j = 0; j < _nVar; j++){
1233                         Ii = i*_nVar +j;
1234                         VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
1235                 }
1236         }
1237         if(_saveConservativeField){
1238                 for (long i = 0; i < _Nmailles; i++){
1239                         // j = 0 : density; j = _nVar - 1 : energy j = 1,..,_nVar-2: momentum
1240                         for (int j = 0; j < _nVar; j++){
1241                                 Ii = i*_nVar +j;
1242                                 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
1243                         }
1244                 }
1245                 _UU.setTime(_time,_nbTimeStep);
1246         }
1247         _VV.setTime(_time,_nbTimeStep);
1248
1249         // create mesh and component info
1250         if (_nbTimeStep ==0){
1251                 string prim_suppress ="rm -rf "+prim+"_*";
1252                 string cons_suppress ="rm -rf "+cons+"_*";
1253
1254                 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
1255                 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
1256
1257                 if(_saveConservativeField){
1258                         _UU.setInfoOnComponent(0,"Density_(kg/m^3)");
1259                         _UU.setInfoOnComponent(1,"Momentum_x");// (kg/m^2/s)
1260                         if (_Ndim>1)
1261                                 _UU.setInfoOnComponent(2,"Momentum_y");// (kg/m^2/s)
1262                         if (_Ndim>2)
1263                                 _UU.setInfoOnComponent(3,"Momentum_z");// (kg/m^2/s)
1264
1265                         _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
1266
1267                         switch(_saveFormat)
1268                         {
1269                         case VTK :
1270                                 _UU.writeVTK(cons);
1271                                 break;
1272                         case MED :
1273                                 _UU.writeMED(cons);
1274                                 break;
1275                         case CSV :
1276                                 _UU.writeCSV(cons);
1277                                 break;
1278                         }
1279                 }
1280                 _VV.setInfoOnComponent(0,"Pressure_(Pa)");
1281                 _VV.setInfoOnComponent(1,"Velocity_x_(m/s)");
1282                 if (_Ndim>1)
1283                         _VV.setInfoOnComponent(2,"Velocity_y_(m/s)");
1284                 if (_Ndim>2)
1285                         _VV.setInfoOnComponent(3,"Velocity_z_(m/s)");
1286                 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
1287
1288                 switch(_saveFormat)
1289                 {
1290                 case VTK :
1291                         _VV.writeVTK(prim);
1292                         break;
1293                 case MED :
1294                         _VV.writeMED(prim);
1295                         break;
1296                 case CSV :
1297                         _VV.writeCSV(prim);
1298                         break;
1299                 }
1300         }
1301         // do not create mesh
1302         else{
1303                 switch(_saveFormat)
1304                 {
1305                 case VTK :
1306                         _VV.writeVTK(prim,false);
1307                         break;
1308                 case MED :
1309                         _VV.writeMED(prim,false);
1310                         break;
1311                 case CSV :
1312                         _VV.writeCSV(prim);
1313                         break;
1314                 }
1315                 if(_saveConservativeField){
1316                         switch(_saveFormat)
1317                         {
1318                         case VTK :
1319                                 _UU.writeVTK(cons,false);
1320                                 break;
1321                         case MED :
1322                                 _UU.writeMED(cons,false);
1323                                 break;
1324                         case CSV :
1325                                 _UU.writeCSV(cons);
1326                                 break;
1327                         }
1328                 }
1329         }
1330         if(_saveVelocity){
1331                 for (long i = 0; i < _Nmailles; i++){
1332                         // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
1333                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
1334                         {
1335                                 int Ii = i*_nVar +1+j;
1336                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
1337                         }
1338                         for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
1339                                 _Vitesse(i,j)=0;
1340                 }
1341                 _Vitesse.setTime(_time,_nbTimeStep);
1342                 if (_nbTimeStep ==0){
1343                         _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
1344                         _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
1345                         _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
1346
1347                         switch(_saveFormat)
1348                         {
1349                         case VTK :
1350                                 _Vitesse.writeVTK(prim+"_Velocity");
1351                                 break;
1352                         case MED :
1353                                 _Vitesse.writeMED(prim+"_Velocity");
1354                                 break;
1355                         case CSV :
1356                                 _Vitesse.writeCSV(prim+"_Velocity");
1357                                 break;
1358                         }
1359                 }
1360                 else{
1361                         switch(_saveFormat)
1362                         {
1363                         case VTK :
1364                                 _Vitesse.writeVTK(prim+"_Velocity",false);
1365                                 break;
1366                         case MED :
1367                                 _Vitesse.writeMED(prim+"_Velocity",false);
1368                                 break;
1369                         case CSV :
1370                                 _Vitesse.writeCSV(prim+"_Velocity");
1371                                 break;
1372                         }
1373                 }
1374         }
1375         if(_isStationary)
1376         {
1377                 prim+="_Stat";
1378                 cons+="_Stat";
1379
1380                 switch(_saveFormat)
1381                 {
1382                 case VTK :
1383                         _VV.writeVTK(prim);
1384                         break;
1385                 case MED :
1386                         _VV.writeMED(prim);
1387                         break;
1388                 case CSV :
1389                         _VV.writeCSV(prim);
1390                         break;
1391                 }
1392
1393                 if(_saveConservativeField){
1394                         switch(_saveFormat)
1395                         {
1396                         case VTK :
1397                                 _UU.writeVTK(cons);
1398                                 break;
1399                         case MED :
1400                                 _UU.writeMED(cons);
1401                                 break;
1402                         case CSV :
1403                                 _UU.writeCSV(cons);
1404                                 break;
1405                         }
1406                 }
1407
1408                 if(_saveVelocity){
1409                         switch(_saveFormat)
1410                         {
1411                         case VTK :
1412                                 _Vitesse.writeVTK(prim+"_Velocity");
1413                                 break;
1414                         case MED :
1415                                 _Vitesse.writeMED(prim+"_Velocity");
1416                                 break;
1417                         case CSV :
1418                                 _Vitesse.writeCSV(prim+"_Velocity");
1419                                 break;
1420                         }
1421                 }
1422         }
1423 }