]> SALOME platform Git repositories - tools/solverlab.git/commitdiff
Salome HOME
Some code factoring
authormichael <michael@localhost.localdomain>
Thu, 13 Jan 2022 17:27:23 +0000 (18:27 +0100)
committermichael <michael@localhost.localdomain>
Thu, 13 Jan 2022 17:27:23 +0000 (18:27 +0100)
CoreFlows/Models/inc/ProblemCoreFlows.hxx
CoreFlows/Models/inc/ProblemFluid.hxx
CoreFlows/Models/inc/StationaryDiffusionEquation.hxx
CoreFlows/Models/src/ProblemFluid.cxx
CoreFlows/Models/src/StationaryDiffusionEquation.cxx

index 5c038cc8ef9ddc01ab568a4983c854aadda5003c..ee2245e6308230a4b484444e59f5c5a6193ac635 100755 (executable)
@@ -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;
index b0cd6255a6d82763eeb74734c12f2def5625147a..7ed239bd262b0bd771d3cf056030841d293dfd5e 100755 (executable)
@@ -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};
index f0c1568618a65ca6179cb942d091071d6685ce3d..1365b2357bf89fdc3f2a4afcc6fe3918490d1845 100755 (executable)
@@ -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);
index bcd96ba7da430d3eb793367d47228a6460ae9b1c..952a5958ea7cce38378c8b95bdbf395d0b1da38c 100755 (executable)
@@ -1,5 +1,6 @@
 #include "ProblemFluid.hxx"
 #include "math.h"
+#include  <numeric>
 
 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)
index 727d75d38f892f373ebb2d8594f6aa497f90af5c..f0d5eb1a91ddf94713599a0abe6b5aaf2ab32a5a 100755 (executable)
@@ -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"<<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 
+}