Merge pull request #160 from streeve/frozen_points

Enable frozen points for boundaries

examples/mechanics/fragmenting_cylinder.cpp

@@ -110,7 +110,7 @@ void fragmentingCylinderExample( const std::string filename )
     using random_type = Kokkos::Random_XorShift64<exec_space>;
     pool_type pool;
     int seed = 456854;
-    pool.init( seed, particles->n_local );
+    pool.init( seed, particles->localOffset() );
     auto init_functor = KOKKOS_LAMBDA( const int pid )
     {
         // Density

src/CabanaPD_BodyTerm.hpp

@@ -26,12 +26,14 @@ struct BodyTerm
 {
     UserFunctor _user_functor;
     bool _force_update;
+    bool _update_frozen;
 
     Timer _timer;
 
-    BodyTerm( UserFunctor user, const bool force )
+    BodyTerm( UserFunctor user, const bool force, const bool update_frozen )
         : _user_functor( user )
         , _force_update( force )
+        , _update_frozen( update_frozen )
     {
     }
 
@@ -41,7 +43,10 @@ struct BodyTerm
     void apply( ExecSpace, ParticleType& particles, const double time )
     {
         _timer.start();
-        Kokkos::RangePolicy<ExecSpace> policy( 0, particles.n_local );
+        std::size_t start = particles.frozenOffset();
+        if ( _update_frozen )
+            start = 0;
+        Kokkos::RangePolicy<ExecSpace> policy( start, particles.localOffset() );
         auto user = _user_functor;
         Kokkos::parallel_for(
             "CabanaPD::BodyTerm::apply", policy,
@@ -56,9 +61,10 @@ struct BodyTerm
 };
 
 template <class UserFunctor>
-auto createBodyTerm( UserFunctor user_functor, const bool force_update )
+auto createBodyTerm( UserFunctor user_functor, const bool force_update,
+                     const bool update_frozen = false )
 {
-    return BodyTerm<UserFunctor>( user_functor, force_update );
+    return BodyTerm<UserFunctor>( user_functor, force_update, update_frozen );
 }
 
 } // namespace CabanaPD

src/CabanaPD_Boundary.hpp

@@ -140,7 +140,7 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary<GeometryType>>
         _timer.start();
 
         _view = index_view_type( "boundary_indices",
-                                 particles.n_local * initial_guess );
+                                 particles.localOffset() * initial_guess );
         _count = index_view_type( "count", 1 );
 
         for ( RegionBoundary plane : planes )
@@ -169,7 +169,8 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary<GeometryType>>
         auto index_space = _view;
         auto count = _count;
         auto x = particles.sliceReferencePosition();
-        Kokkos::RangePolicy<ExecSpace> policy( 0, particles.n_local );
+        // TODO: configure including frozen particles.
+        Kokkos::RangePolicy<ExecSpace> policy( 0, particles.localOffset() );
         auto index_functor = KOKKOS_LAMBDA( const std::size_t pid )
         {
             if ( region.inside( x, pid ) )

src/CabanaPD_Comm.hpp

@@ -270,9 +270,9 @@ class Comm<ParticleType, PMB, TemperatureIndependent>
         halo_ids.rebuild( positions );
 
         // Create the Cabana Halo.
-        halo = std::make_shared<halo_type>( local_grid->globalGrid().comm(),
-                                            particles.n_local, halo_ids._ids,
-                                            halo_ids._destinations, topology );
+        halo = std::make_shared<halo_type>(
+            local_grid->globalGrid().comm(), particles.localOffset(),
+            halo_ids._ids, halo_ids._destinations, topology );
 
         particles.resize( halo->numLocal(), halo->numGhost() );
 

src/CabanaPD_DisplacementProfile.hpp

@@ -56,8 +56,9 @@ void createOutputProfile( MPI_Comm comm, const int num_cell,
             profile( c, 1 ) = user( pid );
         }
     };
+    // TODO: enable ignoring frozen particles.
     Kokkos::RangePolicy<typename memory_space::execution_space> policy(
-        0, particles.n_local );
+        0, particles.localOffset() );
     Kokkos::parallel_for( "displacement_profile", policy, measure_profile );
     auto count_host =
         Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace{}, count );

src/CabanaPD_Force.hpp

@@ -149,23 +149,24 @@ class Force<MemorySpace, BaseForceModel>
     Force( const bool half_neigh, const double delta,
            const ParticleType& particles, const double tol = 1e-14 )
         : _half_neigh( half_neigh )
-        , _neigh_list( neighbor_list_type( particles.sliceReferencePosition(),
-                                           0, particles.n_local, delta + tol,
-                                           1.0, particles.ghost_mesh_lo,
-                                           particles.ghost_mesh_hi ) )
+        , _neigh_list( neighbor_list_type(
+              particles.sliceReferencePosition(), particles.frozenOffset(),
+              particles.localOffset(), delta + tol, 1.0,
+              particles.ghost_mesh_lo, particles.ghost_mesh_hi ) )
     {
     }
 
     // General constructor (necessary for contact, but could be used by any
     // force routine).
     template <class PositionType>
     Force( const bool half_neigh, const double delta,
-           const PositionType& positions, const std::size_t num_local,
-           const double mesh_min[3], const double mesh_max[3],
-           const double tol = 1e-14 )
+           const PositionType& positions, const std::size_t frozen_offset,
+           const std::size_t local_offset, const double mesh_min[3],
+           const double mesh_max[3], const double tol = 1e-14 )
         : _half_neigh( half_neigh )
-        , _neigh_list( neighbor_list_type( positions, 0, num_local, delta + tol,
-                                           1.0, mesh_min, mesh_max ) )
+        , _neigh_list( neighbor_list_type( positions, frozen_offset,
+                                           local_offset, delta + tol, 1.0,
+                                           mesh_min, mesh_max ) )
     {
     }
 
@@ -239,7 +240,6 @@ template <class ForceType, class ParticleType, class ParallelType>
 void computeForce( ForceType& force, ParticleType& particles,
                    const ParallelType& neigh_op_tag, const bool reset = true )
 {
-    auto n_local = particles.n_local;
     auto x = particles.sliceReferencePosition();
     auto u = particles.sliceDisplacement();
     auto f = particles.sliceForce();
@@ -251,14 +251,14 @@ void computeForce( ForceType& force, ParticleType& particles,
 
     // if ( half_neigh )
     // Forces must be atomic for half list
-    // computeForce_half( f_a, x, u, n_local,
+    // computeForce_half( f_a, x, u,
     //                    neigh_op_tag );
 
     // Forces only atomic if using team threading.
     if ( std::is_same<decltype( neigh_op_tag ), Cabana::TeamOpTag>::value )
-        force.computeForceFull( f_a, x, u, particles, n_local, neigh_op_tag );
+        force.computeForceFull( f_a, x, u, particles, neigh_op_tag );
     else
-        force.computeForceFull( f, x, u, particles, n_local, neigh_op_tag );
+        force.computeForceFull( f, x, u, particles, neigh_op_tag );
     Kokkos::fence();
 }
 
@@ -269,7 +269,6 @@ double computeEnergy( ForceType& force, ParticleType& particles,
     double energy = 0.0;
     if constexpr ( is_energy_output<typename ParticleType::output_type>::value )
     {
-        auto n_local = particles.n_local;
         auto x = particles.sliceReferencePosition();
         auto u = particles.sliceDisplacement();
         auto f = particles.sliceForce();
@@ -281,10 +280,9 @@ double computeEnergy( ForceType& force, ParticleType& particles,
 
         // if ( _half_neigh )
         //    energy = computeEnergy_half( force, x, u,
-        //                                  n_local, neigh_op_tag );
+        //                                   neigh_op_tag );
         // else
-        energy = force.computeEnergyFull( W, x, u, particles, n_local,
-                                          neigh_op_tag );
+        energy = force.computeEnergyFull( W, x, u, particles, neigh_op_tag );
         Kokkos::fence();
     }
     return energy;
@@ -296,7 +294,6 @@ template <class ForceType, class ParticleType, class NeighborView,
 void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu,
                    const ParallelType& neigh_op_tag, const bool reset = true )
 {
-    auto n_local = particles.n_local;
     auto x = particles.sliceReferencePosition();
     auto u = particles.sliceDisplacement();
     auto f = particles.sliceForce();
@@ -308,15 +305,14 @@ void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu,
 
     // if ( half_neigh )
     // Forces must be atomic for half list
-    // computeForce_half( f_a, x, u, n_local,
+    // computeForce_half( f_a, x, u,
     //                    neigh_op_tag );
 
     // Forces only atomic if using team threading.
     if ( std::is_same<decltype( neigh_op_tag ), Cabana::TeamOpTag>::value )
-        force.computeForceFull( f_a, x, u, particles, mu, n_local,
-                                neigh_op_tag );
+        force.computeForceFull( f_a, x, u, particles, mu, neigh_op_tag );
     else
-        force.computeForceFull( f, x, u, particles, mu, n_local, neigh_op_tag );
+        force.computeForceFull( f, x, u, particles, mu, neigh_op_tag );
     Kokkos::fence();
 }
 
@@ -328,7 +324,6 @@ double computeEnergy( ForceType& force, ParticleType& particles,
     double energy = 0.0;
     if constexpr ( is_energy_output<typename ParticleType::output_type>::value )
     {
-        auto n_local = particles.n_local;
         auto x = particles.sliceReferencePosition();
         auto u = particles.sliceDisplacement();
         auto f = particles.sliceForce();
@@ -340,9 +335,9 @@ double computeEnergy( ForceType& force, ParticleType& particles,
 
         // if ( _half_neigh )
         //    energy = computeEnergy_half( force, x, u,
-        //                                  n_local, neigh_op_tag );
+        //                                   neigh_op_tag );
         // else
-        energy = force.computeEnergyFull( W, x, u, phi, particles, mu, n_local,
+        energy = force.computeEnergyFull( W, x, u, phi, particles, mu,
                                           neigh_op_tag );
         Kokkos::fence();
     }

src/CabanaPD_HeatTransfer.hpp

@@ -49,10 +49,10 @@ class HeatTransfer : public Force<MemorySpace, BaseForceModel>
 
     template <class TemperatureType, class PosType, class ParticleType,
               class ParallelType>
-    void
-    computeHeatTransferFull( TemperatureType& conduction, const PosType& x,
-                             const PosType& u, const ParticleType& particles,
-                             const int n_local, ParallelType& neigh_op_tag )
+    void computeHeatTransferFull( TemperatureType& conduction, const PosType& x,
+                                  const PosType& u,
+                                  const ParticleType& particles,
+                                  ParallelType& neigh_op_tag )
     {
         _timer.start();
 
@@ -70,7 +70,8 @@ class HeatTransfer : public Force<MemorySpace, BaseForceModel>
                 coeff * ( temp( j ) - temp( i ) ) / xi / xi * vol( j );
         };
 
-        Kokkos::RangePolicy<exec_space> policy( 0, n_local );
+        Kokkos::RangePolicy<exec_space> policy( particles.frozenOffset(),
+                                                particles.localOffset() );
         Cabana::neighbor_parallel_for(
             policy, temp_func, _neigh_list, Cabana::FirstNeighborsTag(),
             neigh_op_tag, "CabanaPD::HeatTransfer::computeFull" );
@@ -86,12 +87,12 @@ class HeatTransfer : public Force<MemorySpace, BaseForceModel>
         const auto rho = particles.sliceDensity();
         const auto conduction = particles.sliceTemperatureConduction();
         auto temp = particles.sliceTemperature();
-        auto n_local = particles.n_local;
         auto euler_func = KOKKOS_LAMBDA( const int i )
         {
             temp( i ) += dt / rho( i ) / model.cp * conduction( i );
         };
-        Kokkos::RangePolicy<exec_space> policy( 0, n_local );
+        Kokkos::RangePolicy<exec_space> policy( particles.frozenOffset(),
+                                                particles.localOffset() );
         Kokkos::parallel_for( "CabanaPD::HeatTransfer::forwardEuler", policy,
                               euler_func );
         _euler_timer.stop();
@@ -104,7 +105,6 @@ void computeHeatTransfer( HeatTransferType& heat_transfer,
                           ParticleType& particles,
                           const ParallelType& neigh_op_tag, const double dt )
 {
-    auto n_local = particles.n_local;
     auto x = particles.sliceReferencePosition();
     auto u = particles.sliceDisplacement();
     auto conduction = particles.sliceTemperatureConduction();
@@ -116,10 +116,10 @@ void computeHeatTransfer( HeatTransferType& heat_transfer,
     // Temperature only needs to be atomic if using team threading.
     if ( std::is_same<decltype( neigh_op_tag ), Cabana::TeamOpTag>::value )
         heat_transfer.computeHeatTransferFull( conduction_a, x, u, particles,
-                                               n_local, neigh_op_tag );
+                                               neigh_op_tag );
     else
         heat_transfer.computeHeatTransferFull( conduction, x, u, particles,
-                                               n_local, neigh_op_tag );
+                                               neigh_op_tag );
     Kokkos::fence();
 
     heat_transfer.forwardEuler( particles, dt );

src/CabanaPD_Integrate.hpp

@@ -106,7 +106,8 @@ class Integrator
             u( i, 1 ) += dt * v( i, 1 );
             u( i, 2 ) += dt * v( i, 2 );
         };
-        Kokkos::RangePolicy<exec_space> policy( 0, v.size() );
+        Kokkos::RangePolicy<exec_space> policy( p.frozenOffset(),
+                                                p.localOffset() );
         Kokkos::parallel_for( "CabanaPD::Integrator::Initial", policy,
                               init_func );
 
@@ -130,7 +131,8 @@ class Integrator
             v( i, 1 ) += half_dt_m * f( i, 1 );
             v( i, 2 ) += half_dt_m * f( i, 2 );
         };
-        Kokkos::RangePolicy<exec_space> policy( 0, v.size() );
+        Kokkos::RangePolicy<exec_space> policy( p.frozenOffset(),
+                                                p.localOffset() );
         Kokkos::parallel_for( "CabanaPD::Integrator::Final", policy,
                               final_func );