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                         Polynoms Poly;
 654                         for( int i=0 ; i<3 ; i++)
 655                                 y[i] = Poly.abs_generalise(vp_dist[i])+_entropicShift[i];
 656                         Poly.abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
 657
 658                         if( _spaceScheme ==pressureCorrection)
 659                                 for( int i=0 ; i<_Ndim ; i++)
 660                                         for( int j=0 ; j<_Ndim ; j++)
 661                                                 _absAroe[(1+i)*_nVar+1+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 662                         else if( _spaceScheme ==lowMach){
 663                                 double M=sqrt(u_2)/c;
 664                                 for( int i=0 ; i<_Ndim ; i++)
 665                                         for( int j=0 ; j<_Ndim ; j++)
 666                                                 _absAroe[(1+i)*_nVar+1+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
 667                         }
 668                 }
 669                 else if( _spaceScheme ==staggered ){
 670                         if(_entropicCorrection)//To do: study entropic correction for staggered
 671                         {
 672                                 *_runLogFile<<"SinglePhaseStaggered::convectionMatrices: entropy scheme not available for staggered scheme"<<endl;
 673                                 _runLogFile->close();
 674                                 throw CdmathException("SinglePhaseStaggered::convectionMatrices: entropy scheme not available for staggered scheme");
 675                         }
 676
 677                         staggeredRoeUpwindingMatrixConservativeVariables( u_n, H, vitesse, k, K);
 678                 }
 679                 else
 680                 {
 681                         *_runLogFile<<"SinglePhaseStaggered::convectionMatrices: scheme not treated"<<endl;
 682                         _runLogFile->close();
 683                         throw CdmathException("SinglePhaseStaggered::convectionMatrices: scheme not treated");
 684                 }
 685
 686                 for(int i=0; i<_nVar*_nVar;i++)
 687                 {
 688                         _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
 689                         _AroePlus[i]  = (_Aroe[i]+_absAroe[i])/2;
 690                 }
 691                 if(_timeScheme==Implicit)
 692                 {
 693                         if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
 694                         {
 695                                 _Vij[0]=_fluides[0]->getPressureFromEnthalpy(_Uroe[_nVar-1]-u_2/2, _Uroe[0]);//pressure
 696                                 _Vij[_nVar-1]=_fluides[0]->getTemperatureFromPressure( _Vij[0], _Uroe[0]);//Temperature
 697                                 for(int idim=0;idim<_Ndim; idim++)
 698                                         _Vij[1+idim]=_Uroe[1+idim];
 699                                 primToConsJacobianMatrix(_Vij);
 700                                 Polynoms Poly;
 701                                 Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
 702                                 Poly.matrixProduct(_AroePlus,  _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
 703                         }
 704                         else
 705                                 for(int i=0; i<_nVar*_nVar;i++)
 706                                 {
 707                                         _AroeMinusImplicit[i] = _AroeMinus[i];
 708                                         _AroePlusImplicit[i]  = _AroePlus[i];
 709                                 }
 710                 }
 711                 if(_verbose && _nbTimeStep%_freqSave ==0)
 712                 {
 713                         displayMatrix(_Aroe, _nVar,"Matrice de Roe");
 714                         displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
 715                         displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
 716                         displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
 717                 }
 718         }
 719
 720         if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
 721         {
 722                 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
 723                 displayMatrix(_AroePlusImplicit,  _nVar,"Matrice _AroePlusImplicit");
 724         }
 725
 726         /*********Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source*****/
 727         if(_entropicCorrection)
 728         {
 729                 InvMatriceRoe( vp_dist);
 730                 Polynoms Poly;
 731                 Poly.matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
 732         }
 733         else if (_spaceScheme==upwind || (_spaceScheme ==pressureCorrection ) || (_spaceScheme ==lowMach ))//upwind sans entropic
 734                 SigneMatriceRoe( vp_dist);
 735         else if (_spaceScheme==centered)//centre  sans entropic
 736                 for(int i=0; i<_nVar*_nVar;i++)
 737                         _signAroe[i] = 0;
 738         else if( _spaceScheme ==staggered )//à tester
 739         {
 740                 double signu;
 741                 if(u_n>0)
 742                         signu=1;
 743                 else if (u_n<0)
 744                         signu=-1;
 745                 else
 746                         signu=0;
 747                 for(int i=0; i<_nVar*_nVar;i++)
 748                         _signAroe[i] = 0;
 749                 _signAroe[0] = signu;
 750                 for(int i=1; i<_nVar-1;i++)
 751                         _signAroe[i*_nVar+i] = -signu;
 752                 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
 753         }
 754         else
 755         {
 756                 *_runLogFile<<"SinglePhaseStaggered::convectionMatrices: well balanced option not treated"<<endl;
 757                 _runLogFile->close();
 758                 throw CdmathException("SinglePhaseStaggered::convectionMatrices: well balanced option not treated");
 759         }
 760 }
 761 void SinglePhaseStaggered::computeScaling(double maxvp)
 762 {
 763         _blockDiag[0]=1;
 764         _invBlockDiag[0]=1;
 765         for(int q=1; q<_nVar-1; q++)
 766         {
 767                 _blockDiag[q]=1./maxvp;//
 768                 _invBlockDiag[q]= maxvp;//1.;//
 769         }
 770         _blockDiag[_nVar - 1]=(_fluides[0]->constante("gamma")-1)/(maxvp*maxvp);//1
 771         _invBlockDiag[_nVar - 1]=  1./_blockDiag[_nVar - 1] ;// 1.;//
 772 }
 773
 774
 775 void SinglePhaseStaggered::sourceVector(PetscScalar * Si, PetscScalar * Ui, PetscScalar * Vi, int i)
 776 {
 777         double phirho=Ui[0], T=Vi[_nVar-1];
 778         double norm_u=0;
 779         for(int k=0; k<_Ndim; k++)
 780                 norm_u+=Vi[1+k]*Vi[1+k];
 781         norm_u=sqrt(norm_u);
 782         if(T>_Tsat)
 783                 Si[0]=_heatPowerField(i)/_latentHeat;
 784         else
 785                 Si[0]=0;
 786         for(int k=1; k<_nVar-1; k++)
 787                 Si[k]  =(_gravite[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*phirho;
 788
 789         Si[_nVar-1]=_heatPowerField(i);
 790
 791         for(int k=0; k<_Ndim; k++)
 792                 Si[_nVar-1] +=(_GravityField3d[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*Vi[1+k]*phirho;
 793
 794         if(_timeScheme==Implicit)
 795         {
 796                 for(int k=0; k<_nVar*_nVar;k++)
 797                         _GravityImplicitationMatrix[k] = 0;
 798                 if(!_usePrimitiveVarsInNewton)
 799                         for(int k=0; k<_nVar;k++)
 800                                 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
 801                 else
 802                 {
 803                         double pression=Vi[0];
 804                         getDensityDerivatives( pression, T, norm_u*norm_u);
 805                         for(int k=0; k<_nVar;k++)
 806                         {
 807                                 _GravityImplicitationMatrix[k*_nVar+0]      =-_gravite[k]*_drho_sur_dp;
 808                                 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
 809                         }
 810                 }
 811         }
 812         if(_verbose && _nbTimeStep%_freqSave ==0)
 813         {
 814                 cout<<"SinglePhaseStaggered::sourceVector"<<endl;
 815                 cout<<"Ui="<<endl;
 816                 for(int k=0;k<_nVar;k++)
 817                         cout<<Ui[k]<<", ";
 818                 cout<<endl;
 819                 cout<<"Vi="<<endl;
 820                 for(int k=0;k<_nVar;k++)
 821                         cout<<Vi[k]<<", ";
 822                 cout<<endl;
 823                 cout<<"Si="<<endl;
 824                 for(int k=0;k<_nVar;k++)
 825                         cout<<Si[k]<<", ";
 826                 cout<<endl;
 827                 if(_timeScheme==Implicit)
 828                         displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
 829         }
 830 }
 831
 832 void SinglePhaseStaggered::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
 833 {
 834         double norm_u=0, u_n=0, rho;
 835         for(int i=0;i<_Ndim;i++)
 836                 u_n += _Uroe[1+i]*_vec_normal[i];
 837
 838         pressureLoss[0]=0;
 839         if(u_n>0){
 840                 for(int i=0;i<_Ndim;i++)
 841                         norm_u += Vi[1+i]*Vi[1+i];
 842                 norm_u=sqrt(norm_u);
 843                 rho=Ui[0];
 844                 for(int i=0;i<_Ndim;i++)
 845                         pressureLoss[1+i]=-1/2*K*rho*norm_u*Vi[1+i];
 846         }
 847         else{
 848                 for(int i=0;i<_Ndim;i++)
 849                         norm_u += Vj[1+i]*Vj[1+i];
 850                 norm_u=sqrt(norm_u);
 851                 rho=Uj[0];
 852                 for(int i=0;i<_Ndim;i++)
 853                         pressureLoss[1+i]=-1/2*K*rho*norm_u*Vj[1+i];
 854         }
 855         pressureLoss[_nVar-1]=-1/2*K*rho*norm_u*norm_u*norm_u;
 856
 857         if(_verbose && _nbTimeStep%_freqSave ==0)
 858         {
 859                 cout<<"SinglePhaseStaggered::pressureLossVector K= "<<K<<endl;
 860                 cout<<"Ui="<<endl;
 861                 for(int k=0;k<_nVar;k++)
 862                         cout<<Ui[k]<<", ";
 863                 cout<<endl;
 864                 cout<<"Vi="<<endl;
 865                 for(int k=0;k<_nVar;k++)
 866                         cout<<Vi[k]<<", ";
 867                 cout<<endl;
 868                 cout<<"Uj="<<endl;
 869                 for(int k=0;k<_nVar;k++)
 870                         cout<<Uj[k]<<", ";
 871                 cout<<endl;
 872                 cout<<"Vj="<<endl;
 873                 for(int k=0;k<_nVar;k++)
 874                         cout<<Vj[k]<<", ";
 875                 cout<<endl;
 876                 cout<<"pressureLoss="<<endl;
 877                 for(int k=0;k<_nVar;k++)
 878                         cout<<pressureLoss[k]<<", ";
 879                 cout<<endl;
 880         }
 881 }
 882
 883 void SinglePhaseStaggered::porosityGradientSourceVector()
 884 {
 885         double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[0], pj=_Vj[0],pij;
 886         for(int i=0;i<_Ndim;i++) {
 887                 u_ni += _Vi[1+i]*_vec_normal[i];
 888                 u_nj += _Vj[1+i]*_vec_normal[i];
 889         }
 890         _porosityGradientSourceVector[0]=0;
 891         rhoj=_Uj[0]/_porosityj;
 892         rhoi=_Ui[0]/_porosityi;
 893         pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
 894         for(int i=0;i<_Ndim;i++)
 895                 _porosityGradientSourceVector[1+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
 896         _porosityGradientSourceVector[_nVar-1]=0;
 897 }
 898
 899 void SinglePhaseStaggered::jacobian(const int &j, string nameOfGroup,double * normale)
 900 {
 901         if(_verbose && _nbTimeStep%_freqSave ==0)
 902                 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
 903
 904         int k;
 905         for(k=0; k<_nVar*_nVar;k++)
 906                 _Jcb[k] = 0;//No implicitation at this stage
 907
 908         _idm[0] = _nVar*j;
 909         for(k=1; k<_nVar; k++)
 910                 _idm[k] = _idm[k-1] + 1;
 911         VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
 912         double q_n=0;//quantité de mouvement normale à la paroi
 913         for(k=0; k<_Ndim; k++)
 914                 q_n+=_externalStates[(k+1)]*normale[k];
 915
 916         // loop of boundary types
 917         if (_limitField[nameOfGroup].bcType==Wall)
 918         {
 919                 for(k=0; k<_nVar;k++)
 920                         _Jcb[k*_nVar + k] = 1;
 921                 for(k=1; k<_nVar-1;k++)
 922                         for(int l=1; l<_nVar-1;l++)
 923                                 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
 924         }
 925         else if (_limitField[nameOfGroup].bcType==Inlet)
 926         {
 927                 return;
 928                 if(q_n<0){
 929                         double v[_Ndim], ve[_Ndim], v2, ve2;
 930
 931                         _idm[0] = _nVar*j;
 932                         for(k=1; k<_nVar; k++)
 933                                 _idm[k] = _idm[k-1] + 1;
 934                         VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 935                         VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
 936
 937                         ve[0] = _limitField[nameOfGroup].v_x[0];
 938                         v[0]=_Vj[1];
 939                         ve2 = ve[0]*ve[0];
 940                         v2 = v[0]*v[0];
 941                         if (_Ndim >1){
 942                                 ve[1] = _limitField[nameOfGroup].v_y[0];
 943                                 v[1]=_Vj[2];
 944                                 ve2 += ve[1]*ve[1];
 945                                 v2 = v[1]*v[1];
 946                         }
 947                         if (_Ndim >2){
 948                                 ve[2] = _limitField[nameOfGroup].v_z[0];
 949                                 v[2]=_Vj[3];
 950                                 ve2 += ve[2]*ve[2];
 951                                 v2 = v[2]*v[2];
 952                         }
 953                         double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,_Uj[0]);
 954                         double total_energy=internal_energy+ve2/2;
 955
 956                         //Mass line
 957                         _Jcb[0]=v2/(2*internal_energy);
 958                         for(k=0; k<_Ndim;k++)
 959                                 _Jcb[1+k]=-v[k]/internal_energy;
 960                         _Jcb[_nVar-1]=1/internal_energy;
 961                         //Momentum lines
 962                         for(int l =1;l<1+_Ndim;l++){
 963                                 _Jcb[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
 964                                 for(k=0; k<_Ndim;k++)
 965                                         _Jcb[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
 966                                 _Jcb[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
 967                         }
 968                         //Energy line
 969                         _Jcb[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
 970                         for(k=0; k<_Ndim;k++)
 971                                 _Jcb[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
 972                         _Jcb[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
 973                 }
 974                 else
 975                         for(k=0;k<_nVar;k++)
 976                                 _Jcb[k*_nVar+k]=1;
 977                 //Old jacobian
 978                 /*
 979                  _Jcb[0] = 1;
 980                 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];//Kieu
 981                 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
 982                 if(_Ndim>1)
 983                 {
 984                         _Jcb[2*_nVar]= _limitField[nameOfGroup].v_y[0];
 985                         v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
 986                         if(_Ndim==3){
 987                                 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
 988                                 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
 989                         }
 990                 }
 991                 _Jcb[(_nVar-1)*_nVar]=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2;
 992                  */
 993         }
 994         else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0){
 995                 return;
 996                 double v[_Ndim], v2=0;
 997                 _idm[0] = _nVar*j;
 998                 for(k=1; k<_nVar; k++)
 999                         _idm[k] = _idm[k-1] + 1;
1000                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1001
1002                 for(k=0; k<_Ndim;k++){
1003                         v[k]=_Vj[1+k];
1004                         v2+=v[k]*v[k];
1005                 }
1006
1007                 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _limitField[nameOfGroup].T);
1008                 double rho_int = _externalStates[0];
1009                 double density_ratio=rho_ext/rho_int;
1010                 //Momentum lines
1011                 for(int l =1;l<1+_Ndim;l++){
1012                         _Jcb[l*_nVar]=-density_ratio*v[l-1];
1013                         _Jcb[l*_nVar+l]=density_ratio;
1014                 }
1015                 //Energy lines
1016                 _Jcb[(_nVar-1)*_nVar]=-v2*density_ratio;
1017                 for(k=0; k<_Ndim;k++)
1018                         _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k];
1019         }
1020         // not wall, not inlet, not inletPressure
1021         else if(_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0))
1022         {
1023                 return;
1024                 double v[_Ndim], v2=0;
1025                 _idm[0] = _nVar*j;
1026                 for(k=1; k<_nVar; k++)
1027                         _idm[k] = _idm[k-1] + 1;
1028                 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1029
1030                 for(k=0; k<_Ndim;k++){
1031                         v[k]=_Vj[1+k];
1032                         v2+=v[k]*v[k];
1033                 }
1034
1035                 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _externalStates[_nVar-1]);
1036                 double rho_int = _externalStates[0];
1037                 double density_ratio=rho_ext/rho_int;
1038                 double internal_energy=_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho_int);
1039                 double total_energy=internal_energy+v2/2;
1040
1041                 //Mass line
1042                 _Jcb[0]=density_ratio*(1-v2/(2*internal_energy));
1043                 for(k=0; k<_Ndim;k++)
1044                         _Jcb[1+k]=density_ratio*v[k]/internal_energy;
1045                 _Jcb[_nVar-1]=-density_ratio/internal_energy;
1046                 //Momentum lines
1047                 for(int l =1;l<1+_Ndim;l++){
1048                         _Jcb[l*_nVar]=density_ratio*v2*v[l-1]/(2*internal_energy);
1049                         for(k=0; k<_Ndim;k++)
1050                                 _Jcb[l*_nVar+1+k]=density_ratio*v[k]*v[l-1]/internal_energy;
1051                         _Jcb[l*_nVar+1+k]-=density_ratio;
1052                         _Jcb[l*_nVar+_nVar-1]=-density_ratio*v[l-1]/internal_energy;
1053                 }
1054                 //Energy line
1055                 _Jcb[(_nVar-1)*_nVar]=density_ratio*v2*total_energy/(2*internal_energy);
1056                 for(k=0; k<_Ndim;k++)
1057                         _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k]*total_energy/internal_energy;
1058                 _Jcb[(_nVar-1)*_nVar+_nVar-1]=density_ratio*(1-total_energy/internal_energy);
1059                 //Old jacobian
1060                 /*
1061                 int idim,jdim;
1062                 double cd = 1,cn=0,p0, gamma;
1063                 _idm[0] = j*_nVar;// Kieu
1064                 for(k=1; k<_nVar;k++)
1065                         _idm[k] = _idm[k-1] + 1;
1066                 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1067                 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1068
1069                 // compute the common numerator and common denominator
1070                 p0=_fluides[0]->constante("p0");
1071                 gamma =_fluides[0]->constante("gamma");
1072                 cn =_limitField[nameOfGroup].p +p0;
1073                 cd = _phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0;
1074                 cd*=cd;
1075                 cd*=(gamma-1);
1076                 //compute the v2
1077                 for(k=1; k<_nVar-1;k++)
1078                         v2+=_externalStates[k]*_externalStates[k];
1079                 // drho_ext/dU
1080                 _JcbDiff[0] = cn*(_phi[_nVar-1] -v2 -p0)/cd;
1081                 for(k=1; k<_nVar-1;k++)
1082                         _JcbDiff[k]=cn*_phi[k]/cd;
1083                 _JcbDiff[_nVar-1]= -cn*_phi[0]/cd;
1084                 //dq_ext/dU
1085                 for(idim=0; idim<_Ndim;idim++)
1086                 {
1087                         //premiere colonne
1088                         _JcbDiff[(1+idim)*_nVar]=-(v2*cn*_phi[idim+1])/(2*cd);
1089                         //colonnes intermediaire
1090                         for(jdim=0; jdim<_Ndim;jdim++)
1091                         {
1092                                 _JcbDiff[(1+idim)*_nVar + jdim + 1] =_externalStates[idim+1]*_phi[jdim+1];
1093                                 _JcbDiff[(1+idim)*_nVar + jdim + 1]*=cn/cd;
1094                         }
1095                         //matrice identite*cn*(rhoe- p0)
1096                         _JcbDiff[(1+idim)*_nVar + idim + 1] +=( cn*(_phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0))/cd;
1097
1098                         //derniere colonne
1099                         _JcbDiff[(1+idim)*_nVar + _nVar-1]=-_phi[idim+1]*cn/cd;
1100                 }
1101                 //drhoE/dU
1102                 _JcbDiff[_nVar*(_nVar-1)] = -(v2*_phi[_nVar -1]*cn)/(2*cd);
1103                 for(int idim=0; idim<_Ndim;idim++)
1104                         _JcbDiff[_nVar*(_nVar-1)+idim+1]=_externalStates[idim+1]*_phi[_nVar -1]*cn/cd;
1105                 _JcbDiff[_nVar*_nVar -1] = -(v2/2+p0)*cn/cd;
1106                  */
1107         }
1108         else  if (_limitField[nameOfGroup].bcType!=Neumann)// not wall, not inlet, not outlet
1109         {
1110                 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1111                 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1112                 _runLogFile->close();
1113                 throw CdmathException("SinglePhaseStaggered::jacobian: This boundary condition is not treated");
1114         }
1115 }
1116
1117
1118 Vector SinglePhaseStaggered::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
1119         if(_verbose && _nbTimeStep%_freqSave ==0)
1120         {
1121                 cout<<"SinglePhaseStaggered::convectionFlux start"<<endl;
1122                 cout<<"Ucons"<<endl;
1123                 cout<<U<<endl;
1124                 cout<<"Vprim"<<endl;
1125                 cout<<V<<endl;
1126         }
1127
1128         double phirho=U(0);//phi rho
1129         Vector phiq(_Ndim);//phi rho u
1130         for(int i=0;i<_Ndim;i++)
1131                 phiq(i)=U(1+i);
1132
1133         double pression=V(0);
1134         Vector vitesse(_Ndim);
1135         for(int i=0;i<_Ndim;i++)
1136                 vitesse(i)=V(1+i);
1137         double Temperature= V(1+_Ndim);
1138
1139         double vitessen=vitesse*normale;
1140         double rho=phirho/porosity;
1141         double e_int=_fluides[0]->getInternalEnergy(Temperature,rho);
1142
1143         Vector F(_nVar);
1144         F(0)=phirho*vitessen;
1145         for(int i=0;i<_Ndim;i++)
1146                 F(1+i)=phirho*vitessen*vitesse(i)+pression*normale(i)*porosity;
1147         F(1+_Ndim)=phirho*(e_int+0.5*vitesse*vitesse+pression/rho)*vitessen;
1148
1149         if(_verbose && _nbTimeStep%_freqSave ==0)
1150         {
1151                 cout<<"SinglePhaseStaggered::convectionFlux end"<<endl;
1152                 cout<<"Flux F(U,V)"<<endl;
1153                 cout<<F<<endl;
1154         }
1155
1156         return F;
1157 }
1158
1159 void SinglePhaseStaggered::convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector vitesse)
1160 {
1161         //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1162         //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
1163         //EOS is more involved with primitive variables
1164         // call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
1165         _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1166         for(int i=0;i<_Ndim;i++)
1167                 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1168         _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1169         for(int i=0;i<_Ndim;i++)
1170         {
1171                 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]+_vec_normal[i];
1172                 for(int j=0;j<_Ndim;j++)
1173                         _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1174                 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1175                 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1176         }
1177         _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1178         for(int i=0;i<_Ndim;i++)
1179                 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1180         _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1181 }
1182
1183 void SinglePhaseStaggered::getDensityDerivatives( double pressure, double temperature, double v2)
1184 {
1185         double rho=_fluides[0]->getDensity(pressure,temperature);
1186         double gamma=_fluides[0]->constante("gamma");
1187         double q=_fluides[0]->constante("q");
1188
1189         if(     !_useDellacherieEOS)
1190         {
1191                 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1192                 double e = fluide0->getInternalEnergy(temperature);
1193                 double cv=fluide0->constante("cv");
1194                 double E=e+0.5*v2;
1195
1196                 _drho_sur_dp=1/((gamma-1)*(e-q));
1197                 _drho_sur_dT=-rho*cv/(e-q);
1198                 _drhoE_sur_dp=E/((gamma-1)*(e-q));
1199                 _drhoE_sur_dT=rho*cv*(1-E/(e-q));
1200         }
1201         else if(_useDellacherieEOS )
1202         {
1203                 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1204                 double h=fluide0->getEnthalpy(temperature);
1205                 double H=h+0.5*v2;
1206                 double cp=fluide0->constante("cp");
1207
1208                 _drho_sur_dp=gamma/((gamma-1)*(h-q));
1209                 _drho_sur_dT=-rho*cp/(h-q);
1210                 _drhoE_sur_dp=gamma*H/((gamma-1)*(h-q))-1;
1211                 _drhoE_sur_dT=rho*cp*(1-H/(h-q));
1212         }
1213         else
1214         {
1215                 *_runLogFile<< "SinglePhaseStaggered::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
1216                 _runLogFile->close();
1217                 throw CdmathException("SinglePhaseStaggered::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
1218         }
1219
1220         if(_verbose && _nbTimeStep%_freqSave ==0)
1221         {
1222                 cout<<"_drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
1223                 cout<<"_drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
1224         }
1225 }
1226 void SinglePhaseStaggered::save(){
1227         string prim(_path+"/SinglePhaseStaggeredPrim_");///Results
1228         string cons(_path+"/SinglePhaseStaggeredCons_");
1229         prim+=_fileName;
1230         cons+=_fileName;
1231
1232         PetscInt Ii;
1233         for (long i = 0; i < _Nmailles; i++){
1234                 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
1235                 for (int j = 0; j < _nVar; j++){
1236                         Ii = i*_nVar +j;
1237                         VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
1238                 }
1239         }
1240         if(_saveConservativeField){
1241                 for (long i = 0; i < _Nmailles; i++){
1242                         // j = 0 : density; j = _nVar - 1 : energy j = 1,..,_nVar-2: momentum
1243                         for (int j = 0; j < _nVar; j++){
1244                                 Ii = i*_nVar +j;
1245                                 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
1246                         }
1247                 }
1248                 _UU.setTime(_time,_nbTimeStep);
1249         }
1250         _VV.setTime(_time,_nbTimeStep);
1251
1252         // create mesh and component info
1253         if (_nbTimeStep ==0){
1254                 string prim_suppress ="rm -rf "+prim+"_*";
1255                 string cons_suppress ="rm -rf "+cons+"_*";
1256
1257                 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
1258                 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
1259
1260                 if(_saveConservativeField){
1261                         _UU.setInfoOnComponent(0,"Density_(kg/m^3)");
1262                         _UU.setInfoOnComponent(1,"Momentum_x");// (kg/m^2/s)
1263                         if (_Ndim>1)
1264                                 _UU.setInfoOnComponent(2,"Momentum_y");// (kg/m^2/s)
1265                         if (_Ndim>2)
1266                                 _UU.setInfoOnComponent(3,"Momentum_z");// (kg/m^2/s)
1267
1268                         _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
1269
1270                         switch(_saveFormat)
1271                         {
1272                         case VTK :
1273                                 _UU.writeVTK(cons);
1274                                 break;
1275                         case MED :
1276                                 _UU.writeMED(cons);
1277                                 break;
1278                         case CSV :
1279                                 _UU.writeCSV(cons);
1280                                 break;
1281                         }
1282                 }
1283                 _VV.setInfoOnComponent(0,"Pressure_(Pa)");
1284                 _VV.setInfoOnComponent(1,"Velocity_x_(m/s)");
1285                 if (_Ndim>1)
1286                         _VV.setInfoOnComponent(2,"Velocity_y_(m/s)");
1287                 if (_Ndim>2)
1288                         _VV.setInfoOnComponent(3,"Velocity_z_(m/s)");
1289                 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
1290
1291                 switch(_saveFormat)
1292                 {
1293                 case VTK :
1294                         _VV.writeVTK(prim);
1295                         break;
1296                 case MED :
1297                         _VV.writeMED(prim);
1298                         break;
1299                 case CSV :
1300                         _VV.writeCSV(prim);
1301                         break;
1302                 }
1303         }
1304         // do not create mesh
1305         else{
1306                 switch(_saveFormat)
1307                 {
1308                 case VTK :
1309                         _VV.writeVTK(prim,false);
1310                         break;
1311                 case MED :
1312                         _VV.writeMED(prim,false);
1313                         break;
1314                 case CSV :
1315                         _VV.writeCSV(prim);
1316                         break;
1317                 }
1318                 if(_saveConservativeField){
1319                         switch(_saveFormat)
1320                         {
1321                         case VTK :
1322                                 _UU.writeVTK(cons,false);
1323                                 break;
1324                         case MED :
1325                                 _UU.writeMED(cons,false);
1326                                 break;
1327                         case CSV :
1328                                 _UU.writeCSV(cons);
1329                                 break;
1330                         }
1331                 }
1332         }
1333         if(_saveVelocity){
1334                 for (long i = 0; i < _Nmailles; i++){
1335                         // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
1336                         for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
1337                         {
1338                                 int Ii = i*_nVar +1+j;
1339                                 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
1340                         }
1341                         for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
1342                                 _Vitesse(i,j)=0;
1343                 }
1344                 _Vitesse.setTime(_time,_nbTimeStep);
1345                 if (_nbTimeStep ==0){
1346                         _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
1347                         _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
1348                         _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
1349
1350                         switch(_saveFormat)
1351                         {
1352                         case VTK :
1353                                 _Vitesse.writeVTK(prim+"_Velocity");
1354                                 break;
1355                         case MED :
1356                                 _Vitesse.writeMED(prim+"_Velocity");
1357                                 break;
1358                         case CSV :
1359                                 _Vitesse.writeCSV(prim+"_Velocity");
1360                                 break;
1361                         }
1362                 }
1363                 else{
1364                         switch(_saveFormat)
1365                         {
1366                         case VTK :
1367                                 _Vitesse.writeVTK(prim+"_Velocity",false);
1368                                 break;
1369                         case MED :
1370                                 _Vitesse.writeMED(prim+"_Velocity",false);
1371                                 break;
1372                         case CSV :
1373                                 _Vitesse.writeCSV(prim+"_Velocity");
1374                                 break;
1375                         }
1376                 }
1377         }
1378         if(_isStationary)
1379         {
1380                 prim+="_Stat";
1381                 cons+="_Stat";
1382
1383                 switch(_saveFormat)
1384                 {
1385                 case VTK :
1386                         _VV.writeVTK(prim);
1387                         break;
1388                 case MED :
1389                         _VV.writeMED(prim);
1390                         break;
1391                 case CSV :
1392                         _VV.writeCSV(prim);
1393                         break;
1394                 }
1395
1396                 if(_saveConservativeField){
1397                         switch(_saveFormat)
1398                         {
1399                         case VTK :
1400                                 _UU.writeVTK(cons);
1401                                 break;
1402                         case MED :
1403                                 _UU.writeMED(cons);
1404                                 break;
1405                         case CSV :
1406                                 _UU.writeCSV(cons);
1407                                 break;
1408                         }
1409                 }
1410
1411                 if(_saveVelocity){
1412                         switch(_saveFormat)
1413                         {
1414                         case VTK :
1415                                 _Vitesse.writeVTK(prim+"_Velocity");
1416                                 break;
1417                         case MED :
1418                                 _Vitesse.writeMED(prim+"_Velocity");
1419                                 break;
1420                         case CSV :
1421                                 _Vitesse.writeCSV(prim+"_Velocity");
1422                                 break;
1423                         }
1424                 }
1425         }
1426 }