Skip to content

Commit

Permalink
Merge pull request #150 from streeve/neighbors_in_force
Browse files Browse the repository at this point in the history 
Move neighbor list into force
  • Loading branch information
@streeve
streeve authored Nov 19, 2024
2 parents 84adbc0 + 4e74a2a commit 5e66739
Show file tree
Hide file tree
Showing 7 changed files with 327 additions and 313 deletions.
129 changes: 90 additions & 39 deletions src/CabanaPD_Force.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,11 +126,11 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i,
}

// Forward declaration.
template <class ExecutionSpace, class ForceType>
template <class MemorySpace, class ForceType>
class Force;

template <class ExecutionSpace>
class Force<ExecutionSpace, BaseForceModel>
template <class MemorySpace>
class Force<MemorySpace, BaseForceModel>
{
protected:
bool _half_neigh;
Expand All @@ -139,19 +139,79 @@ class Force<ExecutionSpace, BaseForceModel>
Timer _energy_timer;

public:
Force( const bool half_neigh )
using neighbor_list_type =
Cabana::VerletList<MemorySpace, Cabana::FullNeighborTag,
Cabana::VerletLayout2D, Cabana::TeamOpTag>;
neighbor_list_type _neigh_list;

// Primary constructor: use positions and construct neighbors.
template <class ParticleType>
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 ) )
{
}

template <class ParticleType, class NeighListType, class ParallelType>
void computeWeightedVolume( ParticleType&, const NeighListType&,
const ParallelType ) const
// Primary constructor: use positions and construct neighbors.
template <class NeighborListType>
Force( const bool half_neigh, const NeighborListType& neighbors )
: _half_neigh( half_neigh )
, _neigh_list( neighbors )
{
}

unsigned getMaxLocalNeighbors()
{
auto neigh = _neigh_list;
unsigned local_max_neighbors;
auto neigh_max = KOKKOS_LAMBDA( const int, unsigned& max_n )
{
max_n =
Cabana::NeighborList<neighbor_list_type>::maxNeighbor( neigh );
};
using exec_space = typename MemorySpace::execution_space;
Kokkos::RangePolicy<exec_space> policy( 0, 1 );
Kokkos::parallel_reduce( policy, neigh_max, local_max_neighbors );
Kokkos::fence();
return local_max_neighbors;
}
template <class ParticleType, class NeighListType, class ParallelType>
void computeDilatation( ParticleType&, const NeighListType&,
const ParallelType ) const

void getNeighborStatistics( unsigned& max_neighbors,
unsigned long long& total_neighbors )
{
auto neigh = _neigh_list;
unsigned local_max_neighbors;
unsigned long long local_total_neighbors;
auto neigh_stats = KOKKOS_LAMBDA( const int, unsigned& max_n,
unsigned long long& total_n )
{
max_n =
Cabana::NeighborList<neighbor_list_type>::maxNeighbor( neigh );
total_n = Cabana::NeighborList<neighbor_list_type>::totalNeighbor(
neigh );
};
using exec_space = typename MemorySpace::execution_space;
Kokkos::RangePolicy<exec_space> policy( 0, 1 );
Kokkos::parallel_reduce( policy, neigh_stats, local_max_neighbors,
local_total_neighbors );
Kokkos::fence();
MPI_Reduce( &local_max_neighbors, &max_neighbors, 1, MPI_UNSIGNED,
MPI_MAX, 0, MPI_COMM_WORLD );
MPI_Reduce( &local_total_neighbors, &total_neighbors, 1,
MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD );
}

// Default to no-op for pair models.
template <class ParticleType, class ParallelType>
void computeWeightedVolume( ParticleType&, const ParallelType ) const
{
}
template <class ParticleType, class ParallelType>
void computeDilatation( ParticleType&, const ParallelType ) const
{
}

Expand All @@ -162,10 +222,8 @@ class Force<ExecutionSpace, BaseForceModel>
/******************************************************************************
Force free functions.
******************************************************************************/
template <class ForceType, class ParticleType, class NeighListType,
class ParallelType>
template <class ForceType, class ParticleType, class ParallelType>
void computeForce( ForceType& force, ParticleType& particles,
const NeighListType& neigh_list,
const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
Expand All @@ -179,23 +237,19 @@ void computeForce( ForceType& force, ParticleType& particles,

// if ( half_neigh )
// Forces must be atomic for half list
// computeForce_half( f_a, x, u, neigh_list, n_local,
// computeForce_half( f_a, x, u, n_local,
// 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, neigh_list, n_local,
neigh_op_tag );
force.computeForceFull( f_a, x, u, particles, n_local, neigh_op_tag );
else
force.computeForceFull( f, x, u, particles, neigh_list, n_local,
neigh_op_tag );
force.computeForceFull( f, x, u, particles, n_local, neigh_op_tag );
Kokkos::fence();
}

template <class ForceType, class ParticleType, class NeighListType,
class ParallelType>
template <class ForceType, class ParticleType, class ParallelType>
double computeEnergy( ForceType& force, ParticleType& particles,
const NeighListType& neigh_list,
const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
Expand All @@ -210,21 +264,20 @@ double computeEnergy( ForceType& force, ParticleType& particles,

double energy;
// if ( _half_neigh )
// energy = computeEnergy_half( force, x, u, neigh_list,
// energy = computeEnergy_half( force, x, u,
// n_local, neigh_op_tag );
// else
energy = force.computeEnergyFull( W, x, u, particles, neigh_list, n_local,
neigh_op_tag );
energy =
force.computeEnergyFull( W, x, u, particles, n_local, neigh_op_tag );
Kokkos::fence();

return energy;
}

// Forces with bond breaking.
template <class ForceType, class ParticleType, class NeighListType,
class NeighborView, class ParallelType>
void computeForce( ForceType& force, ParticleType& particles,
const NeighListType& neigh_list, NeighborView& mu,
template <class ForceType, class ParticleType, class NeighborView,
class ParallelType>
void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu,
const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
Expand All @@ -238,25 +291,23 @@ void computeForce( ForceType& force, ParticleType& particles,

// if ( half_neigh )
// Forces must be atomic for half list
// computeForce_half( f_a, x, u, neigh_list, n_local,
// computeForce_half( f_a, x, u, n_local,
// 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, neigh_list, mu, n_local,
force.computeForceFull( f_a, x, u, particles, mu, n_local,
neigh_op_tag );
else
force.computeForceFull( f, x, u, particles, neigh_list, mu, n_local,
neigh_op_tag );
force.computeForceFull( f, x, u, particles, mu, n_local, neigh_op_tag );
Kokkos::fence();
}

// Energy and damage.
template <class ForceType, class ParticleType, class NeighListType,
class NeighborView, class ParallelType>
template <class ForceType, class ParticleType, class NeighborView,
class ParallelType>
double computeEnergy( ForceType& force, ParticleType& particles,
const NeighListType& neigh_list, NeighborView& mu,
const ParallelType& neigh_op_tag )
NeighborView& mu, const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
auto x = particles.sliceReferencePosition();
Expand All @@ -270,11 +321,11 @@ double computeEnergy( ForceType& force, ParticleType& particles,

double energy;
// if ( _half_neigh )
// energy = computeEnergy_half( force, x, u, neigh_list,
// energy = computeEnergy_half( force, x, u,
// n_local, neigh_op_tag );
// else
energy = force.computeEnergyFull( W, x, u, phi, particles, neigh_list, mu,
n_local, neigh_op_tag );
energy = force.computeEnergyFull( W, x, u, phi, particles, mu, n_local,
neigh_op_tag );
Kokkos::fence();

return energy;
Expand Down
46 changes: 25 additions & 21 deletions src/CabanaPD_HeatTransfer.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,34 +19,40 @@ namespace CabanaPD

// Peridynamic heat transfer with forward-Euler time integration.
// Inherits only because this is a similar neighbor-based kernel.
template <class ExecutionSpace, class ModelType>
class HeatTransfer : public Force<ExecutionSpace, BaseForceModel>
template <class MemorySpace, class ModelType>
class HeatTransfer : public Force<MemorySpace, BaseForceModel>
{
protected:
using base_type = Force<ExecutionSpace, BaseForceModel>;
using base_type = Force<MemorySpace, BaseForceModel>;
using base_type::_half_neigh;
using base_type::_timer;
using model_type = ModelType;
static_assert(
std::is_same_v<typename model_type::fracture_type, Elastic> );

Timer _euler_timer = base_type::_energy_timer;
ModelType _model;
// Using the default exec_space.
using exec_space = typename MemorySpace::execution_space;

public:
HeatTransfer( const bool half_neigh, const ModelType model )
: base_type( half_neigh )
using base_type::_neigh_list;
using model_type = ModelType;
static_assert(
std::is_same_v<typename model_type::fracture_type, Elastic> );

// Running with mechanics as well; no reason to rebuild neighbors.
template <class NeighborType>
HeatTransfer( const bool half_neigh, const NeighborType& neighbors,
const ModelType model )
: base_type( half_neigh, neighbors )
, _model( model )
{
}

template <class TemperatureType, class PosType, class ParticleType,
class NeighListType, class ParallelType>
class ParallelType>
void
computeHeatTransferFull( TemperatureType& conduction, const PosType& x,
const PosType& u, const ParticleType& particles,
const NeighListType& neigh_list, const int n_local,
ParallelType& neigh_op_tag )
const int n_local, ParallelType& neigh_op_tag )
{
_timer.start();

Expand All @@ -64,9 +70,9 @@ class HeatTransfer : public Force<ExecutionSpace, BaseForceModel>
coeff * ( temp( j ) - temp( i ) ) / xi / xi * vol( j );
};

Kokkos::RangePolicy<ExecutionSpace> policy( 0, n_local );
Kokkos::RangePolicy<exec_space> policy( 0, n_local );
Cabana::neighbor_parallel_for(
policy, temp_func, neigh_list, Cabana::FirstNeighborsTag(),
policy, temp_func, _neigh_list, Cabana::FirstNeighborsTag(),
neigh_op_tag, "CabanaPD::HeatTransfer::computeFull" );

_timer.stop();
Expand All @@ -85,19 +91,17 @@ class HeatTransfer : public Force<ExecutionSpace, BaseForceModel>
{
temp( i ) += dt / rho( i ) / model.cp * conduction( i );
};
Kokkos::RangePolicy<ExecutionSpace> policy( 0, n_local );
Kokkos::RangePolicy<exec_space> policy( 0, n_local );
Kokkos::parallel_for( "CabanaPD::HeatTransfer::forwardEuler", policy,
euler_func );
_euler_timer.stop();
}
};

// Heat transfer free function.
template <class HeatTransferType, class ParticleType, class NeighListType,
class ParallelType>
template <class HeatTransferType, class ParticleType, class ParallelType>
void computeHeatTransfer( HeatTransferType& heat_transfer,
ParticleType& particles,
const NeighListType& neigh_list,
const ParallelType& neigh_op_tag, const double dt )
{
auto n_local = particles.n_local;
Expand All @@ -111,11 +115,11 @@ 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, neigh_list, n_local, neigh_op_tag );
heat_transfer.computeHeatTransferFull( conduction_a, x, u, particles,
n_local, neigh_op_tag );
else
heat_transfer.computeHeatTransferFull(
conduction, x, u, particles, neigh_list, n_local, neigh_op_tag );
heat_transfer.computeHeatTransferFull( conduction, x, u, particles,
n_local, neigh_op_tag );
Kokkos::fence();

heat_transfer.forwardEuler( particles, dt );
Expand Down
11 changes: 6 additions & 5 deletions src/CabanaPD_Particles.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,8 +134,9 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
std::array<double, dim> local_mesh_ext;
std::array<double, dim> local_mesh_lo;
std::array<double, dim> local_mesh_hi;
std::array<double, dim> ghost_mesh_lo;
std::array<double, dim> ghost_mesh_hi;
// FIXME: this is for neighborlist construction.
double ghost_mesh_lo[dim];
double ghost_mesh_hi[dim];
std::shared_ptr<
Cabana::Grid::LocalGrid<Cabana::Grid::UniformMesh<double, dim>>>
local_grid;
Expand Down Expand Up @@ -399,9 +400,9 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
auto getForce() { return _plist_f; }
auto getReferencePosition() { return _plist_x; }

void updateCurrentPosition()
void updateCurrentPosition() const
{
_timer.start();
//_timer.start();
// Not using slice function because this is called inside.
auto y = Cabana::slice<0>( _aosoa_y, "current_positions" );
auto x = sliceReferencePosition();
Expand All @@ -414,7 +415,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
};
Kokkos::parallel_for( "CabanaPD::CalculateCurrentPositions", policy,
sum_x_u );
_timer.stop();
//_timer.stop();
}

void resize( int new_local, int new_ghost )
Expand Down
Loading

0 comments on commit 5e66739

Please sign in to comment.