Tupek/field store interface extended

src/smith/differentiable_numerics/CMakeLists.txt

@@ -6,17 +6,20 @@
 
 set(differentiable_numerics_sources
     field_state.cpp
+    field_store.cpp
     differentiable_physics.cpp
     lumped_mass_explicit_newmark_state_advancer.cpp
-    solid_mechanics_state_advancer.cpp
+    solid_mechanics_time_integrator.cpp
     differentiable_solver.cpp
     nonlinear_solve.cpp
     evaluate_objective.cpp
     dirichlet_boundary_conditions.cpp
+    multiphysics_time_integrator.cpp
     )
 
 set(differentiable_numerics_headers
     field_state.hpp
+    field_store.hpp
     state_advancer.hpp
     reaction.hpp
     differentiable_solver.hpp
@@ -25,15 +28,17 @@ set(differentiable_numerics_headers
     explicit_dynamic_solve.hpp
     lumped_mass_weak_form.hpp
     lumped_mass_explicit_newmark_state_advancer.hpp
-    solid_mechanics_state_advancer.hpp
+    solid_mechanics_time_integrator.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
+    multiphysics_time_integrator.hpp
+    solid_dynamics_system.hpp
+    thermo_mechanics_system.hpp
+    system_base.hpp
     )
 
 set(differentiable_numerics_depends

src/smith/differentiable_numerics/differentiable_physics.cpp

@@ -120,12 +120,14 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& dual_nam
       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_states_.empty() && !reaction_names_.empty(),
+                "Reactions were not computed during advanceState, but were requested.");
+
   SLIC_ERROR_IF(
-      reaction_index >= reaction_names_.size(),
+      reaction_index >= reaction_states_.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();
 }
 
@@ -225,7 +227,10 @@ void DifferentiablePhysics::advanceTimestep(double dt)
   dt_prev_ = dt;
 
   TimeInfo time_info(time_, dt, static_cast<size_t>(cycle_));
-  field_states_ = advancer_->advanceState(time_info, *field_shape_displacement_, field_states_, field_params_);
+  auto [states, reactions] =
+      advancer_->advanceState(time_info, *field_shape_displacement_, field_states_, field_params_);
+  field_states_ = states;
+  reaction_states_ = reactions;
 
   cycle_++;
   time_ += dt;

src/smith/differentiable_numerics/differentiable_solver.cpp

@@ -232,6 +232,13 @@ std::vector<DifferentiableBlockSolver::FieldPtr> LinearDifferentiableBlockSolver
         u_bars[static_cast<size_t>(row_i)]->space(), "u_dual_" + std::to_string(row_i));
   }
 
+  // If it's a 1x1 system, pass the operator directly to avoid potential BlockOperator issues with smoothers
+  if (num_rows == 1) {
+    mfem_solver->SetOperator(*jacobian_transposed[0][0]);
+    mfem_solver->Mult(*u_bars[0], *u_duals[0]);
+    return u_duals;
+  }
+
   mfem::Array<int> block_offsets;
   block_offsets.SetSize(num_rows + 1);
   block_offsets[0] = 0;
@@ -303,7 +310,7 @@ std::vector<DifferentiableBlockSolver::FieldPtr> NonlinearDifferentiableBlockSol
 
   auto residual_op_ = std::make_unique<mfem_ext::StdFunctionOperator>(
       block_u->Size(),
-      [&residual_funcs, num_rows, &u_guesses, &block_r](const mfem::Vector& u_, mfem::Vector& r_) {
+      [&residual_funcs, num_rows, &u_guesses, &block_r, &block_offsets](const mfem::Vector& u_, mfem::Vector& r_) {
         const mfem::BlockVector* u = dynamic_cast<const mfem::BlockVector*>(&u_);
         SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a block vector");
         for (int row_i = 0; row_i < num_rows; ++row_i) {
@@ -319,7 +326,11 @@ std::vector<DifferentiableBlockSolver::FieldPtr> NonlinearDifferentiableBlockSol
       },
       [this, &block_offsets, &u_guesses, jacobian_funcs, num_rows](const mfem::Vector& u_) -> mfem::Operator& {
         const mfem::BlockVector* u = dynamic_cast<const mfem::BlockVector*>(&u_);
-        SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a block vector");
+        mfem::BlockVector u_block_wrapper;
+        if (!u) {
+          u_block_wrapper.Update(const_cast<double*>(u_.GetData()), block_offsets);
+          u = &u_block_wrapper;
+        }
         for (int row_i = 0; row_i < num_rows; ++row_i) {
           *u_guesses[static_cast<size_t>(row_i)] = u->GetBlock(row_i);
         }
@@ -359,6 +370,15 @@ std::vector<DifferentiableBlockSolver::FieldPtr> NonlinearDifferentiableBlockSol
         u_bars[static_cast<size_t>(row_i)]->space(), "u_dual_" + std::to_string(row_i));
   }
 
+  auto& linear_solver = nonlinear_solver_->linearSolver();
+
+  // If it's a 1x1 system, pass the operator directly to avoid potential BlockOperator issues with smoothers
+  if (num_rows == 1) {
+    linear_solver.SetOperator(*jacobian_transposed[0][0]);
+    linear_solver.Mult(*u_bars[0], *u_duals[0]);
+    return u_duals;
+  }
+
   mfem::Array<int> block_offsets;
   block_offsets.SetSize(num_rows + 1);
   block_offsets[0] = 0;
@@ -382,7 +402,6 @@ std::vector<DifferentiableBlockSolver::FieldPtr> NonlinearDifferentiableBlockSol
     }
   }
 
-  auto& linear_solver = nonlinear_solver_->linearSolver();
   linear_solver.SetOperator(*block_jac);
   linear_solver.Mult(*block_r, *block_ds);
 

src/smith/differentiable_numerics/differentiable_test_utils.hpp

@@ -15,9 +15,28 @@
 #include "gretl/double_state.hpp"
 #include "smith/differentiable_numerics/field_state.hpp"
 #include "smith/physics/scalar_objective.hpp"
+#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp"
 
 namespace smith {
 
+/**
+ * @brief Verify that reaction forces are zero at non-Dirichlet DOFs.
+ * @param reaction The reaction field to check.
+ * @param bc_manager Boundary condition manager to identify Dirichlet DOFs.
+ * @param tolerance Absolute tolerance for zero check.
+ */
+inline void checkUnconstrainedReactions(const FiniteElementDual& reaction, const BoundaryConditionManager& bc_manager,
+                                        double tolerance = 1e-8)
+{
+  FiniteElementState unconstrained_reactions(reaction.space(), "unconstrained_reactions");
+  unconstrained_reactions = reaction;
+  unconstrained_reactions.SetSubVector(bc_manager.allEssentialTrueDofs(), 0.0);
+
+  double max_unconstrained = unconstrained_reactions.Normlinf();
+  EXPECT_LT(max_unconstrained, tolerance)
+      << "Reaction forces should be zero at non-Dirichlet DOFs. Max violation: " << max_unconstrained;
+}
+
 /// @brief Utility function to construct a smith::functional which evaluates the total kinetic energy
 template <typename DispSpace, typename DensitySpace>
 auto createKineticEnergyIntegrator(smith::Domain& domain, const mfem::ParFiniteElementSpace& velocity_space,