From 18622066d9a049fb43ddec9d03e2b40fefee2ab7 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 30 Dec 2025 16:11:00 -0800 Subject: [PATCH 01/17] Adding a block nonlinear solver. --- .../differentiable_solver.cpp | 125 +++++++++++++++++- .../differentiable_solver.hpp | 30 +++++ 2 files changed, 150 insertions(+), 5 deletions(-) diff --git a/src/smith/differentiable_numerics/differentiable_solver.cpp b/src/smith/differentiable_numerics/differentiable_solver.cpp index 9cd3798776..3ce4860799 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.cpp +++ b/src/smith/differentiable_numerics/differentiable_solver.cpp @@ -141,8 +141,6 @@ std::shared_ptr NonlinearDifferentiableSolver::solve( nonlinear_solver_->setOperator(*residual_op_); nonlinear_solver_->solve(*u); - // std::cout << "solution norm = " << u->Norml2() << std::endl; - return u; } @@ -269,6 +267,116 @@ std::vector LinearDifferentiableBlockSolver return u_duals; } + +NonlinearDifferentiableBlockSolver::NonlinearDifferentiableBlockSolver(std::unique_ptr s) + : nonlinear_solver_(std::move(s)) +{ +} + +void NonlinearDifferentiableBlockSolver::completeSetup(const std::vector&) +{ + //initializeSolver(&nonlinear_solver_->preconditioner(), u); +} + +std::vector NonlinearDifferentiableBlockSolver::solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residual_funcs, + std::function>(const std::vector&)> jacobian_funcs) const +{ + SMITH_MARK_FUNCTION; + + int num_rows = static_cast(u_guesses.size()); + SLIC_ERROR_IF(num_rows < 0, "Number of residual rows must be non-negative"); + + mfem::Array block_offsets; + block_offsets.SetSize(num_rows + 1); + block_offsets[0] = 0; + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_offsets[row_i + 1] = u_guesses[static_cast(row_i)]->space().TrueVSize(); + } + block_offsets.PartialSum(); + + auto block_du = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_du->GetBlock(row_i) = *u_guesses[static_cast(row_i)]; + } + + auto residuals = residual_funcs(u_guesses); + auto block_r = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_r->GetBlock(row_i) = residuals[static_cast(row_i)]; + } + + auto jacs = jacobian_funcs(u_guesses); + auto block_jac = std::make_unique(block_offsets); + for (int i = 0; i < num_rows; ++i) { + for (int j = 0; j < num_rows; ++j) { + block_jac->SetBlock(i, j, jacs[static_cast(i)][static_cast(j)].get()); + } + } + + auto& linear_solver = nonlinear_solver_->linearSolver(); + linear_solver.SetOperator(*block_jac); + linear_solver.Mult(*block_r, *block_du); + + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_guesses[static_cast(row_i)] -= block_du->GetBlock(row_i); + } + *block_du = 0.0; + + return u_guesses; +} + +std::vector NonlinearDifferentiableBlockSolver::solveAdjoint( + const std::vector& u_bars, std::vector>& jacobian_transposed) const +{ + SMITH_MARK_FUNCTION; + + int num_rows = static_cast(u_bars.size()); + SLIC_ERROR_IF(num_rows < 0, "Number of residual rows must be non-negative"); + + std::vector u_duals(static_cast(num_rows)); + for (int row_i = 0; row_i < num_rows; ++row_i) { + u_duals[static_cast(row_i)] = std::make_shared( + u_bars[static_cast(row_i)]->space(), "u_dual_" + std::to_string(row_i)); + } + + mfem::Array block_offsets; + block_offsets.SetSize(num_rows + 1); + block_offsets[0] = 0; + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_offsets[row_i + 1] = u_bars[static_cast(row_i)]->space().TrueVSize(); + } + block_offsets.PartialSum(); + + auto block_ds = std::make_unique(block_offsets); + *block_ds = 0.0; + + auto block_r = std::make_unique(block_offsets); + for (int row_i = 0; row_i < num_rows; ++row_i) { + block_r->GetBlock(row_i) = *u_bars[static_cast(row_i)]; + } + + auto block_jac = std::make_unique(block_offsets); + for (int i = 0; i < num_rows; ++i) { + for (int j = 0; j < num_rows; ++j) { + block_jac->SetBlock(i, j, jacobian_transposed[static_cast(i)][static_cast(j)].get()); + } + } + + auto& linear_solver = nonlinear_solver_->linearSolver(); + linear_solver.SetOperator(*block_jac); + linear_solver.Mult(*block_r, *block_ds); + + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_duals[static_cast(row_i)] = block_ds->GetBlock(row_i); + } + + return u_duals; +} + + + std::shared_ptr buildDifferentiableLinearSolver(LinearSolverOptions linear_opts, const smith::Mesh& mesh) { @@ -277,10 +385,17 @@ std::shared_ptr buildDifferentiableLinearSolver(Line } std::shared_ptr buildDifferentiableNonlinearSolver( - smith::NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh) + NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh) +{ + auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); + return std::make_shared(std::move(solid_solver)); +} + +std::shared_ptr buildDifferentiableNonlinearBlockSolver( + NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh) { - auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); - return std::make_shared(std::move(solid_solver)); + auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); + return std::make_shared(std::move(solid_solver)); } } // namespace smith diff --git a/src/smith/differentiable_numerics/differentiable_solver.hpp b/src/smith/differentiable_numerics/differentiable_solver.hpp index a58540918a..bee488b5d9 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.hpp +++ b/src/smith/differentiable_numerics/differentiable_solver.hpp @@ -175,6 +175,30 @@ class LinearDifferentiableBlockSolver : public DifferentiableBlockSolver { mutable std::unique_ptr mfem_preconditioner; ///< stored mfem block preconditioner }; +/// @brief Implementation of the DifferentiableBlockSolver interface for the special case of nonlinear solves with linear +/// adjoint solves +class NonlinearDifferentiableBlockSolver : public DifferentiableBlockSolver { + public: + /// @brief Construct from a linear solver and linear block precondition which may be used by the linear solver + NonlinearDifferentiableBlockSolver(std::unique_ptr s); + + /// @overload + void completeSetup(const std::vector& us) override; + + /// @overload + std::vector solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residuals, + std::function>(const std::vector&)> jacobians) const override; + + /// @overload + std::vector solveAdjoint(const std::vector& u_bars, + std::vector>& jacobian_transposed) const override; + + mutable std::unique_ptr + nonlinear_solver_; ///< the nonlinear equation solver used for the forward pass +}; + /// @brief Create a differentiable linear solver /// @param linear_opts linear options struct /// @param mesh mesh @@ -189,4 +213,10 @@ std::shared_ptr buildDifferentiableNonlinearSolve LinearSolverOptions linear_opts, const smith::Mesh& mesh); +/// @brief Create a differentiable nonlinear block solver +/// @param nonlinear_opts nonlinear options struct +/// @param linear_opts linear options struct +/// @param mesh mesh +std::shared_ptr buildDifferentiableNonlinearBlockSolver(NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh); + } // namespace smith From 8299a9808f48f1c8a06692ae0071a33923ea8e64 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 8 Jan 2026 13:37:01 -0800 Subject: [PATCH 02/17] Fix style, adding differentiable nonlinear block solver. --- .../differentiable_solver.cpp | 69 ++++++++++++------- .../differentiable_solver.hpp | 11 ++- .../nonlinear_solve.cpp | 54 +++++++++------ .../nonlinear_solve.hpp | 5 +- 4 files changed, 90 insertions(+), 49 deletions(-) diff --git a/src/smith/differentiable_numerics/differentiable_solver.cpp b/src/smith/differentiable_numerics/differentiable_solver.cpp index 3ce4860799..61c984e85b 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.cpp +++ b/src/smith/differentiable_numerics/differentiable_solver.cpp @@ -267,7 +267,6 @@ std::vector LinearDifferentiableBlockSolver return u_duals; } - NonlinearDifferentiableBlockSolver::NonlinearDifferentiableBlockSolver(std::unique_ptr s) : nonlinear_solver_(std::move(s)) { @@ -275,7 +274,7 @@ NonlinearDifferentiableBlockSolver::NonlinearDifferentiableBlockSolver(std::uniq void NonlinearDifferentiableBlockSolver::completeSetup(const std::vector&) { - //initializeSolver(&nonlinear_solver_->preconditioner(), u); + // initializeSolver(&nonlinear_solver_->preconditioner(), u); } std::vector NonlinearDifferentiableBlockSolver::solve( @@ -296,33 +295,56 @@ std::vector NonlinearDifferentiableBlockSol } block_offsets.PartialSum(); - auto block_du = std::make_unique(block_offsets); + auto block_u = std::make_unique(block_offsets); for (int row_i = 0; row_i < num_rows; ++row_i) { - block_du->GetBlock(row_i) = *u_guesses[static_cast(row_i)]; + block_u->GetBlock(row_i) = *u_guesses[static_cast(row_i)]; } - auto residuals = residual_funcs(u_guesses); auto block_r = std::make_unique(block_offsets); - for (int row_i = 0; row_i < num_rows; ++row_i) { - block_r->GetBlock(row_i) = residuals[static_cast(row_i)]; - } - auto jacs = jacobian_funcs(u_guesses); - auto block_jac = std::make_unique(block_offsets); - for (int i = 0; i < num_rows; ++i) { - for (int j = 0; j < num_rows; ++j) { - block_jac->SetBlock(i, j, jacs[static_cast(i)][static_cast(j)].get()); - } - } + auto residual_op_ = std::make_unique( + block_u->Size(), + [residual_funcs, num_rows, &u_guesses, &block_r](const mfem::Vector& u_, mfem::Vector& r_) { + const mfem::BlockVector* u = dynamic_cast(&u_); + SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a blocl vector"); + for (int row_i = 0; row_i < num_rows; ++row_i) { + *u_guesses[static_cast(row_i)] = u->GetBlock(row_i); + } + auto residuals = residual_funcs(u_guesses); + // auto block_r = std::make_unique(block_offsets); + // auto block_r = dynamic_cast(&r_); + SLIC_ERROR_IF(!block_r, "Invalid r cast in block differentiable solver to a block vector"); + for (int row_i = 0; row_i < num_rows; ++row_i) { + auto r = residuals[static_cast(row_i)]; + block_r->GetBlock(row_i) = r; + } + r_ = *block_r; + }, + [this, &block_offsets, &u_guesses, jacobian_funcs, num_rows](const mfem::Vector& u_) -> mfem::Operator& { + const mfem::BlockVector* u = dynamic_cast(&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) { + *u_guesses[static_cast(row_i)] = u->GetBlock(row_i); + } + block_jac_ = std::make_unique(block_offsets); + matrix_of_jacs_ = jacobian_funcs(u_guesses); + for (int i = 0; i < num_rows; ++i) { + for (int j = 0; j < num_rows; ++j) { + auto& J = matrix_of_jacs_[static_cast(i)][static_cast(j)]; + if (J) { + block_jac_->SetBlock(i, j, J.get()); + } + } + } + return *block_jac_; + }); - auto& linear_solver = nonlinear_solver_->linearSolver(); - linear_solver.SetOperator(*block_jac); - linear_solver.Mult(*block_r, *block_du); + nonlinear_solver_->setOperator(*residual_op_); + nonlinear_solver_->solve(*block_u); for (int row_i = 0; row_i < num_rows; ++row_i) { - *u_guesses[static_cast(row_i)] -= block_du->GetBlock(row_i); + *u_guesses[static_cast(row_i)] = block_u->GetBlock(row_i); } - *block_du = 0.0; return u_guesses; } @@ -375,8 +397,6 @@ std::vector NonlinearDifferentiableBlockSol return u_duals; } - - std::shared_ptr buildDifferentiableLinearSolver(LinearSolverOptions linear_opts, const smith::Mesh& mesh) { @@ -384,8 +404,9 @@ std::shared_ptr buildDifferentiableLinearSolver(Line return std::make_shared(std::move(linear_solver), std::move(precond)); } -std::shared_ptr buildDifferentiableNonlinearSolver( - NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh) +std::shared_ptr buildDifferentiableNonlinearSolver(NonlinearSolverOptions nonlinear_opts, + LinearSolverOptions linear_opts, + const smith::Mesh& mesh) { auto solid_solver = std::make_unique(nonlinear_opts, linear_opts, mesh.getComm()); return std::make_shared(std::move(solid_solver)); diff --git a/src/smith/differentiable_numerics/differentiable_solver.hpp b/src/smith/differentiable_numerics/differentiable_solver.hpp index bee488b5d9..745e0509c2 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.hpp +++ b/src/smith/differentiable_numerics/differentiable_solver.hpp @@ -19,6 +19,7 @@ namespace mfem { class Solver; class Vector; class HypreParMatrix; +class BlockOperator; } // namespace mfem namespace smith { @@ -175,8 +176,8 @@ class LinearDifferentiableBlockSolver : public DifferentiableBlockSolver { mutable std::unique_ptr mfem_preconditioner; ///< stored mfem block preconditioner }; -/// @brief Implementation of the DifferentiableBlockSolver interface for the special case of nonlinear solves with linear -/// adjoint solves +/// @brief Implementation of the DifferentiableBlockSolver interface for the special case of nonlinear solves with +/// linear adjoint solves class NonlinearDifferentiableBlockSolver : public DifferentiableBlockSolver { public: /// @brief Construct from a linear solver and linear block precondition which may be used by the linear solver @@ -195,6 +196,9 @@ class NonlinearDifferentiableBlockSolver : public DifferentiableBlockSolver { std::vector solveAdjoint(const std::vector& u_bars, std::vector>& jacobian_transposed) const override; + mutable std::unique_ptr block_jac_; + mutable std::vector> matrix_of_jacs_; + mutable std::unique_ptr nonlinear_solver_; ///< the nonlinear equation solver used for the forward pass }; @@ -217,6 +221,7 @@ std::shared_ptr buildDifferentiableNonlinearSolve /// @param nonlinear_opts nonlinear options struct /// @param linear_opts linear options struct /// @param mesh mesh -std::shared_ptr buildDifferentiableNonlinearBlockSolver(NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh); +std::shared_ptr buildDifferentiableNonlinearBlockSolver( + NonlinearSolverOptions nonlinear_opts, LinearSolverOptions linear_opts, const smith::Mesh& mesh); } // namespace smith diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 7d2ed606a9..c73d0b8d84 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -274,7 +274,7 @@ std::vector block_solve(const std::vector& residual_evals const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const DifferentiableBlockSolver* solver, - const std::vector bc_managers) + const std::vector bc_managers) { SMITH_MARK_FUNCTION; size_t num_rows_ = residual_evals.size(); @@ -326,10 +326,10 @@ std::vector block_solve(const std::vector& residual_evals SMITH_MARK_BEGIN("solve forward"); const size_t num_rows = num_state_inputs.size(); std::vector> input_fields(num_rows); - SLIC_ERROR_IF(num_rows != num_param_inputs.size(), "row count for params and columns are inconsistent"); + SLIC_ERROR_IF(num_rows != num_param_inputs.size(), "row count for params and states are inconsistent"); - // The order of inputs in upstreams seems to be: - // states of residual 0 -> states of residual 1 -> params of residual 0 -> params of residual 1 + // The order of inputs in upstreams is: + // states of residual 0, states of residual 1, ... , params of residual 0, params of residual 1, ... size_t field_count = 0; for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t state_i = 0; state_i < num_state_inputs[row_i]; ++state_i) { @@ -344,12 +344,18 @@ std::vector block_solve(const std::vector& residual_evals std::vector diagonal_fields(num_rows); for (size_t row_i = 0; row_i < num_rows; ++row_i) { - diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][block_indices[row_i][row_i]]); + size_t prime_unknown_row_i = block_indices[row_i][row_i]; + SLIC_ERROR_IF(prime_unknown_row_i == invalid_block_index, + "The primary unknown field (field index for block_index[n][n], must not be invalid)"); + diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][prime_unknown_row_i]); } for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t col_j = 0; col_j < num_rows; ++col_j) { - input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + size_t prime_unknown_ij = block_indices[row_i][col_j]; + if (prime_unknown_ij != invalid_block_index) { + input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + } } } @@ -371,7 +377,6 @@ std::vector block_solve(const std::vector& residual_evals getConstFieldPointers(input_fields[row_i])); residuals[row_i].SetSubVector(bc_managers[row_i]->allEssentialTrueDofs(), 0.0); } - return residuals; }; @@ -391,32 +396,39 @@ std::vector block_solve(const std::vector& residual_evals std::vector tangent_weights(row_field_inputs.size(), 0.0); for (size_t col_j = 0; col_j < num_rows; ++col_j) { size_t field_index_to_diff = block_indices[row_i][col_j]; - tangent_weights[field_index_to_diff] = 1.0; - auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), - getConstFieldPointers(row_field_inputs), tangent_weights); - jacobians[row_i].emplace_back(std::move(jac_ij)); - tangent_weights[field_index_to_diff] = 0.0; + if (field_index_to_diff != invalid_block_index) { + tangent_weights[field_index_to_diff] = 1.0; + auto jac_ij = residual_evals[row_i]->jacobian(time_info, shape_disp_ptr.get(), + getConstFieldPointers(row_field_inputs), tangent_weights); + jacobians[row_i].emplace_back(std::move(jac_ij)); + tangent_weights[field_index_to_diff] = 0.0; + } else { + jacobians[row_i].emplace_back(nullptr); + } } } // Apply BCs to the block system for (size_t row_i = 0; row_i < num_rows; ++row_i) { - mfem::HypreParMatrix* Jii = - jacobians[row_i][row_i]->EliminateRowsCols(bc_managers[row_i]->allEssentialTrueDofs()); - delete Jii; + if (jacobians[row_i][row_i]) { + jacobians[row_i][row_i]->EliminateBC(bc_managers[row_i]->allEssentialTrueDofs(), + mfem::Operator::DiagonalPolicy::DIAG_ONE); + } for (size_t col_j = 0; col_j < num_rows; ++col_j) { if (col_j != row_i) { - jacobians[row_i][col_j]->EliminateRows(bc_managers[row_i]->allEssentialTrueDofs()); - mfem::HypreParMatrix* Jji = - jacobians[col_j][row_i]->EliminateCols(bc_managers[row_i]->allEssentialTrueDofs()); - delete Jji; + if (jacobians[row_i][col_j]) { + jacobians[row_i][col_j]->EliminateRows(bc_managers[row_i]->allEssentialTrueDofs()); + } + if (jacobians[col_j][row_i]) { + mfem::HypreParMatrix* Jji = + jacobians[col_j][row_i]->EliminateCols(bc_managers[row_i]->allEssentialTrueDofs()); + delete Jji; + } } } } - return jacobians; }; - diagonal_fields = solver->solve(diagonal_fields, eval_residuals, eval_jacobians); downstream.set, std::vector>(diagonal_fields); diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp index 07027b274e..df94c229f5 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.hpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -24,6 +24,9 @@ class DifferentiableBlockSolver; class BoundaryConditionManager; class DirichletBoundaryConditions; +/// @brief magic number for represending a field which is not an argument of the weak form. +static constexpr size_t invalid_block_index = std::numeric_limits::max() - 1; + /// @brief Solve a nonlinear system of equations as defined by the weak form, assuming that the field indexed by /// unknown_index is the unknown field /// @param residual_eval The weak form which defines the equations to be solved @@ -61,6 +64,6 @@ std::vector block_solve(const std::vector& residual_evals const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const DifferentiableBlockSolver* solver, - const std::vector bc_managers); + const std::vector bc_managers); } // namespace smith From a94e2fb2ddcf15e82e0e08e96d4febd8a698b08d Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 8 Jan 2026 13:42:47 -0800 Subject: [PATCH 03/17] Fix style. --- .../differentiable_numerics/differentiable_solver.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/smith/differentiable_numerics/differentiable_solver.hpp b/src/smith/differentiable_numerics/differentiable_solver.hpp index 745e0509c2..1eb2b877be 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.hpp +++ b/src/smith/differentiable_numerics/differentiable_solver.hpp @@ -196,8 +196,11 @@ class NonlinearDifferentiableBlockSolver : public DifferentiableBlockSolver { std::vector solveAdjoint(const std::vector& u_bars, std::vector>& jacobian_transposed) const override; - mutable std::unique_ptr block_jac_; - mutable std::vector> matrix_of_jacs_; + mutable std::unique_ptr + block_jac_; ///< Need to hold an instance of a block operator to work with the mfem solver interface + mutable std::vector> + matrix_of_jacs_; ///< Holding vectors of block matrices to that do not going out of scope before the mfem solver + ///< is done with using them in the block_jac_ mutable std::unique_ptr nonlinear_solver_; ///< the nonlinear equation solver used for the forward pass From 6aa3047fac45842457b2e2c31c5ba54b18bdad7a Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 8 Jan 2026 14:19:38 -0800 Subject: [PATCH 04/17] Remove boundary condition option that causes some template priority order issues. --- .../dirichlet_boundary_conditions.hpp | 22 ++----------------- 1 file changed, 2 insertions(+), 20 deletions(-) diff --git a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp index c148350df2..d79fc72b2f 100644 --- a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp +++ b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp @@ -72,39 +72,21 @@ class DirichletBoundaryConditions { /// @brief Specify time and space varying Dirichlet boundary conditions over a domain. /// @param domain All dofs in this domain have boundary conditions applied to it. - /// @param component component direction to apply boundary condition to if the underlying field is a vector-field. /// @param applied_displacement applied_displacement is a functor which takes time, and a /// smith::tensor corresponding to the spatial coordinate. The functor must return a double. For /// example: [](double t, smith::tensor X) { return 1.0; } template - void setVectorBCs(const Domain& domain, int component, AppliedDisplacementFunction applied_displacement) + void setScalarBCs(const Domain& domain, AppliedDisplacementFunction applied_displacement) { - const int field_dim = space_.GetVDim(); - SLIC_ERROR_IF(component >= field_dim, - axom::fmt::format("Trying to set boundary conditions on a field with dim {}, using component {}", - field_dim, component)); auto mfem_coefficient_function = [applied_displacement](const mfem::Vector& X_mfem, double t) { auto X = make_tensor([&X_mfem](int k) { return X_mfem[k]; }); return applied_displacement(t, X); }; auto dof_list = domain.dof_list(&space_); - // scalar ldofs -> vector ldofs - space_.DofsToVDofs(component, dof_list); auto component_disp_bdr_coef_ = std::make_shared(mfem_coefficient_function); - bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_, component); - } - - /// @brief Specify time and space varying Dirichlet boundary conditions over a domain. - /// @param domain All dofs in this domain have boundary conditions applied to it. - /// @param applied_displacement applied_displacement is a functor which takes time, and a - /// smith::tensor corresponding to the spatial coordinate. The functor must return a double. For - /// example: [](double t, smith::tensor X) { return 1.0; } - template - void setScalarBCs(const Domain& domain, AppliedDisplacementFunction applied_displacement) - { - setScalarBCs(domain, 0, applied_displacement); + bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_); } /// @brief Constrain the dofs of a scalar field over a domain From 218b75ec8dadbb55c4b94a31b4b935657a3f1259 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 12 Jan 2026 14:10:10 -0800 Subject: [PATCH 05/17] Still working on block nonlinear solver. --- .../differentiable_numerics/differentiable_solver.cpp | 3 +-- .../dirichlet_boundary_conditions.hpp | 3 ++- src/smith/differentiable_numerics/nonlinear_solve.cpp | 8 ++++++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/smith/differentiable_numerics/differentiable_solver.cpp b/src/smith/differentiable_numerics/differentiable_solver.cpp index ebcd7941a5..7a3b256002 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.cpp +++ b/src/smith/differentiable_numerics/differentiable_solver.cpp @@ -303,7 +303,7 @@ std::vector NonlinearDifferentiableBlockSol auto residual_op_ = std::make_unique( 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](const mfem::Vector& u_, mfem::Vector& r_) { const mfem::BlockVector* u = dynamic_cast(&u_); SLIC_ERROR_IF(!u, "Invalid u cast in block differentiable solver to a blocl vector"); for (int row_i = 0; row_i < num_rows; ++row_i) { @@ -337,7 +337,6 @@ std::vector NonlinearDifferentiableBlockSol } return *block_jac_; }); - nonlinear_solver_->setOperator(*residual_op_); nonlinear_solver_->solve(*block_u); diff --git a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp index d79fc72b2f..d7ddec9a55 100644 --- a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp +++ b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp @@ -84,9 +84,10 @@ class DirichletBoundaryConditions { }; auto dof_list = domain.dof_list(&space_); + space_.DofsToVDofs(static_cast(0), dof_list); auto component_disp_bdr_coef_ = std::make_shared(mfem_coefficient_function); - bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_); + bcs_.addEssential(dof_list, component_disp_bdr_coef_, space_, 0); } /// @brief Constrain the dofs of a scalar field over a domain diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 4bf4643d88..d28af788e5 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -305,7 +305,6 @@ std::vector block_solve(const std::vector& residual_evals } } allFields.push_back(shape_disp); - struct ZeroDualVectors { std::vector operator()(const std::vector& fs) { @@ -319,7 +318,6 @@ std::vector block_solve(const std::vector& residual_evals FieldVecState sol = shape_disp.create_state, std::vector>(allFields, ZeroDualVectors()); - sol.set_eval([=](const gretl::UpstreamStates& upstreams, gretl::DownstreamState& downstream) { SMITH_MARK_BEGIN("solve forward"); const size_t num_rows = num_state_inputs.size(); @@ -371,6 +369,11 @@ std::vector block_solve(const std::vector& residual_evals } for (size_t row_i = 0; row_i < num_rows; ++row_i) { + //std::cout << "r " << row_i << " inputs = "; + //for (auto& f : input_fields[row_i]) { + // std::cout << f->name() << " " << f->Size() << " "; + //} + //std::cout << std::endl; residuals[row_i] = residual_evals[row_i]->residual(time_info, shape_disp_ptr.get(), getConstFieldPointers(input_fields[row_i])); residuals[row_i].SetSubVector(bc_managers[row_i]->allEssentialTrueDofs(), 0.0); @@ -427,6 +430,7 @@ std::vector block_solve(const std::vector& residual_evals } return jacobians; }; + diagonal_fields = solver->solve(diagonal_fields, eval_residuals, eval_jacobians); downstream.set, std::vector>(diagonal_fields); From b3e83493513aba7ceba83b6c8ea1a953fc50bda9 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Fri, 16 Jan 2026 09:23:41 -0800 Subject: [PATCH 06/17] SOme minor updates. --- .../nonlinear_solve.cpp | 17 ++++++----- .../solid_mechanics_state_advancer.cpp | 4 +-- .../time_discretized_weak_form.hpp | 29 ++++++++++--------- 3 files changed, 27 insertions(+), 23 deletions(-) diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index d28af788e5..1f1f34a558 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -346,6 +346,11 @@ std::vector block_solve(const std::vector& residual_evals diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][prime_unknown_row_i]); } + for (size_t row_i = 0; row_i < num_rows; ++row_i) { + FEFieldPtr primal_field_row_i = diagonal_fields[row_i]; + applyBoundaryConditions(time_info.time(), bc_managers[row_i], primal_field_row_i, nullptr); + } + for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t col_j = 0; col_j < num_rows; ++col_j) { size_t prime_unknown_ij = block_indices[row_i][col_j]; @@ -367,13 +372,7 @@ std::vector block_solve(const std::vector& residual_evals *primal_field_row_i = *unknowns[row_i]; applyBoundaryConditions(time_info.time(), bc_managers[row_i], primal_field_row_i, nullptr); } - for (size_t row_i = 0; row_i < num_rows; ++row_i) { - //std::cout << "r " << row_i << " inputs = "; - //for (auto& f : input_fields[row_i]) { - // std::cout << f->name() << " " << f->Size() << " "; - //} - //std::cout << std::endl; residuals[row_i] = residual_evals[row_i]->residual(time_info, shape_disp_ptr.get(), getConstFieldPointers(input_fields[row_i])); residuals[row_i].SetSubVector(bc_managers[row_i]->allEssentialTrueDofs(), 0.0); @@ -574,7 +573,11 @@ std::vector block_solve(const std::vector& residual_evals for (size_t i = 0; i < num_rows_; ++i) { FieldState s = gretl::create_state( zero_dual_from_state(), - [i](const std::vector& sols) { return std::make_shared(*sols[i]); }, + [i](const std::vector& sols) { + auto state_copy = std::make_shared(sols[i]->space(), sols[i]->name()); + *state_copy = *sols[i]; + return state_copy; + }, [i](const std::vector&, const FEFieldPtr&, std::vector& sols_, const FEDualPtr& output_) { *sols_[i] += *output_; }, sol); diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp index 7a9b888e4d..d60e9dfd1e 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp @@ -31,7 +31,7 @@ std::vector SolidMechanicsStateAdvancer::advanceState(const TimeInfo { std::vector states_old = states_old_; if (time_info.cycle() == 0) { - states_old[ACCELERATION] = solve(*solid_dynamic_weak_forms_->quasi_static_weak_form, shape_disp, states_old_, + states_old[ACCELERATION] = solve(*solid_dynamic_weak_forms_->final_reaction_weak_form, shape_disp, states_old_, params, time_info, *solver_, *vector_bcs_, ACCELERATION); } @@ -61,7 +61,7 @@ std::vector SolidMechanicsStateAdvancer::computeReactions(const T { std::vector solid_inputs{states[DISPLACEMENT], states[VELOCITY], states[ACCELERATION]}; solid_inputs.insert(solid_inputs.end(), params.begin(), params.end()); - return {evaluateWeakForm(solid_dynamic_weak_forms_->quasi_static_weak_form, time_info, shape_disp, solid_inputs, + return {evaluateWeakForm(solid_dynamic_weak_forms_->final_reaction_weak_form, time_info, shape_disp, solid_inputs, states[DISPLACEMENT])}; } diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index bde1c31fec..40e4774fd7 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -67,10 +67,12 @@ class TimeDiscretizedWeakForm time_discretized_weak_form; ///< this publically available abstract weak form is a + std::shared_ptr time_discretized_weak_form; ///< this publicaly available abstract weak form is a ///< functions of the current u, u_old, v_old, and a_old, - std::shared_ptr quasi_static_weak_form; ///< this publically available abstract weak form is structly a - ///< function of the current u, v, and a (no time discretization) + std::shared_ptr + final_reaction_weak_form; ///< this publicaly available abstract weak form is structly a + ///< function of the current u, v, and a (no time discretization) + ///< its main purpose is to compute reaction forces after the solve is completed }; template > @@ -93,9 +95,8 @@ class SecondOrderTimeDiscretizedWeakForm>; ///< using - using QuasiStaticWeakFormT = - TimeDiscretizedWeakForm>; ///< using + TimeDiscretizedWeakForm>; ///< using + using FinalReactionFormT = TimeDiscretizedWeakForm>; ///< using /// @brief Constructor SecondOrderTimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr mesh, @@ -110,9 +111,9 @@ class SecondOrderTimeDiscretizedWeakForm(physics_name, mesh, output_mfem_space, input_mfem_spaces_trial_removed); - quasi_static_weak_form = quasi_static_weak_form_; + final_reaction_weak_form_ = + std::make_shared(physics_name, mesh, output_mfem_space, input_mfem_spaces_trial_removed); + final_reaction_weak_form = final_reaction_weak_form_; } /// @overload @@ -129,8 +130,8 @@ class SecondOrderTimeDiscretizedWeakFormaddBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, body_name, - integrand); + final_reaction_weak_form_->addBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, + body_name, integrand); } /// @overload @@ -144,9 +145,9 @@ class SecondOrderTimeDiscretizedWeakForm time_discretized_weak_form_; ///< fully templated time discretized weak form (with time integration rule injected ///< into the q-function) - std::shared_ptr - quasi_static_weak_form_; ///< fully template underlying weak form (no time integration included, a function of - ///< current u, v, and a) + std::shared_ptr + final_reaction_weak_form_; ///< fully template underlying weak form (no time integration included, a function of + ///< current u, v, and a) SecondOrderTimeIntegrationRule time_rule_; ///< encodes the time integration rule }; From c2ced3f6a2e7850a35c087b97e8df76bbe3746eb Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 20 Jan 2026 14:08:37 -0800 Subject: [PATCH 07/17] Trying to simplify adding new physics even more. --- .../differentiable_numerics/field_state.hpp | 10 +- .../tests/CMakeLists.txt | 1 + .../tests/test_thermo_mechanics.cpp | 233 ++++++++++++++++++ .../time_discretized_weak_form.hpp | 3 +- 4 files changed, 244 insertions(+), 3 deletions(-) create mode 100644 src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp diff --git a/src/smith/differentiable_numerics/field_state.hpp b/src/smith/differentiable_numerics/field_state.hpp index e9528faf89..a2937cad27 100644 --- a/src/smith/differentiable_numerics/field_state.hpp +++ b/src/smith/differentiable_numerics/field_state.hpp @@ -212,22 +212,28 @@ inline std::vector spaces(const std::vector< }; /// @brief Get a vector of FieldPtr or DualFieldPtr from a vector of FieldState -inline std::vector getFieldPointers(std::vector& states) +inline std::vector getFieldPointers(std::vector& states, std::vector params = {}) { std::vector pointers; for (auto& t : states) { pointers.push_back(t.get().get()); } + for (auto& p : params) { + pointers.push_back(p.get().get()); + } return pointers; } /// @brief Get a vector of ConstFieldPtr or ConstDualFieldPtr from a vector of FieldState -inline std::vector getConstFieldPointers(const std::vector& states) +inline std::vector getConstFieldPointers(const std::vector& states, const std::vector& params = {}) { std::vector pointers; for (auto& t : states) { pointers.push_back(t.get().get()); } + for (auto& p : params) { + pointers.push_back(p.get().get()); + } return pointers; } diff --git a/src/smith/differentiable_numerics/tests/CMakeLists.txt b/src/smith/differentiable_numerics/tests/CMakeLists.txt index dc5925a793..06681ec856 100644 --- a/src/smith/differentiable_numerics/tests/CMakeLists.txt +++ b/src/smith/differentiable_numerics/tests/CMakeLists.txt @@ -10,6 +10,7 @@ set(differentiable_numerics_test_source test_field_state.cpp test_explicit_dynamics.cpp test_solid_mechanics_state_advancer.cpp + test_thermo_mechanics.cpp ) smith_add_tests( SOURCES ${differentiable_numerics_test_source} diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp new file mode 100644 index 0000000000..a344054348 --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -0,0 +1,233 @@ +// 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 "gtest/gtest.h" + +#include "smith/smith_config.hpp" +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/solver_config.hpp" +#include "smith/mesh_utils/mesh_utils.hpp" + +#include "smith/physics/state/state_manager.hpp" +#include "smith/physics/functional_weak_form.hpp" +#include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/differentiable_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/paraview_writer.hpp" +#include "smith/differentiable_numerics/differentiable_test_utils.hpp" + +namespace smith { + +smith::LinearSolverOptions solid_linear_options{.linear_solver = smith::LinearSolver::CG, + .preconditioner = smith::Preconditioner::HypreJacobi, + .relative_tol = 1e-11, + .absolute_tol = 1e-11, + .max_iterations = 10000, + .print_level = 0}; + +smith::NonlinearSolverOptions solid_nonlinear_opts{.nonlin_solver = NonlinearSolver::TrustRegion, + .relative_tol = 1.0e-10, + .absolute_tol = 1.0e-10, + .max_iterations = 500, + .print_level = 0}; + +static constexpr int dim = 3; +static constexpr int order = 1; + +struct SolidMechanicsMeshFixture : public testing::Test { + double length = 1.0; + double width = 0.04; + int num_elements_x = 12; + int num_elements_y = 2; + int num_elements_z = 2; + double elem_size = length / num_elements_x; + + void SetUp() + { + smith::StateManager::initialize(datastore, "solid"); + auto mfem_shape = mfem::Element::QUADRILATERAL; + mesh = std::make_shared( + mfem::Mesh::MakeCartesian3D(num_elements_x, num_elements_y, num_elements_z, mfem_shape, length, width, width), + "mesh", 0, 0); + mesh->addDomainOfBoundaryElements("left", smith::by_attr(3)); + mesh->addDomainOfBoundaryElements("right", smith::by_attr(5)); + } + + static constexpr double total_simulation_time_ = 1.1; + static constexpr size_t num_steps_ = 4; + static constexpr double dt_ = total_simulation_time_ / num_steps_; + + axom::sidre::DataStore datastore; + std::shared_ptr mesh; +}; + +template +struct FieldType { + FieldType(std::string n) : name(n) {} + std::string name; +}; + +enum class SecondOrderTimeDiscretization +{ + QuasiStatic, + ImplicitNewmark +}; + +enum class FirstOrderTimeDiscretization +{ + QuasiStatic, + BackwardEuler +}; + +struct FieldStore { + FieldStore(std::shared_ptr mesh, size_t storage_size = 50) + : mesh_(mesh), data_store_(std::make_shared(storage_size)) + { + } + + template + void addShapeDisp(FieldType type) + { + shape_disp_.push_back(smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag())); + } + + template + void addUnknown(FieldType type) + { + to_unknowns_[type.name] = unknowns_.size(); + unknowns_.push_back(smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag())); + } + + template + void addDerived(FieldType type, SecondOrderTimeDiscretization time_discretization, std::vector custom_names = {}) + { + if (!custom_names.size()) { + custom_names.push_back(type.name + "_dot"); + custom_names.push_back(type.name + "_ddot"); + } + SLIC_ERROR_IF(custom_names.size()!=2, "Second order time discretized fields must add two derived fields f_dot and f_ddot"); + to_derived_[custom_names[0]] = derived_.size(); + derived_.push_back(smith::createFieldState(*data_store_, Space{}, custom_names[0], mesh_->tag())); + + to_derived_[custom_names[1]] = derived_.size(); + derived_.push_back(smith::createFieldState(*data_store_, Space{}, custom_names[1], mesh_->tag())); + } + + template + void addDerived(FieldType type, FirstOrderTimeDiscretization time_discretization, std::vector custom_names = {}) + { + if (!custom_names.size()) { + custom_names.push_back(type.name + "_dot"); + } + SLIC_ERROR_IF(custom_names.size()!=1, "First order time discretized fields must add exactly one derived field: f_dot"); + to_derived_[custom_names[0]] = derived_.size(); + derived_.push_back(smith::createFieldState(*data_store_, Space{}, custom_names[0], mesh_->tag())); + } + + std::shared_ptr mesh_; + std::shared_ptr data_store_; + + std::vector shape_disp_; + + std::vector unknowns_; + std::vector derived_; + std::vector params_; + + std::map to_unknowns_; + std::map to_derived_; + std::map to_params_; +}; + +struct NewmarkDot +{ + template + SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return (2.0 / t.dt()) * (field_new - field_old) - velo_old; + } +}; + + +template +void createSpaces(const FieldStore& field_store, std::vector& spaces, FirstType type) +{ + const size_t test_index = field_store.to_unknowns_.at(type.name); + spaces.push_back(&field_store.unknowns_[test_index].get()->space()); +} + + +template +void createSpaces(const FieldStore& field_store, std::vector& spaces, FirstType type, Types... types) +{ + const size_t test_index = field_store.to_unknowns_.at(type.name); + spaces.insert(spaces.begin(), &field_store.unknowns_[test_index].get()->space()); + createSpaces(field_store, spaces, types...); +} + + +template +auto createWeakForm(const FieldStore& field_store, std::string name, FieldType test_type, FieldType... field_types) +{ + const size_t test_index = field_store.to_unknowns_.at(test_type.name); + const mfem::ParFiniteElementSpace& test_space = field_store.unknowns_[test_index].get()->space(); + std::vector input_spaces; + createSpaces(field_store, input_spaces, field_types...); + auto weak_form = std::make_shared > >(name, field_store.mesh_, test_space, input_spaces); + return weak_form; +} + + +TEST_F(SolidMechanicsMeshFixture, A) +{ + SMITH_MARK_FUNCTION; + + FieldType> shape_disp_type("shape_displacement"); + FieldType> disp_type("displacement"); + FieldType> temperature_type("temperature"); + + FieldStore field_store(mesh, 100); + + field_store.addShapeDisp(shape_disp_type); + + field_store.addUnknown(disp_type); + field_store.addUnknown(temperature_type); + + field_store.addDerived(disp_type, SecondOrderTimeDiscretization::ImplicitNewmark, {"velocity", "acceleration"}); + field_store.addDerived(temperature_type, FirstOrderTimeDiscretization::BackwardEuler); + + auto weak_form = createWeakForm(field_store, "solid", disp_type, disp_type, temperature_type); + + weak_form->addBodyIntegral(smith::DependsOn<0,1>{}, mesh->entireBodyName(), [](auto /*t*/, auto /*X*/, auto u, auto temperature) { + auto ones = 0.0*get(u); + ones[0] = 1.0; + return smith::tuple{get(u) + ones, get(u)}; + }); + + auto r = weak_form->residual(TimeInfo(0.0,0.0,0), field_store.shape_disp_[0].get().get(), getConstFieldPointers(field_store.unknowns_, field_store.derived_)); + // std::cout << "norm = " << r.Norml2() << std::endl; + + EXPECT_EQ(0,0); + + // auto derived_types = createDerivedFieldTypes(ddot(disp_type), dot(temperature_type)); + + // FieldType> kappa_type; + // FieldType> mu_type; + // auto params = + + // auto weak_form = + // createWeakForm("solid", smith::DependsOnFields(disp_type, dot(disp_type), ddot(disp_type), temperature_type)); +} + +} // namespace smith + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index 40e4774fd7..8afac770dc 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -56,9 +56,10 @@ class TimeDiscretizedWeakForm + template // int... all_active_parameters> void addBodyIntegral(std::string body_name, BodyForceType body_integral) { + auto input_indices = std::make_integer_sequence{}; addBodyIntegral(DependsOn<>{}, body_name, body_integral); } }; From 8b30cd8ca21a6fee62ca33f834a3c7dd7bed3ace Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Wed, 21 Jan 2026 15:14:55 -0800 Subject: [PATCH 08/17] Working through new interface ideas, change newton linesearch default. --- .../solid_mechanics_state_advancer.cpp | 4 +- .../tests/test_thermo_mechanics.cpp | 92 +++++++++---------- .../time_discretized_weak_form.hpp | 16 +++- .../time_integration_rule.hpp | 6 +- src/smith/numerics/solver_config.hpp | 2 +- 5 files changed, 60 insertions(+), 60 deletions(-) diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp index d60e9dfd1e..330ff9661e 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp @@ -46,9 +46,9 @@ std::vector SolidMechanicsStateAdvancer::advanceState(const TimeInfo std::vector states = states_old; states[DISPLACEMENT] = displacement; - states[VELOCITY] = time_rule_.derivative(final_time_info, displacement, states_old[DISPLACEMENT], + states[VELOCITY] = time_rule_.dot(final_time_info, displacement, states_old[DISPLACEMENT], states_old[VELOCITY], states_old[ACCELERATION]); - states[ACCELERATION] = time_rule_.second_derivative(final_time_info, displacement, states_old[DISPLACEMENT], + states[ACCELERATION] = time_rule_.ddot(final_time_info, displacement, states_old[DISPLACEMENT], states_old[VELOCITY], states_old[ACCELERATION]); return states; diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index a344054348..c7a716b7e7 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -16,6 +16,8 @@ #include "smith/differentiable_numerics/field_state.hpp" #include "smith/differentiable_numerics/differentiable_solver.hpp" #include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" #include "smith/differentiable_numerics/paraview_writer.hpp" #include "smith/differentiable_numerics/differentiable_test_utils.hpp" @@ -64,7 +66,7 @@ struct SolidMechanicsMeshFixture : public testing::Test { std::shared_ptr mesh; }; -template +template struct FieldType { FieldType(std::string n) : name(n) {} std::string name; @@ -97,34 +99,16 @@ struct FieldStore { template void addUnknown(FieldType type) { - to_unknowns_[type.name] = unknowns_.size(); - unknowns_.push_back(smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag())); + to_fields_[type.name] = fields_.size(); + fields_.push_back(smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag())); } template - void addDerived(FieldType type, SecondOrderTimeDiscretization time_discretization, std::vector custom_names = {}) + auto addDerived(FieldType type, std::string name) { - if (!custom_names.size()) { - custom_names.push_back(type.name + "_dot"); - custom_names.push_back(type.name + "_ddot"); - } - SLIC_ERROR_IF(custom_names.size()!=2, "Second order time discretized fields must add two derived fields f_dot and f_ddot"); - to_derived_[custom_names[0]] = derived_.size(); - derived_.push_back(smith::createFieldState(*data_store_, Space{}, custom_names[0], mesh_->tag())); - - to_derived_[custom_names[1]] = derived_.size(); - derived_.push_back(smith::createFieldState(*data_store_, Space{}, custom_names[1], mesh_->tag())); - } - - template - void addDerived(FieldType type, FirstOrderTimeDiscretization time_discretization, std::vector custom_names = {}) - { - if (!custom_names.size()) { - custom_names.push_back(type.name + "_dot"); - } - SLIC_ERROR_IF(custom_names.size()!=1, "First order time discretized fields must add exactly one derived field: f_dot"); - to_derived_[custom_names[0]] = derived_.size(); - derived_.push_back(smith::createFieldState(*data_store_, Space{}, custom_names[0], mesh_->tag())); + to_fields_[name] = fields_.size(); + fields_.push_back(smith::createFieldState(*data_store_, Space{}, name, mesh_->tag())); + return FieldType(name); } std::shared_ptr mesh_; @@ -132,13 +116,8 @@ struct FieldStore { std::vector shape_disp_; - std::vector unknowns_; - std::vector derived_; - std::vector params_; - - std::map to_unknowns_; - std::map to_derived_; - std::map to_params_; + std::vector fields_; + std::map to_fields_; }; struct NewmarkDot @@ -154,30 +133,31 @@ struct NewmarkDot template -void createSpaces(const FieldStore& field_store, std::vector& spaces, FirstType type) +void createSpaces(const FieldStore& field_store, std::vector& spaces, FieldType type) { - const size_t test_index = field_store.to_unknowns_.at(type.name); - spaces.push_back(&field_store.unknowns_[test_index].get()->space()); + const size_t test_index = field_store.to_fields_.at(type.name); + spaces.push_back(&field_store.fields_[test_index].get()->space()); } template void createSpaces(const FieldStore& field_store, std::vector& spaces, FirstType type, Types... types) { - const size_t test_index = field_store.to_unknowns_.at(type.name); - spaces.insert(spaces.begin(), &field_store.unknowns_[test_index].get()->space()); + const size_t test_index = field_store.to_fields_.at(type.name); + spaces.insert(spaces.begin(), &field_store.fields_[test_index].get()->space()); createSpaces(field_store, spaces, types...); } template -auto createWeakForm(const FieldStore& field_store, std::string name, FieldType test_type, FieldType... field_types) +auto createWeakForm(std::string name, FieldStore& field_store, FieldType test_type, FieldType... field_types) { - const size_t test_index = field_store.to_unknowns_.at(test_type.name); - const mfem::ParFiniteElementSpace& test_space = field_store.unknowns_[test_index].get()->space(); + const size_t test_index = field_store.to_fields_.at(test_type.name); + const mfem::ParFiniteElementSpace& test_space = field_store.fields_[test_index].get()->space(); std::vector input_spaces; createSpaces(field_store, input_spaces, field_types...); - auto weak_form = std::make_shared > >(name, field_store.mesh_, test_space, input_spaces); + // std::cout << "num spaces = " << input_spaces.size() << std::endl; + auto weak_form = std::make_shared > >(name, field_store.mesh_, test_space, input_spaces); return weak_form; } @@ -197,18 +177,32 @@ TEST_F(SolidMechanicsMeshFixture, A) field_store.addUnknown(disp_type); field_store.addUnknown(temperature_type); - field_store.addDerived(disp_type, SecondOrderTimeDiscretization::ImplicitNewmark, {"velocity", "acceleration"}); - field_store.addDerived(temperature_type, FirstOrderTimeDiscretization::BackwardEuler); + auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); + auto velo_old_type = field_store.addDerived(disp_type, "velocity_old"); + auto accel_old_type = field_store.addDerived(disp_type, "acceleration_old"); + auto temperature_old_type = field_store.addDerived(temperature_type, "temperature_old"); + + auto disp_time_rule = SecondOrderTimeIntegrationRule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + auto temperature_time_rule = BackwardEulerFirstOrderTimeIntegrationRule(); - auto weak_form = createWeakForm(field_store, "solid", disp_type, disp_type, temperature_type); + auto solid_weak_form = createWeakForm("solid", field_store, disp_type, disp_type, temperature_type, disp_old_type, velo_old_type, accel_old_type, temperature_old_type); - weak_form->addBodyIntegral(smith::DependsOn<0,1>{}, mesh->entireBodyName(), [](auto /*t*/, auto /*X*/, auto u, auto temperature) { - auto ones = 0.0*get(u); - ones[0] = 1.0; - return smith::tuple{get(u) + ones, get(u)}; + solid_weak_form->addBodyIntegral(mesh->entireBodyName(), [=](auto time_info, auto /*X*/, auto disp, auto temperature, auto disp_old, auto velo_old, auto accel_old, auto temperature_old) { + auto u = disp_time_rule.value(time_info, disp, disp_old, velo_old, accel_old); + auto a = disp_time_rule.ddot(time_info, disp, disp_old, velo_old, accel_old); + auto theta_dot = temperature_time_rule.dot(time_info, temperature, temperature_old); + return smith::tuple{2.0 * get(a), get(u) * get(theta_dot)}; }); - auto r = weak_form->residual(TimeInfo(0.0,0.0,0), field_store.shape_disp_[0].get().get(), getConstFieldPointers(field_store.unknowns_, field_store.derived_)); + // auto thermal_weak_form = createWeakForm("thermal", field_store, UnknownArgs{}, DerivedArgs<>{}, ParameterArgs<>{}); + // thermal_weak_form->addBodyIntegral(mesh->entireBodyName(), [](auto /*t*/, auto /*X*/, auto temperature, auto u) { + // auto ones = 0.0*get(u); + // ones[0] = 1.0; + // return smith::tuple{get(u) + ones, get(u)}; + // }); + + // solve({wf_1, wf_2}, field_store, solver, bcs, ...); + //auto r = weak_form->residual(TimeInfo(0.0,0.0,0), field_store.shape_disp_[0].get().get(), getConstFieldPointers(field_store.unknowns_, field_store.derived_)); // std::cout << "norm = " << r.Norml2() << std::endl; EXPECT_EQ(0,0); diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index 8afac770dc..0d1a44a272 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -56,11 +56,17 @@ class TimeDiscretizedWeakForm // int... all_active_parameters> + template + void addBodyIntegralImpl(std::string body_name, BodyForceType body_integral, std::integer_sequence) + { + addBodyIntegral(DependsOn{}, body_name, body_integral); + } + + /// @overload + template void addBodyIntegral(std::string body_name, BodyForceType body_integral) { - auto input_indices = std::make_integer_sequence{}; - addBodyIntegral(DependsOn<>{}, body_name, body_integral); + addBodyIntegralImpl(body_name, body_integral, std::make_integer_sequence{}); } }; @@ -128,8 +134,8 @@ class SecondOrderTimeDiscretizedWeakFormaddBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, body_name, integrand); diff --git a/src/smith/differentiable_numerics/time_integration_rule.hpp b/src/smith/differentiable_numerics/time_integration_rule.hpp index dc4e0cd41b..9902d4eef9 100644 --- a/src/smith/differentiable_numerics/time_integration_rule.hpp +++ b/src/smith/differentiable_numerics/time_integration_rule.hpp @@ -34,7 +34,7 @@ struct BackwardEulerFirstOrderTimeIntegrationRule { /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template - SMITH_HOST_DEVICE auto derivative(const TimeInfo& t, const T1& field_new, const T2& field_old) const + SMITH_HOST_DEVICE auto dot(const TimeInfo& t, const T1& field_new, const T2& field_old) const { return (1.0 / t.dt()) * (field_new - field_old); } @@ -67,7 +67,7 @@ struct SecondOrderTimeIntegrationRule { /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template - SMITH_HOST_DEVICE auto derivative([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, [[maybe_unused]] const T4& accel_old) const { @@ -76,7 +76,7 @@ struct SecondOrderTimeIntegrationRule { /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template - SMITH_HOST_DEVICE auto second_derivative([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + SMITH_HOST_DEVICE auto ddot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, [[maybe_unused]] const T4& accel_old) const { diff --git a/src/smith/numerics/solver_config.hpp b/src/smith/numerics/solver_config.hpp index 3ab439046c..56890a4b37 100644 --- a/src/smith/numerics/solver_config.hpp +++ b/src/smith/numerics/solver_config.hpp @@ -431,7 +431,7 @@ struct NonlinearSolverOptions { int max_iterations = 20; /// Maximum line search cutbacks - int max_line_search_iterations = 0; + int max_line_search_iterations = 8; /// Debug print level int print_level = 0; From 39819b30f6c90fa7e0ce42c6ed3450acd1554f02 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Fri, 23 Jan 2026 12:55:16 -0800 Subject: [PATCH 09/17] Some minor adjustments to interface, change default for NewtonLinearsearch to have no-zero line searches. --- .../differentiable_numerics/field_state.hpp | 6 +- .../nonlinear_solve.cpp | 2 +- .../nonlinear_solve.hpp | 2 +- .../solid_mechanics_state_advancer.cpp | 12 +- .../tests/test_thermo_mechanics.cpp | 354 +++++++++++++----- .../time_discretized_weak_form.hpp | 3 +- .../time_integration_rule.hpp | 8 +- src/smith/numerics/equation_solver.cpp | 6 +- src/smith/numerics/equation_solver.hpp | 2 +- src/smith/physics/common.hpp | 11 +- 10 files changed, 289 insertions(+), 117 deletions(-) diff --git a/src/smith/differentiable_numerics/field_state.hpp b/src/smith/differentiable_numerics/field_state.hpp index a2937cad27..4f2322c304 100644 --- a/src/smith/differentiable_numerics/field_state.hpp +++ b/src/smith/differentiable_numerics/field_state.hpp @@ -212,7 +212,8 @@ inline std::vector spaces(const std::vector< }; /// @brief Get a vector of FieldPtr or DualFieldPtr from a vector of FieldState -inline std::vector getFieldPointers(std::vector& states, std::vector params = {}) +inline std::vector getFieldPointers(std::vector& states, + std::vector params = {}) { std::vector pointers; for (auto& t : states) { @@ -225,7 +226,8 @@ inline std::vector getFieldPointers(std::vector } /// @brief Get a vector of ConstFieldPtr or ConstDualFieldPtr from a vector of FieldState -inline std::vector getConstFieldPointers(const std::vector& states, const std::vector& params = {}) +inline std::vector getConstFieldPointers(const std::vector& states, + const std::vector& params = {}) { std::vector pointers; for (auto& t : states) { diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 1f1f34a558..37777a1734 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -272,7 +272,7 @@ std::vector block_solve(const std::vector& residual_evals const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const DifferentiableBlockSolver* solver, - const std::vector bc_managers) + const std::vector& bc_managers) { SMITH_MARK_FUNCTION; size_t num_rows_ = residual_evals.size(); diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp index b89ef78579..eebfbdd40a 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.hpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -64,6 +64,6 @@ std::vector block_solve(const std::vector& residual_evals const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const DifferentiableBlockSolver* solver, - const std::vector bc_managers); + const std::vector& bc_managers); } // namespace smith diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp index 330ff9661e..b40d58e84c 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp @@ -35,21 +35,19 @@ std::vector SolidMechanicsStateAdvancer::advanceState(const TimeInfo params, time_info, *solver_, *vector_bcs_, ACCELERATION); } - TimeInfo final_time_info = time_info.endTimeInfo(); - std::vector solid_inputs{states_old[DISPLACEMENT], states_old[DISPLACEMENT], states_old[VELOCITY], states_old[ACCELERATION]}; auto displacement = solve(*solid_dynamic_weak_forms_->time_discretized_weak_form, shape_disp, solid_inputs, params, - final_time_info, *solver_, *vector_bcs_); + time_info, *solver_, *vector_bcs_); std::vector states = states_old; states[DISPLACEMENT] = displacement; - states[VELOCITY] = time_rule_.dot(final_time_info, displacement, states_old[DISPLACEMENT], - states_old[VELOCITY], states_old[ACCELERATION]); - states[ACCELERATION] = time_rule_.ddot(final_time_info, displacement, states_old[DISPLACEMENT], - states_old[VELOCITY], states_old[ACCELERATION]); + states[VELOCITY] = + time_rule_.dot(time_info, displacement, states_old[DISPLACEMENT], states_old[VELOCITY], states_old[ACCELERATION]); + states[ACCELERATION] = time_rule_.ddot(time_info, displacement, states_old[DISPLACEMENT], states_old[VELOCITY], + states_old[ACCELERATION]); return states; } diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index c7a716b7e7..63579559b7 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -9,10 +9,9 @@ #include "smith/smith_config.hpp" #include "smith/infrastructure/application_manager.hpp" #include "smith/numerics/solver_config.hpp" -#include "smith/mesh_utils/mesh_utils.hpp" #include "smith/physics/state/state_manager.hpp" -#include "smith/physics/functional_weak_form.hpp" + #include "smith/differentiable_numerics/field_state.hpp" #include "smith/differentiable_numerics/differentiable_solver.hpp" #include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" @@ -20,21 +19,79 @@ #include "smith/differentiable_numerics/time_discretized_weak_form.hpp" #include "smith/differentiable_numerics/paraview_writer.hpp" #include "smith/differentiable_numerics/differentiable_test_utils.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" namespace smith { -smith::LinearSolverOptions solid_linear_options{.linear_solver = smith::LinearSolver::CG, - .preconditioner = smith::Preconditioner::HypreJacobi, - .relative_tol = 1e-11, - .absolute_tol = 1e-11, - .max_iterations = 10000, - .print_level = 0}; +/// @brief Green-Saint Venant isotropic thermoelastic model +struct GreenSaintVenantThermoelasticMaterial { + double density; ///< density + double E; ///< Young's modulus + double nu; ///< Poisson's ratio + double C_v; ///< volumetric heat capacity + double alpha; ///< thermal expansion coefficient + double theta_ref; ///< datum temperature for thermal expansion + double kappa; ///< thermal conductivity + + using State = Empty; + + /** + * @brief Evaluate constitutive variables for thermomechanics + * + * @tparam T1 Type of the displacement gradient components (number-like) + * @tparam T2 Type of the temperature (number-like) + * @tparam T3 Type of the temperature gradient components (number-like) + * + * @param[in] grad_u Displacement gradient + * @param[in] theta Temperature + * @param[in] grad_theta Temperature gradient + * @param[in,out] state State variables for this material + * + * @return[out] tuple of constitutive outputs. Contains the + * First Piola stress, the volumetric heat capacity in the reference + * configuration, the heat generated per unit volume during the time + * step (units of energy), and the referential heat flux (units of + * energy per unit time and per unit area). + */ + template + auto operator()(double, State&, const tensor& grad_u, const tensor& grad_v, T3 theta, + const tensor& grad_theta) const + { + const double K = E / (3.0 * (1.0 - 2.0 * nu)); + const double G = 0.5 * E / (1.0 + nu); + const auto Eg = greenStrain(grad_u); + const auto trEg = tr(Eg); + + // stress + static constexpr auto I = Identity(); + const auto S = 2.0 * G * dev(Eg) + K * (trEg - dim * alpha * (theta - theta_ref)) * I; + auto F = grad_u + I; + const auto Piola = dot(F, S); + + // internal heat power + auto greenStrainRate = + 0.5 * (grad_v + transpose(grad_v) + dot(transpose(grad_v), grad_u) + dot(transpose(grad_u), grad_v)); + const auto s0 = -dim * K * alpha * (theta + 273.1) * tr(greenStrainRate); + + // heat flux + const auto q0 = -kappa * grad_theta; + + return smith::tuple{Piola, C_v, s0, q0}; + } +}; + +smith::LinearSolverOptions linear_options{.linear_solver = smith::LinearSolver::Strumpack, + .relative_tol = 1e-8, + .absolute_tol = 1e-8, + .max_iterations = 200, + .print_level = 0}; -smith::NonlinearSolverOptions solid_nonlinear_opts{.nonlin_solver = NonlinearSolver::TrustRegion, - .relative_tol = 1.0e-10, - .absolute_tol = 1.0e-10, - .max_iterations = 500, - .print_level = 0}; +smith::NonlinearSolverOptions nonlinear_opts{.nonlin_solver = NonlinearSolver::NewtonLineSearch, + .relative_tol = 1.9e-6, + .absolute_tol = 1.0e-10, + .max_iterations = 500, + .max_line_search_iterations = 50, + .print_level = 2}; static constexpr int dim = 3; static constexpr int order = 1; @@ -49,39 +106,28 @@ struct SolidMechanicsMeshFixture : public testing::Test { void SetUp() { - smith::StateManager::initialize(datastore, "solid"); + smith::StateManager::initialize(datastore_, "solid"); auto mfem_shape = mfem::Element::QUADRILATERAL; - mesh = std::make_shared( + mesh_ = std::make_shared( mfem::Mesh::MakeCartesian3D(num_elements_x, num_elements_y, num_elements_z, mfem_shape, length, width, width), "mesh", 0, 0); - mesh->addDomainOfBoundaryElements("left", smith::by_attr(3)); - mesh->addDomainOfBoundaryElements("right", smith::by_attr(5)); + mesh_->addDomainOfBoundaryElements("left", smith::by_attr(3)); + mesh_->addDomainOfBoundaryElements("right", smith::by_attr(5)); } static constexpr double total_simulation_time_ = 1.1; static constexpr size_t num_steps_ = 4; static constexpr double dt_ = total_simulation_time_ / num_steps_; - axom::sidre::DataStore datastore; - std::shared_ptr mesh; + axom::sidre::DataStore datastore_; + std::shared_ptr mesh_; }; -template +template struct FieldType { - FieldType(std::string n) : name(n) {} + FieldType(std::string n, int unknown_index_ = -1) : name(n), unknown_index(unknown_index_) {} std::string name; -}; - -enum class SecondOrderTimeDiscretization -{ - QuasiStatic, - ImplicitNewmark -}; - -enum class FirstOrderTimeDiscretization -{ - QuasiStatic, - BackwardEuler + int unknown_index; }; struct FieldStore { @@ -97,124 +143,248 @@ struct FieldStore { } template - void addUnknown(FieldType type) + std::shared_ptr& addUnknown(FieldType& type) { - to_fields_[type.name] = fields_.size(); - fields_.push_back(smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag())); + type.unknown_index = static_cast(num_unknowns_); + to_fields_index_[type.name] = fields_.size(); + to_unknown_index_[type.name] = num_unknowns_; + FieldState new_field = smith::createFieldState(*data_store_, Space{}, type.name, mesh_->tag()); + fields_.push_back(new_field); + ++num_unknowns_; + boundary_conditions_.push_back( + std::make_shared(mesh_->mfemParMesh(), new_field.get()->space())); + SLIC_ASSERT(num_unknowns_ == boundary_conditions.size()); + return boundary_conditions_.back(); } template - auto addDerived(FieldType type, std::string name) + auto addDerived(FieldType, std::string name) { - to_fields_[name] = fields_.size(); + to_fields_index_[name] = fields_.size(); fields_.push_back(smith::createFieldState(*data_store_, Space{}, name, mesh_->tag())); return FieldType(name); } + void 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{argument_name_and_index}; + } + } + + void addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index) + { + size_t field_index = to_fields_index_.at(argument_name); + if (weak_form_name_to_field_indices.count(weak_form_name)) { + weak_form_name_to_field_indices.at(weak_form_name).push_back(field_index); + } else { + weak_form_name_to_field_indices[weak_form_name] = std::vector{field_index}; + } + SLIC_ERROR_IF(argument_index + 1 != weak_form_name_to_field_indices.at(weak_form_name).size(), + "Invalid order for adding weak form arguments."); + } + + void printMap() + { + for (auto& keyval : weak_form_name_to_unknown_name_index_) { + 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 << ", "; + } + std::cout << std::endl; + } + } + + std::vector> indexMap(const std::vector& residual_names) const + { + std::vector> block_indices(residual_names.size()); + + for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { + std::vector& res_indices = block_indices[res_i]; + res_indices = std::vector(num_unknowns_, 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); + + for (const auto& field_name_and_arg_index : arg_info) { + const std::string field_name = field_name_and_arg_index.field_name; + size_t unknown_index = to_unknown_index_.at(field_name); + SLIC_ASSERT(unknown_index < num_unknowns_); + res_indices[unknown_index] = field_name_and_arg_index.field_index_in_residual; + } + } + + return block_indices; + } + + std::vector getBoundaryConditionManagers() const + { + std::vector bcs; + for (auto& bc : boundary_conditions_) { + bcs.push_back(&bc->getBoundaryConditionManager()); + } + return bcs; + } + + size_t getFieldIndex(const std::string& field_name) const { return to_fields_index_.at(field_name); } + + const FieldState& getField(size_t field_index) const { return fields_[field_index]; } + + const FieldState& getShapeDisp() const { return shape_disp_[0]; } + + const std::vector& getFields() const { return fields_; } + + const std::shared_ptr& getMesh() const { return mesh_; } + + private: std::shared_ptr mesh_; std::shared_ptr data_store_; std::vector shape_disp_; - std::vector fields_; - std::map to_fields_; -}; + std::map to_fields_index_; -struct NewmarkDot -{ - template - SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, - [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, - [[maybe_unused]] const T4& accel_old) const - { - return (2.0 / t.dt()) * (field_new - field_old) - velo_old; - } -}; + size_t num_unknowns_ = 0; + std::map to_unknown_index_; + std::vector> boundary_conditions_; + // MRT, do this for readability + struct FieldLabel { + std::string field_name; + size_t field_index_in_residual; + }; -template -void createSpaces(const FieldStore& field_store, std::vector& spaces, FieldType type) -{ - const size_t test_index = field_store.to_fields_.at(type.name); - spaces.push_back(&field_store.fields_[test_index].get()->space()); -} + std::map> weak_form_name_to_unknown_name_index_; + std::map> weak_form_name_to_field_indices; +}; template -void createSpaces(const FieldStore& field_store, std::vector& spaces, FirstType type, Types... types) +void createSpaces(const std::string& weak_form_name, FieldStore& field_store, + std::vector& spaces, size_t arg_num, FirstType type, + Types... types) { - const size_t test_index = field_store.to_fields_.at(type.name); - spaces.insert(spaces.begin(), &field_store.fields_[test_index].get()->space()); - createSpaces(field_store, spaces, types...); + const size_t test_index = field_store.getFieldIndex(type.name); + SLIC_ERROR_IF(spaces.size() != arg_num, "Error creating spaces recursively"); + spaces.push_back(&field_store.getField(test_index).get()->space()); + field_store.addWeakFormArg(weak_form_name, type.name, arg_num); + if (type.unknown_index >= 0) { + field_store.addWeakFormUnknownArg(weak_form_name, type.name, arg_num); + } + if constexpr (sizeof...(types) > 0) { + createSpaces(weak_form_name, field_store, spaces, arg_num + 1, types...); + } } - template -auto createWeakForm(std::string name, FieldStore& field_store, FieldType test_type, FieldType... field_types) +auto createWeakForm(std::string name, FieldType test_type, FieldStore& field_store, + FieldType... field_types) { - const size_t test_index = field_store.to_fields_.at(test_type.name); - const mfem::ParFiniteElementSpace& test_space = field_store.fields_[test_index].get()->space(); + const size_t test_index = field_store.getFieldIndex(test_type.name); + const mfem::ParFiniteElementSpace& test_space = field_store.getField(test_index).get()->space(); std::vector input_spaces; - createSpaces(field_store, input_spaces, field_types...); - // std::cout << "num spaces = " << input_spaces.size() << std::endl; - auto weak_form = std::make_shared > >(name, field_store.mesh_, test_space, input_spaces); - return weak_form; + createSpaces(name, field_store, input_spaces, 0, field_types...); + return std::make_shared>>( + name, field_store.getMesh(), test_space, input_spaces); } +std::vector solve(const std::vector& weak_forms, const FieldStore& field_store, + const DifferentiableBlockSolver* solver, const TimeInfo& time_info) +{ + std::vector weak_form_names; + for (const auto& wf : weak_forms) { + weak_form_names.push_back(wf->name()); + } + std::vector> index_map = field_store.indexMap(weak_form_names); + + for (auto is : index_map) { + for (auto i : is) { + std::cout << i << " "; + } + std::cout << std::endl; + } + + std::vector> params(weak_forms.size()); + std::vector> inputs(weak_forms.size()); + for (auto& s : inputs) { + s = field_store.getFields(); + } + + return block_solve(weak_forms, index_map, field_store.getShapeDisp(), inputs, params, time_info, solver, + field_store.getBoundaryConditionManagers()); +} TEST_F(SolidMechanicsMeshFixture, A) { SMITH_MARK_FUNCTION; FieldType> shape_disp_type("shape_displacement"); + FieldType> disp_type("displacement"); FieldType> temperature_type("temperature"); - FieldStore field_store(mesh, 100); + std::shared_ptr d_nonlinear_solver = + buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); + + FieldStore field_store(mesh_, 100); field_store.addShapeDisp(shape_disp_type); - field_store.addUnknown(disp_type); - field_store.addUnknown(temperature_type); + // addIndependentField, addDependantField + std::shared_ptr& disp_bc = field_store.addUnknown(disp_type); + disp_bc->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { + auto bc = 0.0 * X; + bc[0] = 0.01 * t; + return bc; + }); + vector_bcs->setFixedVectorBCs(mesh->domain("right")); auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); auto velo_old_type = field_store.addDerived(disp_type, "velocity_old"); auto accel_old_type = field_store.addDerived(disp_type, "acceleration_old"); + + std::shared_ptr& temperature_bc = field_store.addUnknown(temperature_type); + auto temperature_old_type = field_store.addDerived(temperature_type, "temperature_old"); auto disp_time_rule = SecondOrderTimeIntegrationRule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); auto temperature_time_rule = BackwardEulerFirstOrderTimeIntegrationRule(); - auto solid_weak_form = createWeakForm("solid", field_store, disp_type, disp_type, temperature_type, disp_old_type, velo_old_type, accel_old_type, temperature_old_type); + auto solid_weak_form = createWeakForm("solid_force", disp_type, field_store, disp_type, disp_old_type, + velo_old_type, accel_old_type, temperature_type, temperature_old_type); - solid_weak_form->addBodyIntegral(mesh->entireBodyName(), [=](auto time_info, auto /*X*/, auto disp, auto temperature, auto disp_old, auto velo_old, auto accel_old, auto temperature_old) { - auto u = disp_time_rule.value(time_info, disp, disp_old, velo_old, accel_old); - auto a = disp_time_rule.ddot(time_info, disp, disp_old, velo_old, accel_old); - auto theta_dot = temperature_time_rule.dot(time_info, temperature, temperature_old); - return smith::tuple{2.0 * get(a), get(u) * get(theta_dot)}; - }); + solid_weak_form->addBodyIntegral( + mesh_->entireBodyName(), [=](auto time_info, auto /*X*/, auto disp, auto disp_old, auto velo_old, auto accel_old, + auto temperature, auto temperature_old) { + auto u = disp_time_rule.value(time_info, disp, disp_old, velo_old, accel_old); + auto a = disp_time_rule.ddot(time_info, disp, disp_old, velo_old, accel_old); + auto theta_dot = temperature_time_rule.dot(time_info, temperature, temperature_old); + return smith::tuple{2.0 * get(a), get(u) * get(theta_dot)}; + }); - // auto thermal_weak_form = createWeakForm("thermal", field_store, UnknownArgs{}, DerivedArgs<>{}, ParameterArgs<>{}); - // thermal_weak_form->addBodyIntegral(mesh->entireBodyName(), [](auto /*t*/, auto /*X*/, auto temperature, auto u) { - // auto ones = 0.0*get(u); - // ones[0] = 1.0; - // return smith::tuple{get(u) + ones, get(u)}; - // }); + auto thermal_weak_form = createWeakForm("thermal_flux", temperature_type, field_store, temperature_type, + temperature_old_type, disp_type); - // solve({wf_1, wf_2}, field_store, solver, bcs, ...); - //auto r = weak_form->residual(TimeInfo(0.0,0.0,0), field_store.shape_disp_[0].get().get(), getConstFieldPointers(field_store.unknowns_, field_store.derived_)); - // std::cout << "norm = " << r.Norml2() << std::endl; + thermal_weak_form->addBodyIntegral( + mesh_->entireBodyName(), [=](auto time_info, auto /*X*/, auto temperature, auto temperature_old, auto) { + auto theta = temperature_time_rule.value(time_info, temperature, temperature_old); + auto theta_dot = temperature_time_rule.dot(time_info, temperature, temperature_old); + return smith::tuple{get(theta_dot), get(theta)}; + }); - EXPECT_EQ(0,0); + std::vector weak_forms{solid_weak_form.get(), thermal_weak_form.get()}; + std::vector disp_temp = solve(weak_forms, field_store, d_nonlinear_solver.get(), TimeInfo(0.0, 1.0)); - // auto derived_types = createDerivedFieldTypes(ddot(disp_type), dot(temperature_type)); + EXPECT_EQ(0, 0); + // auto derived_types = createDerivedFieldTypes(ddot(disp_type), dot(temperature_type)); // FieldType> kappa_type; // FieldType> mu_type; // auto params = - // auto weak_form = - // createWeakForm("solid", smith::DependsOnFields(disp_type, dot(disp_type), ddot(disp_type), temperature_type)); + // createWeakForm("solid", smith::DependsOnFields(disp_type, dot(disp_type), ddot(disp_type), temperature_type)); } } // namespace smith diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index 0d1a44a272..8cd5de4692 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -57,7 +57,8 @@ class TimeDiscretizedWeakForm - void addBodyIntegralImpl(std::string body_name, BodyForceType body_integral, std::integer_sequence) + void addBodyIntegralImpl(std::string body_name, BodyForceType body_integral, + std::integer_sequence) { addBodyIntegral(DependsOn{}, body_name, body_integral); } diff --git a/src/smith/differentiable_numerics/time_integration_rule.hpp b/src/smith/differentiable_numerics/time_integration_rule.hpp index 9902d4eef9..e88d51c488 100644 --- a/src/smith/differentiable_numerics/time_integration_rule.hpp +++ b/src/smith/differentiable_numerics/time_integration_rule.hpp @@ -68,8 +68,8 @@ struct SecondOrderTimeIntegrationRule { /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, - [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, - [[maybe_unused]] const T4& accel_old) const + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const { return (2.0 / t.dt()) * (field_new - field_old) - velo_old; } @@ -77,8 +77,8 @@ struct SecondOrderTimeIntegrationRule { /// @brief evaluate time derivative discretization of the ode state as used by the integration rule template SMITH_HOST_DEVICE auto ddot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, - [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, - [[maybe_unused]] const T4& accel_old) const + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const { auto dt = t.dt(); return (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old; diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index 6c24646dca..540ea62716 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -1019,15 +1019,15 @@ void StrumpackSolver::SetOperator(const mfem::Operator& op) #endif -std::unique_ptr buildNonlinearSolver(const NonlinearSolverOptions& nonlinear_opts, +std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions nonlinear_opts, const LinearSolverOptions& linear_opts, mfem::Solver& prec, MPI_Comm comm) { std::unique_ptr nonlinear_solver; if (nonlinear_opts.nonlin_solver == NonlinearSolver::Newton) { - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0 || nonlinear_opts.max_line_search_iterations != 0, - "Newton's method does not support nonzero min_iterations or max_line_search_iterations"); + nonlinear_opts.max_line_search_iterations = 0; + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, "Newton's method does not support nonzero min_iterations"); nonlinear_solver = std::make_unique(comm, nonlinear_opts); // nonlinear_solver = std::make_unique(comm); } else if (nonlinear_opts.nonlin_solver == NonlinearSolver::LBFGS) { diff --git a/src/smith/numerics/equation_solver.hpp b/src/smith/numerics/equation_solver.hpp index 83bc4673f5..33dcc42621 100644 --- a/src/smith/numerics/equation_solver.hpp +++ b/src/smith/numerics/equation_solver.hpp @@ -271,7 +271,7 @@ std::unique_ptr buildMonolithicMatrix(const mfem::BlockOpe * @param comm The MPI communicator for the supplied nonlinear operators and HypreParVectors * @return The constructed nonlinear solver */ -std::unique_ptr buildNonlinearSolver(const NonlinearSolverOptions& nonlinear_opts, +std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions nonlinear_opts, const LinearSolverOptions& linear_opts, mfem::Solver& preconditioner, MPI_Comm comm = MPI_COMM_WORLD); diff --git a/src/smith/physics/common.hpp b/src/smith/physics/common.hpp index d65f78b6d7..cde66b33d7 100644 --- a/src/smith/physics/common.hpp +++ b/src/smith/physics/common.hpp @@ -22,8 +22,11 @@ struct TimeInfo { { } - /// @brief accessor for time - double time() const { return time_.first; } + /// @brief accessor for time previous time + double time_old() const { return time_.first; } + + /// @brief accessor for the current time + double time() const { return time_.first + dt_.first; } /// @brief accessor for dt double dt() const { return dt_.first; } @@ -31,9 +34,7 @@ struct TimeInfo { /// @brief accessor for cycle size_t cycle() const { return cycle_; } - /// @brief return the time info corresponding to the end of this cycle (time + dt, dt, cycle) - TimeInfo endTimeInfo() const { return TimeInfo(time() + dt(), dt(), cycle()); } - + private: std::pair time_; ///< time and its dual std::pair dt_; ///< timestep and its dual size_t cycle_; ///< cycle, step, iteration count From c2684e82408bb5b3829ea364458926bf8b50e6f7 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 26 Jan 2026 15:39:23 -0800 Subject: [PATCH 10/17] Get an apparently working version of automating multiphysics solves. --- .../dirichlet_boundary_conditions.hpp | 28 ++++++------- .../tests/test_thermo_mechanics.cpp | 42 ++++++++++--------- 2 files changed, 37 insertions(+), 33 deletions(-) diff --git a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp index d7ddec9a55..6f65791e1b 100644 --- a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp +++ b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp @@ -94,35 +94,35 @@ class DirichletBoundaryConditions { template void setFixedScalarBCs(const Domain& domain) { - setVectorBCs(domain, [](auto, auto) { return 0.0; }); + setScalarBCs(domain, [](auto, auto) { return 0.0; }); } - /// @brief Constrain the dofs of a scalar field over a domain - template - void setFixedVectorBCs(const Domain& domain, int component) + /// @brief Constrain the vector dofs over a domain corresponding to a subset of the vector components + template + void setFixedVectorBCs(const Domain& domain, std::vector components) { - setScalarBCs(domain, component, [](auto, auto) { return 0.0; }); + setVectorBCs(domain, components, [](auto, auto) { return smith::tensor{}; }); } - /// @brief Constrain the vector dofs over a domain corresponding to a subset of the vector components - template - void setFixedVectorBCs(const Domain& domain, std::vector components) + /// @brief Constrain the dofs of a scalar field over a domain + template + void setFixedVectorBCs(const Domain& domain, int component) { - setVectorBCs(domain, components, [](auto, auto) { return smith::tensor{}; }); + std::vector components{component}; + setVectorBCs(domain, components); } /// @brief Constrain all the vector dofs over a domain - template + template void setFixedVectorBCs(const Domain& domain) { - const int field_dim = space_.GetVDim(); - SLIC_ERROR_IF(field_dim != spatial_dim, - "Vector boundary conditions current only work if they match the spatial dimension"); + SLIC_ERROR_IF(field_dim != space_.GetVDim(), + "Vector boundary condition field_dim does not match the fields vdim"); std::vector components(static_cast(field_dim)); for (int component = 0; component < field_dim; ++component) { components[static_cast(component)] = component; } - setFixedVectorBCs(domain, components); + setFixedVectorBCs(domain, components); } /// @brief Return the smith BoundaryConditionManager diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index 63579559b7..632760ddb7 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -178,12 +178,12 @@ struct FieldStore { void addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index) { size_t field_index = to_fields_index_.at(argument_name); - if (weak_form_name_to_field_indices.count(weak_form_name)) { - weak_form_name_to_field_indices.at(weak_form_name).push_back(field_index); + if (weak_form_name_to_field_indices_.count(weak_form_name)) { + weak_form_name_to_field_indices_.at(weak_form_name).push_back(field_index); } else { - weak_form_name_to_field_indices[weak_form_name] = std::vector{field_index}; + weak_form_name_to_field_indices_[weak_form_name] = std::vector{field_index}; } - SLIC_ERROR_IF(argument_index + 1 != weak_form_name_to_field_indices.at(weak_form_name).size(), + SLIC_ERROR_IF(argument_index + 1 != weak_form_name_to_field_indices_.at(weak_form_name).size(), "Invalid order for adding weak form arguments."); } @@ -234,7 +234,15 @@ struct FieldStore { const FieldState& getShapeDisp() const { return shape_disp_[0]; } - const std::vector& getFields() const { return fields_; } + std::vector getFields(const std::string& weak_form_name) const + { + auto unknown_field_indices = weak_form_name_to_field_indices_.at(weak_form_name); + std::vector fields_for_residual; + for (auto& i : unknown_field_indices) { + fields_for_residual.push_back(fields_[i]); + } + return fields_for_residual; + } const std::shared_ptr& getMesh() const { return mesh_; } @@ -258,7 +266,7 @@ struct FieldStore { std::map> weak_form_name_to_unknown_name_index_; - std::map> weak_form_name_to_field_indices; + std::map> weak_form_name_to_field_indices_; }; template @@ -299,18 +307,13 @@ std::vector solve(const std::vector& weak_forms, const Fi } std::vector> index_map = field_store.indexMap(weak_form_names); - for (auto is : index_map) { - for (auto i : is) { - std::cout << i << " "; - } - std::cout << std::endl; + std::vector> inputs; + for (size_t i = 0; i < weak_forms.size(); ++i) { + std::string wf_name = weak_forms[i]->name(); + std::vector fields_for_wk = field_store.getFields(wf_name); + inputs.push_back(fields_for_wk); } - std::vector> params(weak_forms.size()); - std::vector> inputs(weak_forms.size()); - for (auto& s : inputs) { - s = field_store.getFields(); - } return block_solve(weak_forms, index_map, field_store.getShapeDisp(), inputs, params, time_info, solver, field_store.getBoundaryConditionManagers()); @@ -332,20 +335,21 @@ TEST_F(SolidMechanicsMeshFixture, A) field_store.addShapeDisp(shape_disp_type); - // addIndependentField, addDependantField std::shared_ptr& disp_bc = field_store.addUnknown(disp_type); - disp_bc->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { + disp_bc->setVectorBCs(mesh_->domain("left"), [](double t, smith::tensor X) { auto bc = 0.0 * X; bc[0] = 0.01 * t; return bc; }); - vector_bcs->setFixedVectorBCs(mesh->domain("right")); + disp_bc->setFixedVectorBCs(mesh_->domain("right")); auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); auto velo_old_type = field_store.addDerived(disp_type, "velocity_old"); auto accel_old_type = field_store.addDerived(disp_type, "acceleration_old"); std::shared_ptr& temperature_bc = field_store.addUnknown(temperature_type); + temperature_bc->setFixedScalarBCs(mesh_->domain("left")); + temperature_bc->setFixedScalarBCs(mesh_->domain("right")); auto temperature_old_type = field_store.addDerived(temperature_type, "temperature_old"); From 8e093aef11e692169166432377be4990bf1542e9 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 26 Jan 2026 15:41:37 -0800 Subject: [PATCH 11/17] Fix style. --- .../dirichlet_boundary_conditions.hpp | 5 ++--- .../differentiable_numerics/tests/test_thermo_mechanics.cpp | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp index 6f65791e1b..a838191769 100644 --- a/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp +++ b/src/smith/differentiable_numerics/dirichlet_boundary_conditions.hpp @@ -113,11 +113,10 @@ class DirichletBoundaryConditions { } /// @brief Constrain all the vector dofs over a domain - template + template void setFixedVectorBCs(const Domain& domain) { - SLIC_ERROR_IF(field_dim != space_.GetVDim(), - "Vector boundary condition field_dim does not match the fields vdim"); + SLIC_ERROR_IF(field_dim != space_.GetVDim(), "Vector boundary condition field_dim does not match the fields vdim"); std::vector components(static_cast(field_dim)); for (int component = 0; component < field_dim; ++component) { components[static_cast(component)] = component; diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index 632760ddb7..9f91b08beb 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -341,7 +341,7 @@ TEST_F(SolidMechanicsMeshFixture, A) bc[0] = 0.01 * t; return bc; }); - disp_bc->setFixedVectorBCs(mesh_->domain("right")); + disp_bc->setFixedVectorBCs(mesh_->domain("right")); auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); auto velo_old_type = field_store.addDerived(disp_type, "velocity_old"); From 9631ff4faec455fbe6bdd5c32892aed032019158 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 27 Jan 2026 08:50:37 -0800 Subject: [PATCH 12/17] Small adjustment to equation solver defaults to ignore linesearch iterations when appropriate. --- src/smith/numerics/equation_solver.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index 540ea62716..35d9fee10e 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -1029,10 +1029,10 @@ std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions nonlinear_opts.max_line_search_iterations = 0; SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, "Newton's method does not support nonzero min_iterations"); nonlinear_solver = std::make_unique(comm, nonlinear_opts); - // nonlinear_solver = std::make_unique(comm); } else if (nonlinear_opts.nonlin_solver == NonlinearSolver::LBFGS) { - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0 || nonlinear_opts.max_line_search_iterations != 0, - "LBFGS does not support nonzero min_iterations or max_line_search_iterations"); + nonlinear_opts.max_line_search_iterations = 0; + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, + "LBFGS does not support nonzero min_iterations"); nonlinear_solver = std::make_unique(comm); } else if (nonlinear_opts.nonlin_solver == NonlinearSolver::NewtonLineSearch) { nonlinear_solver = std::make_unique(comm, nonlinear_opts); @@ -1052,9 +1052,9 @@ std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions // KINSOL else { #ifdef SMITH_USE_SUNDIALS - - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0 || nonlinear_opts.max_line_search_iterations != 0, - "kinsol solvers do not support min_iterations or max_line_search_iterations"); + nonlinear_opts.max_line_search_iterations = 0; + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, + "kinsol solvers do not support min_iterations"); int kinsol_strat = KIN_NONE; From 14fa86e0000a4307add0f32834037f9cd05caf5f Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 27 Jan 2026 08:53:14 -0800 Subject: [PATCH 13/17] Style. --- src/smith/numerics/equation_solver.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index 35d9fee10e..977c4ddaf0 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -1031,8 +1031,7 @@ std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions nonlinear_solver = std::make_unique(comm, nonlinear_opts); } else if (nonlinear_opts.nonlin_solver == NonlinearSolver::LBFGS) { nonlinear_opts.max_line_search_iterations = 0; - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, - "LBFGS does not support nonzero min_iterations"); + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, "LBFGS does not support nonzero min_iterations"); nonlinear_solver = std::make_unique(comm); } else if (nonlinear_opts.nonlin_solver == NonlinearSolver::NewtonLineSearch) { nonlinear_solver = std::make_unique(comm, nonlinear_opts); @@ -1053,8 +1052,7 @@ std::unique_ptr buildNonlinearSolver(NonlinearSolverOptions else { #ifdef SMITH_USE_SUNDIALS nonlinear_opts.max_line_search_iterations = 0; - SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, - "kinsol solvers do not support min_iterations"); + SLIC_ERROR_ROOT_IF(nonlinear_opts.min_iterations != 0, "kinsol solvers do not support min_iterations"); int kinsol_strat = KIN_NONE; From e4a2c733967d8f07fb82b817a63b8ee8977382d6 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 27 Jan 2026 09:54:40 -0800 Subject: [PATCH 14/17] Address an assert and remove or improve some comments. --- src/smith/differentiable_numerics/differentiable_solver.cpp | 4 +--- .../differentiable_numerics/tests/test_thermo_mechanics.cpp | 3 ++- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/smith/differentiable_numerics/differentiable_solver.cpp b/src/smith/differentiable_numerics/differentiable_solver.cpp index 7a3b256002..009ef1396a 100644 --- a/src/smith/differentiable_numerics/differentiable_solver.cpp +++ b/src/smith/differentiable_numerics/differentiable_solver.cpp @@ -273,7 +273,7 @@ NonlinearDifferentiableBlockSolver::NonlinearDifferentiableBlockSolver(std::uniq void NonlinearDifferentiableBlockSolver::completeSetup(const std::vector&) { - // initializeSolver(&nonlinear_solver_->preconditioner(), u); + // eventually may need something like: initializeSolver(&nonlinear_solver_->preconditioner(), u); } std::vector NonlinearDifferentiableBlockSolver::solve( @@ -310,8 +310,6 @@ std::vector NonlinearDifferentiableBlockSol *u_guesses[static_cast(row_i)] = u->GetBlock(row_i); } auto residuals = residual_funcs(u_guesses); - // auto block_r = std::make_unique(block_offsets); - // auto block_r = dynamic_cast(&r_); SLIC_ERROR_IF(!block_r, "Invalid r cast in block differentiable solver to a block vector"); for (int row_i = 0; row_i < num_rows; ++row_i) { auto r = residuals[static_cast(row_i)]; diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index 9f91b08beb..bdbe1ab47a 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -153,7 +153,8 @@ struct FieldStore { ++num_unknowns_; boundary_conditions_.push_back( std::make_shared(mesh_->mfemParMesh(), new_field.get()->space())); - SLIC_ASSERT(num_unknowns_ == boundary_conditions.size()); + SLIC_ERROR_IF(num_unknowns_ != boundary_conditions_.size(), + "Inconcistency between num unknowns and boundary condition size"); return boundary_conditions_.back(); } From 2dd9dfdc28aea524de2b2014992b37536ef8e166 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 27 Jan 2026 12:32:47 -0800 Subject: [PATCH 15/17] Working through thermal mechanics details. --- .../differentiable_solid_mechanics.hpp | 2 +- .../solid_mechanics_state_advancer.cpp | 2 +- .../solid_mechanics_state_advancer.hpp | 6 +- .../test_solid_mechanics_state_advancer.cpp | 6 +- .../tests/test_thermo_mechanics.cpp | 67 +++++++++++++------ .../time_discretized_weak_form.hpp | 6 +- .../time_integration_rule.hpp | 16 ++--- 7 files changed, 64 insertions(+), 41 deletions(-) diff --git a/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp b/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp index c6161392a5..e038af13e8 100644 --- a/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp +++ b/src/smith/differentiable_numerics/differentiable_solid_mechanics.hpp @@ -35,7 +35,7 @@ namespace smith { template auto buildSolidMechanics(std::shared_ptr mesh, std::shared_ptr d_solid_nonlinear_solver, - smith::SecondOrderTimeIntegrationRule time_rule, size_t num_checkpoints, + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, size_t num_checkpoints, std::string physics_name, const std::vector& param_names = {}) { auto graph = std::make_shared(num_checkpoints); diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp index b40d58e84c..67e2b66a97 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.cpp @@ -16,7 +16,7 @@ namespace smith { SolidMechanicsStateAdvancer::SolidMechanicsStateAdvancer( std::shared_ptr solver, std::shared_ptr vector_bcs, std::shared_ptr solid_dynamic_weak_forms, - smith::SecondOrderTimeIntegrationRule time_rule) + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule) : solver_(solver), vector_bcs_(vector_bcs), solid_dynamic_weak_forms_(solid_dynamic_weak_forms), diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp index b81c514d80..622dc5a917 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp @@ -40,7 +40,7 @@ class SolidMechanicsStateAdvancer : public StateAdvancer { SolidMechanicsStateAdvancer(std::shared_ptr solid_solver, std::shared_ptr vector_bcs, std::shared_ptr solid_dynamic_weak_forms, - SecondOrderTimeIntegrationRule time_rule); + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule); /// State enum for indexing convenience enum STATE @@ -70,7 +70,7 @@ class SolidMechanicsStateAdvancer : public StateAdvancer { /// application you will get back: shape_disp, states, params, time, and solid_mechanics_weak_form template static auto buildWeakFormAndStates(const std::shared_ptr& mesh, const std::shared_ptr& graph, - SecondOrderTimeIntegrationRule time_rule, std::string physics_name, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, std::string physics_name, const std::vector& param_names, double initial_time = 0.0) { auto shape_disp = createFieldState(*graph, ShapeDispSpace{}, physics_name + "_shape_displacement", mesh->tag()); @@ -111,7 +111,7 @@ class SolidMechanicsStateAdvancer : public StateAdvancer { std::shared_ptr solid_dynamic_weak_forms_; ///< Solid mechanics time discretized weak forms, user must setup the appropriate ///< integrals. Has both the time discretized and the undiscretized weak forms. - SecondOrderTimeIntegrationRule time_rule_; ///< second order time integration rule. Can compute u, u_dot, u_dot_dot, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule_; ///< second order time integration rule. Can compute u, u_dot, u_dot_dot, ///< given the current predicted u and the previous u, u_dot, u_dot_dot }; diff --git a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp index a49e2a7377..c2a1436060 100644 --- a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp @@ -90,7 +90,7 @@ TEST_F(SolidMechanicsMeshFixture, TransientConstantGravity) std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - smith::SecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); auto [physics, solid_weak_form, bcs] = buildSolidMechanics( @@ -190,7 +190,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - smith::SecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); auto [physics, solid_weak_form, bcs] = buildSolidMechanics( @@ -317,7 +317,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - smith::SecondOrderTimeIntegrationRule time_rule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); auto [physics, solid_weak_form, bcs] = buildSolidMechanics( diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index bdbe1ab47a..c918436333 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -23,6 +23,15 @@ namespace smith { +/** + * @brief Compute Green's strain from the displacement gradient + */ +template +auto greenStrain(const tensor& grad_u) +{ + return 0.5 * (grad_u + transpose(grad_u) + dot(transpose(grad_u), grad_u)); +} + /// @brief Green-Saint Venant isotropic thermoelastic model struct GreenSaintVenantThermoelasticMaterial { double density; ///< density @@ -259,7 +268,6 @@ struct FieldStore { std::map to_unknown_index_; std::vector> boundary_conditions_; - // MRT, do this for readability struct FieldLabel { std::string field_name; size_t field_index_in_residual; @@ -320,18 +328,27 @@ std::vector solve(const std::vector& weak_forms, const Fi field_store.getBoundaryConditionManagers()); } -TEST_F(SolidMechanicsMeshFixture, A) +TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) { SMITH_MARK_FUNCTION; FieldType> shape_disp_type("shape_displacement"); - FieldType> disp_type("displacement"); FieldType> temperature_type("temperature"); std::shared_ptr d_nonlinear_solver = buildDifferentiableNonlinearBlockSolver(nonlinear_opts, linear_options, *mesh_); + double rho = 1.0; + double E = 100.0; + double nu = 0.25; + double c = 1.0; + double alpha = 1.0e-3; + double theta_ref = 0.0; + double k = 1.0; + GreenSaintVenantThermoelasticMaterial material{rho, E, nu, c, alpha, theta_ref, k}; + + FieldStore field_store(mesh_, 100); field_store.addShapeDisp(shape_disp_type); @@ -345,8 +362,8 @@ TEST_F(SolidMechanicsMeshFixture, A) disp_bc->setFixedVectorBCs(mesh_->domain("right")); auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); - auto velo_old_type = field_store.addDerived(disp_type, "velocity_old"); - auto accel_old_type = field_store.addDerived(disp_type, "acceleration_old"); + //auto velo_old_type = field_store.addDerived(disp_type, "velocity_old"); + //auto accel_old_type = field_store.addDerived(disp_type, "acceleration_old"); std::shared_ptr& temperature_bc = field_store.addUnknown(temperature_type); temperature_bc->setFixedScalarBCs(mesh_->domain("left")); @@ -354,30 +371,42 @@ TEST_F(SolidMechanicsMeshFixture, A) auto temperature_old_type = field_store.addDerived(temperature_type, "temperature_old"); - auto disp_time_rule = SecondOrderTimeIntegrationRule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + auto disp_time_rule = QuasiStaticFirstOrderTimeIntegrationRule(); auto temperature_time_rule = BackwardEulerFirstOrderTimeIntegrationRule(); auto solid_weak_form = createWeakForm("solid_force", disp_type, field_store, disp_type, disp_old_type, - velo_old_type, accel_old_type, temperature_type, temperature_old_type); + temperature_type, temperature_old_type); solid_weak_form->addBodyIntegral( - mesh_->entireBodyName(), [=](auto time_info, auto /*X*/, auto disp, auto disp_old, auto velo_old, auto accel_old, + mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto disp, auto disp_old, auto temperature, auto temperature_old) { - auto u = disp_time_rule.value(time_info, disp, disp_old, velo_old, accel_old); - auto a = disp_time_rule.ddot(time_info, disp, disp_old, velo_old, accel_old); - auto theta_dot = temperature_time_rule.dot(time_info, temperature, temperature_old); - return smith::tuple{2.0 * get(a), get(u) * get(theta_dot)}; + auto u = disp_time_rule.value(t_info, disp, disp_old); + auto v = disp_time_rule.dot(t_info, disp, disp_old); + auto T = temperature_time_rule.value(t_info, temperature, temperature_old); + GreenSaintVenantThermoelasticMaterial::State state; + auto [pk, C_v, s0, q0] = + material(t_info.dt(), state, get(u), get(v), get(T), get(T)); + return smith::tuple{smith::zero{}, pk}; }); auto thermal_weak_form = createWeakForm("thermal_flux", temperature_type, field_store, temperature_type, - temperature_old_type, disp_type); + temperature_old_type, disp_type, disp_old_type); + + thermal_weak_form->addBodyIntegral(mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto temperature, + auto temperature_old, auto disp, auto disp_old) { + GreenSaintVenantThermoelasticMaterial::State state; + auto u = disp_time_rule.value(t_info, disp, disp_old); + auto v = disp_time_rule.dot(t_info, disp, disp_old); + auto T = temperature_time_rule.value(t_info, temperature, temperature_old); + auto T_dot = temperature_time_rule.dot(t_info, temperature, temperature_old); + auto [pk, C_v, s0, q0] = + material(t_info.dt(), state, get(u), get(v), get(T), get(T)); + auto dT_dt = get(T_dot); + return smith::tuple{C_v * dT_dt - s0, -q0}; + }); - thermal_weak_form->addBodyIntegral( - mesh_->entireBodyName(), [=](auto time_info, auto /*X*/, auto temperature, auto temperature_old, auto) { - auto theta = temperature_time_rule.value(time_info, temperature, temperature_old); - auto theta_dot = temperature_time_rule.dot(time_info, temperature, temperature_old); - return smith::tuple{get(theta_dot), get(theta)}; - }); + thermal_weak_form->addBodySource(smith::DependsOn<>(), mesh_->entireBodyName(), + [](auto /*t*/, auto /* x */) { return 100.0; }); std::vector weak_forms{solid_weak_form.get(), thermal_weak_form.get()}; std::vector disp_temp = solve(weak_forms, field_store, d_nonlinear_solver.get(), TimeInfo(0.0, 1.0)); diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index 8cd5de4692..cc5496bfad 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -88,7 +88,7 @@ class SecondOrderTimeDiscretizedWeakForm; /// @brief Useful for time-discretized PDEs of second order (involves for first and second derivatives of time). Users /// write q-functions in terns of u, u_dot, u_dot_dot, and the weak form is transformed by the -/// SecondOrderTimeIntegrationRule so that is it globally a function of u, u_old, u_dot_old, u_dot_dot_old, with u as +/// ImplicitNewmarkSecondOrderTimeIntegrationRule so that is it globally a function of u, u_old, u_dot_old, u_dot_dot_old, with u as /// the distinct unknown for the time discretized system. /// @tparam spatial_dim Spatial dimension, 2 or 3. /// @tparam OutputSpace The space corresponding to the output residual for the weak form (test-space). @@ -108,7 +108,7 @@ class SecondOrderTimeDiscretizedWeakForm mesh, - SecondOrderTimeIntegrationRule time_rule, + ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule, const mfem::ParFiniteElementSpace& output_mfem_space, const typename TimeDiscretizedWeakFormT::SpacesT& input_mfem_spaces) : time_rule_(time_rule) @@ -157,7 +157,7 @@ class SecondOrderTimeDiscretizedWeakForm @@ -83,8 +79,6 @@ struct SecondOrderTimeIntegrationRule { auto dt = t.dt(); return (4.0 / (dt * dt)) * (field_new - field_old) - (4.0 / dt) * velo_old - accel_old; } - - SecondOrderTimeIntegrationMethod method_; ///< method specifying time integration rule to inject into the q-function. }; } // namespace smith From 59af3febaf403e1e74d6016db87be6012b2a919e Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 27 Jan 2026 14:50:50 -0800 Subject: [PATCH 16/17] Remove some duplication. --- .../test_solid_mechanics_state_advancer.cpp | 75 ++++--------------- .../tests/test_thermo_mechanics.cpp | 15 +--- .../time_integration_rule.hpp | 37 +++++++++ 3 files changed, 57 insertions(+), 70 deletions(-) diff --git a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp index c2a1436060..5c172179a2 100644 --- a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp @@ -90,7 +90,7 @@ TEST_F(SolidMechanicsMeshFixture, TransientConstantGravity) std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); + smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule; auto [physics, solid_weak_form, bcs] = buildSolidMechanics( @@ -181,20 +181,15 @@ TEST_F(SolidMechanicsMeshFixture, TransientConstantGravity) EXPECT_NEAR(0.0, u_err, 1e-14); } -TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) -{ - SMITH_MARK_FUNCTION; - - std::string physics_name = "solid"; +auto createSolidMechanicsBasePhysics(std::string physics_name, std::shared_ptr mesh) { + size_t num_checkpoints = 200; std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule(smith::SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); - auto [physics, solid_weak_form, bcs] = buildSolidMechanics( - mesh, d_solid_nonlinear_solver, time_rule, 200, physics_name, {"bulk", "shear"}); + mesh, d_solid_nonlinear_solver, ImplicitNewmarkSecondOrderTimeIntegrationRule(), num_checkpoints, physics_name, {"bulk", "shear"}); bcs->setFixedVectorBCs(mesh->domain("right")); bcs->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { @@ -221,7 +216,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) auto shape_disp = physics->getShapeDispFieldState(); auto params = physics->getFieldParams(); - auto initial_states = physics->getInitialFieldStates(); + auto states = physics->getInitialFieldStates(); params[0].get()->setFromFieldFunction([=](smith::tensor) { double scaling = 1.0; @@ -235,6 +230,16 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) physics->resetStates(); + return std::make_tuple(physics, shape_disp, states, params); +} + + +TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) +{ + SMITH_MARK_FUNCTION; + std::string physics_name = "solid"; + auto [physics, shape_disp, initial_states, params] = createSolidMechanicsBasePhysics(physics_name, mesh); + auto pv_writer = smith::createParaviewWriter(*mesh, physics->getFieldStatesAndParamStates(), physics_name); pv_writer.write(0, physics->time(), physics->getFieldStatesAndParamStates()); for (size_t m = 0; m < num_steps_; ++m) { @@ -311,56 +316,8 @@ void adjointBackward(std::shared_ptr physics, smith::FiniteElementD TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) { SMITH_MARK_FUNCTION; - std::string physics_name = "solid"; - - std::shared_ptr d_solid_nonlinear_solver = - buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); - - smith::ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule(SecondOrderTimeIntegrationMethod::IMPLICIT_NEWMARK); - - auto [physics, solid_weak_form, bcs] = - buildSolidMechanics( - mesh, d_solid_nonlinear_solver, time_rule, 200, physics_name, {"bulk", "shear"}); - - bcs->setFixedVectorBCs(mesh->domain("right")); - bcs->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { - auto bc = 0.0 * X; - bc[0] = 0.01 * t; - bc[1] = -0.05 * t; - return bc; - }); - - double E = 100.0; - double nu = 0.25; - auto K = E / (3.0 * (1.0 - 2 * nu)); - auto G = E / (2.0 * (1.0 + nu)); - using MaterialType = solid_mechanics::ParameterizedNeoHookeanSolid; - MaterialType material{.density = 10.0, .K0 = K, .G0 = G}; - - solid_weak_form->addBodyIntegral( - smith::DependsOn<0, 1>{}, mesh->entireBodyName(), - [material](const auto& /*time_info*/, auto /*X*/, auto u, auto /*v*/, auto a, auto bulk, auto shear) { - MaterialType::State state; - auto pk_stress = material(state, get(u), bulk, shear); - return smith::tuple{get(a) * material.density, pk_stress}; - }); - - auto shape_disp = physics->getShapeDispFieldState(); - auto params = physics->getFieldParams(); - auto states = physics->getInitialFieldStates(); - - params[0].get()->setFromFieldFunction([=](smith::tensor) { - double scaling = 1.0; - return scaling * material.K0; - }); - - params[1].get()->setFromFieldFunction([=](smith::tensor) { - double scaling = 1.0; - return scaling * material.G0; - }); - - physics->resetStates(); + auto [physics, shape_disp, initial_states, params] = createSolidMechanicsBasePhysics(physics_name, mesh); double qoi = integrateForward(physics, num_steps_, dt_); SLIC_INFO_ROOT(axom::fmt::format("{}", qoi)); diff --git a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp index c918436333..56eeb7407d 100644 --- a/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_thermo_mechanics.cpp @@ -348,7 +348,6 @@ TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) double k = 1.0; GreenSaintVenantThermoelasticMaterial material{rho, E, nu, c, alpha, theta_ref, k}; - FieldStore field_store(mesh_, 100); field_store.addShapeDisp(shape_disp_type); @@ -362,8 +361,6 @@ TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) disp_bc->setFixedVectorBCs(mesh_->domain("right")); auto disp_old_type = field_store.addDerived(disp_type, "displacement_old"); - //auto velo_old_type = field_store.addDerived(disp_type, "velocity_old"); - //auto accel_old_type = field_store.addDerived(disp_type, "acceleration_old"); std::shared_ptr& temperature_bc = field_store.addUnknown(temperature_type); temperature_bc->setFixedScalarBCs(mesh_->domain("left")); @@ -371,8 +368,8 @@ TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) auto temperature_old_type = field_store.addDerived(temperature_type, "temperature_old"); - auto disp_time_rule = QuasiStaticFirstOrderTimeIntegrationRule(); - auto temperature_time_rule = BackwardEulerFirstOrderTimeIntegrationRule(); + QuasiStaticFirstOrderTimeIntegrationRule disp_time_rule; + BackwardEulerFirstOrderTimeIntegrationRule temperature_time_rule; auto solid_weak_form = createWeakForm("solid_force", disp_type, field_store, disp_type, disp_old_type, temperature_type, temperature_old_type); @@ -411,14 +408,10 @@ TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) std::vector weak_forms{solid_weak_form.get(), thermal_weak_form.get()}; std::vector disp_temp = solve(weak_forms, field_store, d_nonlinear_solver.get(), TimeInfo(0.0, 1.0)); + + EXPECT_EQ(0, 0); - // auto derived_types = createDerivedFieldTypes(ddot(disp_type), dot(temperature_type)); - // FieldType> kappa_type; - // FieldType> mu_type; - // auto params = - // auto weak_form = - // createWeakForm("solid", smith::DependsOnFields(disp_type, dot(disp_type), ddot(disp_type), temperature_type)); } } // namespace smith diff --git a/src/smith/differentiable_numerics/time_integration_rule.hpp b/src/smith/differentiable_numerics/time_integration_rule.hpp index 1f4f97de19..2a48b47d85 100644 --- a/src/smith/differentiable_numerics/time_integration_rule.hpp +++ b/src/smith/differentiable_numerics/time_integration_rule.hpp @@ -81,4 +81,41 @@ struct ImplicitNewmarkSecondOrderTimeIntegrationRule { } }; + +/// @brief encodes rules for time discretizing second order odes (involving first and second time derivatives). +/// When solving f(u, u_dot, u_dot_dot, t) = 0 +/// this class provides the current discrete approximation for u, u_dot, and u_dot_dot as a function of +/// (u^{n+1},u^n,u_dot^n,u_dot_dot^n). +struct QuasiStaticSecondOrderTimeIntegrationRule { + /// @brief Constructor + QuasiStaticSecondOrderTimeIntegrationRule() {} + + /// @brief evaluate value of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto value([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return field_new; + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto dot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return (1.0 / t.dt()) * (field_new - field_old); + } + + /// @brief evaluate time derivative discretization of the ode state as used by the integration rule + template + SMITH_HOST_DEVICE auto ddot([[maybe_unused]] const TimeInfo& t, [[maybe_unused]] const T1& field_new, + [[maybe_unused]] const T2& field_old, [[maybe_unused]] const T3& velo_old, + [[maybe_unused]] const T4& accel_old) const + { + return accel_old; + } +}; + } // namespace smith From 956d6df1fa689de9449ba60c4a7e16cd099a7a4a Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Thu, 29 Jan 2026 12:08:08 -0800 Subject: [PATCH 17/17] Fix style. --- .../solid_mechanics_state_advancer.hpp | 5 ++- .../test_solid_mechanics_state_advancer.cpp | 10 ++--- .../tests/test_thermo_mechanics.cpp | 44 +++++++++++-------- .../time_discretized_weak_form.hpp | 4 +- .../time_integration_rule.hpp | 5 ++- src/smith/physics/solid_mechanics.hpp | 9 +++- 6 files changed, 47 insertions(+), 30 deletions(-) diff --git a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp index 622dc5a917..1218724912 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp +++ b/src/smith/differentiable_numerics/solid_mechanics_state_advancer.hpp @@ -111,8 +111,9 @@ class SolidMechanicsStateAdvancer : public StateAdvancer { std::shared_ptr solid_dynamic_weak_forms_; ///< Solid mechanics time discretized weak forms, user must setup the appropriate ///< integrals. Has both the time discretized and the undiscretized weak forms. - ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule_; ///< second order time integration rule. Can compute u, u_dot, u_dot_dot, - ///< given the current predicted u and the previous u, u_dot, u_dot_dot + ImplicitNewmarkSecondOrderTimeIntegrationRule + time_rule_; ///< second order time integration rule. Can compute u, u_dot, u_dot_dot, + ///< given the current predicted u and the previous u, u_dot, u_dot_dot }; } // namespace smith diff --git a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp index 5c172179a2..e22e123fbb 100644 --- a/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp +++ b/src/smith/differentiable_numerics/tests/test_solid_mechanics_state_advancer.cpp @@ -181,15 +181,16 @@ TEST_F(SolidMechanicsMeshFixture, TransientConstantGravity) EXPECT_NEAR(0.0, u_err, 1e-14); } - -auto createSolidMechanicsBasePhysics(std::string physics_name, std::shared_ptr mesh) { +auto createSolidMechanicsBasePhysics(std::string physics_name, std::shared_ptr mesh) +{ size_t num_checkpoints = 200; std::shared_ptr d_solid_nonlinear_solver = buildDifferentiableNonlinearSolver(solid_nonlinear_opts, solid_linear_options, *mesh); auto [physics, solid_weak_form, bcs] = buildSolidMechanics( - mesh, d_solid_nonlinear_solver, ImplicitNewmarkSecondOrderTimeIntegrationRule(), num_checkpoints, physics_name, {"bulk", "shear"}); + mesh, d_solid_nonlinear_solver, ImplicitNewmarkSecondOrderTimeIntegrationRule(), num_checkpoints, + physics_name, {"bulk", "shear"}); bcs->setFixedVectorBCs(mesh->domain("right")); bcs->setVectorBCs(mesh->domain("left"), [](double t, smith::tensor X) { @@ -233,7 +234,6 @@ auto createSolidMechanicsBasePhysics(std::string physics_name, std::shared_ptr& getAllFields() const { return fields_; } + std::vector getFields(const std::string& weak_form_name) const { auto unknown_field_indices = weak_form_name_to_field_indices_.at(weak_form_name); @@ -283,9 +295,8 @@ void createSpaces(const std::string& weak_form_name, FieldStore& field_store, std::vector& spaces, size_t arg_num, FirstType type, Types... types) { - const size_t test_index = field_store.getFieldIndex(type.name); SLIC_ERROR_IF(spaces.size() != arg_num, "Error creating spaces recursively"); - spaces.push_back(&field_store.getField(test_index).get()->space()); + spaces.push_back(&field_store.getField(type.name).get()->space()); field_store.addWeakFormArg(weak_form_name, type.name, arg_num); if (type.unknown_index >= 0) { field_store.addWeakFormUnknownArg(weak_form_name, type.name, arg_num); @@ -299,8 +310,7 @@ template auto createWeakForm(std::string name, FieldType test_type, FieldStore& field_store, FieldType... field_types) { - const size_t test_index = field_store.getFieldIndex(test_type.name); - const mfem::ParFiniteElementSpace& test_space = field_store.getField(test_index).get()->space(); + const mfem::ParFiniteElementSpace& test_space = field_store.getField(test_type.name).get()->space(); std::vector input_spaces; createSpaces(name, field_store, input_spaces, 0, field_types...); return std::make_shared>>( @@ -374,17 +384,16 @@ TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) auto solid_weak_form = createWeakForm("solid_force", disp_type, field_store, disp_type, disp_old_type, temperature_type, temperature_old_type); - solid_weak_form->addBodyIntegral( - mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto disp, auto disp_old, - auto temperature, auto temperature_old) { - auto u = disp_time_rule.value(t_info, disp, disp_old); - auto v = disp_time_rule.dot(t_info, disp, disp_old); - auto T = temperature_time_rule.value(t_info, temperature, temperature_old); - GreenSaintVenantThermoelasticMaterial::State state; - auto [pk, C_v, s0, q0] = - material(t_info.dt(), state, get(u), get(v), get(T), get(T)); - return smith::tuple{smith::zero{}, pk}; - }); + solid_weak_form->addBodyIntegral(mesh_->entireBodyName(), [=](auto t_info, auto /*X*/, auto disp, auto disp_old, + auto temperature, auto temperature_old) { + auto u = disp_time_rule.value(t_info, disp, disp_old); + auto v = disp_time_rule.dot(t_info, disp, disp_old); + auto T = temperature_time_rule.value(t_info, temperature, temperature_old); + GreenSaintVenantThermoelasticMaterial::State state; + auto [pk, C_v, s0, q0] = + material(t_info.dt(), state, get(u), get(v), get(T), get(T)); + return smith::tuple{smith::zero{}, pk}; + }); auto thermal_weak_form = createWeakForm("thermal_flux", temperature_type, field_store, temperature_type, temperature_old_type, disp_type, disp_old_type); @@ -408,10 +417,9 @@ TEST_F(SolidMechanicsMeshFixture, RunThermoMechanicalCoupled) std::vector weak_forms{solid_weak_form.get(), thermal_weak_form.get()}; std::vector disp_temp = solve(weak_forms, field_store, d_nonlinear_solver.get(), TimeInfo(0.0, 1.0)); - + // auto states = field_store.getFields(); EXPECT_EQ(0, 0); - } } // namespace smith diff --git a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp index cc5496bfad..a98c6d3074 100644 --- a/src/smith/differentiable_numerics/time_discretized_weak_form.hpp +++ b/src/smith/differentiable_numerics/time_discretized_weak_form.hpp @@ -88,8 +88,8 @@ class SecondOrderTimeDiscretizedWeakForm; /// @brief Useful for time-discretized PDEs of second order (involves for first and second derivatives of time). Users /// write q-functions in terns of u, u_dot, u_dot_dot, and the weak form is transformed by the -/// ImplicitNewmarkSecondOrderTimeIntegrationRule so that is it globally a function of u, u_old, u_dot_old, u_dot_dot_old, with u as -/// the distinct unknown for the time discretized system. +/// ImplicitNewmarkSecondOrderTimeIntegrationRule so that is it globally a function of u, u_old, u_dot_old, +/// u_dot_dot_old, with u as the distinct unknown for the time discretized system. /// @tparam spatial_dim Spatial dimension, 2 or 3. /// @tparam OutputSpace The space corresponding to the output residual for the weak form (test-space). /// @tparam TrialInputSpace The space corresponding to the predicted solution u, i.e., the trial solution, the unique diff --git a/src/smith/differentiable_numerics/time_integration_rule.hpp b/src/smith/differentiable_numerics/time_integration_rule.hpp index 2a48b47d85..79a78a67ef 100644 --- a/src/smith/differentiable_numerics/time_integration_rule.hpp +++ b/src/smith/differentiable_numerics/time_integration_rule.hpp @@ -41,7 +41,9 @@ struct BackwardEulerFirstOrderTimeIntegrationRule { }; /// Altenative name for Backward Euler which makes sense when restricting what are typically second order odes, -/// for example transient solid mechanics, to the quasi-static approximation. It happens that the implementation is identical to backward-Euler applied to first order systems as we want to be able to capture current velocity dependencies. +/// for example transient solid mechanics, to the quasi-static approximation. It happens that the implementation is +/// identical to backward-Euler applied to first order systems as we want to be able to capture current velocity +/// dependencies. using QuasiStaticFirstOrderTimeIntegrationRule = BackwardEulerFirstOrderTimeIntegrationRule; /// @brief encodes rules for time discretizing second order odes (involving first and second time derivatives). @@ -81,7 +83,6 @@ struct ImplicitNewmarkSecondOrderTimeIntegrationRule { } }; - /// @brief encodes rules for time discretizing second order odes (involving first and second time derivatives). /// When solving f(u, u_dot, u_dot_dot, t) = 0 /// this class provides the current discrete approximation for u, u_dot, and u_dot_dot as a function of diff --git a/src/smith/physics/solid_mechanics.hpp b/src/smith/physics/solid_mechanics.hpp index 988b1cceec..6ff77a66c6 100644 --- a/src/smith/physics/solid_mechanics.hpp +++ b/src/smith/physics/solid_mechanics.hpp @@ -1647,9 +1647,16 @@ class SolidMechanics, std::integer_se } if (use_warm_start_) { - // Update external forcing + // We want to solve using the forcing from the actual end of step, but, there are cases where the initial + // configuration is already out of equilibrium. In that case, we really want to consider the change in the + // forcing. For most runs, the previous state will be solved for equilibrium, so evaluating with displacement_, + // at time_ would give ~zero for unconstrained nodes. auto r = (*residual_)(time_ + dt, shapeDisplacement(), displacement_, acceleration_, *parameters_[parameter_indices].state...); + auto r_previous = (*residual_)(time_, shapeDisplacement(), displacement_, acceleration_, + *parameters_[parameter_indices].state...); + + r -= r_previous; // use the most recently evaluated Jacobian auto [_, drdu] = (*residual_)(time_, shapeDisplacement(), differentiate_wrt(displacement_), acceleration_,