using namespace std;
+//! enumeration SpaceScheme
+/*! Several numerical schemes are available */
+enum SpaceScheme
+{
+ upwind,/**< classical full upwinding scheme (first order in space) */
+ centered,/**< centered scheme (second order in space) */
+ pressureCorrection,/**< include a pressure correction in the upwind scheme to increase precision at low Mach numbers */
+ lowMach,/**< include an upwinding proportional to the Mach numer scheme to increase precision at low Mach numbers */
+ staggered,/**< scheme inspired by staggered discretisations */
+};
+
//! enumeration NonLinearFormulation
/*! the formulation used to compute the non viscous fluxes */
enum NonLinearFormulation
reducedRoe,/**< compacted formulation of Roe scheme without computation of the fluxes */
};
+//! enumeration phaseType
+/*! The material phase can be Gas or liquid */
+enum phaseType
+{
+ Liquid,/**< Material considered is Liquid */
+ Gas/**< Material considered is Gas */
+};
+
+//! enumeration BoundaryType
+/*! Boundary condition type */
+enum BoundaryType {Wall, InnerWall, Inlet, InletPressure, InletRotationVelocity, InletEnthalpy, Outlet, Neumann, NoTypeSpecified};
+/** \struct LimitField
+ * \brief value of some fields on the boundary */
+struct LimitField{
+ 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;}
+ LimitField(BoundaryType _bcType, double _p, vector<double> _v_x, vector<double> _v_y, vector<double> _v_z,
+ double _T, double _h, double _alpha, double _conc){
+ bcType=_bcType; p=_p; v_x=_v_x; v_y=_v_y; v_z=_v_z; T=_T; h=_h; alpha=_alpha; conc=_conc;
+ }
+
+ BoundaryType bcType;
+ double p;//For outlet (fluid models)
+ vector<double> v_x; vector<double> v_y; vector<double> v_z;//For wall and inlet (fluid models)
+ double T; //for wall and inlet (DriftModel and FiveEqsTwoFluid) and for Dirichlet (DiffusionEquation)
+ double h; //for inlet (TransportEquation)
+ double alpha; //For inlet (IsothermalTwoFluid and FiveEqsTwoFluid)
+ double conc;//For inlet (DriftModel)
+};
+
class ProblemFluid: public ProblemCoreFlows
{
/**\fn
* \brief constructeur de la classe ProblemFluid
*/
- ProblemFluid(void);
+ ProblemFluid(MPI_Comm comm = MPI_COMM_WORLD);
//Gestion du calcul (interface ICoCo)
* */
virtual void validateTimeStep();
+ /** \fn solveTimeStep
+ * \brief calcule les valeurs inconnues au pas de temps +1 .
+ * \details c'est une fonction virtuelle, qui surcharge celle de la classe ProblemCoreFlows
+ * @param void
+ * \return Renvoie false en cas de problème durant le calcul (valeurs non physiques..)
+ * */
+ virtual bool solveTimeStep();//
+
+ /* Boundary conditions */
+ /** \fn setNeumannBoundaryCondition
+ * \brief adds a new boundary condition of type Neumann
+ * \details
+ * \param [in] string the name of the boundary
+ * \param [out] void
+ * */
+ void setNeumannBoundaryCondition(string groupName){
+ _limitField[groupName]=LimitField(Neumann,-1,vector<double>(_Ndim,0),vector<double>(_Ndim,0),vector<double>(_Ndim,0),-1,-1,-1,-1);
+ };
+
/** \fn setOutletBoundaryCondition
* \brief Adds a new boundary condition of type Outlet
* \details
_limitField[groupName]=LimitField(Outlet,referencePressure,vector<double>(_nbPhases,0),vector<double>(_nbPhases,0),vector<double>(_nbPhases,0),-1,-1,-1,-1);
};
+ /** \fn setBoundaryFields
+ * \brief met à jour _limitField ( le type de condition limite )
+ * \details
+ * \param [in] string
+ * \param [out] void
+ * */
+ void setBoundaryFields(map<string, LimitField> boundaryFields){
+ _limitField = boundaryFields;
+ };
+
/** \fn setViscosity
* \brief sets the vector of viscosity coefficients
* @param viscosite is a vector of size equal to the number of phases and containing the viscosity of each phase
* @param void
* @return the number of phases considered in the model
* */
- int getNumberOfPhases(){
+ int getNumberOfPhases() const {
return _nbPhases;
};
- /** \fn computeNewtonVariation
- * \brief Builds and solves the linear system to obtain the variation Ukp1-Uk in a Newton scheme
- * @param void
- * */
- virtual void computeNewtonVariation();
-
/** \fn testConservation
* \brief Teste et affiche la conservation de masse et de la quantité de mouvement
* \Details la fonction est virtuelle pure, on la surcharge dans chacun des modèles
* @param pcType preconditioner (ILU,LU or NONE)
* @param scaling performs a bancing of the system matrix before calling th preconditioner
*/
- void setLinearSolver(linearSolver kspType, preconditioner pcType, bool scaling=false)
+ void setLinearSolver(linearSolver kspType, preconditioner pcType, double maxIts=50, bool scaling=false)
{
- ProblemCoreFlows::setLinearSolver(kspType, pcType);
+ ProblemCoreFlows::setLinearSolver(kspType, pcType, maxIts);
_isScaling= scaling;
};
* \brief returns the value of the latent heat
* @param double L, the value of the latent heat
* */
- double getLatentHeat(){
+ double getLatentHeat() const{
return _latentHeat;
}
* \brief sets the saturation temperature
* @param Tsat double corresponds to saturation temperature
* */
- double getSatTemp(){
+ double getSatTemp() const {
return _Tsat;
}
* \brief set the porosity field;
* @param [in] Field porosity field (field on CELLS)
* */
- void setPorosityField(Field Porosity){
- _porosityField=Porosity;
- _porosityFieldSet=true;
- }
+ void setPorosityField(Field Porosity);
/** \fn getPorosityField
* \brief returns the porosity field;
* @param
* */
- Field getPorosityField(){
+ Field getPorosityField() const {
return _porosityField;
}
* \param [in] string fieldName
* \param [out] void
* */
- void setPorosityField(string fileName, string fieldName){
- _porosityField=Field(fileName, CELLS,fieldName);
- _porosityFieldSet=true;
- }
+ void setPorosityField(string fileName, string fieldName);
/** \fn setPressureLossField
* \brief set the pressure loss coefficients field;
* @param [in] Field pressure loss field (field on FACES)
* */
- void setPressureLossField(Field PressureLoss){
- _pressureLossField=PressureLoss;
- _pressureLossFieldSet=true;
- }
+ void setPressureLossField(Field PressureLoss);
+
/** \fn setPressureLossField
* \brief set the pressure loss coefficient field
* \details localised friction force
* \param [in] string fieldName
* \param [out] void
* */
- void setPressureLossField(string fileName, string fieldName){
- _pressureLossField=Field(fileName, FACES,fieldName);
- _pressureLossFieldSet=true;
- }
+ void setPressureLossField(string fileName, string fieldName);
/** \fn setSectionField
* \brief set the cross section field;
* @param [in] Field cross section field (field on CELLS)
* */
- void setSectionField(Field sectionField){
- _sectionField=sectionField;
- _sectionFieldSet=true;
- }
+ void setSectionField(Field sectionField);
+
/** \fn setSectionField
* \brief set the cross section field
* \details for variable cross section pipe network
* \param [in] string fieldName
* \param [out] void
* */
- void setSectionField(string fileName, string fieldName){
- _sectionField=Field(fileName, CELLS,fieldName);
- _sectionFieldSet=true;
- }
+ void setSectionField(string fileName, string fieldName);
/** \fn setNonLinearFormulation
* \brief sets the formulation used for the computation of non viscous fluxes
* \param [in] void
* \param [out] enum NonLinearFormulation
* */
- NonLinearFormulation getNonLinearFormulation(){
+ NonLinearFormulation getNonLinearFormulation() const{
return _nonLinearFormulation;
}
_usePrimitiveVarsInNewton=usePrimitiveVarsInNewton;
}
- //données initiales
+ /** \fn getSpaceScheme
+ * \brief returns the space scheme name
+ * \param [in] void
+ * \param [out] enum SpaceScheme(upwind, centred, pressureCorrection, pressureCorrection, staggered)
+ * */
+ SpaceScheme getSpaceScheme() const;
+
+ /** \fn setNumericalScheme
+ * \brief sets the numerical method (upwind vs centered and explicit vs implicit)
+ * \details
+ * \param [in] SpaceScheme
+ * \param [in] TimeScheme
+ * \param [out] void
+ * */
+ void setNumericalScheme(SpaceScheme scheme, TimeScheme method=Explicit);
+
+ /* ICoCo code coupling interface */
/*
- virtual vector<string> getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales)
virtual Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc)
- virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale
virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement
*/
+ /** Set input fields to prepare the simulation or coupling **/
+ vector<string> getInputFieldsNames();
+ void setInputField(const string& nameField, Field& inputField );//supply of a required input field
+
+ Field getConservativeField() const ;
+ Field getPrimitiveField() const;
protected :
/** Number of phases in the fluid **/
Field _UU;
/** Field of interfacial states of the VFRoe scheme **/
Field _UUstar, _VVstar;
+
+ SpaceScheme _spaceScheme;
/** the formulation used to compute the non viscous fluxes **/
NonLinearFormulation _nonLinearFormulation;
+ /** PETSc nonlinear solver and line search **/
+ SNES _snes;
+ SNESLineSearch _linesearch;
+ PetscViewer _monitorLineSearch;
+
+ map<string, LimitField> _limitField;
+
/** boolean used to specify that an entropic correction should be used **/
bool _entropicCorrection;
/** Vector containing the eigenvalue jumps for the entropic correction **/
bool _isBoundary;// la face courante est elle une face de bord ?
double _maxvploc;
+ bool solveNewtonPETSc();//Use PETSc Newton methods to solve time step
+
+ /** \fn computeNewtonVariation
+ * \brief Builds and solves the linear system to obtain the variation Ukp1-Uk in a Newton scheme
+ * @param void
+ * */
+ virtual void computeNewtonVariation();
+
+ /** \fn computeNewtonRHS
+ * \brief Builds the right hand side F_X(X) of the linear system in the Newton method to obtain the variation Ukp1-Uk
+ * @param void
+ * */
+ void computeNewtonRHS( Vec X, Vec F_X);
+
+ /** \fn computeSnesRHS
+ * \brief Static function calling computeNewtonRHS to match PETSc nonlinear solver (SNES) structure
+ * @param void
+ * */
+ static int computeSnesRHS(SNES snes, Vec X, Vec F_X, void *ctx);
+
+ /** \fn computeNewtonJacobian
+ * \brief Static function calling computeNewtonJacobian to match PETSc nonlinear solver (SNES) structure
+ * @param void
+ * */
+ void computeNewtonJacobian( Vec X, Mat A);
+
+ /** \fn computeSnesJacobian
+ * \brief Builds the matrix A(X) of the linear system in the Newton method to obtain the variation Ukp1-Uk
+ * @param void
+ * */
+ static int computeSnesJacobian(SNES snes, Vec X, Mat A, Mat Aapprox, void *ctx);
+
/** \fn convectionState
* \brief calcule l'etat de Roe entre deux cellules voisinnes
* \Details c'ets une fonction virtuelle, on la surcharge dans chacun des modèles
* */
virtual void porosityGradientSourceVector()=0;
- //matrice de gravite
-
/** \fn jacobian
* \brief Calcule la jacobienne de la ConditionLimite convection
* \Details est une fonction virtuelle pure, qu'on surcharge dans chacun des modèles
