Material interface

Dockerfile

@@ -127,7 +127,8 @@ RUN . $SPACK_ROOT/share/spack/setup-env.sh && \
     spack -e aperi-mech add [email protected] ~amesos ~amesos2 ~anasazi ~aztec ~belos ~cuda ~epetra ~epetraext ~ifpack ~ifpack2 ~ml ~muelu ~sacado ~shared +exodus +gtest +hdf5 +mpi +stk +zoltan +zoltan2 cxxstd=17 && \
     spack -e aperi-mech add [email protected] && \
     spack -e aperi-mech add [email protected] && \
-    spack -e aperi-mech add [email protected]
+    spack -e aperi-mech add eigen@master
+# eigen@master is used as [email protected] has a lot of GPU related warnings
 
 # Install Packages, aperi-mech
 RUN . $SPACK_ROOT/share/spack/setup-env.sh && \

Dockerfile_Nvidia_H100_GPU

@@ -170,7 +170,8 @@ RUN . $SPACK_ROOT/share/spack/setup-env.sh && \
     spack -e aperi-mech add [email protected] ~amesos ~amesos2 ~anasazi ~aztec ~belos ~epetra ~epetraext ~ifpack ~ifpack2 ~ml ~muelu ~sacado ~shared +cuda +cuda_rdc +exodus +gtest +hdf5 +stk +zoltan +zoltan2 cxxstd=17 cuda_arch=${CUDA_ARCH} && \
     spack -e aperi-mech add [email protected] && \
     spack -e aperi-mech add [email protected] && \
-    spack -e aperi-mech add [email protected]
+    spack -e aperi-mech add eigen@master
+# eigen@master is used as [email protected] has a lot of GPU related warnings
 
 # Install Packages, aperi-mech
 RUN . $SPACK_ROOT/share/spack/setup-env.sh && \

Dockerfile_Nvidia_T4_GPU

@@ -163,7 +163,8 @@ RUN . $SPACK_ROOT/share/spack/setup-env.sh && \
     spack -e aperi-mech add [email protected] ~amesos ~amesos2 ~anasazi ~aztec ~belos ~epetra ~epetraext ~ifpack ~ifpack2 ~ml ~muelu ~sacado ~shared +cuda +cuda_rdc +exodus +gtest +hdf5 +stk +zoltan +zoltan2 cxxstd=17 cuda_arch=${CUDA_ARCH} && \
     spack -e aperi-mech add [email protected] && \
     spack -e aperi-mech add [email protected] && \
-    spack -e aperi-mech add [email protected]
+    spack -e aperi-mech add eigen@master
+# eigen@master is used as [email protected] has a lot of GPU related warnings
 
 # Install Packages, aperi-mech
 RUN . $SPACK_ROOT/share/spack/setup-env.sh && \

do_configure

@@ -60,7 +60,14 @@ cmake_command+=" -D USE_PROTEGO_MECH:BOOL=${USE_PROTEGO_MECH}"
 cmake_command+=" -D CMAKE_BUILD_TYPE:STRING=\"${BUILD_TYPE}\""
 cmake_command+=" -D TRILINOS_PATH:FILEPATH=$(spack location -i trilinos)"
 cmake_command+=" -D Kokkos_ROOT:FILEPATH=$(spack location -i kokkos)"
-cmake_command+=" -D EIGEN_PATH:FILEPATH=$(spack location -i eigen)"
+# Try setting EIGEN_PATH using spack, fallback to a different path if it fails
+spack uninstall -y eigen || true
+if spack location -i eigen &>/dev/null; then
+	eigen_path=$(spack location -i eigen)
+else
+	eigen_path="/home/aperi-mech_docker/eigen_install"
+fi
+cmake_command+=" -D EIGEN_PATH:FILEPATH=${eigen_path}"
 cmake_command+=" -D GTEST_PATH:FILEPATH=$(spack location -i googletest)"
 cmake_command+=" -D YAML-CPP_PATH:FILEPATH=$(spack location -i yaml-cpp)"
 cmake_command+=" -D OPENMPI_PATH:FILEPATH=$(spack location -i openmpi)"

docker-compose.yml

@@ -15,5 +15,7 @@ services:
       dockerfile: Dockerfile
     volumes:
       - .:/home/aperi-mech_docker/aperi-mech
+      - /home/azureuser/eigen:/home/aperi-mech_docker/eigen
+      - /home/azureuser/eigen_install:/home/aperi-mech_docker/eigen_install
       - /home/azureuser/aperifileshare1:/home/aperi-mech_docker/aperifileshare1
     working_dir: /home/aperi-mech_docker/aperi-mech

docker-compose_nvidia_h100_gpu.yml

@@ -15,6 +15,8 @@ services:
       dockerfile: Dockerfile_Nvidia_H100_GPU
     volumes:
       - .:/home/aperi-mech_docker/aperi-mech
+      - /home/azureuser/eigen:/home/aperi-mech_docker/eigen
+      - /home/azureuser/eigen_install:/home/aperi-mech_docker/eigen_install
     working_dir: /home/aperi-mech_docker/aperi-mech
     runtime: nvidia
     deploy:

docker-compose_nvidia_t4_gpu.yml

@@ -16,6 +16,8 @@ services:
       dockerfile: Dockerfile_Nvidia_T4_GPU
     volumes:
       - .:/home/aperi-mech_docker/aperi-mech
+      - /home/azureuser/eigen:/home/aperi-mech_docker/eigen
+      - /home/azureuser/eigen_install_gpu:/home/aperi-mech_docker/eigen_install
       - /home/azureuser/aperifileshare1:/home/aperi-mech_docker/aperifileshare1
     working_dir: /home/aperi-mech_docker/aperi-mech
     runtime: nvidia

include/ComputeInternalForceFunctors.h

@@ -30,12 +30,22 @@ struct ComputeInternalForceFromIntegrationPointFunctor {
             // Compute displacement gradient
             const Eigen::Matrix3d displacement_gradient = node_displacements.transpose() * b_matrix_and_weight.first;
 
+            // Create a map around the state_old and state_new pointers
+            Eigen::InnerStride<Eigen::Dynamic> state_stride(state_bucket_offset);
+            auto state_old_map = Eigen::Map<const Eigen::VectorXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(state_old, m_stress_functor.NumberOfStateVariables(), state_stride);
+            auto state_new_map = Eigen::Map<Eigen::VectorXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(state_new, m_stress_functor.NumberOfStateVariables(), state_stride);
+
             // Compute the 1st pk stress and internal force of the element.
-            const Eigen::Matrix<double, 3, 3> pk1_stress_transpose_volume = m_stress_functor(displacement_gradient, state_old, state_new, state_bucket_offset).transpose() * b_matrix_and_weight.second;
+            Eigen::Matrix<double, 3, 3> pk1_stress;
+            Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> stride(3, 1);
+            auto pk1_stress_map = Eigen::Map<Eigen::Matrix<double, 3, 3>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(pk1_stress.data(), stride);
+            auto displacement_gradient_map = Eigen::Map<const Eigen::Matrix<double, 3, 3>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(displacement_gradient.data(), stride);
+            double timestep = 1.0;  // TODO(jake): This should be passed in
+            m_stress_functor.GetStress(&displacement_gradient_map, nullptr, &state_old_map, &state_new_map, timestep, pk1_stress_map);
 
             // Compute the internal force
             for (size_t i = 0; i < actual_num_neighbors; ++i) {
-                force.row(i).noalias() -= b_matrix_and_weight.first.row(i) * pk1_stress_transpose_volume;
+                force.row(i).noalias() -= b_matrix_and_weight.first.row(i) * pk1_stress.transpose() * b_matrix_and_weight.second;
             }
         }
     }
@@ -66,12 +76,22 @@ struct ComputeInternalForceFromSmoothingCellFunctor {
     KOKKOS_INLINE_FUNCTION void operator()(const Kokkos::Array<Eigen::Matrix<double, 3, 3>, 1> &gathered_node_data_gradient, Eigen::Matrix<double, MaxNumNodes, 3> &force, const Eigen::Matrix<double, MaxNumNodes, 3> &b_matrix, double volume, size_t actual_num_neighbors, const double *state_old = nullptr, double *state_new = nullptr, size_t state_bucket_offset = 1) const {
         const Eigen::Matrix3d &displacement_gradient = gathered_node_data_gradient[0];
 
+        // Create a map around the state_old and state_new pointers
+        Eigen::InnerStride<Eigen::Dynamic> state_stride(state_bucket_offset);
+        auto state_old_map = Eigen::Map<const Eigen::VectorXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(state_old, m_stress_functor.NumberOfStateVariables(), state_stride);
+        auto state_new_map = Eigen::Map<Eigen::VectorXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(state_new, m_stress_functor.NumberOfStateVariables(), state_stride);
+
         // Compute the stress and internal force of the element.
-        const Eigen::Matrix<double, 3, 3> pk1_stress_transpose_neg_volume = m_stress_functor(displacement_gradient, state_old, state_new, state_bucket_offset).transpose() * -volume;
+        Eigen::Matrix<double, 3, 3> pk1_stress;
+        Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> stride(3, 1);
+        auto pk1_stress_map = Eigen::Map<Eigen::Matrix<double, 3, 3>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(pk1_stress.data(), stride);
+        auto displacement_gradient_map = Eigen::Map<const Eigen::Matrix<double, 3, 3>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(displacement_gradient.data(), stride);
+        double timestep = 1.0;  // TODO(jake): This should be passed in
+        m_stress_functor.GetStress(&displacement_gradient_map, nullptr, &state_old_map, &state_new_map, timestep, pk1_stress_map);
 
         // Compute the internal force
         for (size_t i = 0; i < actual_num_neighbors; ++i) {
-            force.row(i).noalias() = b_matrix.row(i) * pk1_stress_transpose_neg_volume;
+            force.row(i).noalias() = -volume * b_matrix.row(i) * pk1_stress.transpose();
         }
     }
     StressFunctor &m_stress_functor;  ///< Functor for computing the stress of the material
@@ -83,7 +103,24 @@ struct ComputeStressOnSmoothingCellFunctor {
         : m_stress_functor(stress_functor) {}
 
     KOKKOS_INLINE_FUNCTION const Eigen::Matrix<double, 3, 3> operator()(const Kokkos::Array<Eigen::Matrix<double, 3, 3>, 1> &gathered_node_data_gradient, const double *state_old = nullptr, double *state_new = nullptr, size_t state_bucket_offset = 1) const {
-        return m_stress_functor(gathered_node_data_gradient[0], state_old, state_new, state_bucket_offset);  // Compute the stress of the element
+        // PK1 stress
+        Eigen::Matrix<double, 3, 3> pk1_stress;
+
+        // Create a map around the state_old and state_new pointers
+        Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic> stride(3, 1);
+        auto pk1_stress_map = Eigen::Map<Eigen::Matrix<double, 3, 3>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(pk1_stress.data(), stride);
+        auto displacement_gradient_map = Eigen::Map<const Eigen::Matrix<double, 3, 3>, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>>(gathered_node_data_gradient[0].data(), stride);
+
+        // Create a map around the state_old and state_new pointers
+        Eigen::InnerStride<Eigen::Dynamic> state_stride(state_bucket_offset);
+        auto state_old_map = Eigen::Map<const Eigen::VectorXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(state_old, m_stress_functor.NumberOfStateVariables(), state_stride);
+        auto state_new_map = Eigen::Map<Eigen::VectorXd, 0, Eigen::InnerStride<Eigen::Dynamic>>(state_new, m_stress_functor.NumberOfStateVariables(), state_stride);
+
+        // Compute the stress of the element.
+        double timestep = 1.0;  // TODO(jake): This should be passed in
+        m_stress_functor.GetStress(&displacement_gradient_map, nullptr, &state_old_map, &state_new_map, timestep, pk1_stress_map);
+
+        return pk1_stress;
     }
     StressFunctor &m_stress_functor;  ///< Functor for computing the stress of the material
 };

include/ElementForceProcessor.h

@@ -110,7 +110,7 @@ class ElementForceProcessor {
                 size_t num_nodes = nodes.size();
 
                 // Set up the field data to gather
-                Eigen::Matrix<double, NumNodes, 3> coordinate_field_data;
+                Eigen::Matrix<double, NumNodes, 3> coordinate_field_data = Eigen::Matrix<double, NumNodes, 3>::Zero();
 
                 // Gather the field data for each node
                 for (size_t i = 0; i < num_nodes; ++i) {
@@ -162,7 +162,7 @@ class ElementForceProcessor {
                 size_t num_nodes = nodes.size();
 
                 // Set up the field data to gather
-                Kokkos::Array<Eigen::Matrix<double, NumNodes, 3>, NumFields> field_data_to_gather;
+                Kokkos::Array<Eigen::Matrix<double, NumNodes, 3>, NumFields> field_data_to_gather = {};
 
                 // Set up the results matrix
                 Eigen::Matrix<double, NumNodes, 3> results_to_scatter;

include/ElementNodeProcessor.h

@@ -16,6 +16,7 @@
 #include <stk_topology/topology.hpp>
 
 #include "AperiStkUtils.h"
+#include "Field.h"
 #include "FieldData.h"
 #include "LogUtils.h"
 #include "MeshData.h"
@@ -80,15 +81,15 @@ class ElementNodeProcessor {
                 size_t num_nodes = nodes.size();
 
                 // Allocate for the node indices
-                Kokkos::Array<stk::mesh::FastMeshIndex, HEX8_NUM_NODES> node_indices;
+                Kokkos::Array<aperi::Index, HEX8_NUM_NODES> node_indices;
 
                 // Get the node indices
                 for (size_t i = 0; i < num_nodes; ++i) {
-                    node_indices[i] = ngp_mesh.fast_mesh_index(nodes[i]);
+                    node_indices[i] = aperi::Index(ngp_mesh.fast_mesh_index(nodes[i]));
                 }
 
                 // Call the action function
-                func(elem_index, node_indices, num_nodes);
+                func(aperi::Index(elem_index), node_indices, num_nodes);
             });
     }
 
@@ -109,44 +110,4 @@ class ElementNodeProcessor {
     stk::mesh::NgpMesh m_ngp_mesh;                 // The ngp mesh object.
 };
 
-// Define the gather kernel
-template <typename T, size_t N>
-struct GatherKernel {
-    GatherKernel(const aperi::MeshData &mesh_data, const aperi::FieldQueryData<T> &field_query_data)
-        : field(stk::mesh::get_updated_ngp_field<T>(*StkGetField(field_query_data, mesh_data.GetMetaData()))) {
-    }
-
-    KOKKOS_FUNCTION void operator()(const stk::mesh::FastMeshIndex &index, Eigen::Vector<T, N> &data) const {
-        for (size_t i = 0; i < N; ++i) {
-            data(i) = field(index, i);
-        }
-    }
-
-   private:
-    stk::mesh::NgpField<T> field;
-};
-
-// Define a scatter kernel
-template <typename T, size_t N>
-struct ScatterKernel {
-    ScatterKernel(const aperi::MeshData &mesh_data, const aperi::FieldQueryData<T> &field_query_data)
-        : field(stk::mesh::get_updated_ngp_field<T>(*StkGetField(field_query_data, mesh_data.GetMetaData()))) {
-    }
-
-    KOKKOS_FUNCTION void AtomicAdd(const stk::mesh::FastMeshIndex &index, const Eigen::Vector<T, N> &data) const {
-        for (size_t i = 0; i < N; ++i) {
-            Kokkos::atomic_add(&field(index, i), data(i));
-        }
-    }
-
-    KOKKOS_FUNCTION void Overwrite(const stk::mesh::FastMeshIndex &index, const Eigen::Vector<T, N> &data) const {
-        for (size_t i = 0; i < N; ++i) {
-            field(index, i) = data(i);
-        }
-    }
-
-   private:
-    stk::mesh::NgpField<T> field;
-};
-
 }  // namespace aperi

include/ExternalForceContribution.h

@@ -6,7 +6,6 @@
 #include "EntityProcessor.h"
 #include "ForceContribution.h"
 #include "IoInputFile.h"
-#include "Material.h"
 #include "MeshData.h"
 
 namespace YAML {