From: michael Date: Thu, 13 Jan 2022 17:27:23 +0000 (+0100) Subject: Some code factoring X-Git-Url: http://git.salome-platform.org/gitweb/?p=tools%2Fsolverlab.git;a=commitdiff_plain;h=2fade9da045cc8c30d1d3fafba05bf970cb5a03a Some code factoring --- diff --git a/CoreFlows/Models/inc/ProblemCoreFlows.hxx b/CoreFlows/Models/inc/ProblemCoreFlows.hxx index 5c038cc..ee2245e 100755 --- a/CoreFlows/Models/inc/ProblemCoreFlows.hxx +++ b/CoreFlows/Models/inc/ProblemCoreFlows.hxx @@ -86,6 +86,14 @@ enum TimeScheme Implicit/**< Implicit numerical scheme */ }; +//! 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) */ +}; + class ProblemCoreFlows { public : @@ -789,6 +797,10 @@ protected : Vec _b;//Linear system right hand side double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps bool _conditionNumber;//computes an estimate of the condition number + /** \fn createKSP + * \brief Create PETSc solver and preconditioner structures + * */ + void createKSP(); //simulation monitoring variables bool _isStationary; diff --git a/CoreFlows/Models/inc/ProblemFluid.hxx b/CoreFlows/Models/inc/ProblemFluid.hxx index b0cd625..7ed239b 100755 --- a/CoreFlows/Models/inc/ProblemFluid.hxx +++ b/CoreFlows/Models/inc/ProblemFluid.hxx @@ -30,22 +30,6 @@ enum SpaceScheme 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 @@ -56,6 +40,14 @@ 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}; diff --git a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx index f0c1568..1365b23 100755 --- a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx +++ b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx @@ -87,6 +87,20 @@ public : 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 @@ -98,6 +112,20 @@ public : 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); diff --git a/CoreFlows/Models/src/ProblemFluid.cxx b/CoreFlows/Models/src/ProblemFluid.cxx index bcd96ba..952a595 100755 --- a/CoreFlows/Models/src/ProblemFluid.cxx +++ b/CoreFlows/Models/src/ProblemFluid.cxx @@ -1,5 +1,6 @@ #include "ProblemFluid.hxx" #include "math.h" +#include using namespace std; @@ -173,8 +174,7 @@ void ProblemFluid::initialize() int *indices = new int[_Nmailles]; - for(int i=0; i<_Nmailles; i++) - indices[i] = i; + std::iota(indices, indices +_Nmailles, 0); VecSetValuesBlocked(_conservativeVars, _Nmailles, indices, initialFieldCons, INSERT_VALUES); VecAssemblyBegin(_conservativeVars); VecAssemblyEnd(_conservativeVars); @@ -197,13 +197,7 @@ void ProblemFluid::initialize() delete[] initialFieldCons; delete[] indices; - //Linear solver - KSPCreate(PETSC_COMM_SELF, &_ksp); - KSPSetType(_ksp, _ksptype); - // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000); - KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts); - KSPGetPC(_ksp, &_pc); - PCSetType(_pc, _pctype); + createKSP(); // Creation du solveur de Newton de PETSc if( _timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB) diff --git a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx index 727d75d..f0d5eb1 100755 --- a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx +++ b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx @@ -71,10 +71,7 @@ StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalcula _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; @@ -224,13 +221,39 @@ void StationaryDiffusionEquation::initialize() _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; @@ -290,7 +313,7 @@ double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop) 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 @@ -405,6 +428,8 @@ double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){ 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) { @@ -528,6 +553,8 @@ double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){ 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) { @@ -743,7 +770,7 @@ void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, precondi 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; @@ -1091,3 +1118,61 @@ StationaryDiffusionEquation::setHeatPowerField(string fileName, string fieldName _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"<