CoreFlows/Models/inc/DiffusionEquation.hxx

   1 //============================================================================
   2 /**
   3  * \file DiffusionEquation.hxx
   4  * \author Michael NDJINGA
   5  * \version 1.0
   6  * \date 24 March 2015
   7  * \brief Heat diffusion equation
   8  * rho*cp*dT/dt - \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 DiffusionEquation DiffusionEquation.hxx "DiffusionEquation.hxx"
  14  *  \brief Scalar heat equation for the Uranium rods temperature
  15  *  \details see \ref DiffusionEqPage for more details
  16  * rho*cp*dT/dt - \lambda\Delta T = \Phi + \lambda_{sf} (T_{fluid}-T)
  17  */
  18 #ifndef DiffusionEquation_HXX_
  19 #define DiffusionEquation_HXX_
  20
  21 #include "ProblemCoreFlows.hxx"
  22 #include "Node.hxx"
  23
  24 using namespace std;
  25
  26 //! enumeration BoundaryType
  27 /*! Boundary condition type  */
  28 enum BoundaryTypeDiffusion      { NeumannDiffusion, DirichletDiffusion, NoneBCDiffusion};
  29
  30 /** \struct LimitField
  31  * \brief value of some fields on the boundary  */
  32 struct LimitFieldDiffusion{
  33         LimitFieldDiffusion(){bcType=NoneBCDiffusion; T=0; normalFlux=0;}
  34         LimitFieldDiffusion(BoundaryTypeDiffusion _bcType, double _T,   double _normalFlux){
  35                 bcType=_bcType; T=_T; normalFlux=_normalFlux;
  36         }
  37
  38         BoundaryTypeDiffusion bcType;
  39         double T; //for Dirichlet
  40         double normalFlux; //for Neumann
  41 };
  42
  43 class DiffusionEquation: public ProblemCoreFlows
  44 {
  45
  46 public :
  47         /** \fn DiffusionEquation
  48                          * \brief Constructor for the temperature diffusion in a solid
  49                          * \param [in] int : space dimension
  50                          * \param [in] double : solid density
  51                          * \param [in] double : solid specific heat at constant pressure
  52                          * \param [in] double : solid conductivity
  53                          *  */
  54
  55         DiffusionEquation( int dim,bool FECalculation=true,double rho=10000,double cp=300,double lambda=5, MPI_Comm comm = MPI_COMM_WORLD);
  56
  57         //Gestion du calcul (ICoCo)
  58         void initialize();
  59         void terminate();//vide la mémoire et enregistre le résultat final
  60         bool initTimeStep(double dt);
  61         double computeTimeStep(bool & stop);//propose un pas de temps pour le calcul. Celà nécessite de discrétiser les opérateur (convection, diffusion, sources) et pour chacun d'employer la condition cfl. En cas de problème durant ce calcul (exemple t=tmax), renvoie stop=true
  62         void abortTimeStep();//efface les inconnues calculées par solveTimeStep() et reinitialise dt à 0
  63         bool iterateTimeStep(bool &ok);
  64         void save();
  65         void validateTimeStep();
  66
  67     /* Boundary conditions */
  68         void setBoundaryFields(map<string, LimitFieldDiffusion> boundaryFields){
  69                 _limitField = boundaryFields;
  70     };
  71         /** \fn setDirichletBoundaryCondition
  72                          * \brief adds a new boundary condition of type Dirichlet
  73                          * \details
  74                          * \param [in] string : the name of the boundary
  75                          * \param [in] double : the value of the temperature at the boundary
  76                          * \param [out] void
  77                          *  */
  78         void setDirichletBoundaryCondition(string groupName,double Temperature=0){
  79                 _limitField[groupName]=LimitFieldDiffusion(DirichletDiffusion,Temperature,-1);
  80         };
  81         /** \fn setDirichletBoundaryCondition
  82                          * \brief adds a new boundary condition of type Dirichlet
  83                          * \details Reads the boundary field in a med file
  84                          * \param [in] string : the name of the boundary
  85                          * \param [in] string : the file name
  86                          * \param [in] string : the field name
  87                          * \param [in] int : the time step number
  88                          * \param [in] int : int corresponding to the enum CELLS or NODES
  89                          * \param [out] void
  90                          *  */
  91         void setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
  92         void setDirichletBoundaryCondition(string groupName, Field bc_field){
  93                 _limitField[groupName]=LimitFieldDiffusion(DirichletDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
  94         }
  95         /** \fn setNeumannBoundaryCondition
  96                          * \brief adds a new boundary condition of type Neumann
  97                          * \details
  98                          * \param [in] string : the name of the boundary
  99                          * \param [out] void
 100                          *  */
 101         void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
 102                 _limitField[groupName]=LimitFieldDiffusion(NeumannDiffusion,-1, normalFlux);
 103         };
 104         /** \fn setNeumannBoundaryCondition
 105                          * \brief adds a new boundary condition of type Neumann
 106                          * \details Reads the boundary field in a med file
 107                          * \param [in] string : the name of the boundary
 108                          * \param [in] string : the file name
 109                          * \param [in] string : the field name
 110                          * \param [in] int : the time step number
 111                          * \param [in] int : int corresponding to the enum CELLS or NODES
 112                          * \param [out] void
 113                          *  */
 114         void setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
 115         void setNeumannBoundaryCondition(string groupName, Field BC_Field){
 116                 _limitField[groupName]=LimitFieldDiffusion(NeumannDiffusion,-1, 0);//This line will be deleted when variable BC are properly treated in solverlab
 117         };
 118         void setDirichletValues(map< int, double> dirichletBoundaryValues);
 119         void setNeumannValues(map< int, double> neumannBoundaryValues);
 120
 121         void setRodDensity(double rho){
 122                 _rho=rho;
 123         };
 124         void setConductivity(double conductivite){
 125                 _conductivity=conductivite;
 126         };
 127         void setDiffusiontensor(Matrix DiffusionTensor){
 128                 _DiffusionTensor=DiffusionTensor;
 129         };
 130
 131
 132         /** Set input fields to prepare the simulation or coupling **/
 133         vector<string> getInputFieldsNames();
 134         void setInputField(const string& nameField, Field& inputField );//supply of a required input field
 135
 136         void setFluidTemperatureField(Field coupledTemperatureField);
 137         void setFluidTemperature(double fluidTemperature){      _fluidTemperature=fluidTemperature;     }
 138         Field& getFluidTemperatureField(){  return _fluidTemperatureField;      }
 139
 140         /*** get output fields names for postprocessing or coupling ***/
 141         vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
 142         Field&         getOutputField(const string& nameField );//Renvoie un champs pour le postraitement
 143
 144         Field& getOutputTemperatureField();//Return the main unknown if present (initialize() should be called first)
 145         Field& getRodTemperatureField();//Return the main unknown if present (initialize() should be called first)
 146
 147     /*********** Generic functions for finite element method ***********/
 148     static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
 149     static int fact(int n);
 150     static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
 151     static int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
 152
 153 protected :
 154         double computeDiffusionMatrix(bool & stop);
 155         double computeDiffusionMatrixFV(bool & stop);
 156         double computeDiffusionMatrixFE(bool & stop);
 157         double computeRHS(bool & stop);
 158
 159         Field _fluidTemperatureField;
 160         bool _fluidTemperatureFieldSet, _diffusionMatrixSet;
 161         double _conductivity,_diffusivity, _fluidTemperature;
 162         double _rho;
 163         double _cp;
 164         Vector _normale;
 165         Matrix _DiffusionTensor;
 166         Vec _Tn, _deltaT, _Tk, _Tkm1, _b0;
 167         Vec _Tn_seq; // Local sequential copy of the parallel vector _Tn, used for saving result files
 168
 169         double _dt_diffusion, _dt_src;
 170         map<string, LimitFieldDiffusion> _limitField;
 171
 172
 173     /************ Data for FE calculation *************/
 174         int _NunknownNodes;/* number of unknown nodes for FE calculation */
 175         int _NboundaryNodes;/* total number of boundary nodes */
 176         int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
 177     std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
 178     std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
 179
 180     /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
 181     bool _dirichletValuesSet;
 182     bool _neumannValuesSet;
 183     std::map< int, double> _dirichletBoundaryValues;
 184     std::map< int, double> _neumannBoundaryValues;
 185 };
 186
 187 #endif /* DiffusionEquation_HXX_ */