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,MPI_Comm comm = MPI_COMM_WORLD);
  54
  55     void setMesh(const Mesh &M);
  56     void setFileName(string fileName){
  57         _fileName = fileName;
  58     }
  59     bool solveStationaryProblem();
  60
  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;
  67
  68         //Gestion du calcul
  69         void initialize();
  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);
  74         void save();
  75
  76     /* Boundary conditions */
  77         void setBoundaryFields(map<string, LimitFieldStationaryDiffusion> boundaryFields){
  78                 _limitField = boundaryFields;
  79     };
  80         /** \fn setDirichletBoundaryCondition
  81                          * \brief adds a new boundary condition of type Dirichlet
  82                          * \details
  83                          * \param [in] string : the name of the boundary
  84                          * \param [in] double : the value of the temperature at the boundary
  85                          * \param [out] void
  86                          *  */
  87         void setDirichletBoundaryCondition(string groupName,double Temperature){
  88                 _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
  89         };
  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
  98                          * \param [out] void
  99                          *  */
 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);
 103         };
 104
 105         /** \fn setNeumannBoundaryCondition
 106                          * \brief adds a new boundary condition of type Neumann
 107                          * \details
 108                          * \param [in] string : the name of the boundary
 109                          * \param [in] double : outward normal flux
 110                          * \param [out] void
 111                          *  */
 112         void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
 113                 _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
 114         };
 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
 123                          * \param [out] void
 124                          *  */
 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);
 128         };
 129
 130         void setDirichletValues(map< int, double> dirichletBoundaryValues);
 131         void setNeumannValues  (map< int, double>   neumannBoundaryValues);
 132
 133         void setConductivity(double conductivite){
 134                 _conductivity=conductivite;
 135         };
 136         void setDiffusiontensor(Matrix DiffusionTensor){
 137                 _DiffusionTensor=DiffusionTensor;
 138         };
 139
 140
 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
 144
 145         void setFluidTemperatureField(Field coupledTemperatureField);
 146         void setFluidTemperature(double fluidTemperature){      _fluidTemperature=fluidTemperature;     }
 147         Field& getFluidTemperatureField(){      return _fluidTemperatureField;  }
 148
 149         /** \fn setHeatPowerField
 150          * \brief set the heat power field (variable in space)
 151          * \details
 152          * \param [in] Field
 153          * \param [out] void
 154          *  */
 155         void setHeatPowerField(Field heatPower);
 156
 157         /** \fn setHeatPowerField
 158          * \brief set the heat power field (variable in space)
 159          * \details
 160          * \param [in] string fileName (including file path)
 161          * \param [in] string fieldName
 162          * \param [out] void
 163          *  */
 164         void setHeatPowerField(string fileName, string fieldName, int iteration = 0, int order = 0, int meshLevel=0);
 165
 166         /** \fn getHeatPowerField
 167          * \brief returns the heat power field
 168          * \details
 169          * \param [in] void
 170          * \param [out] Field
 171          *  */
 172         Field getHeatPowerField(){
 173                 return _heatPowerField;
 174         }
 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
 178
 179     Field& getOutputTemperatureField();
 180         Field& getRodTemperatureField();
 181
 182         /** \fn setVerbose
 183          * \brief Updates display options
 184          * \details
 185          * \param [in] bool
 186          * \param [in] bool
 187          * \param [out] void
 188          *  */
 189         void setVerbose(bool verbose,  bool system=false)
 190         {
 191                 _verbose = verbose;
 192                 _system = system;
 193         };
 194
 195 protected :
 196         //Main unknown field
 197         Field _VV;
 198
 199         int _Ndim;//space dimension
 200         int _nVar;//Number of equations to solve=1
 201
 202     //Mesh data
 203         Mesh _mesh;
 204     bool _meshSet;
 205         bool _initializedMemory;
 206         int _Nmailles;//number of cells for FV calculation
 207         int _neibMaxNbCells;//maximum number of cells around a cell
 208
 209         double _precision;
 210         double _precision_Newton;
 211         double _erreur_rel;//norme(Uk+1-Uk)
 212     bool _computationCompletedSuccessfully;
 213
 214         //Linear solver and petsc
 215         KSP _ksp;
 216         KSPType _ksptype;
 217         PC _pc;
 218         PCType _pctype;
 219         string _pc_hypre;
 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
 223         int _NEWTON_its;
 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
 228
 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
 231
 232         bool _diffusionMatrixSet;
 233         Vector _normale;
 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
 237
 238         double _dt_src;
 239
 240         //Heat transfert variables
 241         double _conductivity, _fluidTemperature;
 242         Field _heatPowerField, _fluidTemperatureField;
 243         bool _heatPowerFieldSet, _fluidTemperatureFieldSet;
 244         double _heatTransfertCoeff, _heatSource;
 245
 246         //Display variables
 247         bool _verbose, _system;
 248         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
 249     //saving parameters
 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
 253
 254         double computeRHS(bool & stop);
 255         double computeDiffusionMatrixFV(bool & stop);
 256         double computeDiffusionMatrixFE(bool & stop);
 257
 258     /************ Data for FE calculation *************/
 259     bool _FECalculation;
 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 */
 267
 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;
 273
 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 */
 280 };
 281
 282 #endif /* StationaryDiffusionEquation_HXX_ */