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