Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add particle property to store current positions #32

Merged
merged 4 commits into from
Oct 4, 2023
Merged

add particle property to store current positions #32

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
c926a2e
Add optional current position output (in addition to reference)
pabloseleson Jul 21, 2023
336f2fa
fixup: allow unused if no particle output enabled
streeve Jul 24, 2023
cf0ab28
Add current position option to solver and examples
streeve Aug 29, 2023
fc61746
Calculate x+u automatically inside slice
streeve Aug 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions examples/crack_branching.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ int main( int argc, char* argv[] )
double t_final = 43e-6;
double dt = 5e-8;
double output_frequency = 5;
bool output_reference = true;

// Material constants
double E = 72e+9; // [Pa]
Expand Down Expand Up @@ -79,7 +80,7 @@ int main( int argc, char* argv[] )
CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
model_type force_model( delta, K, G0 );
CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final,
dt, output_frequency );
dt, output_frequency, output_reference );
inputs.read_args( argc, argv );

// Create particles from mesh.
Expand All @@ -91,7 +92,7 @@ int main( int argc, char* argv[] )
inputs.num_cells, halo_width );

// Define particle initialization.
auto x = particles->sliceRefPosition();
auto x = particles->sliceReferencePosition();
auto v = particles->sliceVelocity();
auto f = particles->sliceForce();
auto rho = particles->sliceDensity();
Expand Down
7 changes: 4 additions & 3 deletions examples/elastic_wave.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ int main( int argc, char* argv[] )
double t_final = 0.6;
double dt = 0.01;
int output_frequency = 5;
bool output_reference = true;
double K = 1.0;
double G = 0.5;
double delta = 0.075;
Expand All @@ -51,7 +52,7 @@ int main( int argc, char* argv[] )
model_type force_model( delta, K, G );

CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final,
dt, output_frequency );
dt, output_frequency, output_reference );
inputs.read_args( argc, argv );

// Create particles from mesh.
Expand All @@ -64,7 +65,7 @@ int main( int argc, char* argv[] )
inputs.num_cells, halo_width );

// Define particle initialization.
auto x = particles->sliceRefPosition();
auto x = particles->sliceReferencePosition();
auto u = particles->sliceDisplacement();
auto v = particles->sliceVelocity();
auto rho = particles->sliceDensity();
Expand Down Expand Up @@ -96,7 +97,7 @@ int main( int argc, char* argv[] )
cabana_pd->init_force();
cabana_pd->run();

x = particles->sliceRefPosition();
x = particles->sliceReferencePosition();
u = particles->sliceDisplacement();
double num_cell_x = inputs.num_cells[0];
auto profile = Kokkos::View<double* [2], memory_space>(
Expand Down
5 changes: 3 additions & 2 deletions examples/kalthoff_winkler.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ int main( int argc, char* argv[] )
double t_final = 70e-6;
double dt = 0.133e-6;
int output_frequency = 10;
bool output_reference = true;

// Material constants
double E = 191e+9; // [Pa]
Expand Down Expand Up @@ -82,7 +83,7 @@ int main( int argc, char* argv[] )
// CabanaPD::ForceModel<CabanaPD::LPS, CabanaPD::Fracture>;
// model_type force_model( delta, K, G, G0 );
CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final,
dt, output_frequency );
dt, output_frequency, output_reference );
inputs.read_args( argc, argv );

// Create particles from mesh.
Expand All @@ -95,7 +96,7 @@ int main( int argc, char* argv[] )
inputs.num_cells, halo_width );

// Define particle initialization.
auto x = particles->sliceRefPosition();
auto x = particles->sliceReferencePosition();
auto v = particles->sliceVelocity();
auto f = particles->sliceForce();
auto rho = particles->sliceDensity();
Expand Down
6 changes: 3 additions & 3 deletions src/CabanaPD_Boundary.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>

auto index_space = _view;
auto count = _count;
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
Kokkos::RangePolicy<ExecSpace> policy( 0, particles.n_local );
auto index_functor = KOKKOS_LAMBDA( const std::size_t pid )
{
Expand Down Expand Up @@ -163,7 +163,7 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
void apply( ExecSpace, ParticleType& particles )
{
auto f = particles.sliceForce();
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto index_space = _index_space._view;
Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
auto value = _value;
Expand Down Expand Up @@ -234,7 +234,7 @@ struct BoundaryCondition<BCIndexSpace, ForceCrackBranchBCTag>
void apply( ExecSpace, ParticleType& particles )
{
auto f = particles.sliceForce();
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto index_space = _index_space._view;
Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
auto value = _value;
Expand Down
2 changes: 1 addition & 1 deletion src/CabanaPD_Comm.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,7 @@ class Comm<ParticleType, PMB>
MPI_Comm_size( local_grid->globalGrid().comm(), &mpi_size );
MPI_Comm_rank( local_grid->globalGrid().comm(), &mpi_rank );

auto positions = particles.sliceRefPosition();
auto positions = particles.sliceReferencePosition();
// Get all 26 neighbor ranks.
auto halo_width = local_grid->haloCellWidth();
auto topology = Cajita::getTopology( *local_grid );
Expand Down
16 changes: 8 additions & 8 deletions src/CabanaPD_Force.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -402,7 +402,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
const ParallelType neigh_op_tag )
{
auto n_local = particles.n_local;
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
const auto vol = particles.sliceVolume();
auto m = particles.sliceWeightedVolume();
Expand Down Expand Up @@ -430,7 +430,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Elastic>>
const ParallelType neigh_op_tag ) const
{
auto n_local = particles.n_local;
const auto x = particles.sliceRefPosition();
const auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
const auto vol = particles.sliceVolume();
auto m = particles.sliceWeightedVolume();
Expand Down Expand Up @@ -567,7 +567,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
const MuView& mu )
{
auto n_local = particles.n_local;
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
const auto vol = particles.sliceVolume();
auto m = particles.sliceWeightedVolume();
Expand Down Expand Up @@ -606,7 +606,7 @@ class Force<ExecutionSpace, ForceModel<LPS, Fracture>>
const MuView& mu ) const
{
auto n_local = particles.n_local;
const auto x = particles.sliceRefPosition();
const auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
const auto vol = particles.sliceVolume();
auto m = particles.sliceWeightedVolume();
Expand Down Expand Up @@ -1210,7 +1210,7 @@ void computeForce( const ForceType& force, ParticleType& particles,
const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
auto f = particles.sliceForce();
auto f_a = particles.sliceForceAtomic();
Expand Down Expand Up @@ -1240,7 +1240,7 @@ double computeEnergy( const ForceType force, ParticleType& particles,
const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
auto f = particles.sliceForce();
auto W = particles.sliceStrainEnergy();
Expand Down Expand Up @@ -1269,7 +1269,7 @@ void computeForce( const ForceType& force, ParticleType& particles,
const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
auto f = particles.sliceForce();
auto f_a = particles.sliceForceAtomic();
Expand Down Expand Up @@ -1300,7 +1300,7 @@ double computeEnergy( const ForceType force, ParticleType& particles,
const ParallelType& neigh_op_tag )
{
auto n_local = particles.n_local;
auto x = particles.sliceRefPosition();
auto x = particles.sliceReferencePosition();
auto u = particles.sliceDisplacement();
auto f = particles.sliceForce();
auto W = particles.sliceStrainEnergy();
Expand Down
4 changes: 3 additions & 1 deletion src/CabanaPD_Input.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,18 +34,20 @@ class Inputs
double final_time;
double timestep;
int output_frequency;
bool output_reference;

bool half_neigh = false;

Inputs( const std::array<int, Dim> nc, std::array<double, Dim> lc,
std::array<double, Dim> hc, const double t_f, const double dt,
const int of )
const int of, const bool output_ref )
: num_cells( nc )
, low_corner( lc )
, high_corner( hc )
, final_time( t_f )
, timestep( dt )
, output_frequency( of )
, output_reference( output_ref )
{
num_steps = final_time / timestep;
}
Expand Down
78 changes: 61 additions & 17 deletions src/CabanaPD_Particles.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,7 @@ class Particles<DeviceType, PMB, Dimension>
using self_type = Particles<DeviceType, PMB, Dimension>;
using device_type = DeviceType;
using memory_space = typename device_type::memory_space;
using execution_space = typename memory_space::execution_space;
static constexpr int dim = Dimension;

// Per particle.
Expand All @@ -111,6 +112,7 @@ class Particles<DeviceType, PMB, Dimension>
// FIXME: enable variable aosoa.
using aosoa_x_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_u_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_y_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_f_type = Cabana::AoSoA<vector_type, memory_space, 1>;
using aosoa_vol_type = Cabana::AoSoA<scalar_type, memory_space, 1>;
using aosoa_nofail_type = Cabana::AoSoA<int_type, memory_space, 1>;
Expand Down Expand Up @@ -213,12 +215,13 @@ class Particles<DeviceType, PMB, Dimension>
int num_particles = particles_per_cell * owned_cells.size();
resize( num_particles, 0 );

auto x = sliceRefPosition();
auto x = sliceReferencePosition();
auto v = sliceVelocity();
auto f = sliceForce();
auto type = sliceType();
auto rho = sliceDensity();
auto u = sliceDisplacement();
auto y = sliceCurrentPosition();
auto vol = sliceVolume();
auto nofail = sliceNoFail();

Expand Down Expand Up @@ -253,6 +256,7 @@ class Particles<DeviceType, PMB, Dimension>
{
x( pid, d ) = cell_coord[d];
u( pid, d ) = 0.0;
y( pid, d ) = 0.0;
v( pid, d ) = 0.0;
f( pid, d ) = 0.0;
}
Expand Down Expand Up @@ -295,13 +299,25 @@ class Particles<DeviceType, PMB, Dimension>
KOKKOS_LAMBDA( const int pid ) { init_functor( pid ); } );
}

auto sliceRefPosition()
auto sliceReferencePosition()
{
streeve marked this conversation as resolved.
Show resolved Hide resolved
return Cabana::slice<0>( _aosoa_x, "positions" );
return Cabana::slice<0>( _aosoa_x, "reference_positions" );
}
auto sliceRefPosition() const
auto sliceReferencePosition() const
{
return Cabana::slice<0>( _aosoa_x, "positions" );
return Cabana::slice<0>( _aosoa_x, "reference_positions" );
}
auto sliceCurrentPosition()
{
// Update before returning data.
updateCurrentPosition();
return Cabana::slice<0>( _aosoa_y, "current_positions" );
}
auto sliceCurrentPosition() const
{
// Update before returning data.
updateCurrentPosition();
return Cabana::slice<0>( _aosoa_y, "current_positions" );
}
auto sliceDisplacement()
{
Expand Down Expand Up @@ -362,36 +378,62 @@ class Particles<DeviceType, PMB, Dimension>
return Cabana::slice<0>( _aosoa_nofail, "no_fail_region" );
}

void updateCurrentPosition()
{
// Not using slice function because this is called inside.
auto y = Cabana::slice<0>( _aosoa_y, "current_positions" );
auto x = sliceReferencePosition();
auto u = sliceDisplacement();
Kokkos::RangePolicy<execution_space> policy( 0, n_local + n_ghost );
auto sum_x_u = KOKKOS_LAMBDA( const std::size_t pid )
{
for ( int d = 0; d < 3; d++ )
y( pid, d ) = x( pid, d ) + u( pid, d );
};
Kokkos::parallel_for( "CabanaPD::CalculateCurrentPositions", policy,
sum_x_u );
}

void resize( int new_local, int new_ghost )
{
n_local = new_local;
n_ghost = new_ghost;

_aosoa_x.resize( new_local + new_ghost );
_aosoa_u.resize( new_local + new_ghost );
_aosoa_y.resize( new_local + new_ghost );
_aosoa_vol.resize( new_local + new_ghost );
_aosoa_f.resize( new_local );
_aosoa_other.resize( new_local );
_aosoa_nofail.resize( new_local + new_ghost );
size = _aosoa_x.size();
};

void output( const int output_step, const double output_time )
auto getPosition( const bool use_reference )
{
if ( use_reference )
return sliceReferencePosition();
else
return sliceCurrentPosition();
}

void output( [[maybe_unused]] const int output_step,
[[maybe_unused]] const double output_time,
[[maybe_unused]] const bool use_reference = true )
{
#ifdef Cabana_ENABLE_HDF5
Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
n_local, sliceRefPosition(), sliceStrainEnergy(), sliceForce(),
sliceDisplacement(), sliceVelocity(), sliceDamage() );
n_local, getPosition( use_reference ), sliceStrainEnergy(),
sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() );
#else
#ifdef Cabana_ENABLE_SILO
Cajita::Experimental::SiloParticleOutput::writePartialRangeTimeStep(
"particles", local_grid->globalGrid(), output_step, output_time, 0,
n_local, sliceRefPosition(), sliceStrainEnergy(), sliceForce(),
sliceDisplacement(), sliceVelocity(), sliceDamage() );
n_local, getPosition( use_reference ), sliceStrainEnergy(),
sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() );
#else
log( std::cout, "No particle output enabled for step ", output_step,
" (", output_time, ")" );
log( std::cout, "No particle output enabled." );
#endif
#endif
}
Expand All @@ -401,6 +443,7 @@ class Particles<DeviceType, PMB, Dimension>
protected:
aosoa_x_type _aosoa_x;
aosoa_u_type _aosoa_u;
aosoa_y_type _aosoa_y;
aosoa_f_type _aosoa_f;
aosoa_vol_type _aosoa_vol;
aosoa_nofail_type _aosoa_nofail;
Expand Down Expand Up @@ -505,12 +548,14 @@ class Particles<DeviceType, LPS, Dimension>
_aosoa_m.resize( new_local + new_ghost );
}

void output( const int output_step, const double output_time )
void output( [[maybe_unused]] const int output_step,
[[maybe_unused]] const double output_time,
[[maybe_unused]] const bool use_reference = true )
{
#ifdef Cabana_ENABLE_HDF5
Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
n_local, base_type::sliceRefPosition(),
n_local, base_type::getPosition( use_reference ),
base_type::sliceStrainEnergy(), base_type::sliceForce(),
base_type::sliceDisplacement(), base_type::sliceVelocity(),
base_type::sliceDamage(), sliceWeightedVolume(),
Expand All @@ -519,14 +564,13 @@ class Particles<DeviceType, LPS, Dimension>
#ifdef Cabana_ENABLE_SILO
Cajita::Experimental::SiloParticleOutput::writePartialRangeTimeStep(
"particles", local_grid->globalGrid(), output_step, output_time, 0,
n_local, base_type::sliceRefPosition(),
n_local, base_type::getPosition( use_reference ),
base_type::sliceStrainEnergy(), base_type::sliceForce(),
base_type::sliceDisplacement(), base_type::sliceVelocity(),
base_type::sliceDamage(), sliceWeightedVolume(),
sliceDilatation() );
#else
log( std::cout, "No particle output enabled for ", output_step, "(",
output_time, ")" );
log( std::cout, "No particle output enabled." );
#endif
#endif
}
Expand Down
Loading
Loading