Gretl solid mechanics

examples/inertia_relief/inertia_relief_example.cpp

@@ -80,37 +80,6 @@ class ParaviewWriter {
   StateVecs dual_states;
 };
 
-auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector<smith::FiniteElementState const*>& states,
-                          std::string output_name)
-{
-  if (output_name == "") {
-    output_name = "default";
-  }
-
-  ParaviewWriter::StateVecs output_states;
-  for (const auto& s : states) {
-    output_states.push_back(std::make_shared<smith::FiniteElementState>(s->space(), s->name()));
-  }
-
-  auto non_const_mesh = const_cast<mfem::ParMesh*>(&mesh);
-  auto paraview_dc = std::make_unique<mfem::ParaViewDataCollection>(output_name, non_const_mesh);
-  int max_order_in_fields = 0;
-
-  // Find the maximum polynomial order in the physics module's states
-  for (const auto& state : output_states) {
-    paraview_dc->RegisterField(state->name(), &state->gridFunction());
-    max_order_in_fields = std::max(max_order_in_fields, state->space().GetOrder(0));
-  }
-
-  // Set the options for the paraview output files
-  paraview_dc->SetLevelsOfDetail(max_order_in_fields);
-  paraview_dc->SetHighOrderOutput(true);
-  paraview_dc->SetDataFormat(mfem::VTKFormat::BINARY);
-  paraview_dc->SetCompression(true);
-
-  return ParaviewWriter(std::move(paraview_dc), output_states, {});
-}
-
 /* Nonlinear problem of the form
  * F(X) = [ r(u) + (dc/du)^T l ] = [ 0 ]
  *        [ -c(u)              ]   [ 0 ]
@@ -302,7 +271,7 @@ int main(int argc, char* argv[])
     return u;
   });
 
-  auto writer = createParaviewOutput(mesh->mfemParMesh(), objective_states, "inertia_relief");
+  auto writer = createParaviewWriter(mesh->mfemParMesh(), objective_states, "inertia_relief");
   if (visualize) {
     writer.write(0, 0.0, objective_states);
   }

src/smith/differentiable_numerics/CMakeLists.txt

@@ -6,20 +6,34 @@
 
 set(differentiable_numerics_sources
     field_state.cpp
-    differentiable_utils.cpp
     differentiable_physics.cpp
     lumped_mass_explicit_newmark_state_advancer.cpp
+    solid_mechanics_state_advancer.cpp
+    differentiable_solver.cpp
+    nonlinear_solve.cpp
+    evaluate_objective.cpp
+    dirichlet_boundary_conditions.cpp
     )
 
 set(differentiable_numerics_headers
     field_state.hpp
     state_advancer.hpp
+    reaction.hpp
+    differentiable_solver.hpp
+    differentiable_physics.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
+    evaluate_objective.hpp
+    dirichlet_boundary_conditions.hpp
+    time_integration_rule.hpp
+    time_discretized_weak_form.hpp
+    differentiable_solid_mechanics.hpp
+    paraview_writer.hpp
+    differentiable_test_utils.hpp
     )
 
 set(differentiable_numerics_depends

src/smith/differentiable_numerics/differentiable_physics.cpp

@@ -4,10 +4,19 @@
 //
 // SPDX-License-Identifier: (BSD-3-Clause)
 
+// 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 +45,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>& reaction_names)
     : BasePhysics(mech_name, mesh, 0, 0.0, false),  // the false is checkpoint_to_disk
       checkpointer_(graph),
-      advancer_(advancer)
+      advancer_(advancer),
+      reaction_names_(reaction_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 +59,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 < reaction_names_.size(); ++i) {
+    reaction_name_to_reaction_index_[reaction_names_[i]] = i;
   }
 
   completeSetup();
@@ -84,21 +99,36 @@ 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 reaction_names_; }
 
 const FiniteElementState& DifferentiablePhysics::state([[maybe_unused]] const std::string& field_name) const
 {
   SLIC_ERROR_IF(
       state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(),
       axom::fmt::format("Could not find field named {0} in mesh with tag \"{1}\" to get", field_name, mesh_->tag()));
   size_t state_index = state_name_to_field_index_.at(field_name);
-  SLIC_ERROR_IF(state_index >= field_states_.size(),
-                "Field states not correctly allocated yet, cannot get state until after initializationStep is called.");
   return *field_states_[state_index].get();
 }
 
+const FiniteElementDual& DifferentiablePhysics::dual(const std::string& dual_name) const
+{
+  SLIC_ERROR_IF(
+      reaction_name_to_reaction_index_.find(dual_name) == reaction_name_to_reaction_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 = reaction_name_to_reaction_index_.at(dual_name);
+  SLIC_ERROR_IF(
+      reaction_index >= reaction_names_.size(),
+      "Dual reactions not correctly allocated yet, cannot get dual until after initializationStep is called.");
+
+  TimeInfo time_info(time_prev_, dt_prev_, static_cast<size_t>(cycle_prev_));
+  reaction_states_ = advancer_->computeReactions(time_info, *field_shape_displacement_, field_states_, field_params_);
+  return *reaction_states_[reaction_index].get();
+}
+
 FiniteElementState DifferentiablePhysics::loadCheckpointedState(const std::string& state_name, int cycle)
 {
   SLIC_ERROR_IF(cycle != cycle_,
@@ -154,14 +184,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(
+        reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_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 = reaction_name_to_reaction_index_.at(reaction_name);
+    *reaction_states_[reaction_index].get_dual() += reaction_dual;
+  }
+}
+
 const FiniteElementState& DifferentiablePhysics::adjoint([[maybe_unused]] const std::string& adjoint_name) const
 {
   // MRT, not implemented
@@ -176,8 +220,12 @@ void DifferentiablePhysics::advanceTimestep(double dt)
     milestones_.push_back(make_milestone(field_states_).step());
   }
 
+  cycle_prev_ = cycle_;
+  time_prev_ = time_;
+  dt_prev_ = dt;
+
   TimeInfo time_info(time_, dt, static_cast<size_t>(cycle_));
-  field_states_ = advancer_->advanceState(*field_shape_displacement_, field_states_, field_params_, time_info);
+  field_states_ = advancer_->advanceState(time_info, *field_shape_displacement_, field_states_, field_params_);
 
   cycle_++;
   time_ += dt;
@@ -232,7 +280,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>& reaction_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 reactionStates
+  std::vector<ReactionState> getReactionStates() const { return reaction_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,24 @@ 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<ReactionState> reaction_states_;             ///< all the reactions registered for the physics
+  std::map<std::string, size_t> reaction_name_to_reaction_index_;  ///< map from reaction names to reaction index
+  std::vector<std::string> reaction_names_;                        ///< names for all the relevant reactions/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_prev_ =
+      0.0;  ///< previous time, saved to reconstruct the start of step time used in computing reaction forces
+  double dt_prev_ =
+      0.0;  ///< previous time increment, saved to reconstruct the start of step time used in computing reaction forces
+  int cycle_prev_ =
+      0;  ///< previous cycle, saved to reconstruct the start of step time used in computing reaction forces
 };
 
 }  // namespace smith