* \date June 2019
* \brief Stationary heat diffusion equation solved with either finite elements or finite volume method.
* -\lambda\Delta T=\Phi + \lambda_{sf} (T_{fluid}-T)
+ * Dirichlet (imposed temperature) or Neumann (imposed normal flux) boundary conditions.
* */
//============================================================================
* \details see \ref StationaryDiffusionEqPage for more details
* -\lambda\Delta T=\Phi(T) + \lambda_{sf} (T_{fluid}-T)
*/
+
#ifndef StationaryDiffusionEquation_HXX_
#define StationaryDiffusionEquation_HXX_
#include "ProblemCoreFlows.hxx"
-#include "Node.hxx"
using namespace std;
+/*! Boundary condition type */
+enum BoundaryTypeStationaryDiffusion { NeumannStationaryDiffusion, DirichletStationaryDiffusion, NoneBCStationaryDiffusion};
+
+/** \struct LimitField
+ * \brief value of some fields on the boundary */
+struct LimitFieldStationaryDiffusion{
+ LimitFieldStationaryDiffusion(){bcType=NoneBCStationaryDiffusion; T=0; normalFlux=0;}
+ LimitFieldStationaryDiffusion(BoundaryTypeStationaryDiffusion _bcType, double _T, double _normalFlux){
+ bcType=_bcType; T=_T; normalFlux=_normalFlux;
+ }
+
+ BoundaryTypeStationaryDiffusion bcType;
+ double T; //for Dirichlet
+ double normalFlux; //for Neumann
+};
+
class StationaryDiffusionEquation
{
/** \fn StationaryDiffusionEquation
* \brief Constructor for the temperature diffusion in a solid
* \param [in] int : space dimension
- * \param [in] double : solid density
- * \param [in] double : solid specific heat at constant pressure
+ * \param [in] bool : numerical method
* \param [in] double : solid conductivity
* */
- StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1);
+ StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1,MPI_Comm comm = MPI_COMM_WORLD);
void setMesh(const Mesh &M);
- void setLinearSolver(linearSolver kspType, preconditioner pcType);
void setFileName(string fileName){
_fileName = fileName;
}
bool solveStationaryProblem();
- Field getOutputTemperatureField();
+ //Linear system and spectrum
+ void setLinearSolver(linearSolver kspType, preconditioner pcType);
+ double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
+ std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
+ std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
+ Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
+
//Gestion du calcul
void initialize();
void terminate();//vide la mémoire et enregistre le résultat final
void save();
/* Boundary conditions */
- void setBoundaryFields(map<string, LimitField> boundaryFields){
+ void setBoundaryFields(map<string, LimitFieldStationaryDiffusion> boundaryFields){
_limitField = boundaryFields;
};
/** \fn setDirichletBoundaryCondition
* \param [out] void
* */
void setDirichletBoundaryCondition(string groupName,double Temperature){
- _limitField[groupName]=LimitField(Dirichlet,-1, vector<double>(_Ndim,0),vector<double>(_Ndim,0),
- vector<double>(_Ndim,0),Temperature,-1,-1,-1);
+ _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
+ };
+ /** \fn setDirichletBoundaryCondition
+ * \brief adds a new boundary condition of type Dirichlet
+ * \details Reads the boundary field in a med file
+ * \param [in] string : the name of the boundary
+ * \param [in] string : the file name
+ * \param [in] string : the field name
+ * \param [in] int : the time step number
+ * \param [in] int : int corresponding to the enum CELLS or NODES
+ * \param [out] void
+ * */
+ void setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
+ void setDirichletBoundaryCondition(string groupName, Field bc_field){
+ _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion, 0, -1);
};
/** \fn setNeumannBoundaryCondition
* \brief adds a new boundary condition of type Neumann
* \details
* \param [in] string : the name of the boundary
+ * \param [in] double : outward normal flux
+ * \param [out] void
+ * */
+ void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
+ _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
+ };
+ /** \fn setNeumannBoundaryCondition
+ * \brief adds a new boundary condition of type Neumann
+ * \details Reads the boundary field in a med file
+ * \param [in] string : the name of the boundary
+ * \param [in] string : the file name
+ * \param [in] string : the field name
+ * \param [in] int : the time step number
+ * \param [in] int : int corresponding to the enum CELLS or NODES
* \param [out] void
* */
- void setNeumannBoundaryCondition(string groupName){
- _limitField[groupName]=LimitField(Neumann,-1, vector<double>(0),vector<double>(0),
- vector<double>(0),-1,-1,-1,-1);
+ void setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
+ void setNeumannBoundaryCondition(string groupName, Field bc_field){
+ _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, 0);
};
void setDirichletValues(map< int, double> dirichletBoundaryValues);
+ void setNeumannValues (map< int, double> neumannBoundaryValues);
void setConductivity(double conductivite){
_conductivity=conductivite;
};
- void setFluidTemperatureField(Field coupledTemperatureField){
- _fluidTemperatureField=coupledTemperatureField;
- _fluidTemperatureFieldSet=true;
- };
- void setFluidTemperature(double fluidTemperature){
- _fluidTemperature=fluidTemperature;
- }
- Field& getRodTemperatureField(){
- return _VV;
- }
- Field& getFluidTemperatureField(){
- return _fluidTemperatureField;
- }
void setDiffusiontensor(Matrix DiffusionTensor){
_DiffusionTensor=DiffusionTensor;
};
+
+ //get input fields to prepare the simulation or coupling
+ vector<string> getInputFieldsNames();
+ void setInputField(const string& nameField, Field& inputField );//supply of a required input field
+
+ void setFluidTemperatureField(Field coupledTemperatureField);
+ void setFluidTemperature(double fluidTemperature){ _fluidTemperature=fluidTemperature; }
+ Field& getFluidTemperatureField(){ return _fluidTemperatureField; }
+
/** \fn setHeatPowerField
* \brief set the heat power field (variable in space)
* \details
* \param [in] Field
* \param [out] void
* */
- void setHeatPowerField(Field heatPower){
- _heatPowerField=heatPower;
- _heatPowerFieldSet=true;
- }
+ void setHeatPowerField(Field heatPower);
/** \fn setHeatPowerField
* \brief set the heat power field (variable in space)
* \param [in] string fieldName
* \param [out] void
* */
- void setHeatPowerField(string fileName, string fieldName){
- _heatPowerField=Field(fileName, CELLS,fieldName);
- _heatPowerFieldSet=true;
+ void setHeatPowerField(string fileName, string fieldName, int iteration = 0, int order = 0, int meshLevel=0);
+
+ /** \fn getHeatPowerField
+ * \brief returns the heat power field
+ * \details
+ * \param [in] void
+ * \param [out] Field
+ * */
+ Field getHeatPowerField(){
+ return _heatPowerField;
}
+ //get output fields names for postprocessing or coupling
+ vector<string> getOutputFieldsNames() ;//liste tous les champs que peut fournir le code pour le postraitement
+ Field& getOutputField(const string& nameField );//Renvoie un champs pour le postraitement
+
+ Field& getOutputTemperatureField();
+ Field& getRodTemperatureField();
+
+ /** \fn setVerbose
+ * \brief Updates display options
+ * \details
+ * \param [in] bool
+ * \param [in] bool
+ * \param [out] void
+ * */
+ void setVerbose(bool verbose, bool system=false)
+ {
+ _verbose = verbose;
+ _system = system;
+ };
protected :
//Main unknown field
double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
bool _conditionNumber;//computes an estimate of the condition number
- map<string, LimitField> _limitField;
+ map<string, LimitFieldStationaryDiffusion> _limitField;
bool _onlyNeumannBC;//if true then the linear system is singular and should be solved up to a constant vector
bool _diffusionMatrixSet;
Vector _normale;
Matrix _DiffusionTensor;
Vec _deltaT, _Tk, _Tkm1, _b0;
+ Vec _Tk_seq; // Local sequential copy of the parallel vector _Tk, used for saving result files
+
double _dt_src;
//Heat transfert variables
double computeRHS(bool & stop);
double computeDiffusionMatrixFV(bool & stop);
+ double computeDiffusionMatrixFE(bool & stop);
/************ Data for FE calculation *************/
bool _FECalculation;
std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
- /*********** Functions for finite element method ***********/
- Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
- double computeDiffusionMatrixFE(bool & stop);
- int fact(int n);
- int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
- int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
-
- /********* Possibility to set a boundary field as Dirichlet boundary condition *********/
+ /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
bool _dirichletValuesSet;
+ bool _neumannValuesSet;
std::map< int, double> _dirichletBoundaryValues;
+ std::map< int, double> _neumannBoundaryValues;
+
+ /**** MPI related variables ***/
+ PetscMPIInt _mpi_size; /* size of communicator */
+ PetscMPIInt _mpi_rank; /* processor rank */
+ VecScatter _scat; /* For the distribution of a local vector */
+ int _globalNbUnknowns, _localNbUnknowns;
+ int _d_nnz, _o_nnz; /* local and "non local" numbers of non zeros corfficients */
};
#endif /* StationaryDiffusionEquation_HXX_ */
