Li/differentiable plastic system

RELEASE-NOTES.md

@@ -31,6 +31,7 @@ The Smith project release numbers follow [Semantic Versioning](http://semver.org
 - Added axisymmetric solid mechanics materials and loads for 2D `(r, z)` meshes.
 - Added host-configs for CZMatrix and RZVector machines
 - Added CI pipelines for CZMatrix machine
+- Added differentiable plastic system with internal variables represented by discontinuous galerkin nodal fields
 
 ### Removed
 

src/smith/differentiable_numerics/CMakeLists.txt

@@ -47,6 +47,7 @@ set(differentiable_numerics_headers
     combined_system.hpp
     system_base.hpp
     differentiable_test_utils.hpp
+    plastic_mechanics_system.hpp
     )
 
 set(differentiable_numerics_depends

src/smith/differentiable_numerics/differentiable_physics.cpp

@@ -67,9 +67,8 @@ DifferentiablePhysics::DifferentiablePhysics(std::shared_ptr<Mesh> mesh, std::sh
 
   reaction_names_.reserve(reaction_infos_.size());
   for (size_t i = 0; i < reaction_infos_.size(); ++i) {
-    SLIC_ERROR_IF(
-        reaction_infos_[i].space == nullptr,
-        axom::fmt::format("Dual '{}' in physics module '{}' has null FE space.", reaction_infos_[i].name, name_));
+    SLIC_ERROR_IF(reaction_infos_[i].space == nullptr,
+                  std::format("Dual '{}' in physics module '{}' has null FE space.", reaction_infos_[i].name, name_));
     reaction_names_.push_back(reaction_infos_[i].name);
     reaction_name_to_reaction_index_[reaction_infos_[i].name] = i;
   }
@@ -135,7 +134,7 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& reaction
   if (reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end()) {
     std::string available;
     for (auto& n : reaction_names_) available += n + " ";
-    SLIC_ERROR(axom::fmt::format(
+    SLIC_ERROR(std::format(
         "Could not find reaction named {0} in mesh with tag \"{1}\" to get. Available reactions (size {2}): {3}",
         reaction_name, mesh_->tag(), reaction_names_.size(), available));
   }
@@ -153,11 +152,10 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& reaction
 
 FiniteElementDual DifferentiablePhysics::loadCheckpointedDual(const std::string& reaction_name, int cycle)
 {
-  SLIC_ERROR_IF(
-      cycle != cycle_,
-      axom::fmt::format("Due to checkpointing restrictions in smith::DifferentiablePhysics, cannot ask for "
-                        "an arbitrary checkpointed reaction cycle, asking for cycle {}, but physics is at cycle {}",
-                        cycle, cycle_));
+  SLIC_ERROR_IF(cycle != cycle_,
+                std::format("Due to checkpointing restrictions in smith::DifferentiablePhysics, cannot ask for "
+                            "an arbitrary checkpointed reaction cycle, asking for cycle {}, but physics is at cycle {}",
+                            cycle, cycle_));
   return dual(reaction_name);
 }
 
@@ -218,9 +216,8 @@ void DifferentiablePhysics::setAdjointLoad(
   for (auto string_reaction_pair : string_to_reaction) {
     std::string field_name = string_reaction_pair.first;
     const smith::FiniteElementDual& reaction = string_reaction_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 reaction named {0} in mesh with tag {1}", field_name, mesh_->tag()));
+    SLIC_ERROR_IF(state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(),
+                  std::format("Could not find reaction 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() += reaction;
   }
@@ -232,10 +229,9 @@ void DifferentiablePhysics::setDualAdjointBcs(
   for (auto string_bc_pair : string_to_bc) {
     std::string reaction_name = string_bc_pair.first;
     const smith::FiniteElementState& reaction_adjoint_load = 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()));
+    SLIC_ERROR_IF(reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end(),
+                  std::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_adjoint_load;
   }

src/smith/differentiable_numerics/differentiable_test_utils.hpp

@@ -79,6 +79,7 @@ inline double checkGradWrt(const gretl::State<double>& objective, smith::FieldSt
 
   std::vector<double> grad_errors;
   auto [grad, grad_fd] = checkGradients(objective, input, *dual_vec, objectiveBase, graph, eps);
+  if (printmore) std::cout << "grad    = " << grad << "\ngrad fd = " << grad_fd << std::endl;
   grad_errors.push_back(std::abs(grad - grad_fd));
 
   for (size_t step = 1; step < num_fd_steps; ++step) {

src/smith/differentiable_numerics/field_store.cpp

@@ -31,16 +31,6 @@ std::shared_ptr<DirichletBoundaryConditions> FieldStore::addBoundaryConditions(F
   return std::make_shared<DirichletBoundaryConditions>(mesh_->mfemParMesh(), field->space());
 }
 
-void FieldStore::addWeakFormUnknownArg(std::string weak_form_name, std::string argument_name, size_t argument_index)
-{
-  FieldLabel argument_name_and_index{.field_name = argument_name, .field_index_in_residual = argument_index};
-  if (weak_form_name_to_unknown_name_index_.count(weak_form_name)) {
-    weak_form_name_to_unknown_name_index_.at(weak_form_name).push_back(argument_name_and_index);
-  } else {
-    weak_form_name_to_unknown_name_index_[weak_form_name] = std::vector<FieldLabel>{argument_name_and_index};
-  }
-}
-
 void FieldStore::addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index)
 {
   // Store the field name instead of index to avoid confusion between states_ and params_ indices
@@ -55,71 +45,43 @@ void FieldStore::addWeakFormArg(std::string weak_form_name, std::string argument
 
 void FieldStore::printMap()
 {
-  for (auto& keyval : weak_form_name_to_unknown_name_index_) {
+  for (const auto& keyval : weak_form_name_to_field_names_) {
     std::cout << "for residual: " << keyval.first << " ";
-    for (auto& name_index : keyval.second) {
-      std::cout << "arg " << name_index.field_name << " at " << name_index.field_index_in_residual << ", ";
+    for (size_t argument_index = 0; argument_index < keyval.second.size(); ++argument_index) {
+      std::cout << "arg " << keyval.second[argument_index] << " at " << argument_index << ", ";
     }
     std::cout << std::endl;
   }
 }
 
-std::vector<std::vector<size_t>> FieldStore::indexMap(const std::vector<std::string>& residual_names) const
+BlockArgumentMap FieldStore::indexMap(const std::vector<std::string>& residual_names) const
 {
-  // Build a local column space: each residual in the subsystem contributes one local column,
-  // corresponding to its "self" diagonal unknown.  The self-unknown is preferably the residual's
-  // reaction (test) field if that field appears in the unknown-arg list for this weak form;
-  // otherwise fall back on the first unknown argument (handles cases like the cycle-zero
-  // acceleration solve, where the reaction field is a dependent/history field).
   std::map<size_t, size_t> global_state_to_local_col;
   for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) {
     const std::string& res_name = residual_names[res_i];
-    size_t global_state_idx = invalid_block_index;
-
-    std::string reaction_name;
-    for (const auto& kv : weak_form_to_test_field_) {
-      if (kv.first == res_name) {
-        reaction_name = kv.second;
-        break;
-      }
-    }
-
-    // Check if the reaction field is one of the registered unknown args for this weak form.
-    bool reaction_is_unknown = false;
-    if (!reaction_name.empty() && weak_form_name_to_unknown_name_index_.count(res_name)) {
-      for (const auto& label : weak_form_name_to_unknown_name_index_.at(res_name)) {
-        if (label.field_name == reaction_name) {
-          reaction_is_unknown = true;
-          break;
-        }
-      }
-    }
-
-    if (reaction_is_unknown) {
-      global_state_idx = to_states_index_.at(reaction_name);
-    } else {
-      const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name);
-      SLIC_ERROR_IF(arg_info.empty(),
-                    "Weak form '" << res_name << "' has no unknown arguments; cannot build index map.");
-      global_state_idx = to_states_index_.at(arg_info.front().field_name);
-    }
-    global_state_to_local_col[global_state_idx] = res_i;
+    size_t global_state_idx = to_states_index_.at(getWeakFormFieldNames(res_name).unknown);
+    auto [owner, inserted] = global_state_to_local_col.emplace(global_state_idx, res_i);
+    SLIC_ERROR_IF(!inserted, std::format("Weak forms '{}' and '{}' both own state index {}",
+                                         residual_names[owner->second], res_name, global_state_idx));
   }
 
-  std::vector<std::vector<size_t>> block_indices(residual_names.size());
+  BlockArgumentMap block_indices(residual_names.size(),
+                                 std::vector<BlockArgumentIndices>(residual_names.size(), BlockArgumentIndices{}));
   for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) {
-    std::vector<size_t>& res_indices = block_indices[res_i];
-    res_indices = std::vector<size_t>(residual_names.size(), invalid_block_index);
     const std::string& res_name = residual_names[res_i];
-    const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name);
+    const auto& field_names = weak_form_name_to_field_names_.at(res_name);
+    size_t state_slot = 0;
+    for (const auto& field_name : field_names) {
+      if (!to_states_index_.count(field_name)) {
+        continue;
+      }
 
-    for (const auto& field_name_and_arg_index : arg_info) {
-      size_t global_state_index = to_states_index_.at(field_name_and_arg_index.field_name);
-      auto it = global_state_to_local_col.find(global_state_index);
-      if (it != global_state_to_local_col.end()) {
-        res_indices[it->second] = field_name_and_arg_index.field_index_in_residual;
+      size_t global_state_index = to_states_index_.at(field_name);
+      auto column = global_state_to_local_col.find(global_state_index);
+      if (column != global_state_to_local_col.end()) {
+        block_indices[res_i][column->second].push_back(state_slot);
       }
-      // else: field belongs to a different subsystem; treat as fixed input here.
+      ++state_slot;
     }
   }
 
@@ -132,7 +94,7 @@ std::vector<const BoundaryConditionManager*> FieldStore::getBoundaryConditionMan
   std::vector<std::string> field_names;
   field_names.reserve(weak_form_names.size());
   for (const auto& wf_name : weak_form_names) {
-    field_names.push_back(getWeakFormReaction(wf_name));
+    field_names.push_back(getWeakFormFieldNames(wf_name).unknown);
   }
   return getBoundaryConditionManagersForFields(field_names);
 }
@@ -277,17 +239,17 @@ const std::vector<FieldState>& FieldStore::getAllFields() const { return states_
 std::vector<FieldState> FieldStore::getStates(const std::string& weak_form_name) const
 {
   // Validate that weak form is registered
-  SLIC_ERROR_ROOT_IF(weak_form_name_to_field_names_.count(weak_form_name) == 0,
-                     axom::fmt::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?",
-                                       weak_form_name));
+  SLIC_ERROR_ROOT_IF(
+      weak_form_name_to_field_names_.count(weak_form_name) == 0,
+      std::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", weak_form_name));
 
   auto field_names = weak_form_name_to_field_names_.at(weak_form_name);
   std::vector<FieldState> fields_for_residual;
   for (auto& name : field_names) {
     // Validate that field exists
     SLIC_ERROR_ROOT_IF(
         to_states_index_.count(name) == 0 && to_params_index_.count(name) == 0,
-        axom::fmt::format("Field '{}' (required by weak form '{}') not found in FieldStore", name, weak_form_name));
+        std::format("Field '{}' (required by weak form '{}') not found in FieldStore", name, weak_form_name));
 
     // Only include state fields, not parameters
     // Parameters are passed separately to avoid duplication in block_solve
@@ -303,31 +265,31 @@ std::vector<FieldState> FieldStore::getStatesFromVectors(const std::string& weak
                                                          const std::vector<FieldState>& param_fields) const
 {
   // Validate that weak form is registered
-  SLIC_ERROR_ROOT_IF(weak_form_name_to_field_names_.count(weak_form_name) == 0,
-                     axom::fmt::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?",
-                                       weak_form_name));
+  SLIC_ERROR_ROOT_IF(
+      weak_form_name_to_field_names_.count(weak_form_name) == 0,
+      std::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", weak_form_name));
 
   auto field_names = weak_form_name_to_field_names_.at(weak_form_name);
   std::vector<FieldState> fields_for_residual;
   for (auto& name : field_names) {
     // Check if it's a state field
     if (to_states_index_.count(name)) {
       size_t idx = to_states_index_.at(name);
-      SLIC_ERROR_ROOT_IF(idx >= state_fields.size(),
-                         axom::fmt::format("State field index {} out of bounds (size={}) for field '{}'", idx,
-                                           state_fields.size(), name));
+      SLIC_ERROR_ROOT_IF(
+          idx >= state_fields.size(),
+          std::format("State field index {} out of bounds (size={}) for field '{}'", idx, state_fields.size(), name));
       fields_for_residual.push_back(state_fields[idx]);
     }
     // Otherwise check if it's a parameter
     else if (to_params_index_.count(name)) {
       size_t idx = to_params_index_.at(name);
       SLIC_ERROR_ROOT_IF(idx >= param_fields.size(),
-                         axom::fmt::format("Parameter field index {} out of bounds (size={}) for field '{}'", idx,
-                                           param_fields.size(), name));
+                         std::format("Parameter field index {} out of bounds (size={}) for field '{}'", idx,
+                                     param_fields.size(), name));
       fields_for_residual.push_back(param_fields[idx]);
     } else {
-      SLIC_ERROR_ROOT(axom::fmt::format("Field or parameter '{}' (required by weak form '{}') not found in FieldStore",
-                                        name, weak_form_name));
+      SLIC_ERROR_ROOT(std::format("Field or parameter '{}' (required by weak form '{}') not found in FieldStore", name,
+                                  weak_form_name));
     }
   }
   return fields_for_residual;
@@ -359,26 +321,27 @@ void FieldStore::markWeakFormInternal(const std::string& weak_form_name)
   internal_weak_forms_.insert(weak_form_name);
 }
 
-void FieldStore::addWeakFormReaction(std::string weak_form_name, std::string field_name)
+void FieldStore::addWeakFormFieldNames(std::string weak_form_name, UnknownAndTestFieldNames field_names)
 {
-  for (auto& kv : weak_form_to_test_field_) {
+  for (auto& kv : weak_form_field_names_) {
     if (kv.first == weak_form_name) {
-      kv.second = field_name;
+      kv.second = std::move(field_names);
       return;
     }
   }
-  weak_form_to_test_field_.push_back({weak_form_name, field_name});
+  weak_form_field_names_.push_back({std::move(weak_form_name), std::move(field_names)});
 }
 
-std::string FieldStore::getWeakFormReaction(const std::string& weak_form_name) const
+const UnknownAndTestFieldNames& FieldStore::getWeakFormFieldNames(const std::string& weak_form_name) const
 {
-  for (const auto& kv : weak_form_to_test_field_) {
+  for (const auto& kv : weak_form_field_names_) {
     if (kv.first == weak_form_name) {
       return kv.second;
     }
   }
-  SLIC_ERROR("Reaction field not found for weak form " << weak_form_name);
-  return "";
+  SLIC_ERROR("Field names not found for weak form " << weak_form_name);
+  static const UnknownAndTestFieldNames empty;
+  return empty;
 }
 
 const std::vector<FieldState>& FieldStore::getParameterFields() const { return params_; }
@@ -405,12 +368,12 @@ std::vector<FieldState> FieldStore::getOutputFieldStates() const
 std::vector<ReactionInfo> FieldStore::getReactionInfos() const
 {
   std::vector<ReactionInfo> infos;
-  for (const auto& kv : weak_form_to_test_field_) {
+  for (const auto& kv : weak_form_field_names_) {
     const std::string& weak_form_name = kv.first;
     if (internal_weak_forms_.count(weak_form_name)) {
       continue;
     }
-    const std::string& field_name = kv.second;
+    const std::string& field_name = kv.second.test;
     infos.push_back({weak_form_name, &getField(field_name).get()->space()});
   }
   return infos;