* \version 1.0
* \date 24 March 2015
* \brief Heat diffusion equation
- * dT/dt - \lambda\Delta T = \Phi + \lambda_{sf} (T_{fluid}-T)
+ * rho*cp*dT/dt - \lambda\Delta T = \Phi + \lambda_{sf} (T_{fluid}-T)
* Dirichlet (imposed temperature) or Neumann (imposed normal flux) boundary conditions.
* */
//============================================================================
/*! \class DiffusionEquation DiffusionEquation.hxx "DiffusionEquation.hxx"
* \brief Scalar heat equation for the Uranium rods temperature
* \details see \ref DiffusionEqPage for more details
- * dT/dt - \lambda\Delta T = \Phi + \lambda_{sf} (T_{fluid}-T)
+ * rho*cp*dT/dt - \lambda\Delta T = \Phi + \lambda_{sf} (T_{fluid}-T)
*/
#ifndef DiffusionEquation_HXX_
#define DiffusionEquation_HXX_
* \param [in] double : solid conductivity
* */
- DiffusionEquation( int dim,bool FECalculation=true,double rho=10000,double cp=300,double lambda=5);
+ DiffusionEquation( int dim,bool FECalculation=true,double rho=10000,double cp=300,double lambda=5, MPI_Comm comm = MPI_COMM_WORLD);
//Gestion du calcul
void initialize();
return _fluidTemperatureField;
}
+ /*********** Generic functions for finite element method ***********/
+ static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
+ static int fact(int n);
+ static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
+ static int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
+
protected :
double computeDiffusionMatrix(bool & stop);
double computeDiffusionMatrixFV(bool & stop);
+ double computeDiffusionMatrixFE(bool & stop);
double computeRHS(bool & stop);
Field _fluidTemperatureField;
double _dt_diffusion, _dt_src;
/************ Data for FE calculation *************/
- bool _FECalculation;
- int _neibMaxNbNodes;/* maximum number of nodes around a node */
int _NunknownNodes;/* number of unknown nodes for FE calculation */
int _NboundaryNodes;/* total number of boundary nodes */
int _NdirichletNodes;/* number of boundary nodes with Dirichlet BC for FE calculation */
std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
- /*********** Functions for finite element method ***********/
- static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
- double computeDiffusionMatrixFE(bool & stop);
- static int fact(int n);
- static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
- static int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
-
TimeScheme _timeScheme;
map<string, LimitFieldDiffusion> _limitField;
* \param [in] double : second Lamé coefficient
* */
- LinearElasticityModel( int dim, bool FECalculation=true, double rho, double lambda, double mu);
+ LinearElasticityModel( int dim, bool FECalculation=true, double rho, double lambda, double mu,MPI_Comm comm = MPI_COMM_WORLD);
void setConstantDensity(double rho) { _rho=rho; }
void setDensityField(Field densityField) { _densityField=densityField; _densityFieldSet=true;}
double computeRHS(bool & stop);
double computeStiffnessMatrixFV(bool & stop);
+ double computeStiffnessMatrixFE(bool & stop);
/************ Data for FE calculation *************/
bool _FECalculation;
std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
- /*********** Functions for finite element method ***********/
- Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
- double computeStiffnessMatrixFE(bool & stop);
- static int fact(int n);
- static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
- static int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
-
/********* Possibility to set a boundary field as Dirichlet boundary condition *********/
bool _dirichletValuesSet;
bool _neumannValuesSet;
std::map< int, double> _dirichletBoundaryValues;
std::map< int, double> _neumannBoundaryValues;
+
+ //MPI variables
+ PetscMPIInt _size; /* size of communicator */
+ PetscMPIInt _rank; /* processor rank */
};
#endif /* LinearElasticityModel_HXX_ */
{
public :
//! Constructeur par défaut
- ProblemCoreFlows();
+ ProblemCoreFlows(MPI_Comm comm = MPI_COMM_WORLD);
virtual ~ProblemCoreFlows();
// -*-*-*- Gestion du calcul (interface ICoCo) -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
* */
void setInitialField(const Field &VV);
+ /** \fn setInitialField
+ * \brief sets the initial field from a field in a med file.
+ * \details This function is added because we have not been able yet to swig properly the enum EntityType. It is replaced by an integer.
+ * \param [in] string : the file name
+ * \param [in] string : the field name
+ * \param [in] int : the time step number
+ * \param [in] int : int corresponding to the enum CELLS, NODES or FACES
+ * \param [out] void
+ * */
+ void setInitialField(string fileName, string fieldName, int timeStepNumber, int field_support_type);
+
/** \fn setInitialField
* \brief sets the initial field from a field in a med file
* \details
* \param [in] string : the file name
* \param [in] string : the field name
* \param [in] int : the time step number
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
- void setInitialField(string fileName, string fieldName, int timeStepNumber);
+ void setInitialField(string fileName, string fieldName, int timeStepNumber, EntityType typeField = CELLS);
/** \fn setInitialFieldConstant
* \brief sets a constant initial field on a mesh stored in a med file
* \details
* \param [in] string : the file name
* \param [in] vector<double> : the value in each cell
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
- void setInitialFieldConstant(string fileName, const vector<double> Vconstant);
+ void setInitialFieldConstant(string fileName, const vector<double> Vconstant, EntityType typeField = CELLS);
/** \fn setInitialFieldConstant
* \brief sets a constant initial field
* \details
* \param [in] Mesh
* \param [in] Vector
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
- void setInitialFieldConstant(const Mesh& M, const Vector Vconstant);
+ void setInitialFieldConstant(const Mesh& M, const Vector Vconstant, EntityType typeField = CELLS);
/** \fn setInitialFieldConstant
* \brief sets a constant initial field
* \details
* \param [in] Mesh
* \param [in] vector<double>
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
- void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant);
+ void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant, EntityType typeField = CELLS);
/** \fn setInitialFieldConstant
* \brief sets a constant initial field
* \param [in] double the highest value in the z direction
* \param [in] string name of the bottom boundary
* \param [in] string name of the top boundary
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
- double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
+ double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="", EntityType typeField = CELLS);
+
+ /** \fn setInitialFieldConstant
+ * \brief sets a constant initial field
+ * \details This function is added because we have not been able yet to swig roperly the enum EntityType. It is replaced by an integer.
+ * \param [in] int the space dimension
+ * \param [in] vector<double> the value in each cell
+ * \param [in] double the lowest value in the x direction
+ * \param [in] double the highest value in the x direction
+ * \param [in] string name of the left boundary
+ * \param [in] string name of the right boundary
+ * \param [in] double the lowest value in the y direction
+ * \param [in] double the highest value in the y direction
+ * \param [in] string name of the back boundary
+ * \param [in] string name of the front boundary
+ * \param [in] double the lowest value in the z direction
+ * \param [in] double the highest value in the z direction
+ * \param [in] string name of the bottom boundary
+ * \param [in] string name of the top boundary
+ * \param [in] integer corresponding to the field support enum : CELLS, NODES or FACES
+ * \param [out] void
+ * */
+ void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
+ double ymin, double ymax, int ny, string backSide, string frontSide,
+ double zmin, double zmax, int nz, string bottomSide, string topSide, int type_of_field );
/** \fn setInitialFieldStepFunction
* \brief sets a step function initial field (Riemann problem)
* \param [in] Vector
* \param [in] double position of the discontinuity on one of the three axis
* \param [in] int direction (axis carrying the discontinuity) : 0 for x, 1 for y, 2 for z
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
- void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0);
+ void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0, EntityType typeField = CELLS);
/** \fn setInitialFieldStepFunction
* \brief sets a constant initial field
* \param [in] double the highest value in the z direction
* \param [in] string name of the bottom boundary
* \param [in] string name of the top boundary
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
void setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
double xmin, double xmax,int nx, string leftSide, string rightSide,
double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
- double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
+ double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="", EntityType typeField = CELLS);
/** \fn setInitialFieldSphericalStepFunction
* \brief sets a step function initial field with value Vin inside the ball with radius Radius and Vout outside
* \param [in] Vector Vout, value outside the ball
* \param [in] double radius of the ball
* \param [in] Vector Center, coordinates of the ball center
+ * \param [in] EntityType : CELLS, NODES or FACES
* \param [out] void
* */
- void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center);
+ void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center, EntityType typeField = CELLS);
/** \fn getTime
* \brief renvoie _time (le temps courant du calcul)
_heatTransfertCoeff=heatTransfertCoeff;
}
- /** \fn setDISPLAY
- * \brief met à jour les paramètres de l'affichage
+ /** \fn setVerbose
+ * \brief Updates display options
* \details
* \param [in] bool
* \param [in] bool
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 */
+
+
};
#endif /* PROBLEMCOREFLOWS_HXX_ */
/**\fn
* \brief constructeur de la classe ProblemFluid
*/
- ProblemFluid(void);
+ ProblemFluid(MPI_Comm comm = MPI_COMM_WORLD);
//Gestion du calcul (interface ICoCo)
* \param [in] double : solid conductivity
* */
- StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1);
+ StationaryDiffusionEquation( int dim,bool FECalculation=true,double lambda=1,MPI_Comm comm = MPI_COMM_WORLD);
void setMesh(const Mesh &M);
void setFileName(string fileName){
_heatPowerFieldSet=true;
}
+ /** \fn setVerbose
+ * \brief Updates display options
+ * \details
+ * \param [in] bool
+ * \param [in] bool
+ * \param [out] void
+ * */
+ void setVerbose(bool verbose, bool system=false)
+ {
+ _verbose = verbose;
+ _system = system;
+ };
+
protected :
//Main unknown field
Field _VV;
double computeRHS(bool & stop);
double computeDiffusionMatrixFV(bool & stop);
+ double computeDiffusionMatrixFE(bool & stop);
/************ Data for FE calculation *************/
bool _FECalculation;
std::vector< int > _boundaryNodeIds;/* List of boundary nodes */
std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
- /*********** Functions for finite element method ***********/
- static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
- double computeDiffusionMatrixFE(bool & stop);
- static int fact(int n);
- static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);
- static int globalNodeIndex(int unknownIndex, std::vector< int > dirichletNodes);
-
/********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/
bool _dirichletValuesSet;
bool _neumannValuesSet;
std::map< int, double> _dirichletBoundaryValues;
std::map< int, double> _neumannBoundaryValues;
+
+ //MPI variables
+ PetscMPIInt _size; /* size of communicator */
+ PetscMPIInt _rank; /* processor rank */
};
#endif /* StationaryDiffusionEquation_HXX_ */
* \param [in] pressureMagnitude : \ref around1bar or \ref around155bars
* \param [in] vector<double> : fluid velocity (assumed constant)
* */
- TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport);
+ TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport, MPI_Comm comm = MPI_COMM_WORLD);
//Gestion du calcul
virtual void initialize();
if(! M.isSquare() )
throw CdmathException("DiffusionEquation::gradientNodal Matrix M should be square !!!");
- int Ndim = M.getNumberOfRows();
+ int Ndim = M.getNumberOfRows()-1;
vector< Matrix > matrices(Ndim);
for (int idim=0; idim<Ndim;idim++){
return result;
}
-DiffusionEquation::DiffusionEquation(int dim, bool FECalculation,double rho,double cp, double lambda){
+DiffusionEquation::DiffusionEquation(int dim, bool FECalculation,double rho,double cp, double lambda, MPI_Comm comm ):ProblemCoreFlows(comm)
+{
/* Control input value are acceptable */
if(rho<_precision or cp<_precision)
{
- std::cout<<"rho = "<<rho<<", cp = "<<cp<< ", precision = "<<_precision;
+ PetscPrintf(PETSC_COMM_WORLD,"rho = %.2f, cp = %.2f, precision = %.2f\n",rho,cp,_precision);
throw CdmathException("Error : parameters rho and cp should be strictly positive");
}
if(lambda < 0.)
{
- std::cout<<"Conductivity = "<<lambda<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Conductivity = %.2f\n",lambda);
throw CdmathException("Error : conductivity parameter lambda cannot be negative");
}
if(dim<=0)
{
- std::cout<<"Space dimension = "<<dim<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Space dimension = %.2f\n",dim);
throw CdmathException("Error : parameter dim cannot be negative");
}
- cout<<"Diffusion problem with density "<<rho<<", specific heat "<< cp<<", conductivity "<< lambda;
+ PetscPrintf(PETSC_COMM_WORLD,"Diffusion problem with density %.2f, specific heat %.2f, conductivity %.2f", rho,cp,lambda);
if(FECalculation)
- cout<<" and finite elements method"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD," and finite elements method\n");
else
- cout<<" and finite volumes method"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD," and finite volumes method\n");
_FECalculation=FECalculation;
/* Finite element data */
- _neibMaxNbNodes=0;
_boundaryNodeIds=std::vector< int >(0);
_dirichletNodeIds=std::vector< int >(0);
_NboundaryNodes=0;
_dt_src=0;
_diffusionMatrixSet=false;
- _fileName = "CoreFlowsDiffusionProblem";
+ _fileName = "SolverlabDiffusionProblem";
_runLogFile=new ofstream;
- /* Default diffusion tensor is identity matrix */
+ /* Default diffusion tensor is diagonal */
_DiffusionTensor=Matrix(_Ndim);
for(int idim=0;idim<_Ndim;idim++)
- _DiffusionTensor(idim,idim)=1;
+ _DiffusionTensor(idim,idim)=_diffusivity;
}
void DiffusionEquation::initialize()
if(_Ndim != _mesh.getSpaceDimension() or _Ndim!=_mesh.getMeshDimension())//for the moment we must have space dim=mesh dim
{
- cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<<_mesh.getMeshDimension()<<", and space dimension is "<<_mesh.getSpaceDimension()<<endl;
+ 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() set initial data first");
else
{
- cout<<"Initialising the diffusion of a solid temperature using ";
+ PetscPrintf(PETSC_COMM_WORLD,"Initialising the diffusion of a solid temperature using ");
*_runLogFile<<"Initialising the diffusion of a solid temperature 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;
}
}
if(_FECalculation)
{
if(_Ndim==1 )//The 1D cdmath mesh is necessarily made of segments
- cout<<"1D Finite element method on segments"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"1D Finite element method on segments\n");
else if(_Ndim==2)
{
if( _mesh.isTriangular() )//Mesh dim=2
- cout<<"2D Finite element method on triangles"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"2D Finite element method on triangles\n");
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;
+ 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
{
- cout<<"Error Finite element with Space dimension "<<_Ndim<< ", and mesh dimension "<<_mesh.getMeshDimension()<< ", mesh should be either triangular either 1D network"<<endl;
+ 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
- cout<<"3D Finite element method on tetrahedra"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"3D Finite element method on tetrahedra\n");
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;
+ 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
- cout<<"1D Finite element method on a 3D network : space dimension is "<<_Ndim<< ", mesh dimension is "<<_mesh.getMeshDimension()<<endl;
+ 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
{
- 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;
+ 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");
_NboundaryNodes=_boundaryNodeIds.size();
if(_NboundaryNodes==_Nnodes)
- cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"<<endl;
+ 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;
- cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <<endl<<endl;
+ 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;
}
}
double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){
+
Cell Cj;
string nameOfGroup;
double coeff;//Diffusion coefficients between nodes i and j
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,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
+ 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
valuesBorder[kdim]=0;
}
GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
- coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
+ coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
VecSetValue(_b,i_int,coeff, ADD_VALUES);
}
}
int NboundaryFaces=boundaryFaces.size();
for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
{
- Face Fi = _mesh.getFace(i);
+ 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 an Neumann BC node (not a Dirichlet BC node)
+ 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
+ else
coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalFlux;
VecSetValue(_b, j_int, coeff, ADD_VALUES);
}
_diffusionMatrixSet=true;
stop=false ;
- cout<<"Maximum diffusivity is "<<_maxvp<< " CFL = "<<_cfl<<" Delta x = "<<_minl<<endl;
+ _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
}
//Compute velocity at the face Fj
- dn=_conductivity*(_DiffusionTensor*normale)*normale;
+ dn=(_DiffusionTensor*normale)*normale;
if(fabs(dn)>_maxvp)
_maxvp=fabs(dn);
return _cfl*_minl*_minl/_maxvp;
}
-double DiffusionEquation::computeRHS(bool & stop){
+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++)
- {
- VecSetValue(_b,i,_heatPowerField(i)/(_rho*_cp),ADD_VALUES);//Contribution of the volumic heat power
- //Contribution due to fluid/solide heat exchange
+ //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/(_rho*_cp)*(_fluidTemperatureField(i)-Ti),ADD_VALUES);
+ VecSetValue(_b,i,(_heatTransfertCoeff*(_fluidTemperatureField(i)-Ti)+_heatPowerField(i))/(_rho*_cp),ADD_VALUES);
}
else//Implicit scheme
- VecSetValue(_b,i,_heatTransfertCoeff/(_rho*_cp)* _fluidTemperatureField(i) ,ADD_VALUES);
- }
+ VecSetValue(_b,i,(_heatTransfertCoeff* _fluidTemperatureField(i) +_heatPowerField(i))/(_rho*_cp) ,ADD_VALUES);
else
{
Cell Ci;
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]);
+ 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)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Right hand side of the linear system\n");
VecView(_b,PETSC_VIEWER_STDOUT_SELF);
-
+ }
stop=false ;
if(_heatTransfertCoeff>_precision)
return _rho*_cp/_heatTransfertCoeff;
}
else//dt<=0
{
- cout<<"DiffusionEquation::initTimeStep dt= "<<dt<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"DiffusionEquation::initTimeStep %.2f = \n",dt);
throw CdmathException("Error DiffusionEquation::initTimeStep : cannot set time step to zero");
}
//At this stage _b contains _b0 + power + heat exchange
_dt = dt;
if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Matrix of the linear system\n");
MatView(_A,PETSC_VIEWER_STDOUT_SELF);
-
+ }
+
return _dt>0;
}
_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 \n",_maxPetscIts);
*_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
converged=false;
return false;
VecCopy(_Tk, _Tkm1);
if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
- cout <<"Valeur propre locale max: " << _maxvp << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Valeur propre locale max: %.2f\n", _maxvp);
_time+=_dt;
_nbTimeStep++;
}
void DiffusionEquation::save(){
- cout<< "Saving numerical results"<<endl<<endl;
+ PetscPrintf(PETSC_COMM_SELF,"Saving numerical results\n\n");
*_runLogFile<< "Saving numerical results"<< endl<<endl;
string resultFile(_path+"/DiffusionEquation");//Results
resultFile+=_fileName;
if(_verbose or _system)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Unknown of the linear system :\n");
VecView(_Tk,PETSC_VIEWER_STDOUT_SELF);
-
+ }
//On remplit le champ
double Ti;
if(!_FECalculation)
using namespace std;
-int LinearElasticityModel::fact(int n)
-{
- return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
-}
-int LinearElasticityModel::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes)
-{//assumes Dirichlet node numbering is strictly increasing
- int j=0;//indice de parcours des noeuds frontière avec CL Dirichlet
- int boundarySize=dirichletNodes.size();
- while(j<boundarySize and dirichletNodes[j]<globalIndex)
- j++;
- if(j==boundarySize)
- return globalIndex-boundarySize;
- else if (dirichletNodes[j]>globalIndex)
- return globalIndex-j;
- else
- throw CdmathException("LinearElasticityModel::unknownNodeIndex : Error : node is a Dirichlet boundary node");
-}
-
-int LinearElasticityModel::globalNodeIndex(int unknownNodeIndex, std::vector< int > dirichletNodes)
-{//assumes Dirichlet boundary node numbering is strictly increasing
- int boundarySize=dirichletNodes.size();
- /* trivial case where all boundary nodes are Neumann BC */
- if(boundarySize==0)
- return unknownNodeIndex;
-
- double unknownNodeMax=-1;//max unknown node number in the interval between jth and (j+1)th Dirichlet boundary nodes
- int j=0;//indice de parcours des noeuds frontière
- //On cherche l'intervale [j,j+1] qui contient le noeud de numéro interieur unknownNodeIndex
- while(j+1<boundarySize and unknownNodeMax<unknownNodeIndex)
- {
- unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
- j++;
- }
-
- if(j+1==boundarySize)
- return unknownNodeIndex+boundarySize;
- else //unknownNodeMax>=unknownNodeIndex) hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j]
- return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1;
-}
-
-LinearElasticityModel::LinearElasticityModel(int dim, bool FECalculation, double rho, double lambda, double mu){
+LinearElasticityModel::LinearElasticityModel(int dim, bool FECalculation, double rho, double lambda, double mu,MPI_Comm comm){
+ /* Initialisation of PETSC */
+ //check if PETSC is already initialised
PetscBool petscInitialized;
PetscInitialized(&petscInitialized);
if(!petscInitialized)
- PetscInitialize(NULL,NULL,0,0);
+ {//check if MPI is already initialised
+ int mpiInitialized;
+ MPI_Initialized(&mpiInitialized);
+ if(mpiInitialized)
+ 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.
+ PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
if(lambda < 0.)
{
}
-Vector LinearElasticityModel::gradientNodal(Matrix M, vector< double > values){
- vector< Matrix > matrices(_Ndim);
-
- for (int idim=0; idim<_Ndim;idim++){
- matrices[idim]=M.deepCopy();
- for (int jdim=0; jdim<_Ndim+1;jdim++)
- matrices[idim](jdim,idim) = values[jdim] ;
- }
-
- Vector result(_Ndim);
- for (int idim=0; idim<_Ndim;idim++)
- result[idim] = matrices[idim].determinant();
-
- return result;
-}
-
double LinearElasticityModel::computeStiffnessMatrix(bool & stop)
{
double result;
M(idim,_Ndim)=1;
}
for (int idim=0; idim<_Ndim+1;idim++)
- GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
+ 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=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ 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= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ j_int= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
}
else if (!dirichletCell_treated)
else
valuesBorder[kdim]=Vector(_Ndim);
}
- GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
+ GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
double coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
VecSetValue(_b,i_int,coeff, ADD_VALUES);
}
int NboundaryFaces=boundaryFaces.size();
for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
{
- Face Fi = _mesh.getFace(i);
+ 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 an Neumann BC node (not a Dirichlet BC node)
{
- j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
+ j_int=DiffusionEquation::unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
if( _neumannValuesSet )
coeff =Fi.getMeasure()/(_Ndim+1)*_neumannBoundaryValues[Fi.getNodeId(j)];
else
for (int j=0; j<nodesId.size();j++)
if(!_mesh.isBorderNode(nodesId[j]))
for (int k=0; k<_Ndim; k++)
- VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds)*nVar+k, _gravity(k)*_densityField(j)*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
+ VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds)*nVar+k, _gravity(k)*_densityField(j)*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
}
}
int globalIndex;
for(int i=0; i<_NunknownNodes; i++)
{
- globalIndex = globalNodeIndex(i, _dirichletNodeIds);
+ globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
for(int j=0; j<_nVar; j++)
{
int k=i*_nVar+j;
using namespace std;
-ProblemCoreFlows::ProblemCoreFlows()
+ProblemCoreFlows::ProblemCoreFlows(MPI_Comm comm)
{
+ /* Initialisation of PETSC */
+ //check if PETSC is already initialised
PetscBool petscInitialized;
PetscInitialized(&petscInitialized);
if(!petscInitialized)
- PetscInitialize(NULL,NULL,0,0);
+ {//check if MPI is already initialised
+ int mpiInitialized;
+ MPI_Initialized(&mpiInitialized);
+ if(mpiInitialized)
+ 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.
+ PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
+
+ if(_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);
+ }
+
+ /* Numerical parameter */
_dt = 0;
_time = 0;
_nbTimeStep=0;
_precision_Newton=_precision;
_erreur_rel= 0;
_isStationary=false;
+ _timeScheme=Explicit;
+ _wellBalancedCorrection=false;
+ _FECalculation=false;
+ _nonLinearSolver = Newton_SOLVERLAB;
+ _conditionNumber=false;
+ _maxvp=0;
+ _runLogFile=new ofstream;
+
+ /* Monitoring of simulation */
+ _restartWithNewTimeScheme=false;
+ _restartWithNewFileName=false;
+ _fileName = "myCoreFlowsProblem";
+ _freqSave = 1;
+ _verbose = false;
+ _system = false;
+
+ /* Mesh parameters */
_Ndim=0;
_minl=0;
_neibMaxNb=0;
- _fileName = "myCoreFlowsProblem";
- _freqSave = 1;
+
+ /* Memory and restart */
_initialDataSet=false;
_initializedMemory=false;
- _restartWithNewTimeScheme=false;
- _restartWithNewFileName=false;
- _timeScheme=Explicit;
- _wellBalancedCorrection=false;
- _FECalculation=false;
- _maxPetscIts=50;
+
+ /* PETSc and linear solver parameters */
_MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
- _maxNewtonIts=50;
_NEWTON_its=0;
- int _PetscIts=0;//the number of iterations of the linear solver
+ _maxPetscIts=50;
+ _maxNewtonIts=50;
+ _PetscIts=0;//the number of iterations of the linear solver
_ksptype = (char*)&KSPGMRES;
_pctype = (char*)&PCLU;
- _nonLinearSolver = Newton_SOLVERLAB;
+
+ /* Physical parameter */
_heatPowerFieldSet=false;
_heatTransfertCoeff=0;
_rodTemperatureFieldSet=false;
_heatSource=0;
- _verbose = false;
- _system = false;
- _conditionNumber=false;
- _maxvp=0;
- _runLogFile=new ofstream;
//extracting current directory
char result[ PATH_MAX ];
}
if(_FECalculation && VV.getTypeOfField()!=NODES)
{
- *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Element method"<<endl;
+ *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method"<<endl;
+ _runLogFile->close();
+ throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method");
+ }
+ else if(!_FECalculation && VV.getTypeOfField()==NODES)
+ {
+ *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method"<<endl;
_runLogFile->close();
- throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Element method");
+ throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method");
}
_VV=VV;
_VV.setName("SOLVERLAB results");
_time=_VV.getTime();
_mesh=_VV.getMesh();
+ _initialDataSet=true;
+
+ //Mesh data
_Nmailles = _mesh.getNumberOfCells();
_Nnodes = _mesh.getNumberOfNodes();
_Nfaces = _mesh.getNumberOfFaces();
_perimeters=Field("Perimeters", CELLS, _mesh,1);
- // find _minl and maximum nb of neibourghs
+ // find _minl (delta x) and maximum nb of neibourghs
_minl = INFINITY;
int nbNeib,indexFace;
Cell Ci;
Face Fk;
if(_verbose)
- cout<<"Computing cell perimeters and mesh minimal diameter"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Computing cell perimeters and mesh minimal diameter\n");
+ //Compute the maximum number of neighbours for nodes or cells
if(VV.getTypeOfField()==NODES)
- {
- _minl = _mesh.getMaxNbNeighbours(NODES);
_neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
- }
else
- for (int i=0; i<_mesh.getNumberOfCells(); i++){
- Ci = _mesh.getCell(i);
- //Detect mesh with junction
- nbNeib=0;
- for(int j=0; j<Ci.getNumberOfFaces(); j++){
- Fk=_mesh.getFace(Ci.getFacesId()[j]);
- nbNeib+=Fk.getNumberOfCells()-1;
- }
- if(nbNeib>_neibMaxNb)
- _neibMaxNb=nbNeib;
- //Compute mesh data
- 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();
- }
- }
- _initialDataSet=true;
-
+ _neibMaxNb=_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();
+ }
+ }else{
+ _minl = min(_minl,Ci.getMeasure());
+ _perimeters(i)=Ci.getNumberOfFaces();
+ }
+ }
if(_verbose)
cout<<_perimeters<<endl;
}
-void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber)
+//Function needed because swig of enum EntityType fails
+void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int field_support_type)
{
- Field VV(fileName,CELLS,fieldName,timeStepNumber,0);
+ if(_FECalculation && field_support_type!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
+
+ Field VV;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ VV = Field(fileName, CELLS, fieldName, timeStepNumber, 0);
+ break;
+ case NODES:
+ VV = Field(fileName, NODES, fieldName, timeStepNumber, 0);
+ break;
+ case FACES:
+ VV = Field(fileName, FACES, fieldName, timeStepNumber, 0);
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+
+ setInitialField(VV);
+}
+//Function needed because swig of enum EntityType fails
+void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
+ double ymin, double ymax, int ny, string backSide, string frontSide,
+ double zmin, double zmax, int nz, string bottomSide, string topSide, int field_support_type)
+{
+ if(_FECalculation && field_support_type!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
+
+ EntityType typeField;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ typeField=CELLS;
+ break;
+ case NODES:
+ typeField=NODES;
+ break;
+ case FACES:
+ typeField=FACES;
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+
+ setInitialFieldConstant( nDim, Vconstant, xmin, xmax, nx, leftSide, rightSide, ymin, ymax, ny, backSide, frontSide, zmin, zmax, nz, bottomSide, topSide, typeField);
+}
+void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, EntityType typeField)
+{
+ if(_FECalculation && typeField!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
+
+ Field VV(fileName, typeField, fieldName, timeStepNumber, 0);
+
setInitialField(VV);
}
-void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant)
+void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant, EntityType typeField)
{
+ if(_FECalculation && typeField!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
Mesh M(fileName);
- Field VV("SOLVERLAB results", CELLS, M, Vconstant.size());
+ Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
- for (int j = 0; j < M.getNumberOfCells(); j++) {
+ for (int j = 0; j < VV.getNumberOfElements(); j++) {
for (int i=0; i< VV.getNumberOfComponents(); i++)
VV(j,i) = Vconstant[i];
}
setInitialField(VV);
}
-void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant)
+void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant, EntityType typeField)
{
- Field VV("SOLVERLAB results", CELLS, M, Vconstant.getNumberOfRows());
+ if(_FECalculation && typeField!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
+
+ Field VV("SOLVERLAB results", typeField, M, Vconstant.getNumberOfRows());
- for (int j = 0; j < M.getNumberOfCells(); j++) {
+ for (int j = 0; j < VV.getNumberOfElements(); j++) {
for (int i=0; i< VV.getNumberOfComponents(); i++)
VV(j,i) = Vconstant(i);
}
setInitialField(VV);
}
-void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant)
+void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant, EntityType typeField)
{
- Field VV("SOLVERLAB results", CELLS, M, Vconstant.size());
+ if(_FECalculation && typeField!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
- for (int j = 0; j < M.getNumberOfCells(); j++) {
+ Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
+
+ for (int j = 0; j < VV.getNumberOfElements(); j++) {
for (int i=0; i< VV.getNumberOfComponents(); i++)
VV(j,i) = Vconstant[i];
}
}
void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
double ymin, double ymax, int ny, string backSide, string frontSide,
- double zmin, double zmax, int nz, string bottomSide, string topSide)
+ double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
{
Mesh M;
- if(nDim==1){
- //cout<<"coucou1 xmin="<<xmin<<", xmax= "<<xmax<< ", nx= "<<nx<<endl;
+ if(nDim==1)
M=Mesh(xmin,xmax,nx);
- //cout<<"coucou2"<<endl;
- }
else if(nDim==2)
M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
else if(nDim==3)
M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
else{
- cout<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3\n");
*_runLogFile<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
_runLogFile->close();
throw CdmathException("Space dimension nDim should be between 1 and 3");
M.setGroupAtPlan(zmin,2,_precision,bottomSide);
}
- setInitialFieldConstant(M, Vconstant);
+ setInitialFieldConstant(M, Vconstant, typeField);
}
-void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction)
+void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction, EntityType typeField)
{
+ if(_FECalculation && typeField!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
+
if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
{
*_runLogFile<<"ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes"<<endl;
_runLogFile->close();
throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes");
}
- Field VV("SOLVERLAB results", CELLS, M, VV_Left.getNumberOfRows());
+ Field VV("SOLVERLAB results", typeField, M, VV_Left.getNumberOfRows());
double component_value;
- for (int j = 0; j < M.getNumberOfCells(); j++) {
- if(direction==0)
- component_value=M.getCell(j).x();
- else if(direction==1)
- component_value=M.getCell(j).y();
- else if(direction==2)
- component_value=M.getCell(j).z();
- else{
+ for (int j = 0; j < VV.getNumberOfElements(); j++)
+ {
+ if(direction<=2)
+ component_value=VV.getElementComponent(j, direction);
+ else
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"Error : space dimension is %d, direction asked is \%d \n",M.getSpaceDimension(),direction);
_runLogFile->close();
throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: direction should be an integer between 0 and 2");
}
void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
double xmin, double xmax, int nx, string leftSide, string rightSide,
double ymin, double ymax, int ny, string backSide, string frontSide,
- double zmin, double zmax, int nz, string bottomSide, string topSide)
+ double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
{
Mesh M;
if(nDim==1)
V_Left(i)=VV_Left[i];
V_Right(i)=VV_Right[i];
}
- setInitialFieldStepFunction(M, V_Left, V_Right, xstep);
+ setInitialFieldStepFunction(M, V_Left, V_Right, xstep, typeField);
}
-void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center)
+void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center, EntityType typeField)
{
+ if(_FECalculation && typeField!= NODES)
+ cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
+
if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
{
- cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Vout.size() = %d, Vin.size()= %d, Center.size() = %d, M.getSpaceDim = %d \n",Vout.size(),Vin.size(),Center.size(), M.getSpaceDimension());
throw CdmathException("ProblemCoreFlows::setInitialFieldSphericalStepFunction : Vector size error");
}
int nVar=Vout.size();
int spaceDim=M.getSpaceDimension();
- Field VV("Primitive variables for spherical step function", CELLS, M, nVar);
+ Field VV("Primitive variables for spherical step function", typeField, M, nVar);
Vector currentPoint(spaceDim);
- for(int i=0;i<M.getNumberOfCells();i++)
+ for(int i=0;i<VV.getNumberOfElements();i++)
{
- currentPoint(0)=M.getCell(i).x();
+ currentPoint(0)=VV.getElementComponent(i,0);
if(spaceDim>1)
{
- currentPoint(1)=M.getCell(i).y();
+ currentPoint(1)=VV.getElementComponent(i,1);
if(spaceDim>2)
- currentPoint(2)=M.getCell(i).z();
+ currentPoint(2)=VV.getElementComponent(i,2);
}
if((currentPoint-Center).norm()<radius)
for(int j=0;j<nVar;j++)
else if (kspType==BCGS)
_ksptype = (char*)&KSPBCGS;
else {
- cout << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!\n");
*_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
_runLogFile->close();
throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
}
// set preconditioner
- if (pcType == NONE)
+ if (pcType == NOPC)
_pctype = (char*)&PCNONE;
else if (pcType ==LU)
_pctype = (char*)&PCLU;
else if (pcType == ICC)
_pctype = (char*)&PCICC;
else {
- cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
- *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!\n");
+ *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
_runLogFile->close();
- throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
+ throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
}
}
bool ok; // Is the time interval successfully solved ?
_isStationary=false;//in case of a second run with a different physics or cfl
- cout<< "Running test case "<< _fileName<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Running test case %d\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;
// Guess the next time step length
_dt=computeTimeStep(stop);
if (stop){
- cout << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Failed computing time step %d, time = %.2f, dt= %.2f, stopping calculation",_nbTimeStep,_time,_dt);
*_runLogFile << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
break;
}
stop=!initTimeStep(_dt);
// Prepare the next time step
if (stop){
- cout << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Failed initializing time step %d, time = %.2f, dt= %.2f, stopping calculation",_nbTimeStep,_time,_dt);
*_runLogFile << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
break;
}
//_dt=computeTimeStep(stop);
}
else{*/
- cout << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Failed solving time step %d, time = %.2f, dt= %.2f, cfl = %.2f, stopping calculation \n",_nbTimeStep,_time,_dt,_cfl);
*_runLogFile << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
stop=true; // Impossible to solve the next time step, the Problem has given up
break;
{
validateTimeStep();
if ((_nbTimeStep-1)%_freqSave ==0){
- cout << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Time step = %d, dt = %.2f, time = %.2f, ||Un+1-Un||= %.2f\n\n",_nbTimeStep,_dt,_time,_erreur_rel);
*_runLogFile << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl<<endl;
}
}
}
}
if(_isStationary){
- cout << "Stationary state reached" <<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Stationary state reached\n");
*_runLogFile << "Stationary state reached" <<endl;
}
else if(_time>=_timeMax){
- cout<<"Maximum time "<<_timeMax<<" reached"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Maximum time %.2f reached\n",_timeMax);
*_runLogFile<<"Maximum time "<<_timeMax<<" reached"<<endl;
}
else if(_nbTimeStep>=_maxNbOfTimeStep){
- cout<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Maximum number of time steps %d reached\n",_maxNbOfTimeStep);
*_runLogFile<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
}
else{
- cout<<"Error problem wants to stop!"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Error problem wants to stop!\n");
*_runLogFile<<"Error problem wants to stop!"<<endl;
}
- cout << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"End of calculation at time t = %.2f and time step number %d\n",_time,_nbTimeStep);
*_runLogFile << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
_runLogFile->close();
if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)//To monitor the convergence of the newton scheme
{
- cout << " Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
+ PetscPrintf(PETSC_COMM_WORLD," Newton iteration %d, %s iterations : %d maximum variation ||Uk+1-Uk||: %.2f\n",_NEWTON_its,_ksptype,_PetscIts,_erreur_rel);
*_runLogFile<< " Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
if(_conditionNumber)
{
PetscReal sv_max, sv_min;
KSPComputeExtremeSingularValues(_ksp, &sv_max, &sv_min);
- cout<<" Singular value max = " << sv_max <<", singular value min = " << sv_min <<", condition number = " << sv_max/sv_min <<endl;
+ PetscPrintf(PETSC_COMM_WORLD," Singular value max = %.2f, singular value min = %.2f, condition number = %.2f\n",sv_max,sv_min,sv_max/sv_min);
*_runLogFile<<" Singular value max = " << sv_max <<", singular value min = " << sv_min <<", condition number = " << sv_max/sv_min <<endl;
}
}
}
if(!converged){
if(_NEWTON_its >= _maxNewtonIts){
- cout << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Maximum number of Newton iterations %d reached\n",_maxNewtonIts);
*_runLogFile << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
}
else if(!ok){
- cout<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
+ PetscPrintf(PETSC_COMM_WORLD,"iterateTimeStep: solving Newton iteration %d Failed\n",_NEWTON_its);
*_runLogFile<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
}
}
else if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)
{
- cout << "Nombre d'iterations de Newton "<< _NEWTON_its << ", Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl << endl;
+ PetscPrintf(PETSC_COMM_WORLD,"Nombre d'iterations de Newton %d, Nombre max d'iterations %s : %d\n\n",_NEWTON_its, _ksptype, _MaxIterLinearSolver);
*_runLogFile <<endl;
*_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its << "Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl << endl;
_MaxIterLinearSolver = 0;
using namespace std;
-ProblemFluid::ProblemFluid(void){
+ProblemFluid::ProblemFluid(MPI_Comm comm):ProblemCoreFlows(comm)
+{
_latentHeat=1e30;
_Tsat=1e30;
_Psat=-1e30;
#include "StationaryDiffusionEquation.hxx"
+#include "DiffusionEquation.hxx"
#include "Node.hxx"
#include "SparseMatrixPetsc.hxx"
#include "math.h"
using namespace std;
-int StationaryDiffusionEquation::fact(int n)
-{
- return (n == 1 || n == 0) ? 1 : fact(n - 1) * n;
-}
-int StationaryDiffusionEquation::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes)
-{//assumes Dirichlet node numbering is strictly increasing
- int j=0;//indice de parcours des noeuds frontière avec CL Dirichlet
- int boundarySize=dirichletNodes.size();
- while(j<boundarySize and dirichletNodes[j]<globalIndex)
- j++;
- if(j==boundarySize)
- return globalIndex-boundarySize;
- else if (dirichletNodes[j]>globalIndex)
- return globalIndex-j;
- else
- throw CdmathException("StationaryDiffusionEquation::unknownNodeIndex : Error : node is a Dirichlet boundary node");
-}
-
-int StationaryDiffusionEquation::globalNodeIndex(int unknownNodeIndex, std::vector< int > dirichletNodes)
-{//assumes Dirichlet boundary node numbering is strictly increasing
- int boundarySize=dirichletNodes.size();
- /* trivial case where all boundary nodes are Neumann BC */
- if(boundarySize==0)
- return unknownNodeIndex;
-
- double unknownNodeMax=-1;//max unknown node number in the interval between jth and (j+1)th Dirichlet boundary nodes
- int j=0;//indice de parcours des noeuds frontière
- //On cherche l'intervale [j,j+1] qui contient le noeud de numéro interieur unknownNodeIndex
- while(j+1<boundarySize and unknownNodeMax<unknownNodeIndex)
- {
- unknownNodeMax += dirichletNodes[j+1]-dirichletNodes[j]-1;
- j++;
- }
-
- if(j+1==boundarySize)
- return unknownNodeIndex+boundarySize;
- else //unknownNodeMax>=unknownNodeIndex, hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j]
- return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1;
-}
-
-StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda){
+StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalculation, double lambda,MPI_Comm comm){
+ /* Initialisation of PETSC */
+ //check if PETSC is already initialised
PetscBool petscInitialized;
PetscInitialized(&petscInitialized);
if(!petscInitialized)
- PetscInitialize(NULL,NULL,0,0);
+ {//check if MPI is already initialised
+ int mpiInitialized;
+ MPI_Initialized(&mpiInitialized);
+ if(mpiInitialized)
+ 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.
+ PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
if(lambda < 0.)
{
_heatPowerFieldSet=false;
_heatTransfertCoeff=0;
_heatSource=0;
+
+ /* Default diffusion tensor is diagonal */
+ _DiffusionTensor=Matrix(_Ndim);
+ for(int idim=0;idim<_Ndim;idim++)
+ _DiffusionTensor(idim,idim)=_conductivity;
}
void StationaryDiffusionEquation::initialize()
}
}
- _DiffusionTensor=Matrix(_Ndim);
- for(int idim=0;idim<_Ndim;idim++)
- _DiffusionTensor(idim,idim)=1;
/**************** Field creation *********************/
if(!_heatPowerFieldSet){
return _dt_src;
}
-Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values)
-{
- if(! M.isSquare() )
- throw CdmathException("DiffusionEquation::gradientNodal Matrix M should be square !!!");
-
- int Ndim = M.getNumberOfRows();
- vector< Matrix > matrices(Ndim);
-
- for (int idim=0; idim<Ndim;idim++){
- matrices[idim]=M.deepCopy();
- for (int jdim=0; jdim<Ndim+1;jdim++)
- matrices[idim](jdim,idim) = values[jdim] ;
- }
-
- Vector result(Ndim);
- for (int idim=0; idim<Ndim;idim++)
- result[idim] = matrices[idim].determinant();
-
- return result;
-}
-
double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
{
double result;
M(idim,_Ndim)=1;
}
for (int idim=0; idim<_Ndim+1;idim++)
- GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim);
+ 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=unknownNodeIndex(nodeIds[idim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
+ 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= unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing
- MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES);
+ 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
else
valuesBorder[kdim]=0;
}
- GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim);
- coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
+ GradShapeFuncBorder=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim);
+ coeff =-1.*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure();
VecSetValue(_b,i_int,coeff, ADD_VALUES);
}
}
int NboundaryFaces=boundaryFaces.size();
for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
{
- Face Fi = _mesh.getFace(i);
+ 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 an Neumann BC node (not a Dirichlet BC node)
+ 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
+ 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
+ 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);
}
//Compute velocity at the face Fj
- dn=_conductivity*(_DiffusionTensor*normale)*normale;
+ dn=(_DiffusionTensor*normale)*normale;
// compute 1/dxi = volume of Ci/area of Fj
inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
}
- else {//Problem for 3D tetrahera
- //cout<< "Fj.getGroupNames().size()= "<<Fj.getGroupNames().size()<<", Fj.x()= "<<Fj.x()<<", Fj.y()= "<<Fj.y()<<", Fj.z()= "<<Fj.z()<<endl;
+ else {
stop=true ;
cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
cout<<"!!!!!! No boundary condition set for boundary named "<<nameOfGroup<< "!!!!!!!!!! _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
return INFINITY;
}
-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)
+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);
+ VecSetValue(_b,i, _heatTransfertCoeff*_fluidTemperatureField(i) + _heatPowerField(i) ,ADD_VALUES);
else
{
Cell Ci;
if(!_mesh.isBorderNode(nodesId[j]))
{
double coeff = _heatTransfertCoeff*_fluidTemperatureField(nodesId[j]) + _heatPowerField(nodesId[j]);
- VecSetValue(_b,unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
+ VecSetValue(_b,DiffusionEquation::unknownNodeIndex(nodesId[j], _dirichletNodeIds), coeff*Ci.getMeasure()/(_Ndim+1),ADD_VALUES);
}
}
}
VecAssemblyBegin(_deltaT);
VecAssemblyEnd( _deltaT);
- for(int i=0; i<_VV.getNumberOfElements(); i++)
+ if(_verbose)
+ cout<<"Début calcul de la variation relative"<<endl;
+
+ 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;
stop=false;
converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
}
throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
}
// set preconditioner
- if (pcType == NONE)
+ if (pcType == NOPC)
_pctype = (char*)&PCNONE;
else if (pcType ==LU)
_pctype = (char*)&PCLU;
else if (pcType == ICC)
_pctype = (char*)&PCICC;
else {
- cout << "!!! Error : only 'NONE', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
- *_runLogFile << "!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
+ cout << "!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl;
+ *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
_runLogFile->close();
- throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
+ throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
}
}
for(int i=0; i<_NunknownNodes; i++)
{
VecGetValues(_Tk, 1, &i, &Ti);
- globalIndex = globalNodeIndex(i, _dirichletNodeIds);
+ globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds);
_VV(globalIndex)=Ti;
}
{
if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Ci.getNodeId(j))==_dirichletNodeIds.end())//node j is an unknown node (not a Dirichlet node)
{
- j_int=unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
+ j_int=DiffusionEquation::unknownNodeIndex(Ci.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
nodal_volumes[j_int]+=Ci.getMeasure()/(_Ndim+1);
}
}
using namespace std;
-TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport){
+TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport, MPI_Comm comm):ProblemCoreFlows(comm)
+{
if(pEstimate==around1bar300KTransport){
_Tref=300;
if(fluid==GasPhase){//Nitrogen pressure 1 bar and temperature 27°C
