diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index df1b8e15..224ef71a 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) -install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(RandomCracks random_cracks.cpp) +target_link_libraries(RandomCracks LINK_PUBLIC CabanaPD) + +install(TARGETS ElasticWave KalthoffWinkler CrackBranching RandomCracks DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/random_cracks.cpp b/examples/random_cracks.cpp new file mode 100644 index 00000000..5d6dda8a --- /dev/null +++ b/examples/random_cracks.cpp @@ -0,0 +1,191 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include +#include +#include + +#include + +#include "mpi.h" + +#include + +#include + +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 = { 400, 160, 8 }; // 400 x 160 x 8 + 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 = 5e-8; + double output_frequency = 5; + bool output_reference = true; + + // 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 + 1e-10; + + // 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; + + // Multiple random pre-notches + + // Number of pre-notches + constexpr int Npn = 200; + + // Minimum and maximum pre-notch length + double minl = 0.001; + double maxl = 0.0025; + + // Initialize pre-notch arrays + Kokkos::Array, Npn> notch_positions; + Kokkos::Array, Npn> notch_v1; + Kokkos::Array, Npn> notch_v2; + + // Reference for random number generator: + // https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution + std::random_device rd; + std::mt19937 gen( rd() ); + std::uniform_real_distribution<> dis( 0.0, 1.0 ); + + // Loop over pre-notches + for ( int n = 0; n < Npn; n++ ) + { + // Random numbers for pre-notch position + double random_number_x = dis( gen ); + double random_number_y = dis( gen ); + + // Coordinates of one endpoint of the pre-notch (random) + // Note: the addition and subtraction of "maxl" ensures the prenotch + // does not extend outside the domain + double Xc1 = + ( low_x + maxl ) + + ( ( high_x - maxl ) - ( low_x + maxl ) ) * random_number_x; + double Yc1 = + ( low_y + maxl ) + + ( ( high_y - maxl ) - ( low_y + maxl ) ) * random_number_y; + Kokkos::Array p0 = { Xc1, Yc1, low_z }; + + // Assign pre-notch position + notch_positions[n] = p0; + + // Random pre-notch length on XY-plane + double random_number_l = dis( gen ); + double l = minl + ( maxl - minl ) * random_number_l; + + // Random pre-notch orientation on XY-plane + double random_number_theta = dis( gen ); + double theta = M_PI * random_number_theta; + + // Pre-notch v1 vector + Kokkos::Array v1 = { l * cos( theta ), l * sin( theta ), + 0 }; + notch_v1[n] = v1; + + // Pre-notch v2 vector + + // Random number for y-component of v2: the angle of v2 in the + // YZ-plane is between 0 and 45 deg. + double random_number_v2_y = dis( gen ); + Kokkos::Array v2 = { 0, random_number_v2_y * thickness, + thickness }; + notch_v2[n] = v2; + } + + CabanaPD::Prenotch prenotch( notch_v1, notch_v2, notch_positions ); + + // Choose force model type. + using model_type = + CabanaPD::ForceModel; + model_type force_model( delta, K, G0 ); + CabanaPD::Inputs<3> inputs( num_cell, low_corner, high_corner, t_final, + dt, output_frequency, output_reference ); + inputs.read_args( argc, argv ); + + // Create particles from mesh. + // Does not set displacements, velocities, etc. + auto particles = std::make_shared< + CabanaPD::Particles>( + exec_space(), inputs.low_corner, inputs.high_corner, + inputs.num_cells, halo_width ); + + // Define particle initialization. + auto x = particles->sliceReferencePosition(); + auto v = particles->sliceVelocity(); + auto f = particles->sliceForce(); + auto rho = particles->sliceDensity(); + auto nofail = particles->sliceNoFail(); + + // Relying on uniform grid here. + double dy = particles->dx[1]; + double b0 = 2e6 / dy; // Pa/m + + 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->updateParticles( 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_Prenotch.hpp b/src/CabanaPD_Prenotch.hpp index d29aa794..def35bc4 100644 --- a/src/CabanaPD_Prenotch.hpp +++ b/src/CabanaPD_Prenotch.hpp @@ -195,22 +195,43 @@ int bondPrenotchIntersection( const Kokkos::Array v1, return keep_bond; } -template +template struct Prenotch { static constexpr std::size_t num_notch = NumNotch; - Kokkos::Array _v1; - Kokkos::Array _v2; - Kokkos::Array, num_notch> _p0_list; + static constexpr std::size_t num_vector = NumNotch; + Kokkos::Array, num_vector> _v1; + Kokkos::Array, num_vector> _v2; + Kokkos::Array, num_notch> _p0; + bool fixed_orientation; + // Default constructor Prenotch() {} + // Constructor if all pre-notches are oriented the same way (e.g. + // Kalthoff-Winkler). Prenotch( Kokkos::Array v1, Kokkos::Array v2, - Kokkos::Array, num_notch> p0_list ) + Kokkos::Array, num_notch> p0 ) + : _v1( { v1 } ) + , _v2( { v2 } ) + , _p0( p0 ) + { + fixed_orientation = true; + } + + // Constructor for general case of any orientation for any number of + // pre-notches. + Prenotch( Kokkos::Array, num_vector> v1, + Kokkos::Array, num_vector> v2, + Kokkos::Array, num_notch> p0 ) : _v1( v1 ) , _v2( v2 ) - , _p0_list( p0_list ) + , _p0( p0 ) { + static_assert( + num_vector == num_notch, + "Number of orientation vectors must match number of pre-notches." ); + fixed_orientation = false; } template policy( 0, particles.n_local ); - auto v1 = _v1; - auto v2 = _v2; - for ( std::size_t p = 0; p < _p0_list.size(); p++ ) + for ( std::size_t p = 0; p < _p0.size(); p++ ) { - auto p0 = _p0_list[p]; + // These will always be different positions. + auto p0 = _p0[p]; + // These may all have the same orientation or all different + // orientations. + auto v1 = getV1( p ); + auto v2 = getV2( p ); auto notch_functor = KOKKOS_LAMBDA( const int i ) { std::size_t num_neighbors = @@ -252,6 +276,22 @@ struct Prenotch Kokkos::parallel_for( "CabanaPD::Prenotch", policy, notch_functor ); } } + + auto getV1( const int p ) + { + if ( fixed_orientation ) + return _v1[0]; + else + return _v1[p]; + } + + auto getV2( const int p ) + { + if ( fixed_orientation ) + return _v2[0]; + else + return _v2[p]; + } }; } // namespace CabanaPD