1 //============================================================================
2 // Name : ProblemCoreFlows
5 // Copyright : CEA Saclay 2014
6 // Description : Generic class for thermal hydraulics problems
7 //============================================================================
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
13 #ifndef PROBLEMFLUID_HXX_
14 #define PROBLEMFLUID_HXX_
17 #include "ProblemCoreFlows.hxx"
18 #include "utilitaire_algebre.h"
22 //! enumeration SpaceScheme
23 /*! Several numerical schemes are available */
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 */
33 //! enumeration pressureEstimate
34 /*! the pressure estimate needed to fit physical parameters */
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) */
41 //! enumeration phaseType
42 /*! The fluid type can be Gas or water */
45 Liquid,/**< Fluid considered is water */
46 Gas/**< Fluid considered is Gas */
49 //! enumeration NonLinearFormulation
50 /*! the formulation used to compute the non viscous fluxes */
51 enum NonLinearFormulation
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 */
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 */
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;
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)
80 class ProblemFluid: public ProblemCoreFlows
85 * \brief constructeur de la classe ProblemFluid
87 ProblemFluid(MPI_Comm comm = MPI_COMM_WORLD);
89 //Gestion du calcul (interface ICoCo)
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.
96 virtual void initialize();
99 * \brief empties the memory
102 virtual void terminate();
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
109 bool initTimeStep(double dt);
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)
116 double computeTimeStep(bool & stop);
118 /** \fn abortTimeStep
119 * \brief Reset the time step dt to 0
121 void abortTimeStep();
123 /** \fn iterateTimeStep
124 * \brief calls computeNewtonVariation to perform one Newton iteration and tests the convergence
126 * @return boolean ok is true is the newton iteration gave a physically acceptable result
128 virtual bool iterateTimeStep(bool &ok);
131 * \brief saves the current results in MED or VTK files
132 * \details It is a pure virtual function overloaded in each model
135 virtual void save()=0;
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
143 virtual void validateTimeStep();
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
149 * \return Renvoie false en cas de problème durant le calcul (valeurs non physiques..)
151 virtual bool solveTimeStep();//
153 /* Boundary conditions */
154 /** \fn setNeumannBoundaryCondition
155 * \brief adds a new boundary condition of type Neumann
157 * \param [in] string the name of the boundary
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);
164 /** \fn setOutletBoundaryCondition
165 * \brief Adds a new boundary condition of type Outlet
167 * \param [in] string : the name of the boundary
168 * \param [in] double : the value of the pressure at the boundary
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);
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
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);
189 /** \fn setBoundaryFields
190 * \brief met à jour _limitField ( le type de condition limite )
195 void setBoundaryFields(map<string, LimitField> boundaryFields){
196 _limitField = boundaryFields;
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
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]);
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
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]);
224 * \brief sets the gravity force in the model
225 * @param gravite is a vector of size equal to the space dimension
227 void setGravity(vector<double> gravite){
228 _GravityField3d = gravite;
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
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++)
241 _fluides[i]->setDragCoeffs(dragCoeffs[i]);
242 _dragCoeffs[i]=dragCoeffs[i];
246 /** \fn getNumberOfPhases
247 * \brief The numer of phase (one or two) depending on the model considered
249 * @return the number of phases considered in the model
251 int getNumberOfPhases() const {
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
260 virtual void testConservation()=0;
263 * \brief saves the velocity field in a separate 3D file so that paraview can display the streamlines
266 void saveVelocity(bool save_v=true){
267 _saveVelocity=save_v;
270 /** \fn saveConservativeField
271 * \brief saves the conservative fields (density, momentum etc...)
273 void saveConservativeField(bool save=true){
274 _saveConservativeField=save;
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
280 void setEntropicCorrection(bool entropyCorr){
281 _entropicCorrection=entropyCorr;
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)
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");
299 _pressureCorrectionOrder=order;
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
314 void setLinearSolver(linearSolver kspType, preconditioner pcType, double maxIts=50, bool scaling=false)
316 ProblemCoreFlows::setLinearSolver(kspType, pcType, maxIts);
320 /** \fn setLatentHeat
321 * \brief Sets the value of the latent heat
322 * @param double L, the value of the latent heat
324 void setLatentHeat(double L){
328 /** \fn getLatentHeat
329 * \brief returns the value of the latent heat
330 * @param double L, the value of the latent heat
332 double getLatentHeat() const{
337 * \brief sets the saturation temperature
338 * @param Tsat double corresponds to saturation temperature
340 void setSatTemp(double Tsat){
345 * \brief sets the saturation temperature
346 * @param Tsat double corresponds to saturation temperature
348 double getSatTemp() const {
352 /** \fn setSatPressure
353 * \brief sets the saturation pressure
354 * @param Psat double corresponds to saturation pressure
356 void setSatPressure(double Psat, double dHsatl_over_dp=0.05){
358 _dHsatl_over_dp=dHsatl_over_dp;
361 /** \fn setPorosityField
362 * \brief set the porosity field;
363 * @param [in] Field porosity field (field on CELLS)
365 void setPorosityField(Field Porosity){
366 Porosity.getMesh().checkFastEquivalWith(_mesh);
367 _porosityField=Porosity;
368 _porosityFieldSet=true;
371 /** \fn getPorosityField
372 * \brief returns the porosity field;
375 Field getPorosityField() const {
376 return _porosityField;
379 /** \fn setPorosityFieldFile
380 * \brief set the porosity field
382 * \param [in] string fileName (including file path)
383 * \param [in] string fieldName
386 void setPorosityField(string fileName, string fieldName){
387 _porosityField=Field(fileName, CELLS,fieldName);
388 _porosityField.getMesh().checkFastEquivalWith(_mesh);
389 _porosityFieldSet=true;
392 /** \fn setPressureLossField
393 * \brief set the pressure loss coefficients field;
394 * @param [in] Field pressure loss field (field on FACES)
396 void setPressureLossField(Field PressureLoss){
397 PressureLoss.getMesh().checkFastEquivalWith(_mesh);
398 _pressureLossField=PressureLoss;
399 _pressureLossFieldSet=true;
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
408 void setPressureLossField(string fileName, string fieldName){
409 _pressureLossField=Field(fileName, FACES,fieldName);
410 _pressureLossField.getMesh().checkFastEquivalWith(_mesh);
411 _pressureLossFieldSet=true;
414 /** \fn setSectionField
415 * \brief set the cross section field;
416 * @param [in] Field cross section field (field on CELLS)
418 void setSectionField(Field sectionField){
419 sectionField.getMesh().checkFastEquivalWith(_mesh);
420 _sectionField=sectionField;
421 _sectionFieldSet=true;
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
430 void setSectionField(string fileName, string fieldName){
431 _sectionField=Field(fileName, CELLS,fieldName);
432 _sectionField.getMesh().checkFastEquivalWith(_mesh);
433 _sectionFieldSet=true;
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
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;
448 /** \fn getNonLinearFormulation
449 * \brief returns the formulation used for the computation of non viscous fluxes
450 * \details Roe, VFRoe, VFFC
452 * \param [out] enum NonLinearFormulation
454 NonLinearFormulation getNonLinearFormulation() const{
455 return _nonLinearFormulation;
458 /** \fn usePrimitiveVarsInNewton
459 * \brief use Primitive Vars instead of conservative vars In Newton scheme for implicit schemes
462 void usePrimitiveVarsInNewton(bool usePrimitiveVarsInNewton){
463 _usePrimitiveVarsInNewton=usePrimitiveVarsInNewton;
466 /** \fn getSpaceScheme
467 * \brief returns the space scheme name
469 * \param [out] enum SpaceScheme(upwind, centred, pressureCorrection, pressureCorrection, staggered)
471 SpaceScheme getSpaceScheme() const;
473 /** \fn setNumericalScheme
474 * \brief sets the numerical method (upwind vs centered and explicit vs implicit)
476 * \param [in] SpaceScheme
477 * \param [in] TimeScheme
480 void setNumericalScheme(SpaceScheme scheme, TimeScheme method=Explicit);
482 /* ICoCo code coupling interface */
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
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
492 Field getConservativeField() const ;
493 Field getPrimitiveField() const;
496 /** Number of phases in the fluid **/
498 /** Field of conservative variables (the primitive variables are defined in the mother class ProblemCoreFlows **/
500 /** Field of interfacial states of the VFRoe scheme **/
501 Field _UUstar, _VVstar;
503 SpaceScheme _spaceScheme;
504 /** the formulation used to compute the non viscous fluxes **/
505 NonLinearFormulation _nonLinearFormulation;
507 /** PETSc nonlinear solver and line search **/
509 SNESLineSearch _linesearch;
510 PetscViewer _monitorLineSearch;
512 map<string, LimitField> _limitField;
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;
521 /** Fluid equation of state **/
522 vector< Fluide* > _fluides;
523 //!Viscosity coefficients
524 vector<double> _viscosite;
525 //!Conductivity coefficients
526 vector<double> _conductivite;
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
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;
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
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 ?
559 bool solveNewtonPETSc();//Use PETSc Newton methods to solve time step
561 /** \fn computeNewtonVariation
562 * \brief Builds and solves the linear system to obtain the variation Ukp1-Uk in a Newton scheme
565 virtual void computeNewtonVariation();
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
571 void computeNewtonRHS( Vec X, Vec F_X);
573 /** \fn computeSnesRHS
574 * \brief Static function calling computeNewtonRHS to match PETSc nonlinear solver (SNES) structure
577 static int computeSnesRHS(SNES snes, Vec X, Vec F_X, void *ctx);
579 /** \fn computeNewtonJacobian
580 * \brief Static function calling computeNewtonJacobian to match PETSc nonlinear solver (SNES) structure
583 void computeNewtonJacobian( Vec X, Mat A);
585 /** \fn computeSnesJacobian
586 * \brief Builds the matrix A(X) of the linear system in the Newton method to obtain the variation Ukp1-Uk
589 static int computeSnesJacobian(SNES snes, Vec X, Mat A, Mat Aapprox, void *ctx);
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
597 virtual void convectionState(const long &i, const long &j, const bool &IsBord)=0;
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
603 virtual void convectionMatrices()=0;
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
612 virtual void diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord)=0;
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
623 void addSourceTermToSecondMember(const int i, int nbNeighboursi,const int j, int nbNeighboursj,bool isBoundary, int ij, double mesureFace);
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
633 virtual void sourceVector(PetscScalar * Si,PetscScalar * Ui,PetscScalar * Vi, int i)=0;
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;
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
655 virtual void pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)=0;
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
671 virtual void porosityGradientSourceVector()=0;
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
679 virtual void jacobian(const int &j, string nameOfGroup,double * normale)=0;
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
687 virtual void jacobianDiff(const int &j, string nameOfGroup)=0;
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
697 virtual void setBoundaryState(string nameOfGroup, const int &j,double *normale)=0;
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
706 virtual void addDiffusionToSecondMember(const int &i,const int &j,bool isBoundary)=0;
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;
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
720 virtual void computeScaling(double offset) =0;
722 // Fonctions utilisant la loi d'etat
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
731 virtual void consToPrim(const double *Ucons, double* Vprim,double porosity=1) = 0;
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
741 virtual void primToCons(const double *V, const int &i, double *U, const int &j)=0;
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
748 //void primToConsJacobianMatrix(double *V)=0;
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)
755 vector< complex<double> > getRacines(vector< double > polynome);
757 // Some supplement functions ---------------------------------------------------------------------------------------------
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
766 virtual void addConvectionToSecondMember(const int &i,const int &j,bool isBord, string groupname="");
768 /** \fn updatePrimitives
769 * \brief updates the global primitive vector from the global conservative vector
772 void updatePrimitives();
774 /** \fn updateConservatives
775 * \brief updates the global conservative vector from the global primitive vector
778 void updateConservatives();
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
784 void AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist);
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
790 void SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist);
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
796 void InvMatriceRoe(vector< complex<double> > valeurs_propres_dist);
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
802 virtual void entropicShift(double* n)=0;
806 #endif /* PROBLEMFLUID_HXX_ */