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)
19 #ifndef StationaryDiffusionEquation_HXX_
20 #define StationaryDiffusionEquation_HXX_
22 #include "ProblemCoreFlows.hxx"
24 /* For the laplacian spectrum */
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] bool : numerical method
54 * \param [in] double : solid conductivity
57 StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1);
59 void setMesh(const Mesh &M);
60 void setFileName(string fileName){
63 bool solveStationaryProblem();
64 Field getOutputTemperatureField();
66 //Linear system and spectrum
67 void setLinearSolver(linearSolver kspType, preconditioner pcType);
68 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
69 std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
70 std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
71 Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
75 void terminate();//vide la mémoire et enregistre le résultat final
76 double computeDiffusionMatrix(bool & stop);
77 double computeTimeStep(bool & stop);//For coupling calculations
78 bool iterateNewtonStep(bool &ok);
81 /* Boundary conditions */
82 void setBoundaryFields(map<string, LimitFieldStationaryDiffusion> boundaryFields){
83 _limitField = boundaryFields;
85 /** \fn setDirichletBoundaryCondition
86 * \brief adds a new boundary condition of type Dirichlet
88 * \param [in] string : the name of the boundary
89 * \param [in] double : the value of the temperature at the boundary
92 void setDirichletBoundaryCondition(string groupName,double Temperature){
93 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
96 /** \fn setNeumannBoundaryCondition
97 * \brief adds a new boundary condition of type Neumann
99 * \param [in] string : the name of the boundary
100 * \param [in] double : outward normal flux
103 void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
104 _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
107 void setDirichletValues(map< int, double> dirichletBoundaryValues);
108 void setNeumannValues (map< int, double> neumannBoundaryValues);
110 void setConductivity(double conductivite){
111 _conductivity=conductivite;
113 void setFluidTemperatureField(Field coupledTemperatureField){
114 _fluidTemperatureField=coupledTemperatureField;
115 _fluidTemperatureFieldSet=true;
117 void setFluidTemperature(double fluidTemperature){
118 _fluidTemperature=fluidTemperature;
120 Field& getRodTemperatureField(){
123 Field& getFluidTemperatureField(){
124 return _fluidTemperatureField;
126 void setDiffusiontensor(Matrix DiffusionTensor){
127 _DiffusionTensor=DiffusionTensor;
130 /** \fn setHeatPowerField
131 * \brief set the heat power field (variable in space)
136 void setHeatPowerField(Field heatPower){
137 _heatPowerField=heatPower;
138 _heatPowerFieldSet=true;
141 /** \fn setHeatPowerField
142 * \brief set the heat power field (variable in space)
144 * \param [in] string fileName (including file path)
145 * \param [in] string fieldName
148 void setHeatPowerField(string fileName, string fieldName){
149 _heatPowerField=Field(fileName, CELLS,fieldName);
150 _heatPowerFieldSet=true;
157 int _Ndim;//space dimension
158 int _nVar;//Number of equations to solve=1
163 bool _initializedMemory;
164 int _Nmailles;//number of cells for FV calculation
165 int _neibMaxNbCells;//maximum number of cells around a cell
168 double _precision_Newton;
169 double _erreur_rel;//norme(Uk+1-Uk)
170 bool _computationCompletedSuccessfully;
172 //Linear solver and petsc
178 int _maxPetscIts;//nombre maximum d'iteration gmres autorisé au cours d'une resolution de système lineaire
179 int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
180 int _PetscIts;//the number of iterations of the linear solver
182 Mat _A;//Linear system matrix
183 Vec _b;//Linear system right hand side
184 double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
185 bool _conditionNumber;//computes an estimate of the condition number
187 map<string, LimitFieldStationaryDiffusion> _limitField;
188 bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
190 bool _diffusionMatrixSet;
192 Matrix _DiffusionTensor;
193 Vec _deltaT, _Tk, _Tkm1, _b0;
196 //Heat transfert variables
197 double _conductivity, _fluidTemperature;
198 Field _heatPowerField, _fluidTemperatureField;
199 bool _heatPowerFieldSet, _fluidTemperatureFieldSet;
200 double _heatTransfertCoeff, _heatSource;
203 bool _verbose, _system;
204 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
206 string _fileName;//name of the calculation
207 string _path;//path to execution directory used for saving results
208 saveFormat _saveFormat;//file saving format : MED, VTK or CSV
210 double computeRHS(bool & stop);
211 double computeDiffusionMatrixFV(bool & stop);
213 /************ Data for FE calculation *************/
215 int _Nnodes;/* number of nodes for FE calculation */
216 int _neibMaxNbNodes;/* maximum number of nodes around a node */
217 int _NunknownNodes;/* number of unknown nodes for FE calculation */
218 int _NboundaryNodes;/* total number of boundary nodes */
219 int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
220 std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
221 std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
223 /*********** Functions for finite element method ***********/
224 Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
225 double computeDiffusionMatrixFE(bool & stop);
227 int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes) const;
228 int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes) const;
230 /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
231 bool _dirichletValuesSet;
232 bool _neumannValuesSet;
233 std::map< int, double> _dirichletBoundaryValues;
234 std::map< int, double> _neumannBoundaryValues;
237 #endif /* StationaryDiffusionEquation_HXX_ */