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