1 //============================================================================
3 * \file StationaryDiffusionEquation.hxx
4 * \author Michael NDJINGA
7 * \brief Stationary heat diffusion equation solved with either finite elements or finite volume method.
8 * -\lambda\Delta T=\Phi + \lambda_{sf} (T_{fluid}-T)
9 * Dirichlet (imposed temperature) or Neumann (imposed normal flux) boundary conditions
11 //============================================================================
13 /*! \class StationaryDiffusionEquation StationaryDiffusionEquation.hxx "StationaryDiffusionEquation.hxx"
14 * \brief Scalar stationary heat equation solved with either finite elements or finite volume method.
15 * \details see \ref StationaryDiffusionEqPage for more details
16 * -\lambda\Delta T=\Phi(T) + \lambda_{sf} (T_{fluid}-T)
18 #ifndef StationaryDiffusionEquation_HXX_
19 #define StationaryDiffusionEquation_HXX_
21 #include "ProblemCoreFlows.hxx"
23 /* for the laplacian spectrum */
29 //! enumeration BoundaryType
30 /*! Boundary condition type */
31 enum BoundaryTypeStationaryDiffusion { NeumannStationaryDiffusion, DirichletStationaryDiffusion, NoneBCStationaryDiffusion};
33 /** \struct LimitField
34 * \brief value of some fields on the boundary */
35 struct LimitFieldStationaryDiffusion{
36 LimitFieldStationaryDiffusion(){bcType=NoneBCStationaryDiffusion; T=0; normalFlux=0;}
37 LimitFieldStationaryDiffusion(BoundaryTypeStationaryDiffusion _bcType, double _T, double _normalFlux){
38 bcType=_bcType; T=_T; normalFlux=_normalFlux;
41 BoundaryTypeStationaryDiffusion bcType;
42 double T; //for Dirichlet
43 double normalFlux; //for Neumann
46 class StationaryDiffusionEquation
50 /** \fn StationaryDiffusionEquation
51 * \brief Constructor for the temperature diffusion in a solid
52 * \param [in] int : space dimension
53 * \param [in] double : solid conductivity
56 StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1);
58 void setMesh(const Mesh &M);
59 void setFileName(string fileName){
62 bool solveStationaryProblem();
63 Field getOutputTemperatureField();
65 //Linear system and spectrum
66 void setLinearSolver(linearSolver kspType, preconditioner pcType);
67 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
68 std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
69 std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
70 Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
74 void terminate();//vide la mémoire et enregistre le résultat final
75 double computeDiffusionMatrix(bool & stop);
76 double computeTimeStep(bool & stop);//For coupling calculations
77 bool iterateNewtonStep(bool &ok);
80 /* Boundary conditions */
81 void setBoundaryFields(map<string, LimitFieldStationaryDiffusion> boundaryFields){
82 _limitField = boundaryFields;
84 /** \fn setDirichletBoundaryCondition
85 * \brief adds a new boundary condition of type Dirichlet
87 * \param [in] string : the name of the boundary
88 * \param [in] double : the value of the temperature at the boundary
91 void setDirichletBoundaryCondition(string groupName,double Temperature){
92 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
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]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
105 void setDirichletValues(map< int, double> dirichletBoundaryValues);
106 void setNeumannValues (map< int, double> neumannBoundaryValues);
108 void setConductivity(double conductivite){
109 _conductivity=conductivite;
111 void setFluidTemperatureField(Field coupledTemperatureField){
112 _fluidTemperatureField=coupledTemperatureField;
113 _fluidTemperatureFieldSet=true;
115 void setFluidTemperature(double fluidTemperature){
116 _fluidTemperature=fluidTemperature;
118 Field& getRodTemperatureField(){
121 Field& getFluidTemperatureField(){
122 return _fluidTemperatureField;
124 void setDiffusiontensor(Matrix DiffusionTensor){
125 _DiffusionTensor=DiffusionTensor;
128 /** \fn setHeatPowerField
129 * \brief set the heat power field (variable in space)
134 void setHeatPowerField(Field heatPower){
135 _heatPowerField=heatPower;
136 _heatPowerFieldSet=true;
139 /** \fn setHeatPowerField
140 * \brief set the heat power field (variable in space)
142 * \param [in] string fileName (including file path)
143 * \param [in] string fieldName
146 void setHeatPowerField(string fileName, string fieldName){
147 _heatPowerField=Field(fileName, CELLS,fieldName);
148 _heatPowerFieldSet=true;
155 int _Ndim;//space dimension
156 int _nVar;//Number of equations to solve=1
161 bool _initializedMemory;
162 int _Nmailles;//number of cells for FV calculation
163 int _neibMaxNbCells;//maximum number of cells around a cell
166 double _precision_Newton;
167 double _erreur_rel;//norme(Uk+1-Uk)
168 bool _computationCompletedSuccessfully;
170 //Linear solver and petsc
176 int _maxPetscIts;//nombre maximum d'iteration gmres autorisé au cours d'une resolution de système lineaire
177 int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
178 int _PetscIts;//the number of iterations of the linear solver
180 Mat _A;//Linear system matrix
181 Vec _b;//Linear system right hand side
182 double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
183 bool _conditionNumber;//computes an estimate of the condition number
185 map<string, LimitFieldStationaryDiffusion> _limitField;
186 bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
188 bool _diffusionMatrixSet;
190 Matrix _DiffusionTensor;
191 Vec _deltaT, _Tk, _Tkm1, _b0;
194 //Heat transfert variables
195 double _conductivity, _fluidTemperature;
196 Field _heatPowerField, _fluidTemperatureField;
197 bool _heatPowerFieldSet, _fluidTemperatureFieldSet;
198 double _heatTransfertCoeff, _heatSource;
201 bool _verbose, _system;
202 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
204 string _fileName;//name of the calculation
205 string _path;//path to execution directory used for saving results
206 saveFormat _saveFormat;//file saving format : MED, VTK or CSV
208 double computeRHS(bool & stop);
209 double computeDiffusionMatrixFV(bool & stop);
211 /************ Data for FE calculation *************/
213 int _Nnodes;/* number of nodes for FE calculation */
214 int _neibMaxNbNodes;/* maximum number of nodes around a node */
215 int _NunknownNodes;/* number of unknown nodes for FE calculation */
216 int _NboundaryNodes;/* total number of boundary nodes */
217 int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
218 std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
219 std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
221 /*********** Functions for finite element method ***********/
222 Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
223 double computeDiffusionMatrixFE(bool & stop);
225 int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes) const;
226 int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes) const;
228 /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
229 bool _dirichletValuesSet;
230 bool _neumannValuesSet;
231 std::map< int, double> _dirichletBoundaryValues;
232 std::map< int, double> _neumannBoundaryValues;
235 #endif /* StationaryDiffusionEquation_HXX_ */