CoreFlows/Models/inc/ProblemFluid.hxx

   1 //============================================================================
   2 // Name        : ProblemCoreFlows
   3 // Author      : M. Ndjinga
   4 // Version     :
   5 // Copyright   : CEA Saclay 2014
   6 // Description : Generic class for thermal hydraulics problems
   7 //============================================================================
   8
   9 /*! \class ProblemFluid ProblemFluid.hxx "ProblemFluid.hxx"
  10  *  \brief Factorises the methods that are common to the non scalar models (fluid models)
  11  *  \details Common functions to fluid models
  12  */
  13 #ifndef PROBLEMFLUID_HXX_
  14 #define PROBLEMFLUID_HXX_
  15
  16 #include "Fluide.h"
  17 #include "ProblemCoreFlows.hxx"
  18 #include "utilitaire_algebre.h"
  19
  20 using namespace std;
  21
  22 //! enumeration NonLinearFormulation
  23 /*! the formulation used to compute the non viscous fluxes */
  24 enum NonLinearFormulation
  25 {
  26         Roe,/**< Ph. Roe non linear formulation is used */
  27         VFRoe,/**< Masella, Faille and Gallouet non linear formulation is used */
  28         VFFC,/**< Ghidaglia, Kumbaro and Le Coq non linear formulation is used */
  29         reducedRoe,/**< compacted formulation of Roe scheme without computation of the fluxes */
  30 };
  31
  32 class ProblemFluid: public ProblemCoreFlows
  33 {
  34
  35 public :
  36         /**\fn
  37          * \brief constructeur de la classe ProblemFluid
  38          */
  39         ProblemFluid(void);
  40
  41         //Gestion du calcul (interface ICoCo)
  42
  43         /** \fn initialize
  44          * \brief Alocates memory and checks that the mesh, the boundary and the intitial data are set
  45          * \Details It is a pure virtual function overloaded y each model.
  46          * @param  void
  47          *  */
  48         virtual void initialize();
  49
  50         /** \fn terminate
  51          * \brief empties the memory
  52          * @param void
  53          *  */
  54         virtual void terminate();
  55
  56         /** \fn initTimeStep
  57          * \brief Sets a new time step dt to be solved later
  58          *  @param  dt is  the value of the time step
  59          *  \return false if dt <0 et True otherwise
  60          *                         */
  61         bool initTimeStep(double dt);
  62
  63         /** \fn computeTimeStep
  64          * \brief Proposes a value for the next time step to be solved using mesh data and cfl coefficient
  65          *  \return  double dt the proposed time step
  66          *  \return  bool stop, true if the calculation should not be continued (stationary state, maximum time or time step numer reached)
  67          *  */
  68         double computeTimeStep(bool & stop);
  69
  70         /** \fn abortTimeStep
  71          * \brief Reset the time step dt to 0
  72          *  */
  73         void abortTimeStep();
  74
  75         /** \fn iterateTimeStep
  76          * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
  77          * @param
  78          * @return boolean ok is true is the newton iteration gave a physically acceptable result
  79          * */
  80         virtual bool iterateTimeStep(bool &ok);
  81
  82         /** \fn save
  83          * \brief saves the current results in MED or VTK files
  84          * \details It is a pure virtual function overloaded in each model
  85          * @param  void
  86          */
  87         virtual void save()=0;
  88
  89         /** \fn validateTimeStep
  90          * \brief Validates the solution computed y solveTimeStep
  91          * \details updates the currens time t=t+dt, save unknown fields, resets the time step dt to 0, tests the stationnarity.
  92          * c It is a pure virtual function overloaded in each model
  93          * @param  void
  94          *  */
  95         virtual void validateTimeStep();
  96
  97         /** \fn setOutletBoundaryCondition
  98          * \brief Adds a new boundary condition of type Outlet
  99          * \details
 100          * \param [in] string : the name of the boundary
 101          * \param [in] double : the value of the pressure at the boundary
 102          * \param [out] void
 103          *  */
 104         void setOutletBoundaryCondition(string groupName,double Pressure){
 105                 _limitField[groupName]=LimitField(Outlet,Pressure,vector<double>(_nbPhases,0),vector<double>(_nbPhases,0),vector<double>(_nbPhases,0),-1,-1,-1,-1);
 106         };
 107
 108         /** \fn setOutletBoundaryCondition
 109          * \brief Adds a new boundary condition of type Outlet taking into account the hydrostatic pressure variations
 110          * \details The pressure is not constant on the boundary but varies linearly with a slope given by the gravity vector
 111          * \param [in] string : the name of the boundary
 112          * \param [in] double : the value of the pressure at the boundary
 113          * \param [in] vector<double> : reference_point position on the boundary where the value Pressure will be imposed
 114          * \param [out] void
 115          *  */
 116         void setOutletBoundaryCondition(string groupName,double referencePressure, vector<double> reference_point){
 117                 /* On the boundary we have P-Pref=rho g\cdot(x-xref) */
 118                 _gravityReferencePoint=reference_point;
 119                 _limitField[groupName]=LimitField(Outlet,referencePressure,vector<double>(_nbPhases,0),vector<double>(_nbPhases,0),vector<double>(_nbPhases,0),-1,-1,-1,-1);
 120         };
 121
 122         /** \fn setViscosity
 123          * \brief sets the vector of viscosity coefficients
 124          * @param viscosite is a vector of size equal to the number of phases and containing the viscosity of each phase
 125          * @return throws an exception if the input vector size is not equal to the number of phases
 126          *       * */
 127         void setViscosity(vector<double> viscosite){
 128                 if(_nbPhases!= viscosite.size())
 129                         throw CdmathException("ProblemFluid::setViscosity: incorrect vector size vs number of phases");
 130                 for(int i=0;i<_nbPhases;i++)
 131                         _fluides[i]->setViscosity(viscosite[i]);
 132         };
 133
 134         /** \fn setConductivity
 135          * \brief sets the vector of conductivity coefficients
 136          * @param conductivite is a vector of size equal to the number of phases and containing the conductivity of each phase
 137          * @return throws an exception if the input vector size is not equal to the number of phases
 138          * */
 139         void setConductivity(vector<double> conductivite){
 140                 if(_nbPhases!= conductivite.size())
 141                         throw CdmathException("ProblemFluid::setConductivity: incorrect vector size vs number of phases");
 142                 for(int i=0;i<_nbPhases;i++)
 143                         _fluides[i]->setConductivity(conductivite[i]);
 144         };
 145
 146         /** \fn setGravity
 147          * \brief sets the gravity force in the model
 148          * @param gravite is a vector of size equal to the space dimension
 149          *                               * */
 150         void setGravity(vector<double> gravite){
 151                 _GravityField3d = gravite;
 152         };
 153
 154         /** \fn setDragCoeffs
 155          * \brief Sets the drag coefficients
 156          * @param dragCoeffs is a  vector of size equal to the number of phases and containing the value of the friction coefficient of each phase
 157          * @return throws an exception if the input vector size is not equal to the numer of phases
 158          * */
 159         void setDragCoeffs(vector<double> dragCoeffs){
 160                 if(_nbPhases!= dragCoeffs.size())
 161                         throw CdmathException("ProblemFluid::setDragCoeffs: incorrect vector size vs number of phases");
 162                 for(int i=0;i<_nbPhases;i++)
 163                 {
 164                         _fluides[i]->setDragCoeffs(dragCoeffs[i]);
 165                         _dragCoeffs[i]=dragCoeffs[i];
 166                 }
 167         };
 168
 169         /** \fn getNumberOfPhases
 170          * \brief The numer of phase (one or two) depending on the model considered
 171          * @param void
 172          * @return the number of phases considered in the model
 173          * */
 174         int getNumberOfPhases(){
 175                 return _nbPhases;
 176         };
 177
 178         /** \fn computeNewtonVariation
 179          * \brief Builds and solves the linear system to obtain the variation Ukp1-Uk in a Newton scheme
 180          * @param void
 181          * */
 182         virtual void computeNewtonVariation();
 183
 184         /** \fn testConservation
 185          * \brief Teste et affiche la conservation de masse et de la quantité de mouvement
 186          * \Details la fonction est virtuelle pure, on la surcharge dans chacun des modèles
 187          * @param void
 188          * */
 189         virtual void testConservation()=0;
 190
 191         /** \fn saveVelocity
 192          * \brief saves the velocity field in a separate 3D file so that paraview can display the streamlines
 193          * @param bool
 194          * */
 195         void saveVelocity(bool save_v=true){
 196                 _saveVelocity=save_v;
 197         }
 198
 199         /** \fn saveConservativeField
 200          * \brief saves the conservative fields (density, momentum etc...)
 201          * */
 202         void saveConservativeField(bool save=true){
 203                 _saveConservativeField=save;
 204         }
 205         /** \fn setEntropicCorrection
 206          * \brief include an entropy correction to avoid non entropic solutions
 207          * @param boolean that is true if entropy correction should be applied
 208          * */
 209         void setEntropicCorrection(bool entropyCorr){
 210                 _entropicCorrection=entropyCorr;
 211         }
 212
 213         /** \fn setPressureCorrectionOrder
 214          * \brief In case a pressure correction scheme is set by a call to setNumericalScheme(pressureCorrection) this function allows the setting of the type of pressure correction to be used
 215          * @param int the order of the pressure correction
 216         * \details The following treatment is applied depending on the value of the input parameter order
 217         * \details 1 -> no pressure correction (pressure correction is applied nowhere in the domain), standard upwind scheme instead is used
 218         * \details 2 -> standard pressure correction is applied everywhere in the domain, even at the boundary faces
 219         * \details 3 -> standard pressure correction is applied only inside the domain (not at the boundary faces)
 220         * \details 4 -> no pressure correction (pressure correction is applied nowhere in the domain), no Riemann problem at wall boundaries (boundary pressure = inner pressure)
 221         * \details 5 -> standard pressure correction is applied everywhere in the domain, no Riemann problem at the boundary (boundary pressure = inner pressure)
 222         * \details 6 -> standard pressure correction is applied inside the domain and a special pressure correction involving gravity is applied at the boundary, no Riemann problem at wall boundaries (boundary pressure = inner pressure+ source term)
 223         * */
 224         void setPressureCorrectionOrder(int order){
 225                 if( order >6 || order <1)
 226                         throw CdmathException("ProblemFluid::setPressureCorrectionOrder Pressure correction order must be an integer between 1 and 4");
 227                 else
 228                         _pressureCorrectionOrder=order;
 229
 230                 if(order==1)
 231                         _spaceScheme=upwind;
 232         }
 233
 234         // Petsc resolution
 235
 236         /** \fn setLinearSolver
 237          * \brief sets the linear solver and preconditioner
 238          * \details virtuelle function overloaded by intanciable classes
 239          * @param kspType linear solver type (GMRES or BICGSTAB)
 240          * @param pcType preconditioner (ILU,LU or NONE)
 241          * @param scaling performs a bancing of the system matrix before calling th preconditioner
 242          */
 243         void setLinearSolver(linearSolver kspType, preconditioner pcType, bool scaling=false)
 244         {
 245                 ProblemCoreFlows::setLinearSolver(kspType, pcType);
 246                 _isScaling= scaling;
 247         };
 248
 249         /** \fn setLatentHeat
 250          * \brief Sets the value of the latent heat
 251          * @param double L, the value of the latent heat
 252          * */
 253         void setLatentHeat(double L){
 254                 _latentHeat=L;
 255         }
 256
 257         /** \fn getLatentHeat
 258          * \brief returns the value of the latent heat
 259          * @param double L, the value of the latent heat
 260          * */
 261         double getLatentHeat(){
 262                 return _latentHeat;
 263         }
 264
 265         /** \fn setSatTemp
 266          * \brief sets the saturation temperature
 267          * @param  Tsat double corresponds to saturation temperature
 268          * */
 269         void setSatTemp(double Tsat){
 270                 _Tsat=Tsat;
 271         }
 272
 273         /** \fn setSatTemp
 274          * \brief sets the saturation temperature
 275          * @param  Tsat double corresponds to saturation temperature
 276          * */
 277         double getSatTemp(){
 278                 return _Tsat;
 279         }
 280
 281         /** \fn setSatPressure
 282          * \brief sets the saturation pressure
 283          * @param  Psat double corresponds to saturation pressure
 284          * */
 285         void setSatPressure(double Psat, double dHsatl_over_dp=0.05){
 286                 _Psat=Psat;
 287                 _dHsatl_over_dp=dHsatl_over_dp;
 288         }
 289
 290         /** \fn setPorosityField
 291          * \brief set the porosity field;
 292          * @param [in] Field porosity field (field on CELLS)
 293          * */
 294         void setPorosityField(Field Porosity){
 295                 _porosityField=Porosity;
 296                 _porosityFieldSet=true;
 297         }
 298
 299         /** \fn getPorosityField
 300          * \brief returns the porosity field;
 301          * @param
 302          * */
 303         Field getPorosityField(){
 304                 return _porosityField;
 305         }
 306
 307         /** \fn setPorosityFieldFile
 308          * \brief set the porosity field
 309          * \details
 310          * \param [in] string fileName (including file path)
 311          * \param [in] string fieldName
 312          * \param [out] void
 313          *  */
 314         void setPorosityField(string fileName, string fieldName){
 315                 _porosityField=Field(fileName, CELLS,fieldName);
 316                 _porosityFieldSet=true;
 317         }
 318
 319         /** \fn setPressureLossField
 320          * \brief set the pressure loss coefficients field;
 321          * @param [in] Field pressure loss field (field on FACES)
 322          * */
 323         void setPressureLossField(Field PressureLoss){
 324                 _pressureLossField=PressureLoss;
 325                 _pressureLossFieldSet=true;
 326         }
 327         /** \fn setPressureLossField
 328          * \brief set the pressure loss coefficient field
 329          * \details localised friction force
 330          * \param [in] string fileName (including file path)
 331          * \param [in] string fieldName
 332          * \param [out] void
 333          *  */
 334         void setPressureLossField(string fileName, string fieldName){
 335                 _pressureLossField=Field(fileName, FACES,fieldName);
 336                 _pressureLossFieldSet=true;
 337         }
 338
 339         /** \fn setSectionField
 340          * \brief set the cross section field;
 341          * @param [in] Field cross section field (field on CELLS)
 342          * */
 343         void setSectionField(Field sectionField){
 344                 _sectionField=sectionField;
 345                 _sectionFieldSet=true;
 346         }
 347         /** \fn setSectionField
 348          * \brief set the cross section field
 349          * \details for variable cross section pipe network
 350          * \param [in] string fileName (including file path)
 351          * \param [in] string fieldName
 352          * \param [out] void
 353          *  */
 354         void setSectionField(string fileName, string fieldName){
 355                 _sectionField=Field(fileName, CELLS,fieldName);
 356                 _sectionFieldSet=true;
 357         }
 358
 359         /** \fn setNonLinearFormulation
 360          * \brief sets the formulation used for the computation of non viscous fluxes
 361          * \details Roe, VFRoe, VFFC
 362          * \param [in] enum NonLinearFormulation
 363          * \param [out] void
 364          *  */
 365         void setNonLinearFormulation(NonLinearFormulation nonLinearFormulation){
 366                 //if(nonLinearFormulation != Roe && nonLinearFormulation != VFRoe && nonLinearFormulation != VFFC && nonLinearFormulation!=reducedRoe)
 367                 //      throw CdmathException("nonLinearFormulation should be either Roe, VFRoe or VFFC");//extra security for swig binding
 368                 _nonLinearFormulation=nonLinearFormulation;
 369         }
 370
 371         /** \fn getNonLinearFormulation
 372          * \brief returns the formulation used for the computation of non viscous fluxes
 373          * \details Roe, VFRoe, VFFC
 374          * \param [in] void
 375          * \param [out] enum NonLinearFormulation
 376          *  */
 377         NonLinearFormulation getNonLinearFormulation(){
 378                 return _nonLinearFormulation;
 379         }
 380
 381         /** \fn usePrimitiveVarsInNewton
 382          * \brief use Primitive Vars instead of conservative vars In Newton scheme for implicit schemes
 383          * \param [in] bool
 384          *  */
 385         void usePrimitiveVarsInNewton(bool usePrimitiveVarsInNewton){
 386                 _usePrimitiveVarsInNewton=usePrimitiveVarsInNewton;
 387         }
 388
 389         //données initiales
 390         /*
 391         virtual vector<string> getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales)
 392         virtual  Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc)
 393         virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale
 394         virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
 395         virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement
 396          */
 397
 398 protected :
 399         /** Number of phases in the fluid **/
 400         int _nbPhases;
 401         /** Field of conservative variables (the primitive variables are defined in the mother class ProblemCoreFlows **/
 402         Field  _UU;
 403         /** Field of interfacial states of the VFRoe scheme **/
 404         Field _UUstar, _VVstar;
 405         /** the formulation used to compute the non viscous fluxes **/
 406         NonLinearFormulation _nonLinearFormulation;
 407
 408         /** boolean used to specify that an entropic correction should be used **/
 409         bool _entropicCorrection;
 410         /** Vector containing the eigenvalue jumps for the entropic correction **/
 411         vector<double> _entropicShift;
 412         /** In case a pressure correction scheme is used some more option regarding the type of pressure correction **/
 413         int _pressureCorrectionOrder;
 414
 415         /** Fluid equation of state **/
 416         vector< Fluide* > _fluides;
 417         //!Viscosity coefficients
 418         vector<double> _viscosite;
 419         //!Conductivity coefficients
 420         vector<double> _conductivite;
 421
 422         /** Source terms **/
 423         vector<double> _gravite, _GravityField3d, _gravityReferencePoint, _dragCoeffs;//_GravityField3d has size _Ndim whereas _gravite has size _Nvar and is usefull for dealing with source term and implicitation of gravity vector
 424         double _latentHeat, _Tsat,_Psat,_dHsatl_over_dp;
 425         Field _porosityField, _pressureLossField, _dp_over_dt, _sectionField;
 426         bool _porosityFieldSet, _pressureLossFieldSet, _sectionFieldSet;
 427         double _porosityi, _porosityj;//porosity of the left and right states around an interface
 428
 429         /** User options **/
 430         bool _isScaling;
 431         bool _saveVelocity;/* In order to display streamlines with paraview */
 432         bool _saveConservativeField;/* Save conservative fields such as density or momentum for instance */
 433         bool _saveInterfacialField;/* Save interfacial fields of the VFRoe scheme */
 434         bool _usePrimitiveVarsInNewton;
 435
 436         // Variables du schema numerique
 437         Vec _conservativeVars, _newtonVariation, _bScaling,_old, _primitiveVars, _Uext,_Uextdiff ,_vecScaling,_invVecScaling, _Vext;
 438         //courant state vector, vector computed at next time step, second member of the equation
 439         PetscScalar *_AroePlus, *_AroeMinus,*_Jcb,*_JcbDiff, *_a, *_blockDiag,  *_invBlockDiag,*_Diffusion, *_GravityImplicitationMatrix;
 440         PetscScalar *_Aroe, *_absAroe, *_signAroe, *_invAroe;
 441         PetscScalar *_AroeMinusImplicit,*_AroePlusImplicit,*_AroeImplicit,*_absAroeImplicit;//negative part of the roe matrix used in the implicit scheme matrix
 442         PetscScalar * _primToConsJacoMat; //Jacobian matrix of the prim-->cons function
 443         PetscScalar *_phi, *_Ui, *_Uj,  *_Vi, *_Vj,  *_Si, *_Sj, * _pressureLossVector, * _porosityGradientSourceVector, *_externalStates;
 444         double *_Uroe, *_Udiff, *_temp, *_l, *_r,  *_vec_normal;
 445         double * _Uij, *_Vij;//Conservative and primitive interfacial states of the VFRoe scheme
 446         int *_idm, *_idn;
 447         double _inv_dxi,_inv_dxj;//diametre des cellules i et j autour d'une face
 448         double _err_press_max,_part_imag_max,_minm1,_minm2;
 449         int _nbMaillesNeg, _nbVpCplx;
 450         bool _isBoundary;// la face courante est elle une face de bord ?
 451         double _maxvploc;
 452
 453         /** \fn convectionState
 454          * \brief calcule l'etat de Roe entre deux cellules voisinnes
 455          * \Details c'ets une fonction virtuelle, on la surcharge dans chacun des modèles
 456          * @param i,j : entiers correspondant aux numéros des cellules à gauche et à droite de l'interface
 457          * @param IsBord : est un booléen qui nous dit si la cellule voisine est sur le bord ou pas
 458          * */
 459         virtual void convectionState(const long &i, const long &j, const bool &IsBord)=0;
 460
 461         /** \fn convectionMatrices
 462          * \brief calcule la matrice de convection de l'etat interfacial entre deux cellules voisinnes
 463          * \Details convectionMatrices est une fonction virtuelle pure, on la surcharge dans chacun des modèles
 464          * */
 465         virtual void convectionMatrices()=0;
 466
 467         /** \fn diffusionStateAndMatrices
 468          * \brief calcule la matrice de diffusion de l'etat interface pour la diffusion
 469          * \Details est une fonction virtuelle pure, on la surcharge dans chacun eds modèles
 470          * @param i,j : entier correspondent aux indices des cellules  à gauche et droite respectivement
 471          * @param IsBord: bollean telling if (i,j) is a boundary face
 472          * @return
 473          * */
 474         virtual void diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord)=0;
 475
 476         /** \fn addSourceTermToSecondMember
 477          * \brief Adds the contribution of source terms to the right hand side of the system: gravity,
 478          *  phase change, heat power and friction
 479          * @param i,j : left and right cell number
 480          * @param nbNeighboursi, integer giving the number of neighbours of cell i
 481          * @param nbNeighboursj, integer giving the number of neighbours of cell j
 482          * @param boolean isBoundary is true for a boundary face (i,j) and false otherwise
 483          * @param double mesureFace the lenght or area of the face
 484          * */
 485         void addSourceTermToSecondMember(const int i, int nbNeighboursi,const int j, int nbNeighboursj,bool isBoundary, int ij, double mesureFace);
 486
 487         /** \fn sourceVector
 488          * \brief Computes the source term (at the exclusion of pressure loss terms)
 489          * \Details pure virtual function, overloaded by each model
 490          * @param Si output vector containing the source  term
 491          * @param Ui, Vi input conservative and primitive vectors
 492          * @param i the cell number. Used to read the power field
 493          * @return
 494          * */
 495         virtual void sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i)=0;
 496
 497         /** \fn convectionFlux
 498          * \brief Computes the convection flux F(U) projected on a vector n
 499          * @param U is the conservative variable vector
 500          * @param V is the primitive variable vector
 501          * @param normal is a unit vector giving the direction where the convection flux matrix F(U) is projected
 502          * @param porosity is the ration of the volume occupied by the fluid in the cell (default value is 1)
 503          * @return The convection flux projected in the direction given by the normal vector: F(U)*normal */
 504         virtual Vector convectionFlux(Vector U,Vector V, Vector normale, double porosity)=0;
 505
 506         /** \fn pressureLossVector
 507          * \brief Computes the contribution of pressure loss terms in the source term
 508          * \Details pure virtual function, overloaded by each model
 509          * @param pressureLoss output vector containing the pressure loss contributions
 510          * @param K, input pressure loss coefficient
 511          * @param Ui input primitive vectors
 512          * @param Vi input conservative vectors
 513          * @param Uj input primitive vectors
 514          * @param Vj input conservative vectors
 515          * @return
 516          * */
 517         virtual void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)=0;
 518
 519         /** \fn porosityGradientSourceVector
 520          * \brief Computes the contribution of the porosity variation in the source term
 521          * \Details pure virtual function, overloaded by each model
 522          * @param porosityGradientVector output vector containing the porosity variation contributions to the source term
 523          * @param Ui input primitive vectors celli
 524          * @param Vi input conservative vectors celli
 525          * @param porosityi input porosity value celli
 526          * @param deltaxi input diameter celli
 527          * @param Uj input primitive vectors cellj
 528          * @param Vj input conservative vectors cellj
 529          * @param porosityj input porosity value cellj
 530          * @param deltaxj input diameter cellj
 531          * @return
 532          * */
 533         virtual void porosityGradientSourceVector()=0;
 534
 535         //matrice de gravite
 536
 537         /** \fn jacobian
 538          * \brief Calcule la jacobienne de la ConditionLimite convection
 539          * \Details est une fonction virtuelle pure, qu'on surcharge dans chacun des modèles
 540          * @param j entier , l'indice de la cellule sur le bord
 541          * @param nameOfGroup  : chaine de caractères, correspond au type de la condition limite
 542          * */
 543         virtual void jacobian(const int &j, string nameOfGroup,double * normale)=0;
 544
 545         /** \fn jacobianDiff
 546          * \brief Calcule la jacobienne de la CL de diffusion
 547          * \Details est une fonction virtuelle pure, qu'on surcharge dans chacun des modèles
 548          * @param j entier , l'indice de la cellule sur le bord
 549          * @param nameOfGroup  : chaine de caractères, correspond au type de la condition limite
 550          * */
 551         virtual void jacobianDiff(const int &j, string nameOfGroup)=0;
 552
 553         /** \fn setBoundaryState
 554          * \brief Calcule l'etat fictif a la frontiere
 555          * \Details est une fonction virtuelle pure, qu'on surcharge dans chacun des modèles
 556          * @param j entier , l'indice de la cellule sur le bord
 557          * @param nameOfGroup  : chaine de caractères, correspond au type de la condition limite
 558          * @param normale est un vecteur de double correspond au vecteur normal à l'interface
 559          * @return
 560          * */
 561         virtual void setBoundaryState(string nameOfGroup, const int &j,double *normale)=0;
 562
 563         /** \fn addDiffusionToSecondMember
 564          * \brief Compute the contribution of the diffusion operator to the right hand side of the system
 565          * \Details this function is pure virtual, and overloaded in each physical model class
 566          * @param i left cell number
 567          * @param j right cell number
 568          * @param boolean isBoundary is true for a boundary face (i,j) and false otherwise
 569          * */
 570         virtual void addDiffusionToSecondMember(const int &i,const int &j,bool isBoundary)=0;
 571
 572         //!Computes the interfacial flux for the VFFC formulation of the staggered upwinding
 573         virtual Vector staggeredVFFCFlux()=0;
 574         //!Compute the corrected interfacial state for lowMach, pressureCorrection and staggered versions of the VFRoe formulation
 575         virtual void applyVFRoeLowMachCorrections(bool isBord, string groupname="")=0;
 576
 577         //remplit les vecteurs de scaling
 578         /** \fn computeScaling
 579          * \brief Special preconditioner based on a matrix scaling strategy
 580          * \Details pure  virtual function overloaded in every model class
 581          * @param offset double , correspond à la plus grande valeur propre
 582          * @return
 583          * */
 584         virtual void computeScaling(double offset) =0;
 585
 586         // Fonctions utilisant la loi d'etat
 587
 588         /** \fn consToPrim
 589          * \brief computes the primitive vector state from a conservative vector state
 590          * \Details ure virtual, implemented by each model
 591          * @param Ucons : conservative variable vector
 592          * @pram Vprim : primitive variable vector
 593          * @param porosity is the porosity coefficient in case of a porous modeling
 594          * */
 595         virtual void consToPrim(const double *Ucons, double* Vprim,double porosity=1) = 0;
 596
 597         /** \fn primToCons
 598          * \brief computes the conservative vector state from a primitive vector state
 599          * \Details pure virtual, implemented by each model
 600          * @param U : conservative variable vector, may contain several states
 601          * @pram V : primitive variable vector, may contain several states
 602          * @param i : index of the conservative state in the vector U
 603          * @param j : index of the primitive state in the vector V
 604          *       */
 605         virtual void primToCons(const double *V, const int &i, double *U, const int &j)=0;
 606
 607         /** \fn primToConsJacobianMatrix
 608          * \brief computes the jacobian matrix of the cons->prim function
 609          * \Details pure virtual, implemented by each model
 610          * @pram V : primitive vector state
 611          *       */
 612         //void primToConsJacobianMatrix(double *V)=0;
 613
 614         /** \fn getRacines
 615          * \brief Computes the roots of a polynomial
 616          * @param polynome is a vector containing the coefficients of the polynom
 617          * @return vector containing the roots (complex numbers)
 618          * */
 619         vector< complex<double> > getRacines(vector< double > polynome);
 620
 621         // Some supplement functions ---------------------------------------------------------------------------------------------
 622
 623         /** \fn addConvectionToSecondMember
 624          * \brief Adds the contribution of the convection to the system right hand side for a face (i,j) inside the domain
 625          * @param i left cell number
 626          * @param j right cell number
 627          * @param isBord is a boolean that is true if the interface (i,j) is a boundary interface
 628          * @param groupname : is a string that may be used when isBord is true to specify which boundary the face (i,j) belongs to
 629          * */
 630         virtual void addConvectionToSecondMember(const int &i,const int &j,bool isBord, string groupname="");
 631
 632         /** \fn updatePrimitives
 633          * \brief updates the global primitive vector from the global conservative vector
 634          * @param void
 635          * */
 636         void updatePrimitives();
 637
 638         /** \fn updateConservatives
 639          * \brief updates the global conservative vector from the global primitive vector
 640          * @param void
 641          * */
 642         void updateConservatives();
 643
 644         /** \fn AbsMatriceRoe
 645          * \brief Computes the absolute value of the Roe matrix
 646          * @param valeurs_propres_dist is the vector of distinct eigenvalues of the Roe matrix
 647          * */
 648         void AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist);
 649
 650         /** \fn SigneMatriceRoe
 651          * \brief Computes the sign of the Roe matrix
 652          * @param valeurs_propres_dist is the vector of distinct eigenvalues of the Roe matrix
 653          * */
 654         void SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist);
 655
 656         /** \fn InvMatriceRoe
 657          * \brief Computes the inverse of the Roe matrix
 658          * @param valeurs_propres_dist is the vector of distinct eigenvalues of the Roe matrix
 659          * */
 660         void InvMatriceRoe(vector< complex<double> > valeurs_propres_dist);
 661
 662         /** \fn entropicShift
 663          * \brief computes the eigenvalue jumps for the entropy correction
 664          * @param normal vector n to the interface between the two cells _l and _r
 665          * */
 666         virtual void entropicShift(double* n)=0;
 667
 668 };
 669
 670 #endif /* PROBLEMFLUID_HXX_ */