From 381dfb91a247806b9b2959852dc517d8f9bba328 Mon Sep 17 00:00:00 2001 From: michael Date: Sat, 30 Oct 2021 22:42:41 +0200 Subject: [PATCH] Improved treatment of boundary conditions --- CoreFlows/Models/inc/DiffusionEquation.hxx | 19 ++- CoreFlows/Models/src/DiffusionEquation.cxx | 179 +++++++++++---------- 2 files changed, 112 insertions(+), 86 deletions(-) diff --git a/CoreFlows/Models/inc/DiffusionEquation.hxx b/CoreFlows/Models/inc/DiffusionEquation.hxx index cb922c6..7b74369 100755 --- a/CoreFlows/Models/inc/DiffusionEquation.hxx +++ b/CoreFlows/Models/inc/DiffusionEquation.hxx @@ -5,12 +5,15 @@ * \version 1.0 * \date 24 March 2015 * \brief Heat diffusion equation + * 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) */ #ifndef DiffusionEquation_HXX_ #define DiffusionEquation_HXX_ @@ -72,7 +75,7 @@ public : * \param [in] double : the value of the temperature at the boundary * \param [out] void * */ - void setDirichletBoundaryCondition(string groupName,double Temperature){ + void setDirichletBoundaryCondition(string groupName,double Temperature=0){ _limitField[groupName]=LimitFieldDiffusion(DirichletDiffusion,Temperature,-1); }; /** \fn setNeumannBoundaryCondition @@ -85,6 +88,9 @@ public : _limitField[groupName]=LimitFieldDiffusion(NeumannDiffusion,-1, normalFlux); }; + void setDirichletValues(map< int, double> dirichletBoundaryValues); + void setNeumannValues(map< int, double> neumannBoundaryValues); + void setRodDensity(double rho){ _rho=rho; }; @@ -132,6 +138,7 @@ protected : /************ 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 */ @@ -139,7 +146,7 @@ protected : 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 + 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); @@ -147,6 +154,12 @@ protected : TimeScheme _timeScheme; map _limitField; -}; + + /********* Possibility to set a boundary field as DirichletNeumann boundary condition *********/ + bool _dirichletValuesSet; + bool _neumannValuesSet; + std::map< int, double> _dirichletBoundaryValues; + std::map< int, double> _neumannBoundaryValues; + }; #endif /* DiffusionEquation_HXX_ */ diff --git a/CoreFlows/Models/src/DiffusionEquation.cxx b/CoreFlows/Models/src/DiffusionEquation.cxx index c93fef8..1ac9cc7 100755 --- a/CoreFlows/Models/src/DiffusionEquation.cxx +++ b/CoreFlows/Models/src/DiffusionEquation.cxx @@ -42,21 +42,26 @@ int DiffusionEquation::globalNodeIndex(int unknownNodeIndex, std::vector< int > if(j+1==boundarySize) return unknownNodeIndex+boundarySize; - else //unknownNodeMax>=unknownNodeIndex) hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j] + else //unknownNodeMax>=unknownNodeIndex, hence our node global number is between dirichletNodes[j-1] and dirichletNodes[j] return unknownNodeIndex - unknownNodeMax + dirichletNodes[j]-1; } -Vector DiffusionEquation::gradientNodal(Matrix M, vector< double > values){ - vector< Matrix > matrices(_Ndim); +Vector DiffusionEquation::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++){ + for (int idim=0; idim_maxvp) - _maxvp=fabs(dij); + MatSetValue(_A,i_int,j_int,_conductivity*(_DiffusionTensor*GradShapeFuncs[idim])*GradShapeFuncs[jdim]/Cj.getMeasure(), ADD_VALUES); } else if (!dirichletCell_treated) {//Second node of the edge is a Dirichlet node dirichletCell_treated=true; for (int kdim=0; kdim<_Ndim+1;kdim++) { - if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),nodeIds[kdim])!=_dirichletNodeIds.end()) - valuesBorder[kdim]=_limitField[nameOfGroup].T; + std::map::iterator it=_dirichletBoundaryValues.find(nodeIds[kdim]); + if( it != _dirichletBoundaryValues.end() ) + { + if( _dirichletValuesSet ) + valuesBorder[kdim]=_dirichletBoundaryValues[it->second]; + else + valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].T; + } else valuesBorder[kdim]=0; } GradShapeFuncBorder=gradientNodal(M,valuesBorder)/fact(_Ndim); - dij =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure(); - VecSetValue(_b,i_int,dij, ADD_VALUES); + coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure(); + VecSetValue(_b,i_int,coeff, ADD_VALUES); } } } } } + //Calcul de la contribution de la condition limite de Neumann au second membre + if( _NdirichletNodes !=_NboundaryNodes) + { + vector< int > boundaryFaces = _mesh.getBoundaryFaceIds(); + int NboundaryFaces=boundaryFaces.size(); + for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord + { + Face Fi = _mesh.getFace(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 + if( _neumannValuesSet ) + coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)]; + 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); @@ -384,11 +399,11 @@ double DiffusionEquation::computeDiffusionMatrixFE(bool & stop){ _diffusionMatrixSet=true; stop=false ; - cout<<"_maxvp= "<<_maxvp<< " cfl= "<<_cfl<<" minl= "<<_minl< division by zero"); else - return _cfl/_maxvp; + return _cfl*_minl*_minl/_maxvp; } double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){ @@ -492,9 +507,9 @@ double DiffusionEquation::computeDiffusionMatrixFV(bool & stop){ _diffusionMatrixSet=true; stop=false ; - cout<<"_maxvp= "<<_maxvp<< " cfl= "<<_cfl<<" minl= "<<_minl< division by zero"); else return _cfl*_minl*_minl/_maxvp; } @@ -632,28 +647,20 @@ bool DiffusionEquation::iterateTimeStep(bool &converged) converged=false; return false; } - else{ + else + { VecCopy(_Tk, _deltaT);//ici on a deltaT=Tk VecAXPY(_deltaT, -1, _Tkm1);//On obtient deltaT=Tk-Tkm1 _erreur_rel= 0; double Ti, dTi; - if(!_FECalculation) - for(int i=0; i<_Nmailles; i++) - { - VecGetValues(_deltaT, 1, &i, &dTi); - VecGetValues(_Tk, 1, &i, &Ti); - if(_erreur_rel < fabs(dTi/Ti)) - _erreur_rel = fabs(dTi/Ti); - } - else - for(int i=0; i<_Nnodes; i++) - { - VecGetValues(_deltaT, 1, &i, &dTi); - VecGetValues(_Tk, 1, &i, &Ti); - if(_erreur_rel < fabs(dTi/Ti)) - _erreur_rel = fabs(dTi/Ti); - } + for(int i=0; i<_VV.getNumberOfElements(); i++) + { + VecGetValues(_deltaT, 1, &i, &dTi); + VecGetValues(_Tk, 1, &i, &Ti); + if(_erreur_rel < fabs(dTi/Ti)) + _erreur_rel = fabs(dTi/Ti); + } converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton } } @@ -670,22 +677,13 @@ void DiffusionEquation::validateTimeStep() _erreur_rel= 0; double Ti, dTi; - if(!_FECalculation) - for(int i=0; i<_Nmailles; i++) - { - VecGetValues(_deltaT, 1, &i, &dTi); - VecGetValues(_Tk, 1, &i, &Ti); - if(_erreur_rel < fabs(dTi/Ti)) - _erreur_rel = fabs(dTi/Ti); - } - else - for(int i=0; i<_Nnodes; i++) - { - VecGetValues(_deltaT, 1, &i, &dTi); - VecGetValues(_Tk, 1, &i, &Ti); - if(_erreur_rel < fabs(dTi/Ti)) - _erreur_rel = fabs(dTi/Ti); - } + for(int i=0; i<_VV.getNumberOfElements(); i++) + { + VecGetValues(_deltaT, 1, &i, &dTi); + VecGetValues(_Tk, 1, &i, &Ti); + if(_erreur_rel < fabs(dTi/Ti)) + _erreur_rel = fabs(dTi/Ti); + } _isStationary =(_erreur_rel <_precision); @@ -806,6 +804,21 @@ void DiffusionEquation::terminate(){ MatDestroy(&_A); } +void +DiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues) +{ + _dirichletValuesSet=true; + _dirichletBoundaryValues=dirichletBoundaryValues; +} + +void +DiffusionEquation::setNeumannValues(map< int, double> neumannBoundaryValues) +{ + _neumannValuesSet=true; + _neumannBoundaryValues=neumannBoundaryValues; +} + + vector DiffusionEquation::getOutputFieldsNames() { vector result(2); -- 2.39.2