Enable particle input from file

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})

examples/read_file.cpp

@@ -0,0 +1,199 @@
+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+template <class ParticleType>
+void read_particles( const std::string filename, ParticleType& particles )
+{
+ // Read particles from csv file.
+ std::vector<double> csv_x;
+ std::vector<double> csv_y;
+ std::vector<double> csv_v;
+
+ std::vector<std::string> row;
+ std::string line, word;
+
+ std::fstream file( filename, 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_v.push_back( std::stod( row[3] ) );
+ }
+ }
+ else
+ throw( "Could not open file." + filename );
+
+ // Create unmanaged Views in order to copy to device.
+ Kokkos::View<double*, Kokkos::HostSpace> x_host( csv_x.data() );
+ Kokkos::View<double*, Kokkos::HostSpace> y_host( csv_y.data() );
+ Kokkos::View<double*, Kokkos::HostSpace> 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<exec_space> 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<int, 3> 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<double, 3> low_corner = { low_x, low_y, low_z };
+ std::array<double, 3> 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<double, 3> p01 = { low_x, y_prenotch1, low_z };
+ Kokkos::Array<double, 3> v1 = { L_prenotch, 0, 0 };
+ Kokkos::Array<double, 3> v2 = { 0, 0, thickness };
+ Kokkos::Array<Kokkos::Array<double, 3>, 1> notch_positions = { p01 };
+ CabanaPD::Prenotch<1> prenotch( v1, v2, notch_positions );
+
+ // Choose force model type.
+ using model_type =
+ CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
+ 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<exec_space, memory_space>;
+ auto particles = std::make_shared<CabanaPD::Particles<
+ 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 );
+
+ // 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<CabanaPD::RegionBoundary> 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<device_type>(
+ inputs, particles, force_model, bc, prenotch );
+ cabana_pd->init_force();
+ cabana_pd->run();
+ }
+
+ MPI_Finalize();
+}

src/CabanaPD_Particles.hpp

@@ -283,6 +283,54 @@ class Particles<DeviceType, PMB, Dimension>
  resize( n_local, 0 );
  size = _aosoa_x.size();
 
+ update_global();
+ }
+
+ // This is necessary after reading in particles from file for a consistent
+ // state.
+ template <class ExecSpace>
+ void update_after_read( const ExecSpace& exec_space, const int hw,
+ const std::array<int, 3> num_cells )
+ {
+ halo_width = hw;
+ auto x = slice_x();
+ 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<ExecSpace>( 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 );
+
+ 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 );