1 //============================================================================
3 * \file DiffusionEquation.hxx
4 * \author Michael NDJINGA
7 * \brief Heat diffusion equation
8 * rho*cp*dT/dt - \lambda\Delta T = \Phi + \lambda_{sf} (T_{fluid}-T)
9 * Dirichlet (imposed temperature) or Neumann (imposed normal flux) boundary conditions.
11 //============================================================================
13 /*! \class DiffusionEquation DiffusionEquation.hxx "DiffusionEquation.hxx"
14 * \brief Scalar heat equation for the Uranium rods temperature
15 * \details see \ref DiffusionEqPage for more details
16 * rho*cp*dT/dt - \lambda\Delta T = \Phi + \lambda_{sf} (T_{fluid}-T)
18 #ifndef DiffusionEquation_HXX_
19 #define DiffusionEquation_HXX_
21 #include "ProblemCoreFlows.hxx"
26 //! enumeration BoundaryType
27 /*! Boundary condition type */
28 enum BoundaryTypeDiffusion { NeumannDiffusion, DirichletDiffusion, NoneBCDiffusion};
30 /** \struct LimitField
31 * \brief value of some fields on the boundary */
32 struct LimitFieldDiffusion{
33 LimitFieldDiffusion(){bcType=NoneBCDiffusion; T=0; normalFlux=0;}
34 LimitFieldDiffusion(BoundaryTypeDiffusion _bcType, double _T, double _normalFlux){
35 bcType=_bcType; T=_T; normalFlux=_normalFlux;
38 BoundaryTypeDiffusion bcType;
39 double T; //for Dirichlet
40 double normalFlux; //for Neumann
43 class DiffusionEquation: public ProblemCoreFlows
47 /** \fn DiffusionEquation
48 * \brief Constructor for the temperature diffusion in a solid
49 * \param [in] int : space dimension
50 * \param [in] double : solid density
51 * \param [in] double : solid specific heat at constant pressure
52 * \param [in] double : solid conductivity
55 DiffusionEquation( int dim,bool FECalculation=true,double rho=10000,double cp=300,double lambda=5, MPI_Comm comm = MPI_COMM_WORLD);
57 //Gestion du calcul (ICoCo)
59 void terminate();//vide la mémoire et enregistre le résultat final
60 bool initTimeStep(double dt);
61 double computeTimeStep(bool & stop);//propose un pas de temps pour le calcul. Celà nécessite de discrétiser les opérateur (convection, diffusion, sources) et pour chacun d'employer la condition cfl. En cas de problème durant ce calcul (exemple t=tmax), renvoie stop=true
62 void abortTimeStep();//efface les inconnues calculées par solveTimeStep() et reinitialise dt à 0
63 bool iterateTimeStep(bool &ok);
65 void validateTimeStep();
67 /* Boundary conditions */
68 void setBoundaryFields(map<string, LimitFieldDiffusion> boundaryFields){
69 _limitField = boundaryFields;
71 /** \fn setDirichletBoundaryCondition
72 * \brief adds a new boundary condition of type Dirichlet
74 * \param [in] string : the name of the boundary
75 * \param [in] double : the value of the temperature at the boundary
78 void setDirichletBoundaryCondition(string groupName,double Temperature=0){
79 _limitField[groupName]=LimitFieldDiffusion(DirichletDiffusion,Temperature,-1);
81 /** \fn setDirichletBoundaryCondition
82 * \brief adds a new boundary condition of type Dirichlet
83 * \details Reads the boundary field in a med file
84 * \param [in] string : the name of the boundary
85 * \param [in] string : the file name
86 * \param [in] string : the field name
87 * \param [in] int : the time step number
88 * \param [in] int : int corresponding to the enum CELLS or NODES
91 void setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
92 void setDirichletBoundaryCondition(string groupName, Field bc_field){
93 _limitField[groupName]=LimitFieldDiffusion(DirichletDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
95 /** \fn setNeumannBoundaryCondition
96 * \brief adds a new boundary condition of type Neumann
98 * \param [in] string : the name of the boundary
101 void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
102 _limitField[groupName]=LimitFieldDiffusion(NeumannDiffusion,-1, normalFlux);
104 /** \fn setNeumannBoundaryCondition
105 * \brief adds a new boundary condition of type Neumann
106 * \details Reads the boundary field in a med file
107 * \param [in] string : the name of the boundary
108 * \param [in] string : the file name
109 * \param [in] string : the field name
110 * \param [in] int : the time step number
111 * \param [in] int : int corresponding to the enum CELLS or NODES
114 void setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
115 void setNeumannBoundaryCondition(string groupName, Field BC_Field){
116 _limitField[groupName]=LimitFieldDiffusion(NeumannDiffusion,-1, 0);//This line will be deleted when variable BC are properly treated in solverlab
118 void setDirichletValues(map< int, double> dirichletBoundaryValues);
119 void setNeumannValues(map< int, double> neumannBoundaryValues);
121 void setRodDensity(double rho){
124 void setConductivity(double conductivite){
125 _conductivity=conductivite;
127 void setDiffusiontensor(Matrix DiffusionTensor){
128 _DiffusionTensor=DiffusionTensor;
132 /** Set input fields to prepare the simulation or coupling **/
133 vector<string> getInputFieldsNames();
134 void setInputField(const string& nameField, Field& inputField );//supply of a required input field
136 void setFluidTemperatureField(Field coupledTemperatureField);
137 void setFluidTemperature(double fluidTemperature){ _fluidTemperature=fluidTemperature; }
138 Field& getFluidTemperatureField(){ return _fluidTemperatureField; }
140 /*** get output fields names for postprocessing or coupling ***/
141 vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
142 Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement
144 Field& getOutputTemperatureField();//Return the main unknown if present (initialize() should be called first)
145 Field& getRodTemperatureField();//Return the main unknown if present (initialize() should be called first)
147 /*********** Generic functions for finite element method ***********/
148 static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
149 static int fact(int n);
150 static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
151 static int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
154 double computeDiffusionMatrix(bool & stop);
155 double computeDiffusionMatrixFV(bool & stop);
156 double computeDiffusionMatrixFE(bool & stop);
157 double computeRHS(bool & stop);
159 Field _fluidTemperatureField;
160 bool _fluidTemperatureFieldSet, _diffusionMatrixSet;
161 double _conductivity,_diffusivity, _fluidTemperature;
165 Matrix _DiffusionTensor;
166 Vec _Tn, _deltaT, _Tk, _Tkm1, _b0;
167 Vec _Tn_seq; // Local sequential copy of the parallel vector _Tn, used for saving result files
169 double _dt_diffusion, _dt_src;
170 map<string, LimitFieldDiffusion> _limitField;
173 /************ Data for FE calculation *************/
174 int _NunknownNodes;/* number of unknown nodes for FE calculation */
175 int _NboundaryNodes;/* total number of boundary nodes */
176 int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
177 std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
178 std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
180 /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
181 bool _dirichletValuesSet;
182 bool _neumannValuesSet;
183 std::map< int, double> _dirichletBoundaryValues;
184 std::map< int, double> _neumannBoundaryValues;
187 #endif /* DiffusionEquation_HXX_ */