Gretl solid mechanics #1501

tupek2 · 2025-12-05T18:57:03Z

Add implicit (differentiable) linear and nonlinear solvers. Add time discretization for second order systems. Add the ability to compute resultant forces. Dynamic patch test and design sensitivity tests are included.

…solid mechanics.

…ction abstraction.

…t residual.

…lid_mechanics.

cmake/thirdparty/SetupSmithThirdParty.cmake

white238 · 2026-01-06T00:17:58Z

src/smith/differentiable_numerics/differentiable_solver.hpp

+/// adjoint solves
+class LinearDifferentiableSolver : public DifferentiableSolver {
+ public:
+  /// @brief Construct from a linear solver and linear precondition which may be used by the linear solver


This sentence needs more linear!

white238 · 2026-01-06T00:18:40Z

src/smith/differentiable_numerics/differentiable_solver.hpp

+      nonlinear_solver_;  ///< the nonlinear equation solver used for the forward pass
+};
+
+/// @brief Abstract interface to DifferentiableBlockSolver inteface. Each dfferenriable block solve should provide


Suggested change

/// @brief Abstract interface to DifferentiableBlockSolver inteface. Each dfferenriable block solve should provide

/// @brief Abstract interface to DifferentiableBlockSolver interface. Each differentiable block solve should provide

white238 · 2026-01-06T00:18:52Z

src/smith/differentiable_numerics/differentiable_solver.hpp

+  using DualPtr = std::shared_ptr<FieldD>;                  ///< using
+  using MatrixPtr = std::unique_ptr<mfem::HypreParMatrix>;  ///< using
+
+  /// @brief required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space


Suggested change

/// @brief required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space

/// @brief Required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space

white238 · 2026-01-06T00:19:46Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.cpp

+{
+}
+
+}  // namespace smith


Suggested change

} // namespace smith

} // namespace smith

I know this is your favorite thing i nitpick....

white238 · 2026-01-06T00:21:05Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp

+  DirichletBoundaryConditions(const Mesh& mesh, mfem::ParFiniteElementSpace& space);
+
+  /// @brief Specify time and space varying Dirichlet boundary conditions over a domain.
+  /// @param domain domain All dofs in this domain have boundary conditions applied to it.


Suggested change

/// @param domain domain All dofs in this domain have boundary conditions applied to it.

/// @param domain All dofs in this domain have boundary conditions applied to it.

pretty sure thats not supposed to be there.

white238 · 2026-01-06T00:24:12Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp

+  }
+
+  /// @brief Specify time and space varying Dirichlet boundary conditions over a domain.
+  /// @param domain domain All dofs in this domain have boundary conditions applied to it.


Suggested change

/// @param domain domain All dofs in this domain have boundary conditions applied to it.

/// @param domain All dofs in this domain have boundary conditions applied to it.

white238 · 2026-01-06T00:25:03Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp

+  mfem::ParFiniteElementSpace& space_;   ///< save the space for the field which will be constrained
+};
+
+}  // namespace smith


Suggested change

} // namespace smith

} // namespace smith

white238 · 2026-01-06T00:28:00Z

src/smith/differentiable_numerics/field_state.hpp

-/// @brief initialize on the gretl::DataStore a ResultantState with values from s
-inline ResultantState create_field_resultant(gretl::DataStore& dataStore, const smith::FEDualPtr& s)
+/// @brief initialize on the gretl::DataStore a ReactionState with values from s
+inline ReactionState create_reaction_state(gretl::DataStore& dataStore, const smith::FEDualPtr& s)


As much as i like pot-hole case, Smith does camel case.

white238 · 2026-01-06T00:28:52Z

src/smith/differentiable_numerics/nonlinear_solve.cpp

@@ -0,0 +1,574 @@
+#include "smith/physics/weak_form.hpp"


Suggested change

#include "smith/physics/weak_form.hpp"

// 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 "smith/physics/weak_form.hpp"

white238 · 2026-01-06T00:30:55Z

src/smith/differentiable_numerics/nonlinear_solve.cpp

+
+    auto eval_residuals = [=](const std::vector<FEFieldPtr>& unknowns) {
+      SLIC_ERROR_IF(unknowns.size() != num_rows,
+                    "block solver unknown size must match the number or residuals in block_solve");


Suggested change

"block solver unknown size must match the number or residuals in block_solve");

"block solver unknowns size must match the number or residuals in block_solve");

.. i think...

white238 · 2026-01-06T00:32:01Z

src/smith/differentiable_numerics/nonlinear_solve.hpp

@@ -0,0 +1,66 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and


Suggested change

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and

// Copyright (c) Lawrence Livermore National Security, LLC and

white238 · 2026-01-06T00:35:15Z

src/smith/differentiable_numerics/reaction.hpp

@@ -0,0 +1,89 @@
+// Copyright (c), Lawrence Livermore National Security, LLC and


Suggested change

// Copyright (c), Lawrence Livermore National Security, LLC and

// Copyright (c) Lawrence Livermore National Security, LLC and

white238 · 2026-01-06T00:35:40Z

src/smith/differentiable_numerics/reaction.hpp

+
+/// @brief gretl-function implementation which evaluates the residual force (which is minus the mechanical force)
+/// given
+/// shape displacement, states and params.  The field_for_residual_space Field is only used to set the approriate size


Suggested change

/// shape displacement, states and params. The field_for_residual_space Field is only used to set the approriate size

/// shape displacement, states and params. The field_for_residual_space Field is only used to set the appropriate size

white238 · 2026-01-06T00:37:21Z

src/smith/differentiable_numerics/reaction.hpp

+    FEDualPtr R = std::make_shared<FiniteElementDual>(field_for_residual_space_->space(),
+                                                      "residual");  // set up output pointer
+    // evaluate the residual with zero acceleration contribution
+    // std::cout << "time info = " << time_info.time() << std::endl;


Delete dead code... does this comment go with the dead code or line below?

white238 · 2026-01-06T00:37:51Z

src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp

@@ -0,0 +1,62 @@
+#include "smith/physics/weak_form.hpp"


Suggested change

#include "smith/physics/weak_form.hpp"

// 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 "smith/physics/weak_form.hpp"

white238 · 2026-01-06T00:38:12Z

src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp

@@ -0,0 +1,118 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and


Suggested change

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and

// Copyright (c) Lawrence Livermore National Security, LLC and

white238 · 2026-01-06T00:38:48Z

src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp

+ * @file solid_mechanics_state_advancer.hpp
+ * .hpp
+ *
+ * @brief Specifies parametrized residuals and various linearized evaluations for arbitrary nonlinear systems of


Suggested change

* @brief Specifies parametrized residuals and various linearized evaluations for arbitrary nonlinear systems of

* @brief Specifies parameterized residuals and various linearized evaluations for arbitrary nonlinear systems of

white238 · 2026-01-06T00:40:00Z

src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp

+  }
+
+  /// @brief Utility function to consistently construct all the weak forms and FieldStates for a solid mechanics
+  /// application You will get back: shape_disp, states, params, time, and solid_mechanics_weak_form


Suggested change

/// application You will get back: shape_disp, states, params, time, and solid_mechanics_weak_form

/// application you will get back: shape_disp, states, params, time, and solid_mechanics_weak_form

white238 · 2026-01-06T00:40:50Z

src/smith/differentiable_numerics/time_discretized_weak_form.hpp

@@ -0,0 +1,154 @@
+// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and


Suggested change

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and

// Copyright (c) Lawrence Livermore National Security, LLC and

white238 · 2026-01-06T00:43:35Z

src/smith/differentiable_numerics/time_integration_rule.hpp

+                                    [[maybe_unused]] const T4& accel_old) const
+  {
+    return (2.0 / t.dt()) * (field_new - field_old) - velo_old;
+    // if (method_ == SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK) {


white238 · 2026-01-06T00:43:45Z

src/smith/differentiable_numerics/time_integration_rule.hpp

+  {
+    auto dt = t.dt();
+    return (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old;
+    // auto accel = (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old;


white238 · 2026-01-06T00:45:23Z

src/smith/physics/common.hpp

  size_t cycle() const { return cycle_; }

+  /// @brief return the time info corresponding to the end of this cycle (time + dt, dt, cycle)
+  TimeInfo end_time_info() const { return TimeInfo(time() + dt(), dt(), cycle()); }


white238 · 2026-01-06T00:47:49Z

The moving of smith/differentiable_numerics/test/differentiable_test_utils.hpp to smith/differentiable_numerics/differentiable_test_utils.hpp confuses me. Should it not be named that?

btalamini · 2026-01-05T18:48:29Z

src/smith/differentiable_numerics/differentiable_physics.cpp

+  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.");


The second error check appears unnecessary to me: the constructor generates the indices in 1:1 correspondence with the names. Also, I don't see a method called initializationStep in this class or the base class. The first error check is not strictly needed either, since the at method raises an exception if the key is not found.

More generally, I prefer having the data validity checking on the data structure itself to keep these kinds of checks from being duplicated at every use.

btalamini · 2026-01-05T18:56:52Z

src/smith/differentiable_numerics/differentiable_physics.cpp

+    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()));


This is the error check duplication I was referring to. There are several more.

btalamini · 2026-01-05T19:13:31Z

src/smith/differentiable_numerics/differentiable_solver.hpp

+  /// @brief destructor
+  virtual ~DifferentiableSolver() {}
+
+  /// @brief required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space


Suggested change

/// @brief required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space

/// @brief required for certain solvers/preconditioners, e.g. when multigrid algorithms want a near null-space

btalamini · 2026-01-05T21:46:10Z

src/smith/differentiable_numerics/nonlinear_solve.cpp

+  // there should be one less input state, as the higher time derivative term (e.g., acceleration), does not have a
+  // predictor


Is this comment up to date? The check below ensures things have the same size, not "one less" as the comment states. If the comment is correct, then I don't understand what it's trying to say. Can you rephrase?

btalamini · 2026-01-05T22:22:55Z

src/smith/differentiable_numerics/differentiable_solver.hpp

+      std::function<std::unique_ptr<mfem::HypreParMatrix>(const smith::FiniteElementState&)> jacobian) const override;
+
+  /// @overload
+  std::shared_ptr<smith::FiniteElementState> solveAdjoint(


This comment is not a point requiring changes in the PR, just a conversation I want to start while it's in mind: Now that we have an abstraction for a solver and adjoint solver, it could be a natural way to remove the memory-intensive copy of the sparse jacobian (for the sparse jacobian transpose). The solveAdjoint method could take the jacobian and use TransposeOperator or something similar to get the action of the transpose. We've been noting and complaining about this memory usage for a while. Should I open an issue or this for follow-up?

I like the idea of an abstraction like adjointSolve taking the responsibility of transposing (instead of putting it on the caller), since that's the main difference from the solve method.

btalamini · 2026-01-05T22:41:25Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp

+  /// smith::tensor<double,spatial_dim> corresponding to the spatial coordinate.  The functor must return a double.  For
+  /// example: [](double t, smith::tensor<double, dim> X) { return 1.0; }
+  template <int spatial_dim, typename AppliedDisplacementFunction>
+  void setScalarBCs(const Domain& domain, int component, AppliedDisplacementFunction applied_displacement)


I find it confusing to have a method named setScalarBCs take an argument for a component. In my view, this is an overload for setVectorBCs for the special case of setting one component, saving the typing of { and } characters around the index.

btalamini · 2026-01-05T22:48:10Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp

+  {
+    int field_dim = space_.GetVDim();
+    for (auto component : components) {
+      SLIC_ERROR_IF(component >= field_dim,


You should probably also check if component < 0.

btalamini · 2026-01-05T22:48:57Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp

It's great that this is factored out for all the physics modules! 👍

btalamini · 2026-01-05T22:55:37Z

src/smith/differentiable_numerics/evaluate_objective.cpp

+DoubleState evaluateObjective(const ScalarObjective& objective, const FieldState& shape_disp,
+                              const std::vector<FieldState>& inputs, const TimeInfo& time_info)
+{
+  const ScalarObjective* objective_ptr = &objective;


Why the indirection? You can't capture objective by reference in the lambdas? I'm not worried about overhead, just trying to clear up my understanding of gretl.

btalamini · 2026-01-06T00:57:08Z

src/smith/differentiable_numerics/nonlinear_solve.cpp

+  for (auto& p : params) allFields.push_back(p);
+  allFields.push_back(shape_disp);
+
+  bool have_bc_field = bc_field;


Just a thought: in the case where BCs are done the classic way, could you use the bc_manager to populate bc_field? Then the code would be unified, no conditionals for have_bc_field.

white238 · 2026-01-06T01:15:14Z

Sorry! wrong button!

lihanyu97 · 2026-01-06T16:00:39Z

src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp

+    }
+    setVectorBCs<spatial_dim>(domain, components, applied_displacement);
+  }
+


For adding dirichlet bc on all components, should we use
bcs_.addEssential(bc_set, coefficient, space)?

lihanyu97 · 2026-01-06T16:16:36Z

src/smith/differentiable_numerics/nonlinear_solve.cpp

+  auto constrained_dofs = bc_manager->allEssentialTrueDofs();
+  if (bc_field_ptr) {
+    for (int i = 0; i < constrained_dofs.Size(); i++) {
+      int j = constrained_dofs[i];
+      (*primal_field)[j] = (*bc_field_ptr)(j);
+    }
+  } else {


Suggested change

auto constrained_dofs = bc_manager->allEssentialTrueDofs();

if (bc_field_ptr) {

for (int i = 0; i < constrained_dofs.Size(); i++) {

int j = constrained_dofs[i];

(*primal_field)[j] = (*bc_field_ptr)(j);

}

} else {

if (bc_field_ptr) {

auto constrained_dofs = bc_manager->allEssentialTrueDofs();

for (int i = 0; i < constrained_dofs.Size(); i++) {

int j = constrained_dofs[i];

(*primal_field)[j] = (*bc_field_ptr)(j);

}

} else {

ebchin · 2026-01-06T22:05:50Z

src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp

+  }
+}
+
+TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics)


Should this be an example instead?

tupek2 added 28 commits November 21, 2025 13:00

Start bringing in additional differentiable utilities, start testing …

00c35f5

…solid mechanics.

Fix equation solver issue, change problem to 2D.

27914af

Fix some print level in trust region solver.

3ccf6f3

Gradient tests in place.

f9fda3c

Start implementing a solid mechanics interface, implement initial rea…

3615eeb

…ction abstraction.

Working toward reactions. Get it compiling, need to test and implemen…

a02d15c

…t residual.

Add solid mechanics and solvers to differentiable physics things.

a5a8f4a

Merge branch 'gretl_submodule' into gretl_solid_mechanics

0bece90

update gretl with new build configs.

05511ff

Working out some design for reaction forces.

d6b2395

Get in an initial design for how to handle resultants/duals.

3a854e2

Merge branch 'develop' into gretl_solid_mechanics

728acbc

Some few comment.

2a05611

Working on testing reaction sensitivities through base-physics.

78d813a

Get more of residual test working.

b12287b

Implement some version of reaction sensitivities in differentiable_so…

103dbfe

…lid_mechanics.

Fix a doc warning.

15b79c3

Merge branch 'gretl_submodule' into gretl_solid_mechanics

1526218

Fix merge mistakes.

fa39a8e

Improve some docs.

86868b0

some more documentation.

f647e9f

More docs.

49d47d3

Continue documenting.

5cdbf01

Updating more docs.

4bad57c

Get all the documentation fixed up.

fcd7baf

Merge branch 'gretl_submodule' into gretl_solid_mechanics

4f3481c

Fixed hidden var, start trying to test dynamics in gravity.

18247a9

Trying to test dynamics.

e735038

white238 reviewed Dec 12, 2025

View reviewed changes

cmake/thirdparty/SetupSmithThirdParty.cmake Outdated Show resolved Hide resolved

Fixed style.

480a6e9

white238 reviewed Jan 6, 2026

View reviewed changes

Move test_utils out of test folder so cz tests can use checkGradWrt

ab05798

white238 reviewed Jan 6, 2026

View reviewed changes

white238 closed this Jan 6, 2026

btalamini approved these changes Jan 6, 2026

View reviewed changes

white238 reopened this Jan 6, 2026

lihanyu97 reviewed Jan 6, 2026

View reviewed changes

ebchin reviewed Jan 6, 2026

View reviewed changes

	/// @brief Abstract interface to DifferentiableBlockSolver inteface. Each dfferenriable block solve should provide
	/// @brief Abstract interface to DifferentiableBlockSolver interface. Each differentiable block solve should provide

	/// @brief required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space
	/// @brief Required for certain solvers/preconditions, e.g. when multigrid algorithms want a near null-space

	/// @param domain domain All dofs in this domain have boundary conditions applied to it.
	/// @param domain All dofs in this domain have boundary conditions applied to it.

-#include "smith/physics/weak_form.hpp"
+// 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 "smith/physics/weak_form.hpp"

	"block solver unknown size must match the number or residuals in block_solve");
	"block solver unknowns size must match the number or residuals in block_solve");

		@@ -0,0 +1,66 @@
		// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and

	// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
	// Copyright (c) Lawrence Livermore National Security, LLC and

		@@ -0,0 +1,89 @@
		// Copyright (c), Lawrence Livermore National Security, LLC and

	/// shape displacement, states and params. The field_for_residual_space Field is only used to set the approriate size
	/// shape displacement, states and params. The field_for_residual_space Field is only used to set the appropriate size

		@@ -0,0 +1,118 @@
		// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and

	* @brief Specifies parametrized residuals and various linearized evaluations for arbitrary nonlinear systems of
	* @brief Specifies parameterized residuals and various linearized evaluations for arbitrary nonlinear systems of

	/// application You will get back: shape_disp, states, params, time, and solid_mechanics_weak_form
	/// application you will get back: shape_disp, states, params, time, and solid_mechanics_weak_form

		@@ -0,0 +1,154 @@
		// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and

		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()));

		// there should be one less input state, as the higher time derivative term (e.g., acceleration), does not have a
		// predictor

Gretl solid mechanics #1501

Are you sure you want to change the base?

Gretl solid mechanics #1501

Uh oh!

Conversation

tupek2 commented Dec 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

white238 commented Jan 6, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

white238 commented Jan 6, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

tupek2 commented Dec 5, 2025 •

edited

Loading