From: michael <michael@localhost.localdomain>
Date: Thu, 13 Jan 2022 17:24:37 +0000 (+0100)
Subject: Updated boundary conditions
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=1d3df271aa134bef0aff333d9ade38755996a01a;p=tools%2Fsolverlab.git

Updated boundary conditions
---

diff --git a/CoreFlows/Models/inc/TransportEquation.hxx b/CoreFlows/Models/inc/TransportEquation.hxx
index 125968f..d1dd186 100755
--- a/CoreFlows/Models/inc/TransportEquation.hxx
+++ b/CoreFlows/Models/inc/TransportEquation.hxx
@@ -20,26 +20,18 @@
 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{
@@ -61,10 +53,10 @@ public :
 	/** \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();
@@ -77,25 +69,74 @@ public :
 	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
diff --git a/CoreFlows/Models/src/TransportEquation.cxx b/CoreFlows/Models/src/TransportEquation.cxx
index aa2a9dc..bd15ace 100755
--- a/CoreFlows/Models/src/TransportEquation.cxx
+++ b/CoreFlows/Models/src/TransportEquation.cxx
@@ -5,39 +5,54 @@
 
 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
@@ -124,14 +139,8 @@ void TransportEquation::initialize()
 		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
 }
@@ -649,3 +658,61 @@ TransportEquation::setRodTemperatureField(Field rodTemperature){
 	_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 
+}