]> SALOME platform Git repositories - tools/solverlab.git/blob - CoreFlows/Models/inc/StationaryDiffusionEquation.hxx
Salome HOME
7810a39d3cd48eaf8d428abedeeb8d8d6ae99d48
[tools/solverlab.git] / CoreFlows / Models / inc / StationaryDiffusionEquation.hxx
1 //============================================================================
2 /**
3  * \file StationaryDiffusionEquation.hxx
4  * \author Michael NDJINGA
5  * \version 1.0
6  * \date June 2019
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.
10  * */
11 //============================================================================
12
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)
17  */
18  
19 #ifndef StationaryDiffusionEquation_HXX_
20 #define StationaryDiffusionEquation_HXX_
21
22 #include "ProblemCoreFlows.hxx"
23
24 using namespace std;
25
26 /*! Boundary condition type  */
27 enum BoundaryTypeStationaryDiffusion    { NeumannStationaryDiffusion, DirichletStationaryDiffusion, NoneBCStationaryDiffusion};
28
29 /** \struct LimitField
30  * \brief value of some fields on the boundary  */
31 struct LimitFieldStationaryDiffusion{
32         LimitFieldStationaryDiffusion(){bcType=NoneBCStationaryDiffusion; T=0; normalFlux=0;}
33         LimitFieldStationaryDiffusion(BoundaryTypeStationaryDiffusion _bcType, double _T,       double _normalFlux){
34                 bcType=_bcType; T=_T; normalFlux=_normalFlux;
35         }
36
37         BoundaryTypeStationaryDiffusion bcType;
38         double T; //for Dirichlet
39         double normalFlux; //for Neumann
40 };
41
42 class StationaryDiffusionEquation
43 {
44
45 public :
46         /** \fn StationaryDiffusionEquation
47                          * \brief Constructor for the temperature diffusion in a solid
48                          * \param [in] int : space dimension
49                          * \param [in] bool : numerical method
50                          * \param [in] double : solid conductivity
51                          *  */
52
53         StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1);
54
55     void setMesh(const Mesh &M);
56     void setFileName(string fileName){
57         _fileName = fileName;
58     }
59     bool solveStationaryProblem();
60     Field getOutputTemperatureField();
61     
62     //Linear system and spectrum
63     void setLinearSolver(linearSolver kspType, preconditioner pcType);
64     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
65     std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
66     std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
67     Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
68
69         //Gestion du calcul
70         void initialize();
71         void terminate();//vide la mémoire et enregistre le résultat final
72         double computeDiffusionMatrix(bool & stop);
73     double computeTimeStep(bool & stop);//For coupling calculations
74         bool iterateNewtonStep(bool &ok);
75         void save();
76
77     /* Boundary conditions */
78         void setBoundaryFields(map<string, LimitFieldStationaryDiffusion> boundaryFields){
79                 _limitField = boundaryFields;
80     };
81         /** \fn setDirichletBoundaryCondition
82                          * \brief adds a new boundary condition of type Dirichlet
83                          * \details
84                          * \param [in] string : the name of the boundary
85                          * \param [in] double : the value of the temperature at the boundary
86                          * \param [out] void
87                          *  */
88         void setDirichletBoundaryCondition(string groupName,double Temperature){
89                 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
90         };
91
92         /** \fn setNeumannBoundaryCondition
93                          * \brief adds a new boundary condition of type Neumann
94                          * \details
95                          * \param [in] string : the name of the boundary
96                          * \param [in] double : outward normal flux
97                          * \param [out] void
98                          *  */
99         void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
100                 _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
101         };
102
103         void setDirichletValues(map< int, double> dirichletBoundaryValues);
104         void setNeumannValues  (map< int, double>   neumannBoundaryValues);
105         
106         void setConductivity(double conductivite){
107                 _conductivity=conductivite;
108         };
109         void setFluidTemperatureField(Field coupledTemperatureField){
110                 _fluidTemperatureField=coupledTemperatureField;
111                 _fluidTemperatureFieldSet=true;
112         };
113         void setFluidTemperature(double fluidTemperature){
114         _fluidTemperature=fluidTemperature;
115         }
116         Field& getRodTemperatureField(){
117                 return _VV;
118         }
119         Field& getFluidTemperatureField(){
120                 return _fluidTemperatureField;
121         }
122         void setDiffusiontensor(Matrix DiffusionTensor){
123                 _DiffusionTensor=DiffusionTensor;
124         };
125
126         /** \fn setHeatPowerField
127          * \brief set the heat power field (variable in space)
128          * \details
129          * \param [in] Field
130          * \param [out] void
131          *  */
132         void setHeatPowerField(Field heatPower){
133                 _heatPowerField=heatPower;
134                 _heatPowerFieldSet=true;
135         }
136
137         /** \fn setHeatPowerField
138          * \brief set the heat power field (variable in space)
139          * \details
140          * \param [in] string fileName (including file path)
141          * \param [in] string fieldName
142          * \param [out] void
143          *  */
144         void setHeatPowerField(string fileName, string fieldName){
145                 _heatPowerField=Field(fileName, CELLS,fieldName);
146                 _heatPowerFieldSet=true;
147         }
148
149 protected :
150         //Main unknown field
151         Field _VV;
152
153         int _Ndim;//space dimension
154         int _nVar;//Number of equations to solve=1
155
156     //Mesh data
157         Mesh _mesh;
158     bool _meshSet;
159         bool _initializedMemory;
160         int _Nmailles;//number of cells for FV calculation
161         int _neibMaxNbCells;//maximum number of cells around a cell
162     
163         double _precision;
164         double _precision_Newton;
165         double _erreur_rel;//norme(Uk+1-Uk)
166     bool _computationCompletedSuccessfully;
167     
168         //Linear solver and petsc
169         KSP _ksp;
170         KSPType _ksptype;
171         PC _pc;
172         PCType _pctype;
173         string _pc_hypre;
174         int _maxPetscIts;//nombre maximum d'iteration gmres autorisé au cours d'une resolution de système lineaire
175         int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
176         int _PetscIts;//the number of iterations of the linear solver
177         int _NEWTON_its;
178         Mat  _A;//Linear system matrix
179         Vec _b;//Linear system right hand side
180         double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
181         bool _conditionNumber;//computes an estimate of the condition number
182
183         map<string, LimitFieldStationaryDiffusion> _limitField;
184     bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
185     
186         bool _diffusionMatrixSet;
187         Vector _normale;
188         Matrix _DiffusionTensor;
189         Vec _deltaT, _Tk, _Tkm1, _b0;
190         double _dt_src;
191     
192         //Heat transfert variables
193         double _conductivity, _fluidTemperature;
194         Field _heatPowerField, _fluidTemperatureField;
195         bool _heatPowerFieldSet, _fluidTemperatureFieldSet;
196         double _heatTransfertCoeff, _heatSource;
197
198         //Display variables
199         bool _verbose, _system;
200         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
201     //saving parameters
202         string _fileName;//name of the calculation
203         string _path;//path to execution directory used for saving results
204         saveFormat _saveFormat;//file saving format : MED, VTK or CSV
205
206         double computeRHS(bool & stop);
207         double computeDiffusionMatrixFV(bool & stop);
208
209     /************ Data for FE calculation *************/
210     bool _FECalculation;
211         int _Nnodes;/* number of nodes for FE calculation */
212         int _neibMaxNbNodes;/* maximum number of nodes around a node */
213         int _NunknownNodes;/* number of unknown nodes for FE calculation */
214         int _NboundaryNodes;/* total number of boundary nodes */
215         int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
216     std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
217     std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
218
219     /*********** Functions for finite element method ***********/
220     Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
221         double computeDiffusionMatrixFE(bool & stop);
222     static int fact(int n);
223     static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
224     static int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
225
226     /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
227     bool _dirichletValuesSet;
228     bool _neumannValuesSet;
229     std::map< int, double> _dirichletBoundaryValues;
230     std::map< int, double> _neumannBoundaryValues;
231 };
232
233 #endif /* StationaryDiffusionEquation_HXX_ */