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_ */