From 0b9138928f5af2e6dc4ea39bef5fc9e7998295c5 Mon Sep 17 00:00:00 2001
From: Sam Reeve <6740307+streeve@users.noreply.github.com>
Date: Mon, 6 May 2024 09:22:05 -0400
Subject: [PATCH] Add thermal elastic example

---
 examples/CMakeLists.txt                       |  2 +
 examples/thermal/CMakeLists.txt               |  3 +
 .../thermal/inputs/thermal_deformation.json   | 13 +++
 examples/thermal/thermal_deformation.cpp      | 82 +++++++++++++++++++
 4 files changed, 100 insertions(+)
 create mode 100644 examples/thermal/CMakeLists.txt
 create mode 100644 examples/thermal/inputs/thermal_deformation.json
 create mode 100644 examples/thermal/thermal_deformation.cpp

diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index df1b8e15..61e33aa7 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -8,3 +8,5 @@ add_executable(CrackBranching crack_branching.cpp)
 target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD)
 
 install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR})
+
+add_subdirectory(thermal)
\ No newline at end of file
diff --git a/examples/thermal/CMakeLists.txt b/examples/thermal/CMakeLists.txt
new file mode 100644
index 00000000..84788a90
--- /dev/null
+++ b/examples/thermal/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})
diff --git a/examples/thermal/inputs/thermal_deformation.json b/examples/thermal/inputs/thermal_deformation.json
new file mode 100644
index 00000000..e6bc5205
--- /dev/null
+++ b/examples/thermal/inputs/thermal_deformation.json
@@ -0,0 +1,13 @@
+{
+    "num_cells"             : {"value": [101, 31, 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": 0.0, "unit": "oC"},
+    "horizon"               : {"value": 0.03, "unit": "m"},  
+    "final_time"            : {"value": 0.0093, "unit": "s"},
+    "timestep"              : {"value": 7.5E-7,  "unit": "s"},
+    "output_frequency"      : {"value": 100},
+    "output_reference"      : {"value": true}
+}
diff --git a/examples/thermal/thermal_deformation.cpp b/examples/thermal/thermal_deformation.cpp
new file mode 100644
index 00000000..576a2dc1
--- /dev/null
+++ b/examples/thermal/thermal_deformation.cpp
@@ -0,0 +1,82 @@
+/****************************************************************************
+ * 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>
+
+void thermalDeformationExample( const std::string filename )
+{
+    using exec_space = Kokkos::DefaultExecutionSpace;
+    using memory_space = typename exec_space::memory_space;
+
+    CabanaPD::Inputs inputs( filename );
+    double E = inputs["elastic_modulus"];
+    double rho0 = inputs["density"];
+    double nu = 0.25; // unitless
+    double K = E / ( 3 * ( 1 - 2 * nu ) );
+    double delta = inputs["horizon"];
+
+    double alpha = inputs["thermal_coefficient"];
+    double temp0 = inputs["reference_temperature"];
+
+    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::PMB;
+
+    // Create particles from mesh.
+    auto particles =
+        std::make_shared<CabanaPD::Particles<memory_space, model_type>>(
+            exec_space(), low_corner, high_corner, num_cells, halo_width );
+
+    auto temp = particles->sliceTemperature();
+    auto x = particles->sliceReferencePosition();
+    auto temp_func = KOKKOS_LAMBDA( const int pid, const double t )
+    {
+        temp( pid ) = 5000.0 * ( x( pid, 1 ) - ( -0.15 ) ) * t;
+    };
+    auto body_term = CabanaPD::createBodyTerm( temp_func );
+
+    auto rho = particles->sliceDensity();
+    auto init_functor = KOKKOS_LAMBDA( const int pid ) { rho( pid ) = rho0; };
+    particles->updateParticles( exec_space{}, init_functor );
+
+    auto force_model =
+        CabanaPD::createForceModel<model_type, CabanaPD::Elastic,
+                                   CabanaPD::TemperatureDependent>(
+            *particles, delta, K, alpha, temp0 );
+    auto cabana_pd = CabanaPD::createSolverElastic<memory_space>(
+        inputs, particles, force_model, body_term );
+    cabana_pd->init_force();
+    cabana_pd->run();
+}
+
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    thermalDeformationExample( argv[1] );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+}