Gretl solid mechanics

src/smith/differentiable_numerics/CMakeLists.txt

@@ -6,20 +6,32 @@
 
 set(differentiable_numerics_sources
     field_state.cpp
-    differentiable_utils.cpp
+    differentiable_test_utils.cpp
     differentiable_physics.cpp
     lumped_mass_explicit_newmark_state_advancer.cpp
+    solid_mechanics_state_advancer.cpp
+    differentiable_solver.cpp
+    nonlinear_solve.cpp
+    dirichlet_boundary_conditions.cpp
     )
 
 set(differentiable_numerics_headers
     field_state.hpp
     state_advancer.hpp
+    reaction.hpp
+    differentiable_solver.hpp
+    differentiable_physics.hpp
+    differentiable_test_utils.hpp
     timestep_estimator.hpp
     explicit_dynamic_solve.hpp
     lumped_mass_weak_form.hpp
-    differentiable_utils.hpp
-    differentiable_physics.hpp
     lumped_mass_explicit_newmark_state_advancer.hpp
+    solid_mechanics_state_advancer.hpp
+    nonlinear_solve.hpp
+    dirichlet_boundary_conditions.hpp
+    time_integration_rule.hpp
+    time_discretized_weak_form.hpp
+    differentiable_solid_mechanics.hpp
     )
 
 set(differentiable_numerics_depends

src/smith/differentiable_numerics/differentiable_physics.cpp

@@ -1,13 +1,10 @@
-// Copyright (c) Lawrence Livermore National Security, LLC and
-// other Smith Project Developers. See the top-level LICENSE file for
-// details.
-//
-// SPDX-License-Identifier: (BSD-3-Clause)
-
+#include "gretl/data_store.hpp"
 #include "smith/differentiable_numerics/differentiable_physics.hpp"
 #include "smith/physics/weak_form.hpp"
 #include "smith/physics/mesh.hpp"
+#include "smith/differentiable_numerics/differentiable_physics.hpp"
 #include "smith/differentiable_numerics/state_advancer.hpp"
+#include "smith/differentiable_numerics/reaction.hpp"
 #include "gretl/data_store.hpp"
 
 namespace smith {
@@ -36,10 +33,12 @@ gretl::State<int> make_milestone(const std::vector<FieldState>& states)
 DifferentiablePhysics::DifferentiablePhysics(std::shared_ptr<Mesh> mesh, std::shared_ptr<gretl::DataStore> graph,
                                              const FieldState& shape_disp, const std::vector<FieldState>& states,
                                              const std::vector<FieldState>& params,
-                                             std::shared_ptr<StateAdvancer> advancer, std::string mech_name)
+                                             std::shared_ptr<StateAdvancer> advancer, std::string mech_name,
+                                             const std::vector<std::string>& resultant_names)
     : BasePhysics(mech_name, mesh, 0, 0.0, false),  // the false is checkpoint_to_disk
       checkpointer_(graph),
-      advancer_(advancer)
+      advancer_(advancer),
+      resultant_names_(resultant_names)
 {
   SLIC_ERROR_IF(states.size() == 0, "Must have a least 1 state for a mechanics.");
   field_shape_displacement_ = std::make_unique<FieldState>(shape_disp);
@@ -48,14 +47,18 @@ DifferentiablePhysics::DifferentiablePhysics(std::shared_ptr<Mesh> mesh, std::sh
     field_states_.push_back(s);
     initial_field_states_.push_back(s);
     state_name_to_field_index_[s.get()->name()] = i;
-    state_names.push_back(s.get()->name());
+    state_names_.push_back(s.get()->name());
   }
 
   for (size_t i = 0; i < params.size(); ++i) {
     const auto& p = params[i];
     field_params_.push_back(p);
     param_name_to_field_index_[p.get()->name()] = i;
-    param_names.push_back(p.get()->name());
+    param_names_.push_back(p.get()->name());
+  }
+
+  for (size_t i = 0; i < resultant_names_.size(); ++i) {
+    resultant_name_to_resultant_index_[resultant_names_[i]] = i;
   }
 
   completeSetup();
@@ -84,9 +87,11 @@ void DifferentiablePhysics::resetAdjointStates()
   gretl_assert(checkpointer_->check_validity());
 }
 
-std::vector<std::string> DifferentiablePhysics::stateNames() const { return state_names; }
+std::vector<std::string> DifferentiablePhysics::stateNames() const { return state_names_; }
 
-std::vector<std::string> DifferentiablePhysics::parameterNames() const { return param_names; }
+std::vector<std::string> DifferentiablePhysics::parameterNames() const { return param_names_; }
+
+std::vector<std::string> DifferentiablePhysics::dualNames() const { return resultant_names_; }
 
 const FiniteElementState& DifferentiablePhysics::state([[maybe_unused]] const std::string& field_name) const
 {
@@ -99,6 +104,21 @@ const FiniteElementState& DifferentiablePhysics::state([[maybe_unused]] const st
   return *field_states_[state_index].get();
 }
 
+const FiniteElementDual& DifferentiablePhysics::dual(const std::string& dual_name) const
+{
+  SLIC_ERROR_IF(
+      resultant_name_to_resultant_index_.find(dual_name) == resultant_name_to_resultant_index_.end(),
+      axom::fmt::format("Could not find dual named {0} in mesh with tag \"{1}\" to get", dual_name, mesh_->tag()));
+  size_t reaction_index = resultant_name_to_resultant_index_.at(dual_name);
+  SLIC_ERROR_IF(
+      reaction_index >= resultant_names_.size(),
+      "Dual reactions not correctly allocated yet, cannot get dual until after initializationStep is called.");
+
+  TimeInfo time_info(time_old_, dt_old_, static_cast<size_t>(cycle_old_));
+  resultant_states_ = advancer_->computeResultants(*field_shape_displacement_, field_states_, field_params_, time_info);
+  return *resultant_states_[reaction_index].get();
+}
+
 FiniteElementState DifferentiablePhysics::loadCheckpointedState(const std::string& state_name, int cycle)
 {
   SLIC_ERROR_IF(cycle != cycle_,
@@ -154,14 +174,28 @@ void DifferentiablePhysics::setAdjointLoad(
   for (auto string_dual_pair : string_to_dual) {
     std::string field_name = string_dual_pair.first;
     const smith::FiniteElementDual& dual = string_dual_pair.second;
-    SLIC_ERROR_IF(
-        state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(),
-        axom::fmt::format("Could not find dual named {0} in mesh with tag {1} to set", field_name, mesh_->tag()));
+    SLIC_ERROR_IF(state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(),
+                  axom::fmt::format("Could not find dual named {0} in mesh with tag {1}", field_name, mesh_->tag()));
     size_t state_index = state_name_to_field_index_.at(field_name);
     *field_states_[state_index].get_dual() += dual;
   }
 }
 
+void DifferentiablePhysics::setDualAdjointBcs(
+    std::unordered_map<std::string, const smith::FiniteElementState&> string_to_bc)
+{
+  for (auto string_bc_pair : string_to_bc) {
+    std::string reaction_name = string_bc_pair.first;
+    const smith::FiniteElementState& reaction_dual = string_bc_pair.second;
+    SLIC_ERROR_IF(
+        resultant_name_to_resultant_index_.find(reaction_name) == resultant_name_to_resultant_index_.end(),
+        axom::fmt::format("When calling setDualAdjointBcs, could not find reaction named {0} in mesh with tag {1}",
+                          reaction_name, mesh_->tag()));
+    size_t reaction_index = resultant_name_to_resultant_index_.at(reaction_name);
+    *resultant_states_[reaction_index].get_dual() += reaction_dual;
+  }
+}
+
 const FiniteElementState& DifferentiablePhysics::adjoint([[maybe_unused]] const std::string& adjoint_name) const
 {
   // MRT, not implemented
@@ -176,6 +210,10 @@ void DifferentiablePhysics::advanceTimestep(double dt)
     milestones_.push_back(make_milestone(field_states_).step());
   }
 
+  cycle_old_ = cycle_;
+  time_old_ = time_;
+  dt_old_ = dt;
+
   TimeInfo time_info(time_, dt, static_cast<size_t>(cycle_));
   field_states_ = advancer_->advanceState(*field_shape_displacement_, field_states_, field_params_, time_info);
 
@@ -232,7 +270,7 @@ DifferentiablePhysics::computeInitialConditionSensitivity() const
   return map;
 }
 
-std::vector<FieldState> DifferentiablePhysics::getAllFieldStates() const
+std::vector<FieldState> DifferentiablePhysics::getFieldStatesAndParamStates() const
 {
   std::vector<FieldState> fields;
   fields.insert(fields.end(), field_states_.begin(), field_states_.end());

src/smith/differentiable_numerics/differentiable_physics.hpp

@@ -26,6 +26,7 @@ class WeakForm;
 class DifferentiableSolver;
 class StateAdvancer;
 class TimestepEstimator;
+class Reaction;
 
 /// @brief Implementation of BasePhysics which uses FieldStates and gretl to track the computational graph, dynamically
 /// checkpoint, and back-propagate sensitivities.
@@ -35,7 +36,7 @@ class DifferentiablePhysics : public BasePhysics {
   DifferentiablePhysics(std::shared_ptr<Mesh> mesh, std::shared_ptr<gretl::DataStore> graph,
                         const FieldState& shape_disp, const std::vector<FieldState>& states,
                         const std::vector<FieldState>& params, std::shared_ptr<StateAdvancer> advancer,
-                        std::string mech_name);
+                        std::string physics_name, const std::vector<std::string>& resultant_names = {});
   /// @brief destructor
   ~DifferentiablePhysics() {}
 
@@ -54,11 +55,14 @@ class DifferentiablePhysics : public BasePhysics {
   /// @overload
   std::vector<std::string> parameterNames() const override;
 
+  /// @overload
+  std::vector<std::string> dualNames() const override;
+
   /// @overload
   const FiniteElementState& state(const std::string& state_name) const override;
 
   /// @overload
-  FiniteElementState loadCheckpointedState(const std::string& state_name, int cycle) override;
+  const FiniteElementDual& dual(const std::string& dual_name) const override;
 
   /// @overload
   const FiniteElementState& shapeDisplacement() const override;
@@ -69,6 +73,9 @@ class DifferentiablePhysics : public BasePhysics {
   /// @overload
   const FiniteElementState& parameter(const std::string& parameter_name) const override;
 
+  /// @overload
+  FiniteElementState loadCheckpointedState(const std::string& state_name, int cycle) override;
+
   /// @overload
   void setState(const std::string& state_name, const FiniteElementState& s) override;
 
@@ -81,6 +88,9 @@ class DifferentiablePhysics : public BasePhysics {
   /// @overload
   void setAdjointLoad(std::unordered_map<std::string, const smith::FiniteElementDual&> string_to_dual) override;
 
+  /// @overload
+  void setDualAdjointBcs(std::unordered_map<std::string, const smith::FiniteElementState&> string_to_bc) override;
+
   /// @overload
   const FiniteElementState& adjoint(const std::string& adjoint_name) const override;
 
@@ -100,12 +110,27 @@ class DifferentiablePhysics : public BasePhysics {
   const std::unordered_map<std::string, const smith::FiniteElementDual&> computeInitialConditionSensitivity()
       const override;
 
+  /// @brief Get all the initial FieldStates
+  std::vector<FieldState> getInitialFieldStates() const { return initial_field_states_; }
+
+  /// @brief Get all the current FieldStates
+  std::vector<FieldState> getFieldStates() const { return field_states_; }
+
+  /// @brief Get all the parameter FieldStates
+  std::vector<FieldState> getFieldParams() const { return field_params_; }
+
   /// @brief Get all the FieldStates... states first, parameters next
-  std::vector<FieldState> getAllFieldStates() const;
+  std::vector<FieldState> getFieldStatesAndParamStates() const;
 
   /// @brief Get the shape displacement FieldState
   FieldState getShapeDispFieldState() const;
 
+  /// @brief Get the current ResultantStates
+  std::vector<ResultantState> getResultantStates() const { return resultant_states_; }
+
+  /// @brief Get state advancer
+  std::shared_ptr<StateAdvancer> getStateAdvancer() const { return advancer_; }
+
  private:
   std::shared_ptr<gretl::DataStore> checkpointer_;  ///< gretl data store manages dynamic checkpointing logic
   std::shared_ptr<StateAdvancer> advancer_;  ///< abstract interface for advancing state from one cycle to the next
@@ -119,13 +144,21 @@ class DifferentiablePhysics : public BasePhysics {
 
   std::map<std::string, size_t> state_name_to_field_index_;  ///< map from state names to field index
   std::map<std::string, size_t> param_name_to_field_index_;  ///< map from param names to param index
-  std::vector<std::string> state_names;                      ///< names of all the states in order
-  std::vector<std::string> param_names;                      ///< names of all the states in order
+  std::vector<std::string> state_names_;                     ///< names of all the states in order
+  std::vector<std::string> param_names_;                     ///< names of all the states in order
+
+  mutable std::vector<ResultantState> resultant_states_;             ///< all the resultants registered for the physics
+  std::map<std::string, size_t> resultant_name_to_resultant_index_;  ///< map from reaction names to resultant index
+  std::vector<std::string> resultant_names_;  ///< names for all the relevant resultants/reactions
 
   std::vector<gretl::Int> milestones_;  ///< a record of the steps in the graph that represent the end conditions of
                                         ///< advanceTimestep(dt). this information is used to halt the gretl graph when
                                         ///< back-propagating to allow users of reverseAdjointTimestep to specify
                                         ///< adjoint loads and to retrieve timestep sensitivity information.
+
+  double time_old_ = 0.0;
+  double dt_old_ = 0.0;
+  int cycle_old_ = 0;
 };
 
 }  // namespace smith

src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp

@@ -0,0 +1,49 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
+// other smith Project Developers. See the top-level LICENSE file for
+// details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+/**
+ * @file differentiable_solid_mechanics.hpp
+ * .hpp
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include "smith/differentiable_numerics/solid_mechanics_state_advancer.hpp"
+#include "smith/differentiable_numerics/time_discretized_weak_form.hpp"
+#include "smith/differentiable_numerics/differentiable_physics.hpp"
+
+namespace smith {
+
+/// Helper function to generate the base-physics for solid mechanics
+/// This will return a tuple of shared pointers to the: BasePhysics, WeakForm, and DirichetBoundaryConditions
+/// Only the BasePhysics needs to stay in scope.  The others are returned to the user so they can define the WeakForm
+/// integrals, and to specify space and time varying boundary conditions
+template <int dim, typename ShapeDispSpace, typename VectorSpace, typename... ParamSpaces>
+auto buildSolidMechanics(std::shared_ptr<smith::Mesh> mesh,
+                         std::shared_ptr<DifferentiableSolver> d_solid_nonlinear_solver,
+                         smith::SecondOrderTimeIntegrationRule time_rule, std::string physics_name,
+                         const std::vector<std::string>& param_names = {})
+{
+  auto graph = std::make_shared<gretl::DataStore>(100);
+  auto [shape_disp, states, params, time, solid_mechanics_weak_form] =
+      SolidMechanicsStateAdvancer::buildWeakFormAndStates<dim, ShapeDispSpace, VectorSpace, ParamSpaces...>(
+          mesh, graph, time_rule, physics_name, param_names);
+
+  auto vector_bcs = std::make_shared<DirichletBoundaryConditions>(
+      mesh->mfemParMesh(), space(states[SolidMechanicsStateAdvancer::DISPLACEMENT]));
+
+  auto state_advancer = std::make_shared<SolidMechanicsStateAdvancer>(d_solid_nonlinear_solver, vector_bcs,
+                                                                      solid_mechanics_weak_form, time_rule);
+
+  auto physics = std::make_shared<DifferentiablePhysics>(mesh, graph, shape_disp, states, params, state_advancer,
+                                                         physics_name, std::vector<std::string>{"reactions"});
+
+  return std::make_tuple(physics, solid_mechanics_weak_form, vector_bcs);
+}
+
+}  // namespace smith