DiffusionEquation( int dim,bool FECalculation=true,double rho=10000,double cp=300,double lambda=5, MPI_Comm comm = MPI_COMM_WORLD);
- //Gestion du calcul
+ //Gestion du calcul (ICoCo)
void initialize();
void terminate();//vide la mémoire et enregistre le résultat final
bool initTimeStep(double dt);
void setConductivity(double conductivite){
_conductivity=conductivite;
};
- void setFluidTemperatureField(Field coupledTemperatureField){
- _fluidTemperatureField=coupledTemperatureField;
- _fluidTemperatureFieldSet=true;
- };
-
void setDiffusiontensor(Matrix DiffusionTensor){
_DiffusionTensor=DiffusionTensor;
};
- void setFluidTemperature(double fluidTemperature){
- _fluidTemperature=fluidTemperature;
- }
- //get output fields 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
+ /** Set input fields to prepare the simulation or coupling **/
+ vector<string> getInputFieldsNames();
+ void setInputField(const string& nameField, Field& inputField );//supply of a required input field
- Field& getRodTemperatureField(){
- return _VV;
+ void setFluidTemperatureField(Field coupledTemperatureField){
+ _fluidTemperatureField=coupledTemperatureField;
+ _fluidTemperatureFieldSet=true;
+ };
+ void setFluidTemperature(double fluidTemperature){
+ _fluidTemperature=fluidTemperature;
}
Field& getFluidTemperatureField(){
return _fluidTemperatureField;
}
+
+ /*** 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();//Return the main unknown if present (initialize() should be called first)
+ Field& getRodTemperatureField();//Return the main unknown if present (initialize() should be called first)
/*********** Generic functions for finite element method ***********/
static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
Vector _normale;
Matrix _DiffusionTensor;
Vec _Tn, _deltaT, _Tk, _Tkm1, _b0;
+ Vec _Tn_seq; // Local sequential copy of the parallel vector _Tn, used for saving result files
+
double _dt_diffusion, _dt_src;
+ TimeScheme _timeScheme;
+ map<string, LimitFieldDiffusion> _limitField;
+
/************ Data for FE calculation *************/
int _NunknownNodes;/* number of unknown nodes for FE calculation */
std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
- TimeScheme _timeScheme;
- map<string, LimitFieldDiffusion> _limitField;
-
/********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
bool _dirichletValuesSet;
bool _neumannValuesSet;
* \param [in] bool, passage par reférence.
* \param [out] bool
* */
- virtual bool iterateTimeStep(bool &converged) = 0; //??
+ virtual bool iterateTimeStep(bool &converged) = 0;
/** \fn isStationary
* \brief vérifie la stationnairité du problème .
* */
virtual double presentTime() const;
+ /** \fn setStationaryMode
+ * \brief Perform the search of a stationary regime
+ * \details
+ * \param [in] bool
+ * \param [out]
+ * */
+ virtual void setStationaryMode(bool stationaryMode){ _stationaryMode=stationaryMode;};
+
+ /** \fn getStationaryMode
+ * \brief Tells if we are seeking a stationary regime
+ * \details
+ * \param [in]
+ * \param [out] bool
+ * */
+ virtual bool getStationaryMode(){return _stationaryMode;};
+
+ /** \fn resetTime
+ * \brief sets the current time (typically to start a new calculation)
+ * \details
+ * \param [in] double
+ * \param [out] void
+ * */
+ void resetTime (double time);
+
/*
//Coupling interface
- virtual vector<string> getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales)
- virtual Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc)
- virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale
+ virtual void getInputMEDDoubleFieldTemplate(const std::string& name, MEDDoubleField& afield) const;//Renvoie le format de champs attendu (maillage, composantes etc)
virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
- virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement
+ virtual void getOutputMEDDoubleField(const std::string& name, MEDDoubleField& afield)//Renvoie un champs pour le postraitement ou le couplage
*/
+ /** Set input fields to prepare the simulation or coupling **/
+ virtual vector<string> getInputFieldsNames()=0;
+ virtual void setInputField(const string& nameField, Field& inputField )=0;//supply of a required input field
+
+ /*! @brief (Optional) Provide the code with a scalar double data.
+ *
+ * See Problem documentation for more details on the time semantic of a scalar value.
+ *
+ * @param[in] name name of the scalar value that is given to the code.
+ * @param[in] val value passed to the code.
+ * @throws ICoCo::WrongArgument exception if the scalar name ('name' parameter) is invalid.
+ */
+ /* virtual void setInputDoubleValue(const std::string& name, const double& val); */
+
+ /*! @brief (Optional) Retrieve a scalar double value from the code.
+ *
+ * See Problem documentation for more details on the time semantic of a scalar value.
+ *
+ * @param[in] name name of the scalar value to be read from the code.
+ * @return the double value read from the code.
+ * @throws ICoCo::WrongArgument exception if the scalar name ('name' parameter) is invalid.
+ */
+ /* virtual double getOutputDoubleValue(const std::string& name) const; */
+
+ /*! @brief (Optional) Similar to setInputDoubleValue() but for an int value.
+ * @sa setInputDoubleValue()
+ */
+ /* virtual void setInputIntValue(const std::string& name, const int& val); */
+
+ /*! @brief (Optional) Similar to getOutputDoubleValue() but for an int value.
+ * @sa getOutputDoubleValue()
+ */
+ /* virtual int getOutputIntValue(const std::string& name) const; */
+
+ /*! @brief (Optional) Similar to setInputDoubleValue() but for an string value.
+ * @sa setInputDoubleValue()
+ */
+ /* virtual void setInputStringValue(const std::string& name, const std::string& val); */
+
+ /*! @brief (Optional) Similar to getOutputDoubleValue() but for an string value.
+ * @sa getOutputDoubleValue()
+ */
+ /* virtual std::string getOutputStringValue(const std::string& name) const; */
+
+ /*! @brief Return ICoCo interface major version number.
+ * @return ICoCo interface major version number (2 at present)
+ */
+ static int GetICoCoMajorVersion() { return 2; }
+
+ /*! @brief (Optional) Get MEDCoupling major version, if the code was built with MEDCoupling support.
+ *
+ * This can be used to assess compatibility between codes when coupling them.
+ *
+ * @return the MEDCoupling major version number (typically 7, 8, 9, ...)
+ */
+ virtual int getMEDCouplingMajorVersion() const{ return MEDCOUPLING_VERSION_MAJOR; };
+
+ /*! @brief (Optional) Indicate whether the code was built with a 64-bits version of MEDCoupling.
+ *
+ * Implemented if the code was built with MEDCoupling support.
+ * This can be used to assess compatibility between codes when coupling them.
+ *
+ * @return the MEDCoupling major version number
+ */
+ virtual bool isMEDCoupling64Bits() const;
+
+ /*! @brief (Optional) Get the list of input scalars accepted by the code.
+ *
+ * @return the list of scalar names that represent inputs of the code
+ * @throws ICoCo::WrongContext exception if called before initialize() or after terminate().
+ */
+ /* virtual std::vector<std::string> getInputValuesNames() const; */
+
+ /*! @brief (Optional) Get the list of output scalars that can be provided by the code.
+ *
+ * @return the list of scalar names that can be returned by the code
+ * @throws ICoCo::WrongContext exception if called before initialize() or after terminate().
+ */
+ /* virtual std::vector<std::string> getOutputValuesNames() const; */
+
Field getUnknownField() const;
//paramètres du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
- /** \fn setPresentTime
- * \brief met à jour _time (le temps courant du calcul)
- * \details
- * \param [in] double
- * \param [out] void
- * */
- void setPresentTime (double time);
-
/** \fn setMaxNbOfTimeStep
* \brief met à jour _maxNbOfTimeStep ( le nombre maximum d'itération du calcul )
* \details
void setHeatPowerField(Field heatPower){
_heatPowerField=heatPower;
_heatPowerFieldSet=true;
+ _isStationary=false;//Source term may be changed after previously reaching a stationary state
}
/** \fn setHeatPowerField
void setHeatPowerField(string fileName, string fieldName){
_heatPowerField=Field(fileName, CELLS,fieldName);
_heatPowerFieldSet=true;
+ _isStationary=false;//Source term may be changed after previously reaching a stationary state
}
/** \fn setHeatSource
- * \brief met à jour la puissance thermique ( _phi )
+ * \brief sets a constant heat power field
* \details
* \param [in] double
* \param [out] void
* */
void setHeatSource(double phi){
_heatSource=phi;
+ _isStationary=false;//Source term may be changed after previously reaching a stationary state
}
/** \fn getHeatPowerField
- * \brief renvoie le champs ?? ( _heatPowerField )
+ * \brief returns the heat power field
* \details
* \param [in] void
* \param [out] Field
return _heatPowerField;
}
- /** \fn setRodTemperatureField ??
- * \brief
- * \details
- * \param [in] Field
- * \param [out] void
- * */
- void setRodTemperatureField(Field rodTemperature){
- _rodTemperatureField=rodTemperature;
- _rodTemperatureFieldSet=true;
- }
-
- /** \fn setRodTemperature ??
- * \brief
- * \details
- * \param [in] double
- * \param [out] void
- * */
- void setRodTemperature(double rodTemp){
- _rodTemperature=rodTemp;
- }
-
- /** \fn getRodTemperatureField
- * \brief
- * \details
- * \param [in] void
- * \param [out] Field
- * */
- virtual Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
- return _rodTemperatureField;
- }
-
/** \fn setHeatTransfertCoeff
* \brief set the heat transfert coefficient for heat exchange between fluid and solid
* \details
* */
void setHeatTransfertCoeff(double heatTransfertCoeff){
_heatTransfertCoeff=heatTransfertCoeff;
+ _isStationary=false;//Source term may be changed after previously reaching a stationary state
}
/** \fn setVerbose
int _Nmailles;//number of cells
int _Nnodes;//number of nodes
int _Nfaces;//number of faces
- int _neibMaxNb;//maximum number of neighbours around a cell
+ int _neibMaxNbCells;//maximum number of neighbours around a cell
int _neibMaxNbNodes;/* maximum number of nodes around a node */
Mesh _mesh;
Field _perimeters;
bool _isStationary;
bool _initialDataSet;
bool _initializedMemory;
+ bool _stationaryMode;//ICoCo V2
bool _restartWithNewTimeScheme;
bool _restartWithNewFileName;
double _timeMax,_time;
ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
//Heat transfert variables
- Field _heatPowerField, _rodTemperatureField;
- bool _heatPowerFieldSet, _rodTemperatureFieldSet;
+ Field _heatPowerField;
+ bool _heatPowerFieldSet;
double _heatTransfertCoeff;
- double _heatSource, _rodTemperature;
+ double _heatSource;
double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this
//Display variables
string _path;//path to execution directory used for saving results
saveFormat _saveFormat;//file saving format : MED, VTK or CSV
- //MPI variables
- PetscMPIInt _size; /* size of communicator */
- PetscMPIInt _rank; /* processor rank */
-
-
+ //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 /* PROBLEMCOREFLOWS_HXX_ */
* */
void setNumericalScheme(SpaceScheme scheme, TimeScheme method=Explicit);
- //données initiales
+ /* ICoCo code coupling interface */
/*
- virtual vector<string> getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales)
virtual Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc)
- virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale
virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement
*/
-
+ /** Set input fields to prepare the simulation or coupling **/
+ vector<string> getInputFieldsNames();
+ void setInputField(const string& nameField, Field& inputField );//supply of a required input field
+
Field getConservativeField() const ;
Field getPrimitiveField() const;
double getReferencePressure() { return _Pref; };
double getReferenceTemperature() { return _Tref; };
- //get output fields for postprocessing or coupling
+ /* get output fields 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& getPressureField();
_fileName = fileName;
}
bool solveStationaryProblem();
- Field getOutputTemperatureField();
//Linear system and spectrum
void setLinearSolver(linearSolver kspType, preconditioner pcType);
void setConductivity(double conductivite){
_conductivity=conductivite;
};
+ 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){
_fluidTemperatureField=coupledTemperatureField;
_fluidTemperatureFieldSet=true;
void setFluidTemperature(double fluidTemperature){
_fluidTemperature=fluidTemperature;
}
- Field& getRodTemperatureField(){
- return _VV;
- }
Field& getFluidTemperatureField(){
return _fluidTemperatureField;
}
- void setDiffusiontensor(Matrix DiffusionTensor){
- _DiffusionTensor=DiffusionTensor;
- };
-
/** \fn setHeatPowerField
* \brief set the heat power field (variable in space)
* \details
_heatPowerFieldSet=true;
}
+ /** \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
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
std::map< int, double> _dirichletBoundaryValues;
std::map< int, double> _neumannBoundaryValues;
- //MPI variables
- PetscMPIInt _size; /* size of communicator */
- PetscMPIInt _rank; /* processor rank */
+ /**** 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_ */
_vitesseTransport=v;
};
- //get output fields for postprocessing or coupling
+ /* set input fields to prepare the simulation */
+ vector<string> getInputFieldsNames();
+ void setInputField(const string& nameField, Field& inputField );//supply of a required input field
+
+ /** \fn setRodTemperatureField
+ * \brief Set the rod temperature field
+ * \details
+ * \param [in] Field
+ * \param [out] void
+ * */
+ void setRodTemperatureField(Field rodTemperature){
+ _rodTemperatureField=rodTemperature;
+ _rodTemperatureFieldSet=true;
+ _isStationary=false;//Source term may be changed after previously reaching a stationary state
+ }
+
+ /** \fn setRodTemperature
+ * \brief Set a constant rod temperature field
+ * \details
+ * \param [in] double
+ * \param [out] void
+ * */
+ void setRodTemperature(double rodTemp){
+ _rodTemperature=rodTemp;
+ _isStationary=false;//Source term may be changed after previously reaching a stationary state
+ }
+
+ /** \fn getRodTemperatureField
+ * \brief
+ * \details
+ * \param [in] void
+ * \param [out] Field
+ * */
+ Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
+ return _rodTemperatureField;
+ }
+
+ /* get output fields 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
return _VV;
}
+ Field& getVoidFractionField(){
+ return _Alpha;
+ }
+
+ Field& getDensityField(){
+ return _Rho;
+ }
+
protected :
double computeTransportMatrix();
double computeRHS();
void updatePrimitives();
+
+ /* Postprocessing fields */
+ Field _TT, _Alpha, _Rho;//Fields of temperature, void fraction, density. Unknown field is enthalpy (_VV)
+ double _rhosatv, _rhosatl;
+ double _Tref, _href, _cpref;
+
double temperature(double h){
return _Tref+(h-_href)/_cpref;
};
return alpha*_rhosatv+(1-alpha)*_rhosatl;
};
- Field _TT, _Alpha, _Rho;//Fields of temperature and coupled temperature
- double _rhosatv, _rhosatl;
- double _Tref, _href, _cpref;
Vector _vitesseTransport, _normale;
bool _transportMatrixSet;
Vec _Hn, _deltaH, _Hk, _Hkm1, _b0;
+ Vec _Hn_seq; // Local sequential copy of the parallel vector _Hn, used for saving result files
double _dt_transport, _dt_src;
map<string, LimitFieldTransport> _limitField;
+
+ /* source terms */
+ bool _rodTemperatureFieldSet;
+ Field _rodTemperatureField;
+ double _rodTemperature;
};
#endif /* TransportEquation_HXX_ */
void DiffusionEquation::initialize()
{
- _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
-
- if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim
+ if(_mpi_rank==0)
{
- PetscPrintf(PETSC_COMM_WORLD,"Problem : dimension defined is %d but mesh dimension= %d, and space dimension is %d",_Ndim,_mesh.getMeshDimension(),_mesh.getSpaceDimension());
- *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<_mesh.getMeshDimension()<<", mesh space dim= "<<_mesh.getSpaceDimension()<<endl;
- *_runLogFile<<"DiffusionEquation::initialize: mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("DiffusionEquation::initialize: mesh has incorrect dimension");
- }
-
- if(!_initialDataSet)
- throw CdmathException("DiffusionEquation::initialize() set initial data first");
- else
- {
- PetscPrintf(PETSC_COMM_WORLD,"Initialising the diffusion of a solid temperature using ");
- *_runLogFile<<"Initialising the diffusion of a solid temperature using ";
- if(!_FECalculation)
- {
- PetscPrintf(PETSC_COMM_WORLD,"Finite volumes method\n\n");
- *_runLogFile<< "Finite volumes method"<<endl<<endl;
- }
- else
- {
- PetscPrintf(PETSC_COMM_WORLD,"Finite elements method\n\n");
- *_runLogFile<< "Finite elements method"<<endl<<endl;
- }
- }
-
- /**************** Field creation *********************/
-
- if(!_heatPowerFieldSet){
- _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
- for(int i =0; i<_VV.getNumberOfElements(); i++)
- _heatPowerField(i) = _heatSource;
- _heatPowerFieldSet=true;
- }
- if(!_fluidTemperatureFieldSet){
- _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
- for(int i =0; i<_VV.getNumberOfElements(); i++)
- _fluidTemperatureField(i) = _fluidTemperature;
- _fluidTemperatureFieldSet=true;
+ _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
+
+ if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim
+ {
+ PetscPrintf(PETSC_COMM_SELF,"Problem : dimension defined is %d but mesh dimension= %d, and space dimension is %d",_Ndim,_mesh.getMeshDimension(),_mesh.getSpaceDimension());
+ *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<_mesh.getMeshDimension()<<", mesh space dim= "<<_mesh.getSpaceDimension()<<endl;
+ *_runLogFile<<"DiffusionEquation::initialize: mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("DiffusionEquation::initialize: mesh has incorrect dimension");
+ }
+
+ if(!_initialDataSet)
+ throw CdmathException("DiffusionEquation::initialize() set initial data first");
+ else
+ {
+ PetscPrintf(PETSC_COMM_SELF,"Initialising the diffusion of a solid temperature using ");
+ *_runLogFile<<"Initialising the diffusion of a solid temperature using ";
+ if(!_FECalculation)
+ {
+ PetscPrintf(PETSC_COMM_SELF,"Finite volumes method\n\n");
+ *_runLogFile<< "Finite volumes method"<<endl<<endl;
+ }
+ else
+ {
+ PetscPrintf(PETSC_COMM_SELF,"Finite elements method\n\n");
+ *_runLogFile<< "Finite elements method"<<endl<<endl;
+ }
+ }
+
+ /**************** Field creation *********************/
+
+ if(!_heatPowerFieldSet){
+ _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
+ for(int i =0; i<_VV.getNumberOfElements(); i++)
+ _heatPowerField(i) = _heatSource;
+ _heatPowerFieldSet=true;
+ }
+ if(!_fluidTemperatureFieldSet){
+ _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
+ for(int i =0; i<_VV.getNumberOfElements(); i++)
+ _fluidTemperatureField(i) = _fluidTemperature;
+ _fluidTemperatureFieldSet=true;
+ }
+
+ /* Détection des noeuds frontière avec une condition limite de Dirichlet */
+ if(_FECalculation)
+ {
+ if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
+ PetscPrintf(PETSC_COMM_SELF,"1D Finite element method on segments\n");
+ else if(_Ndim==2)
+ {
+ if( _mesh.isTriangular() )//Mesh dim=2
+ PetscPrintf(PETSC_COMM_SELF,"2D Finite element method on triangles\n");
+ else if (_mesh.getMeshDimension()==1)//Mesh dim=1
+ PetscPrintf(PETSC_COMM_SELF,"1D Finite element method on a 2D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
+ else
+ {
+ PetscPrintf(PETSC_COMM_SELF,"Error Finite element with space dimension %, and mesh dimension %d, mesh should be either triangular either 1D network\n",_Ndim,_mesh.getMeshDimension());
+ *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
+ }
+ }
+ else if(_Ndim==3)
+ {
+ if( _mesh.isTetrahedral() )//Mesh dim=3
+ PetscPrintf(PETSC_COMM_SELF,"3D Finite element method on tetrahedra\n");
+ else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
+ PetscPrintf(PETSC_COMM_SELF,"2D Finite element method on a 3D surface : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
+ else if (_mesh.getMeshDimension()==1)//Mesh dim=1
+ PetscPrintf(PETSC_COMM_SELF,"1D Finite element method on a 3D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
+ else
+ {
+ PetscPrintf(PETSC_COMM_SELF,"Error Finite element with space dimension %d, and mesh dimension %d, mesh should be either tetrahedral, either a triangularised surface or 1D network",_Ndim,_mesh.getMeshDimension());
+ *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
+ }
+ }
+
+ _boundaryNodeIds = _mesh.getBoundaryNodeIds();
+ _NboundaryNodes=_boundaryNodeIds.size();
+
+ if(_NboundaryNodes==_Nnodes)
+ PetscPrintf(PETSC_COMM_SELF,"!!!!! Warning : all nodes are boundary nodes !!!!!");
+
+ for(int i=0; i<_NboundaryNodes; i++)
+ if(_limitField[(_mesh.getNode(_boundaryNodeIds[i])).getGroupName()].bcType==DirichletDiffusion)
+ _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
+ _NdirichletNodes=_dirichletNodeIds.size();
+ _NunknownNodes=_Nnodes - _NdirichletNodes;
+ PetscPrintf(PETSC_COMM_SELF,"Number of unknown nodes %d, Number of boundary nodes %d, Number of Dirichlet boundary nodes %d\n\n", _NunknownNodes,_NboundaryNodes, _NdirichletNodes);
+ *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
+ }
}
- /* Détection des noeuds frontière avec une condition limite de Dirichlet */
- if(_FECalculation)
- {
- if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
- PetscPrintf(PETSC_COMM_WORLD,"1D Finite element method on segments\n");
- else if(_Ndim==2)
- {
- if( _mesh.isTriangular() )//Mesh dim=2
- PetscPrintf(PETSC_COMM_WORLD,"2D Finite element method on triangles\n");
- else if (_mesh.getMeshDimension()==1)//Mesh dim=1
- PetscPrintf(PETSC_COMM_WORLD,"1D Finite element method on a 2D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
- else
- {
- PetscPrintf(PETSC_COMM_WORLD,"Error Finite element with space dimension %, and mesh dimension %d, mesh should be either triangular either 1D network\n",_Ndim,_mesh.getMeshDimension());
- *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
- }
- }
- else if(_Ndim==3)
- {
- if( _mesh.isTetrahedral() )//Mesh dim=3
- PetscPrintf(PETSC_COMM_WORLD,"3D Finite element method on tetrahedra\n");
- else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
- PetscPrintf(PETSC_COMM_WORLD,"2D Finite element method on a 3D surface : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
- else if (_mesh.getMeshDimension()==1)//Mesh dim=1
- PetscPrintf(PETSC_COMM_WORLD,"1D Finite element method on a 3D network : space dimension is %d, mesh dimension is %d\n",_Ndim,_mesh.getMeshDimension());
- else
- {
- PetscPrintf(PETSC_COMM_WORLD,"Error Finite element with space dimension %d, and mesh dimension %d, mesh should be either tetrahedral, either a triangularised surface or 1D network",_Ndim,_mesh.getMeshDimension());
- *_runLogFile<<"DiffusionEquation::initialize mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types");
- }
- }
-
- _boundaryNodeIds = _mesh.getBoundaryNodeIds();
- _NboundaryNodes=_boundaryNodeIds.size();
-
- if(_NboundaryNodes==_Nnodes)
- PetscPrintf(PETSC_COMM_WORLD,"!!!!! Warning : all nodes are boundary nodes !!!!!");
-
- for(int i=0; i<_NboundaryNodes; i++)
- if(_limitField[(_mesh.getNode(_boundaryNodeIds[i])).getGroupName()].bcType==DirichletDiffusion)
- _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
- _NdirichletNodes=_dirichletNodeIds.size();
- _NunknownNodes=_Nnodes - _NdirichletNodes;
- PetscPrintf(PETSC_COMM_WORLD,"Number of unknown nodes %d, Number of boundary nodes %d, Number of Dirichlet boundary nodes %d\n\n", _NunknownNodes,_NboundaryNodes, _NdirichletNodes);
- *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
- }
-
- //creation de la matrice
- if(!_FECalculation)
- MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
- else
- MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
-
- VecCreate(PETSC_COMM_SELF, &_Tk);
-
if(!_FECalculation)
- VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
- else
- VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
+ _globalNbUnknowns = _Nmailles*_nVar;
+ else{
+ MPI_Bcast(&_NunknownNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ _globalNbUnknowns = _NunknownNodes*_nVar;
+ }
+ /* Vectors creations */
+ VecCreate(PETSC_COMM_WORLD, &_Tk);//main unknown
+ VecSetSizes(_Tk,PETSC_DECIDE,_globalNbUnknowns);
VecSetFromOptions(_Tk);
+ VecGetLocalSize(_Tk, &_localNbUnknowns);
+
VecDuplicate(_Tk, &_Tn);
VecDuplicate(_Tk, &_Tkm1);
VecDuplicate(_Tk, &_deltaT);
VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
VecDuplicate(_Tk, &_b0);//part of the RHS that comes from the boundary conditions. Computed only once at the first time step
- for(int i =0; i<_VV.getNumberOfElements();i++)
- VecSetValue(_Tn,i,_VV(i), INSERT_VALUES);
+ if(_mpi_rank == 0)//Process 0 reads and distributes initial data
+ if(_FECalculation)
+ for(int i = 0; i<_NunknownNodes; i++)
+ {
+ int globalIndex = globalNodeIndex(i, _dirichletNodeIds);
+ VecSetValue(_Tn,i,_VV(globalIndex), INSERT_VALUES);
+ }
+ else
+ for(int i = 0; i<_Nmailles; i++)
+ VecSetValue( _Tn, i, _VV(i), INSERT_VALUES);
+ VecAssemblyBegin(_Tn);
+ VecAssemblyEnd(_Tn);
+
+ /* Matrix creation */
+ MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
+
+ /* Local sequential vector creation */
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecCreateSeq(PETSC_COMM_SELF, _globalNbUnknowns, &_Tn_seq);//For saving results on proc 0
+ VecScatterCreateToZero(_Tn,&_scat,&_Tn_seq);
//Linear solver
- KSPCreate(PETSC_COMM_SELF, &_ksp);
+ KSPCreate(PETSC_COMM_WORLD, &_ksp);
KSPSetType(_ksp, _ksptype);
// if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
_dt_src=computeRHS(stop);
+ VecAssemblyBegin(_b);
+ VecAssemblyEnd( _b);
+
+ if(_verbose or _system)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
+ VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
+ }
+
stop=false;
return min(_dt_diffusion,_dt_src);
}
{
double result;
+ MatZeroEntries(_A);
+ VecZeroEntries(_b0);
+
if(_FECalculation)
result=computeDiffusionMatrixFE(stop);
else
result=computeDiffusionMatrixFV(stop);
+ PetscPrintf(PETSC_COMM_WORLD,"Maximum diffusivity is %.2f, CFL = %.2f, Delta x = %.2f\n",_maxvp,_cfl,_minl);
+
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
+ VecAssemblyBegin(_b0);
+ VecAssemblyEnd( _b0);
+
+
//Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
//update value here if variable heat transfer coefficient
if(_timeScheme == Implicit and _heatTransfertCoeff/(_rho*_cp)>_precision)
MatShift(_A,_heatTransfertCoeff/(_rho*_cp));//Contribution from the liquit/solid heat transfer
if(_verbose or _system)
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
return result;
}
double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){
- Cell Cj;
- string nameOfGroup;
- double coeff;//Diffusion coefficients between nodes i and j
- MatZeroEntries(_A);
- VecZeroEntries(_b);
-
- Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
- std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
- std::vector< int > nodeIds(_Ndim+1);//cell node Ids
- std::vector< Node > nodes(_Ndim+1);//cell nodes
- int i_int, j_int; //index of nodes j and k considered as unknown nodes
- bool dirichletCell_treated;
-
- std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
- for (int idim=0; idim<_Ndim+1;idim++)
- values[idim][idim]=1;
-
- /* parameters for boundary treatment */
- vector< double > valuesBorder(_Ndim+1);
- Vector GradShapeFuncBorder(_Ndim+1);
-
- for (int j=0; j<_Nmailles;j++)
- {
- Cj = _mesh.getCell(j);
-
- for (int idim=0; idim<_Ndim+1;idim++){
- nodeIds[idim]=Cj.getNodeId(idim);
- nodes[idim]=_mesh.getNode(nodeIds[idim]);
- for (int jdim=0; jdim<_Ndim;jdim++)
- M(idim,jdim)=nodes[idim].getPoint()[jdim];
- M(idim,_Ndim)=1;
- }
- for (int idim=0; idim<_Ndim+1;idim++)
- GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
-
- /* Loop on the edges of the cell */
- for (int idim=0; idim<_Ndim+1;idim++)
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
- {//First node of the edge is not Dirichlet node
- i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
- dirichletCell_treated=false;
- for (int jdim=0; jdim<_Ndim+1;jdim++)
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
- {//Second node of the edge is not Dirichlet node
- j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
- MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
- }
- else if (!dirichletCell_treated)
- {//Second node of the edge is a Dirichlet node
- dirichletCell_treated=true;
- for (int kdim=0; kdim<_Ndim+1;kdim++)
- {
- std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
- if( it != _dirichletBoundaryValues.end() )
- {
- if( _dirichletValuesSet )
- valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
- else
- valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
- }
- else
- valuesBorder[kdim]=0;
- }
- GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
- coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
- VecSetValue(_b,i_int,coeff, ADD_VALUES);
- }
- }
- }
- }
+ if(_mpi_rank == 0)
+ {
+ Cell Cj;
+ string nameOfGroup;
+ double coeff;//Diffusion coefficients between nodes i and j
+
+ Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
+ std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
+ std::vector< int > nodeIds(_Ndim+1);//cell node Ids
+ std::vector< Node > nodes(_Ndim+1);//cell nodes
+ int i_int, j_int; //index of nodes j and k considered as unknown nodes
+ bool dirichletCell_treated;
+
+ std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
+ for (int idim=0; idim<_Ndim+1;idim++)
+ values[idim][idim]=1;
+
+ /* parameters for boundary treatment */
+ vector< double > valuesBorder(_Ndim+1);
+ Vector GradShapeFuncBorder(_Ndim+1);
+
+ for (int j=0; j<_Nmailles;j++)
+ {
+ Cj = _mesh.getCell(j);
+
+ for (int idim=0; idim<_Ndim+1;idim++){
+ nodeIds[idim]=Cj.getNodeId(idim);
+ nodes[idim]=_mesh.getNode(nodeIds[idim]);
+ for (int jdim=0; jdim<_Ndim;jdim++)
+ M(idim,jdim)=nodes[idim].getPoint()[jdim];
+ M(idim,_Ndim)=1;
+ }
+ for (int idim=0; idim<_Ndim+1;idim++)
+ GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
+
+ /* Loop on the edges of the cell */
+ for (int idim=0; idim<_Ndim+1;idim++)
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
+ {//First node of the edge is not Dirichlet node
+ i_int=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ dirichletCell_treated=false;
+ for (int jdim=0; jdim<_Ndim+1;jdim++)
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
+ {//Second node of the edge is not Dirichlet node
+ j_int= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
+ }
+ else if (!dirichletCell_treated)
+ {//Second node of the edge is a Dirichlet node
+ dirichletCell_treated=true;
+ for (int kdim=0; kdim<_Ndim+1;kdim++)
+ {
+ std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
+ if( it != _dirichletBoundaryValues.end() )
+ {
+ if( _dirichletValuesSet )
+ valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
+ else
+ valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
+ }
+ else
+ valuesBorder[kdim]=0;
+ }
+ GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
+ coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
+ VecSetValue(_b0,i_int,coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
+ }
+
+ //Calcul de la contribution de la condition limite de Neumann au second membre
+ if( _NdirichletNodes !=_NboundaryNodes)
+ {
+ vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
+ int NboundaryFaces=boundaryFaces.size();
+ for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
+ {
+ Face Fi = _mesh.getFace(boundaryFaces[i]);
+ for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
+ {
+ j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
+ if( _neumannValuesSet )
+ coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
+ else
+ coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
+ VecSetValue(_b0, j_int, coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
}
-
- //Calcul de la contribution de la condition limite de Neumann au second membre
- if( _NdirichletNodes !=_NboundaryNodes)
- {
- vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
- int NboundaryFaces=boundaryFaces.size();
- for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
- {
- Face Fi = _mesh.getFace(boundaryFaces[i]);
- for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
- {
- j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
- if( _neumannValuesSet )
- coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
- else
- coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
- VecSetValue(_b, j_int, coeff, ADD_VALUES);
- }
- }
- }
- }
-
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
_diffusionMatrixSet=true;
- stop=false ;
_maxvp=_diffusivity;//To do : optimise value with the mesh while respecting stability
- PetscPrintf(PETSC_COMM_SELF,"Maximum diffusivity is %.2f, CFL = %.2f, Delta x = %.2f\n",_maxvp,_cfl,_minl);
if(fabs(_maxvp)<_precision)
throw CdmathException("DiffusionEquation::computeDiffusionMatrixFE(): Error computing time step ! Maximum diffusivity is zero => division by zero");
else
}
double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){
- long nbFaces = _mesh.getNumberOfFaces();
- Face Fj;
- Cell Cell1,Cell2;
- string nameOfGroup;
- double inv_dxi, inv_dxj;
- double barycenterDistance;
- Vector normale(_Ndim);
- double dn;
- PetscInt idm, idn;
- std::vector< int > idCells;
- MatZeroEntries(_A);
- VecZeroEntries(_b0);
-
- for (int j=0; j<nbFaces;j++){
- Fj = _mesh.getFace(j);
-
- // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
- idCells = Fj.getCellsId();
- Cell1 = _mesh.getCell(idCells[0]);
- idm = idCells[0];
- for(int l=0; l<Cell1.getNumberOfFaces(); l++){
- if (j == Cell1.getFacesId()[l]){
- for (int idim = 0; idim < _Ndim; ++idim)
- normale[idim] = Cell1.getNormalVector(l,idim);
- break;
- }
- }
-
- //Compute velocity at the face Fj
- dn=(_DiffusionTensor*normale)*normale;
- if(fabs(dn)>_maxvp)
- _maxvp=fabs(dn);
-
- // compute 1/dxi = volume of Ci/area of Fj
- inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
-
- // If Fj is on the boundary
- if (Fj.getNumberOfCells()==1) {
- if(_verbose )
- {
- cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << ") "<<endl;
- }
- nameOfGroup = Fj.getGroupName();
-
- if (_limitField[nameOfGroup].bcType==NeumannDiffusion){
- VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
- }
- else if(_limitField[nameOfGroup].bcType==DirichletDiffusion){
- barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
- MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
- VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
- }
- else {
- stop=true ;
- cout<<"!!!!!!!!!!!!!!!!! Error DiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
- cout<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- cout<<"Accepted boundary conditions are NeumannDiffusion "<<NeumannDiffusion<< " and DirichletDiffusion "<<DirichletDiffusion<<endl;
- *_runLogFile<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- throw CdmathException("Boundary condition not accepted");
- }
-
- // if Fj is inside the domain
- } else if (Fj.getNumberOfCells()==2 ){
- if(_verbose )
- {
- cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
- cout << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << ") "<<endl;
+ if(_mpi_rank == 0)
+ {
+ long nbFaces = _mesh.getNumberOfFaces();
+ Face Fj;
+ Cell Cell1,Cell2;
+ string nameOfGroup;
+ double inv_dxi, inv_dxj;
+ double barycenterDistance;
+ Vector normale(_Ndim);
+ double dn;
+ PetscInt idm, idn;
+ std::vector< int > idCells;
+
+ for (int j=0; j<nbFaces;j++){
+ Fj = _mesh.getFace(j);
+
+ // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
+ idCells = Fj.getCellsId();
+ Cell1 = _mesh.getCell(idCells[0]);
+ idm = idCells[0];
+ for(int l=0; l<Cell1.getNumberOfFaces(); l++){
+ if (j == Cell1.getFacesId()[l]){
+ for (int idim = 0; idim < _Ndim; ++idim)
+ normale[idim] = Cell1.getNormalVector(l,idim);
+ break;
+ }
+ }
+
+ //Compute velocity at the face Fj
+ dn=(_DiffusionTensor*normale)*normale;
+ if(fabs(dn)>_maxvp)
+ _maxvp=fabs(dn);
+
+ // compute 1/dxi = volume of Ci/area of Fj
+ inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
+
+ // If Fj is on the boundary
+ if (Fj.getNumberOfCells()==1) {
+ if(_verbose )
+ {
+ cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << ") "<<endl;
+ }
+ nameOfGroup = Fj.getGroupName();
+
+ if (_limitField[nameOfGroup].bcType==NeumannDiffusion){
+ VecSetValue(_b0,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
+ }
+ else if(_limitField[nameOfGroup].bcType==DirichletDiffusion){
+ barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
+ MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
+ VecSetValue(_b0,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
+ }
+ else {
+ stop=true ;
+ cout<<"!!!!!!!!!!!!!!!!! Error DiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
+ cout<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
+ cout<<"Accepted boundary conditions are NeumannDiffusion "<<NeumannDiffusion<< " and DirichletDiffusion "<<DirichletDiffusion<<endl;
+ *_runLogFile<<"Boundary condition not accepted for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
+ throw CdmathException("Boundary condition not accepted");
+ }
+
+ // if Fj is inside the domain
+ } else if (Fj.getNumberOfCells()==2 ){
+ if(_verbose )
+ {
+ cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
+ cout << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << ") "<<endl;
+ }
+ Cell2 = _mesh.getCell(idCells[1]);
+ idn = idCells[1];
+ if (_Ndim > 1)
+ inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
+ else
+ inv_dxj = 1/Cell2.getMeasure();
+
+ barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
+
+ MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
}
- Cell2 = _mesh.getCell(idCells[1]);
- idn = idCells[1];
- if (_Ndim > 1)
- inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
else
- inv_dxj = 1/Cell2.getMeasure();
-
- barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
-
- MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
+ throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
}
- else
- throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
}
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- VecAssemblyBegin(_b0);
- VecAssemblyEnd(_b0);
-
_diffusionMatrixSet=true;
- stop=false ;
- cout<<"Maximum diffusivity is "<<_maxvp<< " CFL = "<<_cfl<<" Delta x = "<<_minl<<endl;
+ MPI_Bcast(&_maxvp, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);//The determination of _maxvp is optimal, unlike in the FE case
+
if(fabs(_maxvp)<_precision)
throw CdmathException("DiffusionEquation::computeDiffusionMatrixFV(): Error computing time step ! Maximum diffusivity is zero => division by zero");
else
}
double DiffusionEquation::computeRHS(bool & stop){//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeDiffusionMatrix
- VecAssemblyBegin(_b);
- double Ti;
- if(!_FECalculation)
- for (int i=0; i<_Nmailles;i++)
- //Contribution due to fluid/solide heat exchange + Contribution of the volumic heat power
- if(_timeScheme == Explicit)
- {
- VecGetValues(_Tn, 1, &i, &Ti);
- VecSetValue(_b,i,(_heatTransfertCoeff*(_fluidTemperatureField(i)-Ti)+_heatPowerField(i))/(_rho*_cp),ADD_VALUES);
- }
- else//Implicit scheme
- VecSetValue(_b,i,(_heatTransfertCoeff* _fluidTemperatureField(i) +_heatPowerField(i))/(_rho*_cp) ,ADD_VALUES);
- else
- {
- Cell Ci;
- std::vector< int > nodesId;
- for (int i=0; i<_Nmailles;i++)
- {
- Ci=_mesh.getCell(i);
- nodesId=Ci.getNodesId();
- for (int j=0; j<nodesId.size();j++)
- if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
- {
- double coeff = (_heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]))/(_rho*_cp);
- VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
- }
- }
- }
- VecAssemblyEnd(_b);
- if(_verbose or _system)
+ if(_mpi_rank == 0)
{
- PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
- VecView(_b,PETSC_VIEWER_STDOUT_SELF);
+ double Ti;
+ if(!_FECalculation)
+ for (int i=0; i<_Nmailles;i++)
+ //Contribution due to fluid/solide heat exchange + Contribution of the volumic heat power
+ if(_timeScheme == Explicit)
+ {
+ VecGetValues(_Tn, 1, &i, &Ti);
+ VecSetValue(_b,i,(_heatTransfertCoeff*(_fluidTemperatureField(i)-Ti)+_heatPowerField(i))/(_rho*_cp),ADD_VALUES);
+ }
+ else//Implicit scheme
+ VecSetValue(_b,i,(_heatTransfertCoeff* _fluidTemperatureField(i) +_heatPowerField(i))/(_rho*_cp) ,ADD_VALUES);
+ else
+ {
+ Cell Ci;
+ std::vector< int > nodesId;
+ for (int i=0; i<_Nmailles;i++)
+ {
+ Ci=_mesh.getCell(i);
+ nodesId=Ci.getNodesId();
+ for (int j=0; j<nodesId.size();j++)
+ if(!_mesh.isBorderNode(nodesId[j])) //or for better performance nodeIds[idim]>dirichletNodes.upper_bound()
+ {
+ double coeff = (_heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]))/(_rho*_cp);
+ VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
+ }
+ }
+ }
}
- stop=false ;
+
if(_heatTransfertCoeff>_precision)
return _rho*_cp/_heatTransfertCoeff;
else
bool DiffusionEquation::initTimeStep(double dt){
- if(_dt>0 and dt>0)
+ /* tricky because of code coupling */
+ if(_dt>0 and dt>0)//Previous time step was set and used
{
- //Remove the contribution from dt to prepare for new initTimeStep. The diffusion matrix is not recomputed
+ //Remove the contribution from dt to prepare for new time step. The diffusion matrix is not recomputed
if(_timeScheme == Implicit)
MatShift(_A,-1/_dt+1/dt);
//No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
}
else if(dt>0)//_dt==0, first time step
{
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
if(_timeScheme == Implicit)
- MatShift(_A,1/_dt);
+ MatShift(_A,1/dt);
}
else//dt<=0
{
throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero");
}
//At this stage _b contains _b0 + power + heat exchange
- VecAXPY(_b, 1/_dt, _Tn);
+ VecAXPY(_b, 1/dt, _Tn);
_dt = dt;
if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n");
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
}
return _dt>0;
}
void DiffusionEquation::abortTimeStep(){
- //Remove contribution od dt to the RHS
+ //Remove contribution of dt to the RHS
VecAXPY(_b, -1/_dt, _Tn);
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
- //Remove contribution od dt to the matrix
+ //Remove contribution of dt to the matrix
if(_timeScheme == Implicit)
MatShift(_A,-1/_dt);
_dt = 0;
}
else
{
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
#if PETSC_VERSION_GREATER_3_5
KSPSetOperators(_ksp, _A, _A);
{
VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
- _erreur_rel= 0;
- double Ti, dTi;
-
- for(int i=0; i<_VV.getNumberOfElements(); i++)
- {
- VecGetValues(_deltaT, 1, &i, &dTi);
- VecGetValues(_Tk, 1, &i, &Ti);
- if(_erreur_rel < fabs(dTi/Ti))
- _erreur_rel = fabs(dTi/Ti);
- }
+ VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
}
}
{
VecCopy(_Tk, _deltaT);//ici Tk=Tnp1 donc on a deltaT=Tnp1
VecAXPY(_deltaT, -1, _Tn);//On obtient deltaT=Tnp1-Tn
-
- _erreur_rel= 0;
- double Ti, dTi;
-
- for(int i=0; i<_VV.getNumberOfElements(); i++)
- {
- VecGetValues(_deltaT, 1, &i, &dTi);
- VecGetValues(_Tk, 1, &i, &Ti);
- if(_erreur_rel < fabs(dTi/Ti))
- _erreur_rel = fabs(dTi/Ti);
- }
+ VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
_isStationary =(_erreur_rel <_precision);
VecCopy(_Tk, _Tn);
VecCopy(_Tk, _Tkm1);
- if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- PetscPrintf(PETSC_COMM_WORLD,"Valeur propre locale max: %.2f\n", _maxvp);
-
_time+=_dt;
_nbTimeStep++;
+
if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
- save();
+ save();
}
void DiffusionEquation::save(){
- PetscPrintf(PETSC_COMM_SELF,"Saving numerical results\n\n");
+ PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n");
*_runLogFile<< "Saving numerical results"<< endl<<endl;
string resultFile(_path+"/DiffusionEquation");//Results
resultFile+="_";
resultFile+=_fileName;
+ if(_mpi_size>1){
+ VecScatterBegin(_scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD);
+ VecScatterEnd( _scat,_Tn,_Tn_seq,INSERT_VALUES,SCATTER_FORWARD);
+ }
+
if(_verbose or _system)
{
PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n");
- VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
+ VecView(_Tn,PETSC_VIEWER_STDOUT_WORLD);
}
- //On remplit le champ
- double Ti;
- if(!_FECalculation)
- for(int i =0; i<_Nmailles;i++)
- {
- VecGetValues(_Tn, 1, &i, &Ti);
- _VV(i)=Ti;
- }
- else
- {
- int globalIndex;
- for(int i=0; i<_NunknownNodes; i++)
- {
- VecGetValues(_Tk, 1, &i, &Ti);
- globalIndex = globalNodeIndex(i, _dirichletNodeIds);
- _VV(globalIndex)=Ti;//Assumes node numbering starts with border nodes
- }
-
- Node Ni;
- string nameOfGroup;
- for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
- {
- Ni=_mesh.getNode(_dirichletNodeIds[i]);
- nameOfGroup = Ni.getGroupName();
- _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
- }
- }
- _VV.setTime(_time,_nbTimeStep);
-
- // create mesh and component info
- if (_nbTimeStep ==0 || _restartWithNewFileName){
- if (_restartWithNewFileName)
- _restartWithNewFileName=false;
- string suppress ="rm -rf "+resultFile+"_*";
- system(suppress.c_str());//Nettoyage des précédents calculs identiques
-
- _VV.setInfoOnComponent(0,"Temperature_(K)");
- switch(_saveFormat)
- {
- case VTK :
- _VV.writeVTK(resultFile);
- break;
- case MED :
- _VV.writeMED(resultFile);
- break;
- case CSV :
- _VV.writeCSV(resultFile);
- break;
+
+ if(_mpi_rank==0){
+ //On remplit le champ
+ double Ti;
+ if(!_FECalculation)
+ for(int i =0; i<_Nmailles;i++)
+ {
+ if(_mpi_size>1)
+ VecGetValues(_Tn_seq, 1, &i, &Ti);
+ else
+ VecGetValues(_Tn , 1, &i, &Ti);
+
+ _VV(i)=Ti;
+ }
+ else
+ {
+ int globalIndex;
+ for(int i=0; i<_NunknownNodes; i++)
+ {
+ if(_mpi_size>1)
+ VecGetValues(_Tn_seq, 1, &i, &Ti);
+ else
+ VecGetValues(_Tk , 1, &i, &Ti);
+ globalIndex = globalNodeIndex(i, _dirichletNodeIds);
+ _VV(globalIndex)=Ti;
+ }
+
+ Node Ni;
+ string nameOfGroup;
+ for(int i=0; i<_NdirichletNodes; i++)//Assumes node numbering starts with border nodes
+ {
+ Ni=_mesh.getNode(_dirichletNodeIds[i]);
+ nameOfGroup = Ni.getGroupName();
+ _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
+ }
+ }
+ _VV.setTime(_time,_nbTimeStep);
+
+ // create mesh and component info
+ if (_nbTimeStep ==0 || _restartWithNewFileName){
+ if (_restartWithNewFileName)
+ _restartWithNewFileName=false;
+ string suppress ="rm -rf "+resultFile+"_*";
+ system(suppress.c_str());//Nettoyage des précédents calculs identiques
+
+ _VV.setInfoOnComponent(0,"Temperature_(K)");
+ switch(_saveFormat)
+ {
+ case VTK :
+ _VV.writeVTK(resultFile);
+ break;
+ case MED :
+ _VV.writeMED(resultFile);
+ break;
+ case CSV :
+ _VV.writeCSV(resultFile);
+ break;
+ }
}
- }
- else{ // do not create mesh
- switch(_saveFormat)
- {
- case VTK :
- _VV.writeVTK(resultFile,false);
- break;
- case MED :
- _VV.writeMED(resultFile,false);
- break;
- case CSV :
- _VV.writeCSV(resultFile);
- break;
+ else{ // do not create mesh
+ switch(_saveFormat)
+ {
+ case VTK :
+ _VV.writeVTK(resultFile,false);
+ break;
+ case MED :
+ _VV.writeMED(resultFile,false);
+ break;
+ case CSV :
+ _VV.writeCSV(resultFile);
+ break;
+ }
}
+
+ if(_isStationary)
+ {
+ resultFile+="_Stat";
+ switch(_saveFormat)
+ {
+ case VTK :
+ _VV.writeVTK(resultFile);
+ break;
+ case MED :
+ _VV.writeMED(resultFile);
+ break;
+ case CSV :
+ _VV.writeCSV(resultFile);
+ break;
+ }
+ }
}
-
- if(_isStationary)
- {
- resultFile+="_Stat";
- switch(_saveFormat)
- {
- case VTK :
- _VV.writeVTK(resultFile);
- break;
- case MED :
- _VV.writeMED(resultFile);
- break;
- case CSV :
- _VV.writeCSV(resultFile);
- break;
- }
- }
}
void DiffusionEquation::terminate(){
VecDestroy(&_b0);
VecDestroy(&_b);
MatDestroy(&_A);
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecDestroy(&_Tn_seq);
}
void
}
-vector<string> DiffusionEquation::getOutputFieldsNames()
+Field&
+DiffusionEquation::getOutputTemperatureField()
+{
+ if(!_initializedMemory)
+ throw("Computation not initialized. No temperature field available");
+ else
+ return _VV;
+}
+
+Field&
+DiffusionEquation::getRodTemperatureField()
+{
+ return getOutputTemperatureField();
+}
+
+vector<string>
+DiffusionEquation::getInputFieldsNames()
{
vector<string> result(2);
result[0]="FluidTemperature";
- result[1]="RodTemperature";
+ result[1]="HeatPower";
+
+ return result;
+}
+vector<string>
+DiffusionEquation::getOutputFieldsNames()
+{
+ vector<string> result(1);
+
+ result[0]="RodTemperature";
return result;
}
-Field& DiffusionEquation::getOutputField(const string& nameField )
+Field&
+DiffusionEquation::getOutputField(const string& nameField )
{
- if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE" )
- return getFluidTemperatureField();
- else if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
+ if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
return getRodTemperatureField();
else
{
}
}
+void
+DiffusionEquation::setInputField(const string& nameField, Field& inputField )
+{
+ if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
+ return setFluidTemperatureField( inputField) ;
+ else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
+ return setHeatPowerField( inputField );
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
+ throw CdmathException("DiffusionEquation::setInputField error : Unknown Field name");
+ }
+}
PETSC_COMM_WORLD = comm;
PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
}
- MPI_Comm_rank(PETSC_COMM_WORLD,&_rank);
- MPI_Comm_size(PETSC_COMM_WORLD,&_size);
- PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
- PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
+ MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
+ MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
+ PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
+ PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
- if(_size>1)
+ if(_mpi_size>1)
{
PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n");
PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output fields remain sequential ----\n");
- PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrixoperations are done in parallel thanks to PETSc----\n");
- PetscPrintf(PETSC_COMM_WORLD,"---- Processor %d is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n",_rank);
+ PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix operations are done in parallel thanks to PETSc----\n");
+ PetscPrintf(PETSC_COMM_WORLD,"---- Processor %d is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n",_mpi_rank);
}
/* Numerical parameter */
_precision_Newton=_precision;
_erreur_rel= 0;
_isStationary=false;
+ _stationaryMode=false;
_timeScheme=Explicit;
_wellBalancedCorrection=false;
_FECalculation=false;
/* Mesh parameters */
_Ndim=0;
_minl=0;
- _neibMaxNb=0;
+ _neibMaxNbCells=0;
+ _neibMaxNbNodes=0;
/* Memory and restart */
_initialDataSet=false;
_maxNewtonIts=50;
_PetscIts=0;//the number of iterations of the linear solver
_ksptype = (char*)&KSPGMRES;
- _pctype = (char*)&PCLU;
+ if( _mpi_size>1)
+ _pctype = (char*)&PCNONE;
+ else
+ _pctype = (char*)&PCLU;
- /* Physical parameter */
+ /* Physical parameters */
_heatPowerFieldSet=false;
_heatTransfertCoeff=0;
- _rodTemperatureFieldSet=false;
_heatSource=0;
//extracting current directory
void ProblemCoreFlows::setTimeMax(double timeMax){
_timeMax = timeMax;
}
-void ProblemCoreFlows::setPresentTime (double time)
+void ProblemCoreFlows::resetTime (double time)
{
_time=time;
}
}
void ProblemCoreFlows::setInitialField(const Field &VV)
{
-
if(_Ndim != VV.getSpaceDimension()){
*_runLogFile<<"ProblemCoreFlows::setInitialField: mesh has incorrect space dimension"<<endl;
_runLogFile->close();
_VV.setName("SOLVERLAB results");
_time=_VV.getTime();
_mesh=_VV.getMesh();
+
_initialDataSet=true;
//Mesh data
Face Fk;
if(_verbose)
- PetscPrintf(PETSC_COMM_WORLD,"Computing cell perimeters and mesh minimal diameter\n");
+ PetscPrintf(PETSC_COMM_SELF,"Processor %d Computing cell perimeters and mesh minimal diameter\n", _mpi_rank);
//Compute the maximum number of neighbours for nodes or cells
if(VV.getTypeOfField()==NODES)
_neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
else
- _neibMaxNb=_mesh.getMaxNbNeighbours(CELLS);
+ {
+ _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
- //Compute Delta x and the cell perimeters
- for (int i=0; i<_mesh.getNumberOfCells(); i++){
- Ci = _mesh.getCell(i);
- if (_Ndim > 1){
- _perimeters(i)=0;
- for (int k=0 ; k<Ci.getNumberOfFaces() ; k++){
- indexFace=Ci.getFacesId()[k];
- Fk = _mesh.getFace(indexFace);
- _minl = min(_minl,Ci.getMeasure()/Fk.getMeasure());
- _perimeters(i)+=Fk.getMeasure();
+ //Compute Delta x and the cell perimeters
+ for (int i=0; i<_mesh.getNumberOfCells(); i++){
+ Ci = _mesh.getCell(i);
+ if (_Ndim > 1){
+ _perimeters(i)=0;
+ for (int k=0 ; k<Ci.getNumberOfFaces() ; k++){
+ indexFace=Ci.getFacesId()[k];
+ Fk = _mesh.getFace(indexFace);
+ _minl = min(_minl,Ci.getMeasure()/Fk.getMeasure());
+ _perimeters(i)+=Fk.getMeasure();
+ }
+ }else{
+ _minl = min(_minl,Ci.getMeasure());
+ _perimeters(i)=Ci.getNumberOfFaces();
}
- }else{
- _minl = min(_minl,Ci.getMeasure());
- _perimeters(i)=Ci.getNumberOfFaces();
}
+ if(_verbose)
+ cout<<_perimeters<<endl;
}
- if(_verbose)
- cout<<_perimeters<<endl;
+
+ /*** MPI distribution of parameters ***/
+ MPI_Allreduce(&_initialDataSet, &_initialDataSet, 1, MPIU_BOOL, MPI_LOR, PETSC_COMM_WORLD);
+
+ int nbVoisinsMax;//Mettre en attribut ?
+ if(!_FECalculation){
+ MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ nbVoisinsMax = _neibMaxNbCells;
+ }
+ else{
+ MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ nbVoisinsMax = _neibMaxNbNodes;
+ }
+ _d_nnz = (nbVoisinsMax+1)*_nVar;
+ _o_nnz = nbVoisinsMax *_nVar;
}
+
//Function needed because swig of enum EntityType fails
void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int field_support_type)
{
bool ok; // Is the time interval successfully solved ?
_isStationary=false;//in case of a second run with a different physics or cfl
- PetscPrintf(PETSC_COMM_WORLD,"Running test case %d\n",_fileName);
+ PetscPrintf(PETSC_COMM_WORLD,"Running test case %s\n",_fileName);
_runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
*_runLogFile<< "Running test case "<< _fileName<<endl;
ProblemCoreFlows::~ProblemCoreFlows()
{
- /*
PetscBool petscInitialized;
PetscInitialized(&petscInitialized);
if(petscInitialized)
PetscFinalize();
- */
+
delete _runLogFile;
}
}
return _VV;
}
+
+bool
+ProblemCoreFlows::isMEDCoupling64Bits() const
+{
+#ifdef MEDCOUPLING_USE_64BIT_IDS
+ return true;
+#else
+ return false;
+#endif
+};
//creation de la matrice
if(_timeScheme == Implicit)
- MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
+ MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
//creation des vecteurs
VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uext);
if(_restartWithNewTimeScheme)//This is a change of time scheme during a simulation
{
if(_timeScheme == Implicit)
- MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
+ MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
else
MatDestroy(&_A);
_restartWithNewTimeScheme=false;
}
idm = idCells[0];
idn = idCells[1];
- //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNb<<endl;
+ //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNbCells<<endl;
MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
}
else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
+ cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
*_runLogFile<<"Warning: treatment of a junction node"<<endl;
if(!_sectionFieldSet)
}
else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
+ cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
*_runLogFile<<"Warning: treatment of a junction node"<<endl;
if(!_sectionFieldSet)
}
idm = idCells[0];
idn = idCells[1];
- //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNb<<endl;
+ //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNbCells<<endl;
MatSetValuesBlocked(A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
MatSetValuesBlocked(A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
}
else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
+ cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
*_runLogFile<<"Warning: treatment of a junction node"<<endl;
if(!_sectionFieldSet)
if(_timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
SNESDestroy(&_snes);
}
+
+vector<string>
+ProblemFluid::getInputFieldsNames()
+{
+ vector<string> result(1);
+
+ result[0]="HeatPower";
+ result[1]="Porosity";
+ result[2]="PressureLoss";
+ result[3]="Section";
+
+ return result;
+}
+
+void
+ProblemFluid::setInputField(const string& nameField, Field& inputField )
+{
+ if(nameField=="Porosity" || nameField=="POROSITY" || nameField=="Porosité" || nameField=="POROSITE")
+ return setPorosityField( inputField) ;
+ else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
+ return setHeatPowerField( inputField );
+ else if(nameField=="PressureLoss" || nameField=="PRESSURELOSS" || nameField=="PerteDeCharge" || nameField=="PPERTEDECHARGE")
+ return setPressureLossField( inputField) ;
+ else if(nameField=="Section" || nameField=="SECTION" )
+ return setSectionField( inputField );
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
+ throw CdmathException("SinglePhase::setInputField error : Unknown Field name");
+ }
+}
PETSC_COMM_WORLD = comm;
PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
}
- MPI_Comm_rank(PETSC_COMM_WORLD,&_rank);
- MPI_Comm_size(PETSC_COMM_WORLD,&_size);
- PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
- PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
+ MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
+ MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
+ PetscPrintf(PETSC_COMM_WORLD,"Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
+ PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
if(lambda < 0.)
_NEWTON_its=0;
int _PetscIts=0;//the number of iterations of the linear solver
_ksptype = (char*)&KSPGMRES;
- _pctype = (char*)&PCLU;
+ if( _mpi_size>1)
+ _pctype = (char*)&PCNONE;
+ else
+ _pctype = (char*)&PCLU;
+
_conditionNumber=false;
_erreur_rel= 0;
void StationaryDiffusionEquation::initialize()
{
- _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
-
- if(!_meshSet)
- throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first");
- else
- {
- cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
- *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
- if(!_FECalculation)
- {
- cout<< "Finite volumes method"<<endl<<endl;
- *_runLogFile<< "Finite volumes method"<<endl<<endl;
- }
- else
- {
- cout<< "Finite elements method"<<endl<<endl;
- *_runLogFile<< "Finite elements method"<<endl<<endl;
- }
- }
-
- /**************** Field creation *********************/
-
- if(!_heatPowerFieldSet){
- _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
- for(int i =0; i<_VV.getNumberOfElements(); i++)
- _heatPowerField(i) = _heatSource;
- _heatPowerFieldSet=true;
- }
- if(!_fluidTemperatureFieldSet){
- _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
- for(int i =0; i<_VV.getNumberOfElements(); i++)
- _fluidTemperatureField(i) = _fluidTemperature;
- _fluidTemperatureFieldSet=true;
+ if(_mpi_rank==0)
+ {
+ _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
+
+ if(!_meshSet)
+ throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first");
+ else
+ {
+ cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
+ *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using ";
+ if(!_FECalculation)
+ {
+ cout<< "Finite volumes method"<<endl<<endl;
+ *_runLogFile<< "Finite volumes method"<<endl<<endl;
+ }
+ else
+ {
+ cout<< "Finite elements method"<<endl<<endl;
+ *_runLogFile<< "Finite elements method"<<endl<<endl;
+ }
+ }
+
+ /**************** Field creation *********************/
+
+ if(!_heatPowerFieldSet){
+ _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
+ for(int i =0; i<_VV.getNumberOfElements(); i++)
+ _heatPowerField(i) = _heatSource;
+ _heatPowerFieldSet=true;
+ }
+ if(!_fluidTemperatureFieldSet){
+ _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
+ for(int i =0; i<_VV.getNumberOfElements(); i++)
+ _fluidTemperatureField(i) = _fluidTemperature;
+ _fluidTemperatureFieldSet=true;
+ }
+
+ /* Détection des noeuds frontière avec une condition limite de Dirichlet */
+ if(_FECalculation)
+ {
+ if(_NboundaryNodes==_Nnodes)
+ cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
+
+ for(int i=0; i<_NboundaryNodes; i++)
+ {
+ std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
+ if( it != _dirichletBoundaryValues.end() )
+ _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
+ else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
+ {
+ cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
+ *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
+ _runLogFile->close();
+ throw CdmathException("Missing boundary group");
+ }
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
+ {
+ cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
+ cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
+ *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
+ *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
+ _runLogFile->close();
+ throw CdmathException("Missing boundary condition");
+ }
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
+ _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
+ {
+ cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
+ cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
+ *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
+ *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
+ _runLogFile->close();
+ throw CdmathException("Wrong boundary condition");
+ }
+ }
+ _NdirichletNodes=_dirichletNodeIds.size();
+ _NunknownNodes=_Nnodes - _NdirichletNodes;
+ cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
+ *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
+ }
}
-
- /* Détection des noeuds frontière avec une condition limite de Dirichlet */
- if(_FECalculation)
- {
- if(_NboundaryNodes==_Nnodes)
- cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl<<endl;
-
- for(int i=0; i<_NboundaryNodes; i++)
- {
- std::map<int,double>::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]);
- if( it != _dirichletBoundaryValues.end() )
- _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
- else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 )
- {
- cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl;
- *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<<endl;
- _runLogFile->close();
- throw CdmathException("Missing boundary group");
- }
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion)
- {
- cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
- cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
- *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
- *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
- _runLogFile->close();
- throw CdmathException("Missing boundary condition");
- }
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion)
- _dirichletNodeIds.push_back(_boundaryNodeIds[i]);
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion)
- {
- cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
- cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl;
- *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
- *_runLogFile<< "!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion <<endl;
- _runLogFile->close();
- throw CdmathException("Wrong boundary condition");
- }
- }
- _NdirichletNodes=_dirichletNodeIds.size();
- _NunknownNodes=_Nnodes - _NdirichletNodes;
- cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
- *_runLogFile<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
- }
-
- //creation de la matrice
- if(!_FECalculation)
- MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
- else
- MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes, _NunknownNodes, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
-
- VecCreate(PETSC_COMM_SELF, &_Tk);
-
+
if(!_FECalculation)
- VecSetSizes(_Tk,PETSC_DECIDE,_Nmailles);
- else
- VecSetSizes(_Tk,PETSC_DECIDE,_NunknownNodes);
-
+ _globalNbUnknowns = _Nmailles*_nVar;
+ else{
+ MPI_Bcast(&_NunknownNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ _globalNbUnknowns = _NunknownNodes*_nVar;
+ }
+
+ /* Vectors creations */
+ VecCreate(PETSC_COMM_WORLD, &_Tk);//main unknown
+ VecSetSizes(_Tk,PETSC_DECIDE,_globalNbUnknowns);
VecSetFromOptions(_Tk);
+ VecGetLocalSize(_Tk, &_localNbUnknowns);
+
VecDuplicate(_Tk, &_Tkm1);
VecDuplicate(_Tk, &_deltaT);
- VecDuplicate(_Tk, &_b);//RHS of the linear system
+ VecDuplicate(_Tk, &_b);//RHS of the linear system: _b=Tn/dt + _b0 + puisance volumique + couplage thermique avec le fluide
+
+ /* Matrix creation */
+ MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
+
+ /* Local sequential vector creation */
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Tk_seq);//For saving results on proc 0
+ if(_mpi_size==0)
+ _Tk_seq=_Tk;
+ VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq);
//Linear solver
- KSPCreate(PETSC_COMM_SELF, &_ksp);
+ KSPCreate(PETSC_COMM_WORLD, &_ksp);
KSPSetType(_ksp, _ksptype);
// if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
{
double result;
+ MatZeroEntries(_A);
+ VecZeroEntries(_b);
+
if(_FECalculation)
result=computeDiffusionMatrixFE(stop);
else
result=computeDiffusionMatrixFV(stop);
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
+
//Contribution from the solid/fluid heat exchange with assumption of constant heat transfer coefficient
//update value here if variable heat transfer coefficient
if(_heatTransfertCoeff>_precision)
MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer
if(_verbose or _system)
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
return result;
}
double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
- Cell Cj;
- string nameOfGroup;
- double coeff;
- MatZeroEntries(_A);
- VecZeroEntries(_b);
-
- Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
- std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
- std::vector< int > nodeIds(_Ndim+1);//cell node Ids
- std::vector< Node > nodes(_Ndim+1);//cell nodes
- int i_int, j_int; //index of nodes j and k considered as unknown nodes
- bool dirichletCell_treated;
-
- std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
- for (int idim=0; idim<_Ndim+1;idim++)
- values[idim][idim]=1;
-
- /* parameters for boundary treatment */
- vector< double > valuesBorder(_Ndim+1);
- Vector GradShapeFuncBorder(_Ndim+1);
-
- for (int j=0; j<_Nmailles;j++)
- {
- Cj = _mesh.getCell(j);
-
- for (int idim=0; idim<_Ndim+1;idim++){
- nodeIds[idim]=Cj.getNodeId(idim);
- nodes[idim]=_mesh.getNode(nodeIds[idim]);
- for (int jdim=0; jdim<_Ndim;jdim++)
- M(idim,jdim)=nodes[idim].getPoint()[jdim];
- M(idim,_Ndim)=1;
- }
- for (int idim=0; idim<_Ndim+1;idim++)
- GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim);
-
- /* Loop on the edges of the cell */
- for (int idim=0; idim<_Ndim+1;idim++)
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
- {//First node of the edge is not Dirichlet node
- i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
- dirichletCell_treated=false;
- for (int jdim=0; jdim<_Ndim+1;jdim++)
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
- {//Second node of the edge is not Dirichlet node
- j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
- MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
- }
- else if (!dirichletCell_treated)
- {//Second node of the edge is a Dirichlet node
- dirichletCell_treated=true;
- for (int kdim=0; kdim<_Ndim+1;kdim++)
- {
- std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
- if( it != _dirichletBoundaryValues.end() )
- {
- if( _dirichletValuesSet )
- valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
- else
- valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
- }
- else
- valuesBorder[kdim]=0;
- }
- GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
- coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
- VecSetValue(_b,i_int,coeff, ADD_VALUES);
- }
- }
- }
- }
+ if(_mpi_rank == 0)
+ {
+ Cell Cj;
+ string nameOfGroup;
+ double coeff;
+
+ Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix
+ std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes
+ std::vector< int > nodeIds(_Ndim+1);//cell node Ids
+ std::vector< Node > nodes(_Ndim+1);//cell nodes
+ int i_int, j_int; //index of nodes j and k considered as unknown nodes
+ bool dirichletCell_treated;
+
+ std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node
+ for (int idim=0; idim<_Ndim+1;idim++)
+ values[idim][idim]=1;
+
+ /* parameters for boundary treatment */
+ vector< double > valuesBorder(_Ndim+1);
+ Vector GradShapeFuncBorder(_Ndim+1);
+
+ for (int j=0; j<_Nmailles;j++)
+ {
+ Cj = _mesh.getCell(j);
+
+ for (int idim=0; idim<_Ndim+1;idim++){
+ nodeIds[idim]=Cj.getNodeId(idim);
+ nodes[idim]=_mesh.getNode(nodeIds[idim]);
+ for (int jdim=0; jdim<_Ndim;jdim++)
+ M(idim,jdim)=nodes[idim].getPoint()[jdim];
+ M(idim,_Ndim)=1;
+ }
+ for (int idim=0; idim<_Ndim+1;idim++)
+ GradShapeFuncs[idim]=DiffusionEquation::gradientNodal(M,values[idim])/DiffusionEquation::fact(_Ndim);
+
+ /* Loop on the edges of the cell */
+ for (int idim=0; idim<_Ndim+1;idim++)
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[idim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[idim])
+ {//First node of the edge is not Dirichlet node
+ i_int=DiffusionEquation::unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ dirichletCell_treated=false;
+ for (int jdim=0; jdim<_Ndim+1;jdim++)
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[jdim])==_dirichletNodeIds.end())//!_mesh.isBorderNode(nodeIds[jdim])
+ {//Second node of the edge is not Dirichlet node
+ j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ MatSetValue(_A,i_int,j_int,(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
+ }
+ else if (!dirichletCell_treated)
+ {//Second node of the edge is a Dirichlet node
+ dirichletCell_treated=true;
+ for (int kdim=0; kdim<_Ndim+1;kdim++)
+ {
+ std::map<int,double>::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]);
+ if( it != _dirichletBoundaryValues.end() )
+ {
+ if( _dirichletValuesSet )
+ valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
+ else
+ valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T;
+ }
+ else
+ valuesBorder[kdim]=0;
+ }
+ GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
+ coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
+ VecSetValue(_b,i_int,coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
+ }
+
+ //Calcul de la contribution de la condition limite de Neumann au second membre
+ if( _NdirichletNodes !=_NboundaryNodes)
+ {
+ vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
+ int NboundaryFaces=boundaryFaces.size();
+ for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
+ {
+ Face Fi = _mesh.getFace(boundaryFaces[i]);
+ for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
+ {
+ j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
+ if( _neumannValuesSet )
+ coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
+ else
+ coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
+ VecSetValue(_b, j_int, coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
}
-
- //Calcul de la contribution de la condition limite de Neumann au second membre
- if( _NdirichletNodes !=_NboundaryNodes)
- {
- vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
- int NboundaryFaces=boundaryFaces.size();
- for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
- {
- Face Fi = _mesh.getFace(boundaryFaces[i]);
- for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
- {
- if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is a Neumann BC node (not a Dirichlet BC node)
- {
- j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
- if( _neumannValuesSet )
- coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
- else
- coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
- VecSetValue(_b, j_int, coeff, ADD_VALUES);
- }
- }
- }
- }
-
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
if(_onlyNeumannBC) //Check that the matrix is symmetric
{
}
double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){
- long nbFaces = _mesh.getNumberOfFaces();
- Face Fj;
- Cell Cell1,Cell2;
- string nameOfGroup;
- double inv_dxi, inv_dxj;
- double barycenterDistance;
- Vector normale(_Ndim);
- double dn;
- PetscInt idm, idn;
- std::vector< int > idCells;
- MatZeroEntries(_A);
- VecZeroEntries(_b);
-
- for (int j=0; j<nbFaces;j++){
- Fj = _mesh.getFace(j);
-
- // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
- idCells = Fj.getCellsId();
- Cell1 = _mesh.getCell(idCells[0]);
- idm = idCells[0];
- for(int l=0; l<Cell1.getNumberOfFaces(); l++){
- if (j == Cell1.getFacesId()[l]){
- for (int idim = 0; idim < _Ndim; ++idim)
- normale[idim] = Cell1.getNormalVector(l,idim);
- break;
- }
- }
-
- //Compute velocity at the face Fj
- dn=(_DiffusionTensor*normale)*normale;
-
- // compute 1/dxi = volume of Ci/area of Fj
- inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
-
- // If Fj is on the boundary
- if (Fj.getNumberOfCells()==1) {
- if(_verbose )
- {
- cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << ") "<<endl;
- }
-
- std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
- if( it != _dirichletBoundaryValues.end() )
- {
- barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
- MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
- VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
- }
- else
- {
- nameOfGroup = Fj.getGroupName();
-
- if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
- VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
- }
- else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
- barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
- MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
- VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
- }
- else {
- stop=true ;
- cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
- cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
- *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- _runLogFile->close();
- throw CdmathException("Boundary condition not set");
- }
- }
- // if Fj is inside the domain
- } else if (Fj.getNumberOfCells()==2 ){
- if(_verbose )
- {
- cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
- cout << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << ") "<<endl;
+ if(_mpi_rank == 0)
+ {
+ long nbFaces = _mesh.getNumberOfFaces();
+ Face Fj;
+ Cell Cell1,Cell2;
+ string nameOfGroup;
+ double inv_dxi, inv_dxj;
+ double barycenterDistance;
+ Vector normale(_Ndim);
+ double dn;
+ PetscInt idm, idn;
+ std::vector< int > idCells;
+
+ for (int j=0; j<nbFaces;j++){
+ Fj = _mesh.getFace(j);
+
+ // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
+ idCells = Fj.getCellsId();
+ Cell1 = _mesh.getCell(idCells[0]);
+ idm = idCells[0];
+ for(int l=0; l<Cell1.getNumberOfFaces(); l++){
+ if (j == Cell1.getFacesId()[l]){
+ for (int idim = 0; idim < _Ndim; ++idim)
+ normale[idim] = Cell1.getNormalVector(l,idim);
+ break;
+ }
+ }
+
+ //Compute velocity at the face Fj
+ dn=(_DiffusionTensor*normale)*normale;
+
+ // compute 1/dxi = volume of Ci/area of Fj
+ inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
+
+ // If Fj is on the boundary
+ if (Fj.getNumberOfCells()==1) {
+ if(_verbose )
+ {
+ cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << ") "<<endl;
+ }
+
+ std::map<int,double>::iterator it=_dirichletBoundaryValues.find(j);
+ if( it != _dirichletBoundaryValues.end() )
+ {
+ barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
+ MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
+ VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES);
+ }
+ else
+ {
+ nameOfGroup = Fj.getGroupName();
+
+ if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){
+ VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES);
+ }
+ else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){
+ barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
+ MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
+ VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
+ }
+ else {
+ stop=true ;
+ cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
+ cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
+ cout<<"Accepted boundary conditions are NeumannStationaryDiffusion "<<NeumannStationaryDiffusion<< " and DirichletStationaryDiffusion "<<DirichletStationaryDiffusion<<endl;
+ *_runLogFile<<"!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
+ _runLogFile->close();
+ throw CdmathException("Boundary condition not set");
+ }
+ }
+ // if Fj is inside the domain
+ } else if (Fj.getNumberOfCells()==2 ){
+ if(_verbose )
+ {
+ cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
+ cout << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << ") "<<endl;
+ }
+ Cell2 = _mesh.getCell(idCells[1]);
+ idn = idCells[1];
+ if (_Ndim > 1)
+ inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
+ else
+ inv_dxj = 1/Cell2.getMeasure();
+
+ barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
+
+ MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
+ MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
}
- Cell2 = _mesh.getCell(idCells[1]);
- idn = idCells[1];
- if (_Ndim > 1)
- inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
else
- inv_dxj = 1/Cell2.getMeasure();
-
- barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter());
-
- MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES);
- MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES);
+ {
+ *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
+ throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
+ }
}
- else
- {
- *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
- throw CdmathException("StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face");
- }
}
-
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- VecAssemblyBegin(_b);
- VecAssemblyEnd(_b);
if(_onlyNeumannBC) //Check that the matrix is symmetric
{
double StationaryDiffusionEquation::computeRHS(bool & stop)//Contribution of the PDE RHS to the linear systemm RHS (boundary conditions do contribute to the system RHS via the function computeDiffusionMatrix
{
- VecAssemblyBegin(_b);
- if(!_FECalculation)
- for (int i=0; i<_Nmailles;i++)
- VecSetValue(_b,i, _heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i) ,ADD_VALUES);
- else
- {
- Cell Ci;
- std::vector< int > nodesId;
- for (int i=0; i<_Nmailles;i++)
- {
- Ci=_mesh.getCell(i);
- nodesId=Ci.getNodesId();
- for (int j=0; j<nodesId.size();j++)
- if(!_mesh.isBorderNode(nodesId[j]))
- {
- double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
- VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
- }
- }
- }
-
+ if(_mpi_rank == 0)
+ {
+ if(!_FECalculation)
+ for (int i=0; i<_Nmailles;i++)
+ VecSetValue(_b,i, _heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i) ,ADD_VALUES);
+ else
+ {
+ Cell Ci;
+ std::vector< int > nodesId;
+ for (int i=0; i<_Nmailles;i++)
+ {
+ Ci=_mesh.getCell(i);
+ nodesId=Ci.getNodesId();
+ for (int j=0; j<nodesId.size();j++)
+ if(!_mesh.isBorderNode(nodesId[j]))
+ {
+ double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
+ VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
+ }
+ }
+ }
+ }
+ VecAssemblyBegin(_b);
VecAssemblyEnd(_b);
if(_verbose or _system)
- VecView(_b,PETSC_VIEWER_STDOUT_SELF);
+ VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
stop=false ;
return INFINITY;
bool stop;
//Only implicit scheme considered
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
#if PETSC_VERSION_GREATER_3_5
KSPSetOperators(_ksp, _A, _A);
KSPGetResidualNorm(_ksp,&residu);
if (reason!=2 and reason!=3)
{
- cout<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
- cout<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
- cout<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc.");
+ PetscPrintf(PETSC_COMM_WORLD,"!!!!!!!!!!!!! Itérations maximales %d atteintes, résidu = %1.2f, précision demandée= %1.2f",_maxPetscIts,residu,_precision);
+ PetscPrintf(PETSC_COMM_WORLD,"Solver used %s, preconditioner %s, Final number of iteration = %d",_ksptype,_pctype,_PetscIts);
*_runLogFile<<"!!!!!!!!!!!!! Erreur système linéaire : pas de convergence de Petsc."<<endl;
*_runLogFile<<"!!!!!!!!!!!!! Itérations maximales "<<_maxPetscIts<<" atteintes, résidu="<<residu<<", précision demandée= "<<_precision<<endl;
*_runLogFile<<"Solver used "<< _ksptype<<", preconditioner "<<_pctype<<", Final number of iteration= "<<_PetscIts<<endl;
else{
if( _MaxIterLinearSolver < _PetscIts)
_MaxIterLinearSolver = _PetscIts;
- cout<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<< _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"## Système linéaire résolu en %d itérations par le solveur %s et le preconditioneur %s, précision demandée = %1.2f",_PetscIts,_ksptype,_pctype,_precision);
*_runLogFile<<"## Système linéaire résolu en "<<_PetscIts<<" itérations par le solveur "<< _ksptype<<" et le preconditioneur "<<_pctype<<", précision demandée= "<<_precision<<endl<<endl;
VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk
VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1
- _erreur_rel= 0;
- double Ti, dTi;
-
VecAssemblyBegin(_Tk);
VecAssemblyEnd( _Tk);
VecAssemblyBegin(_deltaT);
VecAssemblyEnd( _deltaT);
if(_verbose)
- cout<<"Début calcul de la variation relative"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Début calcul de l'erreur maximale");
+
+ VecNorm(_deltaT,NORM_INFINITY,&_erreur_rel);
- for(int i=0; i<_NunknownNodes; i++)
- {
- VecGetValues(_deltaT, 1, &i, &dTi);
- VecGetValues(_Tk, 1, &i, &Ti);
- if(_erreur_rel < fabs(dTi/Ti))
- _erreur_rel = fabs(dTi/Ti);
- }
if(_verbose)
- cout<<"Fin calcul de la variation relative, erreur relative maximale : " << _erreur_rel <<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Fin calcul de la variation relative, erreur maximale : %1.2f", _erreur_rel );
stop=false;
converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
}
void StationaryDiffusionEquation::setMesh(const Mesh &M)
{
- if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
+ if(_mpi_rank==0)
{
- cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
- *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
- *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
+ if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim
+ {
+ cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<M.getMeshDimension()<<", and space dimension is "<<M.getSpaceDimension()<<endl;
+ *_runLogFile<< "Problem : dim = "<<_Ndim<< " but mesh dim= "<<M.getMeshDimension()<<", mesh space dim= "<<M.getSpaceDimension()<<endl;
+ *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension");
+ }
+
+ _mesh=M;
+ _Nmailles = _mesh.getNumberOfCells();
+ _Nnodes = _mesh.getNumberOfNodes();
+
+ cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
+ *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
+
+ // find maximum nb of neibourghs
+ if(!_FECalculation)
+ {
+ _VV=Field ("Temperature", CELLS, _mesh, 1);
+ _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
+ }
+ else
+ {
+ if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
+ cout<<"1D Finite element method on segments"<<endl;
+ else if(_Ndim==2)
+ {
+ if( _mesh.isTriangular() )//Mesh dim=2
+ cout<<"2D Finite element method on triangles"<<endl;
+ else if (_mesh.getMeshDimension()==1)//Mesh dim=1
+ cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
+ else
+ {
+ cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
+ *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
+ }
+ }
+ else if(_Ndim==3)
+ {
+ if( _mesh.isTetrahedral() )//Mesh dim=3
+ cout<<"3D Finite element method on tetrahedra"<<endl;
+ else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
+ cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
+ else if (_mesh.getMeshDimension()==1)//Mesh dim=1
+ cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
+ else
+ {
+ cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either tetrahedral, either a triangularised surface or 1D network"<<endl;
+ *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
+ _runLogFile->close();
+ throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
+ }
+ }
+
+ _VV=Field ("Temperature", NODES, _mesh, 1);
+
+ _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
+ _boundaryNodeIds = _mesh.getBoundaryNodeIds();
+ _NboundaryNodes=_boundaryNodeIds.size();
+ }
}
-
- _mesh=M;
- _Nmailles = _mesh.getNumberOfCells();
- _Nnodes = _mesh.getNumberOfNodes();
-
- cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;;
- *_runLogFile<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<<endl<<endl;
-
- // find maximum nb of neibourghs
- if(!_FECalculation)
- {
- _VV=Field ("Temperature", CELLS, _mesh, 1);
- _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
- }
- else
- {
- if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
- cout<<"1D Finite element method on segments"<<endl;
- else if(_Ndim==2)
- {
- if( _mesh.isTriangular() )//Mesh dim=2
- cout<<"2D Finite element method on triangles"<<endl;
- else if (_mesh.getMeshDimension()==1)//Mesh dim=1
- cout<<"1D Finite element method on a 2D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
- else
- {
- cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
- *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
- }
- }
- else if(_Ndim==3)
- {
- if( _mesh.isTetrahedral() )//Mesh dim=3
- cout<<"3D Finite element method on tetrahedra"<<endl;
- else if (_mesh.getMeshDimension()==2 and _mesh.isTriangular())//Mesh dim=2
- cout<<"2D Finite element method on a 3D surface : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
- else if (_mesh.getMeshDimension()==1)//Mesh dim=1
- cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
- else
- {
- cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either tetrahedral, either a triangularised surface or 1D network"<<endl;
- *_runLogFile<<"StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"<<endl;
- _runLogFile->close();
- throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types");
- }
- }
-
- _VV=Field ("Temperature", NODES, _mesh, 1);
-
- _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
- _boundaryNodeIds = _mesh.getBoundaryNodeIds();
- _NboundaryNodes=_boundaryNodeIds.size();
- }
+
+ /* MPI distribution parameters */
+ int nbVoisinsMax;//Mettre en attribut ?
+ if(!_FECalculation){
+ MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ nbVoisinsMax = _neibMaxNbCells;
+ }
+ else{
+ MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
+ nbVoisinsMax = _neibMaxNbNodes;
+ }
+ _d_nnz = (nbVoisinsMax+1)*_nVar;
+ _o_nnz = nbVoisinsMax *_nVar;
_meshSet=true;
}
bool stop=false; // Does the Problem want to stop (error) ?
bool converged=false; // has the newton scheme converged (end) ?
- cout<< "!!! Running test case "<< _fileName << " using ";
+ PetscPrintf(PETSC_COMM_WORLD,"!!! Running test case %s using ",_fileName);
*_runLogFile<< "!!! Running test case "<< _fileName<< " using ";
if(!_FECalculation)
{
- cout<< "Finite volumes method"<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Finite volumes method\n\n");
*_runLogFile<< "Finite volumes method"<<endl<<endl;
}
else
{
- cout<< "Finite elements method"<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Finite elements method\n\n");
*_runLogFile<< "Finite elements method"<< endl<<endl;
}
computeDiffusionMatrix( stop);//For the moment the conductivity does not depend on the temperature (linear LHS)
if (stop){
- cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing diffusion matrix, stopping calculation\n");
*_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
_runLogFile->close();
throw CdmathException("Failed computing diffusion matrix");
}
computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS)
if (stop){
- cout << "Error : failed computing right hand side, stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing right hand side, stopping calculation\n");
*_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl;
throw CdmathException("Failed computing right hand side");
}
stop = iterateNewtonStep(converged);
if (stop){
- cout << "Error : failed solving linear system, stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Error : failed solving linear system, stopping calculation\n");
*_runLogFile << "Error : failed linear system, stopping calculation"<< endl;
_runLogFile->close();
throw CdmathException("Failed solving linear system");
}
void StationaryDiffusionEquation::save(){
- cout<< "Saving numerical results"<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n");
*_runLogFile<< "Saving numerical results"<< endl<<endl;
string resultFile(_path+"/StationaryDiffusionEquation");//Results
string suppress ="rm -rf "+resultFile+"_*";
system(suppress.c_str());//Nettoyage des précédents calculs identiques
- if(_verbose or _system)
- VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
-
- double Ti;
- if(!_FECalculation)
- for(int i=0; i<_Nmailles; i++)
- {
- VecGetValues(_Tk, 1, &i, &Ti);
- _VV(i)=Ti;
- }
- else
- {
- int globalIndex;
- for(int i=0; i<_NunknownNodes; i++)
- {
- VecGetValues(_Tk, 1, &i, &Ti);
- globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
- _VV(globalIndex)=Ti;
- }
-
- Node Ni;
- string nameOfGroup;
- for(int i=0; i<_NdirichletNodes; i++)
- {
- Ni=_mesh.getNode(_dirichletNodeIds[i]);
- nameOfGroup = Ni.getGroupName();
- _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
- }
- }
+ if(_mpi_size>1){
+ VecScatterBegin(_scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
+ VecScatterEnd( _scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD);
+ }
- _VV.setInfoOnComponent(0,"Temperature_(K)");
- switch(_saveFormat)
- {
- case VTK :
- _VV.writeVTK(resultFile);
- break;
- case MED :
- _VV.writeMED(resultFile);
- break;
- case CSV :
- _VV.writeCSV(resultFile);
- break;
- }
-}
-Field
-StationaryDiffusionEquation::getOutputTemperatureField()
-{
- if(!_computationCompletedSuccessfully)
- throw("Computation not performed yet or failed. No temperature field available");
- else
- return _VV;
+ if(_verbose or _system)
+ VecView(_Tk,PETSC_VIEWER_STDOUT_WORLD);
+
+ if(_mpi_rank==0){
+ double Ti;
+ if(!_FECalculation)
+ for(int i=0; i<_Nmailles; i++)
+ {
+ if(_mpi_size>1)
+ VecGetValues(_Tk_seq, 1, &i, &Ti);
+ else
+ VecGetValues(_Tk , 1, &i, &Ti);
+ _VV(i)=Ti;
+ }
+ else
+ {
+ int globalIndex;
+ for(int i=0; i<_NunknownNodes; i++)
+ {
+ if(_mpi_size>1)
+ VecGetValues(_Tk_seq, 1, &i, &Ti);
+ else
+ VecGetValues(_Tk , 1, &i, &Ti);
+ globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
+ _VV(globalIndex)=Ti;
+ }
+
+ Node Ni;
+ string nameOfGroup;
+ for(int i=0; i<_NdirichletNodes; i++)
+ {
+ Ni=_mesh.getNode(_dirichletNodeIds[i]);
+ nameOfGroup = Ni.getGroupName();
+ _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T;
+ }
+ }
+
+ _VV.setInfoOnComponent(0,"Temperature_(K)");
+ switch(_saveFormat)
+ {
+ case VTK :
+ _VV.writeVTK(resultFile);
+ break;
+ case MED :
+ _VV.writeMED(resultFile);
+ break;
+ case CSV :
+ _VV.writeCSV(resultFile);
+ break;
+ }
+ }
}
void StationaryDiffusionEquation::terminate()
VecDestroy(&_deltaT);
VecDestroy(&_b);
MatDestroy(&_A);
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecDestroy(&_Tk_seq);
+
+ PetscBool petscInitialized;
+ PetscInitialized(&petscInitialized);
+ if(petscInitialized)
+ PetscFinalize();
+
+ delete _runLogFile;
}
void
StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues)
return my_eigenfield;
}
+
+Field&
+StationaryDiffusionEquation::getOutputTemperatureField()
+{
+ if(!_computationCompletedSuccessfully)
+ throw("Computation not performed yet or failed. No temperature field available");
+ else
+ return _VV;
+}
+
+Field&
+StationaryDiffusionEquation::getRodTemperatureField()
+{
+ return getOutputTemperatureField();
+}
+
+vector<string>
+StationaryDiffusionEquation::getInputFieldsNames()
+{
+ vector<string> result(2);
+
+ result[0]="FluidTemperature";
+ result[1]="HeatPower";
+
+ return result;
+}
+vector<string>
+StationaryDiffusionEquation::getOutputFieldsNames()
+{
+ vector<string> result(1);
+
+ result[0]="RodTemperature";
+
+ return result;
+}
+
+Field&
+StationaryDiffusionEquation::getOutputField(const string& nameField )
+{
+ if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" )
+ return getRodTemperatureField();
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
+ throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name");
+ }
+}
+
+void
+StationaryDiffusionEquation::setInputField(const string& nameField, Field& inputField )
+{
+ if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
+ return setFluidTemperatureField( inputField) ;
+ else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
+ return setHeatPowerField( inputField );
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
+ throw CdmathException("StationaryDiffusionEquation::setInputField error : Unknown Field name");
+ }
+}
_dt_transport=0;
_dt_src=0;
_transportMatrixSet=false;
+ _FECalculation=false;//Only finite volumes available
+ _rodTemperatureFieldSet=false;
+ _rodTemperature=0;
}
void TransportEquation::initialize()
{
- if(!_initialDataSet)
- throw CdmathException("TransportEquation::initialize() set initial data first");
- else if (_VV.getTypeOfField() != CELLS)
- throw CdmathException("TransportEquation::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
- else
- cout<<"Initialising the transport of a fluid enthalpy"<<endl;
- /**************** Field creation *********************/
-
- //post processing fields used only for saving results
- _TT=Field ("Temperature", CELLS, _mesh, 1);
- _Alpha=Field ("Void fraction", CELLS, _mesh, 1);
- _Rho=Field ("Mixture density", CELLS, _mesh, 1);
- //Construction des champs de post-traitement
- VecCreate(PETSC_COMM_SELF, &_Hn);
+ if(_mpi_rank==0)
+ {
+ if(!_initialDataSet)
+ throw CdmathException("TransportEquation::initialize() set initial data first");
+ else if (_VV.getTypeOfField() != CELLS)
+ throw CdmathException("TransportEquation::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
+ else
+ PetscPrintf(PETSC_COMM_SELF,"Initialising the transport of a fluid enthalpy\n");
+
+ /**************** Field creation *********************/
+
+ //post processing fields used only for saving results
+ _TT=Field ("Temperature", CELLS, _mesh, 1);
+ _Alpha=Field ("Void fraction", CELLS, _mesh, 1);
+ _Rho=Field ("Mixture density", CELLS, _mesh, 1);
+ //Construction des champs de post-traitement
+ for(int i =0; i<_Nmailles;i++){
+ _TT(i)=temperature(_VV(i));
+ _Alpha(i)=voidFraction(_VV(i));
+ _Rho(i)=density(_Alpha(i));
+ }
+ if(!_heatPowerFieldSet){
+ _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
+ for(int i =0; i<_Nmailles; i++)
+ _heatPowerField(i) = _heatSource;
+ }
+ if(!_rodTemperatureFieldSet){
+ _rodTemperatureField=Field("Rod temperature",CELLS,_mesh,1);
+ for(int i =0; i<_Nmailles; i++)
+ _rodTemperatureField(i) = _rodTemperature;
+ }
+ }
+
+ _globalNbUnknowns = _Nmailles;
+
+ /* Vectors creations */
+ VecCreate(PETSC_COMM_WORLD, &_Hn);
VecSetSizes(_Hn,PETSC_DECIDE,_Nmailles);
VecSetFromOptions(_Hn);
- for(int i =0; i<_Nmailles;i++){
- _TT(i)=temperature(_VV(i));
- _Alpha(i)=voidFraction(_VV(i));
- _Rho(i)=density(_Alpha(i));
- VecSetValue(_Hn,i,_VV(i), INSERT_VALUES);
- }
- if(!_heatPowerFieldSet){
- _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
- for(int i =0; i<_Nmailles; i++)
- _heatPowerField(i) = _heatSource;
- }
- if(!_rodTemperatureFieldSet){
- _rodTemperatureField=Field("Rod temperature",CELLS,_mesh,1);
- for(int i =0; i<_Nmailles; i++)
- _rodTemperatureField(i) = _rodTemperature;
- }
+ VecGetLocalSize(_Hn, &_localNbUnknowns);
- //creation de la matrice
- MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
VecDuplicate(_Hn, &_Hk);
VecDuplicate(_Hn, &_Hkm1);
VecDuplicate(_Hn, &_deltaH);
VecDuplicate(_Hn, &_b);//RHS of the linear system: _b=Hn/dt + _b0 + puisance
VecDuplicate(_Hn, &_b0);//part of the RHS that comes from the boundary conditions
+ if(_mpi_rank == 0)//Process 0 reads and distributes initial data
+ for(int i =0; i<_Nmailles;i++)
+ VecSetValue(_Hn,i,_VV(i), INSERT_VALUES);
+ VecAssemblyBegin(_Hn);
+ VecAssemblyEnd(_Hn);
+
+ //creation de la matrice
+ MatCreateAIJ(PETSC_COMM_WORLD, _localNbUnknowns, _localNbUnknowns, _globalNbUnknowns, _globalNbUnknowns, _d_nnz, PETSC_NULL, _o_nnz, PETSC_NULL, &_A);
+
+ /* Local sequential vector creation */
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Hn_seq);//For saving results on proc 0
+ VecScatterCreateToZero(_Hn,&_scat,&_Hn_seq);
+
//Linear solver
KSPCreate(PETSC_COMM_SELF, &_ksp);
KSPSetType(_ksp, _ksptype);
double TransportEquation::computeTimeStep(bool & stop){
if(!_transportMatrixSet)
_dt_transport=computeTransportMatrix();
+
_dt_src=computeRHS();
+ if(_verbose or _system)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
+ VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
+ }
+
stop=false;
return min(_dt_transport,_dt_src);
}
double TransportEquation::computeTransportMatrix(){
- long nbFaces = _mesh.getNumberOfFaces();
- Face Fj;
- Cell Cell1,Cell2;
- string nameOfGroup;
- double inv_dxi, inv_dxj;
- Vector normale(_Ndim);
- double un, hk;
- PetscInt idm, idn;
- std::vector< int > idCells;
MatZeroEntries(_A);
VecZeroEntries(_b0);
- for (int j=0; j<nbFaces;j++){
- Fj = _mesh.getFace(j);
-
- // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
- idCells = Fj.getCellsId();
- Cell1 = _mesh.getCell(idCells[0]);
- idm = idCells[0];
- if (_Ndim >1){
- for(int l=0; l<Cell1.getNumberOfFaces(); l++){
- if (j == Cell1.getFacesId()[l]){
- for (int idim = 0; idim < _Ndim; ++idim)
- normale[idim] = Cell1.getNormalVector(l,idim);
- break;
+
+ if(_mpi_rank == 0)
+ {
+ long nbFaces = _mesh.getNumberOfFaces();
+ Face Fj;
+ Cell Cell1,Cell2;
+ string nameOfGroup;
+ double inv_dxi, inv_dxj;
+ Vector normale(_Ndim);
+ double un, hk;
+ PetscInt idm, idn;
+ std::vector< int > idCells;
+ for (int j=0; j<nbFaces;j++){
+ Fj = _mesh.getFace(j);
+
+ // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
+ idCells = Fj.getCellsId();
+ Cell1 = _mesh.getCell(idCells[0]);
+ idm = idCells[0];
+ if (_Ndim >1){
+ for(int l=0; l<Cell1.getNumberOfFaces(); l++){
+ if (j == Cell1.getFacesId()[l]){
+ for (int idim = 0; idim < _Ndim; ++idim)
+ normale[idim] = Cell1.getNormalVector(l,idim);
+ break;
+ }
+ }
+ }else{ // _Ndim = 1 : assume that this is normal mesh : the face index increases in positive direction
+ if (Fj.getNumberOfCells()<2) {
+ if (j==0)
+ normale[0] = -1;
+ else if (j==nbFaces-1)
+ normale[0] = 1;
+ else
+ throw CdmathException("TransportEquation::ComputeTimeStep(): computation of normal vector failed");
+ } else if(Fj.getNumberOfCells()==2){
+ if (idCells[0] < idCells[1])
+ normale[0] = 1;
+ else
+ normale[0] = -1;
}
}
- }else{ // _Ndim = 1 : assume that this is normal mesh : the face index increases in positive direction
- if (Fj.getNumberOfCells()<2) {
- if (j==0)
- normale[0] = -1;
- else if (j==nbFaces-1)
- normale[0] = 1;
- else
- throw CdmathException("TransportEquation::ComputeTimeStep(): computation of normal vector failed");
- } else if(Fj.getNumberOfCells()==2){
- if (idCells[0] < idCells[1])
- normale[0] = 1;
+ //Compute velocity at the face Fj
+ un=normale*_vitesseTransport;
+ if(abs(un)>_maxvp)
+ _maxvp=abs(un);
+
+ // compute 1/dxi = volume of Ci/area of Fj
+ if (_Ndim > 1)
+ inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
+ else
+ inv_dxi = 1/Cell1.getMeasure();
+
+ // If Fj is on the boundary
+ if (Fj.getNumberOfCells()==1) {
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << ") "<<endl;
+ }
+ nameOfGroup = Fj.getGroupName();
+
+ if (_limitField[nameOfGroup].bcType==NeumannTransport || _limitField[nameOfGroup].bcType==OutletTransport ){
+ MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
+ }
+ else if(_limitField[nameOfGroup].bcType==InletTransport || _limitField[nameOfGroup].bcType==DirichletTransport){
+ if(un>0){
+ MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
+ }
+ else{
+ hk=_limitField[nameOfGroup].h;
+ VecSetValue(_b0,idm,-inv_dxi*un*hk, ADD_VALUES);
+ }
+ }
+ else {//bcType=NoneBCTransport
+ cout<<"!!!!!!!!!!!!!!! Error TransportEquation::computeTransportMatrix() !!!!!!!!!!"<<endl;
+ cout<<"!!!!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<" !!!!!!!!!!!!!! "<<endl;
+ cout<<"Accepted boundary conditions are NeumannTransport "<<NeumannTransport<< " and InletTransport "<< InletTransport <<endl;
+ throw CdmathException("Boundary condition not accepted");
+ }
+ // if Fj is inside the domain
+ } else if (Fj.getNumberOfCells()==2 ){
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
+ cout << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << normale[p] << ",";
+ cout << "). "<<endl;
+ }
+ Cell2 = _mesh.getCell(idCells[1]);
+ idn = idCells[1];
+ if (_Ndim > 1)
+ inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
else
- normale[0] = -1;
- }
- }
- //Compute velocity at the face Fj
- un=normale*_vitesseTransport;
- if(abs(un)>_maxvp)
- _maxvp=abs(un);
-
- // compute 1/dxi = volume of Ci/area of Fj
- if (_Ndim > 1)
- inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
- else
- inv_dxi = 1/Cell1.getMeasure();
-
- // If Fj is on the boundary
- if (Fj.getNumberOfCells()==1) {
- if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- {
- cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << ") "<<endl;
- }
- nameOfGroup = Fj.getGroupName();
-
- if (_limitField[nameOfGroup].bcType==NeumannTransport || _limitField[nameOfGroup].bcType==OutletTransport ){
- MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
- }
- else if(_limitField[nameOfGroup].bcType==InletTransport || _limitField[nameOfGroup].bcType==DirichletTransport){
+ inv_dxj = 1/Cell2.getMeasure();
+
if(un>0){
MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
+ MatSetValue(_A,idn,idm,-inv_dxj*un, ADD_VALUES);
}
else{
- hk=_limitField[nameOfGroup].h;
- VecSetValue(_b0,idm,-inv_dxi*un*hk, ADD_VALUES);
+ MatSetValue(_A,idm,idn,inv_dxi*un, ADD_VALUES);
+ MatSetValue(_A,idn,idn,-inv_dxj*un, ADD_VALUES);
}
}
- else {//bcType=NoneBCTransport
- cout<<"!!!!!!!!!!!!!!! Error TransportEquation::computeTransportMatrix() !!!!!!!!!!"<<endl;
- cout<<"!!!!!!!!! Boundary condition not set for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<" !!!!!!!!!!!!!! "<<endl;
- cout<<"Accepted boundary conditions are NeumannTransport "<<NeumannTransport<< " and InletTransport "<< InletTransport <<endl;
- throw CdmathException("Boundary condition not accepted");
- }
- // if Fj is inside the domain
- } else if (Fj.getNumberOfCells()==2 ){
- if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- {
- cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
- cout << " ; vecteur normal=(";
- for(int p=0; p<_Ndim; p++)
- cout << normale[p] << ",";
- cout << "). "<<endl;
- }
- Cell2 = _mesh.getCell(idCells[1]);
- idn = idCells[1];
- if (_Ndim > 1)
- inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
else
- inv_dxj = 1/Cell2.getMeasure();
-
- if(un>0){
- MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
- MatSetValue(_A,idn,idm,-inv_dxj*un, ADD_VALUES);
- }
- else{
- MatSetValue(_A,idm,idn,inv_dxi*un, ADD_VALUES);
- MatSetValue(_A,idn,idn,-inv_dxj*un, ADD_VALUES);
- }
+ throw CdmathException("TransportEquation::ComputeTimeStep(): incompatible number of cells around a face");
}
- else
- throw CdmathException("TransportEquation::ComputeTimeStep(): incompatible number of cells around a face");
}
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- VecAssemblyBegin(_b0);
- VecAssemblyEnd(_b0);
+
_transportMatrixSet=true;
+
+ MPI_Bcast(&_maxvp, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
+ PetscPrintf(PETSC_COMM_WORLD, "Maximum speed is %.2f, CFL = %.2f, Delta x = %.2f\n",_maxvp,_cfl,_minl);
+
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
+ VecAssemblyBegin(_b0);
+ VecAssemblyEnd( _b0);
+
if(abs(_maxvp)<_precision)
throw CdmathException("TransportEquation::computeTransportMatrix(): maximum eigenvalue for time step is zero");
else
}
double TransportEquation::computeRHS(){
double rhomin=INFINITY;
+
VecCopy(_b0,_b);
- VecAssemblyBegin(_b);
- if(_system)
- cout<<"second membre of transport problem"<<endl;
- for (int i=0; i<_Nmailles;i++){
- VecSetValue(_b,i,_heatTransfertCoeff*(_rodTemperatureField(i)-_TT(i))/_Rho(i),ADD_VALUES);
- VecSetValue(_b,i,_heatPowerField(i)/_Rho(i),ADD_VALUES);
+
+ if(_mpi_rank == 0)
+ {
if(_system)
- cout<<_heatPowerField(i)/_Rho(i)<<endl;
- if(_Rho(i)<rhomin)
- rhomin=_Rho(i);
+ cout<<"Second membre of transport problem"<<endl;
+ for (int i=0; i<_Nmailles;i++){
+ VecSetValue(_b,i,_heatTransfertCoeff*(_rodTemperatureField(i)-_TT(i))/_Rho(i),ADD_VALUES);
+ VecSetValue(_b,i,_heatPowerField(i)/_Rho(i),ADD_VALUES);
+ if(_system)
+ cout<<_heatPowerField(i)/_Rho(i)<<endl;
+ if(_Rho(i)<rhomin)
+ rhomin=_Rho(i);
+ }
+ }
+ VecAssemblyBegin(_b);
+ VecAssemblyEnd( _b);
+
+ if(_verbose or _system)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
+ VecView(_b,PETSC_VIEWER_STDOUT_WORLD);
}
- VecAssemblyEnd(_b);
- if(_system)
- VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
return rhomin*_cpref/_heatTransfertCoeff;
}
void TransportEquation::updatePrimitives()
{
- double hi;
- for(int i=0; i<_Nmailles; i++)
+ if(_mpi_size>1){
+ VecScatterBegin(_scat,_Hn,_Hn_seq,INSERT_VALUES,SCATTER_FORWARD);
+ VecScatterEnd( _scat,_Hn,_Hn_seq,INSERT_VALUES,SCATTER_FORWARD);
+ }
+
+ if(_verbose or _system)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n");
+ VecView(_Hn,PETSC_VIEWER_STDOUT_WORLD);
+ }
+
+ if(_mpi_rank == 0)
{
- VecGetValues(_Hk, 1, &i, &hi);
- _VV(i)=hi;
- _TT(i)=temperature(hi);
- _Alpha(i)=voidFraction(hi);
- _Rho(i)=density(_Alpha(i));
+ double hi;
+ for(int i=0; i<_Nmailles; i++)
+ {
+ if(_mpi_size>1)
+ VecGetValues(_Hn_seq, 1, &i, &hi);
+ else
+ VecGetValues(_Hn , 1, &i, &hi);
+ _VV(i)=hi;
+ _TT(i)=temperature(hi);
+ _Alpha(i)=voidFraction(hi);
+ _Rho(i)=density(_Alpha(i));
+ }
}
}
bool TransportEquation::initTimeStep(double dt){
- _dt = dt;
if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
-
- if(_timeScheme == Implicit)
- MatShift(_A,1/_dt);
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n");
+ MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
+ }
+ if(_dt>0 and dt>0)
+ {
+ //Remove the contribution from dt to prepare for new time step. The diffusion matrix is not recomputed
+ if(_timeScheme == Implicit)
+ MatShift(_A,-1/_dt+1/dt);
+ //No need to remove the contribution to the right hand side since it is recomputed from scratch at each time step
+ }
+ else if(dt>0)//_dt==0, first time step
+ {
+ if(_timeScheme == Implicit)
+ MatShift(_A,1/dt);
+ }
+ else//dt<=0
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"TransportEquation::initTimeStep %.2f = \n",dt);
+ throw CdmathException("Error TransportEquation::initTimeStep : cannot set time step to zero");
+ }
+ //At this stage _b contains _b0 + power + heat exchange
+ VecAXPY(_b, 1/dt, _Hn);
+
+ _dt=dt;
+
return _dt>0;
}
void TransportEquation::abortTimeStep(){
+ //Remove contribution of dt to the RHS
VecAXPY(_b, -1/_dt, _Hn);
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
+ //Remove contribution of dt to the matrix
if(_timeScheme == Implicit)
MatShift(_A,-1/_dt);
+
_dt = 0;
}
VecAXPY(_b, 1/_dt, _Hn);
if(_system)
{
- cout << "Vecteur Hn: " << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Vecteur Hn : \n");
VecView(_Hn, PETSC_VIEWER_STDOUT_WORLD);
cout << endl;
- cout<<"Vecteur _b "<<endl;
- VecView(_b, PETSC_VIEWER_STDOUT_SELF);
- cout << "Matrice A "<<endl;
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ PetscPrintf(PETSC_COMM_WORLD,"Vecteur _b : \n");
+ VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
+ PetscPrintf(PETSC_COMM_WORLD,"Matrice A : \n");
+ MatView(_A,PETSC_VIEWER_STDOUT_WORLD);
}
if(_timeScheme == Explicit)
MatMult(_A, _Hn, _Hk);
if(_system)
{
- cout << "Nouveau vecteur Hk: " << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Nouveau vecteur Hk: \n");
VecView(_Hk, PETSC_VIEWER_STDOUT_WORLD);
cout << endl;
}
VecAXPY(_Hk, -1, _b);
if(_system)
{
- cout << "Nouveau vecteur Hk-b: " << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Nouveau vecteur Hk-b: \n");
VecView(_Hk, PETSC_VIEWER_STDOUT_WORLD);
cout << endl;
}
VecScale(_Hk, -_dt);
if(_system)
{
- cout << "Nouveau vecteur dt*(Hk-b): " << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Nouveau vecteur dt*(Hk-b): \n");
VecView(_Hk, PETSC_VIEWER_STDOUT_WORLD);
cout << endl;
}
}
else
{
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
#if PETSC_VERSION_GREATER_3_5
KSPSetOperators(_ksp, _A, _A);
_MaxIterLinearSolver = _PetscIts;
if(_PetscIts>=_maxPetscIts)
{
- cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Systeme lineaire : pas de convergence de Petsc. Itérations maximales %d atteintes", _maxPetscIts);
converged=false;
return false;
}
{
VecCopy(_Hk, _deltaH);//ici on a deltaH=Hk
VecAXPY(_deltaH, -1, _Hkm1);//On obtient deltaH=Hk-Hkm1
- _erreur_rel= 0;
- double hi, dhi;
-
- for(int i=0; i<_Nmailles; i++)
- {
- VecGetValues(_deltaH, 1, &i, &dhi);
- VecGetValues(_Hk, 1, &i, &hi);
- if(_erreur_rel < fabs(dhi/hi))
- _erreur_rel = fabs(dhi/hi);
- }
+ VecNorm(_deltaH,NORM_INFINITY,&_erreur_rel);
+ converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
}
-
- converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
}
- updatePrimitives();
-
VecCopy(_Hk, _Hkm1);
-
return true;
}
void TransportEquation::validateTimeStep()
{
VecCopy(_Hk, _deltaH);//ici Hk=Hnp1 donc on a deltaH=Hnp1
VecAXPY(_deltaH, -1, _Hn);//On obtient deltaH=Hnp1-Hn
+ VecNorm(_deltaH,NORM_INFINITY,&_erreur_rel);
- _erreur_rel= 0;
- double hi, dhi;
-
- for(int i=0; i<_Nmailles; i++)
- {
- VecGetValues(_deltaH, 1, &i, &dhi);
- VecGetValues(_Hk, 1, &i, &hi);
- if(_erreur_rel < fabs(dhi/hi))
- _erreur_rel = fabs(dhi/hi);
- }
_isStationary =(_erreur_rel <_precision);
VecCopy(_Hk, _Hn);
VecCopy(_Hk, _Hkm1);
- if((_nbTimeStep-1)%_freqSave ==0)
- {
- cout <<"Valeur propre locale max: " << _maxvp << endl;
- //Find minimum and maximum void fractions
- double alphamin=INFINITY;
- double alphamax=-INFINITY;
- for(int i=0; i<_Nmailles; i++)
+ updatePrimitives();
+
+ if(_mpi_rank == 0)
+ if((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
{
- if(_Alpha(i)>alphamax)
- alphamax=_Alpha(i);
- if(_Alpha(i)<alphamin)
- alphamin=_Alpha(i);
+ cout <<"Valeur propre locale max: " << _maxvp << endl;
+ //Find minimum and maximum void fractions
+ double alphamin=INFINITY;
+ double alphamax=-INFINITY;
+ for(int i=0; i<_Nmailles; i++)
+ {
+ if(_Alpha(i)>alphamax)
+ alphamax=_Alpha(i);
+ if(_Alpha(i)<alphamin)
+ alphamin=_Alpha(i);
+ }
+ cout<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<endl;
}
- cout<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<endl;
- }
+
_time+=_dt;
_nbTimeStep++;
save();
-
}
void TransportEquation::terminate(){
VecDestroy(&_b0);
VecDestroy(&_b);
MatDestroy(&_A);
+ if(_mpi_size>1 && _mpi_rank == 0)
+ VecDestroy(&_Hn_seq);
}
void TransportEquation::save(){
+ PetscPrintf(PETSC_COMM_WORLD,"Saving numerical results\n\n");
+
string resultFile(_path+"/TransportEquation_");///Results
resultFile+=_fileName;
- _VV.setTime(_time,_nbTimeStep);
- _TT.setTime(_time,_nbTimeStep);
- _Alpha.setTime(_time,_nbTimeStep);
- _Rho.setTime(_time,_nbTimeStep);
-
- // create mesh and component info
- if (_nbTimeStep ==0 || _restartWithNewFileName){
- if (_restartWithNewFileName)
- _restartWithNewFileName=false;
- string suppress ="rm -rf "+resultFile+"_*";
- system(suppress.c_str());//Nettoyage des précédents calculs identiques
-
- switch(_saveFormat)
- {
- case VTK :
- _VV.writeVTK(resultFile+"Enthalpy");
- _TT.writeVTK(resultFile+"Temperature");
- _Alpha.writeVTK(resultFile+"GasFraction");
- _Rho.writeVTK(resultFile+"MixtureDensity");
- break;
- case MED :
- _VV.writeMED(resultFile+"Enthalpy");
- _TT.writeMED(resultFile+"Temperature");
- _Alpha.writeMED(resultFile+"GasFraction");
- _Rho.writeMED(resultFile+"MixtureDensity");
- break;
- case CSV :
- _VV.writeCSV(resultFile+"Enthalpy");
- _TT.writeCSV(resultFile+"Temperature");
- _Alpha.writeCSV(resultFile+"GasFraction");
- _Rho.writeCSV(resultFile+"MixtureDensity");
- break;
+ if(_mpi_rank==0){
+ _VV.setTime(_time,_nbTimeStep);
+ _TT.setTime(_time,_nbTimeStep);
+ _Alpha.setTime(_time,_nbTimeStep);
+ _Rho.setTime(_time,_nbTimeStep);
+
+ // create mesh and component info
+ if (_nbTimeStep ==0 || _restartWithNewFileName){
+ if (_restartWithNewFileName)
+ _restartWithNewFileName=false;
+ string suppress ="rm -rf "+resultFile+"_*";
+ system(suppress.c_str());//Nettoyage des précédents calculs identiques
+
+ switch(_saveFormat)
+ {
+ case VTK :
+ _VV.writeVTK(resultFile+"Enthalpy");
+ _TT.writeVTK(resultFile+"Temperature");
+ _Alpha.writeVTK(resultFile+"GasFraction");
+ _Rho.writeVTK(resultFile+"MixtureDensity");
+ break;
+ case MED :
+ _VV.writeMED(resultFile+"Enthalpy");
+ _TT.writeMED(resultFile+"Temperature");
+ _Alpha.writeMED(resultFile+"GasFraction");
+ _Rho.writeMED(resultFile+"MixtureDensity");
+ break;
+ case CSV :
+ _VV.writeCSV(resultFile+"Enthalpy");
+ _TT.writeCSV(resultFile+"Temperature");
+ _Alpha.writeCSV(resultFile+"GasFraction");
+ _Rho.writeCSV(resultFile+"MixtureDensity");
+ break;
+ }
}
- }
- // do not create mesh
- else{
- switch(_saveFormat)
- {
- case VTK :
- _VV.writeVTK(resultFile+"Enthalpy",false);
- _TT.writeVTK(resultFile+"Temperature",false);
- _Alpha.writeVTK(resultFile+"GasFraction",false);
- _Rho.writeVTK(resultFile+"MixtureDensity",false);
- break;
- case MED :
- _VV.writeMED(resultFile+"Enthalpy",false);
- _TT.writeMED(resultFile+"Temperature",false);
- _Alpha.writeMED(resultFile+"GasFraction",false);
- _Rho.writeMED(resultFile+"MixtureDensity",false);
- break;
- case CSV :
- _VV.writeCSV(resultFile+"Enthalpy");
- _TT.writeCSV(resultFile+"Temperature");
- _Alpha.writeCSV(resultFile+"GasFraction");
- _Rho.writeCSV(resultFile+"MixtureDensity");
- break;
+ // do not create mesh
+ else{
+ switch(_saveFormat)
+ {
+ case VTK :
+ _VV.writeVTK(resultFile+"Enthalpy",false);
+ _TT.writeVTK(resultFile+"Temperature",false);
+ _Alpha.writeVTK(resultFile+"GasFraction",false);
+ _Rho.writeVTK(resultFile+"MixtureDensity",false);
+ break;
+ case MED :
+ _VV.writeMED(resultFile+"Enthalpy",false);
+ _TT.writeMED(resultFile+"Temperature",false);
+ _Alpha.writeMED(resultFile+"GasFraction",false);
+ _Rho.writeMED(resultFile+"MixtureDensity",false);
+ break;
+ case CSV :
+ _VV.writeCSV(resultFile+"Enthalpy");
+ _TT.writeCSV(resultFile+"Temperature");
+ _Alpha.writeCSV(resultFile+"GasFraction");
+ _Rho.writeCSV(resultFile+"MixtureDensity");
+ break;
+ }
}
}
}
-vector<string> TransportEquation::getOutputFieldsNames()
+vector<string> TransportEquation::getInputFieldsNames()
{
vector<string> result(2);
+ result[0]="HeatPower";
+ result[1]="RodTemperature";
+
+ return result;
+}
+
+vector<string> TransportEquation::getOutputFieldsNames()
+{
+ vector<string> result(4);
+
result[0]="Enthalpy";
result[1]="FluidTemperature";
+ result[2]="VoidFraction";
+ result[3]="Density";
return result;
}
Field& TransportEquation::getOutputField(const string& nameField )
{
- if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" )
+ if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE")
return getFluidTemperatureField();
- else if(nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPY" )
+ else if(nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" )
return getEnthalpyField();
- else
+ else if(nameField=="VoidFraction" || nameField=="VOIDFRACTION" || nameField=="TauxDeVide" || nameField=="TAUXDEVIDE")
+ return getVoidFractionField();
+ else if(nameField=="Density" || nameField=="DENSITY" || nameField=="Densité" || nameField=="DENSITE" )
+ return getDensityField();
+ else
{
cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
throw CdmathException("TransportEquation::getOutputField error : Unknown Field name");
}
}
+void
+TransportEquation::setInputField(const string& nameField, Field& inputField )
+{
+ if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TemperatureCombustible" || nameField=="TEMPERATURECOMBUSTIBLE")
+ return setRodTemperatureField( inputField) ;
+ else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
+ return setHeatPowerField( inputField );
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
+ throw CdmathException("TransportEquation::setInputField error : Unknown Field name");
+ }
+}
