Talamini/feature/rate dependence

examples/CMakeLists.txt

@@ -21,32 +21,12 @@ install(
 # Add the simple conduction examples
 add_subdirectory(simple_conduction)
 
-install(
-    FILES
-        simple_conduction/without_input_file.cpp
-    DESTINATION
-        examples/serac/simple_conduction
-)
-
 # Add the contact examples
 add_subdirectory(contact)
 
-install(
-    FILES
-        contact/beam_bending.cpp
-        contact/ironing.cpp
-        contact/sphere.cpp
-        contact/twist.cpp
-    DESTINATION
-        examples/serac/contact
-)
-
 # Add the buckling examples
 add_subdirectory(buckling)
 
-install(
-    FILES
-        buckling/cylinder.cpp
-    DESTINATION
-        examples/serac/buckling
-)
+# Add the uniaxial examples
+add_subdirectory(uniaxial)
+

examples/buckling/CMakeLists.txt

@@ -10,4 +10,12 @@ if(TRIBOL_FOUND AND PETSC_FOUND)
                         OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
                         DEPENDS_ON  serac_physics serac_mesh
                         )
+
+    install(
+        FILES
+            cylinder.cpp
+        DESTINATION
+            examples/serac/buckling
+    )
+
 endif()

examples/contact/CMakeLists.txt

@@ -5,24 +5,27 @@
 # SPDX-License-Identifier: (BSD-3-Clause)
 
 if(TRIBOL_FOUND)
-    blt_add_executable( NAME        contact_beam_bending
-                        SOURCES     beam_bending.cpp
-                        OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                        DEPENDS_ON  serac_physics serac_mesh
-                        )
-    blt_add_executable( NAME        contact_ironing
-                        SOURCES     ironing.cpp
-                        OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                        DEPENDS_ON  serac_physics serac_mesh
-                        )
-    blt_add_executable( NAME        contact_sphere
-                        SOURCES     sphere.cpp
-                        OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                        DEPENDS_ON  serac_physics serac_mesh
-                        )
-    blt_add_executable( NAME        contact_twist
-                        SOURCES     twist.cpp
-                        OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
-                        DEPENDS_ON  serac_physics serac_mesh
-                        )
+    set(CONTACT_EXAMPLES_SOURCES
+        beam_bending.cpp
+        ironing.cpp
+        sphere.cpp
+        twist.cpp
+    )
+
+    foreach(filename ${arg_CONTACT_EXAMPLES_SOURCES})
+        get_filename_component(example_name ${filename} NAME_WE)
+
+        blt_add_executable(NAME        contact_${example_name}
+                           SOURCES     ${filename}
+                           OUTPUT_DIR  ${EXAMPLE_OUTPUT_DIRECTORY}
+                           DEPENDS_ON  serac_physics serac_mesh)
+    endforeach()
+
+    install(
+        FILES
+            ${CONTACT_EXAMPLES_SOURCES}
+        DESTINATION
+            examples/serac/contact
+    )
+
 endif()

examples/simple_conduction/CMakeLists.txt

@@ -10,6 +10,13 @@ blt_add_executable( NAME        simple_conduction_without_input_file
                     DEPENDS_ON  serac_physics serac_mesh serac_functional
                     )
 
+install(
+    FILES
+        without_input_file.cpp
+    DESTINATION
+        examples/serac/simple_conduction
+)
+
 if(SERAC_ENABLE_TESTS)
     blt_add_test(NAME          simple_conduction_without_input_file
                  COMMAND       simple_conduction_without_input_file

examples/uniaxial/CMakeLists.txt

@@ -0,0 +1,20 @@
+# Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
+# other Serac Project Developers. See the top-level LICENSE file for
+# details.
+#
+# SPDX-License-Identifier: (BSD-3-Clause)
+
+blt_add_executable( NAME        uniaxial_tension
+                    SOURCES     uniaxial_tension.cpp
+                    OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY}
+                    DEPENDS_ON  serac_physics serac_mesh
+                    )
+
+# add serac_functional to DEPENDS_ON 
+
+install(
+    FILES
+        uniaxial_tension.cpp
+    DESTINATION
+        examples/serac/uniaxial
+)

examples/uniaxial/uniaxial_tension.cpp

@@ -0,0 +1,168 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
+// other Serac Project Developers. See the top-level LICENSE file for
+// details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#include <string>
+#include <fstream>
+
+#include "axom/inlet.hpp"
+#include "axom/slic/core/SimpleLogger.hpp"
+#include "mfem.hpp"
+
+#include "serac/infrastructure/terminator.hpp"
+#include "serac/mesh/mesh_utils.hpp"
+#include "serac/physics/boundary_conditions/components.hpp"
+#include "serac/physics/materials/solid_material.hpp"
+#include "serac/physics/solid_mechanics.hpp"
+#include "serac/physics/state/state_manager.hpp"
+#include "serac/serac_config.hpp"
+
+#include <cfenv>
+
+template <class Physics>
+void output(double u, double f, const Physics& solid, const std::string& paraview_tag, std::ofstream& file)
+{
+  solid.outputStateToDisk(paraview_tag);
+  file << solid.time() << " " << u << " " << f << std::endl;
+}
+
+int main(int argc, char* argv[])
+{
+  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+  serac::initialize(argc, argv);
+
+  constexpr int p = 2;
+  constexpr int dim = 3;
+  constexpr double x_length = 1.0;
+  constexpr double y_length = 0.1;
+  constexpr double z_length = 0.1;
+  constexpr int elements_in_x = 10;
+  constexpr int elements_in_y = 1;
+  constexpr int elements_in_z = 1;
+
+  int serial_refinements = 0;
+  int parallel_refinements = 0;
+  int time_steps = 100;
+  double strain_rate = 1e-3;
+
+  constexpr double E = 1.0;
+  constexpr double nu = 0.25;
+  constexpr double density = 1.0;
+
+  constexpr double sigma_y = 0.001;
+  constexpr double sigma_sat = 3.0 * sigma_y;
+  constexpr double strain_constant = 10 * sigma_y / E;
+  constexpr double eta = 1e-2;
+
+  constexpr double max_strain = 3 * strain_constant;
+  std::string output_filename = "uniaxial_fd.txt";
+
+  // Handle command line arguments
+  axom::CLI::App app{"Plane strain uniaxial extension of a bar."};
+  // Mesh options
+  app.add_option("--serial-refinements", serial_refinements, "Serial refinement steps", true);
+  app.add_option("--parallel-refinements", parallel_refinements, "Parallel refinement steps", true);
+  app.add_option("--time-steps", time_steps, "Number of time steps to divide simulation", true);
+  app.add_option("--strain-rate", strain_rate, "Nominal strain rate", true);
+  app.add_option("--output-file", output_filename, "Name for force-displacement output file", true);
+  app.set_help_flag("--help");
+
+  CLI11_PARSE(app, argc, argv);
+
+  SLIC_INFO_ROOT(axom::fmt::format("strain rate: {}", strain_rate));
+  SLIC_INFO_ROOT(axom::fmt::format("time_steps: {}", time_steps));
+
+  double max_time = max_strain / strain_rate;
+
+  // Create DataStore
+  const std::string simulation_tag = "uniaxial";
+  const std::string mesh_tag = simulation_tag + "mesh";
+  axom::sidre::DataStore datastore;
+  serac::StateManager::initialize(datastore, simulation_tag + "_data");
+
+  auto mesh = serac::mesh::refineAndDistribute(
+      serac::buildCuboidMesh(elements_in_x, elements_in_y, elements_in_z, x_length, y_length, z_length),
+      serial_refinements, parallel_refinements);
+  auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag);
+
+  // create boundary domains for boundary conditions
+  auto fix_x = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(5));
+  auto fix_y = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(2));
+  auto fix_z = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(1));
+  auto apply_displacement = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr<dim>(3));
+  serac::Domain whole_mesh = serac::EntireDomain(pmesh);
+
+  serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 0};
+
+  serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton,
+                                                  .relative_tol = 1.0e-10,
+                                                  .absolute_tol = 1.0e-12,
+                                                  .max_iterations = 200,
+                                                  .print_level = 1};
+
+  serac::SolidMechanics<p, dim> solid_solver(
+      nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, simulation_tag, mesh_tag);
+
+  using Hardening = serac::solid_mechanics::VoceHardening;
+  using Material = serac::solid_mechanics::J2<Hardening>;
+
+  Hardening hardening{sigma_y, sigma_sat, strain_constant, eta};
+  Material material{E, nu, hardening, density};
+
+  auto internal_states = solid_solver.createQuadratureDataBuffer(Material::State{}, whole_mesh);
+
+  solid_solver.setRateDependentMaterial(material, whole_mesh, internal_states);
+
+  solid_solver.setFixedBCs(fix_x, serac::Component::X);
+  solid_solver.setFixedBCs(fix_y, serac::Component::Y);
+  solid_solver.setFixedBCs(fix_z, serac::Component::Z);
+  auto applied_displacement = [strain_rate](serac::vec3, double t) {
+    return serac::vec3{strain_rate * x_length * t, 0., 0.};
+  };
+  solid_solver.setDisplacementBCs(applied_displacement, apply_displacement, serac::Component::X);
+
+  solid_solver.completeSetup();
+
+  double dt = max_time / (time_steps - 1);
+
+  // get nodes and dofs to compute total force
+  mfem::Array<int> dof_list = apply_displacement.dof_list(&solid_solver.displacement().space());
+  solid_solver.displacement().space().DofsToVDofs(0, dof_list);
+
+  auto compute_net_force = [&dof_list](const serac::FiniteElementDual& reaction) -> double {
+    double R{};
+    for (int i = 0; i < dof_list.Size(); i++) {
+      R += reaction(dof_list[i]);
+    }
+    return R;
+  };
+
+  std::string paraview_tag = simulation_tag + "_paraview";
+  std::ofstream file(output_filename);
+  file << "# time displacement force" << std::endl;
+  {
+    double u = applied_displacement(serac::vec3{}, solid_solver.time())[0];
+    double f = compute_net_force(solid_solver.dual("reactions"));
+    output(u, f, solid_solver, paraview_tag, file);
+  }
+
+  for (int i = 1; i < time_steps; ++i) {
+    SLIC_INFO_ROOT("------------------------------------------");
+    SLIC_INFO_ROOT(axom::fmt::format("TIME STEP {}", i));
+    SLIC_INFO_ROOT(axom::fmt::format("time = {} (out of {})", solid_solver.time() + dt, max_time));
+    serac::logger::flush();
+
+    solid_solver.advanceTimestep(dt);
+
+    double u = applied_displacement(serac::vec3{}, solid_solver.time())[0];
+    double f = compute_net_force(solid_solver.dual("reactions"));
+    output(u, f, solid_solver, paraview_tag, file);
+  }
+
+  file.close();
+  serac::exitGracefully();
+
+  return 0;
+}

src/serac/numerics/functional/tuple_tensor_dual_functions.hpp

@@ -861,10 +861,10 @@ inline SERAC_HOST_DEVICE tuple<vec3, mat3> eig_symm(const mat3& A)
 
     double A11 = dot(s0, A_dev_s0);
     double A12 = dot(s0, A_dev_s1);
-    double A21 = dot(s1, A_dev_s0);
+    double A21 = A12;
     double A22 = dot(s1, A_dev_s1);
 
-    double delta = 0.5 * sgn(A11 - A22) * std::sqrt((A11 - A22) * (A11 - A22) + 4 * A12 * A21);
+    double delta = 0.5 * std::sqrt((A11 - A22) * (A11 - A22) + 4 * A12 * A21);
 
     eta(1) = 0.5 * (A11 + A22) - delta;
     eta(2) = 0.5 * (A11 + A22) + delta;

src/serac/physics/base_physics.cpp

@@ -64,6 +64,7 @@ void BasePhysics::initializeBasePhysicsStates(int cycle, double time)
   timesteps_.clear();
 
   time_ = time;
+  dt_ = 0.0;
   max_time_ = time;
   min_time_ = time;
   cycle_ = cycle;

src/serac/physics/base_physics.hpp

@@ -618,6 +618,11 @@ class BasePhysics {
    */
   double time_;
 
+  /**
+   * @brief Current time step
+   */
+  double dt_;
+
   /**
    * @brief The maximum time reached for the forward solver
    */

src/serac/physics/materials/material_verification_tools.hpp

@@ -43,6 +43,7 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat
                           std::function<double(double)> epsilon_xx, const parameter_types... parameter_functions)
 {
   double t = 0;
+  const double dt = t_max / double(num_steps - 1);
 
   auto state = initial_state;
 
@@ -63,7 +64,6 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat
   output_history.reserve(num_steps);
 
   tensor<double, 3, 3> dudx{};
-  const double dt = t_max / double(num_steps - 1);
   for (size_t i = 0; i < num_steps; i++) {
     auto initial_guess = tensor<double, 2>{dudx[1][1], dudx[2][2]};
     auto epsilon_yy_and_zz = find_root(sigma_yy_and_zz, initial_guess);
@@ -80,6 +80,26 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat
   return output_history;
 }
 
+/**
+ * @brief Drive a rate-dependent material model thorugh a uniaxial tension experiment
+ *
+ * Drives material model through specified axial displacement gradient history.
+ * The time elaspses from 0 up to t_max.
+ * Currently only implemented for isotropic materials (or orthotropic materials with the
+ * principal axes aligned with the coordinate directions).
+ */
+template <typename MaterialType, typename StateType, typename... parameter_types>
+auto uniaxial_stress_test_rate_dependent(double t_max, size_t num_steps, const MaterialType material,
+                                         const StateType initial_state, std::function<double(double)> epsilon_xx,
+                                         const parameter_types... parameter_functions)
+{
+  const double dt = t_max / double(num_steps - 1);
+  auto mat_with_dt = [&material, dt](auto& state, auto du_dx, parameter_types... parameters) {
+    return material(state, dt, du_dx, parameters...);
+  };
+  return uniaxial_stress_test(t_max, num_steps, mat_with_dt, initial_state, epsilon_xx, parameter_functions...);
+}
+
 /**
  * @brief This function takes a material model (and associate state variables),
  *        subjects it to a time history of stimuli, described by `functions ... f`,
@@ -107,12 +127,12 @@ auto single_quadrature_point_test(double t_max, size_t num_steps, const Material
   const double dt = t_max / double(num_steps - 1);
   auto state = initial_state;
 
-  using output_type = decltype(std::tuple{t, state, f(0.0)..., decltype(material(state, f(0.0)...)){}});
+  using output_type = decltype(std::tuple{t, state, f(0.0)..., decltype(material(state, dt, f(0.0)...)){}});
   std::vector<output_type> history;
   history.reserve(num_steps);
 
   for (size_t i = 0; i < num_steps; i++) {
-    auto material_output = material(state, f(t)...);
+    auto material_output = material(state, dt, f(t)...);
     history.push_back(std::tuple{t, state, f(t)..., material_output});
     t += dt;
   }