[dune-pdelab] Problem with UGGrid-GMSH physical-groups in parallel version

Jö Fahlke jorrit at jorrit.de
Mon Oct 24 16:45:55 CEST 2016


There were utility programs in dune-grid which showed how to do this.  The got
removed, but they are still in the history:
https://gitlab.dune-project.org/core/dune-grid/blob/ec8963acda2daa90aea5e48b9b5835c84307da10/src/gmsh-to-alu/main.hh#L155

They do use alugrid instead of UG, but the actual load-balance and
data-transfer part should work identical.  There are differences in how you
create the grid as Christian mentioned and I don't quite remmember which
possibility you need for UG: either you call GmshReader on all ranks and it
will just give you an empty grid on any rank > 0; or you call GmshReader on
rank 0 only and create an empty grid directly with some constructor on the
other ranks.  You'll have to find out which one applies here for yourself, or
maybe someone else can tell you?

Regards,
Jö.

Am Mon, 24. Oct 2016, 14:23:24 +0200 schrieb Shubhangi Gupta:
> Date: Mon, 24 Oct 2016 14:23:24 +0200
> From: Shubhangi Gupta <sgupta at geomar.de>
> To: Christian Engwer <christian.engwer at uni-muenster.de>
> Cc: dune-pdelab mailing list <dune-pdelab at dune-project.org>
> Subject: Re: [dune-pdelab] Problem with UGGrid-GMSH physical-groups in
>  parallel version
> X-No-Auth: unauthenticated sender
> X-No-Relay: not in my network
> User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:45.0) Gecko/20100101
>  Thunderbird/45.3.0
> X-Envelope-From: <sgupta at geomar.de>
> 
> Hi Christian,
> 
> Thanks a lot for your answer.
> 
> Sorry , but I didn't understand exactly what I should do...
> So, something like this?
> 
>         typedef Dune::UGGrid<dim> GridType;
>         GridType grid;
> 
>         // read gmsh file
>         Dune::GridFactory<GridType> factory(&grid);
>         std::vector<int> boundary_index_map;
>         std::vector<int> element_index_map;
>         const std::string grid_file = parameters.getGridName() ;
>         if(helper.rank()==0){
> Dune::GmshReader<GridType>::read(factory,grid_file,true,true);
>             factory.createGrid();
>             grid.loadBalance();
>         }
>         // get view
>         typedef GridType::LeafGridView GV;
>         const GV& gv=grid.leafGridView();
>         using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
> 
> 
> And how can I communicate data along the load balance call?
> Is there an example in pdelab howto that I can refer to?
> 
> Thanks once again,
> warm regards, Shubhangi
> 
> 
> On 24.10.2016 13:47, Christian Engwer wrote:
> >There are two things you have to take care of
> >
> >a) you use UGGrid. UGGrid must be constructed on rank 0 and the you
> >    call load-balance to distribute the grid.
> >b) you collect data from the gmsh file. After load-balance indices
> >    change and thus you have to take case that the correct material
> >    properties and the boundary conditions are available on the remote
> >    rank. This gives two options:
> >    - you can communicate the data along the loadbalance call
> >    - you can store the data initially w.r.t. the global ID and not the
> >      local index, e.g. in a map. Every processor can read the file,
> >      store the association and then we have all required information
> >      locally available and remap onto a local index.
> >
> >Christian
> >
> >On Mon, Oct 24, 2016 at 11:53:46AM +0200, Shubhangi Gupta wrote:
> >>Hi all,
> >>I make grid for my problem using UGGrid with a gmsh file.
> >>I also prescribe the material properties and the boundary conditions using
> >>'physical groups'. I have written the functions below (M and BCType).
> >>These worked fine in the sequential version, but I am now getting
> >>segmentation fault for parallel version.
> >>I use the NonOverlappingEntitySet (similar to the
> >>convection-diffusion/nonoverlappingsinglephaseflow example in
> >>pdelab-how-to.)
> >>Here is how I create the grid:
> >>
> >>         typedef Dune::UGGrid<dim> GridType;
> >>         GridType grid;
> >>         // read gmsh file
> >>         Dune::GridFactory<GridType> factory(&grid);
> >>         std::vector<int> boundary_index_map;
> >>         std::vector<int> element_index_map;
> >>         const std::string grid_file = parameters.getGridName() ;
> >>Dune::GmshReader<GridType>::read(factory,grid_file,true,true);
> >>         factory.createGrid();
> >>         grid.loadBalance();
> >>         std::cout << " after load balance /" << helper.rank() << "/ " <<
> >>grid.size(0) << std::endl;
> >>         //grid.globalRefine(1);
> >>         // std::cout << " after refinement /" << helper.rank() << "/ " <<
> >>grid.size(0) << std::endl;
> >>         // get view
> >>         typedef GridType::LeafGridView GV;
> >>         const GV& gv=grid.leafGridView();
> >>         using ES = Dune::PDELab::NonOverlappingEntitySet<GV>;
> >>         ES es(gv);
> >>         problem_driver(es);
> >>
> >>I am not sure why I am getting the seg-fault.
> >>I will really appreciate if anyone can help me with this issue!
> >>Thanks a lot in advance,
> >>warm regards,
> >>Shubhangi
> >>
> >>PS: Here are the functions M and BCType for reference:
> >>___________________________________________________________________________________________________________________
> >>/** \brief A function for material properties */
> >>template<typename GV, typename RF, typename PGMap>
> >>class M : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1>
> >>>, M<GV,RF,PGMap> >
> >>{
> >>private:
> >>       const GV& gv;
> >>       const Dune::MultipleCodimMultipleGeomTypeMapper<GV, P0Layout> mapper ;
> >>       const PGMap& pg ;
> >>public:
> >>   typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> >
> >>Traits;
> >>   typedef Dune::PDELab::GridFunctionBase < Dune::PDELab::GridFunctionTraits<
> >>GV, RF, 1 , Dune::FieldVector< RF,1 > >, M< GV, RF, PGMap > > BaseT ;
> >>
> >>   //! construct from grid view
> >>   M ( const GV& gv_ , const PGMap& pg_ ) : gv( gv_ ) , mapper( gv_ ) , pg(
> >>pg_ ) {}
> >>
> >>   //! evaluate extended function on element
> >>   inline void evaluate (const typename Traits::ElementType& e,
> >>                                     const typename Traits::DomainType&
> >>xlocal,
> >>                                     typename Traits::RangeType& y) const {
> >>
> >>       // retrieve element index and corresponding material index on level 0
> >>           typename GV::template Codim<0>::EntityPointer ep(e) ;
> >>           while ( ep->level() != 0 ) ep = ep->father() ;
> >>           const int ei = mapper.map(ep) ;
> >>           const int physgroup_index = pg[ ei ] ;
> >>
> >>       // evaluate physical group map and set values accordingly
> >>           switch ( physgroup_index ){
> >>                 case 1     : y = y=1.     ; break ;
> >>                 case 2     : y = y=0.     ; break ;
> >>                 default     : y = 1.         ; break ; // only one material
> >>here
> >>           }
> >>           return ;
> >>   }
> >>
> >>   //! get a reference to the grid view
> >>   inline const GV& getGridView () {return gv;}
> >>};
> >>___________________________________________________________________________________________________________________
> >>/** \brief A function for boundary types */
> >>template<typename GV, typename PGMap >
> >>class BCType : public Dune::PDELab::BoundaryGridFunctionBase<Dune::PDELab::BoundaryGridFunctionTraits<GV,int,1,Dune::FieldVector<int,1>
> >>>,BCType<GV,PGMap> >
> >>{
> >>public :
> >>     typedef Dune::PDELab::BoundaryGridFunctionTraits<
> >>GV,int,1,Dune::FieldVector<int,1> > Traits ;
> >>
> >>     // ! construct from gridview
> >>     BCType ( const GV& gv_ , const PGMap& pgmap_ ) : gv ( gv_ ) , pgmap (
> >>pgmap_ ) {}
> >>     // ! return bctype at point on intersection
> >>     template<typename I>
> >>     inline void evaluate ( I& i , const typename Traits::DomainType& xlocal
> >>, typename Traits::RangeType& y ) const {
> >>
> >>          // evaluate with maps
> >>          int physgroup_index = pgmap[
> >>i.intersection().boundarySegmentIndex() ] ;
> >>          switch ( physgroup_index )
> >>          {
> >>                  case 3    : y = Indices::BCId_neumann;    break ;
> >>                  case 4    : y = Indices::BCId_dirichlet;      break ;
> >>              default : y = Indices::BCId_neumann ; break ; // Neumann
> >>          }
> >>          return ;
> >>     }
> >>
> >>     // ! get a reference to the gridview
> >>     inline const GV& getGridView () { return gv ; }
> >>
> >>private :
> >>     const GV& gv ;
> >>     const PGMap& pgmap ;
> >>
> >>} ;
> >>_________________________________________________________________________________________________________________
> >>
> >>_______________________________________________
> >>dune-pdelab mailing list
> >>dune-pdelab at dune-project.org
> >>http://lists.dune-project.org/mailman/listinfo/dune-pdelab
> >
> 

> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-pdelab


-- 
Jorrit (Jö) Fahlke, Institute for Computational und Applied Mathematics,
University of Münster, Orleans-Ring 10, D-48149 Münster
Tel: +49 251 83 35146 Fax: +49 251 83 32729

Of all the things I've lost, I miss my mind the most.
-- Ozzy Osbourne
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 811 bytes
Desc: Digital signature
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20161024/b7497dc7/attachment.sig>


More information about the dune-pdelab mailing list