using namespace std;
-//! enumeration phase
-/*! The fluid type can be LiquidPhase or water */
-enum phase
-{
- LiquidPhase,/**< Fluid considered is GasPhase */
- GasPhase/**< Fluid considered is Gas */
-};
-
-//! enumeration pressureEstimate
-/*! the pressure estimate needed to fit physical parameters */
-enum pressureMagnitude
-{
- around1bar300KTransport,/**< pressure is around 1 bar and temperature around 300K (for TransportEquation, SinglePhase and IsothermalTwoFluid) or 373 K (saturation for DriftModel and FiveEqsTwoFluid) */
- around155bars600KTransport/**< pressure is around 155 bars and temperature around 618 K (saturation) */
-};
-
//! enumeration BoundaryType
/*! Boundary condition type */
enum BoundaryTypeTransport {InletTransport, OutletTransport, NeumannTransport, DirichletTransport, NoneBCTransport};//Actually Inlet=Dirichlet and Outlet=Neumann
+//! enumeration phaseType
+/*! The material phase can be Solid, Gas or liquid */
+enum FluidMaterial
+{
+ Air,/**< Material considered is air */
+ Water/**< Material considered is water */
+};
+
/** \struct LimitField
* \brief value of some fields on the boundary */
struct LimitFieldTransport{
/** \fn TransportEquation
* \brief Constructor for the enthalpy transport in a fluid
* \param [in] phase : \ref Liquid or \ref Gas
- * \param [in] pressureMagnitude : \ref around1bar or \ref around155bars
+ * \param [in] pressureEstimate : \ref around1bar or \ref around155bars
* \param [in] vector<double> : fluid velocity (assumed constant)
* */
- TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport, MPI_Comm comm = MPI_COMM_WORLD);
+ TransportEquation(FluidMaterial fluid, pressureEstimate pEstimate,vector<double> vitesseTransport, MPI_Comm comm = MPI_COMM_WORLD);
//Gestion du calcul
virtual void initialize();
virtual void validateTimeStep();
/* Boundary conditions */
- /** \fn setIntletBoundaryCondition
+ /** \fn setIntletBoundaryCondition
* \brief adds a new boundary condition of type Inlet
- * \details
+ * \details same as setDirichletBoundaryCondition
* \param [in] string : the name of the boundary
- * \param [in] double : the value of the temperature at the boundary
+ * \param [in] double : the value of the enthalpy at the boundary
* \param [out] void
* */
void setInletBoundaryCondition(string groupName,double enthalpy){
_limitField[groupName]=LimitFieldTransport(InletTransport,-1,enthalpy,-1);
};
+ /** \fn setDirichletBoundaryCondition
+ * \brief adds a new boundary condition of type Dirichlet
+ * \details same as setInletBoundaryCondition
+ * \param [in] string : the name of the boundary
+ * \param [in] double : the value of the enthalpy at the boundary
+ * \param [out] void
+ * */
+ void setDirichletBoundaryCondition(string groupName,double enthalpy){
+ _limitField[groupName]=LimitFieldTransport(DirichletTransport,-1,enthalpy,-1);
+ };
+ /** \fn setDirichletBoundaryCondition
+ * \brief adds a new boundary condition of type Inlet
+ * \details Reads the boundary field in a med file
+ * \param [in] string : the name of the boundary
+ * \param [in] string : the file name
+ * \param [in] string : the field name
+ * \param [in] int : the time step number
+ * \param [in] int : int corresponding to the enum CELLS or NODES
+ * \param [out] void
+ * */
+ void setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
+ void setDirichletBoundaryCondition(string groupName, Field bc_field){
+ _limitField[groupName]=LimitFieldTransport(DirichletTransport, -1, 0, -1);
+ };
- /** \fn setNeumannBoundaryCondition
+ /** \fn setNeumannBoundaryCondition
* \brief adds a new boundary condition of type Neumann
- * \details
+ * \details same as setOutletBoundaryCondition
* \param [in] string the name of the boundary
+ * \param [in] double : the value of the enthalpy flux at the boundary
* \param [out] void
* */
void setNeumannBoundaryCondition(string groupName, double flux=0){
- _limitField[groupName]=LimitFieldTransport(NeumannTransport,-1,flux,-1);
+ _limitField[groupName]=LimitFieldTransport(NeumannTransport,-1,-1,flux);
+ };
+ /** \fn setOutletBoundaryCondition
+ * \brief adds a new boundary condition of type Outlet
+ * \details same as setNeumannBoundaryCondition
+ * \param [in] string the name of the boundary
+ * \param [in] double : the value of the enthalpy flux at the boundary
+ * \param [out] void
+ * */
+ void setOutletBoundaryCondition(string groupName, double flux=0){
+ _limitField[groupName]=LimitFieldTransport(OutletTransport,-1,-1,flux);
+ };
+ /** \fn setNeumannBoundaryCondition
+ * \brief adds a new boundary condition of type Neumann
+ * \details Reads the boundary field in a med file
+ * \param [in] string : the name of the boundary
+ * \param [in] string : the file name
+ * \param [in] string : the field name
+ * \param [in] int : the time step number
+ * \param [in] int : int corresponding to the enum CELLS or NODES
+ * \param [out] void
+ * */
+ void setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
+ void setNeumannBoundaryCondition(string groupName, Field bc_field){
+ _limitField[groupName]=LimitFieldTransport(NeumannTransport,-1,-1, 0);
};
/** \fn setBoundaryFields
using namespace std;
-TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport, MPI_Comm comm):ProblemCoreFlows(comm)
+TransportEquation::TransportEquation(FluidMaterial fluid, pressureEstimate pEstimate,vector<double> vitesseTransport, MPI_Comm comm):ProblemCoreFlows(comm)
{
- if(pEstimate==around1bar300KTransport){
+ if(pEstimate==around1bar300K){
_Tref=300;
- if(fluid==GasPhase){//Nitrogen pressure 1 bar and temperature 27°C
- _href=3.11e5; //nitrogen enthalpy at 1 bar and 300K
- _cpref=1041;//nitrogen specific heat at constant pressure 1 bar and 300K
- //saturation data for nitrogen at 1 bar and 77K
- _hsatv=0.77e5;//nitrogen vapour enthalpy at saturation at 1 bar
- _hsatl=-1.22e5;//nitrogen liquid enthalpy at saturation at 1 bar
- _rhosatv=4.556;//nitrogen vapour density at saturation at 1 bar
- _rhosatl=806.6;//nitrogen liquid density at saturation at 1 bar
- }
- else{
- //Water at pressure 1 bar and temperature 27°C
- _href=1.127e5; //water enthalpy at 1 bar and 300K
- _cpref=4181;//water specific heat at 1 bar and 300K
- //saturation data for water at 1 bar and 373K
- _hsatv=2.675e6;//Gas enthalpy at saturation at 1 bar
- _hsatl=4.175e5;//water enthalpy at saturation at 1 bar
- _rhosatv=0.6;//Gas density at saturation at 1 bar
- _rhosatl=958;//water density at saturation at 1 bar
+ switch (fluid)
+ {
+ case Air://Nitrogen pressure 1 bar and temperature 27°C
+ _href=3.11e5; //nitrogen enthalpy at 1 bar and 300K
+ _cpref=1041;//nitrogen specific heat at constant pressure 1 bar and 300K
+ //saturation data for nitrogen at 1 bar and 77K
+ _hsatv=0.77e5;//nitrogen vapour enthalpy at saturation at 1 bar
+ _hsatl=-1.22e5;//nitrogen liquid enthalpy at saturation at 1 bar
+ _rhosatv=4.556;//nitrogen vapour density at saturation at 1 bar
+ _rhosatl=806.6;//nitrogen liquid density at saturation at 1 bar
+ cout<<"Air at around 1 bar and 300 Kelvin"<<endl;
+ break;
+ case Water://Water at pressure 1 bar and temperature 27°C
+ _href=1.127e5; //water enthalpy at 1 bar and 300K
+ _cpref=4181;//water specific heat at 1 bar and 300K
+ //saturation data for water at 1 bar and 373K
+ _hsatv=2.675e6;//Gas enthalpy at saturation at 1 bar
+ _hsatl=4.175e5;//water enthalpy at saturation at 1 bar
+ _rhosatv=0.6;//Gas density at saturation at 1 bar
+ _rhosatl=958;//water density at saturation at 1 bar
+ cout<<"Water at around 1 bar and 300 Kelvin"<<endl;
+ break;
+ default://Solid phase
+ cout<<"Error TransportEquation::TransportEquation : Only Gas and Liquid phases accepted"<<endl;
+ throw CdmathException("TransportEquation::TransportEquation : Wrong material phase requested");
}
}
else{//around155bars600K
_Tref=618;//=Tsat
- if(fluid==GasPhase){
- _href=2.675e6; //Gas enthalpy at 155 bars and 618K
- _cpref=14001;//Gas specific heat at 155 bar and 618K
- }
- else{//Liquid
- _href=4.175e5;//water enthalpy at 155 bars and 618K
- _cpref=8950;//water specific heat at 155 bar and 618K
+ switch (fluid)
+ {
+ case Air://Nitrogen pressure 155 bars and temperature 327°C
+ _href=2.675e6; //Gas enthalpy at 155 bars and 618K
+ _cpref=14001;//Gas specific heat at 155 bar and 618K
+ cout<<"Air at around 155 bar and 600 Kelvin"<<endl;
+ break;
+ case Water://Water at pressure 155 bars and temperature 327°C
+ _href=4.175e5;//water enthalpy at 155 bars and 618K
+ _cpref=8950;//water specific heat at 155 bar and 618K
+ cout<<"Water at around 155 bar and 600 Kelvin"<<endl;
+ break;
+ default://Solid phase
+ cout<<"Error TransportEquation::TransportEquation : Only Gas and Liquid phases accepted"<<endl;
+ throw CdmathException("TransportEquation::TransportEquation : Wrong material phase requested");
}
//saturation data for water at 155 bars and 618K
_hsatv=2.6e6;//Gas enthalpy at saturation at 155 bars
VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Hn_seq);//For saving results on proc 0
VecScatterCreateToZero(_Hn,&_scat,&_Hn_seq);
- //Linear solver
- KSPCreate(PETSC_COMM_SELF, &_ksp);
- KSPSetType(_ksp, _ksptype);
- // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
- KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
- KSPGetPC(_ksp, &_pc);
- PCSetType(_pc, _pctype);
-
+ createKSP();
+
_initializedMemory=true;
save();//save initial data
}
_rodTemperatureFieldSet=true;
_isStationary=false;//Source term may be changed after previously reaching a stationary state
}
+
+void
+TransportEquation::setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
+ if(_FECalculation && field_support_type != NODES)
+ cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+ else if(!_FECalculation && field_support_type == NODES)
+ cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+
+ Field VV;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case NODES:
+ VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case FACES:
+ VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error TransportEquation::setDirichletBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+ /* For the moment the boundary value is taken constant equal to zero */
+ _limitField[groupName]=LimitFieldTransport(DirichletTransport,-1, 0,-1);//This line will be deleted when variable BC are properly treated in solverlab
+}
+
+void
+TransportEquation::setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
+ if(_FECalculation && field_support_type != NODES)
+ cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+ else if(!_FECalculation && field_support_type == NODES)
+ cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+
+ Field VV;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case NODES:
+ VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case FACES:
+ VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error TransportEquation::setNeumannBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+ /* For the moment the boundary value is taken constant equal to zero */
+ _limitField[groupName]=LimitFieldTransport(NeumannTransport,-1, 0,-1);//This line will be deleted when variable BC are properly treated in solverlab
+}
