Add temperature field

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(ThermalDeformation thermal_deformation.cpp)
+target_link_libraries(ThermalDeformation LINK_PUBLIC CabanaPD)
+
+install(TARGETS ElasticWave KalthoffWinkler CrackBranching ThermalDeformation DESTINATION ${CMAKE_INSTALL_BINDIR})

examples/crack_branching.cpp

@@ -63,9 +63,12 @@ int main( int argc, char* argv[] )
         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 );
         using model_type =
             CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
-        model_type force_model( delta, K, G0 );
+        model_type force_model( delta, K, 0.0, G0 );
 
         // Create particles from mesh.
         // Does not set displacements, velocities, etc.

examples/inputs/kalthoff_winkler_m3.json

@@ -0,0 +1,13 @@
+{
+    "num_cells"         : {"value": [151, 301, 14]},
+    "system_size"       : {"value": [0.1, 0.2, 0.009], "unit": "m"},
+    "density"           : {"value": 8000,  "unit": "kg/m^3"},
+    "elastic_modulus"   : {"value": 191e+9, "unit": "Pa"},
+    "fracture_energy"   : {"value": 42408,   "unit": "J/m^2"},
+    "horizon"           : {"value": 0.002, "unit": "m"}, 
+    "impactor_velocity" : {"value": 16, "unit": "m/sec"}, 
+    "final_time"        : {"value": 70e-6, "unit": "s"},
+    "timestep"          : {"value": 0.133e-6,  "unit": "s"},
+    "output_frequency"  : {"value": 500},
+    "output_reference"  : {"value": true}
+}

examples/inputs/kalthoff_winkler_m4.json

@@ -0,0 +1,13 @@
+{
+    "num_cells"         : {"value": [200, 400, 18]},
+    "system_size"       : {"value": [0.1, 0.2, 0.009], "unit": "m"},
+    "density"           : {"value": 8000,  "unit": "kg/m^3"},
+    "elastic_modulus"   : {"value": 191e+9, "unit": "Pa"},
+    "fracture_energy"   : {"value": 42408,   "unit": "J/m^2"},
+    "horizon"           : {"value": 0.002, "unit": "m"}, 
+    "impactor_velocity" : {"value": 16, "unit": "m/sec"}, 
+    "final_time"        : {"value": 70e-6, "unit": "s"},
+    "timestep"          : {"value": 0.133e-6,  "unit": "s"},
+    "output_frequency"  : {"value": 500},
+    "output_reference"  : {"value": true}
+}

examples/inputs/kalthoff_winkler_m6.json

@@ -0,0 +1,13 @@
+{
+    "num_cells"         : {"value": [300, 600, 27]},
+    "system_size"       : {"value": [0.1, 0.2, 0.009], "unit": "m"},
+    "density"           : {"value": 8000,  "unit": "kg/m^3"},
+    "elastic_modulus"   : {"value": 191e+9, "unit": "Pa"},
+    "fracture_energy"   : {"value": 42408,   "unit": "J/m^2"},
+    "horizon"           : {"value": 0.002, "unit": "m"}, 
+    "impactor_velocity" : {"value": 16, "unit": "m/sec"}, 
+    "final_time"        : {"value": 70e-6, "unit": "s"},
+    "timestep"          : {"value": 0.133e-6,  "unit": "s"},
+    "output_frequency"  : {"value": 500},
+    "output_reference"  : {"value": true}
+}

examples/inputs/thermal.json

@@ -0,0 +1,14 @@
+{
+    "num_cells"       : {"value": [100, 100, 43]},
+    "system_size"     : {"value": [0.028, 0.028, 0.012], "unit": "m"},
+    "radius"          : {"value": 0.0085, "unit": "m"},
+    "density"         : {"value": 19.25e+3,  "unit": ""},
+    "elastic_modulus" : {"value": 340e+9, "unit": ""},
+    "shear_modulus"   : {"value": 0.5, "unit": ""},
+    "thermal_coeff"   : {"value": 4.2E-6, "unit": ""},
+    "horizon"         : {"value": 8.4E-4, "unit": ""},  
+    "final_time"      : {"value": 0.0093, "unit": ""},
+    "timestep"        : {"value": 7.5E-7,  "unit": ""},
+    "output_frequency": {"value": 10},
+    "output_reference": {"value": true}
+}

examples/kalthoff_winkler.cpp

@@ -36,7 +36,7 @@ int main( int argc, char* argv[] )
         double K = E / ( 3.0 * ( 1.0 - 2.0 * nu ) );
         double rho0 = inputs["density"];
         double G0 = inputs["fracture_energy"];
-        // double G = E / ( 2.0 * ( 1.0 + nu ) ); // Only for LPS.
+        double G = E / ( 2.0 * ( 1.0 + nu ) ); // Only for LPS.
 
         // PD horizon
         double delta = inputs["horizon"];
@@ -73,11 +73,8 @@ int main( int argc, char* argv[] )
 
         // Choose force model type.
         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 );
+            CabanaPD::ForceModel<CabanaPD::LPS, CabanaPD::Fracture>;
+        model_type force_model( delta, K, G, G0 );
 
         // Create particles from mesh.
         // Does not set displacements, velocities, etc.

examples/thermal_deformation.cpp

@@ -0,0 +1,117 @@
+/****************************************************************************
+ * 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>
+
+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;
+
+        CabanaPD::Inputs inputs( argv[1] );
+        double E = inputs["elastic_modulus"];
+        double rho0 = inputs["density"];
+        double nu = 0.25;                      // unitless
+        double K = E / ( 3 * ( 1 - 2 * nu ) ); // [Pa]
+        double delta = inputs["horizon"];
+
+        double alpha = inputs["thermal_coeff"]; // [1/oC]
+        // Reference temperature
+        // double temp0 = 0.0;
+
+        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.
+
+        // Choose force model type.
+        using model_type =
+            CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Elastic>;
+        // model_type force_model( delta, K );
+        model_type force_model( delta, K, alpha );
+        // using model_type =
+        //     CabanaPD::ForceModel<CabanaPD::LinearLPS, CabanaPD::Elastic>;
+        // model_type force_model( delta, K, G );
+
+        // Create particles from mesh.
+        // Does not set displacements, velocities, etc.
+        auto particles = std::make_shared<
+            CabanaPD::Particles<memory_space, typename model_type::base_model>>(
+            exec_space(), low_corner, high_corner, num_cells, halo_width );
+        // Do not create particles in the center.
+        double x_center = 0.0;
+        double y_center = -0.0005;
+        double radius = inputs["radius"];
+        auto init_op = KOKKOS_LAMBDA( const int, const double x[3] )
+        {
+            if ( ( ( x[0] - x_center ) * ( x[0] - x_center ) +
+                   ( x[1] - y_center ) * ( x[1] - y_center ) ) <
+                 radius * radius )
+                return false;
+            return true;
+        };
+        particles->createParticles( exec_space(), init_op );
+
+        // Define particle initialization.
+        auto x = particles->sliceReferencePosition();
+        auto u = particles->sliceDisplacement();
+        auto f = particles->sliceForce();
+        auto v = particles->sliceVelocity();
+        auto rho = particles->sliceDensity();
+        // auto temp = particles->sliceTemperature();
+
+        // Domain to apply b.c.
+        CabanaPD::RegionBoundary domain1( low_corner[0], high_corner[0],
+                                          low_corner[1], high_corner[1],
+                                          low_corner[2], high_corner[2] );
+        std::vector<CabanaPD::RegionBoundary> domain = { domain1 };
+
+        auto bc = createBoundaryCondition( CabanaPD::TempBCTag{}, 5000.0,
+                                           exec_space{}, *particles, domain );
+
+        auto init_functor = KOKKOS_LAMBDA( const int pid )
+        {
+            // temp( pid ) = 5000 * x( pid, 1 ) * t_final;
+            rho( pid ) = rho0;
+        };
+        particles->updateParticles( exec_space{}, init_functor );
+
+        auto cabana_pd = CabanaPD::createSolverElastic<memory_space>(
+            inputs, particles, force_model, bc );
+        cabana_pd->init_force();
+        cabana_pd->run();
+
+        double num_cell_x = num_cells[0];
+        auto profile = Kokkos::View<double* [2], memory_space>(
+            Kokkos::ViewAllocateWithoutInitializing( "displacement_profile" ),
+            num_cell_x );
+        int mpi_rank;
+        MPI_Comm_rank( MPI_COMM_WORLD, &mpi_rank );
+        Kokkos::View<int*, memory_space> count( "c", 1 );
+    }
+
+    MPI_Finalize();
+}

src/CabanaPD_Boundary.hpp

@@ -131,6 +131,10 @@ struct ForceUpdateBCTag
 {
 };
 
+struct TempBCTag
+{
+};
+
 struct ZeroBCTag
 {
 };
@@ -155,7 +159,7 @@ struct BoundaryCondition
     }
 
     template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType& )
+    void apply( ExecSpace, ParticleType&, double )
     {
         auto user = _user_functor;
         auto index_space = _index_space._view;
@@ -168,6 +172,42 @@ struct BoundaryCondition
     }
 };
 
+template <class BCIndexSpace>
+struct BoundaryCondition<BCIndexSpace, TempBCTag>
+{
+    double _value;
+    BCIndexSpace _index_space;
+
+    BoundaryCondition( const double value, BCIndexSpace bc_index_space )
+        : _value( value )
+        , _index_space( bc_index_space )
+    {
+    }
+
+    template <class ExecSpace, class Particles>
+    void update( ExecSpace exec_space, Particles particles,
+                 RegionBoundary plane )
+    {
+        _index_space.update( exec_space, particles, plane );
+    }
+
+    template <class ExecSpace, class ParticleType>
+    void apply( ExecSpace, ParticleType& particles, double t )
+    {
+        auto temp = particles.sliceTemperature();
+        auto x = particles.sliceReferencePosition();
+        auto value = _value;
+        auto index_space = _index_space._view;
+        Kokkos::RangePolicy<ExecSpace> policy( 0, index_space.size() );
+        Kokkos::parallel_for(
+            "CabanaPD::BC::apply", policy, KOKKOS_LAMBDA( const int b ) {
+                auto pid = index_space( b );
+                // This is specifically for the thermal deformation problem
+                temp( pid ) = value * x( pid, 1 ) * t;
+            } );
+    }
+};
+
 template <class BCIndexSpace>
 struct BoundaryCondition<BCIndexSpace, ZeroBCTag>
 {
@@ -177,7 +217,7 @@ struct BoundaryCondition<BCIndexSpace, ZeroBCTag>
     }
 
     template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType )
+    void apply( ExecSpace, ParticleType, double )
     {
     }
 };
@@ -202,7 +242,7 @@ struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
     }
 
     template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType& particles )
+    void apply( ExecSpace, ParticleType& particles, double )
     {
         auto f = particles.sliceForce();
         auto index_space = _index_space._view;
@@ -237,7 +277,7 @@ struct BoundaryCondition<BCIndexSpace, ForceUpdateBCTag>
     }
 
     template <class ExecSpace, class ParticleType>
-    void apply( ExecSpace, ParticleType& particles )
+    void apply( ExecSpace, ParticleType& particles, double )
     {
         auto f = particles.sliceForce();
         auto index_space = _index_space._view;
@@ -257,7 +297,7 @@ template <class BoundaryType, class BCTag, class ExecSpace, class Particles>
 auto createBoundaryCondition( BCTag, const double value, ExecSpace exec_space,
                               Particles particles,
                               std::vector<BoundaryType> planes,
-                              const double initial_guess = 0.5 )
+                              const double initial_guess = 1.1 )
 {
     using memory_space = typename Particles::memory_space;
     using bc_index_type = BoundaryIndexSpace<memory_space, BoundaryType>;