From 33799534eb93aa697c1af564a0b9a1039a8c83ea Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 19 May 2026 08:24:27 -0700 Subject: [PATCH 01/28] initial implementation --- .../differentiable_numerics/CMakeLists.txt | 1 + .../plastic_mechanics_system.hpp | 269 ++++++++++++++++++ .../tests/CMakeLists.txt | 1 + .../tests/test_differentiable_plasticity.cpp | 150 ++++++++++ 4 files changed, 421 insertions(+) create mode 100644 src/smith/differentiable_numerics/plastic_mechanics_system.hpp create mode 100644 src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp diff --git a/src/smith/differentiable_numerics/CMakeLists.txt b/src/smith/differentiable_numerics/CMakeLists.txt index 91dc826899..864732bc5a 100644 --- a/src/smith/differentiable_numerics/CMakeLists.txt +++ b/src/smith/differentiable_numerics/CMakeLists.txt @@ -42,6 +42,7 @@ set(differentiable_numerics_headers thermo_mechanics_system.hpp system_base.hpp differentiable_test_utils.hpp + plastic_mechanics_system.hpp ) set(differentiable_numerics_depends diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp new file mode 100644 index 0000000000..36af17c2ae --- /dev/null +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -0,0 +1,269 @@ +// 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) + +/** + * @file electro_mechanics_system.hpp + * @brief Defines the ElectroMechanicsystem struct and its factory function + */ + +#pragma once + +#include "smith/differentiable_numerics/field_store.hpp" +#include "smith/differentiable_numerics/nonlinear_block_solver.hpp" +#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" +#include "smith/differentiable_numerics/multiphysics_time_integrator.hpp" +#include "smith/differentiable_numerics/time_integration_rule.hpp" +#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" +#include "smith/differentiable_numerics/differentiable_physics.hpp" +#include "smith/physics/weak_form.hpp" +#include "smith/differentiable_numerics/system_base.hpp" + +namespace smith { +/** + * @brief Container for a plastic solid mechanics system. + * @tparam dim Spatial dimension. + * @tparam disp_order Order of the displacement basis. + * @tparam parameter_space Finite element spaces for optional parameters. + * + * This system rely on a fixed-point iteration to solve plasticity equivalent to the traditional + * inner-outer Newton method. + * n \in [0, N] denotes the timestep and k \in [0, L] denotes the fixed-point iteration count. + * At every new timestep n+1, first predict new displacement by [u^{n+1,k}, Fp^{n+1,k}, epsilon_p^{n+1,k}] ---> u^{n+1,k+1}. + * Then update internal state variables [u^{n+1,k+1}, Fp^n, epsilon_p^n] ---> [Fp^{n+1,k+1}, epsilon_p^{n+1,k+1}]. + */ +template +struct PlasticMechanicsSystem : public SystemBase { + // Primary weak forms to solve for displacement + /// @brief using for SolidWeakFormType with inputs [u^{n+1, k+1}, Fp^{n+1, k}, epsilon_p^{n+1, k}] + using SolidWeakFormType = TimeDiscretizedWeakForm< + dim, H1, + Parameters, L2, L2, parameter_space...>>; + + /// @brief using for PlasticDeformWeakFormType with inputs [Fp^{n+1, k+1}, Fp^n, epsilon_p^{n+1, k+1}, epsilon_p^n, u^{n+1, k+1}] + using PlasticDeformWeakFormType = TimeDiscretizedWeakForm< + dim, L2, + Parameters, L2, L2, + L2, H1, parameter_space...>>; + + /// @brief using for PlasticStrainWeakFormType with inputs [epsilon_p^{n+1, k+1}, epsilon_p^n, Fp^{n+1, k+1}, Fp^n, u^{n+1, k+1}] + using PlasticStrainWeakFormType = TimeDiscretizedWeakForm< + dim, L2, + Parameters, L2, L2, L2 + H1, parameter_space...>>; + + // Primary weak forms + std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. + + // Internal variable weak forms to update in staggered solve + std::shared_ptr plastic_deform_weak_form; ///< Plastic deformation weak form. + std::shared_ptr plastic_strain_weak_form; ///< Plastic deformation weak form. + + // Primary variable bcs + std::shared_ptr disp_bc; ///< Displacement boundary conditions. + + // Internal variable bcs + std::shared_ptr plastic_deform_bc; + std::shared_ptr plastic_strain_bc; + + std::shared_ptr qusistatic_time_rule; + std::shared_ptr backward_euler_time_rule; + + // Helper functions for plastic deformation gradient bookkeeping + template + SMITH_HOST_DEVICE tensor recoverTensor(const tensor F_state) + { + tensor F_state_tensor{}; + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) { + F_state_tensor(i, j) = F_state(i * dim + j); + } + } + + return F_state_tensor; + } + + template + SMITH_HOST_DEVICE tensor flattenTensor(const tensor F_state) + { + tensor F_state_flattened{}; + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j) { + F_state_flattened(i * dim + j) = F_state(i, j); + } + } + + return F_state_flattened; + } + + /** + * @brief Get the list of all state fields. + * @return std::vector List of state fields. + */ + std::vector getStateFields() const + { + std::vector states; + states.push_back(field_store->getField(prefix("displacement"))); + + // Internal states + states.push_back(field_store->getField(prefix("plastic_defgrad"))); + states.push_back(field_store->getField(prefix("plastic_defgrad_old"))); + states.push_back(field_store->getField(prefix("plastic_strain"))); + states.push_back(field_store->getField(prefix("plastic_strain_old"))); + + return states; + } + + /** + * @brief Create a DifferentiablePhysics object for this system. + * @param physics_name The name of the physics. + * @return std::shared_ptr The differentiable physics object. + */ + std::shared_ptr createDifferentiablePhysics(std::string physics_name) + { + return std::make_shared( + field_store->getMesh(), field_store->graph(), field_store->getShapeDisp(), getStateFields(), + getParameterFields(), advancer, physics_name, + std::vector{prefix("solid_force")}); + } + + /** + * @brief Set the material model for a domain, defining integrals for weak form. + * @tparam MaterialType The material model type. + * @param material The material model instance. + * @param domain_name The name of the domain to apply the material to. + */ + template + void setMaterial(const std::string& domain_name, const MaterialType& mat) + { + solid_weak_form->addBodyIntegral(domain_name, + [&](auto /* t_info */, auto /* X */, auto u, auto Fp, auto epsilon_p, auto... params) { + auto du_dX = get(u); + auto Fp_tensor = this->recoverTensor(get(Fp)); + auto epsilon_pval = get(epsilon_p); + + auto P = mat.firstPiolaStress(du_dX, Fp_tensor, epsilon_pval, params...); + + return tuple{zero{}, P}; + }); + } + + /** + * @brief Set the plasticity update model for a domain, defining integrals for weak form. + * @tparam MaterialType The material model type. + * @param material The material model instance. + * @param domain_name The name of the domain to apply the material to. + */ + template + void setPlasticity(const std::string& domain_name, const MaterialType& mat) + { + auto captured_strain_rule = backward_euler_time_rule; + + plastic_deform_weak_form->addBodyIntegral(domain_name, + [=, this](auto t_info, auto /* X */, auto Fp, auto Fp_old, auto epsilon_p, auto epsilon_p_old, auto u, auto... params) { + auto du_dX = get(u); + auto Fp_tensor = this->recoverTensor(get(Fp)); + auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); + auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); + auto dt = t_info.dt(); + + auto Fp_predict_tensor = mat.plasticDeformGrad(dt, Fp_old_tensor, epsilon_dot, du_dX, params...); + auto Fp_predict = this->flattenTensor(Fp_predict_tensor); + + return tuple(get(Fp) - Fp_predict, zero{}); + }); + + plastic_strain_weak_form->addBodyIntegral(domain_name, + [=, this](auto t_info, auto /* X */, auto epsilon_p, auto epsilon_p_old, auto /* Fp */, auto Fp_old, auto u) { + auto du_dX = get(u); + auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); + auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); + auto dt = t_info.dt(); + + return tuple(epsilon_p_dot - mat.plasticStrain(dt, epsilon_current, epsilon_dot, Fp_old_tensor, du_dX), zero{}); + }); + } +}; + +/** + * @brief Factory function to build a chemomechanics system with L2 state variables. + */ +template +PlasticMechanicsSystem buildPlasticMechanicsSystem( + std::shared_ptr mesh, std::shared_ptr solver, + std::string prepend_name = "", FieldType... parameter_types) +{ + auto field_store = std::make_shared(mesh, 100); + + auto prefix = [&](const std::string& name) { + if (prepend_name.empty()) { + return name; + } + return prepend_name + "_" + name; + }; + + FieldType> shape_disp_type(prefix("shape_displacement")); + field_store->addShapeDisp(shape_disp_type); + + auto quasistatic_time_rule = std::make_shared(); + auto backward_euler_time_rule = std::make_shared(); + + // Displacement with quasi-static rule + FieldType> disp_type(prefix("displacement")); + auto disp_bc = field_store->addIndependent(disp_type, quasistatic_time_rule); + + // State variable fields + FieldType> plastic_defgrad_type(prefix("plastic_defgrad")); + auto plastic_defgrad_bc = field_store->addIndependent(plastic_defgrad_type, backward_euler_time_rule); + auto plastic_defgrad_old_type = + field_store->addDependent(plastic_defgrad_type, FieldStore::TimeDerivative::VAL, prefix("plastic_defgrad_old")); + + FieldType> plastic_strain_type(prefix("plastic_strain")); + auto plastic_strain_bc = field_store->addIndependent(plastic_strain_type, backward_euler_time_rule); + auto plastic_strain_old_type = + field_store->addDependent(plastic_strain_type, FieldStore::TimeDerivative::VAL, prefix("plastic_strain_old")); + + // Parameters + std::vector parameter_fields; + (field_store->addParameter(FieldType(prefix("param_" + parameter_types.name))), ...); + (parameter_fields.push_back(field_store->getField(prefix("param_" + parameter_types.name))), ...); + + using SystemType = PlasticMechanicsSystem; + + // Main weak forms for mechanics + std::string solid_res_name = prefix("solid_residual"); + auto solid_weak_form = std::make_shared( + solid_res_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), + field_store->createSpaces(solid_res_name, disp_type.name, disp_type, + plastic_defgrad_type, plastic_strain_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + std::string plastic_defgrad_res_name = prefix("plastic_defgrad_residual"); + auto plastic_defgrad_weak_form = std::make_shared( + plastic_defgrad_res_name, field_store->getMesh(), field_store->getField(plastic_defgrad_type.name).get()->space(), + field_store->createSpaces(plastic_defgrad_res_name, plastic_defgrad_type.name, plastic_defgrad_type, + plastic_defgrad_old_type, plastic_strain_type, plastic_strain_old_type, disp_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + std::string plastic_strain_res_name = prefix("plastic_defgrad_residual"); + auto plastic_strain_weak_form = std::make_shared( + plastic_strain_res_name, field_store->getMesh(), field_store->getField(plastic_strain_type.name).get()->space(), + field_store->createSpaces(plastic_strain_res_name, plastic_strain_type.name, plastic_strain_type, + plastic_strain_old_type, plastic_defgrad_type, disp_type, + FieldType(prefix("param_" + parameter_types.name))...)); + + // Build solver and advancer + std::vector> weak_forms{ + solid_weak_form, plastic_defgrad_weak_form, plastic_strain_weak_form}; + auto advancer = std::make_shared(field_store, weak_forms, solver); + + return SystemType{{field_store, solver, advancer, parameter_fields, prepend_name}, + solid_weak_form, plastic_defgrad_weak_form, plastic_strain_weak_form, + disp_bc, plastic_defgrad_bc, plastic_strain_bc, quasistatic_time_rule, + backward_euler_time_rule}; +} +} // namespace smith diff --git a/src/smith/differentiable_numerics/tests/CMakeLists.txt b/src/smith/differentiable_numerics/tests/CMakeLists.txt index bbe95cbd55..7c6008bda3 100644 --- a/src/smith/differentiable_numerics/tests/CMakeLists.txt +++ b/src/smith/differentiable_numerics/tests/CMakeLists.txt @@ -16,6 +16,7 @@ set(differentiable_numerics_test_source test_thermal_static.cpp test_solid_static_with_internal_vars.cpp test_multiphysics_time_integrator.cpp + test_differentiable_plasticity.cpp ) smith_add_tests( SOURCES ${differentiable_numerics_test_source} diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp new file mode 100644 index 0000000000..c60a27bb02 --- /dev/null +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -0,0 +1,150 @@ +// 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 +#include "gtest/gtest.h" + +#include "smith/smith_config.hpp" +#include "smith/smith.hpp" + +namespace smith { + +/// @brief Differentiable J2 material with nonlinear isotropic hardening and linear kinematic hardening +template +class DifferentiableJ2SmallStrain { + public: + DifferentiableJ2SmallStrain(HardeningType hardeningModel, double youngsModulus, double poissonsRatio, + double hardeningModulus, double density) + : hardening(hardeningModel), E(youngsModulus), nu(poissonsRatio), Hk(hardeningModulus), rho(density) + {} + + /** @brief calculate the first Piola stress, given the displacement gradient and previous staggered solve material state */ + template + SMITH_HOST_DEVICE auto firstPiolaStress(const tensor& du_dX, const tensor& Fp, + const T3& /* epsilon_p */) const + { + constexpr auto I = Identity(); + const double K = E / (3.0 * (1.0 - 2.0 * nu)); + const double G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp; + auto p = K * tr(el_strain); + auto s = 2.0 * G * dev(el_strain); + + return s + p * I; + } + + /** @brief calculate the plastic deformation gradient */ + // Hanyu : combine plasticDeformGrad and plasticStrain + template + SMITH_HOST_DEVICE auto plasticDeformGrad(double dt, const tensor& Fp_old, + const T2& epsilon_dot, const tensor& du_dX) const + { + using std::sqrt; + constexpr auto I = Identity(); + const double K = E / (3.0 * (1.0 - 2.0 * nu)); + const double G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp_old; + auto p = K * tr(el_strain); + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + auto Np = 1.5 * eta / q; + + return Fp_old + epsilon_dot * dt * Np; + } + + /** @brief calculate the plastic strain */ + template + SMITH_HOST_DEVICE auto plasticStrain(double dt, const T1& epsilon_p, const T2& epsilon_dot, const tensor& Fp_old, + const tensor& du_dX) const + { + using std::sqrt; + constexpr auto I = Identity(); + const double K = E / (3.0 * (1.0 - 2.0 * nu)); + const double G = 0.5 * E / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp_old; + auto p = K * tr(el_strain); + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + return q / dt - (3.0 * G + Hk) * epsilon_dot - this->hardening(epsilon_p, epsilon_dot) / dt; + } + + private: + HardeningType hardening; ///< Flow stress hardening model + double E; ///< Young's modulus + double nu; ///< Poisson's ratio + double Hk; ///< Kinematic hardening modulus + double rho; ///< Mass density +}; + +smith::LinearSolverOptions primal_lin_opts{.linear_solver = smith::LinearSolver::CG, + .preconditioner = smith::Preconditioner::HypreAMG, + .relative_tol = 1e-6, + .absolute_tol = 1e-10, + .max_iterations = 200, + .print_level = 0}; + +smith::NonlinearSolverOptions primal_nonlin_opts{.nonlin_solver = NonlinearSolver::TrustRegion, + .relative_tol = 1.0e-6, + .absolute_tol = 1.0e-8, + .max_iterations = 25, + .print_level = 1}; + +smith::LinearSolverOptions state_lin_opts{.linear_solver = smith::LinearSolver::GMRES, + .preconditioner = smith::Preconditioner::HypreAMG, + .relative_tol = 1e-6, + .absolute_tol = 1e-10, + .max_iterations = 200, + .print_level = 0}; + +smith::NonlinearSolverOptions state_nonlin_opts{.nonlin_solver = NonlinearSolver::NewtonLineSearch, + .relative_tol = 1.0e-6, + .absolute_tol = 1.0e-8, + .max_iterations = 25, + .print_level = 1}; + +TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) +{ + MPI_Barrier(MPI_COMM_WORLD); + + int serial_refinement = 1; + int parallel_refinement = 1; + + static constexpr int dim = 3; + static constexpr int order = 1; + + axom::sidre::DataStore datastore; + smith::StateManager::initialize(datastore, "plasticity_small_strain"); + + std::string filename = SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"; + const std::string meshtag = "mesh"; + auto mesh = std::make_shared(smith::buildMeshFromFile(filename), meshtag, serial_refinement, + parallel_refinement); + + auto staggered_coupled_solver = std::make_shared(10); + auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); + auto state_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + staggered_coupled_solver->addSubsystemSolver({0}, primal_block_solver); + staggered_coupled_solver->addSubsystemSolver({1, 2}, state_block_solver); + + auto plastic_mechanics_system = buildPlasticMechanicsSystem(mesh, staggered_coupled_solver); + + using Hardening = solid_mechanics::LinearHardening; + Hardening hardening{.sigma_y = 50.0, .Hi = 50.0, .eta = 0.0}; + DifferentiableJ2SmallStrain mat(hardening, 1e+4, 0.25, 5.0, 1.0); + + plastic_mechanics_system.setMaterial(mesh->entireBodyName(), mat); +} + +} // namespace smith From 4b4aff65957ab16d8b2457883d3440a83b961f4f Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 19 May 2026 13:47:59 -0700 Subject: [PATCH 02/28] Solver runs --- .../coupled_system_solver.cpp | 6 ++ .../plastic_mechanics_system.hpp | 13 ++-- .../tests/test_differentiable_plasticity.cpp | 77 +++++++++++++------ 3 files changed, 67 insertions(+), 29 deletions(-) diff --git a/src/smith/differentiable_numerics/coupled_system_solver.cpp b/src/smith/differentiable_numerics/coupled_system_solver.cpp index 85599fa745..2b899a5ad0 100644 --- a/src/smith/differentiable_numerics/coupled_system_solver.cpp +++ b/src/smith/differentiable_numerics/coupled_system_solver.cpp @@ -102,8 +102,14 @@ std::vector CoupledSystemSolver::solve( } for (int iter = 0; iter < max_staggered_iterations_; ++iter) { + SLIC_INFO_ROOT("Solving staggered iteration " << iter); + smith::logger::flush(); + // --- Run each stage --- for (size_t stage_idx = 0; stage_idx < active_stages.size(); ++stage_idx) { + SLIC_INFO_ROOT("Solving stage " << stage_idx); + smith::logger::flush(); + const auto& stage = active_stages[stage_idx]; size_t num_stage_blocks = stage.block_indices.size(); diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 36af17c2ae..97f10e9770 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -51,7 +51,7 @@ struct PlasticMechanicsSystem : public SystemBase { /// @brief using for PlasticStrainWeakFormType with inputs [epsilon_p^{n+1, k+1}, epsilon_p^n, Fp^{n+1, k+1}, Fp^n, u^{n+1, k+1}] using PlasticStrainWeakFormType = TimeDiscretizedWeakForm< dim, L2, - Parameters, L2, L2, L2 + Parameters, L2, L2, L2, H1, parameter_space...>>; // Primary weak forms @@ -166,12 +166,11 @@ struct PlasticMechanicsSystem : public SystemBase { plastic_deform_weak_form->addBodyIntegral(domain_name, [=, this](auto t_info, auto /* X */, auto Fp, auto Fp_old, auto epsilon_p, auto epsilon_p_old, auto u, auto... params) { auto du_dX = get(u); - auto Fp_tensor = this->recoverTensor(get(Fp)); auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); auto dt = t_info.dt(); - auto Fp_predict_tensor = mat.plasticDeformGrad(dt, Fp_old_tensor, epsilon_dot, du_dX, params...); + auto Fp_predict_tensor = mat.plasticDeformGrad(dt, Fp_old_tensor, get(epsilon_dot), du_dX, params...); auto Fp_predict = this->flattenTensor(Fp_predict_tensor); return tuple(get(Fp) - Fp_predict, zero{}); @@ -184,7 +183,7 @@ struct PlasticMechanicsSystem : public SystemBase { auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); auto dt = t_info.dt(); - return tuple(epsilon_p_dot - mat.plasticStrain(dt, epsilon_current, epsilon_dot, Fp_old_tensor, du_dX), zero{}); + return tuple(get(epsilon_dot) - mat.plasticStrain(dt, get(epsilon_current), get(epsilon_dot), Fp_old_tensor, du_dX), zero{}); }); } }; @@ -222,7 +221,7 @@ PlasticMechanicsSystem buildPlasticMechanic auto plastic_defgrad_old_type = field_store->addDependent(plastic_defgrad_type, FieldStore::TimeDerivative::VAL, prefix("plastic_defgrad_old")); - FieldType> plastic_strain_type(prefix("plastic_strain")); + FieldType> plastic_strain_type(prefix("plastic_strain")); auto plastic_strain_bc = field_store->addIndependent(plastic_strain_type, backward_euler_time_rule); auto plastic_strain_old_type = field_store->addDependent(plastic_strain_type, FieldStore::TimeDerivative::VAL, prefix("plastic_strain_old")); @@ -249,11 +248,11 @@ PlasticMechanicsSystem buildPlasticMechanic plastic_defgrad_old_type, plastic_strain_type, plastic_strain_old_type, disp_type, FieldType(prefix("param_" + parameter_types.name))...)); - std::string plastic_strain_res_name = prefix("plastic_defgrad_residual"); + std::string plastic_strain_res_name = prefix("plastic_strain_residual"); auto plastic_strain_weak_form = std::make_shared( plastic_strain_res_name, field_store->getMesh(), field_store->getField(plastic_strain_type.name).get()->space(), field_store->createSpaces(plastic_strain_res_name, plastic_strain_type.name, plastic_strain_type, - plastic_strain_old_type, plastic_defgrad_type, disp_type, + plastic_strain_old_type, plastic_defgrad_type, plastic_defgrad_old_type, disp_type, FieldType(prefix("param_" + parameter_types.name))...)); // Build solver and advancer diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index c60a27bb02..66fd9bb487 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -44,12 +44,9 @@ class DifferentiableJ2SmallStrain { const T2& epsilon_dot, const tensor& du_dX) const { using std::sqrt; - constexpr auto I = Identity(); - const double K = E / (3.0 * (1.0 - 2.0 * nu)); const double G = 0.5 * E / (1.0 + nu); auto el_strain = sym(du_dX) - Fp_old; - auto p = K * tr(el_strain); auto s = 2.0 * G * dev(el_strain); auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; auto eta = s - sigma_b; @@ -66,12 +63,9 @@ class DifferentiableJ2SmallStrain { const tensor& du_dX) const { using std::sqrt; - constexpr auto I = Identity(); - const double K = E / (3.0 * (1.0 - 2.0 * nu)); const double G = 0.5 * E / (1.0 + nu); auto el_strain = sym(du_dX) - Fp_old; - auto p = K * tr(el_strain); auto s = 2.0 * G * dev(el_strain); auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; auto eta = s - sigma_b; @@ -90,44 +84,44 @@ class DifferentiableJ2SmallStrain { smith::LinearSolverOptions primal_lin_opts{.linear_solver = smith::LinearSolver::CG, .preconditioner = smith::Preconditioner::HypreAMG, - .relative_tol = 1e-6, - .absolute_tol = 1e-10, + .relative_tol = 1e-12, + .absolute_tol = 1e-12, .max_iterations = 200, .print_level = 0}; -smith::NonlinearSolverOptions primal_nonlin_opts{.nonlin_solver = NonlinearSolver::TrustRegion, - .relative_tol = 1.0e-6, - .absolute_tol = 1.0e-8, +smith::NonlinearSolverOptions primal_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1.0e-12, + .absolute_tol = 1.0e-12, .max_iterations = 25, - .print_level = 1}; + .print_level = 2}; smith::LinearSolverOptions state_lin_opts{.linear_solver = smith::LinearSolver::GMRES, .preconditioner = smith::Preconditioner::HypreAMG, - .relative_tol = 1e-6, + .relative_tol = 1e-10, .absolute_tol = 1e-10, .max_iterations = 200, .print_level = 0}; -smith::NonlinearSolverOptions state_nonlin_opts{.nonlin_solver = NonlinearSolver::NewtonLineSearch, - .relative_tol = 1.0e-6, - .absolute_tol = 1.0e-8, +smith::NonlinearSolverOptions state_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1.0e-10, + .absolute_tol = 1.0e-10, .max_iterations = 25, - .print_level = 1}; + .print_level = 2}; TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) { MPI_Barrier(MPI_COMM_WORLD); - int serial_refinement = 1; - int parallel_refinement = 1; + int serial_refinement = 0; + int parallel_refinement = 0; static constexpr int dim = 3; - static constexpr int order = 1; + static constexpr int order = 2; axom::sidre::DataStore datastore; smith::StateManager::initialize(datastore, "plasticity_small_strain"); - std::string filename = SMITH_REPO_DIR "/data/meshes/patch3D_tets.mesh"; + std::string filename = SMITH_REPO_DIR "/data/meshes/beam-hex.mesh"; const std::string meshtag = "mesh"; auto mesh = std::make_shared(smith::buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); @@ -141,10 +135,49 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) auto plastic_mechanics_system = buildPlasticMechanicsSystem(mesh, staggered_coupled_solver); using Hardening = solid_mechanics::LinearHardening; - Hardening hardening{.sigma_y = 50.0, .Hi = 50.0, .eta = 0.0}; + Hardening hardening{.sigma_y = 500.0, .Hi = 0.0, .eta = 0.5}; DifferentiableJ2SmallStrain mat(hardening, 1e+4, 0.25, 5.0, 1.0); plastic_mechanics_system.setMaterial(mesh->entireBodyName(), mat); + plastic_mechanics_system.setPlasticity(mesh->entireBodyName(), mat); + + // prescribe zero displacement at the supported end of the beam, + mesh->addDomainOfBoundaryElements("support", by_attr(1)); + plastic_mechanics_system.disp_bc->setFixedVectorBCs(mesh->domain("support")); + + // apply a displacement along z to the the tip of the beam + mesh->addDomainOfBoundaryElements("tip", by_attr(2)); + auto translated_in_z = [](double t, tensor) { + tensor u{}; + u[2] = t * (t - 1); + return u; + }; + plastic_mechanics_system.disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); + + double dt = 0.1; + double time = 0.0; + auto shape_disp = plastic_mechanics_system.field_store->getShapeDisp(); + auto states = plastic_mechanics_system.getStateFields(); + auto params = plastic_mechanics_system.getParameterFields(); + + auto pv_writer = createParaviewWriter(*mesh, states, "J2_small_strain"); + pv_writer.write(0, 0.0, states); + + std::vector reactions; + for (size_t step = 0; step < 1; ++step) { + std::tie(states, reactions) = + plastic_mechanics_system.advancer->advanceState(smith::TimeInfo(time, dt, step), shape_disp, states, params); + time += dt; + pv_writer.write(step + 1, time, states); + } + } } // namespace smith + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + smith::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} From dee021c115f47b21b2cb64f58d11426c6d23e179 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 27 May 2026 05:12:42 -0700 Subject: [PATCH 03/28] plastic solve with staggered solver --- .../plastic_mechanics_system.hpp | 27 ++++++------ .../tests/test_differentiable_plasticity.cpp | 43 +++++++++++-------- 2 files changed, 38 insertions(+), 32 deletions(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 97f10e9770..b4b8e0d8a6 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -40,18 +40,18 @@ struct PlasticMechanicsSystem : public SystemBase { /// @brief using for SolidWeakFormType with inputs [u^{n+1, k+1}, Fp^{n+1, k}, epsilon_p^{n+1, k}] using SolidWeakFormType = TimeDiscretizedWeakForm< dim, H1, - Parameters, L2, L2, parameter_space...>>; + Parameters, L2, L2, parameter_space...>>; /// @brief using for PlasticDeformWeakFormType with inputs [Fp^{n+1, k+1}, Fp^n, epsilon_p^{n+1, k+1}, epsilon_p^n, u^{n+1, k+1}] using PlasticDeformWeakFormType = TimeDiscretizedWeakForm< - dim, L2, - Parameters, L2, L2, - L2, H1, parameter_space...>>; + dim, L2, + Parameters, L2, L2, + L2, H1, parameter_space...>>; /// @brief using for PlasticStrainWeakFormType with inputs [epsilon_p^{n+1, k+1}, epsilon_p^n, Fp^{n+1, k+1}, Fp^n, u^{n+1, k+1}] using PlasticStrainWeakFormType = TimeDiscretizedWeakForm< - dim, L2, - Parameters, L2, L2, L2, + dim, L2, + Parameters, L2, L2, L2, H1, parameter_space...>>; // Primary weak forms @@ -141,12 +141,12 @@ struct PlasticMechanicsSystem : public SystemBase { void setMaterial(const std::string& domain_name, const MaterialType& mat) { solid_weak_form->addBodyIntegral(domain_name, - [&](auto /* t_info */, auto /* X */, auto u, auto Fp, auto epsilon_p, auto... params) { + [&](auto t_info, auto /* X */, auto u, auto Fp, auto epsilon_p, auto... params) { auto du_dX = get(u); auto Fp_tensor = this->recoverTensor(get(Fp)); - auto epsilon_pval = get(epsilon_p); + auto dt = t_info.dt(); - auto P = mat.firstPiolaStress(du_dX, Fp_tensor, epsilon_pval, params...); + auto P = mat.firstPiolaStress(dt, du_dX, Fp_tensor, get(epsilon_p), params...); return tuple{zero{}, P}; }); @@ -177,13 +177,14 @@ struct PlasticMechanicsSystem : public SystemBase { }); plastic_strain_weak_form->addBodyIntegral(domain_name, - [=, this](auto t_info, auto /* X */, auto epsilon_p, auto epsilon_p_old, auto /* Fp */, auto Fp_old, auto u) { + [=, this](auto t_info, auto /* X */, auto epsilon_p, auto epsilon_p_old, auto /* Fp */, auto Fp_old, auto u, auto... params) { auto du_dX = get(u); auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); auto dt = t_info.dt(); - return tuple(get(epsilon_dot) - mat.plasticStrain(dt, get(epsilon_current), get(epsilon_dot), Fp_old_tensor, du_dX), zero{}); + auto epsilon_dot_predict = mat.plasticStrain(dt, get(epsilon_current), get(epsilon_dot), Fp_old_tensor, du_dX); + return tuple(get(epsilon_dot) - epsilon_dot_predict, zero{}); }); } }; @@ -216,12 +217,12 @@ PlasticMechanicsSystem buildPlasticMechanic auto disp_bc = field_store->addIndependent(disp_type, quasistatic_time_rule); // State variable fields - FieldType> plastic_defgrad_type(prefix("plastic_defgrad")); + FieldType> plastic_defgrad_type(prefix("plastic_defgrad")); auto plastic_defgrad_bc = field_store->addIndependent(plastic_defgrad_type, backward_euler_time_rule); auto plastic_defgrad_old_type = field_store->addDependent(plastic_defgrad_type, FieldStore::TimeDerivative::VAL, prefix("plastic_defgrad_old")); - FieldType> plastic_strain_type(prefix("plastic_strain")); + FieldType> plastic_strain_type(prefix("plastic_strain")); auto plastic_strain_bc = field_store->addIndependent(plastic_strain_type, backward_euler_time_rule); auto plastic_strain_old_type = field_store->addDependent(plastic_strain_type, FieldStore::TimeDerivative::VAL, prefix("plastic_strain_old")); diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index 66fd9bb487..59653f230a 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -17,13 +17,13 @@ template class DifferentiableJ2SmallStrain { public: DifferentiableJ2SmallStrain(HardeningType hardeningModel, double youngsModulus, double poissonsRatio, - double hardeningModulus, double density) - : hardening(hardeningModel), E(youngsModulus), nu(poissonsRatio), Hk(hardeningModulus), rho(density) + double hardeningModulus, double density, double returnMapTol = 1e-10) + : hardening(hardeningModel), E(youngsModulus), nu(poissonsRatio), Hk(hardeningModulus), rho(density), tol(returnMapTol) {} /** @brief calculate the first Piola stress, given the displacement gradient and previous staggered solve material state */ template - SMITH_HOST_DEVICE auto firstPiolaStress(const tensor& du_dX, const tensor& Fp, + SMITH_HOST_DEVICE auto firstPiolaStress(double /* dt */, const tensor& du_dX, const tensor& Fp, const T3& /* epsilon_p */) const { constexpr auto I = Identity(); @@ -38,7 +38,6 @@ class DifferentiableJ2SmallStrain { } /** @brief calculate the plastic deformation gradient */ - // Hanyu : combine plasticDeformGrad and plasticStrain template SMITH_HOST_DEVICE auto plasticDeformGrad(double dt, const tensor& Fp_old, const T2& epsilon_dot, const tensor& du_dX) const @@ -71,7 +70,9 @@ class DifferentiableJ2SmallStrain { auto eta = s - sigma_b; auto q = sqrt(1.5) * norm(eta); - return q / dt - (3.0 * G + Hk) * epsilon_dot - this->hardening(epsilon_p, epsilon_dot) / dt; + auto epsilon_dot_predict = q / dt - (3.0 * G + Hk) * epsilon_dot - this->hardening(epsilon_p, epsilon_dot) / dt; + auto non_negativity = epsilon_dot_predict >= tol * hardening.sigma_y; + return epsilon_dot_predict * non_negativity; } private: @@ -80,31 +81,32 @@ class DifferentiableJ2SmallStrain { double nu; ///< Poisson's ratio double Hk; ///< Kinematic hardening modulus double rho; ///< Mass density + double tol; }; smith::LinearSolverOptions primal_lin_opts{.linear_solver = smith::LinearSolver::CG, .preconditioner = smith::Preconditioner::HypreAMG, - .relative_tol = 1e-12, - .absolute_tol = 1e-12, + .relative_tol = 1e-15, + .absolute_tol = 1e-15, .max_iterations = 200, .print_level = 0}; smith::NonlinearSolverOptions primal_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, - .relative_tol = 1.0e-12, - .absolute_tol = 1.0e-12, + .relative_tol = 1e-9, + .absolute_tol = 1e-10, .max_iterations = 25, .print_level = 2}; -smith::LinearSolverOptions state_lin_opts{.linear_solver = smith::LinearSolver::GMRES, - .preconditioner = smith::Preconditioner::HypreAMG, +smith::LinearSolverOptions state_lin_opts{.linear_solver = smith::LinearSolver::CG, + .preconditioner = smith::Preconditioner::HypreJacobi, .relative_tol = 1e-10, .absolute_tol = 1e-10, .max_iterations = 200, .print_level = 0}; smith::NonlinearSolverOptions state_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, - .relative_tol = 1.0e-10, - .absolute_tol = 1.0e-10, + .relative_tol = 1e-7, + .absolute_tol = 1e-8, .max_iterations = 25, .print_level = 2}; @@ -112,7 +114,7 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) { MPI_Barrier(MPI_COMM_WORLD); - int serial_refinement = 0; + int serial_refinement = 1; int parallel_refinement = 0; static constexpr int dim = 3; @@ -126,16 +128,19 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) auto mesh = std::make_shared(smith::buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); - auto staggered_coupled_solver = std::make_shared(10); + auto staggered_coupled_solver = std::make_shared(50); auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); - auto state_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + auto strain_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + auto defgrad_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + staggered_coupled_solver->addSubsystemSolver({0}, primal_block_solver); - staggered_coupled_solver->addSubsystemSolver({1, 2}, state_block_solver); + staggered_coupled_solver->addSubsystemSolver({2}, strain_block_solver); + staggered_coupled_solver->addSubsystemSolver({1}, defgrad_block_solver); auto plastic_mechanics_system = buildPlasticMechanicsSystem(mesh, staggered_coupled_solver); using Hardening = solid_mechanics::LinearHardening; - Hardening hardening{.sigma_y = 500.0, .Hi = 0.0, .eta = 0.5}; + Hardening hardening{.sigma_y = 40.0, .Hi = 50.0, .eta = 0.0}; DifferentiableJ2SmallStrain mat(hardening, 1e+4, 0.25, 5.0, 1.0); plastic_mechanics_system.setMaterial(mesh->entireBodyName(), mat); @@ -164,7 +169,7 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) pv_writer.write(0, 0.0, states); std::vector reactions; - for (size_t step = 0; step < 1; ++step) { + for (size_t step = 0; step < 10; ++step) { std::tie(states, reactions) = plastic_mechanics_system.advancer->advanceState(smith::TimeInfo(time, dt, step), shape_disp, states, params); time += dt; From 0b9292dda2891394c8b4e211921fc758709faa51 Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 14 Jul 2026 15:43:52 -0700 Subject: [PATCH 04/28] add finite difference test to differentiable plasticity --- .../coupled_system_solver.cpp | 16 +- .../differentiable_test_utils.hpp | 1 + .../nonlinear_block_solver.hpp | 6 + .../plastic_mechanics_system.hpp | 154 ++++++----- .../tests/test_differentiable_plasticity.cpp | 259 +++++++++++++++--- 5 files changed, 325 insertions(+), 111 deletions(-) diff --git a/src/smith/differentiable_numerics/coupled_system_solver.cpp b/src/smith/differentiable_numerics/coupled_system_solver.cpp index 2b899a5ad0..9ca0537dfe 100644 --- a/src/smith/differentiable_numerics/coupled_system_solver.cpp +++ b/src/smith/differentiable_numerics/coupled_system_solver.cpp @@ -65,6 +65,10 @@ std::vector CoupledSystemSolver::solve( stage.solver->setInnerToleranceMultiplier(inner_tol_factor); } + const bool print_staggered_progress = std::all_of(active_stages.begin(), active_stages.end(), [](const auto& stage) { + return stage.solver && stage.solver->printLevel() > 0; + }); + // Reset each stage solver's convergence tracking (e.g. initial residual norm for rel-tol) for (const auto& stage : active_stages) { stage.solver->resetConvergenceState(); @@ -102,13 +106,17 @@ std::vector CoupledSystemSolver::solve( } for (int iter = 0; iter < max_staggered_iterations_; ++iter) { - SLIC_INFO_ROOT("Solving staggered iteration " << iter); - smith::logger::flush(); + if (print_staggered_progress) { + SLIC_INFO_ROOT("Solving staggered iteration " << iter); + smith::logger::flush(); + } // --- Run each stage --- for (size_t stage_idx = 0; stage_idx < active_stages.size(); ++stage_idx) { - SLIC_INFO_ROOT("Solving stage " << stage_idx); - smith::logger::flush(); + if (print_staggered_progress) { + SLIC_INFO_ROOT("Solving stage " << stage_idx); + smith::logger::flush(); + } const auto& stage = active_stages[stage_idx]; size_t num_stage_blocks = stage.block_indices.size(); diff --git a/src/smith/differentiable_numerics/differentiable_test_utils.hpp b/src/smith/differentiable_numerics/differentiable_test_utils.hpp index 76b4e87700..e547945d97 100644 --- a/src/smith/differentiable_numerics/differentiable_test_utils.hpp +++ b/src/smith/differentiable_numerics/differentiable_test_utils.hpp @@ -136,6 +136,7 @@ inline double checkGradWrt(const gretl::State& objective, smith::FieldSt std::vector grad_errors; auto [grad, grad_fd] = checkGradients(objective, input, *dual_vec, objectiveBase, graph, eps); + if (printmore) std::cout << "grad = " << grad << "\ngrad fd = " << grad_fd << std::endl; grad_errors.push_back(std::abs(grad - grad_fd)); for (size_t step = 1; step < num_fd_steps; ++step) { diff --git a/src/smith/differentiable_numerics/nonlinear_block_solver.hpp b/src/smith/differentiable_numerics/nonlinear_block_solver.hpp index d75b2b1f0a..fd8dab0095 100644 --- a/src/smith/differentiable_numerics/nonlinear_block_solver.hpp +++ b/src/smith/differentiable_numerics/nonlinear_block_solver.hpp @@ -98,6 +98,9 @@ class NonlinearBlockSolverBase { /// @brief Set an inner-solve tolerance multiplier, e.g. for staggered solves. virtual void setInnerToleranceMultiplier(double multiplier) = 0; + /// @brief Return the configured solver print level. + virtual int printLevel() const { return 0; } + protected: mutable bool is_setup_ = false; ///< Records if this block solver has its preconditioner initialized. mutable NonlinearConvergenceContext convergence_context_ = {}; ///< Solver-owned convergence state for one solve. @@ -140,6 +143,9 @@ class NonlinearBlockSolver : public NonlinearBlockSolverBase { /// @brief Set the inner tolerance multiplier. void setInnerToleranceMultiplier(double multiplier) override { inner_tol_multiplier_ = multiplier; } + /// @brief Return the configured solver print level. + int printLevel() const override { return retained_nonlinear_options_ ? retained_nonlinear_options_->print_level : 0; } + /// @brief Build a fresh solver instance from retained config. std::shared_ptr cloneFresh() const; diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index b4b8e0d8a6..74ae9c4179 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -27,12 +27,13 @@ namespace smith { * @tparam dim Spatial dimension. * @tparam disp_order Order of the displacement basis. * @tparam parameter_space Finite element spaces for optional parameters. - * + * * This system rely on a fixed-point iteration to solve plasticity equivalent to the traditional * inner-outer Newton method. * n \in [0, N] denotes the timestep and k \in [0, L] denotes the fixed-point iteration count. - * At every new timestep n+1, first predict new displacement by [u^{n+1,k}, Fp^{n+1,k}, epsilon_p^{n+1,k}] ---> u^{n+1,k+1}. - * Then update internal state variables [u^{n+1,k+1}, Fp^n, epsilon_p^n] ---> [Fp^{n+1,k+1}, epsilon_p^{n+1,k+1}]. + * At every new timestep n+1, first predict new displacement by [u^{n+1,k}, Fp^{n+1,k}, epsilon_p^{n+1,k}] ---> + * u^{n+1,k+1}. Then update internal state variables [u^{n+1,k+1}, Fp^n, epsilon_p^n] ---> [Fp^{n+1,k+1}, + * epsilon_p^{n+1,k+1}]. */ template struct PlasticMechanicsSystem : public SystemBase { @@ -42,33 +43,35 @@ struct PlasticMechanicsSystem : public SystemBase { dim, H1, Parameters, L2, L2, parameter_space...>>; - /// @brief using for PlasticDeformWeakFormType with inputs [Fp^{n+1, k+1}, Fp^n, epsilon_p^{n+1, k+1}, epsilon_p^n, u^{n+1, k+1}] - using PlasticDeformWeakFormType = TimeDiscretizedWeakForm< - dim, L2, - Parameters, L2, L2, - L2, H1, parameter_space...>>; + /// @brief using for PlasticDeformWeakFormType with inputs [Fp^{n+1, k+1}, Fp^n, epsilon_p^{n+1, k+1}, epsilon_p^n, + /// u^{n+1, k+1}] + using PlasticDeformWeakFormType = + TimeDiscretizedWeakForm, + Parameters, L2, L2, + L2, H1, parameter_space...>>; - /// @brief using for PlasticStrainWeakFormType with inputs [epsilon_p^{n+1, k+1}, epsilon_p^n, Fp^{n+1, k+1}, Fp^n, u^{n+1, k+1}] - using PlasticStrainWeakFormType = TimeDiscretizedWeakForm< - dim, L2, - Parameters, L2, L2, L2, - H1, parameter_space...>>; + /// @brief using for PlasticStrainWeakFormType with inputs [epsilon_p^{n+1, k+1}, epsilon_p^n, Fp^{n+1, k+1}, Fp^n, + /// u^{n+1, k+1}] + using PlasticStrainWeakFormType = + TimeDiscretizedWeakForm, + Parameters, L2, L2, + L2, H1, parameter_space...>>; // Primary weak forms - std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. + std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. // Internal variable weak forms to update in staggered solve - std::shared_ptr plastic_deform_weak_form; ///< Plastic deformation weak form. - std::shared_ptr plastic_strain_weak_form; ///< Plastic deformation weak form. + std::shared_ptr plastic_deform_weak_form; ///< Plastic deformation weak form. + std::shared_ptr plastic_strain_weak_form; ///< Plastic deformation weak form. // Primary variable bcs - std::shared_ptr disp_bc; ///< Displacement boundary conditions. + std::shared_ptr disp_bc; ///< Displacement boundary conditions. // Internal variable bcs std::shared_ptr plastic_deform_bc; std::shared_ptr plastic_strain_bc; - std::shared_ptr qusistatic_time_rule; + std::shared_ptr qusistatic_time_rule; std::shared_ptr backward_euler_time_rule; // Helper functions for plastic deformation gradient bookkeeping @@ -118,6 +121,25 @@ struct PlasticMechanicsSystem : public SystemBase { return states; } + /** + * @brief Get the list of physical state fields for visualization. + * @return std::vector List of physical fields suitable for output. + */ + std::vector getOutputFieldStates() const + { + return {field_store->getField(prefix("displacement")), field_store->getField(prefix("plastic_defgrad")), + field_store->getField(prefix("plastic_strain"))}; + } + + /** + * @brief Get information about reaction fields for this system. + * @return List of ReactionInfo structures. + */ + std::vector getReactionInfos() const + { + return {{prefix("solid_force"), &field_store->getField(prefix("displacement")).get()->space()}}; + } + /** * @brief Create a DifferentiablePhysics object for this system. * @param physics_name The name of the physics. @@ -125,10 +147,9 @@ struct PlasticMechanicsSystem : public SystemBase { */ std::shared_ptr createDifferentiablePhysics(std::string physics_name) { - return std::make_shared( - field_store->getMesh(), field_store->graph(), field_store->getShapeDisp(), getStateFields(), - getParameterFields(), advancer, physics_name, - std::vector{prefix("solid_force")}); + return std::make_shared(field_store->getMesh(), field_store->graph(), + field_store->getShapeDisp(), getStateFields(), getParameterFields(), + advancer, physics_name, getReactionInfos()); } /** @@ -140,16 +161,16 @@ struct PlasticMechanicsSystem : public SystemBase { template void setMaterial(const std::string& domain_name, const MaterialType& mat) { - solid_weak_form->addBodyIntegral(domain_name, - [&](auto t_info, auto /* X */, auto u, auto Fp, auto epsilon_p, auto... params) { - auto du_dX = get(u); - auto Fp_tensor = this->recoverTensor(get(Fp)); - auto dt = t_info.dt(); + solid_weak_form->addBodyIntegral( + domain_name, [&](auto t_info, auto /* X */, auto u, auto Fp, auto epsilon_p, auto... params) { + auto du_dX = get(u); + auto Fp_tensor = this->recoverTensor(get(Fp)); + auto dt = t_info.dt(); - auto P = mat.firstPiolaStress(dt, du_dX, Fp_tensor, get(epsilon_p), params...); + auto P = mat.firstPiolaStress(dt, du_dX, Fp_tensor, get(epsilon_p), params...); - return tuple{zero{}, P}; - }); + return tuple{zero{}, P}; + }); } /** @@ -163,29 +184,32 @@ struct PlasticMechanicsSystem : public SystemBase { { auto captured_strain_rule = backward_euler_time_rule; - plastic_deform_weak_form->addBodyIntegral(domain_name, - [=, this](auto t_info, auto /* X */, auto Fp, auto Fp_old, auto epsilon_p, auto epsilon_p_old, auto u, auto... params) { - auto du_dX = get(u); - auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); - auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); - auto dt = t_info.dt(); - - auto Fp_predict_tensor = mat.plasticDeformGrad(dt, Fp_old_tensor, get(epsilon_dot), du_dX, params...); - auto Fp_predict = this->flattenTensor(Fp_predict_tensor); - - return tuple(get(Fp) - Fp_predict, zero{}); - }); - - plastic_strain_weak_form->addBodyIntegral(domain_name, - [=, this](auto t_info, auto /* X */, auto epsilon_p, auto epsilon_p_old, auto /* Fp */, auto Fp_old, auto u, auto... params) { - auto du_dX = get(u); - auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); - auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); - auto dt = t_info.dt(); - - auto epsilon_dot_predict = mat.plasticStrain(dt, get(epsilon_current), get(epsilon_dot), Fp_old_tensor, du_dX); - return tuple(get(epsilon_dot) - epsilon_dot_predict, zero{}); - }); + plastic_deform_weak_form->addBodyIntegral( + domain_name, [=, this](auto t_info, auto /* X */, auto Fp, auto Fp_old, auto epsilon_p, auto epsilon_p_old, + auto u, auto... params) { + auto du_dX = get(u); + auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); + auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); + auto dt = t_info.dt(); + + auto Fp_predict_tensor = mat.plasticDeformGrad(dt, Fp_old_tensor, get(epsilon_dot), du_dX, params...); + auto Fp_predict = this->flattenTensor(Fp_predict_tensor); + + return tuple(get(Fp) - Fp_predict, zero{}); + }); + + plastic_strain_weak_form->addBodyIntegral(domain_name, [=, this](auto t_info, auto /* X */, auto epsilon_p, + auto epsilon_p_old, auto /* Fp */, auto Fp_old, + auto u, auto... params) { + auto du_dX = get(u); + auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); + auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); + auto dt = t_info.dt(); + + auto epsilon_dot_predict = + mat.plasticStrain(dt, get(epsilon_current), get(epsilon_dot), Fp_old_tensor, du_dX, params...); + return tuple(get(epsilon_dot) - epsilon_dot_predict, zero{}); + }); } }; @@ -194,8 +218,8 @@ struct PlasticMechanicsSystem : public SystemBase { */ template PlasticMechanicsSystem buildPlasticMechanicsSystem( - std::shared_ptr mesh, std::shared_ptr solver, - std::string prepend_name = "", FieldType... parameter_types) + std::shared_ptr mesh, std::shared_ptr solver, std::string prepend_name = "", + FieldType... parameter_types) { auto field_store = std::make_shared(mesh, 100); @@ -238,8 +262,7 @@ PlasticMechanicsSystem buildPlasticMechanic std::string solid_res_name = prefix("solid_residual"); auto solid_weak_form = std::make_shared( solid_res_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), - field_store->createSpaces(solid_res_name, disp_type.name, disp_type, - plastic_defgrad_type, plastic_strain_type, + field_store->createSpaces(solid_res_name, disp_type.name, disp_type, plastic_defgrad_type, plastic_strain_type, FieldType(prefix("param_" + parameter_types.name))...)); std::string plastic_defgrad_res_name = prefix("plastic_defgrad_residual"); @@ -257,13 +280,18 @@ PlasticMechanicsSystem buildPlasticMechanic FieldType(prefix("param_" + parameter_types.name))...)); // Build solver and advancer - std::vector> weak_forms{ - solid_weak_form, plastic_defgrad_weak_form, plastic_strain_weak_form}; + std::vector> weak_forms{solid_weak_form, plastic_defgrad_weak_form, + plastic_strain_weak_form}; auto advancer = std::make_shared(field_store, weak_forms, solver); return SystemType{{field_store, solver, advancer, parameter_fields, prepend_name}, - solid_weak_form, plastic_defgrad_weak_form, plastic_strain_weak_form, - disp_bc, plastic_defgrad_bc, plastic_strain_bc, quasistatic_time_rule, - backward_euler_time_rule}; + solid_weak_form, + plastic_defgrad_weak_form, + plastic_strain_weak_form, + disp_bc, + plastic_defgrad_bc, + plastic_strain_bc, + quasistatic_time_rule, + backward_euler_time_rule}; } -} // namespace smith +} // namespace smith diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index 59653f230a..b53e106391 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -12,19 +12,47 @@ namespace smith { +LinearSolverOptions primal_lin_opts{.linear_solver = LinearSolver::CG, + .preconditioner = Preconditioner::HypreAMG, + .relative_tol = 1e-15, + .absolute_tol = 1e-15, + .max_iterations = 200, + .print_level = 0}; + +NonlinearSolverOptions primal_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1e-9, + .absolute_tol = 1e-10, + .max_iterations = 25, + .print_level = 2}; + +LinearSolverOptions state_lin_opts{.linear_solver = LinearSolver::CG, + .preconditioner = Preconditioner::HypreJacobi, + .relative_tol = 1e-10, + .absolute_tol = 1e-10, + .max_iterations = 200, + .print_level = 0}; + +NonlinearSolverOptions state_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1e-7, + .absolute_tol = 1e-8, + .max_iterations = 25, + .print_level = 2}; + /// @brief Differentiable J2 material with nonlinear isotropic hardening and linear kinematic hardening template class DifferentiableJ2SmallStrain { public: DifferentiableJ2SmallStrain(HardeningType hardeningModel, double youngsModulus, double poissonsRatio, - double hardeningModulus, double density, double returnMapTol = 1e-10) - : hardening(hardeningModel), E(youngsModulus), nu(poissonsRatio), Hk(hardeningModulus), rho(density), tol(returnMapTol) - {} + double hardeningModulus, double density) + : hardening(hardeningModel), E(youngsModulus), nu(poissonsRatio), Hk(hardeningModulus), rho(density) + { + } - /** @brief calculate the first Piola stress, given the displacement gradient and previous staggered solve material state */ + /** @brief calculate the first Piola stress, given the displacement gradient and previous staggered solve material + * state */ template - SMITH_HOST_DEVICE auto firstPiolaStress(double /* dt */, const tensor& du_dX, const tensor& Fp, - const T3& /* epsilon_p */) const + SMITH_HOST_DEVICE auto firstPiolaStress(double /* dt */, const tensor& du_dX, + const tensor& Fp, const T3& /* epsilon_p */) const { constexpr auto I = Identity(); const double K = E / (3.0 * (1.0 - 2.0 * nu)); @@ -39,8 +67,8 @@ class DifferentiableJ2SmallStrain { /** @brief calculate the plastic deformation gradient */ template - SMITH_HOST_DEVICE auto plasticDeformGrad(double dt, const tensor& Fp_old, - const T2& epsilon_dot, const tensor& du_dX) const + SMITH_HOST_DEVICE auto plasticDeformGrad(double dt, const tensor& Fp_old, const T2& epsilon_dot, + const tensor& du_dX) const { using std::sqrt; const double G = 0.5 * E / (1.0 + nu); @@ -58,8 +86,8 @@ class DifferentiableJ2SmallStrain { /** @brief calculate the plastic strain */ template - SMITH_HOST_DEVICE auto plasticStrain(double dt, const T1& epsilon_p, const T2& epsilon_dot, const tensor& Fp_old, - const tensor& du_dX) const + SMITH_HOST_DEVICE auto plasticStrain(double dt, const T1& epsilon_p, const T2& epsilon_dot, + const tensor& Fp_old, const tensor& du_dX) const { using std::sqrt; const double G = 0.5 * E / (1.0 + nu); @@ -71,7 +99,7 @@ class DifferentiableJ2SmallStrain { auto q = sqrt(1.5) * norm(eta); auto epsilon_dot_predict = q / dt - (3.0 * G + Hk) * epsilon_dot - this->hardening(epsilon_p, epsilon_dot) / dt; - auto non_negativity = epsilon_dot_predict >= tol * hardening.sigma_y; + auto non_negativity = epsilon_dot_predict >= 0.0; return epsilon_dot_predict * non_negativity; } @@ -81,52 +109,24 @@ class DifferentiableJ2SmallStrain { double nu; ///< Poisson's ratio double Hk; ///< Kinematic hardening modulus double rho; ///< Mass density - double tol; }; -smith::LinearSolverOptions primal_lin_opts{.linear_solver = smith::LinearSolver::CG, - .preconditioner = smith::Preconditioner::HypreAMG, - .relative_tol = 1e-15, - .absolute_tol = 1e-15, - .max_iterations = 200, - .print_level = 0}; - -smith::NonlinearSolverOptions primal_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, - .relative_tol = 1e-9, - .absolute_tol = 1e-10, - .max_iterations = 25, - .print_level = 2}; - -smith::LinearSolverOptions state_lin_opts{.linear_solver = smith::LinearSolver::CG, - .preconditioner = smith::Preconditioner::HypreJacobi, - .relative_tol = 1e-10, - .absolute_tol = 1e-10, - .max_iterations = 200, - .print_level = 0}; - -smith::NonlinearSolverOptions state_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, - .relative_tol = 1e-7, - .absolute_tol = 1e-8, - .max_iterations = 25, - .print_level = 2}; - TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) { MPI_Barrier(MPI_COMM_WORLD); - int serial_refinement = 1; + int serial_refinement = 0; int parallel_refinement = 0; static constexpr int dim = 3; static constexpr int order = 2; axom::sidre::DataStore datastore; - smith::StateManager::initialize(datastore, "plasticity_small_strain"); + StateManager::initialize(datastore, "plasticity_small_strain"); std::string filename = SMITH_REPO_DIR "/data/meshes/beam-hex.mesh"; const std::string meshtag = "mesh"; - auto mesh = std::make_shared(smith::buildMeshFromFile(filename), meshtag, serial_refinement, - parallel_refinement); + auto mesh = std::make_shared(buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); auto staggered_coupled_solver = std::make_shared(50); auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); @@ -171,14 +171,185 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) std::vector reactions; for (size_t step = 0; step < 10; ++step) { std::tie(states, reactions) = - plastic_mechanics_system.advancer->advanceState(smith::TimeInfo(time, dt, step), shape_disp, states, params); + plastic_mechanics_system.advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); time += dt; pv_writer.write(step + 1, time, states); } + auto displacement = plastic_mechanics_system.field_store->getField(plastic_mechanics_system.prefix("displacement")); + double final_disp_l2 = norm(*displacement.get(), 2.0); + EXPECT_GT(final_disp_l2, 0.0); +} + +/// @brief Parameterized differentiable J2 material with nonlinear isotropic hardening and linear kinematic hardening +template +class ParameterizedDifferentiableJ2SmallStrain { + public: + ParameterizedDifferentiableJ2SmallStrain(HardeningType hardeningModel, double baseYoungsModulus, double poissonsRatio, + double baseHardeningModulus, double density) + : hardening(hardeningModel), E0(baseYoungsModulus), nu(poissonsRatio), Hk0(baseHardeningModulus), rho(density) + { + } + + /** @brief calculate the first Piola stress, given the displacement gradient and previous staggered solve material + * state */ + template + SMITH_HOST_DEVICE auto firstPiolaStress(double /* dt */, const tensor& du_dX, + const tensor& Fp, const T3& /* epsilon_p */, const T4& deltaE, + const T5& /* deltaHk */) const + { + constexpr auto I = Identity(); + auto K = (E0 + get(deltaE)) / (3.0 * (1.0 - 2.0 * nu)); + auto G = 0.5 * (E0 + get(deltaE)) / (1.0 + nu); + + auto el_strain = sym(du_dX) - Fp; + auto p = K * tr(el_strain); + auto s = 2.0 * G * dev(el_strain); + + return s + p * I; + } + + /** @brief calculate the plastic deformation gradient */ + template + SMITH_HOST_DEVICE auto plasticDeformGrad(double dt, const tensor& Fp_old, const T2& epsilon_dot, + const tensor& du_dX, const T4& deltaE, const T5& deltaHk) const + { + using std::sqrt; + auto G = 0.5 * (E0 + get(deltaE)) / (1.0 + nu); + auto Hk = Hk0 + get(deltaHk); + + auto el_strain = sym(du_dX) - Fp_old; + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + auto Np = 1.5 * eta / q; + + return Fp_old + epsilon_dot * dt * Np; + } + + /** @brief calculate the plastic strain */ + template + SMITH_HOST_DEVICE auto plasticStrain(double dt, const T1& epsilon_p, const T2& epsilon_dot, + const tensor& Fp_old, const tensor& du_dX, + const T5& deltaE, const T6& deltaHk) const + { + using std::sqrt; + auto G = 0.5 * (E0 + get(deltaE)) / (1.0 + nu); + auto Hk = Hk0 + get(deltaHk); + + auto el_strain = sym(du_dX) - Fp_old; + auto s = 2.0 * G * dev(el_strain); + auto sigma_b = 2.0 / 3.0 * Hk * Fp_old; + auto eta = s - sigma_b; + auto q = sqrt(1.5) * norm(eta); + + auto epsilon_dot_predict = q / dt - (3.0 * G + Hk) * epsilon_dot - this->hardening(epsilon_p, epsilon_dot) / dt; + auto non_negativity = epsilon_dot_predict >= 0.0; + return epsilon_dot_predict * non_negativity; + } + + private: + HardeningType hardening; ///< Flow stress hardening model + double E0; ///< Base Young's modulus + double nu; ///< Poisson's ratio + double Hk0; ///< Base kinematic hardening modulus + double rho; ///< Mass density +}; + +TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) +{ + MPI_Barrier(MPI_COMM_WORLD); + + std::string physics_name = "differentiable_plasticity"; + + state_nonlin_opts.print_level = 0; + primal_nonlin_opts.print_level = 0; + + int serial_refinement = 0; + int parallel_refinement = 0; + + static constexpr int dim = 3; + static constexpr int order = 2; + + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "plasticity_small_strain"); + + std::string filename = SMITH_REPO_DIR "/data/meshes/beam-hex.mesh"; + const std::string meshtag = "mesh"; + auto mesh = std::make_shared(buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); + + auto staggered_coupled_solver = std::make_shared(100); + auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); + auto strain_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + auto defgrad_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); + + staggered_coupled_solver->addSubsystemSolver({0}, primal_block_solver); + staggered_coupled_solver->addSubsystemSolver({2}, strain_block_solver); + staggered_coupled_solver->addSubsystemSolver({1}, defgrad_block_solver); + + auto plastic_mechanics_system = buildPlasticMechanicsSystem( + mesh, staggered_coupled_solver, physics_name, FieldType>("deltaE"), FieldType>("deltaHk")); + + using Hardening = solid_mechanics::LinearHardening; + Hardening hardening{.sigma_y = 25.0, .Hi = 50.0, .eta = 0.0}; + ParameterizedDifferentiableJ2SmallStrain mat(hardening, 1.5e+4, 0.25, 500.0, 1.0); + + plastic_mechanics_system.setMaterial(mesh->entireBodyName(), mat); + plastic_mechanics_system.setPlasticity(mesh->entireBodyName(), mat); + + // prescribe zero displacement at the supported end of the beam, + mesh->addDomainOfBoundaryElements("support", by_attr(1)); + plastic_mechanics_system.disp_bc->setFixedVectorBCs(mesh->domain("support")); + + // apply a displacement along z to the the tip of the beam + mesh->addDomainOfBoundaryElements("tip", by_attr(2)); + auto translated_in_z = [](double t, tensor) { + tensor u{}; + u[2] = t; + return u; + }; + plastic_mechanics_system.disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); + + std::shared_ptr differentiable_plasticity = + plastic_mechanics_system.createDifferentiablePhysics(physics_name); + + auto pv_writer = createParaviewWriter(*mesh, differentiable_plasticity->getFieldStatesAndParamStates(), physics_name, + smith::ParaviewWriter::Options{.write_duals = false}); + pv_writer.write(0, differentiable_plasticity->time(), differentiable_plasticity->getFieldStatesAndParamStates()); + + double dt = 0.1; + for (size_t m = 0; m < 2; ++m) { + differentiable_plasticity->advanceTimestep(dt); + pv_writer.write(m + 1, differentiable_plasticity->time(), + differentiable_plasticity->getFieldStatesAndParamStates()); + } + + auto shape_disp = plastic_mechanics_system.field_store->getShapeDisp(); + auto params = plastic_mechanics_system.getParameterFields(); + + auto displacement = plastic_mechanics_system.field_store->getField(plastic_mechanics_system.prefix("displacement")); + auto disp_squared = 0.5 * innerProduct(displacement, displacement); + gretl::set_as_objective(disp_squared); + + EXPECT_GT(checkGradWrt(disp_squared, shape_disp, 2e-4, 4, true), 0.95); + EXPECT_GT(checkGradWrt(disp_squared, params[0], 100.0, 4, true), 0.95); + EXPECT_GT(checkGradWrt(disp_squared, params[1], 100.0, 4, true), 0.95); + + std::vector sensitivity_fields{shape_disp, params[0], params[1]}; + auto sens_writer = + createParaviewWriter(*mesh, sensitivity_fields, physics_name + "_final_sensitivities"); + + auto& graph = disp_squared.data_store(); + graph.reset(); + gretl::set_as_objective(disp_squared); + graph.back_prop(); + + sens_writer.write(0, differentiable_plasticity->time(), sensitivity_fields); } -} // namespace smith +} // namespace smith int main(int argc, char* argv[]) { From 625867a854843f9ce82f2882f1a2453da8831068 Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 14 Jul 2026 15:46:28 -0700 Subject: [PATCH 05/28] fix style --- .../tests/test_differentiable_plasticity.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index b53e106391..4d54725f52 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -338,8 +338,7 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) EXPECT_GT(checkGradWrt(disp_squared, params[1], 100.0, 4, true), 0.95); std::vector sensitivity_fields{shape_disp, params[0], params[1]}; - auto sens_writer = - createParaviewWriter(*mesh, sensitivity_fields, physics_name + "_final_sensitivities"); + auto sens_writer = createParaviewWriter(*mesh, sensitivity_fields, physics_name + "_final_sensitivities"); auto& graph = disp_squared.data_store(); graph.reset(); From 62cf27462350e4e7c50f39c2ca6797a7f8be1530 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 15 Jul 2026 12:51:20 -0700 Subject: [PATCH 06/28] fix doc --- .github/workflows/ci-tests.yml | 2 +- .../plastic_mechanics_system.hpp | 26 ++++++++++++------- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/.github/workflows/ci-tests.yml b/.github/workflows/ci-tests.yml index a5b1d5fb0d..4f98dd03ed 100644 --- a/.github/workflows/ci-tests.yml +++ b/.github/workflows/ci-tests.yml @@ -79,7 +79,7 @@ jobs: echo "compiler_image ${{ matrix.config.compiler_image }}" echo "host_config ${{ matrix.config.host_config }}" - name: Build and Test ${{ matrix.config.host_config }} ${{ matrix.build_type }} - timeout-minutes: 60 + timeout-minutes: 80 run: | ./scripts/llnl/build_src.py ${{ matrix.config.build_src_opts }} \ --verbose \ diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 74ae9c4179..b382ef1300 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -5,8 +5,8 @@ // SPDX-License-Identifier: (BSD-3-Clause) /** - * @file electro_mechanics_system.hpp - * @brief Defines the ElectroMechanicsystem struct and its factory function + * @file plastic_mechanics_system.hpp + * @brief Defines the PlasticMechanicsystem struct and its factory function */ #pragma once @@ -30,7 +30,7 @@ namespace smith { * * This system rely on a fixed-point iteration to solve plasticity equivalent to the traditional * inner-outer Newton method. - * n \in [0, N] denotes the timestep and k \in [0, L] denotes the fixed-point iteration count. + * $$n \in [0, N]$$ denotes the timestep and $$k \in [0, L]$$ denotes the fixed-point iteration count. * At every new timestep n+1, first predict new displacement by [u^{n+1,k}, Fp^{n+1,k}, epsilon_p^{n+1,k}] ---> * u^{n+1,k+1}. Then update internal state variables [u^{n+1,k+1}, Fp^n, epsilon_p^n] ---> [Fp^{n+1,k+1}, * epsilon_p^{n+1,k+1}]. @@ -68,13 +68,16 @@ struct PlasticMechanicsSystem : public SystemBase { std::shared_ptr disp_bc; ///< Displacement boundary conditions. // Internal variable bcs - std::shared_ptr plastic_deform_bc; - std::shared_ptr plastic_strain_bc; + std::shared_ptr plastic_deform_bc; ///< Plastic deformation gradient conditions. + std::shared_ptr plastic_strain_bc; ///< Plastic strain conditions - std::shared_ptr qusistatic_time_rule; - std::shared_ptr backward_euler_time_rule; + std::shared_ptr qusistatic_time_rule; ///< Quasistatic time integration rule + std::shared_ptr + backward_euler_time_rule; ///< Backward euler time integration rule - // Helper functions for plastic deformation gradient bookkeeping + /** + * @brief Transform tensor variable stored as [dim * dim, 1] to [dim, dim] + */ template SMITH_HOST_DEVICE tensor recoverTensor(const tensor F_state) { @@ -89,6 +92,9 @@ struct PlasticMechanicsSystem : public SystemBase { return F_state_tensor; } + /** + * @brief Transform tensor variable stored as [dim, dim] to [dim * dim, 1] + */ template SMITH_HOST_DEVICE tensor flattenTensor(const tensor F_state) { @@ -155,8 +161,8 @@ struct PlasticMechanicsSystem : public SystemBase { /** * @brief Set the material model for a domain, defining integrals for weak form. * @tparam MaterialType The material model type. - * @param material The material model instance. * @param domain_name The name of the domain to apply the material to. + * @param mat The material model instance. */ template void setMaterial(const std::string& domain_name, const MaterialType& mat) @@ -176,8 +182,8 @@ struct PlasticMechanicsSystem : public SystemBase { /** * @brief Set the plasticity update model for a domain, defining integrals for weak form. * @tparam MaterialType The material model type. - * @param material The material model instance. * @param domain_name The name of the domain to apply the material to. + * @param mat The material model instance. */ template void setPlasticity(const std::string& domain_name, const MaterialType& mat) From bcaa6a654e24db194e5d061bb40427c9e33fac19 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 15 Jul 2026 14:57:59 -0700 Subject: [PATCH 07/28] fix docs --- src/smith/differentiable_numerics/plastic_mechanics_system.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index b382ef1300..0509ac31fe 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -30,7 +30,7 @@ namespace smith { * * This system rely on a fixed-point iteration to solve plasticity equivalent to the traditional * inner-outer Newton method. - * $$n \in [0, N]$$ denotes the timestep and $$k \in [0, L]$$ denotes the fixed-point iteration count. + * n in [0, N] denotes the timestep and k in [0, L] denotes the fixed-point iteration count. * At every new timestep n+1, first predict new displacement by [u^{n+1,k}, Fp^{n+1,k}, epsilon_p^{n+1,k}] ---> * u^{n+1,k+1}. Then update internal state variables [u^{n+1,k+1}, Fp^n, epsilon_p^n] ---> [Fp^{n+1,k+1}, * epsilon_p^{n+1,k+1}]. From ee497588b78d93de8a623760d7f69f292fe7a7e1 Mon Sep 17 00:00:00 2001 From: li97 Date: Thu, 20 Aug 2026 16:21:19 -0700 Subject: [PATCH 08/28] update API --- .../plastic_mechanics_system.hpp | 103 +++++++----------- .../differentiable_numerics/system_solver.cpp | 38 +++++++ .../tests/test_differentiable_plasticity.cpp | 45 ++++---- 3 files changed, 104 insertions(+), 82 deletions(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 0509ac31fe..658ddc8dff 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -16,9 +16,8 @@ #include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp" #include "smith/differentiable_numerics/multiphysics_time_integrator.hpp" #include "smith/differentiable_numerics/time_integration_rule.hpp" -#include "smith/differentiable_numerics/time_discretized_weak_form.hpp" #include "smith/differentiable_numerics/differentiable_physics.hpp" -#include "smith/physics/weak_form.hpp" +#include "smith/physics/functional_weak_form.hpp" #include "smith/differentiable_numerics/system_base.hpp" namespace smith { @@ -37,25 +36,26 @@ namespace smith { */ template struct PlasticMechanicsSystem : public SystemBase { + using SystemBase::SystemBase; // Primary weak forms to solve for displacement /// @brief using for SolidWeakFormType with inputs [u^{n+1, k+1}, Fp^{n+1, k}, epsilon_p^{n+1, k}] - using SolidWeakFormType = TimeDiscretizedWeakForm< + using SolidWeakFormType = FunctionalWeakForm< dim, H1, Parameters, L2, L2, parameter_space...>>; /// @brief using for PlasticDeformWeakFormType with inputs [Fp^{n+1, k+1}, Fp^n, epsilon_p^{n+1, k+1}, epsilon_p^n, /// u^{n+1, k+1}] using PlasticDeformWeakFormType = - TimeDiscretizedWeakForm, - Parameters, L2, L2, - L2, H1, parameter_space...>>; + FunctionalWeakForm, + Parameters, L2, L2, + L2, H1, parameter_space...>>; /// @brief using for PlasticStrainWeakFormType with inputs [epsilon_p^{n+1, k+1}, epsilon_p^n, Fp^{n+1, k+1}, Fp^n, /// u^{n+1, k+1}] using PlasticStrainWeakFormType = - TimeDiscretizedWeakForm, - Parameters, L2, L2, - L2, H1, parameter_space...>>; + FunctionalWeakForm, + Parameters, L2, L2, + L2, H1, parameter_space...>>; // Primary weak forms std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. @@ -71,7 +71,7 @@ struct PlasticMechanicsSystem : public SystemBase { std::shared_ptr plastic_deform_bc; ///< Plastic deformation gradient conditions. std::shared_ptr plastic_strain_bc; ///< Plastic strain conditions - std::shared_ptr qusistatic_time_rule; ///< Quasistatic time integration rule + std::shared_ptr quasistatic_time_rule; ///< Quasistatic time integration rule std::shared_ptr backward_euler_time_rule; ///< Backward euler time integration rule @@ -109,32 +109,15 @@ struct PlasticMechanicsSystem : public SystemBase { return F_state_flattened; } - /** - * @brief Get the list of all state fields. - * @return std::vector List of state fields. - */ - std::vector getStateFields() const - { - std::vector states; - states.push_back(field_store->getField(prefix("displacement"))); - - // Internal states - states.push_back(field_store->getField(prefix("plastic_defgrad"))); - states.push_back(field_store->getField(prefix("plastic_defgrad_old"))); - states.push_back(field_store->getField(prefix("plastic_strain"))); - states.push_back(field_store->getField(prefix("plastic_strain_old"))); - - return states; - } - /** * @brief Get the list of physical state fields for visualization. * @return std::vector List of physical fields suitable for output. */ std::vector getOutputFieldStates() const { - return {field_store->getField(prefix("displacement")), field_store->getField(prefix("plastic_defgrad")), - field_store->getField(prefix("plastic_strain"))}; + return {field_store->getField(field_store->prefix("displacement")), + field_store->getField(field_store->prefix("plastic_defgrad")), + field_store->getField(field_store->prefix("plastic_strain"))}; } /** @@ -143,19 +126,7 @@ struct PlasticMechanicsSystem : public SystemBase { */ std::vector getReactionInfos() const { - return {{prefix("solid_force"), &field_store->getField(prefix("displacement")).get()->space()}}; - } - - /** - * @brief Create a DifferentiablePhysics object for this system. - * @param physics_name The name of the physics. - * @return std::shared_ptr The differentiable physics object. - */ - std::shared_ptr createDifferentiablePhysics(std::string physics_name) - { - return std::make_shared(field_store->getMesh(), field_store->graph(), - field_store->getShapeDisp(), getStateFields(), getParameterFields(), - advancer, physics_name, getReactionInfos()); + return {{field_store->prefix("solid_force"), &field_store->getField(field_store->prefix("displacement")).get()->space()}}; } /** @@ -223,8 +194,8 @@ struct PlasticMechanicsSystem : public SystemBase { * @brief Factory function to build a chemomechanics system with L2 state variables. */ template -PlasticMechanicsSystem buildPlasticMechanicsSystem( - std::shared_ptr mesh, std::shared_ptr solver, std::string prepend_name = "", +std::shared_ptr> buildPlasticMechanicsSystem( + std::shared_ptr mesh, std::shared_ptr solver, std::string prepend_name = "", FieldType... parameter_types) { auto field_store = std::make_shared(mesh, 100); @@ -258,9 +229,14 @@ PlasticMechanicsSystem buildPlasticMechanic field_store->addDependent(plastic_strain_type, FieldStore::TimeDerivative::VAL, prefix("plastic_strain_old")); // Parameters - std::vector parameter_fields; - (field_store->addParameter(FieldType(prefix("param_" + parameter_types.name))), ...); - (parameter_fields.push_back(field_store->getField(prefix("param_" + parameter_types.name))), ...); + auto register_parameter = [&](auto& parameter_type) { + parameter_type.name = prefix("param_" + parameter_type.name); + field_store->addParameter(parameter_type); + }; + + if constexpr (sizeof...(parameter_space) > 0) { + (register_parameter(parameter_types), ...); + } using SystemType = PlasticMechanicsSystem; @@ -269,35 +245,38 @@ PlasticMechanicsSystem buildPlasticMechanic auto solid_weak_form = std::make_shared( solid_res_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), field_store->createSpaces(solid_res_name, disp_type.name, disp_type, plastic_defgrad_type, plastic_strain_type, - FieldType(prefix("param_" + parameter_types.name))...)); + parameter_types...)); std::string plastic_defgrad_res_name = prefix("plastic_defgrad_residual"); auto plastic_defgrad_weak_form = std::make_shared( plastic_defgrad_res_name, field_store->getMesh(), field_store->getField(plastic_defgrad_type.name).get()->space(), field_store->createSpaces(plastic_defgrad_res_name, plastic_defgrad_type.name, plastic_defgrad_type, plastic_defgrad_old_type, plastic_strain_type, plastic_strain_old_type, disp_type, - FieldType(prefix("param_" + parameter_types.name))...)); + parameter_types...)); std::string plastic_strain_res_name = prefix("plastic_strain_residual"); auto plastic_strain_weak_form = std::make_shared( plastic_strain_res_name, field_store->getMesh(), field_store->getField(plastic_strain_type.name).get()->space(), field_store->createSpaces(plastic_strain_res_name, plastic_strain_type.name, plastic_strain_type, plastic_strain_old_type, plastic_defgrad_type, plastic_defgrad_old_type, disp_type, - FieldType(prefix("param_" + parameter_types.name))...)); + parameter_types...)); // Build solver and advancer std::vector> weak_forms{solid_weak_form, plastic_defgrad_weak_form, plastic_strain_weak_form}; - auto advancer = std::make_shared(field_store, weak_forms, solver); - - return SystemType{{field_store, solver, advancer, parameter_fields, prepend_name}, - solid_weak_form, - plastic_defgrad_weak_form, - plastic_strain_weak_form, - disp_bc, - plastic_defgrad_bc, - plastic_strain_bc, - quasistatic_time_rule, - backward_euler_time_rule}; + + auto sys = + std::make_shared(field_store, solver, weak_forms); + + sys->solid_weak_form = solid_weak_form; + sys->plastic_deform_weak_form = plastic_defgrad_weak_form; + sys->plastic_strain_weak_form = plastic_strain_weak_form; + sys->disp_bc = disp_bc; + sys->plastic_deform_bc = plastic_defgrad_bc; + sys->plastic_strain_bc = plastic_strain_bc; + sys->quasistatic_time_rule = quasistatic_time_rule; + sys->backward_euler_time_rule = backward_euler_time_rule; + + return sys; } } // namespace smith diff --git a/src/smith/differentiable_numerics/system_solver.cpp b/src/smith/differentiable_numerics/system_solver.cpp index ca017c511a..338bcfedf6 100644 --- a/src/smith/differentiable_numerics/system_solver.cpp +++ b/src/smith/differentiable_numerics/system_solver.cpp @@ -75,6 +75,8 @@ std::vector SystemSolver::solve(const std::vector& residu SLIC_ERROR_IF(stages_.empty(), "SystemSolver has no stages defined."); size_t num_residuals = residual_evals.size(); + // const bool trace_plasticity = + // num_residuals == 3 && residual_evals[1]->name().find("plastic_defgrad_residual") != std::string::npos; std::vector active_stages = stages_; for (auto& stage : active_stages) { if (stage.block_indices.empty()) { @@ -108,6 +110,14 @@ std::vector SystemSolver::solve(const std::vector& residu // Working copy of states, updated in-place as stages solve std::vector> current_states = states; + // auto trace_row = [&](const char* location, size_t row) { + // for (size_t slot = 0; slot < current_states[row].size(); ++slot) { + // const auto& state = current_states[row][slot]; + // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] {} row={} slot={} name={} norm={:.16e}", location, row, + // slot, state.get()->name(), state.get()->Norml2())); + // } + // }; + // Pre-compute name -> (row, slot) routing so the propagation loop avoids O(N*M) string compares // on every staggered iteration. Field-name identity within current_states is invariant across // the iteration loop: only values are replaced, never the underlying name. @@ -182,6 +192,16 @@ std::vector SystemSolver::solve(const std::vector& residu stage_block_indices.push_back(row_indices); } + // // The plastic mechanics system orders its blocks as displacement, plastic deformation + // // gradient, and plastic strain. Trace the deformation-gradient inputs after the strain + // // stage should have propagated its updated state. + // if (trace_plasticity && stage.block_indices.size() == 1 && stage.block_indices[0] == 1) { + // trace_row("before Fp solve", 1); + // auto residual = eval_residual_and_zero_bcs(1); + // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] Fp residual norm before solve={:.16e}", + // residual.Norml2())); + // } + std::vector stage_solutions = block_solve(stage_residuals, stage_block_indices, shape_disp, stage_states, stage_params, time_info, stage.solver.get(), stage_bc_managers); @@ -200,11 +220,29 @@ std::vector SystemSolver::solve(const std::vector& residu } auto it = field_routing.find(new_state.get()->name()); + // if (trace_plasticity && global_col == 2) { + // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] solved strain name={} norm={:.16e}", + // new_state.get()->name(), new_state.get()->Norml2())); + // if (it == field_routing.end()) { + // SLIC_INFO_ROOT("[plastic-debug] no routing entry for solved strain"); + // } else { + // for (const auto& [row, slot] : it->second) { + // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] strain routes to row={} slot={}", row, slot)); + // } + // } + // } if (it != field_routing.end()) { for (const auto& [r, slot] : it->second) { current_states[r][slot] = new_state; } } + // if (trace_plasticity && global_col == 2) { + // trace_row("after strain routing", 1); + // } + // if (trace_plasticity && global_col == 1) { + // SLIC_INFO_ROOT( + // axom::fmt::format("[plastic-debug] solved Fp norm={:.16e}", new_state.get()->Norml2())); + // } } } diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index 4d54725f52..057968f19f 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -128,7 +128,7 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) const std::string meshtag = "mesh"; auto mesh = std::make_shared(buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); - auto staggered_coupled_solver = std::make_shared(50); + auto staggered_coupled_solver = std::make_shared(50); auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); auto strain_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); auto defgrad_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); @@ -143,12 +143,12 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) Hardening hardening{.sigma_y = 40.0, .Hi = 50.0, .eta = 0.0}; DifferentiableJ2SmallStrain mat(hardening, 1e+4, 0.25, 5.0, 1.0); - plastic_mechanics_system.setMaterial(mesh->entireBodyName(), mat); - plastic_mechanics_system.setPlasticity(mesh->entireBodyName(), mat); + plastic_mechanics_system->setMaterial(mesh->entireBodyName(), mat); + plastic_mechanics_system->setPlasticity(mesh->entireBodyName(), mat); // prescribe zero displacement at the supported end of the beam, mesh->addDomainOfBoundaryElements("support", by_attr(1)); - plastic_mechanics_system.disp_bc->setFixedVectorBCs(mesh->domain("support")); + plastic_mechanics_system->disp_bc->setFixedVectorBCs(mesh->domain("support")); // apply a displacement along z to the the tip of the beam mesh->addDomainOfBoundaryElements("tip", by_attr(2)); @@ -157,13 +157,15 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) u[2] = t * (t - 1); return u; }; - plastic_mechanics_system.disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); + plastic_mechanics_system->disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); double dt = 0.1; double time = 0.0; - auto shape_disp = plastic_mechanics_system.field_store->getShapeDisp(); - auto states = plastic_mechanics_system.getStateFields(); - auto params = plastic_mechanics_system.getParameterFields(); + auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); + auto states = plastic_mechanics_system->field_store->getStateFields(); + auto params = plastic_mechanics_system->field_store->getParameterFields(); + + auto advancer = makeAdvancer(plastic_mechanics_system); auto pv_writer = createParaviewWriter(*mesh, states, "J2_small_strain"); pv_writer.write(0, 0.0, states); @@ -171,12 +173,13 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) std::vector reactions; for (size_t step = 0; step < 10; ++step) { std::tie(states, reactions) = - plastic_mechanics_system.advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); + advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); time += dt; pv_writer.write(step + 1, time, states); } - auto displacement = plastic_mechanics_system.field_store->getField(plastic_mechanics_system.prefix("displacement")); + auto displacement = plastic_mechanics_system->field_store->getField( + plastic_mechanics_system->field_store->prefix("displacement")); double final_disp_l2 = norm(*displacement.get(), 2.0); EXPECT_GT(final_disp_l2, 0.0); } @@ -280,7 +283,7 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) const std::string meshtag = "mesh"; auto mesh = std::make_shared(buildMeshFromFile(filename), meshtag, serial_refinement, parallel_refinement); - auto staggered_coupled_solver = std::make_shared(100); + auto staggered_coupled_solver = std::make_shared(100); auto primal_block_solver = buildNonlinearBlockSolver(primal_nonlin_opts, primal_lin_opts, *mesh); auto strain_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); auto defgrad_block_solver = buildNonlinearBlockSolver(state_nonlin_opts, state_lin_opts, *mesh); @@ -296,12 +299,12 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) Hardening hardening{.sigma_y = 25.0, .Hi = 50.0, .eta = 0.0}; ParameterizedDifferentiableJ2SmallStrain mat(hardening, 1.5e+4, 0.25, 500.0, 1.0); - plastic_mechanics_system.setMaterial(mesh->entireBodyName(), mat); - plastic_mechanics_system.setPlasticity(mesh->entireBodyName(), mat); + plastic_mechanics_system->setMaterial(mesh->entireBodyName(), mat); + plastic_mechanics_system->setPlasticity(mesh->entireBodyName(), mat); // prescribe zero displacement at the supported end of the beam, mesh->addDomainOfBoundaryElements("support", by_attr(1)); - plastic_mechanics_system.disp_bc->setFixedVectorBCs(mesh->domain("support")); + plastic_mechanics_system->disp_bc->setFixedVectorBCs(mesh->domain("support")); // apply a displacement along z to the the tip of the beam mesh->addDomainOfBoundaryElements("tip", by_attr(2)); @@ -310,10 +313,11 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) u[2] = t; return u; }; - plastic_mechanics_system.disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); + plastic_mechanics_system->disp_bc->setVectorBCs(mesh->domain("tip"), {2}, translated_in_z); - std::shared_ptr differentiable_plasticity = - plastic_mechanics_system.createDifferentiablePhysics(physics_name); + auto advancer = makeAdvancer(plastic_mechanics_system); + std::unique_ptr differentiable_plasticity = + makeDifferentiablePhysics(plastic_mechanics_system, advancer, physics_name); auto pv_writer = createParaviewWriter(*mesh, differentiable_plasticity->getFieldStatesAndParamStates(), physics_name, smith::ParaviewWriter::Options{.write_duals = false}); @@ -326,10 +330,11 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) differentiable_plasticity->getFieldStatesAndParamStates()); } - auto shape_disp = plastic_mechanics_system.field_store->getShapeDisp(); - auto params = plastic_mechanics_system.getParameterFields(); + auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); + auto params = plastic_mechanics_system->field_store->getParameterFields(); - auto displacement = plastic_mechanics_system.field_store->getField(plastic_mechanics_system.prefix("displacement")); + auto displacement = plastic_mechanics_system->field_store->getField( + plastic_mechanics_system->field_store->prefix("displacement")); auto disp_squared = 0.5 * innerProduct(displacement, displacement); gretl::set_as_objective(disp_squared); From 78b2ff8fc0d9450edcc62410de1f49242899d473 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 24 Aug 2026 15:46:45 -0700 Subject: [PATCH 09/28] Use index logic instead of strings to establish the location of unknowns in weak forms, and how they are updated for staggered solves. --- .../differentiable_numerics/field_store.cpp | 98 ++++++------- .../differentiable_numerics/field_store.hpp | 104 +++++++------- .../multiphysics_time_integrator.cpp | 20 +-- .../nonlinear_solve.cpp | 67 +++++---- .../nonlinear_solve.hpp | 22 +-- .../plastic_mechanics_system.hpp | 55 ++++---- .../solid_mechanics_system.hpp | 4 +- .../state_variable_system.hpp | 2 +- .../differentiable_numerics/system_base.cpp | 49 ++++--- .../differentiable_numerics/system_base.hpp | 10 +- .../differentiable_numerics/system_solver.cpp | 87 ++++-------- .../differentiable_numerics/system_solver.hpp | 8 +- .../tests/test_combined_thermo_mechanics.cpp | 2 +- .../tests/test_differentiable_plasticity.cpp | 25 +++- .../tests/test_mixed_poisson.cpp | 2 +- .../test_multiphysics_time_integrator.cpp | 133 ++++++++++++++++-- .../tests/test_porous_heat_sink.cpp | 2 +- .../thermal_system.hpp | 2 +- src/smith/physics/functional_weak_form.hpp | 3 + .../tests/test_functional_weak_form.cpp | 26 ++++ 20 files changed, 418 insertions(+), 303 deletions(-) diff --git a/src/smith/differentiable_numerics/field_store.cpp b/src/smith/differentiable_numerics/field_store.cpp index 20d0869ade..80cda52ae0 100644 --- a/src/smith/differentiable_numerics/field_store.cpp +++ b/src/smith/differentiable_numerics/field_store.cpp @@ -9,6 +9,8 @@ #include "smith/differentiable_numerics/nonlinear_solve.hpp" #include "gretl/wang_checkpoint_strategy.hpp" +#include + namespace smith { FieldStore::FieldStore(std::shared_ptr mesh, size_t storage_size, std::string prepend_name) @@ -64,62 +66,43 @@ void FieldStore::printMap() } } -std::vector> FieldStore::indexMap(const std::vector& residual_names) const +BlockIndexMap FieldStore::indexMap(const std::vector& residual_names) const { - // Build a local column space: each residual in the subsystem contributes one local column, - // corresponding to its "self" diagonal unknown. The self-unknown is preferably the residual's - // reaction (test) field if that field appears in the unknown-arg list for this weak form; - // otherwise fall back on the first unknown argument (handles cases like the cycle-zero - // acceleration solve, where the reaction field is a dependent/history field). std::map global_state_to_local_col; for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { const std::string& res_name = residual_names[res_i]; - size_t global_state_idx = invalid_block_index; - - std::string reaction_name; - for (const auto& kv : weak_form_to_test_field_) { - if (kv.first == res_name) { - reaction_name = kv.second; - break; - } - } - - // Check if the reaction field is one of the registered unknown args for this weak form. - bool reaction_is_unknown = false; - if (!reaction_name.empty() && weak_form_name_to_unknown_name_index_.count(res_name)) { - for (const auto& label : weak_form_name_to_unknown_name_index_.at(res_name)) { - if (label.field_name == reaction_name) { - reaction_is_unknown = true; - break; - } - } - } - - if (reaction_is_unknown) { - global_state_idx = to_states_index_.at(reaction_name); - } else { - const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name); - SLIC_ERROR_IF(arg_info.empty(), - "Weak form '" << res_name << "' has no unknown arguments; cannot build index map."); - global_state_idx = to_states_index_.at(arg_info.front().field_name); - } - global_state_to_local_col[global_state_idx] = res_i; + size_t global_state_idx = to_states_index_.at(getWeakFormFieldNames(res_name).unknown); + auto [owner, inserted] = global_state_to_local_col.emplace(global_state_idx, res_i); + SLIC_ERROR_IF(!inserted, axom::fmt::format("Weak forms '{}' and '{}' both own state index {}", + residual_names[owner->second], res_name, global_state_idx)); } - std::vector> block_indices(residual_names.size()); + BlockIndexMap block_indices(residual_names.size(), + std::vector(residual_names.size(), BlockArgumentIndices{})); 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(residual_names.size(), invalid_block_index); const std::string& res_name = residual_names[res_i]; - const auto& arg_info = weak_form_name_to_unknown_name_index_.at(res_name); + const auto& field_names = weak_form_name_to_field_names_.at(res_name); + const auto& dependencies = weak_form_name_to_unknown_name_index_.at(res_name); + size_t state_slot = 0; + for (size_t argument_index = 0; argument_index < field_names.size(); ++argument_index) { + const auto& field_name = field_names[argument_index]; + if (!to_states_index_.count(field_name)) { + continue; + } - for (const auto& field_name_and_arg_index : arg_info) { - size_t global_state_index = to_states_index_.at(field_name_and_arg_index.field_name); - auto it = global_state_to_local_col.find(global_state_index); - if (it != global_state_to_local_col.end()) { - res_indices[it->second] = field_name_and_arg_index.field_index_in_residual; + size_t global_state_index = to_states_index_.at(field_name); + auto column = global_state_to_local_col.find(global_state_index); + if (column != global_state_to_local_col.end()) { + bool marked_dependency = std::any_of(dependencies.begin(), dependencies.end(), [&](const auto& dependency) { + return dependency.field_index_in_residual == argument_index && dependency.field_name == field_name; + }); + SLIC_ERROR_IF(!marked_dependency, + axom::fmt::format("Field '{}' in weak form '{}' belongs to solved column {} but is not marked " + "as an unknown dependency", + field_name, res_name, column->second)); + block_indices[res_i][column->second].push_back(state_slot); } - // else: field belongs to a different subsystem; treat as fixed input here. + ++state_slot; } } @@ -132,7 +115,7 @@ std::vector FieldStore::getBoundaryConditionMan std::vector field_names; field_names.reserve(weak_form_names.size()); for (const auto& wf_name : weak_form_names) { - field_names.push_back(getWeakFormReaction(wf_name)); + field_names.push_back(getWeakFormFieldNames(wf_name).unknown); } return getBoundaryConditionManagersForFields(field_names); } @@ -359,26 +342,27 @@ void FieldStore::markWeakFormInternal(const std::string& weak_form_name) internal_weak_forms_.insert(weak_form_name); } -void FieldStore::addWeakFormReaction(std::string weak_form_name, std::string field_name) +void FieldStore::addWeakFormFieldNames(std::string weak_form_name, UnknownAndTestFieldNames field_names) { - for (auto& kv : weak_form_to_test_field_) { + for (auto& kv : weak_form_field_names_) { if (kv.first == weak_form_name) { - kv.second = field_name; + kv.second = std::move(field_names); return; } } - weak_form_to_test_field_.push_back({weak_form_name, field_name}); + weak_form_field_names_.push_back({std::move(weak_form_name), std::move(field_names)}); } -std::string FieldStore::getWeakFormReaction(const std::string& weak_form_name) const +const UnknownAndTestFieldNames& FieldStore::getWeakFormFieldNames(const std::string& weak_form_name) const { - for (const auto& kv : weak_form_to_test_field_) { + for (const auto& kv : weak_form_field_names_) { if (kv.first == weak_form_name) { return kv.second; } } - SLIC_ERROR("Reaction field not found for weak form " << weak_form_name); - return ""; + SLIC_ERROR("Field names not found for weak form " << weak_form_name); + static const UnknownAndTestFieldNames empty; + return empty; } const std::vector& FieldStore::getParameterFields() const { return params_; } @@ -405,12 +389,12 @@ std::vector FieldStore::getOutputFieldStates() const std::vector FieldStore::getReactionInfos() const { std::vector infos; - for (const auto& kv : weak_form_to_test_field_) { + for (const auto& kv : weak_form_field_names_) { const std::string& weak_form_name = kv.first; if (internal_weak_forms_.count(weak_form_name)) { continue; } - const std::string& field_name = kv.second; + const std::string& field_name = kv.second.test; infos.push_back({weak_form_name, &getField(field_name).get()->space()}); } return infos; diff --git a/src/smith/differentiable_numerics/field_store.hpp b/src/smith/differentiable_numerics/field_store.hpp index d85b3a6338..939aef4205 100644 --- a/src/smith/differentiable_numerics/field_store.hpp +++ b/src/smith/differentiable_numerics/field_store.hpp @@ -7,6 +7,7 @@ #pragma once #include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" #include "smith/differentiable_numerics/time_integration_rule.hpp" #include "smith/physics/functional_weak_form.hpp" #include "smith/physics/mesh.hpp" @@ -30,15 +31,18 @@ struct ReactionInfo { const mfem::ParFiniteElementSpace* space = nullptr; ///< The finite element space of the dual field. }; +/// Names identifying a weak form's test field and solver-owned unknown. +struct UnknownAndTestFieldNames { + std::string unknown; ///< Field owned and solved by this residual row. + std::string test; ///< Field defining residual test space and reaction output. +}; + /** - * @brief Representation of a field type with a name and a flag indicating whether it is an - * active Jacobian unknown in the current weak-form context. + * @brief Representation of a field type with a name and dependency flag. * - * @c is_unknown is intentionally a per-instance flag, not a global property. The same field - * may be a Jacobian variable in one weak form (e.g. displacement in the main solid solve) and - * a fixed input in another (e.g. displacement in the stress projection). Callers control this - * by passing the same @c FieldType object (set by @c addIndependent) or a plain copy with the - * default @c is_unknown = false. + * @c is_unknown is intentionally a per-instance flag. It marks every occurrence whose value + * depends on a field solved within the same system. Solver stages determine which dependencies + * are active during each block solve. History and external fields remain fixed. * * @tparam Space The finite element space type. * @tparam Time The time integration type (unused by default). @@ -50,11 +54,11 @@ struct FieldType { /** * @brief Construct a new FieldType object. * @param n Name of the field. - * @param is_unknown_ Whether this field is a Jacobian unknown in the current context. + * @param is_unknown_ Whether this field depends on a system-solved unknown. */ FieldType(std::string n, bool is_unknown_ = false) : name(n), is_unknown(is_unknown_) {} std::string name; ///< Name of the field. - bool is_unknown = false; ///< True if this field is a Jacobian variable in the current weak-form context. + bool is_unknown = false; ///< True when this argument depends on a system-solved unknown. }; /** @@ -129,7 +133,7 @@ struct FieldStore { * * Registers the field as an unknown by setting @c type.is_unknown = true, so the same * @c FieldType object can later be passed to @c createSpaces to mark this argument - * as an active Jacobian variable. Also creates a boundary-condition slot keyed by field + * as a solved dependency. Also creates a boundary-condition slot keyed by field * name that callers can populate after this call returns. * * @tparam Space The finite element space type. @@ -173,9 +177,9 @@ struct FieldStore { * **Return value:** a @c FieldType whose @c name is the name of the newly registered * field and @c is_unknown is @c false. This object is intentionally returned (rather than * discarded) so that callers can pass it directly to @c createSpaces when assembling the - * weak-form argument list. To make it the Jacobian variable for a specific weak form (e.g. + * weak-form argument list. To make it a solved dependency for a specific weak form (e.g. * the cycle-zero acceleration solve), copy the returned object and set @c is_unknown = true - * before passing to @c createSpaces. + * before passing it to @c createSpaces. * * @tparam Space The finite element space type (must match the independent field). * @param independent_field The @c FieldType of the independent (predicted) field. @@ -236,17 +240,8 @@ struct FieldStore { */ void addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index); - /** - * @brief Register the reaction (test) field for a weak form. - * - * The reaction field is the field whose test function space the weak form integrates against. - * It determines which field's degrees of freedom the assembled residual is "returned to" - * (i.e. the field whose force/flux vector is populated). - * - * @param weak_form_name Name of the weak form. - * @param field_name Name of the reaction field. - */ - void addWeakFormReaction(std::string weak_form_name, std::string field_name); + /// Register test-space and solver-owned field names for a weak form. + void addWeakFormFieldNames(std::string weak_form_name, UnknownAndTestFieldNames field_names); /** * @brief Mark a weak form as internal so it is excluded from getReactionInfos(). @@ -256,53 +251,49 @@ struct FieldStore { */ void markWeakFormInternal(const std::string& weak_form_name); - /** - * @brief Get the name of the reaction (test) field for a weak form. - * @param weak_form_name Name of the weak form. - * @return std::string Name of the reaction field. - */ - std::string getWeakFormReaction(const std::string& weak_form_name) const; + /// Get explicit test-space and solver-owned field names for a weak form. + const UnknownAndTestFieldNames& getWeakFormFieldNames(const std::string& weak_form_name) const; /** * @brief Register all input fields for a weak form and return their FE spaces. * * This is the primary setup method for constructing a weak form. It: - * 1. Registers @p reaction_field_name as the reaction/test field via @c addWeakFormReaction. + * 1. Registers explicit test-space and solver-owned field names. * 2. Iterates over every @c FieldType in @p types (in order), registering each as an input - * argument. If @c type.is_unknown is @c true, the field is also registered as an active - * Jacobian unknown for this weak form. + * argument. The solver-owned field and every argument with @c type.is_unknown are registered + * as dependencies on system-solved unknowns. * 3. Returns the ordered vector of finite element spaces. * - * A field may have @c is_unknown = true in one weak form and @c false in another (e.g. - * displacement is a Jacobian variable in the main solid solve but a fixed input in the stress - * projection). Callers control this by passing the @c FieldType returned from @c addIndependent - * (has @c is_unknown = true) or a plain copy constructed from the field name (has @c is_unknown - * = false). Similarly, a dependent field can be made the Jacobian variable for a specific weak - * form (e.g. acceleration in the cycle-zero solve) by copying the returned @c FieldType and - * setting @c is_unknown = true. + * Every occurrence of a field solved by another row in the same system must have + * @c is_unknown = true. Solver stages choose active Jacobian columns. Fixed history fields and + * fields owned by external systems remain false. * * @param weak_form_name Name of the weak form being constructed. - * @param reaction_field_name Name of the test/reaction field (may differ from the first input). + * @param field_names Names of the test-space and solver-owned fields. * @param types Ordered list of @c FieldType descriptors for every input argument. * @return std::vector Ordered input FE spaces. */ template std::vector createSpaces(const std::string& weak_form_name, - const std::string& reaction_field_name, + UnknownAndTestFieldNames field_names, FieldTypes... types) { - addWeakFormReaction(weak_form_name, reaction_field_name); + addWeakFormFieldNames(weak_form_name, field_names); std::vector spaces; size_t arg_num = 0; + bool unknown_found = false; auto register_field = [&](auto type) { spaces.push_back(&getField(type.name).get()->space()); addWeakFormArg(weak_form_name, type.name, arg_num); - if (type.is_unknown) { + if (type.is_unknown || type.name == field_names.unknown) { addWeakFormUnknownArg(weak_form_name, type.name, arg_num); } + unknown_found = unknown_found || type.name == field_names.unknown; ++arg_num; }; (register_field(types), ...); + SLIC_ERROR_IF(!unknown_found, "Unknown field '" << field_names.unknown << "' is not an argument of weak form '" + << weak_form_name << "'"); return spaces; } @@ -332,18 +323,18 @@ struct FieldStore { /** * @brief Generate an index map for the residuals. * @param residual_names Names of the residuals. - * @return std::vector> The index map. + * @return Slot-list matrix mapping residual rows and solved columns. */ - std::vector> indexMap(const std::vector& residual_names) const; + BlockIndexMap indexMap(const std::vector& residual_names) const; /** * @brief Get the boundary condition managers for the given weak forms, one per residual row. * - * For each weak form in @p weak_form_names the reaction (test) field name is looked up. The + * For each weak form in @p weak_form_names the solver-owned field name is looked up. The * returned manager is selected by consulting the registered @c TimeIntegrationMapping s: - * - reaction = primary or history slot -> value-level BC manager - * - reaction = second-derivative (ddot) slot -> second-derivative BC manager - * - reaction has its own DBC entry not tied to a mapping -> that DBC's value manager + * - unknown = primary or history slot -> value-level BC manager + * - unknown = second-derivative (ddot) slot -> second-derivative BC manager + * - unknown has its own DBC entry not tied to a mapping -> that DBC's value manager * - otherwise -> @c nullptr (solver skips null entries) * * The second-derivative manager is rebuilt on each call, so late value-BC additions are reflected. @@ -493,10 +484,9 @@ 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_names_; - std::vector> weak_form_to_test_field_; + std::vector> weak_form_field_names_; std::set internal_weak_forms_; ///< weak forms excluded from getReactionInfos() (subsystem-internal) std::vector, TimeIntegrationMapping>> time_integration_rules_; @@ -509,13 +499,15 @@ struct FieldStore { * Thin convenience wrapper: registers @p test_type as the reaction field, registers all * @p field_types as input arguments, and constructs the weak form in one call. */ -template -auto createWeakForm(std::string name, FieldType test_type, FieldStore& field_store, - FieldType... field_types) +template +auto createWeakForm(std::string name, FieldType test_type, FieldType unknown_type, + FieldStore& field_store, FieldType... field_types) { - return std::make_shared>>( + return std::make_shared< + FunctionalWeakForm>>( name, field_store.getMesh(), field_store.getField(test_type.name).get()->space(), - field_store.createSpaces(name, test_type.name, field_types...)); + field_store.createSpaces(name, {.unknown = unknown_type.name, .test = test_type.name}, unknown_type, + field_types...)); } /** diff --git a/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp index 7979726a51..9f6d1a8623 100644 --- a/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp @@ -25,8 +25,8 @@ MultiphysicsTimeIntegrator::MultiphysicsTimeIntegrator(std::shared_ptrweak_forms.size(); ++i) { const std::string wf_name = system_->weak_forms[i]->name(); - const std::string reaction_name = system_->field_store->getWeakFormReaction(wf_name); - main_unknown_name_to_local_idx_[reaction_name] = i; + const std::string unknown_name = system_->field_store->getWeakFormFieldNames(wf_name).unknown; + main_unknown_name_to_local_idx_[unknown_name] = i; } } @@ -69,14 +69,14 @@ std::pair, std::vector> MultiphysicsTimeI SLIC_ERROR_ROOT_IF(cycle_zero_unknowns.size() != cz_sys->weak_forms.size(), "Cycle zero system result count does not match number of cycle-zero weak forms"); - SLIC_ERROR_ROOT_IF(!cz_sys->solve_result_field_names.empty() && - cz_sys->solve_result_field_names.size() != cz_sys->weak_forms.size(), - "Cycle zero solve_result_field_names size does not match number of weak forms"); + SLIC_ERROR_ROOT_IF( + !cz_sys->solved_field_names.empty() && cz_sys->solved_field_names.size() != cz_sys->weak_forms.size(), + "Cycle zero solved_field_names size does not match number of weak forms"); for (size_t i = 0; i < cz_sys->weak_forms.size(); ++i) { const std::string result_field_name = - cz_sys->solve_result_field_names.empty() - ? system_->field_store->getWeakFormReaction(cz_sys->weak_forms[i]->name()) - : cz_sys->solve_result_field_names[i]; + cz_sys->solved_field_names.empty() + ? system_->field_store->getWeakFormFieldNames(cz_sys->weak_forms[i]->name()).unknown + : cz_sys->solved_field_names[i]; size_t result_field_state_idx = system_->field_store->getFieldIndex(result_field_name); current_states[result_field_state_idx] = cycle_zero_unknowns[i]; system_->field_store->setField(result_field_state_idx, cycle_zero_unknowns[i]); @@ -122,8 +122,8 @@ std::pair, std::vector> MultiphysicsTimeI for (const auto& ps : post_solve_systems_) { auto ps_unknowns = ps->solve(time_info); for (size_t i = 0; i < ps->weak_forms.size(); ++i) { - const std::string reaction_name = ps->field_store->getWeakFormReaction(ps->weak_forms[i]->name()); - size_t u_idx = ps->field_store->getFieldIndex(reaction_name); + const std::string unknown_name = ps->field_store->getWeakFormFieldNames(ps->weak_forms[i]->name()).unknown; + size_t u_idx = ps->field_store->getFieldIndex(unknown_name); ps->field_store->setField(u_idx, ps_unknowns[i]); } } diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index e600f6cec6..7929a8b36d 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -26,9 +26,8 @@ void applyBoundaryConditions(double time, const smith::BoundaryConditionManager* } } -std::vector block_solve(const std::vector& residual_evals, - const std::vector> block_indices, const FieldState& shape_disp, - const std::vector>& states, +std::vector block_solve(const std::vector& residual_evals, const BlockIndexMap block_indices, + const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, const std::vector& bc_managers) @@ -50,18 +49,20 @@ std::vector block_solve(const std::vector& residual_evals // Validate block_indices bounds against states sizes for (size_t row = 0; row < num_rows_; ++row) { + std::vector slot_owners(states[row].size(), invalid_block_index); for (size_t col = 0; col < num_rows_; ++col) { - size_t idx = block_indices[row][col]; - if (idx != invalid_block_index) { + for (size_t idx : block_indices[row][col]) { SLIC_ERROR_IF(idx >= states[row].size(), - axom::fmt::format("block_indices[{}][{}] = {} is out of bounds (states[{}].size() = {})", row, - col, idx, row, states[row].size())); + axom::fmt::format("block_indices[{}][{}] contains {} outside states[{}] size {}", row, col, idx, + row, states[row].size())); + SLIC_ERROR_IF( + slot_owners[idx] != invalid_block_index, + axom::fmt::format("State slot {} in row {} belongs to columns {} and {}", idx, row, slot_owners[idx], col)); + slot_owners[idx] = col; } } - // Validate that diagonal entry is not invalid - SLIC_ERROR_IF( - block_indices[row][row] == invalid_block_index, - axom::fmt::format("block_indices[{}][{}] (diagonal entry) must not be invalid_block_index", row, row)); + SLIC_ERROR_IF(block_indices[row][row].empty(), + axom::fmt::format("block_indices[{}][{}] diagonal entry must not be empty", row, row)); } std::vector num_state_inputs; @@ -115,9 +116,7 @@ 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) { - 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)"); + size_t prime_unknown_row_i = block_indices[row_i][row_i].front(); diagonal_fields[row_i] = std::make_shared(*input_fields[row_i][prime_unknown_row_i]); } @@ -128,9 +127,8 @@ std::vector block_solve(const std::vector& residual_evals 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]; - if (prime_unknown_ij != invalid_block_index) { - input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + for (size_t slot : block_indices[row_i][col_j]) { + input_fields[row_i][slot] = diagonal_fields[col_j]; } } } @@ -172,13 +170,17 @@ std::vector block_solve(const std::vector& residual_evals std::vector row_field_inputs = input_fields[row_i]; 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]; - if (field_index_to_diff != invalid_block_index) { - tangent_weights[field_index_to_diff] = 1.0; + const auto& field_indices_to_diff = block_indices[row_i][col_j]; + if (!field_indices_to_diff.empty()) { + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 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; + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 0.0; + } } else { jacobians[row_i].emplace_back(nullptr); } @@ -244,14 +246,15 @@ std::vector block_solve(const std::vector& residual_evals // make a copy so we don't accidentally override the original copy 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]]); + diagonal_fields[row_i] = + std::make_shared(*input_fields[row_i][block_indices[row_i][row_i].front()]); *diagonal_fields[row_i] = *s[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) { - if (block_indices[row_i][col_j] != invalid_block_index) { - input_fields[row_i][block_indices[row_i][col_j]] = diagonal_fields[col_j]; + for (size_t slot : block_indices[row_i][col_j]) { + input_fields[row_i][slot] = diagonal_fields[col_j]; } } } @@ -272,13 +275,17 @@ std::vector block_solve(const std::vector& residual_evals std::vector row_field_inputs = input_fields[row_i]; 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]; - if (field_index_to_diff != invalid_block_index) { - tangent_weights[field_index_to_diff] = 1.0; + const auto& field_indices_to_diff = block_indices[row_i][col_j]; + if (!field_indices_to_diff.empty()) { + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 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; + for (size_t slot : field_indices_to_diff) { + tangent_weights[slot] = 0.0; + } } else { jacobians[row_i].emplace_back(nullptr); } @@ -351,8 +358,8 @@ std::vector block_solve(const std::vector& residual_evals // No sensitivity needed for primal fields for (size_t row_i = 0; row_i < num_rows; ++row_i) { for (size_t col_j = 0; col_j < num_rows; ++col_j) { - if (block_indices[row_i][col_j] != invalid_block_index) { - field_sensitivities[row_i][block_indices[row_i][col_j]] = nullptr; + for (size_t slot : block_indices[row_i][col_j]) { + field_sensitivities[row_i][slot] = nullptr; } } } diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp index c083bd018b..c3d79fb0a4 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.hpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -28,15 +28,21 @@ class DirichletBoundaryConditions; /// @brief magic number for representing a field which is not an argument of the weak form. static constexpr size_t invalid_block_index = std::numeric_limits::max() - 1; +/// Local weak-form argument slots associated with one solved field. +using BlockArgumentIndices = std::vector; + +/// Matrix mapping residual rows and solved columns to local argument slots. +using BlockIndexMap = std::vector>; + /// @brief Solve a block nonlinear system of equations as defined by the vector of weak form /// @param residual_evals Vector of weak forms which define the equations to be solved -/// @param block_indices Matrix of index arguments specifying where in each WeakForm the unknown fields are passed in. +/// @param block_indices Matrix specifying every argument slot where each solved field is passed. /// Example: for a 2 weak-form system, with weak-forms, r1, r2 /// r1(a,b,c) /// r2(b,d,e,a) // with unknowns (with respect to the solver) being a, and b. -// r1 has unknowns a,b in the ‘slots’ 0, 1 -// r2 has unknowns a,b, in the ‘slots’ 3,0 +// r1 has unknowns a,b in slots {0}, {1} +// r2 has unknowns a,b in slots {3}, {0} /// @param shape_disp The mesh-morphed shape displacement /// @param states The time varying states as inputs to the weak form /// @param params The fixed field parameters as inputs to the weak form @@ -44,9 +50,8 @@ static constexpr size_t invalid_block_index = std::numeric_limits::max() /// @param solver The nonlinear block solver used to solve the system of equations /// @param bc_managers Holds information about which degrees of freedom (DOFS) /// @return Vector of field solutions satisfying the weak forms -std::vector block_solve(const std::vector& residual_evals, - const std::vector> block_indices, const FieldState& shape_disp, - const std::vector>& states, +std::vector block_solve(const std::vector& residual_evals, const BlockIndexMap block_indices, + const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, const std::vector& bc_managers); @@ -67,8 +72,9 @@ inline FieldState solve(const WeakForm& residual_eval, const FieldState& shape_d const DirichletBoundaryConditions& bcs, size_t unknown_state_index = 0) { std::vector bc_managers{&bcs.getBoundaryConditionManager()}; - auto solutions = block_solve({const_cast(&residual_eval)}, {{unknown_state_index}}, shape_disp, {states}, - {params}, time_info, &solver, bc_managers); + BlockIndexMap block_indices{{BlockArgumentIndices{unknown_state_index}}}; + auto solutions = block_solve({const_cast(&residual_eval)}, block_indices, shape_disp, {states}, {params}, + time_info, &solver, bc_managers); return solutions[0]; } diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 658ddc8dff..042a646105 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -126,7 +126,8 @@ struct PlasticMechanicsSystem : public SystemBase { */ std::vector getReactionInfos() const { - return {{field_store->prefix("solid_force"), &field_store->getField(field_store->prefix("displacement")).get()->space()}}; + return {{field_store->prefix("solid_force"), + &field_store->getField(field_store->prefix("displacement")).get()->space()}}; } /** @@ -198,39 +199,32 @@ std::shared_ptr> bui std::shared_ptr mesh, std::shared_ptr solver, std::string prepend_name = "", FieldType... parameter_types) { - auto field_store = std::make_shared(mesh, 100); + auto field_store = std::make_shared(mesh, 100, prepend_name); - auto prefix = [&](const std::string& name) { - if (prepend_name.empty()) { - return name; - } - return prepend_name + "_" + name; - }; - - FieldType> shape_disp_type(prefix("shape_displacement")); + FieldType> shape_disp_type("shape_displacement"); field_store->addShapeDisp(shape_disp_type); auto quasistatic_time_rule = std::make_shared(); auto backward_euler_time_rule = std::make_shared(); // Displacement with quasi-static rule - FieldType> disp_type(prefix("displacement")); + FieldType> disp_type("displacement"); auto disp_bc = field_store->addIndependent(disp_type, quasistatic_time_rule); // State variable fields - FieldType> plastic_defgrad_type(prefix("plastic_defgrad")); + FieldType> plastic_defgrad_type("plastic_defgrad"); auto plastic_defgrad_bc = field_store->addIndependent(plastic_defgrad_type, backward_euler_time_rule); auto plastic_defgrad_old_type = - field_store->addDependent(plastic_defgrad_type, FieldStore::TimeDerivative::VAL, prefix("plastic_defgrad_old")); + field_store->addDependent(plastic_defgrad_type, FieldStore::TimeDerivative::VAL, "plastic_defgrad_old"); - FieldType> plastic_strain_type(prefix("plastic_strain")); + FieldType> plastic_strain_type("plastic_strain"); auto plastic_strain_bc = field_store->addIndependent(plastic_strain_type, backward_euler_time_rule); auto plastic_strain_old_type = - field_store->addDependent(plastic_strain_type, FieldStore::TimeDerivative::VAL, prefix("plastic_strain_old")); + field_store->addDependent(plastic_strain_type, FieldStore::TimeDerivative::VAL, "plastic_strain_old"); // Parameters auto register_parameter = [&](auto& parameter_type) { - parameter_type.name = prefix("param_" + parameter_type.name); + parameter_type.name = "param_" + parameter_type.name; field_store->addParameter(parameter_type); }; @@ -241,33 +235,34 @@ std::shared_ptr> bui using SystemType = PlasticMechanicsSystem; // Main weak forms for mechanics - std::string solid_res_name = prefix("solid_residual"); + std::string solid_res_name = field_store->prefix("solid_residual"); auto solid_weak_form = std::make_shared( solid_res_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), - field_store->createSpaces(solid_res_name, disp_type.name, disp_type, plastic_defgrad_type, plastic_strain_type, - parameter_types...)); + field_store->createSpaces(solid_res_name, {.unknown = disp_type.name, .test = disp_type.name}, disp_type, + plastic_defgrad_type, plastic_strain_type, parameter_types...)); - std::string plastic_defgrad_res_name = prefix("plastic_defgrad_residual"); + std::string plastic_defgrad_res_name = field_store->prefix("plastic_defgrad_residual"); auto plastic_defgrad_weak_form = std::make_shared( plastic_defgrad_res_name, field_store->getMesh(), field_store->getField(plastic_defgrad_type.name).get()->space(), - field_store->createSpaces(plastic_defgrad_res_name, plastic_defgrad_type.name, plastic_defgrad_type, - plastic_defgrad_old_type, plastic_strain_type, plastic_strain_old_type, disp_type, - parameter_types...)); + field_store->createSpaces(plastic_defgrad_res_name, + {.unknown = plastic_defgrad_type.name, .test = plastic_defgrad_type.name}, + plastic_defgrad_type, plastic_defgrad_old_type, plastic_strain_type, + plastic_strain_old_type, disp_type, parameter_types...)); - std::string plastic_strain_res_name = prefix("plastic_strain_residual"); + std::string plastic_strain_res_name = field_store->prefix("plastic_strain_residual"); auto plastic_strain_weak_form = std::make_shared( plastic_strain_res_name, field_store->getMesh(), field_store->getField(plastic_strain_type.name).get()->space(), - field_store->createSpaces(plastic_strain_res_name, plastic_strain_type.name, plastic_strain_type, - plastic_strain_old_type, plastic_defgrad_type, plastic_defgrad_old_type, disp_type, - parameter_types...)); + field_store->createSpaces(plastic_strain_res_name, + {.unknown = plastic_strain_type.name, .test = plastic_strain_type.name}, + plastic_strain_type, plastic_strain_old_type, plastic_defgrad_type, + plastic_defgrad_old_type, disp_type, parameter_types...)); // Build solver and advancer std::vector> weak_forms{solid_weak_form, plastic_defgrad_weak_form, plastic_strain_weak_form}; - auto sys = - std::make_shared(field_store, solver, weak_forms); - + auto sys = std::make_shared(field_store, solver, weak_forms); + sys->solid_weak_form = solid_weak_form; sys->plastic_deform_weak_form = plastic_defgrad_weak_form; sys->plastic_strain_weak_form = plastic_strain_weak_form; diff --git a/src/smith/differentiable_numerics/solid_mechanics_system.hpp b/src/smith/differentiable_numerics/solid_mechanics_system.hpp index cc3a885984..78cf957b91 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/solid_mechanics_system.hpp @@ -246,7 +246,7 @@ auto registerSolidMechanicsFields(std::shared_ptr field_store, auto physics_fields = PhysicsFields, H1, H1, H1>{ - field_store, FieldType>(field_store->prefix("displacement_solve_state")), + field_store, FieldType>(field_store->prefix("displacement_solve_state"), true), FieldType>(field_store->prefix("displacement")), FieldType>(field_store->prefix("velocity")), FieldType>(field_store->prefix("acceleration"))}; @@ -335,7 +335,7 @@ auto buildSolidMechanicsSystemImpl(std::shared_ptr field_store, cons auto cycle_zero_solver = detail::makeCycleZeroSolver(solver, *field_store->getMesh()); auto cycle_zero_system = makeSystem(field_store, cycle_zero_solver, {sys->solid_weak_form}); - cycle_zero_system->solve_result_field_names = {accel_old_type.name}; + cycle_zero_system->solved_field_names = {accel_old_type.name}; auto cycle_zero_inputs = std::vector{disp_old_type.name, disp_old_type.name, velo_old_type.name, accel_old_type.name}; auto append_if_state = [&](const auto& field) { diff --git a/src/smith/differentiable_numerics/state_variable_system.hpp b/src/smith/differentiable_numerics/state_variable_system.hpp index 25c237d2de..1d1be1e680 100644 --- a/src/smith/differentiable_numerics/state_variable_system.hpp +++ b/src/smith/differentiable_numerics/state_variable_system.hpp @@ -130,7 +130,7 @@ auto registerInternalVariableFields(std::shared_ptr field_store, } return PhysicsFields{ - field_store, FieldType(field_store->prefix("state_solve_state")), + field_store, FieldType(field_store->prefix("state_solve_state"), true), FieldType(field_store->prefix("state"))}; } diff --git a/src/smith/differentiable_numerics/system_base.cpp b/src/smith/differentiable_numerics/system_base.cpp index e13589eb90..ba1f5917e8 100644 --- a/src/smith/differentiable_numerics/system_base.cpp +++ b/src/smith/differentiable_numerics/system_base.cpp @@ -8,6 +8,8 @@ #include "smith/differentiable_numerics/reaction.hpp" #include "smith/differentiable_numerics/nonlinear_solve.hpp" +#include + namespace smith { std::vector SystemBase::solve(const TimeInfo& time_info) const @@ -16,44 +18,56 @@ std::vector SystemBase::solve(const TimeInfo& time_info) const for (const auto& wf : weak_forms) { weak_form_names.push_back(wf->name()); } - std::vector> index_map = field_store->indexMap(weak_form_names); + BlockIndexMap index_map; + if (solved_field_names.empty()) { + index_map = field_store->indexMap(weak_form_names); + } std::vector> inputs; + std::vector> input_field_indices; if (!solve_input_field_names.empty()) { SLIC_ERROR_IF(solve_input_field_names.size() != weak_forms.size(), "solve_input_field_names size must match weak_forms size"); } for (size_t i = 0; i < weak_forms.size(); ++i) { std::vector fields_for_wk; + std::vector indices_for_wk; if (solve_input_field_names.empty()) { fields_for_wk = field_store->getStates(weak_forms[i]->name()); } else { for (const auto& field_name : solve_input_field_names[i]) { fields_for_wk.push_back(field_store->getField(field_name)); + indices_for_wk.push_back(field_store->getFieldIndex(field_name)); } } inputs.push_back(fields_for_wk); + input_field_indices.push_back(indices_for_wk); } std::vector bc_field_names; - if (!solve_result_field_names.empty()) { - SLIC_ERROR_IF(solve_result_field_names.size() != weak_forms.size(), - "solve_result_field_names size must match weak_forms size"); - bc_field_names = solve_result_field_names; + if (!solved_field_names.empty()) { + SLIC_ERROR_IF(solved_field_names.size() != weak_forms.size(), "solved_field_names size must match weak_forms size"); + SLIC_ERROR_IF(solve_input_field_names.empty(), "solve_input_field_names must accompany solved_field_names"); + bc_field_names = solved_field_names; + index_map = BlockIndexMap(weak_forms.size(), std::vector(weak_forms.size())); + std::set result_field_indices; + for (const auto& field_name : solved_field_names) { + size_t field_index = field_store->getFieldIndex(field_name); + SLIC_ERROR_IF(!result_field_indices.insert(field_index).second, + "Each solve result field must have exactly one owning row"); + } for (size_t row = 0; row < weak_forms.size(); ++row) { for (size_t col = 0; col < weak_forms.size(); ++col) { - index_map[row][col] = invalid_block_index; - for (size_t arg = 0; arg < inputs[row].size(); ++arg) { - if (inputs[row][arg].get()->name() == solve_result_field_names[col]) { - index_map[row][col] = arg; - break; + size_t result_field_index = field_store->getFieldIndex(solved_field_names[col]); + for (size_t arg = 0; arg < input_field_indices[row].size(); ++arg) { + if (input_field_indices[row][arg] == result_field_index) { + index_map[row][col].push_back(arg); } } } - SLIC_ERROR_IF(index_map[row][row] == invalid_block_index, "Requested solve result field '" - << solve_result_field_names[row] - << "' is not an argument of weak form '" - << weak_form_names[row] << "'"); + SLIC_ERROR_IF(index_map[row][row].empty(), "Requested solve result field '" + << solved_field_names[row] << "' is not an argument of weak form '" + << weak_form_names[row] << "'"); } } @@ -64,9 +78,8 @@ std::vector SystemBase::solve(const TimeInfo& time_info) const for (auto& p : weak_forms) { weak_form_ptrs.push_back(p.get()); } - auto bc_managers = solve_result_field_names.empty() - ? field_store->getBoundaryConditionManagers(weak_form_names) - : field_store->getBoundaryConditionManagersForFields(bc_field_names); + auto bc_managers = solved_field_names.empty() ? field_store->getBoundaryConditionManagers(weak_form_names) + : field_store->getBoundaryConditionManagersForFields(bc_field_names); return solver->solve(weak_form_ptrs, index_map, field_store->getShapeDisp(), inputs, wk_params, time_info, bc_managers); } @@ -78,7 +91,7 @@ std::vector SystemBase::computeReactions(const TimeInfo& time_inf auto params = field_store->getParameterFields(); for (const auto& wf : weak_forms) { std::vector wf_fields = field_store->getStatesFromVectors(wf->name(), states_for_reactions, params); - std::string test_field_name = field_store->getWeakFormReaction(wf->name()); + std::string test_field_name = field_store->getWeakFormFieldNames(wf->name()).test; size_t test_field_idx = field_store->getFieldIndex(test_field_name); FieldState test_field = states_for_reactions[test_field_idx]; reactions.push_back(smith::evaluateWeakForm(wf, time_info, field_store->getShapeDisp(), wf_fields, test_field)); diff --git a/src/smith/differentiable_numerics/system_base.hpp b/src/smith/differentiable_numerics/system_base.hpp index 4a4d9e7ec1..395f26b0f6 100644 --- a/src/smith/differentiable_numerics/system_base.hpp +++ b/src/smith/differentiable_numerics/system_base.hpp @@ -61,7 +61,7 @@ struct SystemBase { ///< independently; cycle-zero solves do not couple across subsystems. std::vector> post_solve_systems; ///< Optional systems solved after main state update. std::vector - solve_result_field_names; ///< Optional per-weak-form fields to solve/update instead of reaction fields. + solved_field_names; ///< Optional per-weak-form fields overriding registered solver-owned fields. std::vector> solve_input_field_names; ///< Optional per-weak-form input field ordering used during solve. @@ -121,10 +121,10 @@ auto buildWeakFormWithCouplingImpl(const FieldStorePtr& field_store, const std:: auto coupling_fields_tuple = std::get(args_tuple); return std::apply( [&](const auto&... coupling_fields) { - return std::make_shared(weak_form_name, field_store->getMesh(), - field_store->getField(unknown_field_name).get()->space(), - field_store->createSpaces(weak_form_name, unknown_field_name, - std::get(args_tuple)..., coupling_fields...)); + return std::make_shared( + weak_form_name, field_store->getMesh(), field_store->getField(unknown_field_name).get()->space(), + field_store->createSpaces(weak_form_name, {.unknown = unknown_field_name, .test = unknown_field_name}, + std::get(args_tuple)..., coupling_fields...)); }, coupling_fields_tuple); } diff --git a/src/smith/differentiable_numerics/system_solver.cpp b/src/smith/differentiable_numerics/system_solver.cpp index 338bcfedf6..883cfd0370 100644 --- a/src/smith/differentiable_numerics/system_solver.cpp +++ b/src/smith/differentiable_numerics/system_solver.cpp @@ -13,8 +13,6 @@ #include #include -#include -#include #include #include @@ -65,8 +63,7 @@ void SystemSolver::appendStagesWithBlockMapping(const SystemSolver& subsystem_so } std::vector SystemSolver::solve(const std::vector& residual_evals, - const std::vector>& block_indices, - const FieldState& shape_disp, + const BlockIndexMap& block_indices, const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, @@ -75,8 +72,21 @@ std::vector SystemSolver::solve(const std::vector& residu SLIC_ERROR_IF(stages_.empty(), "SystemSolver has no stages defined."); size_t num_residuals = residual_evals.size(); - // const bool trace_plasticity = - // num_residuals == 3 && residual_evals[1]->name().find("plastic_defgrad_residual") != std::string::npos; + SLIC_ERROR_IF(states.size() != num_residuals, "State rows must match residual count"); + SLIC_ERROR_IF(block_indices.size() != num_residuals, "Block index rows must match residual count"); + for (size_t row = 0; row < num_residuals; ++row) { + SLIC_ERROR_IF(block_indices[row].size() != num_residuals, + axom::fmt::format("Block index columns for row {} must match residual count", row)); + for (size_t col = 0; col < num_residuals; ++col) { + for (size_t slot : block_indices[row][col]) { + SLIC_ERROR_IF(slot >= states[row].size(), + axom::fmt::format("Block slot {} for row {}, column {} exceeds state count {}", slot, row, col, + states[row].size())); + } + } + SLIC_ERROR_IF(block_indices[row][row].empty(), + axom::fmt::format("Block index diagonal for row {} must not be empty", row)); + } std::vector active_stages = stages_; for (auto& stage : active_stages) { if (stage.block_indices.empty()) { @@ -110,24 +120,6 @@ std::vector SystemSolver::solve(const std::vector& residu // Working copy of states, updated in-place as stages solve std::vector> current_states = states; - // auto trace_row = [&](const char* location, size_t row) { - // for (size_t slot = 0; slot < current_states[row].size(); ++slot) { - // const auto& state = current_states[row][slot]; - // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] {} row={} slot={} name={} norm={:.16e}", location, row, - // slot, state.get()->name(), state.get()->Norml2())); - // } - // }; - - // Pre-compute name -> (row, slot) routing so the propagation loop avoids O(N*M) string compares - // on every staggered iteration. Field-name identity within current_states is invariant across - // the iteration loop: only values are replaced, never the underlying name. - std::unordered_map>> field_routing; - for (size_t r = 0; r < num_residuals; ++r) { - for (size_t slot = 0; slot < current_states[r].size(); ++slot) { - field_routing[current_states[r][slot].get()->name()].emplace_back(r, slot); - } - } - // Helper lambda to assemble input pointers, evaluate residual, and zero essential BCs auto eval_residual_and_zero_bcs = [&](size_t global_row) { std::vector input_ptrs; @@ -172,7 +164,7 @@ std::vector SystemSolver::solve(const std::vector& residu size_t num_stage_blocks = stage.block_indices.size(); std::vector stage_residuals; - std::vector> stage_block_indices; + BlockIndexMap stage_block_indices; std::vector> stage_states; std::vector> stage_params; std::vector stage_bc_managers; @@ -184,7 +176,7 @@ std::vector SystemSolver::solve(const std::vector& residu stage_states.push_back(current_states[global_row]); stage_params.push_back(params[global_row]); - std::vector row_indices(num_stage_blocks, invalid_block_index); + std::vector row_indices(num_stage_blocks); for (size_t col_idx = 0; col_idx < num_stage_blocks; ++col_idx) { size_t global_col = stage.block_indices[col_idx]; row_indices[col_idx] = block_indices[global_row][global_col]; @@ -192,57 +184,26 @@ std::vector SystemSolver::solve(const std::vector& residu stage_block_indices.push_back(row_indices); } - // // The plastic mechanics system orders its blocks as displacement, plastic deformation - // // gradient, and plastic strain. Trace the deformation-gradient inputs after the strain - // // stage should have propagated its updated state. - // if (trace_plasticity && stage.block_indices.size() == 1 && stage.block_indices[0] == 1) { - // trace_row("before Fp solve", 1); - // auto residual = eval_residual_and_zero_bcs(1); - // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] Fp residual norm before solve={:.16e}", - // residual.Norml2())); - // } - std::vector stage_solutions = block_solve(stage_residuals, stage_block_indices, shape_disp, stage_states, stage_params, time_info, stage.solver.get(), stage_bc_managers); - // Propagate updated fields to every residual input that references the solved field. - // Match by field name (looked up via the pre-computed routing map): coupling fields appear - // as fixed inputs in other rows and therefore do not have a valid unknown-block entry there. + // Propagate updated fields to every residual input sharing the solved FieldStore index. // Apply relaxation: x_new = omega * x_solved + (1 - omega) * x_k. for (size_t i = 0; i < num_stage_blocks; ++i) { size_t global_col = stage.block_indices[i]; FieldState new_state = stage_solutions[i]; if (stage.relaxation_factor != 1.0) { - FieldState old_state = current_states[global_col][block_indices[global_col][global_col]]; + FieldState old_state = current_states[global_col][block_indices[global_col][global_col].front()]; new_state = weighted_average(new_state, old_state, stage.relaxation_factor); } - auto it = field_routing.find(new_state.get()->name()); - // if (trace_plasticity && global_col == 2) { - // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] solved strain name={} norm={:.16e}", - // new_state.get()->name(), new_state.get()->Norml2())); - // if (it == field_routing.end()) { - // SLIC_INFO_ROOT("[plastic-debug] no routing entry for solved strain"); - // } else { - // for (const auto& [row, slot] : it->second) { - // SLIC_INFO_ROOT(axom::fmt::format("[plastic-debug] strain routes to row={} slot={}", row, slot)); - // } - // } - // } - if (it != field_routing.end()) { - for (const auto& [r, slot] : it->second) { - current_states[r][slot] = new_state; + for (size_t row = 0; row < num_residuals; ++row) { + for (size_t slot : block_indices[row][global_col]) { + current_states[row][slot] = new_state; } } - // if (trace_plasticity && global_col == 2) { - // trace_row("after strain routing", 1); - // } - // if (trace_plasticity && global_col == 1) { - // SLIC_INFO_ROOT( - // axom::fmt::format("[plastic-debug] solved Fp norm={:.16e}", new_state.get()->Norml2())); - // } } } @@ -275,7 +236,7 @@ std::vector SystemSolver::solve(const std::vector& residu std::vector final_solutions; final_solutions.reserve(num_residuals); for (size_t r = 0; r < num_residuals; ++r) { - size_t s_idx = block_indices[r][r]; + size_t s_idx = block_indices[r][r].front(); final_solutions.push_back(current_states[r][s_idx]); } diff --git a/src/smith/differentiable_numerics/system_solver.hpp b/src/smith/differentiable_numerics/system_solver.hpp index 71351f1c53..fabca974d2 100644 --- a/src/smith/differentiable_numerics/system_solver.hpp +++ b/src/smith/differentiable_numerics/system_solver.hpp @@ -10,6 +10,7 @@ #include #include #include "smith/differentiable_numerics/field_state.hpp" +#include "smith/differentiable_numerics/nonlinear_solve.hpp" #include "smith/numerics/solver_config.hpp" #include "smith/physics/common.hpp" @@ -59,16 +60,15 @@ class SystemSolver { /// @brief Solves the multiphysics system using staggered iterations. /// @param residual_evals Vector of WeakForm evaluations for each block. - /// @param block_indices Block indices for each residual evaluation. + /// @param block_indices Argument slot lists for each residual and solved field. /// @param shape_disp Current shape displacement. /// @param states Nested vector of field states. /// @param params Nested vector of parameters. /// @param time_info Current time information. /// @param bc_managers Managers for boundary conditions. /// @return Updated field states. - std::vector solve(const std::vector& residual_evals, - const std::vector>& block_indices, const FieldState& shape_disp, - const std::vector>& states, + std::vector solve(const std::vector& residual_evals, const BlockIndexMap& block_indices, + const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const std::vector& bc_managers) const; diff --git a/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp index 856d5e42c0..0a2575753e 100644 --- a/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp @@ -514,7 +514,7 @@ TEST_F(ThermoMechanicsMeshFixture, CombinedSystemCarriesCycleZeroSystems) ASSERT_EQ(solid_system->cycle_zero_systems.size(), 1u); EXPECT_EQ(combined_system->cycle_zero_systems.size(), solid_system->cycle_zero_systems.size()); - EXPECT_EQ(combined_system->cycle_zero_systems[0]->solve_result_field_names, std::vector{"acceleration"}); + EXPECT_EQ(combined_system->cycle_zero_systems[0]->solved_field_names, std::vector{"acceleration"}); } } // namespace smith diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index 057968f19f..e97bf43027 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -164,6 +164,14 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); auto states = plastic_mechanics_system->field_store->getStateFields(); auto params = plastic_mechanics_system->field_store->getParameterFields(); + const double initial_plastic_defgrad_l2 = norm( + *plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_defgrad")) + .get(), + 2.0); + const double initial_plastic_strain_l2 = norm( + *plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_strain")) + .get(), + 2.0); auto advancer = makeAdvancer(plastic_mechanics_system); @@ -172,16 +180,21 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) std::vector reactions; for (size_t step = 0; step < 10; ++step) { - std::tie(states, reactions) = - advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); + std::tie(states, reactions) = advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); time += dt; pv_writer.write(step + 1, time, states); } - auto displacement = plastic_mechanics_system->field_store->getField( - plastic_mechanics_system->field_store->prefix("displacement")); + auto displacement = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); + auto plastic_defgrad = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_defgrad")); + auto plastic_strain = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_strain")); double final_disp_l2 = norm(*displacement.get(), 2.0); EXPECT_GT(final_disp_l2, 0.0); + EXPECT_GT(norm(*plastic_defgrad.get(), 2.0), initial_plastic_defgrad_l2 + 1.0e-10); + EXPECT_GT(norm(*plastic_strain.get(), 2.0), initial_plastic_strain_l2 + 1.0e-10); } /// @brief Parameterized differentiable J2 material with nonlinear isotropic hardening and linear kinematic hardening @@ -333,8 +346,8 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); auto params = plastic_mechanics_system->field_store->getParameterFields(); - auto displacement = plastic_mechanics_system->field_store->getField( - plastic_mechanics_system->field_store->prefix("displacement")); + auto displacement = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); auto disp_squared = 0.5 * innerProduct(displacement, displacement); gretl::set_as_objective(disp_squared); diff --git a/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp b/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp index 9c21313fd3..e9f8e5b61c 100644 --- a/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp +++ b/src/smith/differentiable_numerics/tests/test_mixed_poisson.cpp @@ -216,7 +216,7 @@ TEST_P(BlockPreconditionerTest, BlockSolve) auto nonlinear_block_solver = smith::buildNonlinearBlockSolver(nonlin_opts, linear_options, *mesh); - auto sols = block_solve({&con_form, &bal_form}, {{0, 1}, {0, 1}}, shape_disp, {con_arguments, bal_arguments}, + auto sols = block_solve({&con_form, &bal_form}, {{{0}, {1}}, {{0}, {1}}}, shape_disp, {con_arguments, bal_arguments}, {flux_params, potential_params}, smith::TimeInfo(time.get(), dt.get(), cycle), nonlinear_block_solver.get(), {flux_bc_manager.get(), potential_bc_manager.get()}); diff --git a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp index c2cb034247..fdfb0d838e 100644 --- a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp @@ -121,6 +121,24 @@ class ConstantNonlinearBlockSolver : public NonlinearBlockSolverBase { mutable int solve_calls_ = 0; }; +class ResidualRecordingNonlinearBlockSolver : public NoOpNonlinearBlockSolver { + public: + std::vector solve( + const std::vector& u_guesses, + std::function(const std::vector&)> residual, + std::function>(const std::vector&)> jacobian) const override + { + auto residuals = residual(u_guesses); + residual_norm_ = residuals.front().Norml2(); + return NoOpNonlinearBlockSolver::solve(u_guesses, std::move(residual), std::move(jacobian)); + } + + double residualNorm() const { return residual_norm_; } + + private: + mutable double residual_norm_ = 0.0; +}; + class NeedsInitialSolveRule : public QuasiStaticRule { public: bool requiresInitialAccelerationSolve() const override { return true; } @@ -131,8 +149,9 @@ auto buildScalarDiffusionWeakForm(const std::string& name, std::shared_ptr FieldTypeT field_type) { using WeakFormType = FunctionalWeakForm<2, H1<1>, Parameters>>; - auto weak_form = std::make_shared(name, mesh, fs->getField(field_type.name).get()->space(), - fs->createSpaces(name, field_type.name, field_type)); + auto weak_form = std::make_shared( + name, mesh, fs->getField(field_type.name).get()->space(), + fs->createSpaces(name, {.unknown = field_type.name, .test = field_type.name}, field_type)); weak_form->addBodyIntegral(mesh->entireBodyName(), [](auto, auto, auto u) { return smith::tuple{0.0 * get(u), get(u)}; }); return weak_form; @@ -144,10 +163,10 @@ auto buildSecondOrderMainWeakForm(const std::string& name, std::shared_ptr VelocityFieldType velocity_type, AccelerationFieldType acceleration_type) { using WeakFormType = FunctionalWeakForm<2, H1<1>, Parameters, H1<1>, H1<1>, H1<1>>>; - auto weak_form = - std::make_shared(name, mesh, fs->getField(displacement_type.name).get()->space(), - fs->createSpaces(name, displacement_type.name, displacement_type, - displacement_old_type, velocity_type, acceleration_type)); + auto weak_form = std::make_shared( + name, mesh, fs->getField(displacement_type.name).get()->space(), + fs->createSpaces(name, {.unknown = displacement_type.name, .test = displacement_type.name}, displacement_type, + displacement_old_type, velocity_type, acceleration_type)); weak_form->addBodySource(mesh->entireBodyName(), [](auto, auto, auto...) { return 0.0; }); return weak_form; } @@ -160,7 +179,8 @@ auto buildSecondOrderCycleZeroWeakForm(const std::string& name, std::shared_ptr< using WeakFormType = FunctionalWeakForm<2, H1<1>, Parameters, H1<1>, H1<1>>>; auto weak_form = std::make_shared( name, mesh, fs->getField(acceleration_type.name).get()->space(), - fs->createSpaces(name, acceleration_type.name, displacement_type, velocity_type, acceleration_type)); + fs->createSpaces(name, {.unknown = acceleration_type.name, .test = acceleration_type.name}, displacement_type, + velocity_type, acceleration_type)); weak_form->addBodySource(mesh->entireBodyName(), [](auto, auto, auto...) { return 0.0; }); return weak_form; } @@ -173,7 +193,7 @@ bool allNodesOnBoundary(const std::vector& nodes, double x_target) } // namespace -TEST(MultiphysicsTimeIntegrator, CycleZeroUsesBcsForReactionFieldNotUnknownZero) +TEST(MultiphysicsTimeIntegrator, CycleZeroUsesOwnedFieldBoundaryConditions) { axom::sidre::DataStore datastore; StateManager::initialize(datastore, "multiphysics_time_integrator"); @@ -335,7 +355,7 @@ TEST(MultiphysicsTimeIntegrator, CycleZeroSolveResultUpdatesAccelerationField) cycle_zero_system->field_store = field_store; cycle_zero_system->weak_forms = {cycle_zero_wf}; cycle_zero_system->solver = cycle_zero_solver; - cycle_zero_system->solve_result_field_names = {acceleration_type.name}; + cycle_zero_system->solved_field_names = {acceleration_type.name}; cycle_zero_system->solve_input_field_names = { {displacement_old_type.name, displacement_old_type.name, velocity_type.name, acceleration_type.name}}; @@ -514,6 +534,101 @@ TEST(SystemSolver, AppendsStagesWithBlockMappingForCombinedSubsystems) StateManager::reset(); } +TEST(SystemSolver, RoutesRepeatedDependenciesWithoutUpdatingHistory) +{ + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "system_solver_repeated_dependency_routing"); + + auto mesh = std::make_shared(mfem::Mesh::MakeCartesian2D(2, 2, mfem::Element::QUADRILATERAL, true, 1.0, 1.0), + "system_solver_repeated_dependency_routing_mesh"); + auto field_store = std::make_shared(mesh, 20); + FieldType> shape_disp_type("shape_displacement"); + field_store->addShapeDisp(shape_disp_type); + + auto static_rule = std::make_shared(); + auto backward_euler_rule = std::make_shared(); + FieldType> field_a_type("field_a"); + FieldType> field_b_type("field_b"); + field_store->addIndependent(field_a_type, backward_euler_rule); + auto field_a_old_type = field_store->addDependent(field_a_type, FieldStore::TimeDerivative::VAL, "field_a_old"); + field_store->addIndependent(field_b_type, static_rule); + + auto field_a_weak_form = buildScalarDiffusionWeakForm("field_a_residual", mesh, field_store, field_a_type); + using FieldBWeakForm = FunctionalWeakForm<2, H1<1>, Parameters, H1<1>, H1<1>, H1<1>>>; + auto field_b_weak_form = std::make_shared( + "field_b_residual", mesh, field_store->getField(field_b_type.name).get()->space(), + field_store->createSpaces("field_b_residual", {.unknown = field_b_type.name, .test = field_b_type.name}, + field_b_type, field_a_type, field_a_type, field_a_old_type)); + field_b_weak_form->addBodySource(mesh->entireBodyName(), + [](auto, auto, auto, auto field_a_first, auto field_a_second, auto field_a_old) { + return field_a_first + field_a_second - 2.0 * field_a_old; + }); + + field_store->setField(field_store->getFieldIndex(field_a_old_type.name), field_store->getField(field_a_type.name)); + + const std::vector residual_names{field_a_weak_form->name(), field_b_weak_form->name()}; + const auto block_indices = field_store->indexMap(residual_names); + EXPECT_EQ(block_indices[0][0], BlockArgumentIndices({0})); + EXPECT_TRUE(block_indices[0][1].empty()); + EXPECT_EQ(block_indices[1][0], BlockArgumentIndices({1, 2})); + EXPECT_EQ(block_indices[1][1], BlockArgumentIndices({0})); + + auto field_a_solver = std::make_shared(3.0); + auto field_b_solver = std::make_shared(); + SystemSolver solver(1, true); + solver.addSubsystemSolver({0}, field_a_solver, 0.5); + solver.addSubsystemSolver({1}, field_b_solver); + + const std::vector residuals{field_a_weak_form.get(), field_b_weak_form.get()}; + const std::vector> states{field_store->getStates(field_a_weak_form->name()), + field_store->getStates(field_b_weak_form->name())}; + const std::vector> params(residuals.size()); + const auto bc_managers = field_store->getBoundaryConditionManagers(residual_names); + const TimeInfo time_info(0.0, 1.0, 0); + auto solutions = + solver.solve(residuals, block_indices, field_store->getShapeDisp(), states, params, time_info, bc_managers); + + EXPECT_NEAR((*solutions[0].get())[0], 1.5, 1.0e-12); + + auto expected_field_b_states = states[1]; + expected_field_b_states[1] = solutions[0]; + expected_field_b_states[2] = solutions[0]; + auto expected_residual = field_b_weak_form->residual(time_info, field_store->getShapeDisp().get().get(), + getConstFieldPointers(expected_field_b_states)); + + EXPECT_NEAR(field_b_solver->residualNorm(), expected_residual.Norml2(), 1.0e-12); + EXPECT_GT(expected_residual.Norml2(), 1.0e-8); + + StateManager::reset(); +} + +TEST(FieldStore, UsesExplicitUnknownFieldInsteadOfTestField) +{ + axom::sidre::DataStore datastore; + StateManager::initialize(datastore, "field_store_explicit_unknown_field"); + + auto mesh = std::make_shared(mfem::Mesh::MakeCartesian2D(2, 2, mfem::Element::QUADRILATERAL, true, 1.0, 1.0), + "field_store_explicit_unknown_field_mesh"); + auto field_store = std::make_shared(mesh, 10); + FieldType> test_type("test_field"); + FieldType> unknown_type("unknown_field"); + auto static_rule = std::make_shared(); + field_store->addIndependent(test_type, static_rule); + auto unknown_bcs = field_store->addIndependent(unknown_type, static_rule); + + field_store->createSpaces("residual", {.unknown = unknown_type.name, .test = test_type.name}, test_type, + unknown_type); + + const auto block_indices = field_store->indexMap({"residual"}); + EXPECT_EQ(field_store->getWeakFormFieldNames("residual").test, test_type.name); + EXPECT_EQ(field_store->getWeakFormFieldNames("residual").unknown, unknown_type.name); + EXPECT_EQ(block_indices[0][0], BlockArgumentIndices({1})); + EXPECT_EQ(field_store->getBoundaryConditionManagers({"residual"})[0], &unknown_bcs->getBoundaryConditionManager()); + EXPECT_EQ(field_store->getReactionInfos()[0].space, &field_store->getField(test_type.name).get()->space()); + + StateManager::reset(); +} + } // namespace smith int main(int argc, char* argv[]) diff --git a/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp b/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp index af2b9429d2..b52a69b763 100644 --- a/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp +++ b/src/smith/differentiable_numerics/tests/test_porous_heat_sink.cpp @@ -303,7 +303,7 @@ TEST_P(BlockPreconditionerTest, BlockSolve) auto nonlinear_block_solver = smith::buildNonlinearBlockSolver(nonlin_opts, linear_options, *mesh); - auto sols = block_solve({&T1_form, &T2_form}, {{0, 1}, {0, 1}}, shape_disp, {T1_arguments, T2_arguments}, + auto sols = block_solve({&T1_form, &T2_form}, {{{0}, {1}}, {{0}, {1}}}, shape_disp, {T1_arguments, T2_arguments}, {T1_params, T2_params}, smith::TimeInfo(time.get(), dt.get(), cycle), nonlinear_block_solver.get(), {T1_bc_manager.get(), T2_bc_manager.get()}); diff --git a/src/smith/differentiable_numerics/thermal_system.hpp b/src/smith/differentiable_numerics/thermal_system.hpp index 1a7f54baee..d0845857c9 100644 --- a/src/smith/differentiable_numerics/thermal_system.hpp +++ b/src/smith/differentiable_numerics/thermal_system.hpp @@ -194,7 +194,7 @@ auto registerThermalFields(std::shared_ptr field_store, field_store->addDependent(temperature_type, FieldStore::TimeDerivative::VAL, "temperature"); return PhysicsFields, H1>{ - field_store, FieldType>(field_store->prefix("temperature_solve_state")), + field_store, FieldType>(field_store->prefix("temperature_solve_state"), true), FieldType>(field_store->prefix("temperature"))}; } diff --git a/src/smith/physics/functional_weak_form.hpp b/src/smith/physics/functional_weak_form.hpp index 926660a66d..f00e3245e1 100644 --- a/src/smith/physics/functional_weak_form.hpp +++ b/src/smith/physics/functional_weak_form.hpp @@ -333,6 +333,9 @@ class FunctionalWeakForm, [[maybe_unused]] const std::vector& quad_fields = {}) const override { validateFields(fields, "jacobian"); + SLIC_ERROR_IF(jacobian_weights.size() != fields.size(), + axom::fmt::format("jacobian(): jacobian_weights.size()={} but fields.size()={}", + jacobian_weights.size(), fields.size())); SetCurrentTimeInfoRAII clear_current_time_info_on_exit(current_time_info_, time_info); std::unique_ptr J; diff --git a/src/smith/physics/tests/test_functional_weak_form.cpp b/src/smith/physics/tests/test_functional_weak_form.cpp index 758ee9e911..ba14f854d1 100644 --- a/src/smith/physics/tests/test_functional_weak_form.cpp +++ b/src/smith/physics/tests/test_functional_weak_form.cpp @@ -103,6 +103,8 @@ struct WeakFormFixture : public testing::Test { [](auto /*t_info*/, auto /*x*/, auto n, auto... /*args*/) { return 1.0 * n; }); f_weak_form->addBodySource(mesh->entireBodyName(), [](auto /*t_info*/, auto /*x*/, auto u, auto... /*args*/) { return u; }); + f_weak_form->addBodySource(mesh->entireBodyName(), [](auto /*t_info*/, auto /*x*/, auto /*u*/, auto velocity, + auto /*density*/) { return 0.25 * velocity; }); f_weak_form->addBodySource(mesh->entireBodyName(), [](auto /*t_info*/, auto x, auto... /*args*/) { return 0.5 * x; }); @@ -250,6 +252,30 @@ TEST_F(WeakFormFixture, JvpConsistency) } } +TEST_F(WeakFormFixture, MultipleJacobianWeightsSumArgumentDerivatives) +{ + auto input_fields = getConstFieldPointers(states, params); + const std::vector displacement_weights = {1.0, 0.0, 0.0}; + const std::vector velocity_weights = {0.0, 1.0, 0.0}; + const std::vector combined_weights = {1.0, 1.0, 0.0}; + + auto displacement_jacobian = weak_form->jacobian(time_info, shape_disp.get(), input_fields, displacement_weights); + auto velocity_jacobian = weak_form->jacobian(time_info, shape_disp.get(), input_fields, velocity_weights); + auto combined_jacobian = weak_form->jacobian(time_info, shape_disp.get(), input_fields, combined_weights); + + smith::FiniteElementDual expected(states[DISP].space(), "expected"); + smith::FiniteElementDual velocity_contribution(states[DISP].space(), "velocity_contribution"); + smith::FiniteElementDual actual(states[DISP].space(), "actual"); + displacement_jacobian->Mult(state_tangents[DISP], expected); + velocity_jacobian->Mult(state_tangents[DISP], velocity_contribution); + expected += velocity_contribution; + combined_jacobian->Mult(state_tangents[DISP], actual); + + EXPECT_GT(velocity_contribution.Norml2(), 1.0e-12); + actual -= expected; + EXPECT_NEAR(actual.Norml2(), 0.0, 1.0e-12); +} + TEST_F(WeakFormFixture, ForwardsOriginalTimeInfoToIntegrands) { using TrialSpace = VectorSpace; From 6891a60032a107bfcdbf8db12742996a26454cbd Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 24 Aug 2026 17:02:27 -0700 Subject: [PATCH 10/28] Renames, remove is unknown, and other simplifications. --- .../differentiable_numerics/field_store.cpp | 35 +++-------- .../differentiable_numerics/field_store.hpp | 62 +++---------------- .../multiphysics_time_integrator.cpp | 10 +-- .../nonlinear_solve.cpp | 5 +- .../nonlinear_solve.hpp | 7 ++- .../solid_mechanics_system.hpp | 8 +-- .../state_variable_system.hpp | 4 +- .../differentiable_numerics/system_base.cpp | 29 +++++---- .../differentiable_numerics/system_base.hpp | 2 +- .../differentiable_numerics/system_solver.cpp | 4 +- .../differentiable_numerics/system_solver.hpp | 2 +- .../tests/test_combined_thermo_mechanics.cpp | 2 +- .../tests/test_differentiable_plasticity.cpp | 11 ++-- .../test_multiphysics_time_integrator.cpp | 5 +- .../thermal_system.hpp | 4 +- 15 files changed, 65 insertions(+), 125 deletions(-) diff --git a/src/smith/differentiable_numerics/field_store.cpp b/src/smith/differentiable_numerics/field_store.cpp index 80cda52ae0..a975cce84c 100644 --- a/src/smith/differentiable_numerics/field_store.cpp +++ b/src/smith/differentiable_numerics/field_store.cpp @@ -9,8 +9,6 @@ #include "smith/differentiable_numerics/nonlinear_solve.hpp" #include "gretl/wang_checkpoint_strategy.hpp" -#include - namespace smith { FieldStore::FieldStore(std::shared_ptr mesh, size_t storage_size, std::string prepend_name) @@ -33,16 +31,6 @@ std::shared_ptr FieldStore::addBoundaryConditions(F return std::make_shared(mesh_->mfemParMesh(), field->space()); } -void FieldStore::addWeakFormUnknownArg(std::string weak_form_name, std::string argument_name, size_t argument_index) -{ - FieldLabel argument_name_and_index{.field_name = argument_name, .field_index_in_residual = argument_index}; - if (weak_form_name_to_unknown_name_index_.count(weak_form_name)) { - weak_form_name_to_unknown_name_index_.at(weak_form_name).push_back(argument_name_and_index); - } else { - weak_form_name_to_unknown_name_index_[weak_form_name] = std::vector{argument_name_and_index}; - } -} - void FieldStore::addWeakFormArg(std::string weak_form_name, std::string argument_name, size_t argument_index) { // Store the field name instead of index to avoid confusion between states_ and params_ indices @@ -57,16 +45,16 @@ void FieldStore::addWeakFormArg(std::string weak_form_name, std::string argument void FieldStore::printMap() { - for (auto& keyval : weak_form_name_to_unknown_name_index_) { + for (const auto& keyval : weak_form_name_to_field_names_) { std::cout << "for residual: " << keyval.first << " "; - for (auto& name_index : keyval.second) { - std::cout << "arg " << name_index.field_name << " at " << name_index.field_index_in_residual << ", "; + for (size_t argument_index = 0; argument_index < keyval.second.size(); ++argument_index) { + std::cout << "arg " << keyval.second[argument_index] << " at " << argument_index << ", "; } std::cout << std::endl; } } -BlockIndexMap FieldStore::indexMap(const std::vector& residual_names) const +BlockArgumentMap FieldStore::indexMap(const std::vector& residual_names) const { std::map global_state_to_local_col; for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { @@ -77,15 +65,13 @@ BlockIndexMap FieldStore::indexMap(const std::vector& residual_name residual_names[owner->second], res_name, global_state_idx)); } - BlockIndexMap block_indices(residual_names.size(), - std::vector(residual_names.size(), BlockArgumentIndices{})); + BlockArgumentMap block_indices(residual_names.size(), + std::vector(residual_names.size(), BlockArgumentIndices{})); for (size_t res_i = 0; res_i < residual_names.size(); ++res_i) { const std::string& res_name = residual_names[res_i]; const auto& field_names = weak_form_name_to_field_names_.at(res_name); - const auto& dependencies = weak_form_name_to_unknown_name_index_.at(res_name); size_t state_slot = 0; - for (size_t argument_index = 0; argument_index < field_names.size(); ++argument_index) { - const auto& field_name = field_names[argument_index]; + for (const auto& field_name : field_names) { if (!to_states_index_.count(field_name)) { continue; } @@ -93,13 +79,6 @@ BlockIndexMap FieldStore::indexMap(const std::vector& residual_name size_t global_state_index = to_states_index_.at(field_name); auto column = global_state_to_local_col.find(global_state_index); if (column != global_state_to_local_col.end()) { - bool marked_dependency = std::any_of(dependencies.begin(), dependencies.end(), [&](const auto& dependency) { - return dependency.field_index_in_residual == argument_index && dependency.field_name == field_name; - }); - SLIC_ERROR_IF(!marked_dependency, - axom::fmt::format("Field '{}' in weak form '{}' belongs to solved column {} but is not marked " - "as an unknown dependency", - field_name, res_name, column->second)); block_indices[res_i][column->second].push_back(state_slot); } ++state_slot; diff --git a/src/smith/differentiable_numerics/field_store.hpp b/src/smith/differentiable_numerics/field_store.hpp index 939aef4205..a76adfb41b 100644 --- a/src/smith/differentiable_numerics/field_store.hpp +++ b/src/smith/differentiable_numerics/field_store.hpp @@ -38,12 +38,7 @@ struct UnknownAndTestFieldNames { }; /** - * @brief Representation of a field type with a name and dependency flag. - * - * @c is_unknown is intentionally a per-instance flag. It marks every occurrence whose value - * depends on a field solved within the same system. Solver stages determine which dependencies - * are active during each block solve. History and external fields remain fixed. - * + * @brief Representation of a named field type. * @tparam Space The finite element space type. * @tparam Time The time integration type (unused by default). */ @@ -51,14 +46,9 @@ template struct FieldType { using space_type = Space; ///< The finite element space type. - /** - * @brief Construct a new FieldType object. - * @param n Name of the field. - * @param is_unknown_ Whether this field depends on a system-solved unknown. - */ - FieldType(std::string n, bool is_unknown_ = false) : name(n), is_unknown(is_unknown_) {} - std::string name; ///< Name of the field. - bool is_unknown = false; ///< True when this argument depends on a system-solved unknown. + /// Construct a field type with the given name. + FieldType(std::string n) : name(n) {} + std::string name; ///< Name of the field. }; /** @@ -131,13 +121,10 @@ struct FieldStore { /** * @brief Add an independent field (a solver unknown) to the store. * - * Registers the field as an unknown by setting @c type.is_unknown = true, so the same - * @c FieldType object can later be passed to @c createSpaces to mark this argument - * as a solved dependency. Also creates a boundary-condition slot keyed by field - * name that callers can populate after this call returns. + * Also creates a boundary-condition slot keyed by field name that callers can populate. * * @tparam Space The finite element space type. - * @param type The field type specification; @c type.is_unknown is set to @c true on return. + * @param type The field type specification. * @param time_rule The time integration rule governing how this unknown and its dependents * are related across time steps. * @return std::shared_ptr The boundary conditions for this field. @@ -147,7 +134,6 @@ struct FieldStore { std::shared_ptr time_rule) { type.name = prefix(type.name); - type.is_unknown = true; to_states_index_[type.name] = states_.size(); FieldState new_field = smith::createFieldState(*graph_, Space{}, type.name, mesh_->tag()); states_.push_back(new_field); @@ -174,12 +160,8 @@ struct FieldStore { * time integration rule can reconstruct the current rate from the pair * (predicted_value, stored_old_value). * - * **Return value:** a @c FieldType whose @c name is the name of the newly registered - * field and @c is_unknown is @c false. This object is intentionally returned (rather than - * discarded) so that callers can pass it directly to @c createSpaces when assembling the - * weak-form argument list. To make it a solved dependency for a specific weak form (e.g. - * the cycle-zero acceleration solve), copy the returned object and set @c is_unknown = true - * before passing it to @c createSpaces. + * Returns a descriptor for the newly registered field so callers can pass it directly to + * @c createSpaces when assembling the weak-form argument list. * * @tparam Space The finite element space type (must match the independent field). * @param independent_field The @c FieldType of the independent (predicted) field. @@ -224,14 +206,6 @@ struct FieldStore { return FieldType(name); } - /** - * @brief Register an argument to a weak form as an unknown. - * @param weak_form_name Name of the weak form. - * @param argument_name Name of the argument field. - * @param argument_index Index of the argument in the weak form's argument list. - */ - void addWeakFormUnknownArg(std::string weak_form_name, std::string argument_name, size_t argument_index); - /** * @brief Register an argument to a weak form. * @param weak_form_name Name of the weak form. @@ -259,15 +233,9 @@ struct FieldStore { * * This is the primary setup method for constructing a weak form. It: * 1. Registers explicit test-space and solver-owned field names. - * 2. Iterates over every @c FieldType in @p types (in order), registering each as an input - * argument. The solver-owned field and every argument with @c type.is_unknown are registered - * as dependencies on system-solved unknowns. + * 2. Registers every @c FieldType in @p types as an ordered input argument. * 3. Returns the ordered vector of finite element spaces. * - * Every occurrence of a field solved by another row in the same system must have - * @c is_unknown = true. Solver stages choose active Jacobian columns. Fixed history fields and - * fields owned by external systems remain false. - * * @param weak_form_name Name of the weak form being constructed. * @param field_names Names of the test-space and solver-owned fields. * @param types Ordered list of @c FieldType descriptors for every input argument. @@ -285,9 +253,6 @@ struct FieldStore { auto register_field = [&](auto type) { spaces.push_back(&getField(type.name).get()->space()); addWeakFormArg(weak_form_name, type.name, arg_num); - if (type.is_unknown || type.name == field_names.unknown) { - addWeakFormUnknownArg(weak_form_name, type.name, arg_num); - } unknown_found = unknown_found || type.name == field_names.unknown; ++arg_num; }; @@ -325,7 +290,7 @@ struct FieldStore { * @param residual_names Names of the residuals. * @return Slot-list matrix mapping residual rows and solved columns. */ - BlockIndexMap indexMap(const std::vector& residual_names) const; + BlockArgumentMap indexMap(const std::vector& residual_names) const; /** * @brief Get the boundary condition managers for the given weak forms, one per residual row. @@ -475,15 +440,8 @@ struct FieldStore { /// Boundary conditions keyed by primary unknown field name. std::map> boundary_conditions_; - struct FieldLabel { - std::string field_name; - size_t field_index_in_residual; - }; - std::shared_ptr addBoundaryConditions(FEFieldPtr field); - std::map> weak_form_name_to_unknown_name_index_; - std::map> weak_form_name_to_field_names_; std::vector> weak_form_field_names_; diff --git a/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp index 9f6d1a8623..57f4d75f04 100644 --- a/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/multiphysics_time_integrator.cpp @@ -69,14 +69,14 @@ std::pair, std::vector> MultiphysicsTimeI SLIC_ERROR_ROOT_IF(cycle_zero_unknowns.size() != cz_sys->weak_forms.size(), "Cycle zero system result count does not match number of cycle-zero weak forms"); - SLIC_ERROR_ROOT_IF( - !cz_sys->solved_field_names.empty() && cz_sys->solved_field_names.size() != cz_sys->weak_forms.size(), - "Cycle zero solved_field_names size does not match number of weak forms"); + SLIC_ERROR_ROOT_IF(!cz_sys->solve_result_field_names.empty() && + cz_sys->solve_result_field_names.size() != cz_sys->weak_forms.size(), + "Cycle zero solve_result_field_names size does not match number of weak forms"); for (size_t i = 0; i < cz_sys->weak_forms.size(); ++i) { const std::string result_field_name = - cz_sys->solved_field_names.empty() + cz_sys->solve_result_field_names.empty() ? system_->field_store->getWeakFormFieldNames(cz_sys->weak_forms[i]->name()).unknown - : cz_sys->solved_field_names[i]; + : cz_sys->solve_result_field_names[i]; size_t result_field_state_idx = system_->field_store->getFieldIndex(result_field_name); current_states[result_field_state_idx] = cycle_zero_unknowns[i]; system_->field_store->setField(result_field_state_idx, cycle_zero_unknowns[i]); diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 7929a8b36d..698ee87d59 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -26,8 +26,9 @@ void applyBoundaryConditions(double time, const smith::BoundaryConditionManager* } } -std::vector block_solve(const std::vector& residual_evals, const BlockIndexMap block_indices, - const FieldState& shape_disp, const std::vector>& states, +std::vector block_solve(const std::vector& residual_evals, + const BlockArgumentMap& block_indices, const FieldState& shape_disp, + const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, const std::vector& bc_managers) diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp index c3d79fb0a4..33497b3c02 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.hpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -32,7 +32,7 @@ static constexpr size_t invalid_block_index = std::numeric_limits::max() using BlockArgumentIndices = std::vector; /// Matrix mapping residual rows and solved columns to local argument slots. -using BlockIndexMap = std::vector>; +using BlockArgumentMap = std::vector>; /// @brief Solve a block nonlinear system of equations as defined by the vector of weak form /// @param residual_evals Vector of weak forms which define the equations to be solved @@ -50,7 +50,8 @@ using BlockIndexMap = std::vector>; /// @param solver The nonlinear block solver used to solve the system of equations /// @param bc_managers Holds information about which degrees of freedom (DOFS) /// @return Vector of field solutions satisfying the weak forms -std::vector block_solve(const std::vector& residual_evals, const BlockIndexMap block_indices, +std::vector block_solve(const std::vector& residual_evals, + const BlockArgumentMap& block_indices, const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, @@ -72,7 +73,7 @@ inline FieldState solve(const WeakForm& residual_eval, const FieldState& shape_d const DirichletBoundaryConditions& bcs, size_t unknown_state_index = 0) { std::vector bc_managers{&bcs.getBoundaryConditionManager()}; - BlockIndexMap block_indices{{BlockArgumentIndices{unknown_state_index}}}; + BlockArgumentMap block_indices{{BlockArgumentIndices{unknown_state_index}}}; auto solutions = block_solve({const_cast(&residual_eval)}, block_indices, shape_disp, {states}, {params}, time_info, &solver, bc_managers); return solutions[0]; diff --git a/src/smith/differentiable_numerics/solid_mechanics_system.hpp b/src/smith/differentiable_numerics/solid_mechanics_system.hpp index 78cf957b91..9e9a603239 100644 --- a/src/smith/differentiable_numerics/solid_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/solid_mechanics_system.hpp @@ -246,7 +246,7 @@ auto registerSolidMechanicsFields(std::shared_ptr field_store, auto physics_fields = PhysicsFields, H1, H1, H1>{ - field_store, FieldType>(field_store->prefix("displacement_solve_state"), true), + field_store, FieldType>(field_store->prefix("displacement_solve_state")), FieldType>(field_store->prefix("displacement")), FieldType>(field_store->prefix("velocity")), FieldType>(field_store->prefix("acceleration"))}; @@ -307,7 +307,7 @@ auto buildSolidMechanicsSystemImpl(std::shared_ptr field_store, cons auto disp_time_rule_ptr = std::make_shared(); FieldType> shape_disp_type(field_store->prefix("shape_displacement")); - FieldType> disp_type(field_store->prefix("displacement_solve_state"), true); + FieldType> disp_type(field_store->prefix("displacement_solve_state")); FieldType> disp_old_type(field_store->prefix("displacement")); FieldType> velo_old_type(field_store->prefix("velocity")); FieldType> accel_old_type(field_store->prefix("acceleration")); @@ -335,7 +335,7 @@ auto buildSolidMechanicsSystemImpl(std::shared_ptr field_store, cons auto cycle_zero_solver = detail::makeCycleZeroSolver(solver, *field_store->getMesh()); auto cycle_zero_system = makeSystem(field_store, cycle_zero_solver, {sys->solid_weak_form}); - cycle_zero_system->solved_field_names = {accel_old_type.name}; + cycle_zero_system->solve_result_field_names = {accel_old_type.name}; auto cycle_zero_inputs = std::vector{disp_old_type.name, disp_old_type.name, velo_old_type.name, accel_old_type.name}; auto append_if_state = [&](const auto& field) { @@ -351,7 +351,7 @@ auto buildSolidMechanicsSystemImpl(std::shared_ptr field_store, cons } if (has_stress_output) { - FieldType> stress_type(field_store->prefix("stress"), true); + FieldType> stress_type(field_store->prefix("stress")); FieldType> disp_as_input(disp_type.name); std::string stress_name = field_store->prefix("stress_projection"); sys->stress_weak_form = detail::buildWeakFormWithCoupling( diff --git a/src/smith/differentiable_numerics/state_variable_system.hpp b/src/smith/differentiable_numerics/state_variable_system.hpp index 1d1be1e680..752dd1d768 100644 --- a/src/smith/differentiable_numerics/state_variable_system.hpp +++ b/src/smith/differentiable_numerics/state_variable_system.hpp @@ -130,7 +130,7 @@ auto registerInternalVariableFields(std::shared_ptr field_store, } return PhysicsFields{ - field_store, FieldType(field_store->prefix("state_solve_state"), true), + field_store, FieldType(field_store->prefix("state_solve_state")), FieldType(field_store->prefix("state"))}; } @@ -160,7 +160,7 @@ auto buildInternalVariableSystemImpl(std::shared_ptr field_store, co { auto internal_variable_time_rule = std::make_shared(); - FieldType state_type(field_store->prefix("state_solve_state"), true); + FieldType state_type(field_store->prefix("state_solve_state")); FieldType state_old_type(field_store->prefix("state")); auto internal_variable_bc = field_store->getBoundaryConditions(state_type.name); diff --git a/src/smith/differentiable_numerics/system_base.cpp b/src/smith/differentiable_numerics/system_base.cpp index ba1f5917e8..7908039ce2 100644 --- a/src/smith/differentiable_numerics/system_base.cpp +++ b/src/smith/differentiable_numerics/system_base.cpp @@ -18,8 +18,8 @@ std::vector SystemBase::solve(const TimeInfo& time_info) const for (const auto& wf : weak_forms) { weak_form_names.push_back(wf->name()); } - BlockIndexMap index_map; - if (solved_field_names.empty()) { + BlockArgumentMap index_map; + if (solve_result_field_names.empty()) { index_map = field_store->indexMap(weak_form_names); } @@ -45,20 +45,21 @@ std::vector SystemBase::solve(const TimeInfo& time_info) const } std::vector bc_field_names; - if (!solved_field_names.empty()) { - SLIC_ERROR_IF(solved_field_names.size() != weak_forms.size(), "solved_field_names size must match weak_forms size"); - SLIC_ERROR_IF(solve_input_field_names.empty(), "solve_input_field_names must accompany solved_field_names"); - bc_field_names = solved_field_names; - index_map = BlockIndexMap(weak_forms.size(), std::vector(weak_forms.size())); + if (!solve_result_field_names.empty()) { + SLIC_ERROR_IF(solve_result_field_names.size() != weak_forms.size(), + "solve_result_field_names size must match weak_forms size"); + SLIC_ERROR_IF(solve_input_field_names.empty(), "solve_input_field_names must accompany solve_result_field_names"); + bc_field_names = solve_result_field_names; + index_map = BlockArgumentMap(weak_forms.size(), std::vector(weak_forms.size())); std::set result_field_indices; - for (const auto& field_name : solved_field_names) { + for (const auto& field_name : solve_result_field_names) { size_t field_index = field_store->getFieldIndex(field_name); SLIC_ERROR_IF(!result_field_indices.insert(field_index).second, "Each solve result field must have exactly one owning row"); } for (size_t row = 0; row < weak_forms.size(); ++row) { for (size_t col = 0; col < weak_forms.size(); ++col) { - size_t result_field_index = field_store->getFieldIndex(solved_field_names[col]); + size_t result_field_index = field_store->getFieldIndex(solve_result_field_names[col]); for (size_t arg = 0; arg < input_field_indices[row].size(); ++arg) { if (input_field_indices[row][arg] == result_field_index) { index_map[row][col].push_back(arg); @@ -66,8 +67,9 @@ std::vector SystemBase::solve(const TimeInfo& time_info) const } } SLIC_ERROR_IF(index_map[row][row].empty(), "Requested solve result field '" - << solved_field_names[row] << "' is not an argument of weak form '" - << weak_form_names[row] << "'"); + << solve_result_field_names[row] + << "' is not an argument of weak form '" << weak_form_names[row] + << "'"); } } @@ -78,8 +80,9 @@ std::vector SystemBase::solve(const TimeInfo& time_info) const for (auto& p : weak_forms) { weak_form_ptrs.push_back(p.get()); } - auto bc_managers = solved_field_names.empty() ? field_store->getBoundaryConditionManagers(weak_form_names) - : field_store->getBoundaryConditionManagersForFields(bc_field_names); + auto bc_managers = solve_result_field_names.empty() + ? field_store->getBoundaryConditionManagers(weak_form_names) + : field_store->getBoundaryConditionManagersForFields(bc_field_names); return solver->solve(weak_form_ptrs, index_map, field_store->getShapeDisp(), inputs, wk_params, time_info, bc_managers); } diff --git a/src/smith/differentiable_numerics/system_base.hpp b/src/smith/differentiable_numerics/system_base.hpp index 395f26b0f6..c4a1c31de9 100644 --- a/src/smith/differentiable_numerics/system_base.hpp +++ b/src/smith/differentiable_numerics/system_base.hpp @@ -61,7 +61,7 @@ struct SystemBase { ///< independently; cycle-zero solves do not couple across subsystems. std::vector> post_solve_systems; ///< Optional systems solved after main state update. std::vector - solved_field_names; ///< Optional per-weak-form fields overriding registered solver-owned fields. + solve_result_field_names; ///< Optional per-weak-form fields overriding registered solver-owned fields. std::vector> solve_input_field_names; ///< Optional per-weak-form input field ordering used during solve. diff --git a/src/smith/differentiable_numerics/system_solver.cpp b/src/smith/differentiable_numerics/system_solver.cpp index 883cfd0370..d3878dd2e9 100644 --- a/src/smith/differentiable_numerics/system_solver.cpp +++ b/src/smith/differentiable_numerics/system_solver.cpp @@ -63,7 +63,7 @@ void SystemSolver::appendStagesWithBlockMapping(const SystemSolver& subsystem_so } std::vector SystemSolver::solve(const std::vector& residual_evals, - const BlockIndexMap& block_indices, const FieldState& shape_disp, + const BlockArgumentMap& block_indices, const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, @@ -164,7 +164,7 @@ std::vector SystemSolver::solve(const std::vector& residu size_t num_stage_blocks = stage.block_indices.size(); std::vector stage_residuals; - BlockIndexMap stage_block_indices; + BlockArgumentMap stage_block_indices; std::vector> stage_states; std::vector> stage_params; std::vector stage_bc_managers; diff --git a/src/smith/differentiable_numerics/system_solver.hpp b/src/smith/differentiable_numerics/system_solver.hpp index fabca974d2..1d8157391c 100644 --- a/src/smith/differentiable_numerics/system_solver.hpp +++ b/src/smith/differentiable_numerics/system_solver.hpp @@ -67,7 +67,7 @@ class SystemSolver { /// @param time_info Current time information. /// @param bc_managers Managers for boundary conditions. /// @return Updated field states. - std::vector solve(const std::vector& residual_evals, const BlockIndexMap& block_indices, + std::vector solve(const std::vector& residual_evals, const BlockArgumentMap& block_indices, const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const std::vector& bc_managers) const; diff --git a/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp b/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp index 0a2575753e..856d5e42c0 100644 --- a/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp +++ b/src/smith/differentiable_numerics/tests/test_combined_thermo_mechanics.cpp @@ -514,7 +514,7 @@ TEST_F(ThermoMechanicsMeshFixture, CombinedSystemCarriesCycleZeroSystems) ASSERT_EQ(solid_system->cycle_zero_systems.size(), 1u); EXPECT_EQ(combined_system->cycle_zero_systems.size(), solid_system->cycle_zero_systems.size()); - EXPECT_EQ(combined_system->cycle_zero_systems[0]->solved_field_names, std::vector{"acceleration"}); + EXPECT_EQ(combined_system->cycle_zero_systems[0]->solve_result_field_names, std::vector{"acceleration"}); } } // namespace smith diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index e97bf43027..28ca55567f 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -180,13 +180,14 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) std::vector reactions; for (size_t step = 0; step < 10; ++step) { - std::tie(states, reactions) = advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); + std::tie(states, reactions) = + advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); time += dt; pv_writer.write(step + 1, time, states); } - auto displacement = - plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); + auto displacement = plastic_mechanics_system->field_store->getField( + plastic_mechanics_system->field_store->prefix("displacement")); auto plastic_defgrad = plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_defgrad")); auto plastic_strain = @@ -346,8 +347,8 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); auto params = plastic_mechanics_system->field_store->getParameterFields(); - auto displacement = - plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); + auto displacement = plastic_mechanics_system->field_store->getField( + plastic_mechanics_system->field_store->prefix("displacement")); auto disp_squared = 0.5 * innerProduct(displacement, displacement); gretl::set_as_objective(disp_squared); diff --git a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp index fdfb0d838e..a750ef8c64 100644 --- a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp @@ -215,9 +215,6 @@ TEST(MultiphysicsTimeIntegrator, CycleZeroUsesOwnedFieldBoundaryConditions) FieldType> displacement_type("displacement"); auto displacement_bc = field_store->addIndependent(displacement_type, quasi_static); - ASSERT_TRUE(temperature_type.is_unknown); - ASSERT_TRUE(displacement_type.is_unknown); - temperature_bc->setScalarBCs<2>(mesh->domain("left"), [](double, tensor) { return 0.0; }); displacement_bc->setScalarBCs<2>(mesh->domain("right"), [](double, tensor) { return 1.0; }); @@ -355,7 +352,7 @@ TEST(MultiphysicsTimeIntegrator, CycleZeroSolveResultUpdatesAccelerationField) cycle_zero_system->field_store = field_store; cycle_zero_system->weak_forms = {cycle_zero_wf}; cycle_zero_system->solver = cycle_zero_solver; - cycle_zero_system->solved_field_names = {acceleration_type.name}; + cycle_zero_system->solve_result_field_names = {acceleration_type.name}; cycle_zero_system->solve_input_field_names = { {displacement_old_type.name, displacement_old_type.name, velocity_type.name, acceleration_type.name}}; diff --git a/src/smith/differentiable_numerics/thermal_system.hpp b/src/smith/differentiable_numerics/thermal_system.hpp index d0845857c9..f82ca2ef30 100644 --- a/src/smith/differentiable_numerics/thermal_system.hpp +++ b/src/smith/differentiable_numerics/thermal_system.hpp @@ -194,7 +194,7 @@ auto registerThermalFields(std::shared_ptr field_store, field_store->addDependent(temperature_type, FieldStore::TimeDerivative::VAL, "temperature"); return PhysicsFields, H1>{ - field_store, FieldType>(field_store->prefix("temperature_solve_state"), true), + field_store, FieldType>(field_store->prefix("temperature_solve_state")), FieldType>(field_store->prefix("temperature"))}; } @@ -214,7 +214,7 @@ auto buildThermalSystemImpl(std::shared_ptr field_store, const Coupl auto temperature_time_rule_ptr = std::make_shared(); FieldType> shape_disp_type(field_store->prefix("shape_displacement")); - FieldType> temperature_type(field_store->prefix("temperature_solve_state"), true); + FieldType> temperature_type(field_store->prefix("temperature_solve_state")); FieldType> temperature_old_type(field_store->prefix("temperature")); auto temperature_bc = field_store->getBoundaryConditions(temperature_type.name); From 6943e6344d4be9a77bd00b8d97fbe151799e525c Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Mon, 24 Aug 2026 17:10:21 -0700 Subject: [PATCH 11/28] Fix style. --- src/smith/differentiable_numerics/nonlinear_solve.cpp | 5 ++--- src/smith/differentiable_numerics/nonlinear_solve.hpp | 3 +-- .../tests/test_differentiable_plasticity.cpp | 11 +++++------ 3 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 698ee87d59..10f45c3fc9 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -26,9 +26,8 @@ void applyBoundaryConditions(double time, const smith::BoundaryConditionManager* } } -std::vector block_solve(const std::vector& residual_evals, - const BlockArgumentMap& block_indices, const FieldState& shape_disp, - const std::vector>& states, +std::vector block_solve(const std::vector& residual_evals, const BlockArgumentMap& block_indices, + const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, const std::vector& bc_managers) diff --git a/src/smith/differentiable_numerics/nonlinear_solve.hpp b/src/smith/differentiable_numerics/nonlinear_solve.hpp index 33497b3c02..1825bb7ab5 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.hpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.hpp @@ -50,8 +50,7 @@ using BlockArgumentMap = std::vector>; /// @param solver The nonlinear block solver used to solve the system of equations /// @param bc_managers Holds information about which degrees of freedom (DOFS) /// @return Vector of field solutions satisfying the weak forms -std::vector block_solve(const std::vector& residual_evals, - const BlockArgumentMap& block_indices, +std::vector block_solve(const std::vector& residual_evals, const BlockArgumentMap& block_indices, const FieldState& shape_disp, const std::vector>& states, const std::vector>& params, const TimeInfo& time_info, const NonlinearBlockSolverBase* solver, diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index 28ca55567f..e97bf43027 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -180,14 +180,13 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) std::vector reactions; for (size_t step = 0; step < 10; ++step) { - std::tie(states, reactions) = - advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); + std::tie(states, reactions) = advancer->advanceState(TimeInfo(time, dt, step), shape_disp, states, params); time += dt; pv_writer.write(step + 1, time, states); } - auto displacement = plastic_mechanics_system->field_store->getField( - plastic_mechanics_system->field_store->prefix("displacement")); + auto displacement = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); auto plastic_defgrad = plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("plastic_defgrad")); auto plastic_strain = @@ -347,8 +346,8 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) auto shape_disp = plastic_mechanics_system->field_store->getShapeDisp(); auto params = plastic_mechanics_system->field_store->getParameterFields(); - auto displacement = plastic_mechanics_system->field_store->getField( - plastic_mechanics_system->field_store->prefix("displacement")); + auto displacement = + plastic_mechanics_system->field_store->getField(plastic_mechanics_system->field_store->prefix("displacement")); auto disp_squared = 0.5 * innerProduct(displacement, displacement); gretl::set_as_objective(disp_squared); From 6aeb1486aeddfef15391342d8608bd3205c49047 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 25 Aug 2026 16:50:45 -0700 Subject: [PATCH 12/28] allow backward compatable interface for field store field initialization. --- .../differentiable_numerics/field_store.hpp | 14 +++++++++++++ .../plastic_mechanics_system.hpp | 18 ++++++++--------- .../differentiable_numerics/system_base.hpp | 8 ++++---- .../test_multiphysics_time_integrator.cpp | 20 +++++++++---------- 4 files changed, 35 insertions(+), 25 deletions(-) diff --git a/src/smith/differentiable_numerics/field_store.hpp b/src/smith/differentiable_numerics/field_store.hpp index a76adfb41b..9a22bd2299 100644 --- a/src/smith/differentiable_numerics/field_store.hpp +++ b/src/smith/differentiable_numerics/field_store.hpp @@ -262,6 +262,20 @@ struct FieldStore { return spaces; } + /** + * @brief Register input fields when test and unknown field names match. + * @param weak_form_name Name of the weak form being constructed. + * @param field_name Shared test-space and solver-owned field name. + * @param types Ordered list of @c FieldType descriptors for every input argument. + * @return std::vector Ordered input FE spaces. + */ + template + std::vector createSpaces(const std::string& weak_form_name, + const std::string& field_name, FieldTypes... types) + { + return createSpaces(weak_form_name, {.unknown = field_name, .test = field_name}, types...); + } + /** * @brief Mapping between primary and history/derivative fields for time integration. */ diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 042a646105..71d4ceced9 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -238,24 +238,22 @@ std::shared_ptr> bui std::string solid_res_name = field_store->prefix("solid_residual"); auto solid_weak_form = std::make_shared( solid_res_name, field_store->getMesh(), field_store->getField(disp_type.name).get()->space(), - field_store->createSpaces(solid_res_name, {.unknown = disp_type.name, .test = disp_type.name}, disp_type, - plastic_defgrad_type, plastic_strain_type, parameter_types...)); + field_store->createSpaces(solid_res_name, disp_type.name, disp_type, plastic_defgrad_type, plastic_strain_type, + parameter_types...)); std::string plastic_defgrad_res_name = field_store->prefix("plastic_defgrad_residual"); auto plastic_defgrad_weak_form = std::make_shared( plastic_defgrad_res_name, field_store->getMesh(), field_store->getField(plastic_defgrad_type.name).get()->space(), - field_store->createSpaces(plastic_defgrad_res_name, - {.unknown = plastic_defgrad_type.name, .test = plastic_defgrad_type.name}, - plastic_defgrad_type, plastic_defgrad_old_type, plastic_strain_type, - plastic_strain_old_type, disp_type, parameter_types...)); + field_store->createSpaces(plastic_defgrad_res_name, plastic_defgrad_type.name, plastic_defgrad_type, + plastic_defgrad_old_type, plastic_strain_type, plastic_strain_old_type, disp_type, + parameter_types...)); std::string plastic_strain_res_name = field_store->prefix("plastic_strain_residual"); auto plastic_strain_weak_form = std::make_shared( plastic_strain_res_name, field_store->getMesh(), field_store->getField(plastic_strain_type.name).get()->space(), - field_store->createSpaces(plastic_strain_res_name, - {.unknown = plastic_strain_type.name, .test = plastic_strain_type.name}, - plastic_strain_type, plastic_strain_old_type, plastic_defgrad_type, - plastic_defgrad_old_type, disp_type, parameter_types...)); + field_store->createSpaces(plastic_strain_res_name, plastic_strain_type.name, plastic_strain_type, + plastic_strain_old_type, plastic_defgrad_type, plastic_defgrad_old_type, disp_type, + parameter_types...)); // Build solver and advancer std::vector> weak_forms{solid_weak_form, plastic_defgrad_weak_form, diff --git a/src/smith/differentiable_numerics/system_base.hpp b/src/smith/differentiable_numerics/system_base.hpp index c4a1c31de9..fde91bd22b 100644 --- a/src/smith/differentiable_numerics/system_base.hpp +++ b/src/smith/differentiable_numerics/system_base.hpp @@ -121,10 +121,10 @@ auto buildWeakFormWithCouplingImpl(const FieldStorePtr& field_store, const std:: auto coupling_fields_tuple = std::get(args_tuple); return std::apply( [&](const auto&... coupling_fields) { - return std::make_shared( - weak_form_name, field_store->getMesh(), field_store->getField(unknown_field_name).get()->space(), - field_store->createSpaces(weak_form_name, {.unknown = unknown_field_name, .test = unknown_field_name}, - std::get(args_tuple)..., coupling_fields...)); + return std::make_shared(weak_form_name, field_store->getMesh(), + field_store->getField(unknown_field_name).get()->space(), + field_store->createSpaces(weak_form_name, unknown_field_name, + std::get(args_tuple)..., coupling_fields...)); }, coupling_fields_tuple); } diff --git a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp index a750ef8c64..d9bfd9121a 100644 --- a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp @@ -149,9 +149,8 @@ auto buildScalarDiffusionWeakForm(const std::string& name, std::shared_ptr FieldTypeT field_type) { using WeakFormType = FunctionalWeakForm<2, H1<1>, Parameters>>; - auto weak_form = std::make_shared( - name, mesh, fs->getField(field_type.name).get()->space(), - fs->createSpaces(name, {.unknown = field_type.name, .test = field_type.name}, field_type)); + auto weak_form = std::make_shared(name, mesh, fs->getField(field_type.name).get()->space(), + fs->createSpaces(name, field_type.name, field_type)); weak_form->addBodyIntegral(mesh->entireBodyName(), [](auto, auto, auto u) { return smith::tuple{0.0 * get(u), get(u)}; }); return weak_form; @@ -163,10 +162,10 @@ auto buildSecondOrderMainWeakForm(const std::string& name, std::shared_ptr VelocityFieldType velocity_type, AccelerationFieldType acceleration_type) { using WeakFormType = FunctionalWeakForm<2, H1<1>, Parameters, H1<1>, H1<1>, H1<1>>>; - auto weak_form = std::make_shared( - name, mesh, fs->getField(displacement_type.name).get()->space(), - fs->createSpaces(name, {.unknown = displacement_type.name, .test = displacement_type.name}, displacement_type, - displacement_old_type, velocity_type, acceleration_type)); + auto weak_form = + std::make_shared(name, mesh, fs->getField(displacement_type.name).get()->space(), + fs->createSpaces(name, displacement_type.name, displacement_type, + displacement_old_type, velocity_type, acceleration_type)); weak_form->addBodySource(mesh->entireBodyName(), [](auto, auto, auto...) { return 0.0; }); return weak_form; } @@ -179,8 +178,7 @@ auto buildSecondOrderCycleZeroWeakForm(const std::string& name, std::shared_ptr< using WeakFormType = FunctionalWeakForm<2, H1<1>, Parameters, H1<1>, H1<1>>>; auto weak_form = std::make_shared( name, mesh, fs->getField(acceleration_type.name).get()->space(), - fs->createSpaces(name, {.unknown = acceleration_type.name, .test = acceleration_type.name}, displacement_type, - velocity_type, acceleration_type)); + fs->createSpaces(name, acceleration_type.name, displacement_type, velocity_type, acceleration_type)); weak_form->addBodySource(mesh->entireBodyName(), [](auto, auto, auto...) { return 0.0; }); return weak_form; } @@ -554,8 +552,8 @@ TEST(SystemSolver, RoutesRepeatedDependenciesWithoutUpdatingHistory) using FieldBWeakForm = FunctionalWeakForm<2, H1<1>, Parameters, H1<1>, H1<1>, H1<1>>>; auto field_b_weak_form = std::make_shared( "field_b_residual", mesh, field_store->getField(field_b_type.name).get()->space(), - field_store->createSpaces("field_b_residual", {.unknown = field_b_type.name, .test = field_b_type.name}, - field_b_type, field_a_type, field_a_type, field_a_old_type)); + field_store->createSpaces("field_b_residual", field_b_type.name, field_b_type, field_a_type, field_a_type, + field_a_old_type)); field_b_weak_form->addBodySource(mesh->entireBodyName(), [](auto, auto, auto, auto field_a_first, auto field_a_second, auto field_a_old) { return field_a_first + field_a_second - 2.0 * field_a_old; From 2c8dd882c5a159603b31a139e67e5b7a68ea4368 Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Tue, 25 Aug 2026 18:07:48 -0700 Subject: [PATCH 13/28] Update some docs, named variables, simplify usage a bit. --- .../test_multiphysics_time_integrator.cpp | 42 +++++++++++++++---- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp index d9bfd9121a..2225ef95aa 100644 --- a/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp +++ b/src/smith/differentiable_numerics/tests/test_multiphysics_time_integrator.cpp @@ -531,25 +531,36 @@ TEST(SystemSolver, AppendsStagesWithBlockMappingForCombinedSubsystems) TEST(SystemSolver, RoutesRepeatedDependenciesWithoutUpdatingHistory) { + constexpr int dim = 2; + constexpr int order = 1; + using ShapeSpace = H1; + using FieldBTestSpace = H1; + using FieldBUnknownSpace = H1; + using FieldAFirstOccurrenceSpace = H1; + using FieldASecondOccurrenceSpace = H1; + using FieldAHistorySpace = H1; + using FieldBParameters = + Parameters; + using FieldBWeakForm = FunctionalWeakForm; + axom::sidre::DataStore datastore; StateManager::initialize(datastore, "system_solver_repeated_dependency_routing"); auto mesh = std::make_shared(mfem::Mesh::MakeCartesian2D(2, 2, mfem::Element::QUADRILATERAL, true, 1.0, 1.0), "system_solver_repeated_dependency_routing_mesh"); auto field_store = std::make_shared(mesh, 20); - FieldType> shape_disp_type("shape_displacement"); + FieldType shape_disp_type("shape_displacement"); field_store->addShapeDisp(shape_disp_type); auto static_rule = std::make_shared(); auto backward_euler_rule = std::make_shared(); - FieldType> field_a_type("field_a"); - FieldType> field_b_type("field_b"); + FieldType field_a_type("field_a"); + FieldType field_b_type("field_b"); field_store->addIndependent(field_a_type, backward_euler_rule); auto field_a_old_type = field_store->addDependent(field_a_type, FieldStore::TimeDerivative::VAL, "field_a_old"); field_store->addIndependent(field_b_type, static_rule); auto field_a_weak_form = buildScalarDiffusionWeakForm("field_a_residual", mesh, field_store, field_a_type); - using FieldBWeakForm = FunctionalWeakForm<2, H1<1>, Parameters, H1<1>, H1<1>, H1<1>>>; auto field_b_weak_form = std::make_shared( "field_b_residual", mesh, field_store->getField(field_b_type.name).get()->space(), field_store->createSpaces("field_b_residual", field_b_type.name, field_b_type, field_a_type, field_a_type, @@ -559,7 +570,7 @@ TEST(SystemSolver, RoutesRepeatedDependenciesWithoutUpdatingHistory) return field_a_first + field_a_second - 2.0 * field_a_old; }); - field_store->setField(field_store->getFieldIndex(field_a_old_type.name), field_store->getField(field_a_type.name)); + *field_store->getField(field_a_old_type.name).get() = *field_store->getField(field_a_type.name).get(); const std::vector residual_names{field_a_weak_form->name(), field_b_weak_form->name()}; const auto block_indices = field_store->indexMap(residual_names); @@ -568,10 +579,17 @@ TEST(SystemSolver, RoutesRepeatedDependenciesWithoutUpdatingHistory) EXPECT_EQ(block_indices[1][0], BlockArgumentIndices({1, 2})); EXPECT_EQ(block_indices[1][1], BlockArgumentIndices({0})); - auto field_a_solver = std::make_shared(3.0); + // Force field_a to a known nonzero solution without testing nonlinear convergence. Relaxation produces a + // predictable update, isolating whether SystemSolver propagates every repeated current-field slot while leaving + // history unchanged. + constexpr double field_a_solved_value = 3.0; + constexpr double field_a_relaxation_factor = 0.5; + constexpr int max_staggered_iterations = 1; + constexpr bool exact_staggered_steps = true; + auto field_a_solver = std::make_shared(field_a_solved_value); auto field_b_solver = std::make_shared(); - SystemSolver solver(1, true); - solver.addSubsystemSolver({0}, field_a_solver, 0.5); + SystemSolver solver(max_staggered_iterations, exact_staggered_steps); + solver.addSubsystemSolver({0}, field_a_solver, field_a_relaxation_factor); solver.addSubsystemSolver({1}, field_b_solver); const std::vector residuals{field_a_weak_form.get(), field_b_weak_form.get()}; @@ -583,8 +601,10 @@ TEST(SystemSolver, RoutesRepeatedDependenciesWithoutUpdatingHistory) auto solutions = solver.solve(residuals, block_indices, field_store->getShapeDisp(), states, params, time_info, bc_managers); - EXPECT_NEAR((*solutions[0].get())[0], 1.5, 1.0e-12); + // field_a starts at zero, ConstantNonlinearBlockSolver returns 3.0, and 0.5 relaxation produces 1.5. + EXPECT_NEAR((*solutions[0].get())[0], field_a_relaxation_factor * field_a_solved_value, 1.0e-12); + // Expected field_b inputs contain relaxed field_a in both current-value slots and unchanged field_a_old history. auto expected_field_b_states = states[1]; expected_field_b_states[1] = solutions[0]; expected_field_b_states[2] = solutions[0]; @@ -611,13 +631,17 @@ TEST(FieldStore, UsesExplicitUnknownFieldInsteadOfTestField) field_store->addIndependent(test_type, static_rule); auto unknown_bcs = field_store->addIndependent(unknown_type, static_rule); + // Most weak forms use one field for both solver-owned unknown and test space. This case separates those roles: + // unknown_field owns the solve and boundary conditions, while test_field defines residual and reaction spaces. field_store->createSpaces("residual", {.unknown = unknown_type.name, .test = test_type.name}, test_type, unknown_type); const auto block_indices = field_store->indexMap({"residual"}); EXPECT_EQ(field_store->getWeakFormFieldNames("residual").test, test_type.name); EXPECT_EQ(field_store->getWeakFormFieldNames("residual").unknown, unknown_type.name); + // Unlike the matching-name case, the diagonal unknown occupies argument slot 1 instead of test-field slot 0. EXPECT_EQ(block_indices[0][0], BlockArgumentIndices({1})); + // Boundary conditions follow the solver-owned unknown; reaction output follows the test space. EXPECT_EQ(field_store->getBoundaryConditionManagers({"residual"})[0], &unknown_bcs->getBoundaryConditionManager()); EXPECT_EQ(field_store->getReactionInfos()[0].space, &field_store->getField(test_type.name).get()->space()); From cb1dabea39aee3217c31e16e3463e358e402e4bf Mon Sep 17 00:00:00 2001 From: Michael Tupek Date: Wed, 26 Aug 2026 09:51:58 -0700 Subject: [PATCH 14/28] Lower the order for the plasticity solve so the gradient checks run faster enough in debug. --- .../tests/test_differentiable_plasticity.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index e97bf43027..9086b6da13 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -287,7 +287,7 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) int parallel_refinement = 0; static constexpr int dim = 3; - static constexpr int order = 2; + static constexpr int order = 1; axom::sidre::DataStore datastore; StateManager::initialize(datastore, "plasticity_small_strain"); From 25ce6bb442c31d172b66305e4755575c6074716e Mon Sep 17 00:00:00 2001 From: li97 Date: Fri, 4 Sep 2026 09:52:14 -0700 Subject: [PATCH 15/28] update finite difference test parameterization --- .../tests/test_differentiable_plasticity.cpp | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index 9086b6da13..cc52a25e13 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -215,8 +215,9 @@ class ParameterizedDifferentiableJ2SmallStrain { const T5& /* deltaHk */) const { constexpr auto I = Identity(); - auto K = (E0 + get(deltaE)) / (3.0 * (1.0 - 2.0 * nu)); - auto G = 0.5 * (E0 + get(deltaE)) / (1.0 + nu); + auto E = E0 * (1.0 + get(deltaE)); + auto K = E / (3.0 * (1.0 - 2.0 * nu)); + auto G = 0.5 * E / (1.0 + nu); auto el_strain = sym(du_dX) - Fp; auto p = K * tr(el_strain); @@ -231,8 +232,9 @@ class ParameterizedDifferentiableJ2SmallStrain { const tensor& du_dX, const T4& deltaE, const T5& deltaHk) const { using std::sqrt; - auto G = 0.5 * (E0 + get(deltaE)) / (1.0 + nu); - auto Hk = Hk0 + get(deltaHk); + auto E = E0 * (1.0 + get(deltaE)); + auto Hk = Hk0 * (1.0 + get(deltaHk)); + auto G = 0.5 * E / (1.0 + nu); auto el_strain = sym(du_dX) - Fp_old; auto s = 2.0 * G * dev(el_strain); @@ -252,8 +254,9 @@ class ParameterizedDifferentiableJ2SmallStrain { const T5& deltaE, const T6& deltaHk) const { using std::sqrt; - auto G = 0.5 * (E0 + get(deltaE)) / (1.0 + nu); - auto Hk = Hk0 + get(deltaHk); + auto E = E0 * (1.0 + get(deltaE)); + auto Hk = Hk0 * (1.0 + get(deltaHk)); + auto G = 0.5 * E / (1.0 + nu); auto el_strain = sym(du_dX) - Fp_old; auto s = 2.0 * G * dev(el_strain); @@ -309,7 +312,7 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) mesh, staggered_coupled_solver, physics_name, FieldType>("deltaE"), FieldType>("deltaHk")); using Hardening = solid_mechanics::LinearHardening; - Hardening hardening{.sigma_y = 25.0, .Hi = 50.0, .eta = 0.0}; + Hardening hardening{.sigma_y = 20.0, .Hi = 50.0, .eta = 0.0}; ParameterizedDifferentiableJ2SmallStrain mat(hardening, 1.5e+4, 0.25, 500.0, 1.0); plastic_mechanics_system->setMaterial(mesh->entireBodyName(), mat); @@ -352,8 +355,8 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) gretl::set_as_objective(disp_squared); EXPECT_GT(checkGradWrt(disp_squared, shape_disp, 2e-4, 4, true), 0.95); - EXPECT_GT(checkGradWrt(disp_squared, params[0], 100.0, 4, true), 0.95); - EXPECT_GT(checkGradWrt(disp_squared, params[1], 100.0, 4, true), 0.95); + EXPECT_GT(checkGradWrt(disp_squared, params[0], 1e-2, 4, true), 0.95); + EXPECT_GT(checkGradWrt(disp_squared, params[1], 1e-1, 4, true), 0.95); std::vector sensitivity_fields{shape_disp, params[0], params[1]}; auto sens_writer = createParaviewWriter(*mesh, sensitivity_fields, physics_name + "_final_sensitivities"); From b49926a5142c7da9f930c7f20356fb5673469281 Mon Sep 17 00:00:00 2001 From: Hanyu Li <125513291+lihanyu97@users.noreply.github.com> Date: Wed, 9 Sep 2026 13:19:38 -0500 Subject: [PATCH 16/28] Update .github/workflows/ci-tests.yml Co-authored-by: Chris White --- .github/workflows/ci-tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci-tests.yml b/.github/workflows/ci-tests.yml index 4f98dd03ed..a5b1d5fb0d 100644 --- a/.github/workflows/ci-tests.yml +++ b/.github/workflows/ci-tests.yml @@ -79,7 +79,7 @@ jobs: echo "compiler_image ${{ matrix.config.compiler_image }}" echo "host_config ${{ matrix.config.host_config }}" - name: Build and Test ${{ matrix.config.host_config }} ${{ matrix.build_type }} - timeout-minutes: 80 + timeout-minutes: 60 run: | ./scripts/llnl/build_src.py ${{ matrix.config.build_src_opts }} \ --verbose \ From 14dfaf2d84afab2372f25e5a90ed5dd8ecd7713b Mon Sep 17 00:00:00 2001 From: Hanyu Li <125513291+lihanyu97@users.noreply.github.com> Date: Wed, 9 Sep 2026 15:34:21 -0500 Subject: [PATCH 17/28] pass by reference for recoverTensor Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../differentiable_numerics/plastic_mechanics_system.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 71d4ceced9..d972325e88 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -79,17 +79,17 @@ struct PlasticMechanicsSystem : public SystemBase { * @brief Transform tensor variable stored as [dim * dim, 1] to [dim, dim] */ template - SMITH_HOST_DEVICE tensor recoverTensor(const tensor F_state) + SMITH_HOST_DEVICE tensor recoverTensor(const tensor& f_state) const { - tensor F_state_tensor{}; + tensor f_state_tensor{}; for (int i = 0; i < dim; ++i) { for (int j = 0; j < dim; ++j) { - F_state_tensor(i, j) = F_state(i * dim + j); + f_state_tensor(i, j) = f_state(i * dim + j); } } - return F_state_tensor; + return f_state_tensor; } /** From 86a3789fb7ea358ecab5498ea868e8751b950de2 Mon Sep 17 00:00:00 2001 From: Hanyu Li <125513291+lihanyu97@users.noreply.github.com> Date: Wed, 9 Sep 2026 15:34:49 -0500 Subject: [PATCH 18/28] pass by reference for flettenTensor Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../differentiable_numerics/plastic_mechanics_system.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index d972325e88..133bab6e40 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -96,17 +96,17 @@ struct PlasticMechanicsSystem : public SystemBase { * @brief Transform tensor variable stored as [dim, dim] to [dim * dim, 1] */ template - SMITH_HOST_DEVICE tensor flattenTensor(const tensor F_state) + SMITH_HOST_DEVICE tensor flattenTensor(const tensor& f_state) const { - tensor F_state_flattened{}; + tensor f_state_flattened{}; for (int i = 0; i < dim; ++i) { for (int j = 0; j < dim; ++j) { - F_state_flattened(i * dim + j) = F_state(i, j); + f_state_flattened(i * dim + j) = f_state(i, j); } } - return F_state_flattened; + return f_state_flattened; } /** From 0e5b2bb5cbb23f10faf2042cb9b79bc20db79484 Mon Sep 17 00:00:00 2001 From: Hanyu Li <125513291+lihanyu97@users.noreply.github.com> Date: Thu, 10 Sep 2026 11:53:09 -0500 Subject: [PATCH 19/28] fix comment Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/smith/differentiable_numerics/plastic_mechanics_system.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 133bab6e40..3d804cfaec 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -6,7 +6,7 @@ /** * @file plastic_mechanics_system.hpp - * @brief Defines the PlasticMechanicsystem struct and its factory function + * @brief Defines the PlasticMechanicsSystem struct and its factory function. */ #pragma once From dd3cc13c3a9b12622452ebcbec18cbef4a1d1648 Mon Sep 17 00:00:00 2001 From: Hanyu Li <125513291+lihanyu97@users.noreply.github.com> Date: Thu, 10 Sep 2026 11:53:29 -0500 Subject: [PATCH 20/28] fix comment Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/smith/differentiable_numerics/plastic_mechanics_system.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 3d804cfaec..9908af4f17 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -27,7 +27,7 @@ namespace smith { * @tparam disp_order Order of the displacement basis. * @tparam parameter_space Finite element spaces for optional parameters. * - * This system rely on a fixed-point iteration to solve plasticity equivalent to the traditional + * This system relies on a fixed-point iteration to solve plasticity, equivalent to a traditional * inner-outer Newton method. * n in [0, N] denotes the timestep and k in [0, L] denotes the fixed-point iteration count. * At every new timestep n+1, first predict new displacement by [u^{n+1,k}, Fp^{n+1,k}, epsilon_p^{n+1,k}] ---> From c096c693740672332ed8d58735e7b32e8738cdc9 Mon Sep 17 00:00:00 2001 From: Hanyu Li <125513291+lihanyu97@users.noreply.github.com> Date: Thu, 10 Sep 2026 11:53:52 -0500 Subject: [PATCH 21/28] fix comment Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/smith/differentiable_numerics/plastic_mechanics_system.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 9908af4f17..2b78410893 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -192,7 +192,7 @@ struct PlasticMechanicsSystem : public SystemBase { }; /** - * @brief Factory function to build a chemomechanics system with L2 state variables. + * @brief Factory function to build a plastic mechanics system with L2 internal variables. */ template std::shared_ptr> buildPlasticMechanicsSystem( From c92fc252231369363b98826658a590da0a2ac54e Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 22 Sep 2026 09:53:23 -0700 Subject: [PATCH 22/28] remove unused epsilon_current --- src/smith/differentiable_numerics/plastic_mechanics_system.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 2b78410893..58f4e14727 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -167,7 +167,7 @@ struct PlasticMechanicsSystem : public SystemBase { auto u, auto... params) { auto du_dX = get(u); auto Fp_old_tensor = this->recoverTensor(get(Fp_old)); - auto [epsilon_current, epsilon_dot] = captured_strain_rule->interpolate(t_info, epsilon_p, epsilon_p_old); + auto epsilon_dot = captured_strain_rule->dot(t_info, epsilon_p, epsilon_p_old); auto dt = t_info.dt(); auto Fp_predict_tensor = mat.plasticDeformGrad(dt, Fp_old_tensor, get(epsilon_dot), du_dX, params...); From 47c5ab78e15d161d9ae8cf329bffcc112b70a249 Mon Sep 17 00:00:00 2001 From: li97 Date: Tue, 6 Oct 2026 16:34:23 -0700 Subject: [PATCH 23/28] remove test dependency --- .../tests/test_differentiable_plasticity.cpp | 58 +++++++++++-------- 1 file changed, 33 insertions(+), 25 deletions(-) diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index cc52a25e13..99b4c4b8f0 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -12,31 +12,18 @@ namespace smith { -LinearSolverOptions primal_lin_opts{.linear_solver = LinearSolver::CG, - .preconditioner = Preconditioner::HypreAMG, - .relative_tol = 1e-15, - .absolute_tol = 1e-15, - .max_iterations = 200, - .print_level = 0}; - -NonlinearSolverOptions primal_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, - .relative_tol = 1e-9, - .absolute_tol = 1e-10, - .max_iterations = 25, - .print_level = 2}; - -LinearSolverOptions state_lin_opts{.linear_solver = LinearSolver::CG, - .preconditioner = Preconditioner::HypreJacobi, - .relative_tol = 1e-10, +LinearSolverOptions lin_opts{.linear_solver = LinearSolver::CG, + .preconditioner = Preconditioner::HypreAMG, + .relative_tol = 1e-15, + .absolute_tol = 1e-15, + .max_iterations = 200, + .print_level = 0}; + +NonlinearSolverOptions nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, + .relative_tol = 1e-9, .absolute_tol = 1e-10, - .max_iterations = 200, - .print_level = 0}; - -NonlinearSolverOptions state_nonlin_opts{.nonlin_solver = NonlinearSolver::Newton, - .relative_tol = 1e-7, - .absolute_tol = 1e-8, - .max_iterations = 25, - .print_level = 2}; + .max_iterations = 25, + .print_level = 2}; /// @brief Differentiable J2 material with nonlinear isotropic hardening and linear kinematic hardening template @@ -115,6 +102,17 @@ TEST(DifferentiablePlasticity, J2SmallStrainLinearHardening) { MPI_Barrier(MPI_COMM_WORLD); + NonlinearSolverOptions primal_nonlin_opts = nonlin_opts; + LinearSolverOptions primal_lin_opts = lin_opts; + + NonlinearSolverOptions state_nonlin_opts = nonlin_opts; + state_nonlin_opts.relative_tol = 1e-7; + state_nonlin_opts.absolute_tol = 1e-8; + LinearSolverOptions state_lin_opts = lin_opts; + state_lin_opts.preconditioner = Preconditioner::HypreJacobi; + state_lin_opts.relative_tol = 1e-10; + state_lin_opts.absolute_tol = 1e-10; + int serial_refinement = 0; int parallel_refinement = 0; @@ -283,8 +281,18 @@ TEST(DifferentiablePlasticity, PlasticLoadingFinitDiff) std::string physics_name = "differentiable_plasticity"; - state_nonlin_opts.print_level = 0; + NonlinearSolverOptions primal_nonlin_opts = nonlin_opts; primal_nonlin_opts.print_level = 0; + LinearSolverOptions primal_lin_opts = lin_opts; + + NonlinearSolverOptions state_nonlin_opts = nonlin_opts; + state_nonlin_opts.relative_tol = 1e-7; + state_nonlin_opts.absolute_tol = 1e-8; + state_nonlin_opts.print_level = 0; + LinearSolverOptions state_lin_opts = lin_opts; + state_lin_opts.preconditioner = Preconditioner::HypreJacobi; + state_lin_opts.relative_tol = 1e-10; + state_lin_opts.absolute_tol = 1e-10; int serial_refinement = 0; int parallel_refinement = 0; From 0b403cad167e1c08f6fdc5aa9e67d9fd21adb62f Mon Sep 17 00:00:00 2001 From: li97 Date: Fri, 16 Oct 2026 15:32:31 -0700 Subject: [PATCH 24/28] add release note --- RELEASE-NOTES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 950392b794..4feef4fe35 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -31,6 +31,7 @@ The Smith project release numbers follow [Semantic Versioning](http://semver.org - Added axisymmetric solid mechanics materials and loads for 2D `(r, z)` meshes. - Added host-configs for CZMatrix and RZVector machines - Added CI pipelines for CZMatrix machine +- Added differentiable plastic system with internal variables represented by discontinuous galerkin nodal fields ### Removed From a6c20d176421626c03e6407510a4e728dfd4aeb7 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 21 Oct 2026 10:43:53 -0700 Subject: [PATCH 25/28] switch to std::format --- .../differentiable_physics.cpp | 28 +++++++--------- .../differentiable_numerics/field_store.cpp | 32 +++++++++---------- .../differentiable_numerics/field_store.hpp | 6 ++-- .../nonlinear_solve.cpp | 8 ++--- .../plastic_mechanics_system.hpp | 4 +-- .../differentiable_numerics/system_solver.cpp | 19 ++++++----- .../tests/test_solid_dynamics.cpp | 27 ++++++++-------- src/smith/numerics/block_preconditioner.cpp | 21 ++++++------ src/smith/numerics/equation_solver.cpp | 6 ++-- src/smith/numerics/nonlinear_convergence.cpp | 1 - src/smith/physics/functional_weak_form.hpp | 29 ++++++++--------- 11 files changed, 85 insertions(+), 96 deletions(-) diff --git a/src/smith/differentiable_numerics/differentiable_physics.cpp b/src/smith/differentiable_numerics/differentiable_physics.cpp index 4928db1a2b..58353be536 100644 --- a/src/smith/differentiable_numerics/differentiable_physics.cpp +++ b/src/smith/differentiable_numerics/differentiable_physics.cpp @@ -67,9 +67,8 @@ DifferentiablePhysics::DifferentiablePhysics(std::shared_ptr mesh, std::sh reaction_names_.reserve(reaction_infos_.size()); for (size_t i = 0; i < reaction_infos_.size(); ++i) { - SLIC_ERROR_IF( - reaction_infos_[i].space == nullptr, - axom::fmt::format("Dual '{}' in physics module '{}' has null FE space.", reaction_infos_[i].name, name_)); + SLIC_ERROR_IF(reaction_infos_[i].space == nullptr, + std::format("Dual '{}' in physics module '{}' has null FE space.", reaction_infos_[i].name, name_)); reaction_names_.push_back(reaction_infos_[i].name); reaction_name_to_reaction_index_[reaction_infos_[i].name] = i; } @@ -135,7 +134,7 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& reaction if (reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end()) { std::string available; for (auto& n : reaction_names_) available += n + " "; - SLIC_ERROR(axom::fmt::format( + SLIC_ERROR(std::format( "Could not find reaction named {0} in mesh with tag \"{1}\" to get. Available reactions (size {2}): {3}", reaction_name, mesh_->tag(), reaction_names_.size(), available)); } @@ -153,11 +152,10 @@ const FiniteElementDual& DifferentiablePhysics::dual(const std::string& reaction FiniteElementDual DifferentiablePhysics::loadCheckpointedDual(const std::string& reaction_name, int cycle) { - SLIC_ERROR_IF( - cycle != cycle_, - axom::fmt::format("Due to checkpointing restrictions in smith::DifferentiablePhysics, cannot ask for " - "an arbitrary checkpointed reaction cycle, asking for cycle {}, but physics is at cycle {}", - cycle, cycle_)); + SLIC_ERROR_IF(cycle != cycle_, + std::format("Due to checkpointing restrictions in smith::DifferentiablePhysics, cannot ask for " + "an arbitrary checkpointed reaction cycle, asking for cycle {}, but physics is at cycle {}", + cycle, cycle_)); return dual(reaction_name); } @@ -218,9 +216,8 @@ void DifferentiablePhysics::setAdjointLoad( for (auto string_reaction_pair : string_to_reaction) { std::string field_name = string_reaction_pair.first; const smith::FiniteElementDual& reaction = string_reaction_pair.second; - SLIC_ERROR_IF( - state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(), - axom::fmt::format("Could not find reaction named {0} in mesh with tag {1}", field_name, mesh_->tag())); + SLIC_ERROR_IF(state_name_to_field_index_.find(field_name) == state_name_to_field_index_.end(), + std::format("Could not find reaction named {0} in mesh with tag {1}", field_name, mesh_->tag())); size_t state_index = state_name_to_field_index_.at(field_name); *field_states_[state_index].get_dual() += reaction; } @@ -232,10 +229,9 @@ void DifferentiablePhysics::setDualAdjointBcs( for (auto string_bc_pair : string_to_bc) { std::string reaction_name = string_bc_pair.first; const smith::FiniteElementState& reaction_adjoint_load = string_bc_pair.second; - SLIC_ERROR_IF( - reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end(), - axom::fmt::format("When calling setDualAdjointBcs, could not find reaction named {0} in mesh with tag {1}", - reaction_name, mesh_->tag())); + SLIC_ERROR_IF(reaction_name_to_reaction_index_.find(reaction_name) == reaction_name_to_reaction_index_.end(), + std::format("When calling setDualAdjointBcs, could not find reaction named {0} in mesh with tag {1}", + reaction_name, mesh_->tag())); size_t reaction_index = reaction_name_to_reaction_index_.at(reaction_name); *reaction_states_[reaction_index].get_dual() += reaction_adjoint_load; } diff --git a/src/smith/differentiable_numerics/field_store.cpp b/src/smith/differentiable_numerics/field_store.cpp index a975cce84c..04ff05cf4a 100644 --- a/src/smith/differentiable_numerics/field_store.cpp +++ b/src/smith/differentiable_numerics/field_store.cpp @@ -61,8 +61,8 @@ BlockArgumentMap FieldStore::indexMap(const std::vector& residual_n const std::string& res_name = residual_names[res_i]; size_t global_state_idx = to_states_index_.at(getWeakFormFieldNames(res_name).unknown); auto [owner, inserted] = global_state_to_local_col.emplace(global_state_idx, res_i); - SLIC_ERROR_IF(!inserted, axom::fmt::format("Weak forms '{}' and '{}' both own state index {}", - residual_names[owner->second], res_name, global_state_idx)); + SLIC_ERROR_IF(!inserted, std::format("Weak forms '{}' and '{}' both own state index {}", + residual_names[owner->second], res_name, global_state_idx)); } BlockArgumentMap block_indices(residual_names.size(), @@ -239,9 +239,9 @@ const std::vector& FieldStore::getAllFields() const { return states_ std::vector FieldStore::getStates(const std::string& weak_form_name) const { // Validate that weak form is registered - SLIC_ERROR_ROOT_IF(weak_form_name_to_field_names_.count(weak_form_name) == 0, - axom::fmt::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", - weak_form_name)); + SLIC_ERROR_ROOT_IF( + weak_form_name_to_field_names_.count(weak_form_name) == 0, + std::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", weak_form_name)); auto field_names = weak_form_name_to_field_names_.at(weak_form_name); std::vector fields_for_residual; @@ -249,7 +249,7 @@ std::vector FieldStore::getStates(const std::string& weak_form_name) // Validate that field exists SLIC_ERROR_ROOT_IF( to_states_index_.count(name) == 0 && to_params_index_.count(name) == 0, - axom::fmt::format("Field '{}' (required by weak form '{}') not found in FieldStore", name, weak_form_name)); + std::format("Field '{}' (required by weak form '{}') not found in FieldStore", name, weak_form_name)); // Only include state fields, not parameters // Parameters are passed separately to avoid duplication in block_solve @@ -265,9 +265,9 @@ std::vector FieldStore::getStatesFromVectors(const std::string& weak const std::vector& param_fields) const { // Validate that weak form is registered - SLIC_ERROR_ROOT_IF(weak_form_name_to_field_names_.count(weak_form_name) == 0, - axom::fmt::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", - weak_form_name)); + SLIC_ERROR_ROOT_IF( + weak_form_name_to_field_names_.count(weak_form_name) == 0, + std::format("Weak form '{}' not found in FieldStore. Did you forget to call addReaction()?", weak_form_name)); auto field_names = weak_form_name_to_field_names_.at(weak_form_name); std::vector fields_for_residual; @@ -275,21 +275,21 @@ std::vector FieldStore::getStatesFromVectors(const std::string& weak // Check if it's a state field if (to_states_index_.count(name)) { size_t idx = to_states_index_.at(name); - SLIC_ERROR_ROOT_IF(idx >= state_fields.size(), - axom::fmt::format("State field index {} out of bounds (size={}) for field '{}'", idx, - state_fields.size(), name)); + SLIC_ERROR_ROOT_IF( + idx >= state_fields.size(), + std::format("State field index {} out of bounds (size={}) for field '{}'", idx, state_fields.size(), name)); fields_for_residual.push_back(state_fields[idx]); } // Otherwise check if it's a parameter else if (to_params_index_.count(name)) { size_t idx = to_params_index_.at(name); SLIC_ERROR_ROOT_IF(idx >= param_fields.size(), - axom::fmt::format("Parameter field index {} out of bounds (size={}) for field '{}'", idx, - param_fields.size(), name)); + std::format("Parameter field index {} out of bounds (size={}) for field '{}'", idx, + param_fields.size(), name)); fields_for_residual.push_back(param_fields[idx]); } else { - SLIC_ERROR_ROOT(axom::fmt::format("Field or parameter '{}' (required by weak form '{}') not found in FieldStore", - name, weak_form_name)); + SLIC_ERROR_ROOT(std::format("Field or parameter '{}' (required by weak form '{}') not found in FieldStore", name, + weak_form_name)); } } return fields_for_residual; diff --git a/src/smith/differentiable_numerics/field_store.hpp b/src/smith/differentiable_numerics/field_store.hpp index 9a22bd2299..8dcd9edada 100644 --- a/src/smith/differentiable_numerics/field_store.hpp +++ b/src/smith/differentiable_numerics/field_store.hpp @@ -109,9 +109,9 @@ struct FieldStore { // Verify the space matches by checking vdim (== Space::components). auto& existing = params_[to_params_index_.at(type.name)]; SLIC_ERROR_ROOT_IF(existing.get()->space().GetVDim() != Space::components, - axom::fmt::format("Parameter '{}' re-registered with a different space " - "(existing vdim={}, new vdim={})", - type.name, existing.get()->space().GetVDim(), Space::components)); + std::format("Parameter '{}' re-registered with a different space " + "(existing vdim={}, new vdim={})", + type.name, existing.get()->space().GetVDim(), Space::components)); return; } to_params_index_[type.name] = params_.size(); diff --git a/src/smith/differentiable_numerics/nonlinear_solve.cpp b/src/smith/differentiable_numerics/nonlinear_solve.cpp index 10f45c3fc9..42244f7d20 100644 --- a/src/smith/differentiable_numerics/nonlinear_solve.cpp +++ b/src/smith/differentiable_numerics/nonlinear_solve.cpp @@ -53,16 +53,16 @@ std::vector block_solve(const std::vector& residual_evals for (size_t col = 0; col < num_rows_; ++col) { for (size_t idx : block_indices[row][col]) { SLIC_ERROR_IF(idx >= states[row].size(), - axom::fmt::format("block_indices[{}][{}] contains {} outside states[{}] size {}", row, col, idx, - row, states[row].size())); + std::format("block_indices[{}][{}] contains {} outside states[{}] size {}", row, col, idx, row, + states[row].size())); SLIC_ERROR_IF( slot_owners[idx] != invalid_block_index, - axom::fmt::format("State slot {} in row {} belongs to columns {} and {}", idx, row, slot_owners[idx], col)); + std::format("State slot {} in row {} belongs to columns {} and {}", idx, row, slot_owners[idx], col)); slot_owners[idx] = col; } } SLIC_ERROR_IF(block_indices[row][row].empty(), - axom::fmt::format("block_indices[{}][{}] diagonal entry must not be empty", row, row)); + std::format("block_indices[{}][{}] diagonal entry must not be empty", row, row)); } std::vector num_state_inputs; diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 58f4e14727..808a14713d 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -58,11 +58,11 @@ struct PlasticMechanicsSystem : public SystemBase { L2, H1, parameter_space...>>; // Primary weak forms - std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. + std::shared_ptr solid_weak_form; // Internal variable weak forms to update in staggered solve std::shared_ptr plastic_deform_weak_form; ///< Plastic deformation weak form. - std::shared_ptr plastic_strain_weak_form; ///< Plastic deformation weak form. + std::shared_ptr plastic_strain_weak_form; ///< Plastic strain weak form. // Primary variable bcs std::shared_ptr disp_bc; ///< Displacement boundary conditions. diff --git a/src/smith/differentiable_numerics/system_solver.cpp b/src/smith/differentiable_numerics/system_solver.cpp index d3878dd2e9..b18345181a 100644 --- a/src/smith/differentiable_numerics/system_solver.cpp +++ b/src/smith/differentiable_numerics/system_solver.cpp @@ -14,7 +14,6 @@ #include #include #include -#include namespace smith { @@ -35,7 +34,7 @@ void SystemSolver::addSubsystemSolver(const std::vector& block_indices, { SLIC_ERROR_IF(!solver, "SystemSolver stage solver must be non-null"); SLIC_ERROR_IF(relaxation_factor <= 0.0 || relaxation_factor > 1.0, - axom::fmt::format("Stage relaxation_factor {} must be in (0, 1]", relaxation_factor)); + std::format("Stage relaxation_factor {} must be in (0, 1]", relaxation_factor)); stages_.push_back(Stage{block_indices, std::move(solver), relaxation_factor}); } @@ -53,8 +52,8 @@ void SystemSolver::appendStagesWithBlockMapping(const SystemSolver& subsystem_so remapped_block_indices.reserve(stage.block_indices.size()); for (size_t local_block_index : stage.block_indices) { SLIC_ERROR_IF(local_block_index >= global_block_indices.size(), - axom::fmt::format("Local block index {} exceeds subsystem size {}", local_block_index, - global_block_indices.size())); + std::format("Local block index {} exceeds subsystem size {}", local_block_index, + global_block_indices.size())); remapped_block_indices.push_back(global_block_indices[local_block_index]); } } @@ -76,16 +75,16 @@ std::vector SystemSolver::solve(const std::vector& residu SLIC_ERROR_IF(block_indices.size() != num_residuals, "Block index rows must match residual count"); for (size_t row = 0; row < num_residuals; ++row) { SLIC_ERROR_IF(block_indices[row].size() != num_residuals, - axom::fmt::format("Block index columns for row {} must match residual count", row)); + std::format("Block index columns for row {} must match residual count", row)); for (size_t col = 0; col < num_residuals; ++col) { for (size_t slot : block_indices[row][col]) { SLIC_ERROR_IF(slot >= states[row].size(), - axom::fmt::format("Block slot {} for row {}, column {} exceeds state count {}", slot, row, col, - states[row].size())); + std::format("Block slot {} for row {}, column {} exceeds state count {}", slot, row, col, + states[row].size())); } } SLIC_ERROR_IF(block_indices[row][row].empty(), - axom::fmt::format("Block index diagonal for row {} must not be empty", row)); + std::format("Block index diagonal for row {} must not be empty", row)); } std::vector active_stages = stages_; for (auto& stage : active_stages) { @@ -95,7 +94,7 @@ std::vector SystemSolver::solve(const std::vector& residu } for (size_t block_index : stage.block_indices) { SLIC_ERROR_IF(block_index >= num_residuals, - axom::fmt::format("Stage block index {} exceeds residual count {}", block_index, num_residuals)); + std::format("Stage block index {} exceeds residual count {}", block_index, num_residuals)); } } // Set the inner tolerance factor based on the number of stages. For single-stage @@ -226,7 +225,7 @@ std::vector SystemSolver::solve(const std::vector& residu } } if (all_converged) { - SLIC_INFO_ROOT(axom::fmt::format("Staggered iteration converged after {} iteration(s)", iter + 1)); + SLIC_INFO_ROOT(std::format("Staggered iteration converged after {} iteration(s)", iter + 1)); break; } } diff --git a/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp b/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp index 333e8d39c0..4f0bf3ca80 100644 --- a/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp +++ b/src/smith/differentiable_numerics/tests/test_solid_dynamics.cpp @@ -418,15 +418,15 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesGretl) reaction_squared.data_store().back_prop(); for (auto s : initial_states) { - SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", s.get()->name(), s.get()->Norml2(), s.get_dual()->Norml2())); + SLIC_INFO_ROOT(std::format("{} {} {}", s.get()->name(), s.get()->Norml2(), s.get_dual()->Norml2())); } - SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", shape_disp.get()->name(), shape_disp.get()->Norml2(), - shape_disp.get_dual()->Norml2())); + SLIC_INFO_ROOT( + std::format("{} {} {}", shape_disp.get()->name(), shape_disp.get()->Norml2(), shape_disp.get_dual()->Norml2())); for (size_t p = 0; p < params.size(); ++p) { - SLIC_INFO_ROOT(axom::fmt::format("{} {} {}", params[p].get()->name(), params[p].get()->Norml2(), - params[p].get_dual()->Norml2())); + SLIC_INFO_ROOT( + std::format("{} {} {}", params[p].get()->name(), params[p].get()->Norml2(), params[p].get_dual()->Norml2())); } } @@ -473,7 +473,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) createSolidMechanicsBasePhysics(physics_name, mesh, solid_nonlinear_opts, solid_linear_options); double qoi = integrateForward(*physics, num_steps_, dt_, physics_name + "_reactions"); - SLIC_INFO_ROOT(axom::fmt::format("{}", qoi)); + SLIC_INFO_ROOT(std::format("{}", qoi)); // Check that reaction forces are zero away from Dirichlet DOFs auto reactions = physics->getReactionStates(); @@ -491,14 +491,13 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesBasePhysics) auto state_sensitivities = physics->computeInitialConditionSensitivity(); for (auto name_and_state_sensitivity : state_sensitivities) { - SLIC_INFO_ROOT( - axom::fmt::format("{} {}", name_and_state_sensitivity.first, name_and_state_sensitivity.second.Norml2())); + SLIC_INFO_ROOT(std::format("{} {}", name_and_state_sensitivity.first, name_and_state_sensitivity.second.Norml2())); } - SLIC_INFO_ROOT(axom::fmt::format("{} {}", shape_sensitivity.name(), shape_sensitivity.Norml2())); + SLIC_INFO_ROOT(std::format("{} {}", shape_sensitivity.name(), shape_sensitivity.Norml2())); for (size_t p = 0; p < num_params; ++p) { - SLIC_INFO_ROOT(axom::fmt::format("{} {}", parameter_sensitivities[p].name(), parameter_sensitivities[p].Norml2())); + SLIC_INFO_ROOT(std::format("{} {}", parameter_sensitivities[p].name(), parameter_sensitivities[p].Norml2())); } } @@ -557,13 +556,13 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesComparison) // Compare shape sensitivities double shape_diff = diff_norm(*shape_dispG.get_dual(), shape_sensitivityB); - SLIC_INFO_ROOT(axom::fmt::format("Shape sensitivity difference: {:.6e}", shape_diff)); + SLIC_INFO_ROOT(std::format("Shape sensitivity difference: {:.6e}", shape_diff)); EXPECT_LT(shape_diff, tol); // Compare parameter sensitivities for (size_t p = 0; p < num_params; ++p) { double param_diff = diff_norm(*paramsG[p].get_dual(), parameter_sensitivitiesB[p]); - SLIC_INFO_ROOT(axom::fmt::format("Parameter {} sensitivity difference: {:.6e}", p, param_diff)); + SLIC_INFO_ROOT(std::format("Parameter {} sensitivity difference: {:.6e}", p, param_diff)); EXPECT_LT(param_diff, tol); } @@ -572,7 +571,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesComparison) for (const auto& suffix : state_suffixes) { std::string nameG = physics_name + "_gretl_" + suffix; std::string nameB = physics_name + "_base_" + suffix; - SLIC_INFO_ROOT(axom::fmt::format("Comparing sensitivity for {}: {} vs {}", suffix, nameG, nameB)); + SLIC_INFO_ROOT(std::format("Comparing sensitivity for {}: {} vs {}", suffix, nameG, nameB)); // Find Gretl dual const FiniteElementDual* dualG = nullptr; @@ -585,7 +584,7 @@ TEST_F(SolidMechanicsMeshFixture, SensitivitiesComparison) ASSERT_NE(dualG, nullptr) << "Could not find Gretl dual for " << nameG; double state_diff = diff_norm(*dualG, initial_condition_sensitivitiesB.at(nameB)); - SLIC_INFO_ROOT(axom::fmt::format("Initial state {} sensitivity difference: {:.6e}", suffix, state_diff)); + SLIC_INFO_ROOT(std::format("Initial state {} sensitivity difference: {:.6e}", suffix, state_diff)); EXPECT_LT(state_diff, tol); } } diff --git a/src/smith/numerics/block_preconditioner.cpp b/src/smith/numerics/block_preconditioner.cpp index 38478a85cf..a0a543f499 100644 --- a/src/smith/numerics/block_preconditioner.cpp +++ b/src/smith/numerics/block_preconditioner.cpp @@ -7,7 +7,6 @@ #include "mfem.hpp" #include "axom/slic/core/SimpleLogger.hpp" -#include "axom/fmt.hpp" namespace smith { @@ -64,10 +63,9 @@ void BlockDiagonalPreconditioner::SetOperator(const mfem::Operator& jacobian) block_jacobian_ = dynamic_cast(&jacobian); MFEM_VERIFY(block_jacobian_, "Jacobian must be a BlockOperator"); - SLIC_ERROR_ROOT_IF( - block_jacobian_->NumRowBlocks() != num_blocks_ || block_jacobian_->NumColBlocks() != num_blocks_, - axom::fmt::format("BlockDiagonalPreconditioner solver count ({}) must match block operator size ({}x{})", - num_blocks_, block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); + SLIC_ERROR_ROOT_IF(block_jacobian_->NumRowBlocks() != num_blocks_ || block_jacobian_->NumColBlocks() != num_blocks_, + std::format("BlockDiagonalPreconditioner solver count ({}) must match block operator size ({}x{})", + num_blocks_, block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); block_offsets_.MakeRef(const_cast&>(block_jacobian_->RowOffsets())); solver_diag_ = std::make_unique(block_offsets_); @@ -219,8 +217,8 @@ void BlockTriangularPreconditioner::SetOperator(const mfem::Operator& jacobian) SLIC_ERROR_ROOT_IF( block_jacobian_->NumRowBlocks() != num_blocks_ || block_jacobian_->NumColBlocks() != num_blocks_, - axom::fmt::format("BlockTriangularPreconditioner solver count ({}) must match block operator size ({}x{})", - num_blocks_, block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); + std::format("BlockTriangularPreconditioner solver count ({}) must match block operator size ({}x{})", num_blocks_, + block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); block_offsets_.MakeRef(const_cast&>(block_jacobian_->RowOffsets())); @@ -406,11 +404,10 @@ void BlockSchurPreconditioner::SetOperator(const mfem::Operator& jacobian) MFEM_VERIFY(block_jacobian_, "Jacobian must be a BlockOperator"); SLIC_ERROR_ROOT_IF(block_jacobian_->NumRowBlocks() != 2 || block_jacobian_->NumColBlocks() != 2, - axom::fmt::format("BlockSchurPreconditioner requires a 2x2 block operator, got {}x{}", - block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); - SLIC_ERROR_ROOT_IF( - mfem_solvers_.size() != 2, - axom::fmt::format("BlockSchurPreconditioner requires exactly 2 solvers, got {}", mfem_solvers_.size())); + std::format("BlockSchurPreconditioner requires a 2x2 block operator, got {}x{}", + block_jacobian_->NumRowBlocks(), block_jacobian_->NumColBlocks())); + SLIC_ERROR_ROOT_IF(mfem_solvers_.size() != 2, + std::format("BlockSchurPreconditioner requires exactly 2 solvers, got {}", mfem_solvers_.size())); block_offsets_.MakeRef(const_cast&>(block_jacobian_->RowOffsets())); solver_diag_ = std::make_unique(block_offsets_); diff --git a/src/smith/numerics/equation_solver.cpp b/src/smith/numerics/equation_solver.cpp index d7adc577e6..89a26c9a0f 100644 --- a/src/smith/numerics/equation_solver.cpp +++ b/src/smith/numerics/equation_solver.cpp @@ -123,9 +123,9 @@ bool monolithicizeOperatorIfNeeded(const LinearSolverOptions& linear_options, mf gradient_operator = buildMonolithicMatrix(*block_gradient).release(); SLIC_DEBUG_ROOT( - axom::fmt::format("Automatically monolithicizing block Jacobian for linear solver {} with " - "preconditioner {}", - linearName(linear_options.linear_solver), preconditionerName(linear_options.preconditioner))); + std::format("Automatically monolithicizing block Jacobian for linear solver {} with " + "preconditioner {}", + linearName(linear_options.linear_solver), preconditionerName(linear_options.preconditioner))); return true; } diff --git a/src/smith/numerics/nonlinear_convergence.cpp b/src/smith/numerics/nonlinear_convergence.cpp index 454e6c8118..cd036849d3 100644 --- a/src/smith/numerics/nonlinear_convergence.cpp +++ b/src/smith/numerics/nonlinear_convergence.cpp @@ -9,7 +9,6 @@ #include #include -#include "axom/fmt.hpp" #include "axom/slic.hpp" namespace smith { diff --git a/src/smith/physics/functional_weak_form.hpp b/src/smith/physics/functional_weak_form.hpp index ea7e02787f..2eb7a29dba 100644 --- a/src/smith/physics/functional_weak_form.hpp +++ b/src/smith/physics/functional_weak_form.hpp @@ -334,8 +334,8 @@ class FunctionalWeakForm, { validateFields(fields, "jacobian"); SLIC_ERROR_IF(jacobian_weights.size() != fields.size(), - axom::fmt::format("jacobian(): jacobian_weights.size()={} but fields.size()={}", - jacobian_weights.size(), fields.size())); + std::format("jacobian(): jacobian_weights.size()={} but fields.size()={}", jacobian_weights.size(), + fields.size())); SetCurrentTimeInfoRAII clear_current_time_info_on_exit(current_time_info_, time_info); std::unique_ptr J; @@ -477,7 +477,7 @@ class FunctionalWeakForm, { if constexpr (I < sizeof...(InputSpaces)) { using Space = typename std::tuple_element>::type; - validateSpace(*input_mfem_spaces[I], axom::fmt::format("input[{}]", I)); + validateSpace(*input_mfem_spaces[I], std::format("input[{}]", I)); validateInputSpaces(input_mfem_spaces); } } @@ -489,7 +489,7 @@ class FunctionalWeakForm, if constexpr (I < sizeof...(InputSpaces)) { using Space = typename std::tuple_element>::type; validateSpace(fields[I]->space(), - axom::fmt::format("{}(): field[{}] ('{}')", method_name, I, fields[I]->name())); + std::format("{}(): field[{}] ('{}')", method_name, I, fields[I]->name())); validateFieldsRecursive(fields, method_name); } } @@ -498,8 +498,8 @@ class FunctionalWeakForm, void validateFields(const std::vector& fields, const std::string& method_name) const { SLIC_ERROR_ROOT_IF(fields.size() != sizeof...(InputSpaces), - axom::fmt::format("{}(): fields.size()={} but weak form expects {} InputSpaces", method_name, - fields.size(), sizeof...(InputSpaces))); + std::format("{}(): fields.size()={} but weak form expects {} InputSpaces", method_name, + fields.size(), sizeof...(InputSpaces))); if constexpr (sizeof...(InputSpaces) > 0) { validateFieldsRecursive<0>(fields, method_name); @@ -519,27 +519,26 @@ class FunctionalWeakForm, bool is_h1 = (fec_name.find("H1") != std::string::npos || fec_name.find("ND_") != std::string::npos || // Nedelec elements fec_name.find("Linear") != std::string::npos); - SLIC_ERROR_ROOT_IF(!is_h1, axom::fmt::format("Space '{}': Template specifies H1 family but mfem space uses '{}'", - space_name, fec_name)); + SLIC_ERROR_ROOT_IF(!is_h1, std::format("Space '{}': Template specifies H1 family but mfem space uses '{}'", + space_name, fec_name)); } else if constexpr (Space::family == Family::L2) { std::string fec_name = fec->Name(); bool is_l2 = (fec_name.find("L2") != std::string::npos || fec_name.find("DG") != std::string::npos || fec_name.find("Const") != std::string::npos); - SLIC_ERROR_ROOT_IF( - !is_l2, axom::fmt::format("Space '{}': Template specifies L2/DG family but mfem space uses '{}'", space_name, - fec_name)); + SLIC_ERROR_ROOT_IF(!is_l2, std::format("Space '{}': Template specifies L2/DG family but mfem space uses '{}'", + space_name, fec_name)); } // Validate order SLIC_ERROR_ROOT_IF(fec->GetOrder() != Space::order, - axom::fmt::format("Space '{}': Template specifies order {} but mfem space has order {}", - space_name, Space::order, fec->GetOrder())); + std::format("Space '{}': Template specifies order {} but mfem space has order {}", space_name, + Space::order, fec->GetOrder())); // Validate vector dimension SLIC_ERROR_ROOT_IF( mfem_space.GetVDim() != Space::components, - axom::fmt::format("Space '{}': Template specifies {} components but mfem space has {} components (VDim)", - space_name, Space::components, mfem_space.GetVDim())); + std::format("Space '{}': Template specifies {} components but mfem space has {} components (VDim)", space_name, + Space::components, mfem_space.GetVDim())); } /// @brief Utility to get array of jacobian functions, one for each input field in fs From 323d640e4622400e80db9ee197963547d99ce4af Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 21 Oct 2026 13:14:30 -0700 Subject: [PATCH 26/28] update docs --- src/smith/differentiable_numerics/plastic_mechanics_system.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp index 808a14713d..22b2830d08 100644 --- a/src/smith/differentiable_numerics/plastic_mechanics_system.hpp +++ b/src/smith/differentiable_numerics/plastic_mechanics_system.hpp @@ -58,7 +58,7 @@ struct PlasticMechanicsSystem : public SystemBase { L2, H1, parameter_space...>>; // Primary weak forms - std::shared_ptr solid_weak_form; + std::shared_ptr solid_weak_form; ///< Solid mechanics weak form. // Internal variable weak forms to update in staggered solve std::shared_ptr plastic_deform_weak_form; ///< Plastic deformation weak form. From 50f67e02c912766f91fe3e94a4ea3a06e39fe008 Mon Sep 17 00:00:00 2001 From: li97 Date: Wed, 21 Oct 2026 14:51:36 -0700 Subject: [PATCH 27/28] switch to include specific header --- .../tests/test_differentiable_plasticity.cpp | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp index 99b4c4b8f0..d90cb74b3b 100644 --- a/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp +++ b/src/smith/differentiable_numerics/tests/test_differentiable_plasticity.cpp @@ -4,11 +4,25 @@ // // SPDX-License-Identifier: (BSD-3-Clause) -#include #include "gtest/gtest.h" -#include "smith/smith_config.hpp" -#include "smith/smith.hpp" +#include "gretl/data_store.hpp" + +#include "smith/mesh_utils/mesh_utils.hpp" +#include "smith/infrastructure/application_manager.hpp" +#include "smith/numerics/solver_config.hpp" + +#include "smith/physics/state/state_manager.hpp" +#include "smith/physics/functional_objective.hpp" +#include "smith/physics/boundary_conditions/boundary_condition_manager.hpp" +#include "smith/physics/materials/solid_material.hpp" + +#include "smith/differentiable_numerics/nonlinear_block_solver.hpp" +#include "smith/differentiable_numerics/system_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" +#include "smith/differentiable_numerics/plastic_mechanics_system.hpp" namespace smith { From a8c78c6bd07a599cfb38f4d08b4552be40cef015 Mon Sep 17 00:00:00 2001 From: Hanyu Li <125513291+lihanyu97@users.noreply.github.com> Date: Wed, 21 Oct 2026 16:52:54 -0500 Subject: [PATCH 28/28] Update src/smith/differentiable_numerics/system_solver.cpp Co-authored-by: Chris White --- src/smith/differentiable_numerics/system_solver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smith/differentiable_numerics/system_solver.cpp b/src/smith/differentiable_numerics/system_solver.cpp index b18345181a..284505110e 100644 --- a/src/smith/differentiable_numerics/system_solver.cpp +++ b/src/smith/differentiable_numerics/system_solver.cpp @@ -13,7 +13,7 @@ #include #include -#include +#include "axom/slic.hpp" namespace smith {