]> SALOME platform Git repositories - tools/solverlab.git/blob - CoreFlows/Models/inc/StationaryDiffusionEquation.hxx
Salome HOME
update CoreFlows
[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 #ifndef StationaryDiffusionEquation_HXX_
19 #define StationaryDiffusionEquation_HXX_
20
21 #include "ProblemCoreFlows.hxx"
22
23 /* for the laplacian spectrum */
24 #include <slepceps.h>
25 #include <slepcsvd.h>
26
27 using namespace std;
28
29 //! enumeration BoundaryType
30 /*! Boundary condition type  */
31 enum BoundaryTypeStationaryDiffusion    { NeumannStationaryDiffusion, DirichletStationaryDiffusion, NoneBCStationaryDiffusion};
32
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;
39         }
40
41         BoundaryTypeStationaryDiffusion bcType;
42         double T; //for Dirichlet
43         double normalFlux; //for Neumann
44 };
45
46 class StationaryDiffusionEquation
47 {
48
49 public :
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
54                          *  */
55
56         StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1);
57
58     void setMesh(const Mesh &M);
59     void setFileName(string fileName){
60         _fileName = fileName;
61     }
62     bool solveStationaryProblem();
63     Field getOutputTemperatureField();
64     
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;
71
72         //Gestion du calcul
73         void initialize();
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);
78         void save();
79
80     /* Boundary conditions */
81         void setBoundaryFields(map<string, LimitFieldStationaryDiffusion> boundaryFields){
82                 _limitField = boundaryFields;
83     };
84         /** \fn setDirichletBoundaryCondition
85                          * \brief adds a new boundary condition of type Dirichlet
86                          * \details
87                          * \param [in] string : the name of the boundary
88                          * \param [in] double : the value of the temperature at the boundary
89                          * \param [out] void
90                          *  */
91         void setDirichletBoundaryCondition(string groupName,double Temperature){
92                 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
93         };
94
95         /** \fn setNeumannBoundaryCondition
96                          * \brief adds a new boundary condition of type Neumann
97                          * \details
98                          * \param [in] string : the name of the boundary
99                          * \param [out] void
100                          *  */
101         void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
102                 _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
103         };
104
105         void setDirichletValues(map< int, double> dirichletBoundaryValues);
106         void setNeumannValues  (map< int, double>   neumannBoundaryValues);
107         
108         void setConductivity(double conductivite){
109                 _conductivity=conductivite;
110         };
111         void setFluidTemperatureField(Field coupledTemperatureField){
112                 _fluidTemperatureField=coupledTemperatureField;
113                 _fluidTemperatureFieldSet=true;
114         };
115         void setFluidTemperature(double fluidTemperature){
116         _fluidTemperature=fluidTemperature;
117         }
118         Field& getRodTemperatureField(){
119                 return _VV;
120         }
121         Field& getFluidTemperatureField(){
122                 return _fluidTemperatureField;
123         }
124         void setDiffusiontensor(Matrix DiffusionTensor){
125                 _DiffusionTensor=DiffusionTensor;
126         };
127
128         /** \fn setHeatPowerField
129          * \brief set the heat power field (variable in space)
130          * \details
131          * \param [in] Field
132          * \param [out] void
133          *  */
134         void setHeatPowerField(Field heatPower){
135                 _heatPowerField=heatPower;
136                 _heatPowerFieldSet=true;
137         }
138
139         /** \fn setHeatPowerField
140          * \brief set the heat power field (variable in space)
141          * \details
142          * \param [in] string fileName (including file path)
143          * \param [in] string fieldName
144          * \param [out] void
145          *  */
146         void setHeatPowerField(string fileName, string fieldName){
147                 _heatPowerField=Field(fileName, CELLS,fieldName);
148                 _heatPowerFieldSet=true;
149         }
150
151 protected :
152         //Main unknown field
153         Field _VV;
154
155         int _Ndim;//space dimension
156         int _nVar;//Number of equations to solve=1
157
158     //Mesh data
159         Mesh _mesh;
160     bool _meshSet;
161         bool _initializedMemory;
162         int _Nmailles;//number of cells for FV calculation
163         int _neibMaxNbCells;//maximum number of cells around a cell
164     
165         double _precision;
166         double _precision_Newton;
167         double _erreur_rel;//norme(Uk+1-Uk)
168     bool _computationCompletedSuccessfully;
169     
170         //Linear solver and petsc
171         KSP _ksp;
172         KSPType _ksptype;
173         PC _pc;
174         PCType _pctype;
175         string _pc_hypre;
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
179         int _NEWTON_its;
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
184
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
187     
188         bool _diffusionMatrixSet;
189         Vector _normale;
190         Matrix _DiffusionTensor;
191         Vec _deltaT, _Tk, _Tkm1, _b0;
192         double _dt_src;
193     
194         //Heat transfert variables
195         double _conductivity, _fluidTemperature;
196         Field _heatPowerField, _fluidTemperatureField;
197         bool _heatPowerFieldSet, _fluidTemperatureFieldSet;
198         double _heatTransfertCoeff, _heatSource;
199
200         //Display variables
201         bool _verbose, _system;
202         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
203     //saving parameters
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
207
208         double computeRHS(bool & stop);
209         double computeDiffusionMatrixFV(bool & stop);
210
211     /************ Data for FE calculation *************/
212     bool _FECalculation;
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 */
220
221     /*********** Functions for finite element method ***********/
222     Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
223         double computeDiffusionMatrixFE(bool & stop);
224     int fact(int n);
225     int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes) const;
226     int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes) const;
227
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;
233 };
234
235 #endif /* StationaryDiffusionEquation_HXX_ */