Fragmenting cylinder example

README.md

@@ -117,34 +117,39 @@ Examples which only include mechanics and fracture are with `examples/mechanics`
 The first example is an elastic wave propagating through a cube from an initial Gaussian radial displacement profile from [1]. Assuming the build paths above, the example can be run with:
 
 ```
-./CabanaPD/build/install/bin/ElasticWave CabanaPD/examples/inputs/elastic_wave.json
+./CabanaPD/build/install/bin/ElasticWave CabanaPD/examples/mechanics/inputs/elastic_wave.json
 ```
 
 The second example is the Kalthoff-Winkler experiment [2], where an impactor
 causes crack propagation at an angle from two pre-notches on a steel plate. The
 example can be run with:
 
 ```
-./CabanaPD/build/install/bin/KalthoffWinkler CabanaPD/examples/inputs/kalthoff_winkler.json
+./CabanaPD/build/install/bin/KalthoffWinkler CabanaPD/examples/mechanics/inputs/kalthoff_winkler.json
 ```
 
-The third example is crack branching in a pre-notched soda-lime glass plate due to traction loading [3]. The example can be
-run with:
+The third example is crack branching in a pre-notched soda-lime glass plate due to traction loading [3]. The example can be run with:
 
 ```
-./CabanaPD/build/install/bin/CrackBranching CabanaPD/examples/inputs/crack_branching.json
+./CabanaPD/build/install/bin/CrackBranching CabanaPD/examples/mechanics/inputs/crack_branching.json
+```
+
+The fourth example is a fragmenting cylinder due to internal pressure [4]. The example can be run with:
+
+```
+./CabanaPD/build/install/bin/FragmentingCylinder CabanaPD/examples/mechanics/inputs/fragmenting_cylinder.json
 ```
 
 ### Thermomechanics
 Examples which demonstrate temperature-dependent mechanics and fracture are with `examples/thermomechanics`.
 
-The first example is thermoelastic deformation in a homogeneous plate due to linear thermal loading [4]. The example can be run with:
+The first example is thermoelastic deformation in a homogeneous plate due to linear thermal loading [5]. The example can be run with:
 
 ```
 ./CabanaPD/build/install/bin/ThermalDeformation CabanaPD/examples/thermomechanics/thermal_deformation.json
 ```
 
-The second example is crack initiation and propagation in an alumina ceramic plate due to a thermal shock caused by water quenching [5]. The example can be run with:
+The second example is crack initiation and propagation in an alumina ceramic plate due to a thermal shock caused by water quenching [6]. The example can be run with:
 
 ```
 ./CabanaPD/build/install/bin/ThermalCrack CabanaPD/examples/thermomechanics/thermal_crack.json
@@ -164,9 +169,11 @@ Kunze, and L.W. Meyer, eds., Vol 1, DGM Informationsgesellschaft Verlag (1988)
 
 [3] F. Bobaru and G. Zhang, Why do cracks branch? A peridynamic investigation of dynamic brittle fracture, International Journal of Fracture 196 (2015): 59–98.
 
-[4] D. He, D. Huang, and D. Jiang, Modeling and studies of fracture in functionally graded materials under thermal shock loading using peridynamics, Theoretical and Applied Fracture Mechanics 111 (2021): 102852.
+[4] D.J. Littlewood, M.L. Parks, J.T. Foster, J.A. Mitchell, and P. Diehl, The peridigm meshfree peridynamics code, Journal of Peridynamics and Nonlocal Modeling 6 (2024): 118–148.  
+
+[5] D. He, D. Huang, and D. Jiang, Modeling and studies of fracture in functionally graded materials under thermal shock loading using peridynamics, Theoretical and Applied Fracture Mechanics 111 (2021): 102852.
 
-[5] C.P. Jiang, X.F. Wu, J. Li, F. Song, Y.F. Shao, X.H. Xu, and P. Yan, A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock, Acta Materialia 60 (2012): 4540–4550.
+[6] C.P. Jiang, X.F. Wu, J. Li, F. Song, Y.F. Shao, X.H. Xu, and P. Yan, A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock, Acta Materialia 60 (2012): 4540–4550.
 
 ## Contributing
 

examples/mechanics/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(FragmentingCylinder fragmenting_cylinder.cpp)
+target_link_libraries(FragmentingCylinder LINK_PUBLIC CabanaPD)
+
+install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/mechanics/fragmenting_cylinder.cpp

@@ -0,0 +1,165 @@
+/****************************************************************************
+ * 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 <fstream>
+#include <iostream>
+
+#include "mpi.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <CabanaPD.hpp>
+
+// Simulate an expanding cylinder resulting in fragmentation
+void fragmentingCylinderExample( const std::string filename )
+{
+    // ====================================================
+    //             Use default 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 K = inputs["bulk_modulus"];
+    // double G = inputs["shear_modulus"]; // Only for LPS.
+    double sc = inputs["critical_stretch"];
+    double delta = inputs["horizon"];
+    delta += 1e-10;
+    // For PMB or LPS with influence_type == 1
+    double G0 = 9 * K * delta * ( sc * sc ) / 5;
+    // For LPS with influence_type == 0 (default)
+    // double G0 = 15 * K * delta * ( sc * sc ) / 8;
+
+    // ====================================================
+    //                  Discretization
+    // ====================================================
+    // FIXME: set halo width based on delta
+    std::array<double, 3> low_corner = inputs["low_corner"];
+    std::array<double, 3> high_corner = inputs["high_corner"];
+    std::array<int, 3> num_cells = inputs["num_cells"];
+    int m = std::floor( delta /
+                        ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) );
+    int halo_width = m + 1; // Just to be safe.
+
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    using model_type = CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
+    model_type force_model( delta, K, G0 );
+    // using model_type =
+    //      CabanaPD::ForceModel<CabanaPD::LPS, CabanaPD::Fracture>;
+    // model_type force_model( delta, K, G, G0 );
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    // Does not set displacements, velocities, etc.
+    auto particles = std::make_shared<
+        CabanaPD::Particles<memory_space, typename model_type::base_model,
+                            typename model_type::thermal_type>>(
+        exec_space(), low_corner, high_corner, num_cells, halo_width );
+
+    // ====================================================
+    //            Custom particle initialization
+    // ====================================================
+    double x_center = 0.5 * ( low_corner[0] + high_corner[0] );
+    double y_center = 0.5 * ( low_corner[1] + high_corner[1] );
+    double Rout = inputs["cylinder_outer_radius"];
+    double Rin = inputs["cylinder_inner_radius"];
+
+    // Do not create particles outside given cylindrical region
+    auto init_op = KOKKOS_LAMBDA( const int, const double x[3] )
+    {
+        double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) +
+                     ( x[1] - y_center ) * ( x[1] - y_center );
+        if ( rsq < Rin * Rin || rsq > Rout * Rout )
+            return false;
+        return true;
+    };
+    particles->createParticles( exec_space(), init_op );
+
+    auto rho = particles->sliceDensity();
+    auto x = particles->sliceReferencePosition();
+    auto v = particles->sliceVelocity();
+    auto f = particles->sliceForce();
+
+    double vrmax = inputs["max_radial_velocity"];
+    double vrmin = inputs["min_radial_velocity"];
+    double vzmax = inputs["max_vertical_velocity"];
+    double zmin = low_corner[2];
+
+    auto dx = particles->dx;
+    double height = inputs["system_size"][2];
+    double factor = inputs["grid_perturbation_factor"];
+
+    using pool_type = Kokkos::Random_XorShift64_Pool<exec_space>;
+    using random_type = Kokkos::Random_XorShift64<exec_space>;
+    pool_type pool;
+    int seed = 456854;
+    pool.init( seed, particles->n_local );
+    auto init_functor = KOKKOS_LAMBDA( const int pid )
+    {
+        // Density
+        rho( pid ) = rho0;
+
+        // Perturb particle positions
+        auto gen = pool.get_state();
+        for ( std::size_t d = 0; d < 3; d++ )
+        {
+            auto rand =
+                Kokkos::rand<random_type, double>::draw( gen, 0.0, 1.0 );
+            x( pid, d ) += ( 2.0 * rand - 1.0 ) * factor * dx[d];
+        }
+        pool.free_state( gen );
+
+        // Velocity
+        double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * height ) ) - 1;
+        double vr = vrmax - vrmin * zfactor * zfactor;
+        v( pid, 0 ) =
+            vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
+        v( pid, 1 ) =
+            vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) );
+        v( pid, 2 ) = vzmax * zfactor;
+    };
+    particles->updateParticles( exec_space{}, init_functor );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+
+    // Define empty pre-notch
+    CabanaPD::Prenotch<0> prenotch;
+
+    auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
+        inputs, particles, force_model, prenotch );
+    cabana_pd->init();
+    cabana_pd->run();
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    fragmentingCylinderExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

examples/mechanics/inputs/fragmenting_cylinder.json

@@ -0,0 +1,20 @@
+{
+    "num_cells"                : {"value": [51, 51, 101]},
+    "system_size"              : {"value": [0.05, 0.05, 0.1], "unit": "m"},
+    "grid_perturbation_factor" : {"value": 0.1},
+    "density"                  : {"value": 7800, "unit": "kg/m^3"},
+    "bulk_modulus"             : {"value": 130e+9, "unit": "Pa"},
+    "shear_modulus"            : {"value": 78e+9, "unit": "Pa"},
+    "critical_stretch"         : {"value": 0.02},
+    "horizon"                  : {"value": 0.00417462, "unit": "m"},  
+    "cylinder_outer_radius"    : {"value": 0.025, "unit": "m"},
+    "cylinder_inner_radius"    : {"value": 0.02, "unit": "m"},
+    "max_radial_velocity"      : {"value": 200, "unit": "m/s"}, 
+    "min_radial_velocity"      : {"value": 50, "unit": "m/s"},
+    "max_vertical_velocity"    : {"value": 100, "unit": "m/s"}, 
+    "final_time"               : {"value": 2.5e-4, "unit": "s"},
+    "timestep"                 : {"value": 1.7e-07, "unit": "s"},
+    "timestep_safety_factor"   : {"value": 0.70},
+    "output_frequency"         : {"value": 50},
+    "output_reference"         : {"value": false}
+}

src/CabanaPD_Particles.hpp

@@ -137,7 +137,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     std::shared_ptr<
         Cabana::Grid::LocalGrid<Cabana::Grid::UniformMesh<double, dim>>>
         local_grid;
-    double dx[dim];
+    Kokkos::Array<double, dim> dx;
 
     int halo_width;