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