From 3d359444f5eef4fc0423834fc20b815d4f178ed9 Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Thu, 27 Jul 2023 12:27:01 -0500 Subject: [PATCH 01/34] add function to read csv files --- src/CabanaPD_Particles.hpp | 94 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 6a17b447..7bde9543 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -288,6 +288,100 @@ class Particles MPI_COMM_WORLD ); } + template + void create_particles_from_file( const ExecSpace& exec_space, std::string file_name ) + { + + // read the csv file + std::vector csv_x; + std::vector csv_y; + std::vector csv_vol; + + std::vector row; + std::string line, word; + + std::fstream file (file_name, std::ios::in); + if(file.is_open()) + { + std::getline(file, line); + while(std::getline(file, line)) + { + row.clear(); + + std::stringstream str(line); + + while(std::getline(str, word, ',')) + { + row.push_back(word); + } + csv_x.push_back(std::stod(row[1])); + csv_y.push_back(std::stod(row[2])); + csv_vol.push_back(std::stod(row[3])); + + } + } + + // Create a local mesh and owned space. + auto local_mesh = Cajita::createLocalMesh( *local_grid ); + auto owned_cells = local_grid->indexSpace( + Cajita::Own(), Cajita::Cell(), Cajita::Local() ); + + int particles_per_cell = 1; + int num_particles = particles_per_cell * csv_x.size(); + resize( num_particles, 0 ); + + auto x = slice_x(); + auto v = slice_v(); + auto f = slice_f(); + auto type = slice_type(); + auto rho = slice_rho(); + auto u = slice_u(); + auto vol = slice_vol(); + auto nofail = slice_nofail(); + + auto created = Kokkos::View( + Kokkos::ViewAllocateWithoutInitializing( "particle_created" ), + num_particles ); + + // Initialize particles. + int mpi_rank = -1; + MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ); + + double cell_coord[3]; + + for (size_t i = 0; i < csv_x.size(); i++){ + + cell_coord = {csv_x[i],csv_y[i],0}; + // Set the particle position. + for ( int d = 0; d < 3; d++ ) + { + x( i, d ) = cell_coord[d]; + u( i, d ) = 0.0; + v( i, d ) = 0.0; + f( i, d ) = 0.0; + } + // FIXME: hardcoded + type( i ) = 0; + nofail( i ) = 0; + rho( i ) = 1.0; + + // Get the volume of the cell. + int empty[3]; + vol( i ) = csv_vol[i]; + + } + + n_local = csv_x.size();; + resize( n_local, 0 ); + size = _aosoa_x.size(); + + // Not using Allreduce because global count is only used for printing. + MPI_Reduce( &n_local, &n_global, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, + MPI_COMM_WORLD ); + + } + + template void update_particles( const ExecSpace, const FunctorType init_functor ) { From 5f6984ee464b532f52e982db00ac3f2db8f10a4a Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 27 Jul 2023 17:21:30 -0400 Subject: [PATCH 02/34] format --- src/CabanaPD_Particles.hpp | 105 ++++++++++++++++++------------------- 1 file changed, 52 insertions(+), 53 deletions(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 7bde9543..690dbaa3 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -289,7 +289,8 @@ class Particles } template - void create_particles_from_file( const ExecSpace& exec_space, std::string file_name ) + void create_particles_from_file( const ExecSpace& exec_space, + std::string file_name ) { // read the csv file @@ -298,28 +299,27 @@ class Particles std::vector csv_vol; std::vector row; - std::string line, word; + std::string line, word; - std::fstream file (file_name, std::ios::in); - if(file.is_open()) - { - std::getline(file, line); - while(std::getline(file, line)) - { - row.clear(); + std::fstream file( file_name, std::ios::in ); + if ( file.is_open() ) + { + std::getline( file, line ); + while ( std::getline( file, line ) ) + { + row.clear(); - std::stringstream str(line); + std::stringstream str( line ); - while(std::getline(str, word, ',')) - { - row.push_back(word); + while ( std::getline( str, word, ',' ) ) + { + row.push_back( word ); + } + csv_x.push_back( std::stod( row[1] ) ); + csv_y.push_back( std::stod( row[2] ) ); + csv_vol.push_back( std::stod( row[3] ) ); } - csv_x.push_back(std::stod(row[1])); - csv_y.push_back(std::stod(row[2])); - csv_vol.push_back(std::stod(row[3])); - - } - } + } // Create a local mesh and owned space. auto local_mesh = Cajita::createLocalMesh( *local_grid ); @@ -343,44 +343,43 @@ class Particles Kokkos::ViewAllocateWithoutInitializing( "particle_created" ), num_particles ); - // Initialize particles. - int mpi_rank = -1; - MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ); - - double cell_coord[3]; - - for (size_t i = 0; i < csv_x.size(); i++){ - - cell_coord = {csv_x[i],csv_y[i],0}; - // Set the particle position. - for ( int d = 0; d < 3; d++ ) - { - x( i, d ) = cell_coord[d]; - u( i, d ) = 0.0; - v( i, d ) = 0.0; - f( i, d ) = 0.0; - } - // FIXME: hardcoded - type( i ) = 0; - nofail( i ) = 0; - rho( i ) = 1.0; - - // Get the volume of the cell. - int empty[3]; - vol( i ) = csv_vol[i]; + // Initialize particles. + int mpi_rank = -1; + MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ); - } + double cell_coord[3]; - n_local = csv_x.size();; - resize( n_local, 0 ); - size = _aosoa_x.size(); - - // Not using Allreduce because global count is only used for printing. - MPI_Reduce( &n_local, &n_global, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, - MPI_COMM_WORLD ); + for ( size_t i = 0; i < csv_x.size(); i++ ) + { - } + cell_coord = { csv_x[i], csv_y[i], 0 }; + // Set the particle position. + for ( int d = 0; d < 3; d++ ) + { + x( i, d ) = cell_coord[d]; + u( i, d ) = 0.0; + v( i, d ) = 0.0; + f( i, d ) = 0.0; + } + // FIXME: hardcoded + type( i ) = 0; + nofail( i ) = 0; + rho( i ) = 1.0; + + // Get the volume of the cell. + int empty[3]; + vol( i ) = csv_vol[i]; + } + n_local = csv_x.size(); + ; + resize( n_local, 0 ); + size = _aosoa_x.size(); + + // Not using Allreduce because global count is only used for printing. + MPI_Reduce( &n_local, &n_global, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, + MPI_COMM_WORLD ); + } template void update_particles( const ExecSpace, const FunctorType init_functor ) From 80c1dee6b77be12fc92df54caad0d2ecc2386f6d Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 1 Aug 2023 14:34:22 -0400 Subject: [PATCH 03/34] Move csv read to example file; update after read --- examples/read_file.cpp | 199 +++++++++++++++++++++++++++++++++++++ src/CabanaPD_Particles.hpp | 127 ++++++++--------------- 2 files changed, 240 insertions(+), 86 deletions(-) create mode 100644 examples/read_file.cpp diff --git a/examples/read_file.cpp b/examples/read_file.cpp new file mode 100644 index 00000000..47e43d13 --- /dev/null +++ b/examples/read_file.cpp @@ -0,0 +1,199 @@ +#include +#include + +#include "mpi.h" + +#include + +#include + +template +void read_particles( const std::string filename, ParticleType& particles ) +{ + // Read particles from csv file. + std::vector csv_x; + std::vector csv_y; + std::vector csv_vol; + + std::vector row; + std::string line, word; + + std::fstream file( file_name, std::ios::in ); + if ( file.is_open() ) + { + std::getline( file, line ); + while ( std::getline( file, line ) ) + { + row.clear(); + + std::stringstream str( line ); + + while ( std::getline( str, word, ',' ) ) + { + row.push_back( word ); + } + csv_x.push_back( std::stod( row[1] ) ); + csv_y.push_back( std::stod( row[2] ) ); + csv_vol.push_back( std::stod( row[3] ) ); + } + } + else + throw( "Could not open file." + filename ); + + // Create unmanaged Views in order to copy to device. + Kokkos::View x_host( csv_x.data() ); + Kokkos::View y_host( csv_y.data() ); + Kokkos::View vol_host( csv_v.data() ); + + // Copy to the device. + using memory_space = typename ParticleType::memory_space; + auto x = Kokkos::create_mirror_view_and_copy( memory_space(), x_host ); + auto y = Kokkos::create_mirror_view_and_copy( memory_space(), y_host ); + auto vol = Kokkos::create_mirror_view_and_copy( memory_space(), vol_host ); + + auto px = particles.slice_x(); + auto pvol = particles.slice_vol(); + auto v = slice_v(); + auto f = slice_f(); + auto type = slice_type(); + auto rho = slice_rho(); + auto u = slice_u(); + auto nofail = slice_nofail(); + + using exec_space = typename memory_space::execution_space; + Kokkos::parallel_for( + "copy_to_particles", + Kokkos::RangePolicy policy( 0, x.size() ), + KOKKOS_LAMBDA( const int i ) { + // Set the particle position and volume. + pvol( pid ) = vol( i ); + px( pid, 0 ) = x( i ); + px( pid, 1 ) = y( i ); + + // Initialize everything else to zero. + px( pid, 2 ) = 0.0; + for ( int d = 0; d < 3; d++ ) + { + u( pid, d ) = 0.0; + v( pid, d ) = 0.0; + f( pid, d ) = 0.0; + } + type( pid ) = 0; + nofail( pid ) = 0; + rho( pid ) = 1.0; + } ); +} + +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + { + Kokkos::ScopeGuard scope_guard( argc, argv ); + + // FIXME: change backend at compile time for now. + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // Plate dimensions (m) + double height = 0.1; + double width = 0.04; + double thickness = 0.002; + + // Domain + std::array num_cell = { 300, 121, 6 }; // 300 x 120 x 6 + double low_x = -0.5 * height; + double low_y = -0.5 * width; + double low_z = -0.5 * thickness; + double high_x = 0.5 * height; + double high_y = 0.5 * width; + double high_z = 0.5 * thickness; + std::array low_corner = { low_x, low_y, low_z }; + std::array high_corner = { high_x, high_y, high_z }; + + // Time + double t_final = 43e-6; + double dt = 6.7e-8; + double output_frequency = 5; + + // Material constants + double E = 72e+9; // [Pa] + double nu = 0.25; // unitless + double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa] + double rho0 = 2440; // [kg/m^3] + double G0 = 3.8; // [J/m^2] + + // PD horizon + double delta = 0.001; + + // FIXME: set halo width based on delta + int m = std::floor( + delta / ( ( high_corner[0] - low_corner[0] ) / num_cell[0] ) ); + int halo_width = m + 1; + + // Prenotch + double L_prenotch = height / 2.0; + double y_prenotch1 = 0.0; + Kokkos::Array p01 = { low_x, y_prenotch1, low_z }; + Kokkos::Array v1 = { L_prenotch, 0, 0 }; + Kokkos::Array v2 = { 0, 0, thickness }; + Kokkos::Array, 1> notch_positions = { p01 }; + CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions ); + + // Choose force model type. + using model_type = + CabanaPD::ForceModel; + model_type force_model( delta, K, G0 ); + CabanaPD::Inputs inputs( num_cell, low_corner, high_corner, t_final, dt, + output_frequency ); + inputs.read_args( argc, argv ); + + // Default construct to then read particles. + using device_type = Kokkos::Device; + auto particles = std::make_shared>(); + + // Read particles from file. + std::string file_name = "file.csv"; + read_particles( file_name, particles ); + // Update after reading. + particles.update_after_read( exec_space(), halo_width, num_cell ); + + // Define particle initialization. + auto x = particles->slice_x(); + auto v = particles->slice_v(); + auto f = particles->slice_f(); + auto rho = particles->slice_rho(); + auto nofail = particles->slice_nofail(); + + // Relying on uniform grid here. + double dy = particles->dy; + double b0 = 2e6 / dy; // Pa + + CabanaPD::RegionBoundary plane1( low_x, high_x, low_y - dy, low_y + dy, + low_z, high_z ); + CabanaPD::RegionBoundary plane2( low_x, high_x, high_y - dy, + high_y + dy, low_z, high_z ); + std::vector planes = { plane1, plane2 }; + auto bc = + createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, + exec_space{}, *particles, planes, b0 ); + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + rho( pid ) = rho0; + // Set the no-fail zone. + if ( x( pid, 1 ) <= plane1.low_y + delta + 1e-10 || + x( pid, 1 ) >= plane2.high_y - delta - 1e-10 ) + nofail( pid ) = 1; + }; + particles->update_particles( exec_space{}, init_functor ); + + // FIXME: use createSolver to switch backend at runtime. + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model, bc, prenotch ); + cabana_pd->init_force(); + cabana_pd->run(); + } + + MPI_Finalize(); +} diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 690dbaa3..d2fb849b 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -283,99 +283,54 @@ class Particles resize( n_local, 0 ); size = _aosoa_x.size(); - // Not using Allreduce because global count is only used for printing. - MPI_Reduce( &n_local, &n_global, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, - MPI_COMM_WORLD ); + update_global(); } + // This is necessary after reading in particles from file for a consistent + // state. template - void create_particles_from_file( const ExecSpace& exec_space, - std::string file_name ) + void update_after_read( const ExecSpace& exec_space, const int hw, + const std::array num_cells ) { - - // read the csv file - std::vector csv_x; - std::vector csv_y; - std::vector csv_vol; - - std::vector row; - std::string line, word; - - std::fstream file( file_name, std::ios::in ); - if ( file.is_open() ) - { - std::getline( file, line ); - while ( std::getline( file, line ) ) - { - row.clear(); - - std::stringstream str( line ); - - while ( std::getline( str, word, ',' ) ) - { - row.push_back( word ); - } - csv_x.push_back( std::stod( row[1] ) ); - csv_y.push_back( std::stod( row[2] ) ); - csv_vol.push_back( std::stod( row[3] ) ); - } - } - - // Create a local mesh and owned space. - auto local_mesh = Cajita::createLocalMesh( *local_grid ); - auto owned_cells = local_grid->indexSpace( - Cajita::Own(), Cajita::Cell(), Cajita::Local() ); - - int particles_per_cell = 1; - int num_particles = particles_per_cell * csv_x.size(); - resize( num_particles, 0 ); - + halo_width = hw; auto x = slice_x(); - auto v = slice_v(); - auto f = slice_f(); - auto type = slice_type(); - auto rho = slice_rho(); - auto u = slice_u(); - auto vol = slice_vol(); - auto nofail = slice_nofail(); - - auto created = Kokkos::View( - Kokkos::ViewAllocateWithoutInitializing( "particle_created" ), - num_particles ); - - // Initialize particles. - int mpi_rank = -1; - MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank ); - - double cell_coord[3]; - - for ( size_t i = 0; i < csv_x.size(); i++ ) - { - - cell_coord = { csv_x[i], csv_y[i], 0 }; - // Set the particle position. - for ( int d = 0; d < 3; d++ ) - { - x( i, d ) = cell_coord[d]; - u( i, d ) = 0.0; - v( i, d ) = 0.0; - f( i, d ) = 0.0; - } - // FIXME: hardcoded - type( i ) = 0; - nofail( i ) = 0; - rho( i ) = 1.0; - - // Get the volume of the cell. - int empty[3]; - vol( i ) = csv_vol[i]; - } + n_local = x.size(); + n_ghost = 0; + size = n_local; + update_global(); + + double min_x = 0.0; + double min_y = 0.0; + double min_z = 0.0; + double max_x = 0.0; + double max_y = 0.0; + double max_z = 0.0; + Kokkos::parallel_reduce( + "CabanaPD::Particles::min_max_positions", + Kokkos::RangePolicy( exec_space, 0, n_local ), + KOKKOS_LAMBDA( const int i, int& min_x, int& min_y, int& min_z, + int& max_x, int& max_y, int& max_z ) { + if ( x( i, 0 ) > max_x ) + max_x = x( i, 0 ); + else if ( x( i, 0 ) < min_x ) + min_x = x( i, 0 ); + if ( x( i, 1 ) > max_y ) + max_y = x( i, 1 ); + else if ( x( i, 1 ) < min_y ) + min_y = x( i, 1 ); + if ( x( i, 2 ) > max_z ) + max_z = x( i, 2 ); + else if ( x( i, 2 ) < min_z ) + min_z = x( i, 2 ); + }, + min_x, min_y, min_z, max_x, max_y, max_z ); - n_local = csv_x.size(); - ; - resize( n_local, 0 ); - size = _aosoa_x.size(); + create_domain( { min_x, min_y, min_z }, { max_x, max_y, max_z }, + num_cells ); + } + void update_global() + { // Not using Allreduce because global count is only used for printing. MPI_Reduce( &n_local, &n_global, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); From 2137af1542e3e1eaecabbb7e2d80d57aad874b61 Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Wed, 30 Aug 2023 23:07:32 -0500 Subject: [PATCH 04/34] Add read file example to CMake --- examples/CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index df1b8e15..79a697a6 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -7,4 +7,7 @@ target_link_libraries(KalthoffWinkler LINK_PUBLIC CabanaPD) add_executable(CrackBranching crack_branching.cpp) target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) +add_executable(ReadFile read_file.cpp) +target_link_libraries(ReadFile LINK_PUBLIC CabanaPD) + install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR}) From 181608ac43c81e766af7697e6f7ff01d631e4021 Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Wed, 30 Aug 2023 23:10:34 -0500 Subject: [PATCH 05/34] Fix some C++ compiler errors related to tread the CSV file --- examples/read_file.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 47e43d13..353299a4 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -13,12 +13,12 @@ void read_particles( const std::string filename, ParticleType& particles ) // Read particles from csv file. std::vector csv_x; std::vector csv_y; - std::vector csv_vol; + std::vector csv_v; std::vector row; std::string line, word; - std::fstream file( file_name, std::ios::in ); + std::fstream file( filename, std::ios::in ); if ( file.is_open() ) { std::getline( file, line ); @@ -34,7 +34,7 @@ void read_particles( const std::string filename, ParticleType& particles ) } csv_x.push_back( std::stod( row[1] ) ); csv_y.push_back( std::stod( row[2] ) ); - csv_vol.push_back( std::stod( row[3] ) ); + csv_v.push_back( std::stod( row[3] ) ); } } else From 4eadd5129d6f45d6716d96062379f9f5c3be12d8 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 31 Aug 2023 11:15:16 -0400 Subject: [PATCH 06/34] fixup: type mistake --- src/CabanaPD_Particles.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 254a2789..dcce2565 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -304,8 +304,9 @@ class Particles Kokkos::parallel_reduce( "CabanaPD::Particles::min_max_positions", Kokkos::RangePolicy( exec_space, 0, n_local ), - KOKKOS_LAMBDA( const int i, int& min_x, int& min_y, int& min_z, - int& max_x, int& max_y, int& max_z ) { + KOKKOS_LAMBDA( const int i, double& min_x, double& min_y, + double& min_z, double& max_x, double& max_y, + double& max_z ) { if ( x( i, 0 ) > max_x ) max_x = x( i, 0 ); else if ( x( i, 0 ) < min_x ) From d464f37e37a30871feaed35a8b1e5d46b07c3f60 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 31 Aug 2023 11:15:25 -0400 Subject: [PATCH 07/34] fixup: file read example compilation --- examples/read_file.cpp | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 353299a4..fc6e1af1 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -53,22 +53,21 @@ void read_particles( const std::string filename, ParticleType& particles ) auto px = particles.slice_x(); auto pvol = particles.slice_vol(); - auto v = slice_v(); - auto f = slice_f(); - auto type = slice_type(); - auto rho = slice_rho(); - auto u = slice_u(); - auto nofail = slice_nofail(); + auto v = particles.slice_v(); + auto f = particles.slice_f(); + auto type = particles.slice_type(); + auto rho = particles.slice_rho(); + auto u = particles.slice_u(); + auto nofail = particles.slice_nofail(); using exec_space = typename memory_space::execution_space; Kokkos::parallel_for( - "copy_to_particles", - Kokkos::RangePolicy policy( 0, x.size() ), - KOKKOS_LAMBDA( const int i ) { + "copy_to_particles", Kokkos::RangePolicy( 0, x.size() ), + KOKKOS_LAMBDA( const int pid ) { // Set the particle position and volume. - pvol( pid ) = vol( i ); - px( pid, 0 ) = x( i ); - px( pid, 1 ) = y( i ); + pvol( pid ) = vol( pid ); + px( pid, 0 ) = x( pid ); + px( pid, 1 ) = y( pid ); // Initialize everything else to zero. px( pid, 2 ) = 0.0; @@ -154,9 +153,9 @@ int main( int argc, char* argv[] ) // Read particles from file. std::string file_name = "file.csv"; - read_particles( file_name, particles ); + read_particles( file_name, *particles ); // Update after reading. - particles.update_after_read( exec_space(), halo_width, num_cell ); + particles->update_after_read( exec_space(), halo_width, num_cell ); // Define particle initialization. auto x = particles->slice_x(); From 2d04d49f3d8c9063b2c716bf28768570461d091c Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Sat, 23 Sep 2023 16:18:30 -0400 Subject: [PATCH 08/34] fixup: function rename --- examples/read_file.cpp | 28 ++++++++++++++-------------- src/CabanaPD_Particles.hpp | 6 +++--- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index fc6e1af1..5a293cbf 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -51,14 +51,14 @@ void read_particles( const std::string filename, ParticleType& particles ) auto y = Kokkos::create_mirror_view_and_copy( memory_space(), y_host ); auto vol = Kokkos::create_mirror_view_and_copy( memory_space(), vol_host ); - auto px = particles.slice_x(); - auto pvol = particles.slice_vol(); - auto v = particles.slice_v(); - auto f = particles.slice_f(); - auto type = particles.slice_type(); - auto rho = particles.slice_rho(); - auto u = particles.slice_u(); - auto nofail = particles.slice_nofail(); + auto px = particles.sliceRefPosition(); + auto pvol = particles.sliceVolume(); + auto v = particles.sliceVelocity(); + auto f = particles.sliceForce(); + auto type = particles.sliceType(); + auto rho = particles.sliceDensity(); + auto u = particles.sliceDisplacement(); + auto nofail = particles.sliceNoFail(); using exec_space = typename memory_space::execution_space; Kokkos::parallel_for( @@ -158,11 +158,11 @@ int main( int argc, char* argv[] ) particles->update_after_read( exec_space(), halo_width, num_cell ); // Define particle initialization. - auto x = particles->slice_x(); - auto v = particles->slice_v(); - auto f = particles->slice_f(); - auto rho = particles->slice_rho(); - auto nofail = particles->slice_nofail(); + auto x = particles->sliceRefPosition(); + auto v = particles->sliceVelocity(); + auto f = particles->sliceForce(); + auto rho = particles->sliceDensity(); + auto nofail = particles->sliceNoFail(); // Relying on uniform grid here. double dy = particles->dy; @@ -185,7 +185,7 @@ int main( int argc, char* argv[] ) x( pid, 1 ) >= plane2.high_y - delta - 1e-10 ) nofail( pid ) = 1; }; - particles->update_particles( exec_space{}, init_functor ); + particles->updateParticles( exec_space{}, init_functor ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 3d99b0c7..827bde88 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -288,7 +288,7 @@ class Particles const std::array num_cells ) { halo_width = hw; - auto x = slice_x(); + auto x = sliceRefPosition(); n_local = x.size(); n_ghost = 0; size = n_local; @@ -321,8 +321,8 @@ class Particles }, min_x, min_y, min_z, max_x, max_y, max_z ); - create_domain( { min_x, min_y, min_z }, { max_x, max_y, max_z }, - num_cells ); + createDomain( { min_x, min_y, min_z }, { max_x, max_y, max_z }, + num_cells ); } void update_global() From ade2d1bd6d6584a957ff6afb3098b41dd1d6d562 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 28 Sep 2023 17:33:17 -0400 Subject: [PATCH 09/34] fixup: file read and example details --- examples/read_file.cpp | 93 +++++++++++--------------------------- src/CabanaPD_Input.cpp | 6 +++ src/CabanaPD_Input.hpp | 1 + src/CabanaPD_Particles.hpp | 41 ++++++++++++----- 4 files changed, 64 insertions(+), 77 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 5a293cbf..34ebdf20 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -32,24 +32,29 @@ void read_particles( const std::string filename, ParticleType& particles ) { row.push_back( word ); } - csv_x.push_back( std::stod( row[1] ) ); - csv_y.push_back( std::stod( row[2] ) ); - csv_v.push_back( std::stod( row[3] ) ); + csv_x.push_back( std::stod( row[0] ) ); + csv_y.push_back( std::stod( row[1] ) ); + csv_v.push_back( std::stod( row[2] ) ); } } else throw( "Could not open file." + filename ); // Create unmanaged Views in order to copy to device. - Kokkos::View x_host( csv_x.data() ); - Kokkos::View y_host( csv_y.data() ); - Kokkos::View vol_host( csv_v.data() ); + Kokkos::View x_host( csv_x.data(), + csv_x.size() ); + Kokkos::View y_host( csv_y.data(), + csv_y.size() ); + Kokkos::View vol_host( csv_v.data(), + csv_v.size() ); // Copy to the device. using memory_space = typename ParticleType::memory_space; auto x = Kokkos::create_mirror_view_and_copy( memory_space(), x_host ); auto y = Kokkos::create_mirror_view_and_copy( memory_space(), y_host ); auto vol = Kokkos::create_mirror_view_and_copy( memory_space(), vol_host ); + // Resize internal variables (no ghosts initially). + particles.resize( x.size(), 0 ); auto px = particles.sliceRefPosition(); auto pvol = particles.sliceVolume(); @@ -68,6 +73,7 @@ void read_particles( const std::string filename, ParticleType& particles ) pvol( pid ) = vol( pid ); px( pid, 0 ) = x( pid ); px( pid, 1 ) = y( pid ); + px( pid, 2 ) = 0.0; // Initialize everything else to zero. px( pid, 2 ) = 0.0; @@ -93,56 +99,28 @@ int main( int argc, char* argv[] ) using exec_space = Kokkos::DefaultExecutionSpace; using memory_space = typename exec_space::memory_space; - // Plate dimensions (m) - double height = 0.1; - double width = 0.04; - double thickness = 0.002; - - // Domain - std::array num_cell = { 300, 121, 6 }; // 300 x 120 x 6 - double low_x = -0.5 * height; - double low_y = -0.5 * width; - double low_z = -0.5 * thickness; - double high_x = 0.5 * height; - double high_y = 0.5 * width; - double high_z = 0.5 * thickness; - std::array low_corner = { low_x, low_y, low_z }; - std::array high_corner = { high_x, high_y, high_z }; - // Time - double t_final = 43e-6; - double dt = 6.7e-8; - double output_frequency = 5; + double t_final = 1e-6; + double dt = 1e-8; + double output_frequency = 10; // Material constants double E = 72e+9; // [Pa] double nu = 0.25; // unitless double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa] - double rho0 = 2440; // [kg/m^3] double G0 = 3.8; // [J/m^2] // PD horizon - double delta = 0.001; - - // FIXME: set halo width based on delta - int m = std::floor( - delta / ( ( high_corner[0] - low_corner[0] ) / num_cell[0] ) ); - int halo_width = m + 1; - - // Prenotch - double L_prenotch = height / 2.0; - double y_prenotch1 = 0.0; - Kokkos::Array p01 = { low_x, y_prenotch1, low_z }; - Kokkos::Array v1 = { L_prenotch, 0, 0 }; - Kokkos::Array v2 = { 0, 0, thickness }; - Kokkos::Array, 1> notch_positions = { p01 }; - CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions ); + double delta = 0.1; + int halo_width = 1; // Choose force model type. using model_type = CabanaPD::ForceModel; model_type force_model( delta, K, G0 ); - CabanaPD::Inputs inputs( num_cell, low_corner, high_corner, t_final, dt, + std::array zero = { 0.0, 0.0, 0.0 }; + std::array num_cell = { 2, 20, 1 }; + CabanaPD::Inputs inputs( num_cell, zero, zero, t_final, dt, output_frequency ); inputs.read_args( argc, argv ); @@ -152,10 +130,10 @@ int main( int argc, char* argv[] ) device_type, typename model_type::base_model>>(); // Read particles from file. - std::string file_name = "file.csv"; - read_particles( file_name, *particles ); - // Update after reading. - particles->update_after_read( exec_space(), halo_width, num_cell ); + read_particles( inputs.input_file, *particles ); + // Update after reading. Currently requires fixed cell spacing. + double dx = 0.5; + particles->updateAfterRead( exec_space(), halo_width, dx ); // Define particle initialization. auto x = particles->sliceRefPosition(); @@ -164,28 +142,11 @@ int main( int argc, char* argv[] ) auto rho = particles->sliceDensity(); auto nofail = particles->sliceNoFail(); - // Relying on uniform grid here. - double dy = particles->dy; - double b0 = 2e6 / dy; // Pa - - CabanaPD::RegionBoundary plane1( low_x, high_x, low_y - dy, low_y + dy, - low_z, high_z ); - CabanaPD::RegionBoundary plane2( low_x, high_x, high_y - dy, - high_y + dy, low_z, high_z ); - std::vector planes = { plane1, plane2 }; + CabanaPD::Prenotch<0> prenotch; + std::vector planes = {}; auto bc = createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, - exec_space{}, *particles, planes, b0 ); - - auto init_functor = KOKKOS_LAMBDA( const int pid ) - { - rho( pid ) = rho0; - // Set the no-fail zone. - if ( x( pid, 1 ) <= plane1.low_y + delta + 1e-10 || - x( pid, 1 ) >= plane2.high_y - delta - 1e-10 ) - nofail( pid ) = 1; - }; - particles->updateParticles( exec_space{}, init_functor ); + exec_space{}, *particles, planes, 10.0 ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( diff --git a/src/CabanaPD_Input.cpp b/src/CabanaPD_Input.cpp index 4c9f4c22..6de7dab2 100644 --- a/src/CabanaPD_Input.cpp +++ b/src/CabanaPD_Input.cpp @@ -116,6 +116,12 @@ void Inputs::read_args( int argc, char* argv[] ) output_file = argv[i + 1]; ++i; } + else if ( ( strcmp( argv[i], "-i" ) == 0 ) || + ( strcmp( argv[i], "--input-file" ) == 0 ) ) + { + input_file = argv[i + 1]; + ++i; + } else if ( ( strcmp( argv[i], "-e" ) == 0 ) || ( strcmp( argv[i], "--error-file" ) == 0 ) ) { diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp index 1eea581a..5d99e3e1 100644 --- a/src/CabanaPD_Input.hpp +++ b/src/CabanaPD_Input.hpp @@ -21,6 +21,7 @@ class Inputs public: std::string output_file = "cabanaPD.out"; std::string error_file = "cabanaPD.err"; + std::string input_file = "particles.csv"; std::string device_type = "SERIAL"; std::array num_cells; diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 827bde88..a2fa8e02 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -284,8 +284,7 @@ class Particles // This is necessary after reading in particles from file for a consistent // state. template - void update_after_read( const ExecSpace& exec_space, const int hw, - const std::array num_cells ) + void updateAfterRead( const ExecSpace& exec_space, const int hw, double dx ) { halo_width = hw; auto x = sliceRefPosition(); @@ -294,12 +293,18 @@ class Particles size = n_local; update_global(); - double min_x = 0.0; - double min_y = 0.0; - double min_z = 0.0; - double max_x = 0.0; - double max_y = 0.0; - double max_z = 0.0; + double max_x; + Kokkos::Max max_x_reducer( max_x ); + double min_x; + Kokkos::Min min_x_reducer( min_x ); + double max_y; + Kokkos::Max max_y_reducer( max_y ); + double min_y; + Kokkos::Min min_y_reducer( min_y ); + double max_z; + Kokkos::Max max_z_reducer( max_z ); + double min_z; + Kokkos::Min min_z_reducer( min_z ); Kokkos::parallel_reduce( "CabanaPD::Particles::min_max_positions", Kokkos::RangePolicy( exec_space, 0, n_local ), @@ -319,10 +324,24 @@ class Particles else if ( x( i, 2 ) < min_z ) min_z = x( i, 2 ); }, - min_x, min_y, min_z, max_x, max_y, max_z ); + min_x_reducer, min_y_reducer, min_z_reducer, max_x_reducer, + max_y_reducer, max_z_reducer ); - createDomain( { min_x, min_y, min_z }, { max_x, max_y, max_z }, - num_cells ); + std::array min_corner = { min_x, min_y, min_z }; + std::array max_corner = { max_x, max_y, max_z }; + + std::array num_cells; + for ( int d = 0; d < 3; d++ ) + { + // if ( max_corner[d] - min_corner[d] < dx ) + { + max_corner[d] += dx / 2.0; + min_corner[d] -= dx / 2.0; + } + num_cells[d] = + static_cast( ( max_corner[d] - min_corner[d] ) / dx ); + } + createDomain( min_corner, max_corner, num_cells ); } void update_global() From f7f2ce18324a4f5fb0c0e7e4e84039b97eaa7bc4 Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Thu, 5 Oct 2023 10:03:59 -0500 Subject: [PATCH 10/34] Add load in force --- examples/read_file.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 34ebdf20..06e4772b 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -111,7 +111,7 @@ int main( int argc, char* argv[] ) double G0 = 3.8; // [J/m^2] // PD horizon - double delta = 0.1; + double delta = 1.5; int halo_width = 1; // Choose force model type. @@ -143,10 +143,14 @@ int main( int argc, char* argv[] ) auto nofail = particles->sliceNoFail(); CabanaPD::Prenotch<0> prenotch; - std::vector planes = {}; + + CabanaPD::RegionBoundary plane1( 0 , 1, 9, 10, + -1, 1 ); + + std::vector planes = { plane1}; auto bc = createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, - exec_space{}, *particles, planes, 10.0 ); + exec_space{}, *particles, planes, 1.0 ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( From 02b75f0f29251016aa3fbd6087ce5e8d23d21339 Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Thu, 5 Oct 2023 10:13:58 -0500 Subject: [PATCH 11/34] Clang formart --- examples/read_file.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 06e4772b..bb10643a 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -144,10 +144,9 @@ int main( int argc, char* argv[] ) CabanaPD::Prenotch<0> prenotch; - CabanaPD::RegionBoundary plane1( 0 , 1, 9, 10, - -1, 1 ); + CabanaPD::RegionBoundary plane1( 0, 1, 9, 10, -1, 1 ); - std::vector planes = { plane1}; + std::vector planes = { plane1 }; auto bc = createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, exec_space{}, *particles, planes, 1.0 ); From 456f13550296de876ee46e630ac92105412466e9 Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Thu, 5 Oct 2023 20:26:19 -0500 Subject: [PATCH 12/34] Update force --- examples/read_file.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index bb10643a..d53c964e 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -101,7 +101,7 @@ int main( int argc, char* argv[] ) // Time double t_final = 1e-6; - double dt = 1e-8; + double dt = 1e-9; double output_frequency = 10; // Material constants @@ -111,7 +111,7 @@ int main( int argc, char* argv[] ) double G0 = 3.8; // [J/m^2] // PD horizon - double delta = 1.5; + double delta = 3*0.1299999999999999; int halo_width = 1; // Choose force model type. @@ -119,7 +119,7 @@ int main( int argc, char* argv[] ) CabanaPD::ForceModel; model_type force_model( delta, K, G0 ); std::array zero = { 0.0, 0.0, 0.0 }; - std::array num_cell = { 2, 20, 1 }; + std::array num_cell = { 8, 80, 1 }; CabanaPD::Inputs inputs( num_cell, zero, zero, t_final, dt, output_frequency ); inputs.read_args( argc, argv ); @@ -132,7 +132,7 @@ int main( int argc, char* argv[] ) // Read particles from file. read_particles( inputs.input_file, *particles ); // Update after reading. Currently requires fixed cell spacing. - double dx = 0.5; + double dx =0.1299999999999999; particles->updateAfterRead( exec_space(), halo_width, dx ); // Define particle initialization. @@ -145,11 +145,12 @@ int main( int argc, char* argv[] ) CabanaPD::Prenotch<0> prenotch; CabanaPD::RegionBoundary plane1( 0, 1, 9, 10, -1, 1 ); + std::vector planes = { plane1 }; auto bc = - createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, - exec_space{}, *particles, planes, 1.0 ); + createBoundaryCondition( CabanaPD::ForceUpdateBCTag{}, + exec_space{}, *particles, planes, 2e6 / dx / dx ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( From ce2e627eb3da5d98c089aed1f9c42b1735fbd703 Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Thu, 5 Oct 2023 20:26:38 -0500 Subject: [PATCH 13/34] Clang format --- examples/read_file.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index d53c964e..ad3bc80e 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -111,7 +111,7 @@ int main( int argc, char* argv[] ) double G0 = 3.8; // [J/m^2] // PD horizon - double delta = 3*0.1299999999999999; + double delta = 3 * 0.1299999999999999; int halo_width = 1; // Choose force model type. @@ -132,7 +132,7 @@ int main( int argc, char* argv[] ) // Read particles from file. read_particles( inputs.input_file, *particles ); // Update after reading. Currently requires fixed cell spacing. - double dx =0.1299999999999999; + double dx = 0.1299999999999999; particles->updateAfterRead( exec_space(), halo_width, dx ); // Define particle initialization. @@ -145,12 +145,11 @@ int main( int argc, char* argv[] ) CabanaPD::Prenotch<0> prenotch; CabanaPD::RegionBoundary plane1( 0, 1, 9, 10, -1, 1 ); - std::vector planes = { plane1 }; auto bc = - createBoundaryCondition( CabanaPD::ForceUpdateBCTag{}, - exec_space{}, *particles, planes, 2e6 / dx / dx ); + createBoundaryCondition( CabanaPD::ForceUpdateBCTag{}, exec_space{}, + *particles, planes, 2e6 / dx / dx ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( From 01b342f8b6e6fbf1da0495510f22fc560c955c9d Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 6 Oct 2023 13:11:07 -0400 Subject: [PATCH 14/34] Generalize symmetric 1d boundary condition --- examples/crack_branching.cpp | 2 +- src/CabanaPD_Boundary.hpp | 43 ++++++++++++++++++++++++++++++------ 2 files changed, 37 insertions(+), 8 deletions(-) diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp index e0c80731..1e7d99cd 100644 --- a/examples/crack_branching.cpp +++ b/examples/crack_branching.cpp @@ -107,7 +107,7 @@ int main( int argc, char* argv[] ) high_y + dy, low_z, high_z ); std::vector planes = { plane1, plane2 }; auto bc = - createBoundaryCondition( CabanaPD::ForceCrackBranchBCTag{}, + createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes, b0 ); auto init_functor = KOKKOS_LAMBDA( const int pid ) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index d0cb8e31..08fd3e98 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -118,7 +118,7 @@ struct BoundaryIndexSpace template auto createBoundaryIndexSpace( ExecSpace exec_space, Particles particles, - std::vector planes, + std::vector planes, const double initial_guess ) { using memory_space = typename Particles::memory_space; @@ -133,13 +133,14 @@ struct ForceUpdateBCTag { }; -struct ForceCrackBranchBCTag +struct ForceSymmetric1dBCTag { }; template struct BoundaryCondition; +// Reset BC (all dimensions). template struct BoundaryCondition { @@ -176,6 +177,7 @@ struct BoundaryCondition } }; +// Increment BC (all dimensions). template struct BoundaryCondition { @@ -211,15 +213,22 @@ struct BoundaryCondition } }; +// Symmetric 1d BC applied with opposite sign based on position from the +// midplane in the specified direction. template -struct BoundaryCondition +struct BoundaryCondition { double _value; BCIndexSpace _index_space; + int _dim; + double _center; - BoundaryCondition( const double value, BCIndexSpace bc_index_space ) + BoundaryCondition( const double value, BCIndexSpace bc_index_space, + const int dim = 1, const double center = 0.0 ) : _value( value ) , _index_space( bc_index_space ) + , _dim( dim ) + , _center( center ) { } @@ -238,18 +247,21 @@ struct BoundaryCondition auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; + auto dim = _dim; + auto center = _center; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { auto pid = index_space( b ); - // This is specifically for the crack branching. - auto sign = std::abs( x( pid, 1 ) ) / x( pid, 1 ); + auto relative_x = x( pid, dim ) - center; + auto sign = std::abs( relative_x ) / relative_x; f( pid, 1 ) += value * sign; } ); } }; // FIXME: relatively large initial guess for allocation. -template +template auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles, std::vector planes, const double value, @@ -262,6 +274,23 @@ auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles, return BoundaryCondition( value, bc_indices ); } +// FIXME: relatively large initial guess for allocation. +template +auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space, + Particles particles, + std::vector planes, + const double value, const int dim, + const double center, + const double initial_guess = 0.5 ) +{ + using memory_space = typename Particles::memory_space; + using bc_index_type = BoundaryIndexSpace; + bc_index_type bc_indices = createBoundaryIndexSpace( + exec_space, particles, planes, initial_guess ); + return BoundaryCondition( + value, bc_indices, dim, center ); +} + } // namespace CabanaPD #endif From 5b5f71898b755e6cbc5e1fa724edf2b7d4e7eecc Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 6 Oct 2023 13:11:34 -0400 Subject: [PATCH 15/34] Use new generalized symmetric BC --- examples/read_file.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index ad3bc80e..3bc36a16 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -147,9 +147,12 @@ int main( int argc, char* argv[] ) CabanaPD::RegionBoundary plane1( 0, 1, 9, 10, -1, 1 ); std::vector planes = { plane1 }; - auto bc = - createBoundaryCondition( CabanaPD::ForceUpdateBCTag{}, exec_space{}, - *particles, planes, 2e6 / dx / dx ); + int bc_dim = 2; + double center = particles->local_mesh_ext[bc_dim] / 2.0 + + particles->local_mesh_lo[bc_dim]; + auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, + exec_space{}, *particles, planes, + 2e6 / dx / dx, bc_dim, center ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( From d08830d06a691da405ef4e229eebef7bf2d89c2e Mon Sep 17 00:00:00 2001 From: Patrick Diehl Date: Fri, 6 Oct 2023 15:01:53 -0500 Subject: [PATCH 16/34] Center the geometry and pull on top and bottom --- examples/read_file.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 189cabb7..e329a0ec 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -139,9 +139,11 @@ int main( int argc, char* argv[] ) CabanaPD::Prenotch<0> prenotch; - CabanaPD::RegionBoundary plane1( 0, 1, 9, 10, -1, 1 ); + CabanaPD::RegionBoundary planeUpper( -0.5, 0.5, 4.5, 5, -1, 1 ); + CabanaPD::RegionBoundary planeLower( -0.5, 0.5, -5, -4.5, -1, 1 ); - std::vector planes = { plane1 }; + std::vector planes = { planeUpper, + planeLower }; int bc_dim = 2; double center = particles->local_mesh_ext[bc_dim] / 2.0 + particles->local_mesh_lo[bc_dim]; From 93945c630abf450ba194f0d100288edd5757d559 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 12 Oct 2023 09:44:20 -0700 Subject: [PATCH 17/34] fixup: run as 3d for now --- examples/read_file.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index e329a0ec..6ba411ea 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -116,13 +116,13 @@ int main( int argc, char* argv[] ) using model_type = CabanaPD::ForceModel; model_type force_model( delta, K, G0 ); - CabanaPD::Inputs<2> inputs( t_final, dt, output_frequency, true ); + CabanaPD::Inputs<3> inputs( t_final, dt, output_frequency, true ); inputs.read_args( argc, argv ); // Default construct to then read particles. using device_type = Kokkos::Device; auto particles = std::make_shared>(); + device_type, typename model_type::base_model, 3>>(); // Read particles from file. read_particles( inputs.input_file, *particles ); From 6cc42defe5276eecc1046da27e5cf89e2624c103 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 12 Oct 2023 09:45:43 -0700 Subject: [PATCH 18/34] fixup: clean up file example --- examples/read_file.cpp | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 6ba411ea..ad3624e6 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -109,7 +109,8 @@ int main( int argc, char* argv[] ) double G0 = 3.8; // [J/m^2] // PD horizon - double delta = 3 * 0.1299999999999999; + double dx = 0.1299999999999999; + double delta = 3 * dx; int halo_width = 1; // Choose force model type. @@ -127,16 +128,8 @@ int main( int argc, char* argv[] ) // Read particles from file. read_particles( inputs.input_file, *particles ); // Update after reading. Currently requires fixed cell spacing. - double dx = 0.1299999999999999; particles->updateAfterRead( exec_space(), halo_width, dx ); - // Define particle initialization. - auto x = particles->sliceReferencePosition(); - auto v = particles->sliceVelocity(); - auto f = particles->sliceForce(); - auto rho = particles->sliceDensity(); - auto nofail = particles->sliceNoFail(); - CabanaPD::Prenotch<0> prenotch; CabanaPD::RegionBoundary planeUpper( -0.5, 0.5, 4.5, 5, -1, 1 ); From 8d021e63f9f3fc8398f0a8fdaaea49daae6f328d Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 16 Oct 2023 08:51:33 -0400 Subject: [PATCH 19/34] fixup: init particle fields for file read --- examples/read_file.cpp | 6 +++++- src/CabanaPD_Particles.hpp | 10 ++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index ad3624e6..1b3f0550 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -64,6 +64,7 @@ void read_particles( const std::string filename, ParticleType& particles ) auto rho = particles.sliceDensity(); auto u = particles.sliceDisplacement(); auto nofail = particles.sliceNoFail(); + auto damage = particles.sliceDamage(); using exec_space = typename memory_space::execution_space; Kokkos::parallel_for( @@ -75,7 +76,9 @@ void read_particles( const std::string filename, ParticleType& particles ) px( pid, 1 ) = y( pid ); // Initialize everything else to zero. - for ( int d = 0; d < 2; d++ ) + // Currently psuedo-2d. + px( pid, 2 ) = 0.0; + for ( int d = 0; d < 3; d++ ) { u( pid, d ) = 0.0; v( pid, d ) = 0.0; @@ -84,6 +87,7 @@ void read_particles( const std::string filename, ParticleType& particles ) type( pid ) = 0; nofail( pid ) = 0; rho( pid ) = 1.0; + damage( pid ) = 0; } ); } diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 0dc60304..9119ecdf 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -639,6 +639,16 @@ class Particles _aosoa_theta.resize( 0 ); } + template + void updateAfterRead( const ExecSpace& exec_space, const int hw, double dx ) + { + base_type::updateAfterRead( exec_space, hw, dx ); + + // Only need to resize LPS variables + _aosoa_theta.resize( n_local + n_ghost ); + _aosoa_m.resize( n_local + n_ghost ); + } + auto sliceDilatation() { return Cabana::slice<0>( _aosoa_theta, "dilatation" ); From 5671fac8f3638e51cdb2af9fddbb3e3bd4697711 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 16 Oct 2023 08:51:57 -0400 Subject: [PATCH 20/34] fixup: assert BC plane inputs are reasonable --- src/CabanaPD_Boundary.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 7c1361a6..889700c5 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -42,7 +42,12 @@ struct RegionBoundary , low_y( _low_y ) , high_y( _high_y ) , low_z( _low_z ) - , high_z( _high_z ){}; + , high_z( _high_z ) + { + assert( low_x < high_x ); + assert( low_y < high_y ); + assert( low_z < high_z ); + }; }; template From e005fa1fc15a737e1a1c238ced05f97570542860 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 16 Oct 2023 09:26:10 -0400 Subject: [PATCH 21/34] fixup: file read inputs --- examples/read_file.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 1b3f0550..1f280fc1 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -113,7 +113,7 @@ int main( int argc, char* argv[] ) double G0 = 3.8; // [J/m^2] // PD horizon - double dx = 0.1299999999999999; + double dx = 0.5; double delta = 3 * dx; int halo_width = 1; @@ -136,17 +136,17 @@ int main( int argc, char* argv[] ) CabanaPD::Prenotch<0> prenotch; - CabanaPD::RegionBoundary planeUpper( -0.5, 0.5, 4.5, 5, -1, 1 ); - CabanaPD::RegionBoundary planeLower( -0.5, 0.5, -5, -4.5, -1, 1 ); + CabanaPD::RegionBoundary planeUpper( -0.5, 1.5, 8.5, 10.5, -1, 1 ); + CabanaPD::RegionBoundary planeLower( -0.5, 1.5, -0.5, 1.5, -1, 1 ); std::vector planes = { planeUpper, planeLower }; - int bc_dim = 2; + int bc_dim = 1; double center = particles->local_mesh_ext[bc_dim] / 2.0 + particles->local_mesh_lo[bc_dim]; auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes, - 2e6 / dx / dx, bc_dim, center ); + 2e6, bc_dim, center ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( From fab6caa5b3e7bfed33a640321ebb29d086c768bc Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 30 Oct 2023 16:48:50 -0400 Subject: [PATCH 22/34] fixup: allow displacement BC --- examples/read_file.cpp | 2 +- src/CabanaPD_Boundary.hpp | 17 ++++++----------- src/CabanaPD_Solver.hpp | 37 +++++++++++++++++++++++++++++++------ 3 files changed, 38 insertions(+), 18 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 1f280fc1..b5e9cf49 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -150,7 +150,7 @@ int main( int argc, char* argv[] ) // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, bc, prenotch ); + inputs, particles, force_model, bc, prenotch, false ); cabana_pd->init_force(); cabana_pd->run(); } diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 889700c5..0e98c1c1 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -165,11 +165,9 @@ struct BoundaryCondition _index_space.update( exec_space, particles, plane ); } - template - void apply( ExecSpace, ParticleType& particles ) + template + void apply( ExecSpace, FieldType& f, PositionType& ) { - auto f = particles.sliceForce(); - auto x = particles.sliceReferencePosition(); auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; @@ -202,10 +200,9 @@ struct BoundaryCondition _index_space.update( exec_space, particles, plane ); } - template - void apply( ExecSpace, ParticleType& particles ) + template + void apply( ExecSpace, FieldType& f, PositionType& ) { - auto f = particles.sliceForce(); auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; @@ -244,11 +241,9 @@ struct BoundaryCondition _index_space.update( exec_space, particles, plane ); } - template - void apply( ExecSpace, ParticleType& particles ) + template + void apply( ExecSpace, FieldType& f, PositionType& x ) { - auto f = particles.sliceForce(); - auto x = particles.sliceReferencePosition(); auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 81cf8934..88eef4e1 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -357,9 +357,11 @@ class SolverFracture SolverFracture( input_type _inputs, std::shared_ptr _particles, force_model_type force_model, bc_type bc, - prenotch_type prenotch ) + prenotch_type prenotch, + const bool force_boundary_condition ) : base_type( _inputs, _particles, force_model ) , boundary_condition( bc ) + , force_bc( force_boundary_condition ) { init_timer.reset(); @@ -391,8 +393,18 @@ class SolverFracture computeEnergy( *force, *particles, *neighbors, mu, neigh_iter_tag() ); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles ); - + auto x = particles->sliceReferencePosition(); + // FIXME: this needs to be generalized to any field. + if ( force_bc ) + { + auto f = particles->sliceForce(); + boundary_condition.apply( exec_space(), f, x ); + } + else + { + auto u = particles->sliceDisplacement(); + boundary_condition.apply( exec_space(), u, x ); + } particles->output( 0, 0.0, output_reference ); init_time += init_timer.seconds(); } @@ -436,7 +448,18 @@ class SolverFracture force_time += force_timer.seconds(); // Add boundary condition. - boundary_condition.apply( exec_space{}, *particles ); + auto x = particles->sliceReferencePosition(); + // FIXME: this needs to be generalized to any field. + if ( force_bc ) + { + auto f = particles->sliceForce(); + boundary_condition.apply( exec_space(), f, x ); + } + else + { + auto u = particles->sliceDisplacement(); + boundary_condition.apply( exec_space(), u, x ); + } // Integrate - velocity Verlet second half. integrate_timer.reset(); @@ -473,6 +496,7 @@ class SolverFracture using base_type::neighbors; using base_type::particles; bc_type boundary_condition; + bool force_bc; using NeighborView = typename Kokkos::View; NeighborView mu; @@ -510,11 +534,12 @@ template auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model, BCType bc, PrenotchType prenotch ) + ForceModel model, BCType bc, PrenotchType prenotch, + const bool force_bc = true ) { return std::make_shared>( - inputs, particles, model, bc, prenotch ); + inputs, particles, model, bc, prenotch, force_bc ); } /* From 773c3a5a48e8a6b88e9703dd82692cc9eef380a2 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 31 Oct 2023 11:42:38 -0400 Subject: [PATCH 23/34] fixup: remove dx for file reads --- examples/read_file.cpp | 3 +-- src/CabanaPD_Particles.hpp | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index b5e9cf49..c1ac1c95 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -115,7 +115,6 @@ int main( int argc, char* argv[] ) // PD horizon double dx = 0.5; double delta = 3 * dx; - int halo_width = 1; // Choose force model type. using model_type = @@ -132,7 +131,7 @@ int main( int argc, char* argv[] ) // Read particles from file. read_particles( inputs.input_file, *particles ); // Update after reading. Currently requires fixed cell spacing. - particles->updateAfterRead( exec_space(), halo_width, dx ); + particles->updateAfterRead( exec_space(), delta ); CabanaPD::Prenotch<0> prenotch; diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 9119ecdf..9922555f 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -292,9 +292,11 @@ class Particles // state. template std::enable_if_t<3 == NSD, void> - updateAfterRead( const ExecSpace& exec_space, const int hw, double dx ) + updateAfterRead( const ExecSpace& exec_space, const double cutoff ) { - halo_width = hw; + // Because this may be a non-uniform mesh, build a background mesh with + // dx=cutoff and a halo_width=1 + halo_width = 1; auto x = sliceReferencePosition(); n_local = x.size(); n_ghost = 0; @@ -341,19 +343,19 @@ class Particles std::array num_cells; for ( int d = 0; d < 3; d++ ) { - max_corner[d] += dx / 2.0; - min_corner[d] -= dx / 2.0; num_cells[d] = - static_cast( ( max_corner[d] - min_corner[d] ) / dx ); + static_cast( ( max_corner[d] - min_corner[d] ) / cutoff ); } createDomain( min_corner, max_corner, num_cells ); } template std::enable_if_t<2 == NSD, void> - updateAfterRead( const ExecSpace& exec_space, const int hw, double dx ) + updateAfterRead( const ExecSpace& exec_space, double cutoff ) { - halo_width = hw; + // Because this may be a non-uniform mesh, build a background mesh with + // dx=cutoff and a halo_width=1 + halo_width = 1; auto x = sliceReferencePosition(); n_local = x.size(); n_ghost = 0; @@ -390,10 +392,8 @@ class Particles std::array num_cells; for ( int d = 0; d < 2; d++ ) { - max_corner[d] += dx / 2.0; - min_corner[d] -= dx / 2.0; num_cells[d] = - static_cast( ( max_corner[d] - min_corner[d] ) / dx ); + static_cast( ( max_corner[d] - min_corner[d] ) / cutoff ); } createDomain( min_corner, max_corner, num_cells ); } From 905ae39f16001b27c787e30c9f66ac19cb5331b6 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 1 Nov 2023 14:32:43 -0400 Subject: [PATCH 24/34] fixup: force non-zero extent for pseudo-2d --- src/CabanaPD_Particles.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 9922555f..3cd4fb45 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -295,7 +295,7 @@ class Particles updateAfterRead( const ExecSpace& exec_space, const double cutoff ) { // Because this may be a non-uniform mesh, build a background mesh with - // dx=cutoff and a halo_width=1 + // dx=cutoff and force halo_width=1 halo_width = 1; auto x = sliceReferencePosition(); n_local = x.size(); @@ -345,6 +345,11 @@ class Particles { num_cells[d] = static_cast( ( max_corner[d] - min_corner[d] ) / cutoff ); + if ( num_cells[d] == 0 ) + { + num_cells[d]++; + max_corner[d] += 1e-16; + } } createDomain( min_corner, max_corner, num_cells ); } @@ -354,7 +359,7 @@ class Particles updateAfterRead( const ExecSpace& exec_space, double cutoff ) { // Because this may be a non-uniform mesh, build a background mesh with - // dx=cutoff and a halo_width=1 + // dx=cutoff and force halo_width=1 halo_width = 1; auto x = sliceReferencePosition(); n_local = x.size(); From 317abcd0e22095b7c51f3e7787d01b7b60f12eda Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 2 Nov 2023 12:19:19 -0400 Subject: [PATCH 25/34] fixup: memory space vs device type --- examples/read_file.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index ea0e1e2f..43975d7e 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -125,9 +125,8 @@ int main( int argc, char* argv[] ) inputs.read_args( argc, argv ); // Default construct to then read particles. - using device_type = Kokkos::Device; auto particles = std::make_shared>(); + memory_space, typename model_type::base_model, 3>>(); // Read particles from file. read_particles( inputs.input_file, *particles ); @@ -149,7 +148,7 @@ int main( int argc, char* argv[] ) 2e6, bc_dim, center ); // FIXME: use createSolver to switch backend at runtime. - auto cabana_pd = CabanaPD::createSolverFracture( + auto cabana_pd = CabanaPD::createSolverFracture( inputs, particles, force_model, bc, prenotch, false ); cabana_pd->init_force(); cabana_pd->run(); From 22245d9ead4aa86f2a27a832b57267292c741229 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 6 Nov 2023 17:33:49 -0500 Subject: [PATCH 26/34] fixup: BC regions --- examples/read_file.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 43975d7e..6522af35 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -135,8 +135,14 @@ int main( int argc, char* argv[] ) CabanaPD::Prenotch<0> prenotch; - CabanaPD::RegionBoundary planeUpper( -0.5, 1.5, 8.5, 10.5, -1, 1 ); - CabanaPD::RegionBoundary planeLower( -0.5, 1.5, -0.5, 1.5, -1, 1 ); + CabanaPD::RegionBoundary planeUpper( + particles->local_mesh_lo[0], particles->local_mesh_hi[0], + particles->local_mesh_hi[1] - delta, particles->local_mesh_hi[1], + particles->local_mesh_lo[2], particles->local_mesh_hi[2] ); + CabanaPD::RegionBoundary planeLower( + particles->local_mesh_lo[0], particles->local_mesh_hi[0], + particles->local_mesh_lo[1], particles->local_mesh_lo[1] + delta, + particles->local_mesh_lo[2], particles->local_mesh_hi[2] ); std::vector planes = { planeUpper, planeLower }; From bb1512cb6c08a1743ea72a276087da2bc643bd23 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 6 Nov 2023 17:35:28 -0500 Subject: [PATCH 27/34] Switch back to dx; issues with LinkedCellList (in VerletList) with single bin in a given dimension --- examples/read_file.cpp | 6 +++--- src/CabanaPD_Particles.hpp | 22 ++++++++++++++++------ 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 6522af35..6ecae5f7 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -113,9 +113,9 @@ int main( int argc, char* argv[] ) double G0 = 3.8; // [J/m^2] // PD horizon - double dx = 0.12; int m = 3; - double delta = m * dx; + double delta = 4.0; + double dx = delta / m; // Choose force model type. using model_type = @@ -131,7 +131,7 @@ int main( int argc, char* argv[] ) // Read particles from file. read_particles( inputs.input_file, *particles ); // Update after reading. Currently requires fixed cell spacing. - particles->updateAfterRead( exec_space(), delta ); + particles->updateAfterRead( exec_space(), dx ); CabanaPD::Prenotch<0> prenotch; diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 9983569a..55e33a29 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -288,7 +288,7 @@ class Particles // state. template std::enable_if_t<3 == NSD, void> - updateAfterRead( const ExecSpace& exec_space, const double cutoff ) + updateAfterRead( const ExecSpace& exec_space, const double dx ) { // Because this may be a non-uniform mesh, build a background mesh with // dx=cutoff and force halo_width=1 @@ -340,11 +340,13 @@ class Particles for ( int d = 0; d < 3; d++ ) { num_cells[d] = - static_cast( ( max_corner[d] - min_corner[d] ) / cutoff ); + static_cast( ( max_corner[d] - min_corner[d] ) / dx ); if ( num_cells[d] == 0 ) { num_cells[d]++; - max_corner[d] += 1e-16; + // Issues with pseudo-2d with only one bin in a given direction. + max_corner[d] += dx * 3; + min_corner[d] -= dx * 3; } } createDomain( min_corner, max_corner, num_cells ); @@ -352,7 +354,7 @@ class Particles template std::enable_if_t<2 == NSD, void> - updateAfterRead( const ExecSpace& exec_space, double cutoff ) + updateAfterRead( const ExecSpace& exec_space, double dx ) { // Because this may be a non-uniform mesh, build a background mesh with // dx=cutoff and force halo_width=1 @@ -394,9 +396,17 @@ class Particles for ( int d = 0; d < 2; d++ ) { num_cells[d] = - static_cast( ( max_corner[d] - min_corner[d] ) / cutoff ); + static_cast( ( max_corner[d] - min_corner[d] ) / dx ); + if ( num_cells[d] == 0 ) + { + num_cells[d]++; + // Potentially issues with pseudo-1d with only one bin in a + // given direction. + max_corner[d] += dx * 3; + min_corner[d] -= dx * 3; + } + createDomain( min_corner, max_corner, num_cells ); } - createDomain( min_corner, max_corner, num_cells ); } void update_global() From 45bd8682f06673dd831a376f191dbcb6885a3d56 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 6 Nov 2023 17:39:43 -0500 Subject: [PATCH 28/34] fixup: force BC changed to displacement --- examples/read_file.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 6ecae5f7..921b2829 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -151,7 +151,7 @@ int main( int argc, char* argv[] ) particles->local_mesh_lo[bc_dim]; auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes, - 2e6, bc_dim, center ); + 1e-4, bc_dim, center ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( From 4e6eb1a0e3fd0ff385e841a61f8076f0c993deba Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 6 Nov 2023 19:08:21 -0500 Subject: [PATCH 29/34] Add no-fail region --- examples/read_file.cpp | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 921b2829..5bca3bad 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -133,16 +133,30 @@ int main( int argc, char* argv[] ) // Update after reading. Currently requires fixed cell spacing. particles->updateAfterRead( exec_space(), dx ); + // Do this separately so that the mesh is already set up. + auto x = particles->sliceReferencePosition(); + auto nofail = particles->sliceNoFail(); + auto f = particles->sliceForce(); + double max_bc = particles->local_mesh_hi[1] - delta; + double min_bc = particles->local_mesh_lo[1] + delta; + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + if ( x( pid, 1 ) <= min_bc + 1e-10 || + x( pid, 1 ) >= max_bc - 1e-10 ) + nofail( pid ) = 1; + }; + particles->updateParticles( exec_space{}, init_functor ); + CabanaPD::Prenotch<0> prenotch; CabanaPD::RegionBoundary planeUpper( - particles->local_mesh_lo[0], particles->local_mesh_hi[0], - particles->local_mesh_hi[1] - delta, particles->local_mesh_hi[1], - particles->local_mesh_lo[2], particles->local_mesh_hi[2] ); + particles->local_mesh_lo[0], particles->local_mesh_hi[0], max_bc, + particles->local_mesh_hi[1], particles->local_mesh_lo[2], + particles->local_mesh_hi[2] ); CabanaPD::RegionBoundary planeLower( particles->local_mesh_lo[0], particles->local_mesh_hi[0], - particles->local_mesh_lo[1], particles->local_mesh_lo[1] + delta, - particles->local_mesh_lo[2], particles->local_mesh_hi[2] ); + particles->local_mesh_lo[1], min_bc, particles->local_mesh_lo[2], + particles->local_mesh_hi[2] ); std::vector planes = { planeUpper, planeLower }; @@ -151,7 +165,7 @@ int main( int argc, char* argv[] ) particles->local_mesh_lo[bc_dim]; auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes, - 1e-4, bc_dim, center ); + 1e-9, bc_dim, center ); // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( From 77b21968397af24fa0c1376fdde0afc63e69960a Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 10 Nov 2023 07:04:53 -0500 Subject: [PATCH 30/34] Add split values for symmetric BC --- examples/read_file.cpp | 5 ++--- src/CabanaPD_Boundary.hpp | 43 +++++++++++++++++++++++++++++++++------ 2 files changed, 39 insertions(+), 9 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 5bca3bad..04aeb732 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -165,11 +165,10 @@ int main( int argc, char* argv[] ) particles->local_mesh_lo[bc_dim]; auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes, - 1e-9, bc_dim, center ); + 2e4, 0.0, bc_dim, center ); - // FIXME: use createSolver to switch backend at runtime. auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, bc, prenotch, false ); + inputs, particles, force_model, bc, prenotch ); cabana_pd->init_force(); cabana_pd->run(); } diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 0e98c1c1..975190ce 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -220,14 +220,27 @@ struct BoundaryCondition template struct BoundaryCondition { - double _value; + double _value_top; + double _value_bottom; BCIndexSpace _index_space; int _dim; double _center; BoundaryCondition( const double value, BCIndexSpace bc_index_space, const int dim = 1, const double center = 0.0 ) - : _value( value ) + : _value_top( value ) + , _value_bottom( value ) + , _index_space( bc_index_space ) + , _dim( dim ) + , _center( center ) + { + } + + BoundaryCondition( const double value_top, const double value_bottom, + BCIndexSpace bc_index_space, const int dim = 1, + const double center = 0.0 ) + : _value_top( value_top ) + , _value_bottom( value_bottom ) , _index_space( bc_index_space ) , _dim( dim ) , _center( center ) @@ -246,15 +259,17 @@ struct BoundaryCondition { auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); - auto value = _value; auto dim = _dim; auto center = _center; + auto value_top = _value_top; + auto value_bottom = _value_bottom; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { auto pid = index_space( b ); - auto relative_x = x( pid, dim ) - center; - auto sign = std::abs( relative_x ) / relative_x; - f( pid, 1 ) += value * sign; + if ( x( pid, dim ) > center ) + f( pid, dim ) += value_top; + else + f( pid, dim ) += value_bottom; } ); } }; @@ -291,6 +306,22 @@ auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space, value, bc_indices, dim, center ); } +// FIXME: relatively large initial guess for allocation. +template +auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space, + Particles particles, + std::vector planes, + const double value_top, const double value_bottom, + const int dim, const double center, + const double initial_guess = 0.5 ) +{ + using memory_space = typename Particles::memory_space; + using bc_index_type = BoundaryIndexSpace; + bc_index_type bc_indices = createBoundaryIndexSpace( + exec_space, particles, planes, initial_guess ); + return BoundaryCondition( + value_top, value_bottom, bc_indices, dim, center ); +} } // namespace CabanaPD #endif From d47622edb7af8adef71f51a5fbebebfb414c290b Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 10 Nov 2023 08:09:26 -0500 Subject: [PATCH 31/34] Guard against duplicate points in file reading --- examples/read_file.cpp | 49 ++++++++++++++++++++++++++++++------------ 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 04aeb732..39181a30 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -66,29 +66,50 @@ void read_particles( const std::string filename, ParticleType& particles ) auto nofail = particles.sliceNoFail(); auto damage = particles.sliceDamage(); + // Using atomic memory for the count. + using CountView = + typename Kokkos::View>; + CountView count( "count" ); + using exec_space = typename memory_space::execution_space; Kokkos::parallel_for( "copy_to_particles", Kokkos::RangePolicy( 0, x.size() ), KOKKOS_LAMBDA( const int pid ) { + // Guard against duplicate points. This threaded loop should be + // faster than reading the data on the host. + for ( int p = 0; p < pid; p++ ) + { + if ( ( Kokkos::abs( x( p ) - x( pid ) ) < 1e-14 ) && + ( Kokkos::abs( y( p ) - y( pid ) ) < 1e-14 ) ) + return; + } + const std::size_t c = count()++; + // Set the particle position and volume. - pvol( pid ) = vol( pid ); - px( pid, 0 ) = x( pid ); - px( pid, 1 ) = y( pid ); + pvol( c ) = vol( pid ); + px( c, 0 ) = x( pid ); + px( c, 1 ) = y( pid ); // Initialize everything else to zero. // Currently psuedo-2d. - px( pid, 2 ) = 0.0; + px( c, 2 ) = 0.0; for ( int d = 0; d < 3; d++ ) { - u( pid, d ) = 0.0; - v( pid, d ) = 0.0; - f( pid, d ) = 0.0; + u( c, d ) = 0.0; + v( c, d ) = 0.0; + f( c, d ) = 0.0; } - type( pid ) = 0; - nofail( pid ) = 0; - rho( pid ) = 1.0; - damage( pid ) = 0; + type( c ) = 0; + nofail( c ) = 0; + rho( c ) = 1.0; + damage( c ) = 0; } ); + + int final_count = 0; + Kokkos::deep_copy( final_count, count ); + if ( final_count < static_cast( x.size() ) ) + particles.resize( final_count, 0 ); } int main( int argc, char* argv[] ) @@ -102,9 +123,9 @@ int main( int argc, char* argv[] ) using memory_space = typename exec_space::memory_space; // Time - double t_final = 1e-6; - double dt = 1e-9; - double output_frequency = 10; + double t_final = 1.0; + double dt = 1e-7; + double output_frequency = 100; // Material constants double E = 72e+9; // [Pa] From ab6f3e5cb9ea90e08b63eaafd16e7007d8cbb0dc Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 10 Nov 2023 09:39:48 -0500 Subject: [PATCH 32/34] Add time dependent BC with linear ramp --- examples/crack_branching.cpp | 8 +- examples/read_file.cpp | 9 +- src/CabanaPD_Boundary.hpp | 186 ++++++++++++++++++++++++++--------- src/CabanaPD_Solver.hpp | 9 +- 4 files changed, 154 insertions(+), 58 deletions(-) diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp index d71812a4..cc045cc3 100644 --- a/examples/crack_branching.cpp +++ b/examples/crack_branching.cpp @@ -106,9 +106,11 @@ int main( int argc, char* argv[] ) CabanaPD::RegionBoundary plane2( low_x, high_x, high_y - dy, high_y + dy, low_z, high_z ); std::vector planes = { plane1, plane2 }; - auto bc = - createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, - exec_space{}, *particles, planes, b0 ); + int bc_dim = 1; + double center = 0.0; + auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, + exec_space{}, *particles, planes, b0, + bc_dim, center ); auto init_functor = KOKKOS_LAMBDA( const int pid ) { diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 39181a30..8d7fe8bb 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -123,7 +123,8 @@ int main( int argc, char* argv[] ) using memory_space = typename exec_space::memory_space; // Time - double t_final = 1.0; + double t_final = 1e-1; + double t_ramp = 1e-3; double dt = 1e-7; double output_frequency = 100; @@ -184,9 +185,9 @@ int main( int argc, char* argv[] ) int bc_dim = 1; double center = particles->local_mesh_ext[bc_dim] / 2.0 + particles->local_mesh_lo[bc_dim]; - auto bc = createBoundaryCondition( CabanaPD::ForceSymmetric1dBCTag{}, - exec_space{}, *particles, planes, - 2e4, 0.0, bc_dim, center ); + auto bc = createBoundaryCondition( + CabanaPD::ForceSymmetric1dBCTag{}, exec_space{}, *particles, planes, + 2e4, 0.0, bc_dim, center, 0.0, t_ramp ); auto cabana_pd = CabanaPD::createSolverFracture( inputs, particles, force_model, bc, prenotch ); diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 975190ce..32912c42 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -151,10 +151,19 @@ struct BoundaryCondition { double _value; BCIndexSpace _index_space; - - BoundaryCondition( const double value, BCIndexSpace bc_index_space ) + double _time_start; + double _time_ramp; + double _time_end; + + BoundaryCondition( BCIndexSpace bc_index_space, const double value, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) : _value( value ) , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) { } @@ -166,16 +175,28 @@ struct BoundaryCondition } template - void apply( ExecSpace, FieldType& f, PositionType& ) + void apply( ExecSpace, FieldType& f, PositionType&, const double t ) { auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; + auto start = _time_start; + auto end = _time_end; + auto ramp = _time_ramp; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { - auto pid = index_space( b ); - for ( int d = 0; d < 3; d++ ) - f( pid, d ) = value; + double current = value; + // Initial linear ramp. + if ( t < ramp ) + { + current = value * ( t - start ) / ( ramp - start ); + } + if ( t > start && t < end ) + { + auto pid = index_space( b ); + for ( int d = 0; d < 3; d++ ) + f( pid, d ) = current; + } } ); } }; @@ -186,11 +207,23 @@ struct BoundaryCondition { double _value; BCIndexSpace _index_space; - - BoundaryCondition( const double value, BCIndexSpace bc_index_space ) + double _time_start; + double _time_ramp; + double _time_end; + + BoundaryCondition( BCIndexSpace bc_index_space, const double value, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) : _value( value ) , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) { + assert( _time_ramp >= _time_start ); + assert( _time_end >= _time_ramp ); + assert( _time_end >= _time_start ); } template @@ -201,16 +234,28 @@ struct BoundaryCondition } template - void apply( ExecSpace, FieldType& f, PositionType& ) + void apply( ExecSpace, FieldType& f, PositionType&, const double t ) { auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); auto value = _value; + auto start = _time_start; + auto end = _time_end; + auto ramp = _time_ramp; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { - auto pid = index_space( b ); - for ( int d = 0; d < 3; d++ ) - f( pid, d ) += value; + double current = value; + // Initial linear ramp. + if ( t < ramp ) + { + current = value * ( t - start ) / ( ramp - start ); + } + if ( t > start && t < end ) + { + auto pid = index_space( b ); + for ( int d = 0; d < 3; d++ ) + f( pid, d ) += current; + } } ); } }; @@ -220,31 +265,56 @@ struct BoundaryCondition template struct BoundaryCondition { + BCIndexSpace _index_space; double _value_top; double _value_bottom; - BCIndexSpace _index_space; + double _time_start; + double _time_ramp; + double _time_end; + int _dim; double _center; - BoundaryCondition( const double value, BCIndexSpace bc_index_space, - const int dim = 1, const double center = 0.0 ) - : _value_top( value ) + BoundaryCondition( BCIndexSpace bc_index_space, const double value, + const int dim, const double center, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) + : _index_space( bc_index_space ) + , _value_top( value ) , _value_bottom( value ) - , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) , _dim( dim ) , _center( center ) { + assert( _time_ramp >= _time_start ); + assert( _time_end >= _time_ramp ); + assert( _time_end >= _time_start ); + assert( _dim >= 0 ); + assert( _dim < 3 ); } - BoundaryCondition( const double value_top, const double value_bottom, - BCIndexSpace bc_index_space, const int dim = 1, - const double center = 0.0 ) - : _value_top( value_top ) + BoundaryCondition( BCIndexSpace bc_index_space, const double value_top, + const double value_bottom, const int dim, + const double center, const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max() ) + : _index_space( bc_index_space ) + , _value_top( value_top ) , _value_bottom( value_bottom ) - , _index_space( bc_index_space ) + , _time_start( start ) + , _time_ramp( ramp ) + , _time_end( end ) , _dim( dim ) , _center( center ) { + assert( _time_ramp >= _time_start ); + assert( _time_end >= _time_ramp ); + assert( _time_end >= _time_start ); + assert( _dim >= 0 ); + assert( _dim < 3 ); } template @@ -255,7 +325,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, FieldType& f, PositionType& x ) + void apply( ExecSpace, FieldType& f, PositionType& x, const double t ) { auto index_space = _index_space._view; Kokkos::RangePolicy policy( 0, index_space.size() ); @@ -263,13 +333,28 @@ struct BoundaryCondition auto center = _center; auto value_top = _value_top; auto value_bottom = _value_bottom; + auto start = _time_start; + auto end = _time_end; + auto ramp = _time_ramp; Kokkos::parallel_for( "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) { - auto pid = index_space( b ); - if ( x( pid, dim ) > center ) - f( pid, dim ) += value_top; - else - f( pid, dim ) += value_bottom; + auto current_top = value_top; + auto current_bottom = value_bottom; + // Initial linear ramp. + if ( t < ramp ) + { + auto t_factor = ( t - start ) / ( ramp - start ); + current_top = value_top * t_factor; + current_bottom = value_bottom * t_factor; + } + if ( t > start && t < end ) + { + auto pid = index_space( b ); + if ( x( pid, dim ) > center ) + f( pid, dim ) += current_top; + else + f( pid, dim ) += current_bottom; + } } ); } }; @@ -277,50 +362,57 @@ struct BoundaryCondition // FIXME: relatively large initial guess for allocation. template -auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles, - std::vector planes, - const double value, - const double initial_guess = 0.5 ) +auto createBoundaryCondition( + BCTag, ExecSpace exec_space, Particles particles, + std::vector planes, const double value, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max(), + const double initial_guess = 0.5 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); - return BoundaryCondition( value, bc_indices ); + return BoundaryCondition( bc_indices, value, start, + ramp, end ); } // FIXME: relatively large initial guess for allocation. template -auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space, - Particles particles, - std::vector planes, - const double value, const int dim, - const double center, - const double initial_guess = 0.5 ) +auto createBoundaryCondition( + ForceSymmetric1dBCTag, ExecSpace exec_space, Particles particles, + std::vector planes, const double value, const int dim, + const double center, const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max(), + const double initial_guess = 0.5 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); return BoundaryCondition( - value, bc_indices, dim, center ); + bc_indices, value, dim, center, start, ramp, end ); } // FIXME: relatively large initial guess for allocation. template -auto createBoundaryCondition( ForceSymmetric1dBCTag, ExecSpace exec_space, - Particles particles, - std::vector planes, - const double value_top, const double value_bottom, - const int dim, const double center, - const double initial_guess = 0.5 ) +auto createBoundaryCondition( + ForceSymmetric1dBCTag, ExecSpace exec_space, Particles particles, + std::vector planes, const double value_top, + const double value_bottom, const int dim, const double center, + const double start = 0.0, + const double ramp = std::numeric_limits::max(), + const double end = std::numeric_limits::max(), + const double initial_guess = 0.5 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); return BoundaryCondition( - value_top, value_bottom, bc_indices, dim, center ); + bc_indices, value_top, value_bottom, dim, center, start, ramp, end ); } } // namespace CabanaPD diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 2157db8c..800e78d3 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -398,12 +398,12 @@ class SolverFracture if ( force_bc ) { auto f = particles->sliceForce(); - boundary_condition.apply( exec_space(), f, x ); + boundary_condition.apply( exec_space(), f, x, 0.0 ); } else { auto u = particles->sliceDisplacement(); - boundary_condition.apply( exec_space(), u, x ); + boundary_condition.apply( exec_space(), u, x, 0.0 ); } particles->output( 0, 0.0, output_reference ); init_time += init_timer.seconds(); @@ -450,15 +450,16 @@ class SolverFracture // Add boundary condition. auto x = particles->sliceReferencePosition(); // FIXME: this needs to be generalized to any field. + auto time = step * inputs->timestep; if ( force_bc ) { auto f = particles->sliceForce(); - boundary_condition.apply( exec_space(), f, x ); + boundary_condition.apply( exec_space(), f, x, time ); } else { auto u = particles->sliceDisplacement(); - boundary_condition.apply( exec_space(), u, x ); + boundary_condition.apply( exec_space(), u, x, time ); } // Integrate - velocity Verlet second half. From e4107f8c7ec12f6efae6731150390e24fbb5cd54 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 10 Nov 2023 12:21:57 -0500 Subject: [PATCH 33/34] fixup: file read abs() --- examples/read_file.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/read_file.cpp b/examples/read_file.cpp index 8d7fe8bb..e04f95b6 100644 --- a/examples/read_file.cpp +++ b/examples/read_file.cpp @@ -77,11 +77,13 @@ void read_particles( const std::string filename, ParticleType& particles ) "copy_to_particles", Kokkos::RangePolicy( 0, x.size() ), KOKKOS_LAMBDA( const int pid ) { // Guard against duplicate points. This threaded loop should be - // faster than reading the data on the host. + // faster than checking during the host read. for ( int p = 0; p < pid; p++ ) { - if ( ( Kokkos::abs( x( p ) - x( pid ) ) < 1e-14 ) && - ( Kokkos::abs( y( p ) - y( pid ) ) < 1e-14 ) ) + auto xdiff = x( p ) - x( pid ); + auto ydiff = y( p ) - y( pid ); + if ( ( Kokkos::abs( xdiff ) < 1e-14 ) && + ( Kokkos::abs( ydiff ) < 1e-14 ) ) return; } const std::size_t c = count()++; From 89c9e6c5e7288b8a379baeea13c47e9cfc80a224 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 10 Nov 2023 12:35:16 -0500 Subject: [PATCH 34/34] Update HIP CI version --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 17cc5f92..4faf395d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -118,7 +118,7 @@ jobs: matrix: cxx: ['hipcc'] cmake_build_type: ['Release'] - kokkos_ver: ['3.6.01'] + kokkos_ver: ['4.0.01'] runs-on: ubuntu-20.04 container: ghcr.io/ecp-copa/ci-containers/rocm:latest steps: