diff --git a/examples/mechanics/CMakeLists.txt b/examples/mechanics/CMakeLists.txt index a810f419..ad346f7c 100644 --- a/examples/mechanics/CMakeLists.txt +++ b/examples/mechanics/CMakeLists.txt @@ -22,4 +22,16 @@ target_link_libraries(PlateWithHole LINK_PUBLIC CabanaPD) add_executable(CrackInclusion crack_inclusion.cpp) target_link_libraries(CrackInclusion LINK_PUBLIC CabanaPD) -install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder RandomCracks DogboneTensileTest PlateWithHole CrackInclusion DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(FlatIndenter flat_indenter.cpp) +target_link_libraries(FlatIndenter LINK_PUBLIC CabanaPD) + +add_executable(SphericalIndenter spherical_indenter.cpp) +target_link_libraries(SphericalIndenter LINK_PUBLIC CabanaPD) + +add_executable(ConicalIndenter conical_indenter.cpp) +target_link_libraries(ConicalIndenter LINK_PUBLIC CabanaPD) + +add_executable(BerkovichIndenter berkovich_indenter.cpp) +target_link_libraries(BerkovichIndenter LINK_PUBLIC CabanaPD) + +install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder RandomCracks DogboneTensileTest PlateWithHole CrackInclusion FlatIndenter SphericalIndenter ConicalIndenter BerkovichIndenter DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/mechanics/berkovich_indenter.cpp b/examples/mechanics/berkovich_indenter.cpp new file mode 100644 index 00000000..ee369651 --- /dev/null +++ b/examples/mechanics/berkovich_indenter.cpp @@ -0,0 +1,276 @@ +/**************************************************************************** + * Copyright (c) 2022 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 "mpi.h" + +#include + +#include + +// Simulate a Berkovich indenter on a plate. +void berkovichIndenterExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double horizon = inputs["horizon"]; + horizon += 1e-10; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + + // ==================================================== + // Force model + // ==================================================== + using model_type = CabanaPD::PMB; + // using contact_type = CabanaPD::NormalRepulsionModel; + CabanaPD::ForceModel force_model( model_type{}, horizon, K, G0 ); + + // ==================================================== + // Particle generation + // ==================================================== + CabanaPD::Particles particles( memory_space{}, model_type{} ); + // CabanaPD::Particles particles( memory_space{}, contact_type{} ); + + // Note that individual inputs can be passed instead (see other examples). + particles.domain( inputs ); + particles.create( exec_space{} ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + double dx = particles.dx[0]; + double dy = particles.dx[1]; + double dz = particles.dx[2]; + double alpha = inputs["indenter_angle"]; + // Convert angle to radians. + alpha *= CabanaPD::pi / 180.0; + // Compute tan( alpha ) to use multiple times. + double tan_alpha = Kokkos::tan( alpha ); + + double v0 = inputs["indenter_velocity"]; + double tfinal = inputs["final_time"]; + // Height + double H = v0 * tfinal; + // Base radius + double R = H * tan_alpha; + + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + + // TODO: Adapt region to pyramid base shape for force calculation. + // Define boundary condition regions to constrain displacements. + CabanaPD::Region pressure_region( + low_corner[0], high_corner[0], low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + + // Define boundary condition regions to constrain displacements. + CabanaPD::Region low_x_plane( + low_corner[0] - dx, low_corner[0] + dx, low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + CabanaPD::Region high_x_plane( + high_corner[0] - dx, high_corner[0] + dx, low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + CabanaPD::Region low_y_plane( + low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, + low_corner[2], high_corner[2] ); + CabanaPD::Region high_y_plane( + low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy, + low_corner[2], high_corner[2] ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles.sliceDensity(); + auto nofail = particles.sliceNoFail(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + nofail( pid ) = 1; + }; + particles.update( exec_space{}, init_functor ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + + /* + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx / 2.0; + r_extend *= dx; + + contact_type contact_model( horizon, r_c, r_extend, K ); + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + */ + + // ==================================================== + // Boundary conditions + // ==================================================== + // Create BC last to ensure ghost particles are included. + auto x = solver.particles.sliceReferencePosition(); + auto u = solver.particles.sliceDisplacement(); + // double v0 = inputs["indenter_velocity"]; + + // z-coordinate of top layer of particles + double z_top = high_corner[2] - 0.5 * dz; + // Initial z-coordinate of the indenter tip + double z0_indenter = z_top; + + auto disp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + // z-coordinate of the indenter tip + double z_indenter = z0_indenter - v0 * t; + + // Check if indenter tip touches the plate (top layer of particles) + if ( z_indenter < z_top ) + { + // Height of indenter + double H = z_top - z_indenter; + + // Edge length of base equilateral triangle + double L = 2.0 * std::sqrt( 3.0 ) * tan_alpha * H; + + // Inradius of base equilateral triangle + double R = std::sqrt( 3.0 ) * L / 6.0; + + // Indenter displacement + double x_i = x( pid, 0 ); + double y_i = x( pid, 1 ); + + if ( y_i > -R && y_i <= std::sqrt( 3.0 ) * x_i + 2 * R && + y_i <= -std::sqrt( 3.0 ) * x_i + 2 * R ) + { + // Find distance to edges. + double d1 = std::abs( y_i + R ); + double d2 = + 0.5 * std::abs( -std::sqrt( 3.0 ) * x_i + y_i - 2 * R ); + double d3 = + 0.5 * std::abs( std::sqrt( 3.0 ) * x_i + y_i - 2 * R ); + + // Find distance to nearest edge. + double d = std::min( { d1, d2, d3 } ); + + // Find inradius of triangle for which particle is on its + // perimenter + double r = R - d; + + // Find edge length of triangle for which particle is on its + // perimenter + double l = 6.0 * r / std::sqrt( 3.0 ); + + // Indenter displacement + u( pid, 2 ) = z_indenter + + l / ( 2.0 * std::sqrt( 3.0 ) * tan_alpha ) - + x( pid, 2 ); + } + } + + // Edge constraints + if ( low_x_plane.inside( x, pid ) || high_x_plane.inside( x, pid ) ) + { + u( pid, 0 ) = 0; + u( pid, 1 ) = 0; + u( pid, 2 ) = 0; + } + + if ( low_y_plane.inside( x, pid ) || high_y_plane.inside( x, pid ) ) + { + u( pid, 0 ) = 0; + u( pid, 1 ) = 0; + u( pid, 2 ) = 0; + } + }; + + auto bc = createBoundaryCondition( + disp_func, exec_space{}, solver.particles, true, pressure_region, + low_x_plane, high_x_plane, low_y_plane, high_y_plane ); + + //======================================== + // OUTPUTS + //======================================== + auto f = solver.particles.sliceForce(); + + // Output force in x-direction. + auto force_func_x = KOKKOS_LAMBDA( const int p ) + { + return f( p, 0 ) * dx * dy * dz; + }; + auto output_fx = CabanaPD::createOutputTimeSeries( + "output_force_x.txt", inputs, exec_space{}, solver.particles, + force_func_x, pressure_region ); + + // Output force in y-direction. + auto force_func_y = KOKKOS_LAMBDA( const int p ) + { + return f( p, 1 ) * dx * dy * dz; + }; + auto output_fy = CabanaPD::createOutputTimeSeries( + "output_force_y.txt", inputs, exec_space{}, solver.particles, + force_func_y, pressure_region ); + + // Output force in z-direction. + auto force_func_z = KOKKOS_LAMBDA( const int p ) + { + return f( p, 2 ) * dx * dy * dz; + }; + auto output_fz = CabanaPD::createOutputTimeSeries( + "output_force_z.txt", inputs, exec_space{}, solver.particles, + force_func_z, pressure_region ); + + // ==================================================== + // Simulation run + // ==================================================== + solver.init( bc ); + + solver.updateRegion( output_fx, output_fy, output_fz ); + + solver.run( bc, output_fx, output_fy, output_fz ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + berkovichIndenterExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/mechanics/conical_indenter.cpp b/examples/mechanics/conical_indenter.cpp new file mode 100644 index 00000000..a5c4640d --- /dev/null +++ b/examples/mechanics/conical_indenter.cpp @@ -0,0 +1,247 @@ +/**************************************************************************** + * Copyright (c) 2022 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 "mpi.h" + +#include + +#include + +// Simulate a conical indenter on a plate. +void conicalIndenterExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double horizon = inputs["horizon"]; + horizon += 1e-10; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + + // ==================================================== + // Force model + // ==================================================== + using model_type = CabanaPD::PMB; + // using contact_type = CabanaPD::NormalRepulsionModel; + CabanaPD::ForceModel force_model( model_type{}, horizon, K, G0 ); + + // ==================================================== + // Particle generation + // ==================================================== + CabanaPD::Particles particles( memory_space{}, model_type{} ); + // CabanaPD::Particles particles( memory_space{}, contact_type{} ); + + // Note that individual inputs can be passed instead (see other examples). + particles.domain( inputs ); + particles.create( exec_space{} ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + double dx = particles.dx[0]; + double dy = particles.dx[1]; + double dz = particles.dx[2]; + double alpha = inputs["indenter_angle"]; + // Convert angle to radians. + alpha *= CabanaPD::pi / 180.0; + // Compute tan( alpha ) to use multiple times. + double tan_alpha = Kokkos::tan( alpha ); + + double v0 = inputs["indenter_velocity"]; + double tfinal = inputs["final_time"]; + // Height + double H = v0 * tfinal; + // Base radius + double R = H * tan_alpha; + + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + + // Define cylindrical region to apply pressure. + CabanaPD::Region pressure_region( + 0.0, R, high_corner[2] - dz, high_corner[2] + dz, x_center, y_center ); + + // Define boundary condition regions to constrain displacements. + CabanaPD::Region low_x_plane( + low_corner[0] - dx, low_corner[0] + dx, low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + CabanaPD::Region high_x_plane( + high_corner[0] - dx, high_corner[0] + dx, low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + CabanaPD::Region low_y_plane( + low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, + low_corner[2], high_corner[2] ); + CabanaPD::Region high_y_plane( + low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy, + low_corner[2], high_corner[2] ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles.sliceDensity(); + auto nofail = particles.sliceNoFail(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + nofail( pid ) = 1; + }; + particles.update( exec_space{}, init_functor ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + + /* + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx / 2.0; + r_extend *= dx; + + contact_type contact_model( horizon, r_c, r_extend, K ); + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + */ + + // ==================================================== + // Boundary conditions + // ==================================================== + // Create BC last to ensure ghost particles are included. + auto x = solver.particles.sliceReferencePosition(); + auto u = solver.particles.sliceDisplacement(); + // double v0 = inputs["indenter_velocity"]; + + // z-coordinate of top layer of particles + double z_top = high_corner[2] - 0.5 * dz; + // Initial z-coordinate of the indenter tip + double z0_indenter = z_top; + + auto disp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + // z-coordinate of the indenter tip + double z_indenter = z0_indenter - v0 * t; + + // Check if indenter tip touches the plate (top layer of particles) + if ( z_indenter < z_top ) + { + double r_indenter = ( z_top - z_indenter ) * tan_alpha; + + // Distance squared of particle to center on XY plane (assumes + // indenter is centered on the XY plane) + double r_sq = + ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + + // Indenter displacement + if ( r_sq <= r_indenter * r_indenter ) + u( pid, 2 ) = ( z_indenter + std::sqrt( r_sq ) / tan_alpha ) - + x( pid, 2 ); + } + + // Edge constraints + if ( low_x_plane.inside( x, pid ) || high_x_plane.inside( x, pid ) ) + { + u( pid, 0 ) = 0; + u( pid, 1 ) = 0; + u( pid, 2 ) = 0; + } + + if ( low_y_plane.inside( x, pid ) || high_y_plane.inside( x, pid ) ) + { + u( pid, 0 ) = 0; + u( pid, 1 ) = 0; + u( pid, 2 ) = 0; + } + }; + + auto bc = createBoundaryCondition( + disp_func, exec_space{}, solver.particles, true, pressure_region, + low_x_plane, high_x_plane, low_y_plane, high_y_plane ); + + //======================================== + // OUTPUTS + //======================================== + auto f = solver.particles.sliceForce(); + + // Output force in x-direction. + auto force_func_x = KOKKOS_LAMBDA( const int p ) + { + return f( p, 0 ) * dx * dy * dz; + }; + auto output_fx = CabanaPD::createOutputTimeSeries( + "output_force_x.txt", inputs, exec_space{}, solver.particles, + force_func_x, pressure_region ); + + // Output force in y-direction. + auto force_func_y = KOKKOS_LAMBDA( const int p ) + { + return f( p, 1 ) * dx * dy * dz; + }; + auto output_fy = CabanaPD::createOutputTimeSeries( + "output_force_y.txt", inputs, exec_space{}, solver.particles, + force_func_y, pressure_region ); + + // Output force in z-direction. + auto force_func_z = KOKKOS_LAMBDA( const int p ) + { + return f( p, 2 ) * dx * dy * dz; + }; + auto output_fz = CabanaPD::createOutputTimeSeries( + "output_force_z.txt", inputs, exec_space{}, solver.particles, + force_func_z, pressure_region ); + + // ==================================================== + // Simulation run + // ==================================================== + solver.init( bc ); + + solver.updateRegion( output_fx, output_fy, output_fz ); + + solver.run( bc, output_fx, output_fy, output_fz ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + conicalIndenterExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/mechanics/flat_indenter.cpp b/examples/mechanics/flat_indenter.cpp new file mode 100644 index 00000000..69a41e10 --- /dev/null +++ b/examples/mechanics/flat_indenter.cpp @@ -0,0 +1,164 @@ +/**************************************************************************** + * Copyright (c) 2022 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 "mpi.h" + +#include + +#include + +// Simulate crack branching from an pre-crack. +void flatIndenterExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double horizon = inputs["horizon"]; + horizon += 1e-10; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + + // ==================================================== + // Pre-notch + // ==================================================== +/* + double height = inputs["system_size"][0]; + double thickness = inputs["system_size"][2]; + double L_prenotch = height / 2.0; + double y_prenotch = 0.0; + Kokkos::Array p01 = { low_corner[0], y_prenotch, low_corner[2] }; + 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 ); +*/ + // ==================================================== + // Force model + // ==================================================== + using model_type = CabanaPD::PMB; + CabanaPD::ForceModel force_model( model_type{}, horizon, K, G0 ); + + // ==================================================== + // Particle generation + // ==================================================== + CabanaPD::Particles particles( memory_space{}, model_type{} ); + + // Note that individual inputs can be passed instead (see other examples). + particles.domain( inputs ); + particles.create( exec_space{} ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + + double dz = particles.dx[2]; + CabanaPD::Region square_pressure( + 0.5 * low_corner[0], 0.5 * high_corner[0], 0.5 * low_corner[1], + 0.5 * high_corner[1], high_corner[2] - dz, high_corner[2] + dz ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles.sliceDensity(); + auto x = particles.sliceReferencePosition(); + auto v = particles.sliceVelocity(); + auto f = particles.sliceForce(); + auto nofail = particles.sliceNoFail(); + + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + + // Density + rho( pid ) = rho0; +/* + // No-fail zone + if ( x( pid, 1 ) <= plane1.low[1] + horizon + 1e-10 || + x( pid, 1 ) >= plane2.high[1] - horizon - 1e-10 ) + nofail( pid ) = 1; +*/ + }; + particles.update( exec_space{}, init_functor ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + + // ==================================================== + // Boundary conditions + // ==================================================== + // Create BC last to ensure ghost particles are included. + double sigma0 = inputs["traction"]; + double b0 = -sigma0 / dz; + auto indent_force = 0; + f = solver.particles.sliceForce(); + x = solver.particles.sliceReferencePosition(); + // Create a symmetric force BC in the z-direction. + + auto bc_op = KOKKOS_LAMBDA( const int pid, const double ) + { + double xsq = x(pid,0) * x(pid,0); + double ysq = x(pid,1) * x(pid,1); + double rsq = xsq + ysq; + + if ( xsq + ysq < 9e-6 ) + { + f( pid, 2 ) += b0; + } + + + }; + auto bc = createBoundaryCondition( bc_op, exec_space{}, solver.particles, + true, square_pressure ); + + // ==================================================== + // Simulation run + // ==================================================== + //solver.init( bc, prenotch ); + solver.init( bc ); + solver.run( bc ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + flatIndenterExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/mechanics/indentation_test.cpp b/examples/mechanics/indentation_test.cpp new file mode 100644 index 00000000..29af5745 --- /dev/null +++ b/examples/mechanics/indentation_test.cpp @@ -0,0 +1,164 @@ +/**************************************************************************** + * Copyright (c) 2022 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 "mpi.h" + +#include + +#include + +// Simulate crack branching from an pre-crack. +void crackBranchingExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double horizon = inputs["horizon"]; + horizon += 1e-10; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + + // ==================================================== + // Pre-notch + // ==================================================== + /* + double height = inputs["system_size"][0]; + double thickness = inputs["system_size"][2]; + double L_prenotch = height / 2.0; + double y_prenotch = 0.0; + Kokkos::Array p01 = { low_corner[0], y_prenotch, low_corner[2] }; + 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 ); + */ + + // ==================================================== + // Force model + // ==================================================== + using model_type = CabanaPD::PMB; + CabanaPD::ForceModel force_model( model_type{}, horizon, K, G0 ); + + // ==================================================== + // Particle generation + // ==================================================== + CabanaPD::Particles particles( memory_space{}, model_type{} ); + + // Note that individual inputs can be passed instead (see other examples). + particles.domain( inputs ); + particles.create( exec_space{} ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + // double dy = particles.dx[1]; + // CabanaPD::Region plane1( + // low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, + // low_corner[2], high_corner[2] ); + // CabanaPD::Region plane2( + // low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + + // dy, low_corner[2], high_corner[2] ); + + double dz = particles.dx[2]; + CabanaPD::Region square_pressure( + 0.5 * low_corner[0], 0.5 * high_corner[0], 0.5 * low_corner[1], + 0.5 * high_corner[1], high_corner[2] - dz, high_corner[2] + dz ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles.sliceDensity(); + auto x = particles.sliceReferencePosition(); + auto v = particles.sliceVelocity(); + auto f = particles.sliceForce(); + auto nofail = particles.sliceNoFail(); + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + // No-fail zone + // if ( x( pid, 1 ) <= plane1.low[1] + horizon + 1e-10 || + // x( pid, 1 ) >= plane2.high[1] - horizon - 1e-10 ) + // nofail( pid ) = 1; + }; + particles.update( exec_space{}, init_functor ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + + // ==================================================== + // Boundary conditions + // ==================================================== + // Create BC last to ensure ghost particles are included. + double sigma0 = inputs["traction"]; + // double b0 = sigma0 / dy; + double b0 = -sigma0 / dz; + f = solver.particles.sliceForce(); + x = solver.particles.sliceReferencePosition(); + // Create a symmetric force BC in the y-direction. + auto bc_op = KOKKOS_LAMBDA( const int pid, const double ) + { + // auto ypos = x( pid, 1 ); + // auto sign = std::abs( ypos ) / ypos; + // f( pid, 1 ) += b0 * sign; + f( pid, 2 ) += b0; + }; + // auto bc = createBoundaryCondition( bc_op, exec_space{}, solver.particles, + // true, plane1, plane2 ); + + auto bc = createBoundaryCondition( bc_op, exec_space{}, solver.particles, + true, square_pressure ); + + // ==================================================== + // Simulation run + // ==================================================== + // solver.init( bc, prenotch ); + solver.init( bc ); + solver.run( bc ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + crackBranchingExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/mechanics/inputs/berkovich_indenter.json b/examples/mechanics/inputs/berkovich_indenter.json new file mode 100644 index 00000000..82bfbf92 --- /dev/null +++ b/examples/mechanics/inputs/berkovich_indenter.json @@ -0,0 +1,17 @@ +{ + "num_cells" : {"value": [200, 200, 1]}, + "system_size" : {"value": [0.005, 0.005, 2.5e-05], "unit": "m"}, + "density" : {"value": 2440, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 1e+7, "unit": "Pa"}, + "fracture_energy" : {"value": 3.8, "unit": "J/m^2"}, + "horizon" : {"value": 0.75e-4, "unit": "m"}, + "indenter_angle" : {"value": 65.35, "unit": "degrees"}, + "indenter_velocity" : {"value": 2.0, "unit": "m/s"}, + "contact_horizon_factor" : {"value": 0.9}, + "contact_horizon_extend_factor" : {"value": 0.01}, + "final_time" : {"value": 5e-4, "unit": "s"}, + "timestep" : {"value": 0.35e-8, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": false} +} diff --git a/examples/mechanics/inputs/conical_indenter.json b/examples/mechanics/inputs/conical_indenter.json new file mode 100644 index 00000000..3c50e78c --- /dev/null +++ b/examples/mechanics/inputs/conical_indenter.json @@ -0,0 +1,17 @@ +{ + "num_cells" : {"value": [200, 200, 1]}, + "system_size" : {"value": [0.005, 0.005, 2.5e-05], "unit": "m"}, + "density" : {"value": 2440, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 1e+7, "unit": "Pa"}, + "fracture_energy" : {"value": 3.8, "unit": "J/m^2"}, + "horizon" : {"value": 0.75e-4, "unit": "m"}, + "indenter_angle" : {"value": 45, "unit": "degrees"}, + "indenter_velocity" : {"value": 2.0, "unit": "m/s"}, + "contact_horizon_factor" : {"value": 0.9}, + "contact_horizon_extend_factor" : {"value": 0.01}, + "final_time" : {"value": 5e-4, "unit": "s"}, + "timestep" : {"value": 0.35e-8, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": false} +} diff --git a/examples/mechanics/inputs/flat_indenter.json b/examples/mechanics/inputs/flat_indenter.json new file mode 100644 index 00000000..3b9b3fa5 --- /dev/null +++ b/examples/mechanics/inputs/flat_indenter.json @@ -0,0 +1,14 @@ +{ + "num_cells" : {"value": [400, 160, 8]}, + "system_size" : {"value": [0.1, 0.04, 0.002], "unit": "m"}, + "density" : {"value": 2440, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 72e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 3.8, "unit": "J/m^2"}, + "horizon" : {"value": 0.001, "unit": "m"}, + "traction" : {"value": 2e6, "unit": "Pa"}, + "final_time" : {"value": 43e-6, "unit": "s"}, + "timestep" : {"value": 4.5e-8, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 5}, + "output_reference" : {"value": true} +} diff --git a/examples/mechanics/inputs/indentation_test.json b/examples/mechanics/inputs/indentation_test.json new file mode 100644 index 00000000..f20aa5f9 --- /dev/null +++ b/examples/mechanics/inputs/indentation_test.json @@ -0,0 +1,19 @@ +{ + "num_cells_OLD" : {"value": [400, 400, 80]}, + "num_cells" : {"value": [200, 200, 40]}, + "system_size_OLD" : {"value": [0.1, 0.04, 0.002], "unit": "m"}, + "system_size" : {"value": [0.1, 0.1, 0.02], "unit": "m"}, + "density" : {"value": 2440, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 72e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 3.8, "unit": "J/m^2"}, + "horizon_ref" : {"value": 0.001, "unit": "m"}, + "horizon" : {"value": 0.0015, "unit": "m"}, + "traction" : {"value": 2e6, "unit": "Pa"}, + "final_time" : {"value": 43e-6, "unit": "s"}, + "timestep_OLD" : {"value": 4.5e-8, "unit": "s"}, + "timestep" : {"value": 6.75e-8, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 5}, + "output_reference_OLD" : {"value": true}, + "output_reference" : {"value": false} +} diff --git a/examples/mechanics/inputs/spherical_indenter.json b/examples/mechanics/inputs/spherical_indenter.json new file mode 100644 index 00000000..3da4718b --- /dev/null +++ b/examples/mechanics/inputs/spherical_indenter.json @@ -0,0 +1,17 @@ +{ + "num_cells" : {"value": [200, 200, 1]}, + "system_size" : {"value": [0.005, 0.005, 2.5e-05], "unit": "m"}, + "density" : {"value": 2440, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 1e+8, "unit": "Pa"}, + "fracture_energy" : {"value": 3.8, "unit": "J/m^2"}, + "horizon" : {"value": 0.75e-4, "unit": "m"}, + "indenter_radius" : {"value": 3e-3, "unit": "m"}, + "indenter_velocity" : {"value": 2.0, "unit": "m/s"}, + "contact_horizon_factor" : {"value": 0.9}, + "contact_horizon_extend_factor" : {"value": 0.01}, + "final_time" : {"value": 5e-4, "unit": "s"}, + "timestep" : {"value": 0.35e-8, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": false} +} diff --git a/examples/mechanics/spherical_indenter.cpp b/examples/mechanics/spherical_indenter.cpp new file mode 100644 index 00000000..8fe9a03b --- /dev/null +++ b/examples/mechanics/spherical_indenter.cpp @@ -0,0 +1,247 @@ +/**************************************************************************** + * Copyright (c) 2022 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 "mpi.h" + +#include + +#include + +// Simulate a spherical indenter on a plate. +void sphericalIndenterExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double horizon = inputs["horizon"]; + horizon += 1e-10; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + + // ==================================================== + // Force model + // ==================================================== + using model_type = CabanaPD::PMB; + // using contact_type = CabanaPD::NormalRepulsionModel; + CabanaPD::ForceModel force_model( model_type{}, horizon, K, G0 ); + + // ==================================================== + // Particle generation + // ==================================================== + CabanaPD::Particles particles( memory_space{}, model_type{} ); + // CabanaPD::Particles particles( memory_space{}, contact_type{} ); + + // Note that individual inputs can be passed instead (see other examples). + particles.domain( inputs ); + particles.create( exec_space{} ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + double dx = particles.dx[0]; + double dy = particles.dx[1]; + double dz = particles.dx[2]; + double R = inputs["indenter_radius"]; + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + + // Define cylindrical region to apply pressure. + CabanaPD::Region pressure_region( + 0.0, R, high_corner[2] - dz, high_corner[2] + dz, x_center, y_center ); + + // Define boundary condition regions to constrain displacements. + CabanaPD::Region low_x_plane( + low_corner[0] - dx, low_corner[0] + dx, low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + CabanaPD::Region high_x_plane( + high_corner[0] - dx, high_corner[0] + dx, low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + CabanaPD::Region low_y_plane( + low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, + low_corner[2], high_corner[2] ); + CabanaPD::Region high_y_plane( + low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy, + low_corner[2], high_corner[2] ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles.sliceDensity(); + auto nofail = particles.sliceNoFail(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + nofail( pid ) = 1; + }; + particles.update( exec_space{}, init_functor ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + + /* + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx / 2.0; + r_extend *= dx; + + contact_type contact_model( horizon, r_c, r_extend, K ); + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + */ + + // ==================================================== + // Boundary conditions + // ==================================================== + // Create BC last to ensure ghost particles are included. + auto x = solver.particles.sliceReferencePosition(); + auto u = solver.particles.sliceDisplacement(); + double v0 = inputs["indenter_velocity"]; + + // z-coordinate of top layer of particles + double z_top = high_corner[2] - 0.5 * dz; + // Initial z-coordinate of the indenter center + double z0_indenter = z_top + R; + + auto disp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + // z-coordinate of the indenter center + double z_indenter = z0_indenter - v0 * t; + + // Check if indenter touches the plate (top layer of particles) + if ( z_indenter - R < z_top ) + { + double r_indenter_sq; + + if ( z_indenter > z_top ) + { + // Current radius of indenter squared + r_indenter_sq = + R * R - ( z_indenter - z_top ) * ( z_indenter - z_top ); + } + else + { + // Current radius of indenter squared + r_indenter_sq = R * R; + } + + // Distance squared of particle to center on XY plane (assumes + // indenter is centered on the XY plane) + double r_sq = + ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + + // Indenter displacement + if ( r_sq <= r_indenter_sq ) + u( pid, 2 ) = + ( z_indenter - std::sqrt( R * R - r_sq ) ) - x( pid, 2 ); + } + + // Edge constraints + if ( low_x_plane.inside( x, pid ) || high_x_plane.inside( x, pid ) ) + { + u( pid, 0 ) = 0; + u( pid, 1 ) = 0; + u( pid, 2 ) = 0; + } + + if ( low_y_plane.inside( x, pid ) || high_y_plane.inside( x, pid ) ) + { + u( pid, 0 ) = 0; + u( pid, 1 ) = 0; + u( pid, 2 ) = 0; + } + }; + + auto bc = createBoundaryCondition( + disp_func, exec_space{}, solver.particles, true, pressure_region, + low_x_plane, high_x_plane, low_y_plane, high_y_plane ); + + //======================================== + // OUTPUTS + //======================================== + auto f = solver.particles.sliceForce(); + + // Output force in x-direction. + auto force_func_x = KOKKOS_LAMBDA( const int p ) + { + return f( p, 0 ) * dx * dy * dz; + }; + auto output_fx = CabanaPD::createOutputTimeSeries( + "output_force_x.txt", inputs, exec_space{}, solver.particles, + force_func_x, pressure_region ); + + // Output force in y-direction. + auto force_func_y = KOKKOS_LAMBDA( const int p ) + { + return f( p, 1 ) * dx * dy * dz; + }; + auto output_fy = CabanaPD::createOutputTimeSeries( + "output_force_y.txt", inputs, exec_space{}, solver.particles, + force_func_y, pressure_region ); + + // Output force in z-direction. + auto force_func_z = KOKKOS_LAMBDA( const int p ) + { + return f( p, 2 ) * dx * dy * dz; + }; + auto output_fz = CabanaPD::createOutputTimeSeries( + "output_force_z.txt", inputs, exec_space{}, solver.particles, + force_func_z, pressure_region ); + + // ==================================================== + // Simulation run + // ==================================================== + solver.init( bc ); + + solver.updateRegion( output_fx, output_fy, output_fz ); + + solver.run( bc, output_fx, output_fy, output_fz ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + sphericalIndenterExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +}