From 0d4b425ed75e0cb6e20b3624a31db7569b47e0bc Mon Sep 17 00:00:00 2001 From: michael Date: Fri, 5 Nov 2021 23:16:32 +0100 Subject: [PATCH] Prepare structures for parallel simulation --- CoreFlows/Models/inc/DiffusionEquation.hxx | 22 +- .../Models/inc/LinearElasticityModel.hxx | 14 +- CoreFlows/Models/inc/ProblemCoreFlows.hxx | 70 +++- CoreFlows/Models/inc/ProblemFluid.hxx | 2 +- .../inc/StationaryDiffusionEquation.hxx | 27 +- CoreFlows/Models/inc/TransportEquation.hxx | 2 +- CoreFlows/Models/src/DiffusionEquation.cxx | 102 +++--- .../Models/src/LinearElasticityModel.cxx | 89 ++--- CoreFlows/Models/src/ProblemCoreFlows.cxx | 318 ++++++++++++------ CoreFlows/Models/src/ProblemFluid.cxx | 3 +- .../src/StationaryDiffusionEquation.cxx | 138 +++----- CoreFlows/Models/src/TransportEquation.cxx | 3 +- 12 files changed, 436 insertions(+), 354 deletions(-) diff --git a/CoreFlows/Models/inc/DiffusionEquation.hxx b/CoreFlows/Models/inc/DiffusionEquation.hxx index 7b74369..2a3110b 100755 --- a/CoreFlows/Models/inc/DiffusionEquation.hxx +++ b/CoreFlows/Models/inc/DiffusionEquation.hxx @@ -5,7 +5,7 @@ * \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. * */ //============================================================================ @@ -13,7 +13,7 @@ /*! \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_ @@ -52,7 +52,7 @@ public : * \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(); @@ -121,9 +121,16 @@ public : 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; @@ -137,21 +144,12 @@ protected : 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 _limitField; diff --git a/CoreFlows/Models/inc/LinearElasticityModel.hxx b/CoreFlows/Models/inc/LinearElasticityModel.hxx index bdaee35..da3fcb8 100755 --- a/CoreFlows/Models/inc/LinearElasticityModel.hxx +++ b/CoreFlows/Models/inc/LinearElasticityModel.hxx @@ -55,7 +55,7 @@ public : * \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;} @@ -164,6 +164,7 @@ protected : double computeRHS(bool & stop); double computeStiffnessMatrixFV(bool & stop); + double computeStiffnessMatrixFE(bool & stop); /************ Data for FE calculation *************/ bool _FECalculation; @@ -175,18 +176,15 @@ protected : 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_ */ diff --git a/CoreFlows/Models/inc/ProblemCoreFlows.hxx b/CoreFlows/Models/inc/ProblemCoreFlows.hxx index 64b2d73..ac9e73c 100755 --- a/CoreFlows/Models/inc/ProblemCoreFlows.hxx +++ b/CoreFlows/Models/inc/ProblemCoreFlows.hxx @@ -90,7 +90,7 @@ class ProblemCoreFlows { public : //! Constructeur par défaut - ProblemCoreFlows(); + ProblemCoreFlows(MPI_Comm comm = MPI_COMM_WORLD); virtual ~ProblemCoreFlows(); // -*-*-*- Gestion du calcul (interface ICoCo) -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* @@ -249,42 +249,57 @@ public : * */ 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 : the value in each cell + * \param [in] EntityType : CELLS, NODES or FACES * \param [out] void * */ - void setInitialFieldConstant(string fileName, const vector Vconstant); + void setInitialFieldConstant(string fileName, const vector 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 + * \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 @@ -303,11 +318,36 @@ public : * \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 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 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 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) @@ -317,9 +357,10 @@ public : * \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 @@ -340,12 +381,13 @@ public : * \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 VV_Left, vector 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 @@ -355,9 +397,10 @@ public : * \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) @@ -606,8 +649,8 @@ public : _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 @@ -716,6 +759,11 @@ protected : 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_ */ diff --git a/CoreFlows/Models/inc/ProblemFluid.hxx b/CoreFlows/Models/inc/ProblemFluid.hxx index f0100f4..9989b61 100755 --- a/CoreFlows/Models/inc/ProblemFluid.hxx +++ b/CoreFlows/Models/inc/ProblemFluid.hxx @@ -84,7 +84,7 @@ public : /**\fn * \brief constructeur de la classe ProblemFluid */ - ProblemFluid(void); + ProblemFluid(MPI_Comm comm = MPI_COMM_WORLD); //Gestion du calcul (interface ICoCo) diff --git a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx index 7bd065e..ab0f1c3 100755 --- a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx +++ b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx @@ -50,7 +50,7 @@ public : * \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){ @@ -146,6 +146,19 @@ public : _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; @@ -205,6 +218,7 @@ protected : double computeRHS(bool & stop); double computeDiffusionMatrixFV(bool & stop); + double computeDiffusionMatrixFE(bool & stop); /************ Data for FE calculation *************/ bool _FECalculation; @@ -216,18 +230,15 @@ protected : 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_ */ diff --git a/CoreFlows/Models/inc/TransportEquation.hxx b/CoreFlows/Models/inc/TransportEquation.hxx index 5acc4ed..6f17e7f 100755 --- a/CoreFlows/Models/inc/TransportEquation.hxx +++ b/CoreFlows/Models/inc/TransportEquation.hxx @@ -64,7 +64,7 @@ public : * \param [in] pressureMagnitude : \ref around1bar or \ref around155bars * \param [in] vector : fluid velocity (assumed constant) * */ - TransportEquation(phase fluid, pressureMagnitude pEstimate,vector vitesseTransport); + TransportEquation(phase fluid, pressureMagnitude pEstimate,vector vitesseTransport, MPI_Comm comm = MPI_COMM_WORLD); //Gestion du calcul virtual void initialize(); diff --git a/CoreFlows/Models/src/DiffusionEquation.cxx b/CoreFlows/Models/src/DiffusionEquation.cxx index 1ac9cc7..b8993c6 100755 --- a/CoreFlows/Models/src/DiffusionEquation.cxx +++ b/CoreFlows/Models/src/DiffusionEquation.cxx @@ -51,7 +51,7 @@ Vector DiffusionEquation::gradientNodal(Matrix M, vector< double > values) 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 values) 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 = "<(0); _dirichletNodeIds=std::vector< int >(0); _NboundaryNodes=0; @@ -118,14 +118,14 @@ DiffusionEquation::DiffusionEquation(int dim, bool FECalculation,double rho,doub _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() @@ -134,7 +134,7 @@ 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()<close(); @@ -145,16 +145,16 @@ void DiffusionEquation::initialize() 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"<close(); throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types"); @@ -196,14 +196,14 @@ void DiffusionEquation::initialize() else if(_Ndim==3) { if( _mesh.isTetrahedral() )//Mesh dim=3 - cout<<"3D Finite element method on tetrahedra"<close(); throw CdmathException("DiffusionEquation::initialize: mesh has incorrect cell types"); @@ -214,14 +214,14 @@ void DiffusionEquation::initialize() _NboundaryNodes=_boundaryNodeIds.size(); if(_NboundaryNodes==_Nnodes) - cout<<"!!!!! Warning : all nodes are boundary nodes !!!!!"< division by zero"); else @@ -436,7 +438,7 @@ double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){ } //Compute velocity at the face Fj - dn=_conductivity*(_DiffusionTensor*normale)*normale; + dn=(_DiffusionTensor*normale)*normale; if(fabs(dn)>_maxvp) _maxvp=fabs(dn); @@ -514,22 +516,19 @@ double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){ 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; @@ -541,7 +540,7 @@ double DiffusionEquation::computeRHS(bool & stop){ for (int j=0; jdirichletNodes.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); } } @@ -549,8 +548,10 @@ double DiffusionEquation::computeRHS(bool & stop){ 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; @@ -576,7 +577,7 @@ bool DiffusionEquation::initTimeStep(double dt){ } else//dt<=0 { - cout<<"DiffusionEquation::initTimeStep dt= "<0; } @@ -642,7 +646,7 @@ bool DiffusionEquation::iterateTimeStep(bool &converged) _MaxIterLinearSolver = _PetscIts; if(_PetscIts>=_maxPetscIts) { - cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"< 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(jglobalIndex) - 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=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.) { @@ -263,22 +236,6 @@ void LinearElasticityModel::initialize() } -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; @@ -328,20 +285,20 @@ double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){ 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) @@ -360,7 +317,7 @@ double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){ 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); } @@ -376,12 +333,12 @@ double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){ 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 @@ -540,7 +497,7 @@ double LinearElasticityModel::computeRHS(bool & stop)//Contribution of the PDE R for (int j=0; j1) + { + 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; @@ -34,35 +57,45 @@ ProblemCoreFlows::ProblemCoreFlows() _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 ]; @@ -126,96 +159,163 @@ void ProblemCoreFlows::setInitialField(const Field &VV) } if(_FECalculation && VV.getTypeOfField()!=NODES) { - *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Element method"<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"<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"<_neibMaxNb) - _neibMaxNb=nbNeib; - //Compute mesh data - if (_Ndim > 1){ - _perimeters(i)=0; - for (int k=0 ; k 1){ + _perimeters(i)=0; + for (int k=0 ; k 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!!!"< Vconstant) +void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector Vconstant, EntityType typeField) { + if(_FECalculation && typeField!= NODES) + cout<<"Warning : finite element simulation should have initial field on nodes!!!"< Vconstant) +void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector 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!!!"< 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="<close(); throw CdmathException("Space dimension nDim should be between 1 and 3"); @@ -253,28 +350,30 @@ void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector V 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!!!"<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"); } @@ -290,7 +389,7 @@ void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector VV_Left, vector 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) @@ -320,29 +419,32 @@ void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector1) { - 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()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; @@ -403,10 +505,10 @@ void ProblemCoreFlows::setLinearSolver(linearSolver kspType, preconditioner pcTy 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 !!!" ); } } @@ -435,7 +537,7 @@ bool ProblemCoreFlows::run() 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<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<=_timeMax){ - cout<<"Maximum time "<<_timeMax<<" reached"<=_maxNbOfTimeStep){ - cout<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<close(); @@ -558,14 +660,14 @@ bool ProblemCoreFlows::solveTimeStep(){ 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 < 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(jglobalIndex) - 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=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.) { @@ -121,6 +95,11 @@ StationaryDiffusionEquation::StationaryDiffusionEquation(int dim, bool FECalcula _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() @@ -145,9 +124,6 @@ void StationaryDiffusionEquation::initialize() } } - _DiffusionTensor=Matrix(_Ndim); - for(int idim=0;idim<_Ndim;idim++) - _DiffusionTensor(idim,idim)=1; /**************** Field creation *********************/ if(!_heatPowerFieldSet){ @@ -279,27 +255,6 @@ double StationaryDiffusionEquation::computeTimeStep(bool & stop){ 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; idimclose(); - throw CdmathException("!!! Error : only 'NONE' or 'LU' or 'ILU' preconditioners are acceptable !!!" ); + throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" ); } } @@ -881,7 +841,7 @@ void StationaryDiffusionEquation::save(){ for(int i=0; i<_NunknownNodes; i++) { VecGetValues(_Tk, 1, &i, &Ti); - globalIndex = globalNodeIndex(i, _dirichletNodeIds); + globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds); _VV(globalIndex)=Ti; } @@ -962,7 +922,7 @@ StationaryDiffusionEquation::getEigenvalues(int nev, EPSWhich which, double tol) { 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); } } diff --git a/CoreFlows/Models/src/TransportEquation.cxx b/CoreFlows/Models/src/TransportEquation.cxx index cdd5413..de715fd 100755 --- a/CoreFlows/Models/src/TransportEquation.cxx +++ b/CoreFlows/Models/src/TransportEquation.cxx @@ -5,7 +5,8 @@ using namespace std; -TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,vector vitesseTransport){ +TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,vector vitesseTransport, MPI_Comm comm):ProblemCoreFlows(comm) +{ if(pEstimate==around1bar300KTransport){ _Tref=300; if(fluid==GasPhase){//Nitrogen pressure 1 bar and temperature 27°C -- 2.39.2