X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=CoreFlows%2FModels%2Fsrc%2FStationaryDiffusionEquation.cxx;h=f0d5eb1a91ddf94713599a0abe6b5aaf2ab32a5a;hb=2fade9da045cc8c30d1d3fafba05bf970cb5a03a;hp=eccee5edf9de74306573e40370d238d9eb981e12;hpb=8d75250c8a620dc1bab221405876864c1227b9bb;p=tools%2Fsolverlab.git diff --git a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx index eccee5e..f0d5eb1 100755 --- a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx +++ b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx @@ -1,4 +1,5 @@ #include "StationaryDiffusionEquation.hxx" +#include "DiffusionEquation.hxx" #include "Node.hxx" #include "SparseMatrixPetsc.hxx" #include "math.h" @@ -8,61 +9,34 @@ using namespace std; -int StationaryDiffusionEquation::fact(int n) -{ - return (n == 1 || n == 0) ? 1 : fact(n - 1) * n; -} -int StationaryDiffusionEquation::unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes) const -{//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) const -{//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,&_mpi_rank); + MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size); + PetscPrintf(PETSC_COMM_WORLD,"\n Simulation on %d processors\n",_mpi_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",_mpi_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.) { std::cout<<"Conductivity="<open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation - if(!_meshSet) - throw CdmathException("StationaryDiffusionEquation::initialize() set mesh first"); - else - { - cout<<"!!!! Initialisation of the computation of the temperature diffusion in a solid using "; - *_runLogFile<<"!!!!! Initialisation of the computation of the temperature diffusion in a solid using "; - if(!_FECalculation) - { - cout<< "Finite volumes method"<::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]); - if( it != _dirichletBoundaryValues.end() ) - _dirichletNodeIds.push_back(_boundaryNodeIds[i]); - else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 ) - { - cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl; - *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<close(); - throw CdmathException("Missing boundary group"); - } - else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion) - { - cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl; - cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; - *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<close(); - throw CdmathException("Missing boundary condition"); - } - else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion) - _dirichletNodeIds.push_back(_boundaryNodeIds[i]); - else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion) - { - cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl; - cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; - *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<close(); - throw CdmathException("Wrong boundary condition"); - } - } - _NdirichletNodes=_dirichletNodeIds.size(); - _NunknownNodes=_Nnodes - _NdirichletNodes; - cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation + + if(!_meshSet) + throw CdmathException("!!!!!!!!StationaryDiffusionEquation::initialize() set mesh first"); + else + { + cout<<"\n Initialisation of the computation of the temperature diffusion in a solid using "; + *_runLogFile<<"\n Initialisation of the computation of the temperature diffusion in a solid using "; + if(!_FECalculation) + { + cout<< "Finite volumes method"<::iterator it=_dirichletBoundaryValues.find(_boundaryNodeIds[i]); + if( it != _dirichletBoundaryValues.end() ) + _dirichletNodeIds.push_back(_boundaryNodeIds[i]); + else if( _mesh.getNode(_boundaryNodeIds[i]).getGroupNames().size()==0 ) + { + cout<<"!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<< endl; + *_runLogFile<< "!!! No boundary group set for boundary node" << _boundaryNodeIds[i]<close(); + throw CdmathException("Missing boundary group"); + } + else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCStationaryDiffusion) + { + cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl; + cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; + *_runLogFile<< "!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<close(); + throw CdmathException("Missing boundary condition"); + } + else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletStationaryDiffusion) + _dirichletNodeIds.push_back(_boundaryNodeIds[i]); + else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannStationaryDiffusion) + { + cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl; + cout<<"!!! Accepted boundary conditions are DirichletStationaryDiffusion "<< DirichletStationaryDiffusion <<" and NeumannStationaryDiffusion "<< NeumannStationaryDiffusion << endl; + *_runLogFile<< "!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<close(); + throw CdmathException("Wrong boundary condition"); + } + } + _NdirichletNodes=_dirichletNodeIds.size(); + _NunknownNodes=_Nnodes - _NdirichletNodes; + cout<<"Number of unknown nodes " << _NunknownNodes <<", Number of boundary nodes " << _NboundaryNodes<< ", Number of Dirichlet boundary nodes " << _NdirichletNodes <1 && _mpi_rank == 0) + VecCreateSeq(PETSC_COMM_SELF,_globalNbUnknowns,&_Tk_seq);//For saving results on proc 0 + if(_mpi_size==0) + _Tk_seq=_Tk; + VecScatterCreateToZero(_Tk,&_scat,&_Tk_seq); + + //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; @@ -273,14 +280,6 @@ void StationaryDiffusionEquation::initialize() *_runLogFile<<"### The system matrix being singular, we seek a zero sum solution, and exact (LU and CHOLESKY) and incomplete factorisations (ILU and ICC) may fail."< 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 StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop) { double result; + MatZeroEntries(_A); + VecZeroEntries(_b); + if(_FECalculation) result=computeDiffusionMatrixFE(stop); else result=computeDiffusionMatrixFV(stop); + MatAssemblyBegin(_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(_heatTransfertCoeff>_precision) MatShift(_A,_heatTransfertCoeff);//Contribution from the liquit/solid heat transfer if(_verbose or _system) - MatView(_A,PETSC_VIEWER_STDOUT_SELF); + MatView(_A,PETSC_VIEWER_STDOUT_WORLD); return result; } double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){ - Cell Cj; - string nameOfGroup; - double dn, coeff; - MatZeroEntries(_A); - VecZeroEntries(_b); - - Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix - std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes - std::vector< int > nodeIds(_Ndim+1);//cell node Ids - std::vector< Node > nodes(_Ndim+1);//cell nodes - int i_int, j_int; //index of nodes j and k considered as unknown nodes - bool dirichletCell_treated; - - std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node - for (int idim=0; idim<_Ndim+1;idim++) - values[idim][idim]=1; - - /* parameters for boundary treatment */ - vector< double > valuesBorder(_Ndim+1); - Vector GradShapeFuncBorder(_Ndim+1); - - for (int j=0; j<_Nmailles;j++) - { - Cj = _mesh.getCell(j); - - for (int idim=0; idim<_Ndim+1;idim++){ - nodeIds[idim]=Cj.getNodeId(idim); - nodes[idim]=_mesh.getNode(nodeIds[idim]); - for (int jdim=0; jdim<_Ndim;jdim++) - M(idim,jdim)=nodes[idim].getPoint()[jdim]; - M(idim,_Ndim)=1; - } - for (int idim=0; idim<_Ndim+1;idim++) - GradShapeFuncs[idim]=gradientNodal(M,values[idim])/fact(_Ndim); + if(_mpi_rank == 0) + { + Cell Cj; + string nameOfGroup; + double coeff; + + Matrix M(_Ndim+1,_Ndim+1);//cell geometry matrix + std::vector< Vector > GradShapeFuncs(_Ndim+1);//shape functions of cell nodes + std::vector< int > nodeIds(_Ndim+1);//cell node Ids + std::vector< Node > nodes(_Ndim+1);//cell nodes + int i_int, j_int; //index of nodes j and k considered as unknown nodes + bool dirichletCell_treated; + + std::vector< vector< double > > values(_Ndim+1,vector< double >(_Ndim+1,0));//values of shape functions on cell node + for (int idim=0; idim<_Ndim+1;idim++) + values[idim][idim]=1; + + /* parameters for boundary treatment */ + vector< double > valuesBorder(_Ndim+1); + Vector GradShapeFuncBorder(_Ndim+1); + + for (int j=0; j<_Nmailles;j++) + { + Cj = _mesh.getCell(j); + + for (int idim=0; idim<_Ndim+1;idim++){ + nodeIds[idim]=Cj.getNodeId(idim); + nodes[idim]=_mesh.getNode(nodeIds[idim]); + for (int jdim=0; jdim<_Ndim;jdim++) + M(idim,jdim)=nodes[idim].getPoint()[jdim]; + M(idim,_Ndim)=1; + } + for (int idim=0; idim<_Ndim+1;idim++) + 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=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= DiffusionEquation::unknownNodeIndex(nodeIds[jdim], _dirichletNodeIds);//assumes Dirichlet boundary node numbering is strictly increasing + MatSetValue(_A,i_int,j_int,(_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++) + { + 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=DiffusionEquation::gradientNodal(M,valuesBorder)/DiffusionEquation::fact(_Ndim); + coeff =-1.*(_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(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 a Neumann BC node (not a Dirichlet BC node) + { + j_int=DiffusionEquation::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); + } + } + } + } + } - /* 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 - 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 - 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++) - { - 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); - coeff =-_conductivity*(_DiffusionTensor*GradShapeFuncBorder)*GradShapeFuncs[idim]/Cj.getMeasure(); - VecSetValue(_b,i_int,coeff, ADD_VALUES); - } - } - } - } + 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) + { + cout<<"Singular matrix is not symmetric, tolerance= "<< _precision< 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); - VecAssemblyEnd(_b); _diffusionMatrixSet=true; stop=false ; @@ -448,149 +445,160 @@ double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){ } double StationaryDiffusionEquation::computeDiffusionMatrixFV(bool & stop){ - long nbFaces = _mesh.getNumberOfFaces(); - Face Fj; - Cell Cell1,Cell2; - string nameOfGroup; - double inv_dxi, inv_dxj; - double barycenterDistance; - Vector normale(_Ndim); - double dn; - PetscInt idm, idn; - std::vector< int > idCells; - MatZeroEntries(_A); - VecZeroEntries(_b); - - for (int j=0; j::iterator it=_dirichletBoundaryValues.find(j); - if( it != _dirichletBoundaryValues.end() ) - { - barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); - MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); - VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES); - } - else - { - nameOfGroup = Fj.getGroupName(); - - if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){ - VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES); - } - else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){ - barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); - MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); - VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES); - } - else { - stop=true ; - cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"< idCells; + + for (int j=0; j::iterator it=_dirichletBoundaryValues.find(j); + if( it != _dirichletBoundaryValues.end() ) + { + barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); + MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); + VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*it->second, ADD_VALUES); + } + else + { + nameOfGroup = Fj.getGroupName(); + + if (_limitField[nameOfGroup].bcType==NeumannStationaryDiffusion){ + VecSetValue(_b,idm, -dn*inv_dxi*_limitField[nameOfGroup].normalFlux, ADD_VALUES); + } + else if(_limitField[nameOfGroup].bcType==DirichletStationaryDiffusion){ + barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter()); + MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES); + VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES); + } + else { + stop=true ; + cout<<"!!!!!!!!!!!!!!! Error StationaryDiffusionEquation::computeDiffusionMatrixFV !!!!!!!!!!"< 1) + inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); + else + inv_dxj = 1/Cell2.getMeasure(); + + barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter()); + + MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES); + MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES); } - Cell2 = _mesh.getCell(idCells[1]); - idn = idCells[1]; - if (_Ndim > 1) - inv_dxj = Fj.getMeasure()/Cell2.getMeasure(); else - inv_dxj = 1/Cell2.getMeasure(); - - barycenterDistance=Cell1.getBarryCenter().distance(Cell2.getBarryCenter()); - - MatSetValue(_A,idm,idm, dn*inv_dxi/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idm,idn,-dn*inv_dxi/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idn,idn, dn*inv_dxj/barycenterDistance, ADD_VALUES); - MatSetValue(_A,idn,idm,-dn*inv_dxj/barycenterDistance, ADD_VALUES); + { + *_runLogFile<<"StationaryDiffusionEquation::computeDiffusionMatrixFV(): incompatible number of cells around a face"< nodesId; - for (int i=0; i<_Nmailles;i++) - { - Ci=_mesh.getCell(i); - nodesId=Ci.getNodesId(); - for (int j=0; j nodesId; + for (int i=0; i<_Nmailles;i++) + { + Ci=_mesh.getCell(i); + nodesId=Ci.getNodesId(); + for (int j=0; jclose(); - throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"); + if(_Ndim != M.getSpaceDimension() or _Ndim!=M.getMeshDimension())//for the moment we must have space dim=mesh dim + { + cout<< "Problem : dimension defined is "<<_Ndim<< " but mesh dimension= "<close(); + throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect dimension"); + } + + _mesh=M; + _Nmailles = _mesh.getNumberOfCells(); + _Nnodes = _mesh.getNumberOfNodes(); + + cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<close(); + throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); + } + } + else if(_Ndim==3) + { + if( _mesh.isTetrahedral() )//Mesh dim=3 + cout<<"3D Finite element method on tetrahedra"<close(); + throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); + } + } + + _VV=Field ("Temperature", NODES, _mesh, 1); + + _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES); + _boundaryNodeIds = _mesh.getBoundaryNodeIds(); + _NboundaryNodes=_boundaryNodeIds.size(); + } } - - _mesh=M; - _Nmailles = _mesh.getNumberOfCells(); - _Nnodes = _mesh.getNumberOfNodes(); - - cout<<"Mesh has "<< _Nmailles << " cells and " << _Nnodes << " nodes"<close(); - throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); - } - } - else if(_Ndim==3) - { - if( _mesh.isTetrahedral() )//Mesh dim=3 - cout<<"3D Finite element method on tetrahedra"<close(); - throw CdmathException("StationaryDiffusionEquation::setMesh: mesh has incorrect cell types"); - } - } - - _VV=Field ("Temperature", NODES, _mesh, 1); - - _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES); - _boundaryNodeIds = _mesh.getBoundaryNodeIds(); - _NboundaryNodes=_boundaryNodeIds.size(); - } + + /* MPI distribution parameters */ + int nbVoisinsMax;//Mettre en attribut ? + if(!_FECalculation){ + MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD); + MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD); + nbVoisinsMax = _neibMaxNbCells; + } + else{ + MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD); + MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD); + nbVoisinsMax = _neibMaxNbNodes; + } + _d_nnz = (nbVoisinsMax+1)*_nVar; + _o_nnz = nbVoisinsMax *_nVar; _meshSet=true; } @@ -758,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; @@ -769,10 +781,10 @@ void StationaryDiffusionEquation::setLinearSolver(linearSolver kspType, precondi 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; + cout << "!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!" << endl; + *_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 !!!" ); } } @@ -787,36 +799,36 @@ bool StationaryDiffusionEquation::solveStationaryProblem() bool stop=false; // Does the Problem want to stop (error) ? bool converged=false; // has the newton scheme converged (end) ? - cout<< "!!! Running test case "<< _fileName << " using "; + PetscPrintf(PETSC_COMM_WORLD,"!!! Running test case %s using ",_fileName); *_runLogFile<< "!!! Running test case "<< _fileName<< " using "; if(!_FECalculation) { - cout<< "Finite volumes method"<close(); throw CdmathException("Failed computing diffusion matrix"); } computeRHS(stop);//For the moment the heat power does not depend on the unknown temperature (linear RHS) if (stop){ - cout << "Error : failed computing right hand side, stopping calculation"<< endl; + PetscPrintf(PETSC_COMM_WORLD,"Error : failed computing right hand side, stopping calculation\n"); *_runLogFile << "Error : failed computing right hand side, stopping calculation"<< endl; throw CdmathException("Failed computing right hand side"); } stop = iterateNewtonStep(converged); if (stop){ - cout << "Error : failed solving linear system, stopping calculation"<< endl; + PetscPrintf(PETSC_COMM_WORLD,"Error : failed solving linear system, stopping calculation\n"); *_runLogFile << "Error : failed linear system, stopping calculation"<< endl; _runLogFile->close(); throw CdmathException("Failed solving linear system"); @@ -858,7 +870,7 @@ bool StationaryDiffusionEquation::solveStationaryProblem() } void StationaryDiffusionEquation::save(){ - cout<< "Saving numerical results"<1){ + VecScatterBegin(_scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD); + VecScatterEnd( _scat,_Tk,_Tk_seq,INSERT_VALUES,SCATTER_FORWARD); + } - _VV.setInfoOnComponent(0,"Temperature_(K)"); - switch(_saveFormat) - { - case VTK : - _VV.writeVTK(resultFile); - break; - case MED : - _VV.writeMED(resultFile); - break; - case CSV : - _VV.writeCSV(resultFile); - break; - } -} -Field -StationaryDiffusionEquation::getOutputTemperatureField() -{ - if(!_computationCompletedSuccessfully) - throw("Computation not performed yet or failed. No temperature field available"); - else - return _VV; + if(_verbose or _system) + VecView(_Tk,PETSC_VIEWER_STDOUT_WORLD); + + if(_mpi_rank==0){ + double Ti; + if(!_FECalculation) + for(int i=0; i<_Nmailles; i++) + { + if(_mpi_size>1) + VecGetValues(_Tk_seq, 1, &i, &Ti); + else + VecGetValues(_Tk , 1, &i, &Ti); + _VV(i)=Ti; + } + else + { + int globalIndex; + for(int i=0; i<_NunknownNodes; i++) + { + if(_mpi_size>1) + VecGetValues(_Tk_seq, 1, &i, &Ti); + else + VecGetValues(_Tk , 1, &i, &Ti); + globalIndex = DiffusionEquation::globalNodeIndex(i, _dirichletNodeIds); + _VV(globalIndex)=Ti; + } + + Node Ni; + string nameOfGroup; + for(int i=0; i<_NdirichletNodes; i++) + { + Ni=_mesh.getNode(_dirichletNodeIds[i]); + nameOfGroup = Ni.getGroupName(); + _VV(_dirichletNodeIds[i])=_limitField[nameOfGroup].T; + } + } + + _VV.setInfoOnComponent(0,"Temperature_(K)"); + switch(_saveFormat) + { + case VTK : + _VV.writeVTK(resultFile); + break; + case MED : + _VV.writeMED(resultFile); + break; + case CSV : + _VV.writeCSV(resultFile); + break; + } + } } void StationaryDiffusionEquation::terminate() @@ -930,6 +947,15 @@ void StationaryDiffusionEquation::terminate() VecDestroy(&_deltaT); VecDestroy(&_b); MatDestroy(&_A); + if(_mpi_size>1 && _mpi_rank == 0) + VecDestroy(&_Tk_seq); + + PetscBool petscInitialized; + PetscInitialized(&petscInitialized); + if(petscInitialized) + PetscFinalize(); + + delete _runLogFile; } void StationaryDiffusionEquation::setDirichletValues(map< int, double> dirichletBoundaryValues) @@ -967,7 +993,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); } } @@ -992,12 +1018,161 @@ StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, doubl MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol); Field my_eigenfield; - if(_FECalculation) - my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev); - else - my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev); + my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev); my_eigenfield.setFieldByDataArrayDouble(d); return my_eigenfield; } + +Field& +StationaryDiffusionEquation::getOutputTemperatureField() +{ + if(!_computationCompletedSuccessfully) + throw("Computation not performed yet or failed. No temperature field available"); + else + return _VV; +} + +Field& +StationaryDiffusionEquation::getRodTemperatureField() +{ + return getOutputTemperatureField(); +} + +vector +StationaryDiffusionEquation::getInputFieldsNames() +{ + vector result(2); + + result[0]="FluidTemperature"; + result[1]="HeatPower"; + + return result; +} +vector +StationaryDiffusionEquation::getOutputFieldsNames() +{ + vector result(1); + + result[0]="RodTemperature"; + + return result; +} + +Field& +StationaryDiffusionEquation::getOutputField(const string& nameField ) +{ + if(nameField=="RodTemperature" || nameField=="RODTEMPERATURE" || nameField=="TEMPERATURECOMBUSTIBLE" || nameField=="TemperatureCombustible" ) + return getRodTemperatureField(); + else + { + cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl; + throw CdmathException("DiffusionEquation::getOutputField error : Unknown Field name"); + } +} + +void +StationaryDiffusionEquation::setInputField(const string& nameField, Field& inputField ) +{ + if(!_meshSet) + throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setInputField set the mesh first"); + + if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" || nameField=="TemperatureFluide" || nameField=="TEMPERATUREFLUIDE") + return setFluidTemperatureField( inputField) ; + else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" ) + return setHeatPowerField( inputField ); + else + { + cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl; + throw CdmathException("StationaryDiffusionEquation::setInputField error : Unknown Field name"); + } +} + +void +StationaryDiffusionEquation::setFluidTemperatureField(Field coupledTemperatureField){ + if(!_meshSet) + throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setFluidTemperatureField set initial field first"); + + coupledTemperatureField.getMesh().checkFastEquivalWith(_mesh); + _fluidTemperatureField=coupledTemperatureField; + _fluidTemperatureFieldSet=true; +}; + +void +StationaryDiffusionEquation::setHeatPowerField(Field heatPower){ + if(!_meshSet) + throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first"); + + heatPower.getMesh().checkFastEquivalWith(_mesh); + _heatPowerField=heatPower; + _heatPowerFieldSet=true; +} + +void +StationaryDiffusionEquation::setHeatPowerField(string fileName, string fieldName, int iteration, int order, int meshLevel){ + if(!_meshSet) + throw CdmathException("!!!!!!!! StationaryDiffusionEquation::setHeatPowerField set initial field first"); + + _heatPowerField=Field(fileName, CELLS,fieldName, iteration, order, meshLevel); + _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"<