staggered,/**< scheme inspired by staggered discretisations */
};
-//! enumeration pressureEstimate
-/*! the pressure estimate needed to fit physical parameters */
-enum pressureEstimate
-{
- around1bar300K,/**< pressure is around 1 bar and temperature around 300K (for TransportEquation, SinglePhase and IsothermalTwoFluid) or 373 K (saturation for DriftModel and FiveEqsTwoFluid) */
- around155bars600K/**< pressure is around 155 bars and temperature around 618 K (saturation) */
-};
-
-//! enumeration phaseType
-/*! The fluid type can be Gas or water */
-enum phaseType
-{
- Liquid,/**< Fluid considered is water */
- Gas/**< Fluid considered is Gas */
-};
-
//! enumeration NonLinearFormulation
/*! the formulation used to compute the non viscous fluxes */
enum NonLinearFormulation
reducedRoe,/**< compacted formulation of Roe scheme without computation of the fluxes */
};
+//! enumeration phaseType
+/*! The material phase can be Gas or liquid */
+enum phaseType
+{
+ Liquid,/**< Material considered is Liquid */
+ Gas/**< Material considered is Gas */
+};
+
//! enumeration BoundaryType
/*! Boundary condition type */
enum BoundaryType {Wall, InnerWall, Inlet, InletPressure, InletRotationVelocity, InletEnthalpy, Outlet, Neumann, NoTypeSpecified};
void setDirichletBoundaryCondition(string groupName,double Temperature){
_limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,Temperature,-1);
};
+ /** \fn setDirichletBoundaryCondition
+ * \brief adds a new boundary condition of type Dirichlet
+ * \details Reads the boundary field in a med file
+ * \param [in] string : the name of the boundary
+ * \param [in] string : the file name
+ * \param [in] string : the field name
+ * \param [in] int : the time step number
+ * \param [in] int : int corresponding to the enum CELLS or NODES
+ * \param [out] void
+ * */
+ void setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
+ void setDirichletBoundaryCondition(string groupName, Field bc_field){
+ _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion, 0, -1);
+ };
/** \fn setNeumannBoundaryCondition
* \brief adds a new boundary condition of type Neumann
void setNeumannBoundaryCondition(string groupName, double normalFlux=0){
_limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, normalFlux);
};
+ /** \fn setNeumannBoundaryCondition
+ * \brief adds a new boundary condition of type Neumann
+ * \details Reads the boundary field in a med file
+ * \param [in] string : the name of the boundary
+ * \param [in] string : the file name
+ * \param [in] string : the field name
+ * \param [in] int : the time step number
+ * \param [in] int : int corresponding to the enum CELLS or NODES
+ * \param [out] void
+ * */
+ void setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
+ void setNeumannBoundaryCondition(string groupName, Field bc_field){
+ _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,-1, 0);
+ };
void setDirichletValues(map< int, double> dirichletBoundaryValues);
void setNeumannValues (map< int, double> neumannBoundaryValues);
_NEWTON_its=0;
int _PetscIts=0;//the number of iterations of the linear solver
_ksptype = (char*)&KSPGMRES;
- if( _mpi_size>1)
- _pctype = (char*)&PCNONE;
- else
- _pctype = (char*)&PCLU;
+ _pctype = (char*)&PCILU;
_conditionNumber=false;
_erreur_rel= 0;
_Tk_seq=_Tk;
VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq);
- //Linear solver
+ //PETSc Linear solver
KSPCreate(PETSC_COMM_WORLD, &_ksp);
KSPSetType(_ksp, _ksptype);
- // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
KSPGetPC(_ksp, &_pc);
- PCSetType(_pc, _pctype);
+ //PETSc preconditioner
+ if(_mpi_size==1 )
+ PCSetType(_pc, _pctype);
+ else
+ {
+ PCSetType(_pc, PCBJACOBI);//Global preconditioner is block jacobi
+ if(_pctype != (char*)&PCILU)//Default pc type is ilu
+ {
+ PetscOptionsSetValue(NULL,"-sub_pc_type ",_pctype);
+ PetscOptionsSetValue(NULL,"-sub_ksp_type ","preonly");
+ //If the above setvalue does not work, try the following
+ /*
+ KSPSetUp(_ksp);//to set the block Jacobi data structures (including creation of an internal KSP context for each block)
+ KSP * subKSP;
+ PC subpc;
+ int nlocal;//nb local blocs (should equal 1)
+ PCBJacobiGetSubKSP(_pc,&nlocal,NULL,&subKSP);
+ if(nlocal==1)
+ {
+ KSPSetType(subKSP[0], KSPPREONLY);//local block solver is same as global
+ KSPGetPC(subKSP[0],&subpc);
+ PCSetType(subpc,_pctype);
+ }
+ else
+ throw CdmathException("PC Block Jacobi, more than one block in this processor!!");
+ */
+ }
+ }
//Checking whether all boundary conditions are Neumann boundary condition
//if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
result=computeDiffusionMatrixFV(stop);
MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd( _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(_onlyNeumannBC) //Check that the matrix is symmetric
{
PetscBool isSymetric;
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
MatIsSymmetric(_A,_precision,&isSymetric);
if(!isSymetric)
{
if(_onlyNeumannBC) //Check that the matrix is symmetric
{
PetscBool isSymetric;
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd( _A, MAT_FINAL_ASSEMBLY);
MatIsSymmetric(_A,_precision,&isSymetric);
if(!isSymetric)
{
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;
_heatPowerField.getMesh().checkFastEquivalWith(_mesh);
_heatPowerFieldSet=true;
}
+
+void
+StationaryDiffusionEquation::setDirichletBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
+ if(_FECalculation && field_support_type != NODES)
+ cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+ else if(!_FECalculation && field_support_type == NODES)
+ cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+
+ Field VV;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case NODES:
+ VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case FACES:
+ VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error StationaryDiffusionEquation::setDirichletBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+ /* For the moment the boundary value is taken constant equal to zero */
+ _limitField[groupName]=LimitFieldStationaryDiffusion(DirichletStationaryDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
+}
+
+void
+StationaryDiffusionEquation::setNeumannBoundaryCondition(string groupName, string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type){
+ if(_FECalculation && field_support_type != NODES)
+ cout<<"Warning : finite element simulation should have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+ else if(!_FECalculation && field_support_type == NODES)
+ cout<<"Warning : finite volume simulation should not have boundary field on nodes!!! Change parameter field_support_type"<<endl;
+
+ Field VV;
+
+ switch(field_support_type)
+ {
+ case CELLS:
+ VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case NODES:
+ VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ case FACES:
+ VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
+ break;
+ default:
+ std::ostringstream message;
+ message << "Error StationaryDiffusionEquation::setNeumannBoundaryCondition \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
+ throw CdmathException(message.str().c_str());
+ }
+ /* For the moment the boundary value is taken constant equal to zero */
+ _limitField[groupName]=LimitFieldStationaryDiffusion(NeumannStationaryDiffusion,0,-1);//This line will be deleted when variable BC are properly treated in solverlab
+}
