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"
26 /*! Boundary condition type */
27 enum BoundaryTypeStationaryDiffusion { NeumannStationaryDiffusion, DirichletStationaryDiffusion, NoneBCStationaryDiffusion};
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;
37 BoundaryTypeStationaryDiffusion bcType;
38 double T; //for Dirichlet
39 double normalFlux; //for Neumann
42 class StationaryDiffusionEquation
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
53 StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1,MPI_Comm comm = MPI_COMM_WORLD);
55 void setMesh(const Mesh &M);
56 void setFileName(string fileName){
59 bool solveStationaryProblem();
61 //Linear system and spectrum
62 void setLinearSolver(linearSolver kspType, preconditioner pcType);
63 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
64 std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
65 std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
66 Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
70 void terminate();//vide la mémoire et enregistre le résultat final
71 double computeDiffusionMatrix(bool & stop);
72 double computeTimeStep(bool & stop);//For coupling calculations
73 bool iterateNewtonStep(bool &ok);
76 /* Boundary conditions */
77 void setBoundaryFields(map<string, LimitFieldStationaryDiffusion> boundaryFields){
78 _limitField = boundaryFields;
80 /** \fn setDirichletBoundaryCondition
81 * \brief adds a new boundary condition of type Dirichlet
83 * \param [in] string : the name of the boundary
84 * \param [in] double : the value of the temperature at the boundary
87 void setDirichletBoundaryCondition(string groupName,double Temperature){
88 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
90 /** \fn setDirichletBoundaryCondition
91 * \brief adds a new boundary condition of type Dirichlet
92 * \details Reads the boundary field in a med file
93 * \param [in] string : the name of the boundary
94 * \param [in] string : the file name
95 * \param [in] string : the field name
96 * \param [in] int : the time step number
97 * \param [in] int : int corresponding to the enum CELLS or NODES
100 void setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
101 void setDirichletBoundaryCondition(string groupName, Field bc_field){
102 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion, 0, -1);
105 /** \fn setNeumannBoundaryCondition
106 * \brief adds a new boundary condition of type Neumann
108 * \param [in] string : the name of the boundary
109 * \param [in] double : outward normal flux
112 void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
113 _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
115 /** \fn setNeumannBoundaryCondition
116 * \brief adds a new boundary condition of type Neumann
117 * \details Reads the boundary field in a med file
118 * \param [in] string : the name of the boundary
119 * \param [in] string : the file name
120 * \param [in] string : the field name
121 * \param [in] int : the time step number
122 * \param [in] int : int corresponding to the enum CELLS or NODES
125 void setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
126 void setNeumannBoundaryCondition(string groupName, Field bc_field){
127 _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, 0);
130 void setDirichletValues(map< int, double> dirichletBoundaryValues);
131 void setNeumannValues (map< int, double> neumannBoundaryValues);
133 void setConductivity(double conductivite){
134 _conductivity=conductivite;
136 void setDiffusiontensor(Matrix DiffusionTensor){
137 _DiffusionTensor=DiffusionTensor;
141 //get input fields to prepare the simulation or coupling
142 vector<string> getInputFieldsNames();
143 void setInputField(const string& nameField, Field& inputField );//supply of a required input field
145 void setFluidTemperatureField(Field coupledTemperatureField);
146 void setFluidTemperature(double fluidTemperature){ _fluidTemperature=fluidTemperature; }
147 Field& getFluidTemperatureField(){ return _fluidTemperatureField; }
149 /** \fn setHeatPowerField
150 * \brief set the heat power field (variable in space)
155 void setHeatPowerField(Field heatPower);
157 /** \fn setHeatPowerField
158 * \brief set the heat power field (variable in space)
160 * \param [in] string fileName (including file path)
161 * \param [in] string fieldName
164 void setHeatPowerField(string fileName, string fieldName, int iteration = 0, int order = 0, int meshLevel=0);
166 /** \fn getHeatPowerField
167 * \brief returns the heat power field
172 Field getHeatPowerField(){
173 return _heatPowerField;
175 //get output fields names for postprocessing or coupling
176 vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
177 Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement
179 Field& getOutputTemperatureField();
180 Field& getRodTemperatureField();
183 * \brief Updates display options
189 void setVerbose(bool verbose, bool system=false)
199 int _Ndim;//space dimension
200 int _nVar;//Number of equations to solve=1
205 bool _initializedMemory;
206 int _Nmailles;//number of cells for FV calculation
207 int _neibMaxNbCells;//maximum number of cells around a cell
210 double _precision_Newton;
211 double _erreur_rel;//norme(Uk+1-Uk)
212 bool _computationCompletedSuccessfully;
214 //Linear solver and petsc
220 int _maxPetscIts;//nombre maximum d'iteration gmres autorisé au cours d'une resolution de système lineaire
221 int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
222 int _PetscIts;//the number of iterations of the linear solver
224 Mat _A;//Linear system matrix
225 Vec _b;//Linear system right hand side
226 double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
227 bool _conditionNumber;//computes an estimate of the condition number
229 map<string, LimitFieldStationaryDiffusion> _limitField;
230 bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
232 bool _diffusionMatrixSet;
234 Matrix _DiffusionTensor;
235 Vec _deltaT, _Tk, _Tkm1, _b0;
236 Vec _Tk_seq; // Local sequential copy of the parallel vector _Tk, used for saving result files
240 //Heat transfert variables
241 double _conductivity, _fluidTemperature;
242 Field _heatPowerField, _fluidTemperatureField;
243 bool _heatPowerFieldSet, _fluidTemperatureFieldSet;
244 double _heatTransfertCoeff, _heatSource;
247 bool _verbose, _system;
248 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
250 string _fileName;//name of the calculation
251 string _path;//path to execution directory used for saving results
252 saveFormat _saveFormat;//file saving format : MED, VTK or CSV
254 double computeRHS(bool & stop);
255 double computeDiffusionMatrixFV(bool & stop);
256 double computeDiffusionMatrixFE(bool & stop);
258 /************ Data for FE calculation *************/
260 int _Nnodes;/* number of nodes for FE calculation */
261 int _neibMaxNbNodes;/* maximum number of nodes around a node */
262 int _NunknownNodes;/* number of unknown nodes for FE calculation */
263 int _NboundaryNodes;/* total number of boundary nodes */
264 int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
265 std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
266 std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
268 /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
269 bool _dirichletValuesSet;
270 bool _neumannValuesSet;
271 std::map< int, double> _dirichletBoundaryValues;
272 std::map< int, double> _neumannBoundaryValues;
274 /**** MPI related variables ***/
275 PetscMPIInt _mpi_size; /* size of communicator */
276 PetscMPIInt _mpi_rank; /* processor rank */
277 VecScatter _scat; /* For the distribution of a local vector */
278 int _globalNbUnknowns, _localNbUnknowns;
279 int _d_nnz, _o_nnz; /* local and "non local" numbers of non zeros corfficients */
282 #endif /* StationaryDiffusionEquation_HXX_ */