Add thermal deformation example

README.md

@@ -107,8 +107,14 @@ ctest
 
 ## Examples
 
-Once built and installed, CabanaPD examples can be run. Timing and energy
-information is output to file and particle output is written to files (if enabled within Cabana) that can be visualized with Paraview and similar applications. 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:
+Once built and installed, CabanaPD `examples/` can be run. Timing and energy
+information is output to file and particle output is written to files (if enabled within Cabana) that can be visualized with Paraview and similar applications. 
+New examples can be created by using any of the current cases as a template. All inputs are specified in the example JSON files within the relevant `inputs/` subdirectory.
+
+### Mechanics
+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
@@ -129,7 +135,15 @@ run with:
 ./CabanaPD/build/install/bin/CrackBranching CabanaPD/examples/inputs/crack_branching.json
 ```
 
-New examples can be created by using any of the current cases as a template. All inputs are currently specified in the example source files themselves.
+### Thermomechanics
+Examples which demonstrate temperature-dependent mechanics and fracture are with `examples/thermomechanics`.
+
+The first example demonstrates a thermoelastic problem without fracture with a homogeneous plate
+under linear thermal loading [4]. The example can be run with:
+
+```
+./CabanaPD/build/install/bin/ThermalDeformation CabanaPD/examples/thermomechanics/thermal_deformation.json
+```
 
 ## References
 
@@ -145,6 +159,8 @@ 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.
+
 ## Contributing
 
 We encourage you to contribute to CabanaPD! Please check the

examples/CMakeLists.txt

@@ -1,10 +1,2 @@
-add_executable(ElasticWave elastic_wave.cpp)
-target_link_libraries(ElasticWave LINK_PUBLIC CabanaPD)
-
-add_executable(KalthoffWinkler kalthoff_winkler.cpp)
-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_subdirectory(mechanics)
+add_subdirectory(thermomechanics)

examples/mechanics/CMakeLists.txt

@@ -0,0 +1,10 @@
+add_executable(ElasticWave elastic_wave.cpp)
+target_link_libraries(ElasticWave LINK_PUBLIC CabanaPD)
+
+add_executable(KalthoffWinkler kalthoff_winkler.cpp)
+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})

examples/crack_branching.cpp → examples/mechanics/crack_branching.cpp

@@ -79,7 +79,8 @@ void crackBranchingExample( const std::string filename )
     // ====================================================
     // Does not set displacements, velocities, etc.
     auto particles = std::make_shared<
-        CabanaPD::Particles<memory_space, typename model_type::base_model>>(
+        CabanaPD::Particles<memory_space, typename model_type::base_model,
+                            typename model_type::thermal_type>>(
         exec_space(), low_corner, high_corner, num_cells, halo_width );
 
     // ====================================================
@@ -109,8 +110,8 @@ void crackBranchingExample( const std::string filename )
         auto sign = std::abs( ypos ) / ypos;
         Cabana::get( p_f, CabanaPD::Field::Force(), 1 ) += b0 * sign;
     };
-    auto bc =
-        createBoundaryCondition( bc_op, exec_space{}, *particles, planes );
+    auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes,
+                                       true );
 
     // ====================================================
     //            Custom particle initialization
@@ -136,9 +137,9 @@ void crackBranchingExample( const std::string filename )
     //                   Simulation run
     // ====================================================
     auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
-        inputs, particles, force_model, bc, prenotch );
-    cabana_pd->init_force();
-    cabana_pd->run();
+        inputs, particles, force_model, prenotch );
+    cabana_pd->init();
+    cabana_pd->run( bc );
 }
 
 // Initialize MPI+Kokkos.

examples/elastic_wave.cpp → examples/mechanics/elastic_wave.cpp

@@ -64,15 +64,10 @@ void elasticWaveExample( const std::string filename )
     // ====================================================
     // Does not set displacements, velocities, etc.
     auto particles = std::make_shared<
-        CabanaPD::Particles<memory_space, typename model_type::base_model>>(
+        CabanaPD::Particles<memory_space, typename model_type::base_model,
+                            typename model_type::thermal_type>>(
         exec_space(), low_corner, high_corner, num_cells, halo_width );
 
-    // ====================================================
-    //                Boundary conditions
-    // ====================================================
-    auto bc = CabanaPD::createBoundaryCondition<memory_space>(
-        CabanaPD::ZeroBCTag{} );
-
     // ====================================================
     //            Custom particle initialization
     // ====================================================
@@ -110,8 +105,8 @@ void elasticWaveExample( const std::string filename )
     //                   Simulation run
     // ====================================================
     auto cabana_pd = CabanaPD::createSolverElastic<memory_space>(
-        inputs, particles, force_model, bc );
-    cabana_pd->init_force();
+        inputs, particles, force_model );
+    cabana_pd->init();
     cabana_pd->run();
 
     // ====================================================

examples/kalthoff_winkler.cpp → examples/mechanics/kalthoff_winkler.cpp

@@ -89,7 +89,8 @@ void kalthoffWinklerExample( const std::string filename )
     // ====================================================
     // Does not set displacements, velocities, etc.
     auto particles = std::make_shared<
-        CabanaPD::Particles<memory_space, typename model_type::base_model>>(
+        CabanaPD::Particles<memory_space, typename model_type::base_model,
+                            typename model_type::thermal_type>>(
         exec_space(), low_corner, high_corner, num_cells, halo_width );
 
     // ====================================================
@@ -128,9 +129,9 @@ void kalthoffWinklerExample( const std::string filename )
     //                   Simulation run
     // ====================================================
     auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
-        inputs, particles, force_model, bc, prenotch );
-    cabana_pd->init_force();
-    cabana_pd->run();
+        inputs, particles, force_model, prenotch );
+    cabana_pd->init();
+    cabana_pd->run( bc );
 }
 
 // Initialize MPI+Kokkos.

examples/thermomechanics/CMakeLists.txt

@@ -0,0 +1,3 @@
+add_executable(ThermalDeformation thermal_deformation.cpp)
+target_link_libraries(ThermalDeformation LINK_PUBLIC CabanaPD)
+install(TARGETS ThermalDeformation DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/thermomechanics/inputs/thermal_deformation.json

@@ -0,0 +1,13 @@
+{
+    "num_cells"             : {"value": [100, 30, 3]},
+    "system_size"           : {"value": [1.0, 0.3, 0.03], "unit": "m"},
+    "density"               : {"value": 3980,  "unit": "kg/m^3"},
+    "elastic_modulus"       : {"value": 370e+9, "unit": "Pa"},
+    "thermal_coefficient"   : {"value": 7.5E-6, "unit": "oC^{-1}"},
+    "reference_temperature" : {"value": 20.0, "unit": "oC"},
+    "horizon"               : {"value": 0.03, "unit": "m"},  
+    "final_time"            : {"value": 0.01, "unit": "s"},
+    "timestep"              : {"value": 7.5E-7,  "unit": "s"},
+    "output_frequency"      : {"value": 100},
+    "output_reference"      : {"value": true}
+}

examples/thermomechanics/thermal_deformation.cpp

@@ -0,0 +1,136 @@
+/****************************************************************************
+ * 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 thermally-induced deformation in a rectangular plate.
+void thermalDeformationExample( 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 and problem parameters
+    // ====================================================
+    // Material parameters
+    double rho0 = inputs["density"];
+    double E = inputs["elastic_modulus"];
+    double nu = 0.25;
+    double K = E / ( 3 * ( 1 - 2 * nu ) );
+    double delta = inputs["horizon"];
+    double alpha = inputs["thermal_coefficient"];
+
+    // Problem parameters
+    double temp0 = inputs["reference_temperature"];
+
+    // ====================================================
+    //                  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 type
+    // ====================================================
+    using model_type = CabanaPD::PMB;
+    using thermal_type = CabanaPD::TemperatureDependent;
+
+    // ====================================================
+    //                 Particle generation
+    // ====================================================
+    // Does not set displacements, velocities, etc.
+    auto particles = std::make_shared<
+        CabanaPD::Particles<memory_space, model_type, thermal_type>>(
+        exec_space(), low_corner, high_corner, num_cells, halo_width );
+
+    // ====================================================
+    //            Custom particle initialization
+    // ====================================================
+    auto rho = particles->sliceDensity();
+    auto init_functor = KOKKOS_LAMBDA( const int pid ) { rho( pid ) = rho0; };
+    particles->updateParticles( exec_space{}, init_functor );
+
+    // ====================================================
+    //                    Force model
+    // ====================================================
+    auto force_model =
+        CabanaPD::createForceModel<model_type, CabanaPD::Elastic, thermal_type>(
+            *particles, delta, K, alpha, temp0 );
+
+    // ====================================================
+    //                   Create solver
+    // ====================================================
+    auto cabana_pd = CabanaPD::createSolverElastic<memory_space>(
+        inputs, particles, force_model );
+
+    // ====================================================
+    //                   Imposed field
+    // ====================================================
+    auto x = particles->sliceReferencePosition();
+    auto temp = particles->sliceTemperature();
+    const double low_corner_y = low_corner[1];
+    // This is purposely delayed until after solver init so that ghosted
+    // particles are correctly taken into account for lambda capture here.
+    auto temp_func = KOKKOS_LAMBDA( const int pid, const double t )
+    {
+        temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t;
+    };
+    auto body_term = CabanaPD::createBodyTerm( temp_func, false );
+
+    // ====================================================
+    //                   Simulation run
+    // ====================================================
+    cabana_pd->init( body_term );
+    cabana_pd->run( body_term );
+
+    // ====================================================
+    //                      Outputs
+    // ====================================================
+    // Output displacement along the x-axis
+    createDisplacementProfile( MPI_COMM_WORLD, num_cells[0], 0,
+                               "xdisplacement_profile.txt", *particles );
+
+    // Output displacement along the y-axis
+    createDisplacementProfile( MPI_COMM_WORLD, num_cells[1], 1,
+                               "ydisplacement_profile.txt", *particles );
+}
+
+// Initialize MPI+Kokkos.
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    thermalDeformationExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}

src/CabanaPD_BodyTerm.hpp

@@ -25,11 +25,13 @@ template <class UserFunctor>
 struct BodyTerm
 {
     UserFunctor _user_functor;
+    bool _force_update;
 
     Timer _timer;
 
-    BodyTerm( UserFunctor user )
+    BodyTerm( UserFunctor user, const bool force )
         : _user_functor( user )
+        , _force_update( force )
     {
     }
 
@@ -47,14 +49,16 @@ struct BodyTerm
         _timer.stop();
     }
 
+    auto forceUpdate() { return _force_update; }
+
     auto time() { return _timer.time(); };
     auto timeInit() { return 0.0; };
 };
 
 template <class UserFunctor>
-auto createBodyTerm( UserFunctor user_functor )
+auto createBodyTerm( UserFunctor user_functor, const bool force_update )
 {
-    return BodyTerm<UserFunctor>( user_functor );
+    return BodyTerm<UserFunctor>( user_functor, force_update );
 }
 
 } // namespace CabanaPD