Sierra Toolkit Version of the Day
UseCase_Rebal_4.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 #include <use_cases/UseCase_Rebal_4.hpp>
00010 
00011 #include <stk_util/parallel/Parallel.hpp>
00012 #include <stk_util/parallel/ParallelReduce.hpp>
00013 
00014 #include <stk_mesh/base/FieldData.hpp>
00015 #include <stk_mesh/base/GetEntities.hpp>
00016 #include <stk_mesh/base/MetaData.hpp>
00017 #include <stk_mesh/base/BulkData.hpp>
00018 
00019 #include <stk_mesh/fem/CoordinateSystems.hpp>
00020 #include <stk_mesh/fem/FEMMetaData.hpp>
00021 #include <stk_mesh/fem/FEMHelpers.hpp>
00022 #include <stk_mesh/fem/CreateAdjacentEntities.hpp>
00023 
00024 #include <stk_mesh/fixtures/HexFixture.hpp>
00025 
00026 #include <stk_rebalance/Rebalance.hpp>
00027 #include <stk_rebalance/Partition.hpp>
00028 #include <stk_rebalance/ZoltanPartition.hpp>
00029 
00030 #include <stk_rebalance_utils/RebalanceUtils.hpp>
00031 
00032 //----------------------------------------------------------------------
00033 
00034 using namespace stk::mesh::fixtures;
00035 
00036 typedef stk::mesh::Field<double> ScalarField ;
00037 
00038 namespace stk {
00039 namespace rebalance {
00040 namespace use_cases {
00041 
00042 class GreedySideset : public Partition {
00043   public :
00044   struct MeshInfo {
00045     std::vector<mesh::Entity *>      mesh_entities;
00046     const VectorField * nodal_coord_ref ;
00047     const ScalarField * elem_weight_ref;
00048     std::vector<unsigned>            dest_proc_ids ;
00049 
00051     MeshInfo():
00052       nodal_coord_ref(NULL),
00053       elem_weight_ref(NULL) {}
00054 
00056     ~MeshInfo() {}
00057   };
00058   explicit GreedySideset(ParallelMachine pm,
00059                          const stk::mesh::PartVector & surfaces,
00060                          mesh::BulkData   & bulk_data);
00061   virtual ~GreedySideset();
00062   virtual void reset_dest_proc_data();
00063   virtual void set_mesh_info ( const std::vector<mesh::Entity *> &mesh_entities,
00064                                const VectorField   * nodal_coord_ref,
00065                                const ScalarField   * elem_weight_ref=NULL);
00066   virtual void determine_new_partition(bool &RebalancingNeeded);
00067   virtual unsigned num_elems() const;
00068   virtual int get_new_partition(stk::mesh::EntityProcVec &new_partition);
00069   virtual bool partition_dependents_needed()const;
00070   bool find_mesh_entity(const mesh::Entity * entity, unsigned & moid) const;
00071   unsigned destination_proc(const unsigned moid) const;
00072   void set_destination_proc(const unsigned moid, const unsigned proc );
00073   MeshInfo  mesh_information_;
00074   unsigned  total_number_entities_;
00075   const stk::mesh::PartVector & surfaces_;
00076   mesh::BulkData   & bulk_data_;
00077 };
00078 
00079 GreedySideset::GreedySideset(ParallelMachine pm,
00080                              const stk::mesh::PartVector & surfaces,
00081                              mesh::BulkData   & bulk_data) :
00082   stk::rebalance::Partition(pm),
00083   mesh_information_(),
00084   surfaces_(surfaces),
00085   bulk_data_(bulk_data) {}
00086 GreedySideset::~GreedySideset() {}
00087 void GreedySideset::reset_dest_proc_data() {
00088   const int proc = parallel_machine_rank(comm_);
00089   const unsigned size = mesh_information_.mesh_entities.size();
00090   mesh_information_.dest_proc_ids.assign(size, proc);
00091 }
00092 void GreedySideset::set_mesh_info ( const std::vector<mesh::Entity *> &mesh_entities,
00093                                     const VectorField   * nodal_coord_ref,
00094                                     const ScalarField   * elem_weight_ref){
00095   MeshInfo mesh_info;
00096 
00097   /* Keep track of the total number of elements. */
00098   total_number_entities_ = mesh_entities.size();
00099 
00100   mesh_info.mesh_entities = mesh_entities;
00101   mesh_info.nodal_coord_ref = nodal_coord_ref;
00102   mesh_info.elem_weight_ref = elem_weight_ref;
00103 
00109   mesh_info.dest_proc_ids.assign(mesh_entities.size(), stk::parallel_machine_rank(comm_));
00110 
00111   mesh_information_ = mesh_info;
00112 }
00113 
00114 unsigned GreedySideset::num_elems() const {return total_number_entities_ ;}
00115 int GreedySideset::get_new_partition(stk::mesh::EntityProcVec &new_partition){
00116 std::vector<mesh::Entity*>::iterator i=mesh_information_.mesh_entities.begin();
00117 std::vector<unsigned>     ::iterator j=mesh_information_.dest_proc_ids.begin();
00118   for (;i != mesh_information_.mesh_entities.end(),
00119         j != mesh_information_.dest_proc_ids.end();
00120         ++i,++j) {
00121     mesh::Entity * mesh_entity = *i;
00122     unsigned proc = *j;
00123     mesh::EntityProc et(mesh_entity, proc);
00124     new_partition.push_back(et);
00125   }
00126   return 0;
00127 }
00128 bool GreedySideset::partition_dependents_needed()const{return true;}
00129 
00130 bool GreedySideset::find_mesh_entity(const mesh::Entity * entity, unsigned & moid) const
00131 {
00132   unsigned len = mesh_information_.mesh_entities.size();
00133   for(moid = 0; moid < len; ++moid)
00134   {
00135     if(mesh_information_.mesh_entities[moid] == entity) return true;
00136   }
00137   return false;
00138 }
00139 unsigned GreedySideset::destination_proc(const unsigned moid) const
00140 {
00141   return mesh_information_.dest_proc_ids[ moid ];
00142 }
00143 void GreedySideset::set_destination_proc(const unsigned moid,
00144                                          const unsigned proc )
00145 {
00146   mesh_information_.dest_proc_ids[ moid ] = proc;
00147 }
00148 
00149 
00150 void GreedySideset::determine_new_partition(bool &RebalancingNeeded) {
00151 
00152   reset_dest_proc_data();
00153 
00154   stk::mesh::fem::FEMMetaData & fem_meta = stk::mesh::fem::FEMMetaData::get(bulk_data_);
00155   const stk::mesh::EntityRank side_rank = fem_meta.side_rank();
00156   const stk::mesh::EntityRank elem_rank = fem_meta.element_rank();
00157 
00158   // Select active ghosted side faces.
00159   stk::mesh::Selector selector(!fem_meta.locally_owned_part() &
00160                                 stk::mesh::selectIntersection(surfaces_));
00161 
00162   mesh::EntityVector sides;
00163   mesh::get_selected_entities(selector, bulk_data_.buckets(side_rank), sides);
00164 
00165   const unsigned p_rank = bulk_data_.parallel_rank();
00166   size_t local_changes = 0;
00167   const unsigned nSide = sides.size();
00168   for(unsigned iSide = 0; iSide < nSide; ++iSide)
00169   {
00170     const mesh::Entity & side = *sides[iSide];
00171     const unsigned sideProc = side.owner_rank();
00172     ThrowRequireMsg(sideProc!=p_rank,
00173      "When iterating Non-locally owned sides, found a locally owned side.");
00174 
00175     stk::mesh::PairIterRelation iElem = side.relations(elem_rank);
00176     for ( ; iElem.first != iElem.second; ++iElem.first ) {
00177       const mesh::Entity & elem = *iElem.first->entity();
00178       unsigned moid;
00179       const bool mesh_entity_found = find_mesh_entity(&elem, moid);
00180       if (mesh_entity_found) {
00181         const unsigned elemProc = elem.owner_rank();
00182         ThrowRequireMsg(elemProc==p_rank,
00183           "When iterating locally owned elements, found a non-locally owned element.");
00184         const unsigned destProc = destination_proc(moid);
00185         ThrowRequireMsg(destProc==p_rank || destProc==sideProc,
00186          " Sanity check failed: "
00187          "It's possible that an element is connected to "
00188          "two sides that are owned by different procs.  We don't "
00189          "yet handle that situation here but we can, at least, "
00190          "detect it. ");
00191         if(elemProc != sideProc)
00192         {
00193           ++local_changes;
00194           set_destination_proc(moid, sideProc);
00195         }
00196       }
00197     }
00198   }
00199   size_t global_changes = 0;
00200   stk::all_reduce_sum (comm_, &local_changes, &global_changes, 1);
00201   RebalancingNeeded = global_changes > 0;
00202 }
00203 
00204 enum { nx = 3, ny = 3 };
00205 
00206 bool test_greedy_sideset ( stk::ParallelMachine comm )
00207 {
00208   unsigned spatial_dimension = 2;
00209   std::vector<std::string> rank_names = stk::mesh::fem::entity_rank_names(spatial_dimension);
00210   stk::mesh::fem::FEMMetaData fem_meta;
00211   fem_meta.FEM_initialize(spatial_dimension, rank_names);
00212   stk::mesh::MetaData & meta_data = stk::mesh::fem::FEMMetaData::get_meta_data(fem_meta);
00213   stk::mesh::BulkData bulk_data( meta_data , comm , 100 );
00214   const stk::mesh::EntityRank element_rank    = fem_meta.element_rank();
00215   const stk::mesh::EntityRank node_rank       = fem_meta.node_rank();
00216 
00217   stk::mesh::fem::CellTopology quad_top(shards::getCellTopologyData<shards::Quadrilateral<4> >());
00218   stk::mesh::fem::CellTopology line_top(shards::getCellTopologyData<shards::Line<2> >());
00219   stk::mesh::Part & quad_part( fem_meta.declare_part("quad", quad_top ) );
00220   stk::mesh::Part & side_part( fem_meta.declare_part("line", line_top ) );
00221   VectorField & coord_field( fem_meta.declare_field< VectorField >( "coordinates" ) );
00222   ScalarField & weight_field( fem_meta.declare_field< ScalarField >( "element_weights" ) );
00223 
00224   stk::mesh::put_field( coord_field , node_rank , fem_meta.universal_part() );
00225   stk::mesh::put_field(weight_field , element_rank , fem_meta.universal_part() );
00226 
00227   fem_meta.commit();
00228   const unsigned p_rank = bulk_data.parallel_rank();
00229   bulk_data.modification_begin();
00230 
00231   if ( !p_rank ) {
00232 
00233     std::vector<std::vector<stk::mesh::Entity*> > quads(nx);
00234     for ( unsigned ix = 0 ; ix < nx ; ++ix ) quads[ix].resize(ny);
00235 
00236     const unsigned nnx = nx + 1 ;
00237     for ( unsigned iy = 0 ; iy < ny ; ++iy ) {
00238       for ( unsigned ix = 0 ; ix < nx ; ++ix ) {
00239         stk::mesh::EntityId elem = 1 + ix + iy * nx ;
00240         stk::mesh::EntityId nodes[4] ;
00241         nodes[0] = 1 + ix + iy * nnx ;
00242         nodes[1] = 2 + ix + iy * nnx ;
00243         nodes[2] = 2 + ix + ( iy + 1 ) * nnx ;
00244         nodes[3] = 1 + ix + ( iy + 1 ) * nnx ;
00245 
00246         stk::mesh::Entity &q = stk::mesh::fem::declare_element( bulk_data , quad_part , elem , nodes );
00247         quads[ix][iy] = &q;
00248       }
00249     }
00250 
00251     for ( unsigned iy = 0 ; iy < ny ; ++iy ) {
00252       for ( unsigned ix = 0 ; ix < nx ; ++ix ) {
00253         stk::mesh::EntityId elem = 1 + ix + iy * nx ;
00254         stk::mesh::Entity * e = bulk_data.get_entity( element_rank, elem );
00255         double * const e_weight = stk::mesh::field_data( weight_field , *e );
00256         *e_weight = 1.0;
00257       }
00258     }
00259     for ( unsigned iy = 0 ; iy <= ny ; ++iy ) {
00260       for ( unsigned ix = 0 ; ix <= nx ; ++ix ) {
00261         stk::mesh::EntityId nid = 1 + ix + iy * nnx ;
00262         stk::mesh::Entity * n = bulk_data.get_entity( node_rank, nid );
00263         double * const coord = stk::mesh::field_data( coord_field , *n );
00264         coord[0] = .1*ix;
00265         coord[1] = .1*iy;
00266         coord[2] = 0;
00267       }
00268     }
00269   }
00270 
00271   bulk_data.modification_end();
00272 
00273   // create some sides and faces to rebalance.
00274   stk::mesh::PartVector add_parts;
00275   stk::mesh::create_adjacent_entities(bulk_data, add_parts);
00276 
00277   bulk_data.modification_begin();
00278 
00279   const stk::mesh::PartVector surfaces(1, &side_part);
00280   {
00281     const stk::mesh::PartVector empty_remove_parts;
00282     stk::mesh::fem::FEMMetaData & fmeta = stk::mesh::fem::FEMMetaData::get(bulk_data);
00283     const stk::mesh::EntityRank side_rank = fmeta.side_rank();
00284     stk::mesh::Selector selector2( fmeta.locally_owned_part());
00285     mesh::EntityVector sides;
00286     mesh::get_selected_entities(selector2, bulk_data.buckets(side_rank), sides);
00287 
00288     const unsigned nSide = sides.size();
00289     for(unsigned iSide = 0; iSide < nSide; ++iSide)
00290     {
00291       mesh::Entity & side = *sides[iSide];
00292       if (side.identifier()==7) {
00293         bulk_data.change_entity_parts(side, surfaces, empty_remove_parts);
00294       }
00295     }
00296   }
00297   bulk_data.modification_end();
00298 
00299   // Zoltan partition is specialized form a virtual base class, stk::rebalance::Partition.
00300   // Other specializations are possible.
00301   Teuchos::ParameterList emptyList;
00302   stk::rebalance::Zoltan zoltan_partition(comm, spatial_dimension, emptyList);
00303   stk::mesh::Selector selector3(fem_meta.locally_owned_part());
00304   stk::rebalance::rebalance(bulk_data, selector3, &coord_field, NULL, zoltan_partition);
00305   {
00306     const int  print_stats = 1;
00307     int        nentity     = 0;
00308     double     entity_wgt  = 0;
00309     int        ncuts       = 0;
00310     double     cut_wgt     = 0;
00311     int        nboundary   = 0;
00312     int        nadj        = 0;
00313     const int ierr = zoltan_partition.evaluate (print_stats, &nentity, &entity_wgt, &ncuts, &cut_wgt, &nboundary, &nadj);
00314     std::cout <<" Information returned from the Zoltan evaluate function:"<<std::endl;
00315     std::cout <<" Error Code:             :"<<ierr      <<std::endl;
00316     std::cout <<" Number of entities:     :"<<nentity   <<std::endl;
00317     std::cout <<" Number of cuts:         :"<<ncuts     <<std::endl;
00318     std::cout <<" Cut Weight:             :"<<cut_wgt   <<std::endl;
00319     std::cout <<" Number on Boundary:     :"<<nboundary <<std::endl;
00320     std::cout <<" Number Adjancent:       :"<<nadj      <<std::endl;
00321     {
00322       stk::mesh::fem::FEMMetaData & fmeta = stk::mesh::fem::FEMMetaData::get(bulk_data);
00323       const stk::mesh::EntityRank side_rank = fmeta.side_rank();
00324       const stk::mesh::EntityRank elem_rank = fmeta.element_rank();
00325       const mesh::Entity *s = bulk_data.get_entity(side_rank,7);
00326       if (s) {
00327         const mesh::Entity & side = *s;
00328         if (p_rank == side.owner_rank()) {
00329           stk::mesh::PairIterRelation iElem = side.relations(elem_rank);
00330           for ( ; iElem.first != iElem.second; ++iElem.first ) {
00331             const mesh::Entity & elem = *iElem.first->entity();
00332             const unsigned elemProc = elem.owner_rank();
00333             if (elemProc!=p_rank) {
00334               std::cout <<p_rank<<" Good: Found element of of side 7 not owned."
00335                         <<" Element "<<elemProc
00336                         <<" on processor "<<p_rank
00337                         <<" owned by "<<elemProc<<std::endl;
00338             }
00339           }
00340         }
00341       }
00342     }
00343   }
00344 
00345   const double imbalance_threshold = 1.5;
00346   bool do_rebal = imbalance_threshold < stk::rebalance::check_balance(bulk_data, NULL, element_rank);
00347 
00348   if( !p_rank )
00349     std::cout << std::endl
00350      << "Use Case 4: imbalance_threshold after rebalance 1 = " << imbalance_threshold <<", "<<do_rebal << std::endl;
00351 
00352   {
00353     stk::rebalance::use_cases::GreedySideset greedy_sideset(comm, surfaces, bulk_data);
00354     stk::mesh::Selector selector4(fem_meta.locally_owned_part());
00355     stk::rebalance::rebalance(bulk_data, selector4, &coord_field, NULL, greedy_sideset);
00356   }
00357 
00358   do_rebal = imbalance_threshold < stk::rebalance::check_balance(bulk_data, NULL, element_rank);
00359 
00360   if( !p_rank )
00361     std::cout << std::endl
00362      << "Use Case 4: imbalance_threshold after rebalance 2 = " << imbalance_threshold <<", "<<do_rebal << std::endl;
00363   {
00364     stk::mesh::fem::FEMMetaData & fmeta = stk::mesh::fem::FEMMetaData::get(bulk_data);
00365     const stk::mesh::EntityRank side_rank = fmeta.side_rank();
00366     const stk::mesh::EntityRank elem_rank = fmeta.element_rank();
00367     mesh::Entity *s = bulk_data.get_entity(side_rank,7);
00368     if (s) {
00369       mesh::Entity & side = *s;
00370       if (p_rank == side.owner_rank()) {
00371         stk::mesh::PairIterRelation iElem = side.relations(elem_rank);
00372         for ( ; iElem.first != iElem.second; ++iElem.first ) {
00373           const mesh::Entity & elem = *iElem.first->entity();
00374           const unsigned elemProc = elem.owner_rank();
00375           if (elemProc!=p_rank) {
00376             std::cerr <<p_rank<<" Error: Found element of of side 7 not owned:"<<elemProc<<std::endl;
00377           }
00378           ThrowRequireMsg(elemProc==p_rank,
00379            "Use case 4 error check failed. Found element of of side 7 not owned.");
00380         }
00381       }
00382     }
00383   }
00384   // Check that we satisfy our threshhold
00385   bool result = !do_rebal ;
00386 
00387   // And verify that all dependent entities are on the same proc as their parent element
00388   {
00389     stk::mesh::EntityVector entities;
00390     stk::mesh::Selector selector5 = fem_meta.locally_owned_part();
00391 
00392     get_selected_entities(selector5, bulk_data.buckets(node_rank), entities);
00393     result &= verify_dependent_ownership(element_rank, entities);
00394     get_selected_entities(selector5, bulk_data.buckets(fem_meta.side_rank()), entities);
00395     result &= verify_dependent_ownership(element_rank, entities);
00396   }
00397 
00398 
00399 
00400   return result;
00401 }
00402 
00403 } //namespace use_cases
00404 } //namespace rebalance
00405 } //namespace stk
00406 
00407 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines