From e6e4e35e6a7e774ee13a074dfd38f3b277ef14fb Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Tue, 31 Dec 2024 07:59:49 -0800 Subject: [PATCH 01/32] Add time increment input to plasticity models and get J2 material tests to pass --- .../materials/material_verification_tools.hpp | 10 +++++----- src/serac/physics/materials/solid_material.hpp | 4 ++-- .../physics/materials/tests/J2_material.cpp | 17 +++++++++++------ 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/src/serac/physics/materials/material_verification_tools.hpp b/src/serac/physics/materials/material_verification_tools.hpp index 0e90eda48e..b20d54fbfb 100644 --- a/src/serac/physics/materials/material_verification_tools.hpp +++ b/src/serac/physics/materials/material_verification_tools.hpp @@ -43,6 +43,7 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat std::function epsilon_xx, const parameter_types... parameter_functions) { double t = 0; + const double dt = t_max / double(num_steps - 1); auto state = initial_state; @@ -55,7 +56,7 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat du_dx[1][1] = epsilon_yy; du_dx[2][2] = epsilon_zz; auto state_copy = state; - auto stress = material(state_copy, du_dx, parameter_functions(t)...); + auto stress = material(state_copy, dt, du_dx, parameter_functions(t)...); return tensor{{stress[1][1], stress[2][2]}}; }; @@ -63,7 +64,6 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat output_history.reserve(num_steps); tensor dudx{}; - const double dt = t_max / double(num_steps - 1); for (size_t i = 0; i < num_steps; i++) { auto initial_guess = tensor{dudx[1][1], dudx[2][2]}; auto epsilon_yy_and_zz = find_root(sigma_yy_and_zz, initial_guess); @@ -71,7 +71,7 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat dudx[1][1] = epsilon_yy_and_zz[0]; dudx[2][2] = epsilon_yy_and_zz[1]; - auto stress = material(state, dudx, parameter_functions(t)...); + auto stress = material(state, dt, dudx, parameter_functions(t)...); output_history.push_back(tuple{t, dudx, stress, state}); t += dt; @@ -107,12 +107,12 @@ auto single_quadrature_point_test(double t_max, size_t num_steps, const Material const double dt = t_max / double(num_steps - 1); auto state = initial_state; - using output_type = decltype(std::tuple{t, state, f(0.0)..., decltype(material(state, f(0.0)...)){}}); + using output_type = decltype(std::tuple{t, state, f(0.0)..., decltype(material(state, dt, f(0.0)...)){}}); std::vector history; history.reserve(num_steps); for (size_t i = 0; i < num_steps; i++) { - auto material_output = material(state, f(t)...); + auto material_output = material(state, dt, f(t)...); history.push_back(std::tuple{t, state, f(t)..., material_output}); t += dt; } diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index bacdb15fa9..3d1c7e04f7 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -252,7 +252,7 @@ struct J2SmallStrain { /** @brief calculate the Cauchy stress, given the displacement gradient and previous material state */ template - auto operator()(State& state, const T du_dX) const + auto operator()(State& state, double dt, const tensor& du_dX) const { using std::sqrt; constexpr auto I = Identity(); @@ -313,7 +313,7 @@ struct J2 { /** @brief calculate the Cauchy stress, given the displacement gradient and previous material state */ template - auto operator()(State& state, const T du_dX) const + auto operator()(State& state, double dt, const tensor& du_dX) const { using std::sqrt; constexpr auto I = Identity(); diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index e990db9f8d..fbffa6e38a 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -138,7 +138,8 @@ TEST(J2SmallStrain, SatisfiesConsistency) Hardening hardening_law{.sigma_y = 0.1, .n = 2.0, .eps0 = 0.01}; Material material{.E = 1.0, .nu = 0.25, .hardening = hardening_law, .Hk = 0.0, .density = 1.0}; auto internal_state = Material::State{}; - tensor stress = material(internal_state, du_dx); + const double dt = 1.0; + tensor stress = material(internal_state, dt, du_dx); double mises = std::sqrt(1.5) * norm(dev(stress)); double flow_stress = hardening_law(internal_state.accumulated_plastic_strain); EXPECT_NEAR(mises, flow_stress, 1e-9 * mises); @@ -250,7 +251,9 @@ TEST(J2, DerivativeCorrectness) }}; // clang-format on - auto stress_and_tangent = material(internal_state, make_dual(H)); + const double dt = 1.0; + + auto stress_and_tangent = material(internal_state, dt, make_dual(H)); auto tangent = get_gradient(stress_and_tangent); // make sure that this load case is actually yielding @@ -260,10 +263,10 @@ TEST(J2, DerivativeCorrectness) // finite difference evaluations auto internal_state_old_p = Material::State{}; - auto stress_p = material(internal_state_old_p, H + epsilon * dH); + auto stress_p = material(internal_state_old_p, dt, H + epsilon * dH); auto internal_state_old_m = Material::State{}; - auto stress_m = material(internal_state_old_m, H - epsilon * dH); + auto stress_m = material(internal_state_old_m, dt, H - epsilon * dH); // Make sure the finite difference evaluations all took the same branch (yielding). ASSERT_GT(internal_state_old_p.accumulated_plastic_strain, 1e-3); @@ -298,6 +301,8 @@ TEST(J2, FrameIndifference) }}; //clang-format on + const double dt = 1.0; + // before we check frame indifference, make sure Q is a rotation constexpr auto I = Identity<3>(); ASSERT_LT(norm(dot(transpose(Q), Q) - I), 1e-14); @@ -306,7 +311,7 @@ TEST(J2, FrameIndifference) auto internal_state = Material::State{}; // Piola stress components in original coordinate frame - auto P = material(internal_state, H); + auto P = material(internal_state, dt, H); // make sure that this load case is actually yielding ASSERT_GT(internal_state.accumulated_plastic_strain, 1e-3); @@ -321,7 +326,7 @@ TEST(J2, FrameIndifference) auto internal_state_star = Material::State{}; // stress in second coordinate frame - auto P_star = material(internal_state_star, H_star); + auto P_star = material(internal_state_star, dt, H_star); auto error = P_star - dot(Q, P); ASSERT_LT(norm(error), 1e-13*norm(P)); From 20fdde22fc42994a57788ee6d5f5e9fa32ea977e Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Thu, 2 Jan 2025 13:42:05 -0800 Subject: [PATCH 02/32] Implement rate dependence in J2 and get it to work in material point test --- .../physics/materials/solid_material.hpp | 41 ++++++++----- .../physics/materials/tests/J2_material.cpp | 61 ++++++++++++++++--- 2 files changed, 80 insertions(+), 22 deletions(-) diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index 3d1c7e04f7..c986eaf292 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -163,12 +163,19 @@ struct NeoHookeanAdditiveSplit { double G; ///< shear modulus }; +template +auto overstress(double viscosity, T accumulated_plastic_strain_rate) +{ + return viscosity*accumulated_plastic_strain_rate; +} + /** * @brief Linear isotropic hardening law */ struct LinearHardening { double sigma_y; ///< yield strength double Hi; ///< Isotropic hardening modulus + double eta; /** * @brief Computes the flow stress @@ -177,10 +184,10 @@ struct LinearHardening { * @param accumulated_plastic_strain The uniaxial equivalent accumulated plastic strain * @return Flow stress value */ - template - auto operator()(const T accumulated_plastic_strain) const + template + auto operator()(T1 accumulated_plastic_strain, T2 accumulated_plastic_strain_rate) const { - return sigma_y + Hi * accumulated_plastic_strain; + return sigma_y + Hi * accumulated_plastic_strain + overstress(eta, accumulated_plastic_strain_rate); }; }; @@ -191,6 +198,7 @@ struct PowerLawHardening { double sigma_y; ///< yield strength double n; ///< hardening index in reciprocal form double eps0; ///< reference value of accumulated plastic strain + double eta = 0.0; /** * @brief Computes the flow stress @@ -199,11 +207,11 @@ struct PowerLawHardening { * @param accumulated_plastic_strain The uniaxial equivalent accumulated plastic strain * @return Flow stress value */ - template - auto operator()(const T accumulated_plastic_strain) const + template + auto operator()(T1 accumulated_plastic_strain, T2 accumulated_plastic_strain_rate) const { using std::pow; - return sigma_y * pow(1.0 + accumulated_plastic_strain / eps0, 1.0 / n); + return sigma_y * pow(1.0 + accumulated_plastic_strain / eps0, 1.0 / n) + overstress(eta, accumulated_plastic_strain_rate); }; }; @@ -216,6 +224,7 @@ struct VoceHardening { double sigma_y; ///< yield strength double sigma_sat; ///< saturation value of flow strength double strain_constant; ///< The constant dictating how fast the exponential decays + double eta = 0.0; ///< linear viscosity for rate dependence /** * @brief Computes the flow stress @@ -224,11 +233,13 @@ struct VoceHardening { * @param accumulated_plastic_strain The uniaxial equivalent accumulated plastic strain * @return Flow stress value */ - template - auto operator()(const T accumulated_plastic_strain) const + template + auto operator()(T1 epsilon_p, T2 epsilon_p_dot) const { using std::exp; - return sigma_sat - (sigma_sat - sigma_y) * exp(-accumulated_plastic_strain / strain_constant); + auto Y_eq = sigma_sat - (sigma_sat - sigma_y) * exp(-epsilon_p / strain_constant); + auto Y_vis = overstress(eta, epsilon_p_dot); + return Y_eq + Y_vis; }; }; @@ -269,8 +280,8 @@ struct J2SmallStrain { // (ii) admissibility const double eqps_old = state.accumulated_plastic_strain; - auto residual = [eqps_old, G, *this](auto delta_eqps, auto trial_q) { - return trial_q - (3.0 * G + Hk) * delta_eqps - this->hardening(eqps_old + delta_eqps); + auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_q) { + return trial_q - (3.0 * G + Hk) * delta_eqps - this->hardening(eqps_old + delta_eqps, delta_eqps/dt); }; if (residual(0.0, get_value(q)) > tol * hardening.sigma_y) { // (iii) return mapping @@ -280,7 +291,7 @@ struct J2SmallStrain { // variables, the return map won't be repeated. ScalarSolverOptions opts{.xtol = 0, .rtol = tol * hardening.sigma_y, .max_iter = 25}; double lower_bound = 0.0; - double upper_bound = (get_value(q) - hardening(eqps_old)) / (3.0 * G + Hk); + double upper_bound = (get_value(q) - hardening(eqps_old, 0.0)) / (3.0 * G + Hk); auto [delta_eqps, status] = solve_scalar_equation(residual, 0.0, lower_bound, upper_bound, opts, q); auto Np = 1.5 * eta / q; @@ -332,8 +343,8 @@ struct J2 { // (ii) admissibility const double eqps_old = state.accumulated_plastic_strain; - auto residual = [eqps_old, G, *this](auto delta_eqps, auto trial_mises) { - return trial_mises - 3.0 * G * delta_eqps - this->hardening(eqps_old + delta_eqps); + auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_mises) { + return trial_mises - 3.0 * G * delta_eqps - this->hardening(eqps_old + delta_eqps, delta_eqps/dt); }; if (residual(0.0, get_value(q)) > tol * hardening.sigma_y) { // (iii) return mapping @@ -343,7 +354,7 @@ struct J2 { // variables, the return map won't be repeated. ScalarSolverOptions opts{.xtol = 0, .rtol = tol * hardening.sigma_y, .max_iter = 25}; double lower_bound = 0.0; - double upper_bound = (get_value(q) - hardening(eqps_old)) / (3.0 * G); + double upper_bound = (get_value(q) - hardening(eqps_old, 0.0)) / (3.0 * G); auto [delta_eqps, status] = solve_scalar_equation(residual, 0.0, lower_bound, upper_bound, opts, q); auto Np = 1.5 * s / q; diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index fbffa6e38a..42b0921d7f 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -16,6 +16,9 @@ #include "serac/numerics/functional/tensor.hpp" #include "serac/physics/materials/material_verification_tools.hpp" +#include +#include + namespace serac { using namespace serac; @@ -86,7 +89,7 @@ TEST(J2SmallStrain, Verification) using Hardening = solid_mechanics::LinearHardening; using Material = solid_mechanics::J2SmallStrain; - Hardening hardening{.sigma_y = sigma_y, .Hi = 0.0}; + Hardening hardening{.sigma_y = sigma_y, .Hi = 0.0, .eta = 0.0}; Material material{.E = E, .nu = nu, .hardening = hardening, .Hk = 0.0, .density = 1.0}; Material::State initial_state{}; @@ -115,9 +118,10 @@ TEST(J2SmallStrain, Verification) TEST(J2, PowerLawHardeningWorksWithDuals) { double sigma_y = 1.0; - solid_mechanics::PowerLawHardening hardening_law{.sigma_y = sigma_y, .n = 2.0, .eps0 = 0.01}; + solid_mechanics::PowerLawHardening hardening_law{.sigma_y = sigma_y, .n = 2.0, .eps0 = 0.01, .eta = 0.0}; double eqps = 0.1; - auto flow_stress = hardening_law(make_dual(eqps)); + double eqps_dot = 1.0; + auto flow_stress = hardening_law(make_dual(eqps), eqps_dot); EXPECT_GT(flow_stress.value, sigma_y); EXPECT_GT(flow_stress.gradient, 0.0); }; @@ -135,13 +139,13 @@ TEST(J2SmallStrain, SatisfiesConsistency) using Hardening = solid_mechanics::PowerLawHardening; using Material = solid_mechanics::J2SmallStrain; - Hardening hardening_law{.sigma_y = 0.1, .n = 2.0, .eps0 = 0.01}; + Hardening hardening_law{.sigma_y = 0.1, .n = 2.0, .eps0 = 0.01, .eta = 0.0}; Material material{.E = 1.0, .nu = 0.25, .hardening = hardening_law, .Hk = 0.0, .density = 1.0}; auto internal_state = Material::State{}; const double dt = 1.0; tensor stress = material(internal_state, dt, du_dx); double mises = std::sqrt(1.5) * norm(dev(stress)); - double flow_stress = hardening_law(internal_state.accumulated_plastic_strain); + double flow_stress = hardening_law(internal_state.accumulated_plastic_strain, 0.0); EXPECT_NEAR(mises, flow_stress, 1e-9 * mises); double twoG = material.E / (1 + material.nu); @@ -159,7 +163,7 @@ TEST(J2SmallStrain, Uniaxial) double sigma_y = 0.01; double Hi = E / 100.0; - Hardening hardening{.sigma_y = sigma_y, .Hi = Hi}; + Hardening hardening{.sigma_y = sigma_y, .Hi = Hi, .eta = 0.0}; Material material{.E = E, .nu = nu, .hardening = hardening, .Hk = 0.0, .density = 1.0}; auto internal_state = Material::State{}; @@ -196,7 +200,7 @@ TEST(J2, Uniaxial) double sigma_y = 0.01; double Hi = E / 100.0; - Hardening hardening{.sigma_y = sigma_y, .Hi = Hi}; + Hardening hardening{.sigma_y = sigma_y, .Hi = Hi, .eta = 0.0}; Material material{.E = E, .nu = 0.25, .hardening = hardening, .density = 1.0}; auto internal_state = Material::State{}; @@ -336,6 +340,49 @@ TEST(J2, FrameIndifference) ASSERT_LT(norm(error), 1e-13*norm(internal_state.Fpinv)); } +TEST(J2, PlotOutput) +{ + /* Log strain J2 plasticity has the nice feature that the exact uniaxial stress solution from + small strain plasticity are applicable, if you replace the lineasr strain with log strain + and use the Kirchhoff stress as the output. + */ + + using Hardening = solid_mechanics::VoceHardening; + using Material = solid_mechanics::J2; + + constexpr double E = 1.0; + constexpr double sigma_y = 0.001; + constexpr double sigma_sat = 3*sigma_y; + constexpr double strain_constant = 10*sigma_y/E; + constexpr double eta = 1e-2; + + constexpr double max_strain = 3*strain_constant; + constexpr double strain_rate = 1e-1; + const std::string tag = "m1"; + constexpr double t_max = max_strain/strain_rate; + + + + Hardening hardening{sigma_y, sigma_sat, strain_constant, eta}; + Material material{.E = E, .nu = 0.25, .hardening = hardening, .density = 1.0}; + + auto internal_state = Material::State{}; + auto strain = [=](double t) { return max_strain*t/t_max; }; + auto response_history = uniaxial_stress_test(t_max, 100, material, internal_state, strain); + + std::ofstream file("J2_" + tag + ".txt"); + + for (auto r : response_history) { + auto [t, dudx, P, state] = r; + auto TK = dot(P, transpose(dudx + Identity<3>())); + double e = std::log1p(dudx[0][0]); // log strain + double s = TK[0][0]; // Kirchhoff stress + double pe = -std::log(get<3>(r).Fpinv[0][0]); // plastic strain + file << t << " " << e << " " << s << std::endl; + } + file.close(); +}; + } // namespace serac int main(int argc, char* argv[]) From 4f25d8e67e45de654f6285d9c9590e2208709108 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 3 Jan 2025 11:04:09 -0800 Subject: [PATCH 03/32] Pass time increment to material models in solid mechanics --- src/serac/physics/base_physics.hpp | 2 ++ src/serac/physics/materials/solid_material.hpp | 8 ++++---- src/serac/physics/solid_mechanics.hpp | 14 +++++++++++--- src/serac/physics/tests/solid_statics_patch.cpp | 16 ++++++++++------ 4 files changed, 27 insertions(+), 13 deletions(-) diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index b7925f260c..bfee79de35 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -618,6 +618,8 @@ class BasePhysics { */ double time_; + double time_prev_; + /** * @brief The maximum time reached for the forward solver */ diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index c986eaf292..067ca15af7 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -36,7 +36,7 @@ struct LinearIsotropic { * @return The Cauchy stress */ template - SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor& du_dX) const + SERAC_HOST_DEVICE auto operator()(State& /* state */, double /* dt */, const tensor& du_dX) const { auto I = Identity(); auto lambda = K - (2.0 / 3.0) * G; @@ -72,7 +72,7 @@ struct StVenantKirchhoff { * @return The Cauchy stress */ template - auto operator()(State&, const tensor& grad_u) const + auto operator()(State&, double /* dt */, const tensor& grad_u) const { static constexpr auto I = Identity(); auto F = grad_u + I; @@ -107,7 +107,7 @@ struct NeoHookean { * @return The Cauchy stress */ template - SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor& du_dX) const + SERAC_HOST_DEVICE auto operator()(State& /* state */, double /* dt */, const tensor& du_dX) const { using std::log1p; constexpr auto I = Identity(); @@ -140,7 +140,7 @@ struct NeoHookeanAdditiveSplit { * @return The first Piola stress */ template - SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor& du_dX) const + SERAC_HOST_DEVICE auto operator()(State& /* state */, double /* dt */, const tensor& du_dX) const { using std::pow; constexpr auto I = Identity(); diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 564e79009e..97265ddd4b 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -740,11 +740,16 @@ class SolidMechanics, std::integer_se template struct MaterialStressFunctor { /// @brief Constructor for the functor - MaterialStressFunctor(Material material) : material_(material) {} + MaterialStressFunctor(Material material, SolidMechanics& owner) : material_(material), owner_(&owner) {} /// @brief Material model Material material_; + /// @brief Pointer to the `SolidMechanics` instance that owns the current `MaterialStressFunctor` + /// + /// This is needed to compute the time increment from the enclosing + const SolidMechanics* owner_; + /** * @brief Material stress response call * @@ -766,8 +771,9 @@ class SolidMechanics, std::integer_se { auto du_dX = get(displacement); auto d2u_dt2 = get(acceleration); + const double dt = owner_->time_ - owner_->time_prev_; - auto stress = material_(state, du_dX, params...); + auto stress = material_(state, dt, du_dX, params...); return serac::tuple{material_.density * d2u_dt2, stress}; } @@ -804,7 +810,7 @@ class SolidMechanics, std::integer_se { static_assert(std::is_same_v || std::is_same_v, "invalid quadrature data provided in setMaterial()"); - MaterialStressFunctor material_functor(material); + MaterialStressFunctor material_functor(material, *this); residual_->AddDomainIntegral( Dimension{}, DependsOn<0, 1, @@ -1170,6 +1176,7 @@ class SolidMechanics, std::integer_se if (is_quasistatic_) { quasiStaticSolve(dt); } else { + time_prev_ = time_; // The current ode interface tracks 2 times, one internally which we have a handle to via time_, // and one here via the step interface. // We are ignoring this one, and just using the internal version for now. @@ -1515,6 +1522,7 @@ class SolidMechanics, std::integer_se // warm start must be called prior to the time update so that the previous Jacobians can be used consistently // throughout. warmStartDisplacement(dt); + time_prev_ = time_; time_ += dt; // this method is essentially equivalent to the 1-liner diff --git a/src/serac/physics/tests/solid_statics_patch.cpp b/src/serac/physics/tests/solid_statics_patch.cpp index a8b19c7971..7b6906c4db 100644 --- a/src/serac/physics/tests/solid_statics_patch.cpp +++ b/src/serac/physics/tests/solid_statics_patch.cpp @@ -77,7 +77,7 @@ class AffineSolution { * @param essential_boundaries Boundary attributes on which essential boundary conditions are desired */ template - void applyLoads(const Material& material, SolidMechanics& sf, Domain essential_boundary) const + void applyLoads(const Material& material, SolidMechanics& sf, Domain essential_boundary, double dt) const { // essential BCs auto ebc_func = [*this](tensor X, double) { return this->eval(X); }; @@ -86,7 +86,7 @@ class AffineSolution { // natural BCs typename Material::State state; - tensor P = material(state, A); + tensor P = material(state, dt, A); auto traction = [P](auto, auto n0, auto) { return dot(P, n0); }; Domain entire_boundary = EntireBoundary(sf.mesh()); sf.setTraction(traction, entire_boundary); @@ -216,14 +216,16 @@ double solution_error(PatchBoundaryCondition bc) solid.setMaterial(mat, domain); + constexpr double dt = 1.0; + Domain essential_boundary = Domain::ofBoundaryElements(pmesh, by_attr(essentialBoundaryAttributes(bc))); - exact_displacement.applyLoads(mat, solid, essential_boundary); + exact_displacement.applyLoads(mat, solid, essential_boundary, dt); // Finalize the data structures solid.completeSetup(); // Perform the quasi-static solve - solid.advanceTimestep(1.0); + solid.advanceTimestep(dt); // Output solution for debugging // solid.outputStateToDisk("paraview_output"); @@ -338,7 +340,9 @@ double pressure_error() }); // clang-format on - tensor P = mat(state, H); + constexpr double dt = 1.0; + + tensor P = mat(state, dt, H); auto F = H + Identity(); auto sigma = dot(P, transpose(F)) / det(F); double pressure = -sigma[0][0]; @@ -358,7 +362,7 @@ double pressure_error() solid.completeSetup(); // Perform the quasi-static solve - solid.advanceTimestep(1.0); + solid.advanceTimestep(dt); solid.outputStateToDisk(); From 459fff0806c73bb2099ae1ebe705ee8dcaab6391 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 3 Jan 2025 15:17:22 -0800 Subject: [PATCH 04/32] Implement first draft of uniaxial tension example --- examples/CMakeLists.txt | 10 ++ examples/uniaxial/CMakeLists.txt | 13 +++ examples/uniaxial/uniaxial_tension.cpp | 135 +++++++++++++++++++++++++ 3 files changed, 158 insertions(+) create mode 100644 examples/uniaxial/CMakeLists.txt create mode 100644 examples/uniaxial/uniaxial_tension.cpp diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 2aceea0fff..4d9eabe428 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -50,3 +50,13 @@ install( DESTINATION examples/serac/buckling ) + +# Add the uniaxial examples +add_subdirectory(uniaxial) + +install( + FILES + uniaxial/uniaxial.cpp + DESTINATION + examples/serac/uniaxial +) diff --git a/examples/uniaxial/CMakeLists.txt b/examples/uniaxial/CMakeLists.txt new file mode 100644 index 0000000000..d72d1cd598 --- /dev/null +++ b/examples/uniaxial/CMakeLists.txt @@ -0,0 +1,13 @@ +# Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and +# other Serac Project Developers. See the top-level LICENSE file for +# details. +# +# SPDX-License-Identifier: (BSD-3-Clause) + +blt_add_executable( NAME uniaxial_tension + SOURCES uniaxial_tension.cpp + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON serac_physics serac_mesh + ) + +# add serac_functional to DEPENDS_ON \ No newline at end of file diff --git a/examples/uniaxial/uniaxial_tension.cpp b/examples/uniaxial/uniaxial_tension.cpp new file mode 100644 index 0000000000..0328663687 --- /dev/null +++ b/examples/uniaxial/uniaxial_tension.cpp @@ -0,0 +1,135 @@ +// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include + +#include "axom/inlet.hpp" +#include "axom/slic/core/SimpleLogger.hpp" +#include "mfem.hpp" + +#include "serac/infrastructure/terminator.hpp" +#include "serac/mesh/mesh_utils.hpp" +#include "serac/physics/boundary_conditions/components.hpp" +#include "serac/physics/materials/solid_material.hpp" +#include "serac/physics/solid_mechanics.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/serac_config.hpp" + +int main(int argc, char* argv[]) +{ + serac::initialize(argc, argv); + + constexpr int p = 2; + constexpr int dim = 3; + constexpr double x_length = 1.0; + constexpr double y_length = 0.1; + constexpr double z_length = 0.1; + constexpr int elements_in_x = 10; + constexpr int elements_in_y = 1; + constexpr int elements_in_z = 1; + + int serial_refinements = 0; + int parallel_refinements = 0; + int time_steps = 20; + double strain_rate = 1e-3; + + constexpr double E = 1.0; + constexpr double nu = 0.25; + + constexpr double sigma_y = 0.001; + constexpr double sigma_sat = 3.0*sigma_y; + constexpr double strain_constant = 10*sigma_y/E; + constexpr double eta = 1e-2; + + constexpr double max_strain = 3*strain_constant; + + // Handle command line arguments + axom::CLI::App app{"Plane strain uniaxial extension of a bar."}; + // Mesh options + app.add_option("--serial-refinements", serial_refinements, "Serial refinement steps", true); + app.add_option("--parallel-refinements", parallel_refinements, "Parallel refinement steps", true); + app.add_option("--time-steps", time_steps, + "Number of time steps to divide simulation", true); + app.add_option("--strain_rate", strain_rate, "Nominal strain rate", true); + + double max_time = max_strain/strain_rate; + + // Create DataStore + const std::string simulation_tag = "uniaxial"; + const std::string mesh_tag = simulation_tag + "mesh"; + axom::sidre::DataStore datastore; + serac::StateManager::initialize(datastore, simulation_tag + "_data"); + + auto serial_mesh = std::make_unique( + mfem::Mesh::MakeCartesian3D(elements_in_x, elements_in_y, elements_in_z, + mfem::Element::QUADRILATERAL, true, x_length, + y_length, z_length)); + serial_mesh->Print(); + auto mesh = serac::mesh::refineAndDistribute(std::move(*serial_mesh), serial_refinements, parallel_refinements); + auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); + + // create boundary domains for boundary conditions + auto fix_x = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(1)); + auto fix_y = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(2)); + auto fix_z = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(4)); + auto apply_displacement = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(3)); + serac::Domain whole_mesh = serac::EntireDomain(pmesh); + + serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 1}; +#ifndef MFEM_USE_STRUMPACK + SLIC_INFO_ROOT("Uniaxial app requires MFEM built with strumpack."); + return 1; +#endif + + serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton, + .relative_tol = 1.0e-10, + .absolute_tol = 1.0e-12, + .max_iterations = 200, + .print_level = 1}; + + serac::SolidMechanics solid_solver( + nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, simulation_tag, + mesh_tag); + + using Hardening = serac::solid_mechanics::VoceHardening; + using Material = serac::solid_mechanics::J2; + + Hardening hardening{sigma_y, sigma_sat, strain_constant, eta}; + Material material{E, nu, hardening, .density = 1.0}; + + auto internal_states = solid_solver.createQuadratureDataBuffer(Material::State{}, whole_mesh); + + solid_solver.setMaterial(material, whole_mesh, internal_states); + + solid_solver.setFixedBCs(fix_x, serac::Component::X); + solid_solver.setFixedBCs(fix_y, serac::Component::Y); + solid_solver.setFixedBCs(fix_z, serac::Component::Z); + auto applied_displacement = [strain_rate](serac::vec3, double t){ return serac::vec3{strain_rate*x_length*t, 0., 0.}; }; + solid_solver.setDisplacementBCs(applied_displacement, apply_displacement, serac::Component::X); + + solid_solver.completeSetup(); + + std::string paraview_tag = simulation_tag + "_paraview"; + solid_solver.outputStateToDisk(paraview_tag); + + double dt = max_time/(time_steps - 1); + + for (int i = 0; i < time_steps; ++i) { + SLIC_INFO_ROOT("------------------------------------------"); + SLIC_INFO_ROOT(axom::fmt::format("TIME STEP {}", i)); + SLIC_INFO_ROOT(axom::fmt::format("time = {} (out of {})", solid_solver.time(), max_time)); + serac::logger::flush(); + + solid_solver.advanceTimestep(dt); + + // Output the sidre-based plot files + solid_solver.outputStateToDisk(paraview_tag); + } + + serac::exitGracefully(); + + return 0; +} From b0eab5a74ea22588458f72bccdb41088960223ba Mon Sep 17 00:00:00 2001 From: Chris White Date: Fri, 3 Jan 2025 16:57:06 -0800 Subject: [PATCH 05/32] clean up example logic --- examples/CMakeLists.txt | 30 ---------------- examples/buckling/CMakeLists.txt | 8 +++++ examples/contact/CMakeLists.txt | 43 ++++++++++++----------- examples/simple_conduction/CMakeLists.txt | 7 ++++ examples/uniaxial/CMakeLists.txt | 9 ++++- 5 files changed, 46 insertions(+), 51 deletions(-) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 4d9eabe428..76ca23b338 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -21,42 +21,12 @@ install( # Add the simple conduction examples add_subdirectory(simple_conduction) -install( - FILES - simple_conduction/without_input_file.cpp - DESTINATION - examples/serac/simple_conduction -) - # Add the contact examples add_subdirectory(contact) -install( - FILES - contact/beam_bending.cpp - contact/ironing.cpp - contact/sphere.cpp - contact/twist.cpp - DESTINATION - examples/serac/contact -) - # Add the buckling examples add_subdirectory(buckling) -install( - FILES - buckling/cylinder.cpp - DESTINATION - examples/serac/buckling -) - # Add the uniaxial examples add_subdirectory(uniaxial) -install( - FILES - uniaxial/uniaxial.cpp - DESTINATION - examples/serac/uniaxial -) diff --git a/examples/buckling/CMakeLists.txt b/examples/buckling/CMakeLists.txt index c67b60419a..2f0cfc6d33 100644 --- a/examples/buckling/CMakeLists.txt +++ b/examples/buckling/CMakeLists.txt @@ -10,4 +10,12 @@ if(TRIBOL_FOUND AND PETSC_FOUND) OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} DEPENDS_ON serac_physics serac_mesh ) + + install( + FILES + cylinder.cpp + DESTINATION + examples/serac/buckling + ) + endif() diff --git a/examples/contact/CMakeLists.txt b/examples/contact/CMakeLists.txt index d514654c5d..a67800cef0 100644 --- a/examples/contact/CMakeLists.txt +++ b/examples/contact/CMakeLists.txt @@ -5,24 +5,27 @@ # SPDX-License-Identifier: (BSD-3-Clause) if(TRIBOL_FOUND) - blt_add_executable( NAME contact_beam_bending - SOURCES beam_bending.cpp - OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} - DEPENDS_ON serac_physics serac_mesh - ) - blt_add_executable( NAME contact_ironing - SOURCES ironing.cpp - OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} - DEPENDS_ON serac_physics serac_mesh - ) - blt_add_executable( NAME contact_sphere - SOURCES sphere.cpp - OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} - DEPENDS_ON serac_physics serac_mesh - ) - blt_add_executable( NAME contact_twist - SOURCES twist.cpp - OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} - DEPENDS_ON serac_physics serac_mesh - ) + set(CONTACT_EXAMPLES_SOURCES + beam_bending.cpp + ironing.cpp + sphere.cpp + twist.cpp + ) + + foreach(filename ${arg_CONTACT_EXAMPLES_SOURCES}) + get_filename_component(example_name ${filename} NAME_WE) + + blt_add_executable(NAME contact_${example_name} + SOURCES ${filename} + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON serac_physics serac_mesh) + endforeach() + + install( + FILES + ${CONTACT_EXAMPLES_SOURCES} + DESTINATION + examples/serac/contact + ) + endif() diff --git a/examples/simple_conduction/CMakeLists.txt b/examples/simple_conduction/CMakeLists.txt index 6c58adbf11..91893180df 100644 --- a/examples/simple_conduction/CMakeLists.txt +++ b/examples/simple_conduction/CMakeLists.txt @@ -10,6 +10,13 @@ blt_add_executable( NAME simple_conduction_without_input_file DEPENDS_ON serac_physics serac_mesh serac_functional ) +install( + FILES + without_input_file.cpp + DESTINATION + examples/serac/simple_conduction +) + if(SERAC_ENABLE_TESTS) blt_add_test(NAME simple_conduction_without_input_file COMMAND simple_conduction_without_input_file diff --git a/examples/uniaxial/CMakeLists.txt b/examples/uniaxial/CMakeLists.txt index d72d1cd598..a0195ddae9 100644 --- a/examples/uniaxial/CMakeLists.txt +++ b/examples/uniaxial/CMakeLists.txt @@ -10,4 +10,11 @@ blt_add_executable( NAME uniaxial_tension DEPENDS_ON serac_physics serac_mesh ) -# add serac_functional to DEPENDS_ON \ No newline at end of file +# add serac_functional to DEPENDS_ON + +install( + FILES + uniaxial_tension.cpp + DESTINATION + examples/serac/uniaxial +) From 31473cfd90673301c72e6cecde9f899fdcd47f94 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 3 Jan 2025 17:44:54 -0800 Subject: [PATCH 06/32] Fix app and get working --- examples/uniaxial/uniaxial_tension.cpp | 43 +++++++++++++++++--------- 1 file changed, 28 insertions(+), 15 deletions(-) diff --git a/examples/uniaxial/uniaxial_tension.cpp b/examples/uniaxial/uniaxial_tension.cpp index 0328663687..035207d934 100644 --- a/examples/uniaxial/uniaxial_tension.cpp +++ b/examples/uniaxial/uniaxial_tension.cpp @@ -5,6 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) #include +#include #include "axom/inlet.hpp" #include "axom/slic/core/SimpleLogger.hpp" @@ -33,7 +34,7 @@ int main(int argc, char* argv[]) int serial_refinements = 0; int parallel_refinements = 0; - int time_steps = 20; + int time_steps = 100; double strain_rate = 1e-3; constexpr double E = 1.0; @@ -53,7 +54,7 @@ int main(int argc, char* argv[]) app.add_option("--parallel-refinements", parallel_refinements, "Parallel refinement steps", true); app.add_option("--time-steps", time_steps, "Number of time steps to divide simulation", true); - app.add_option("--strain_rate", strain_rate, "Nominal strain rate", true); + app.add_option("--strain-rate", strain_rate, "Nominal strain rate", true); double max_time = max_strain/strain_rate; @@ -63,26 +64,17 @@ int main(int argc, char* argv[]) axom::sidre::DataStore datastore; serac::StateManager::initialize(datastore, simulation_tag + "_data"); - auto serial_mesh = std::make_unique( - mfem::Mesh::MakeCartesian3D(elements_in_x, elements_in_y, elements_in_z, - mfem::Element::QUADRILATERAL, true, x_length, - y_length, z_length)); - serial_mesh->Print(); - auto mesh = serac::mesh::refineAndDistribute(std::move(*serial_mesh), serial_refinements, parallel_refinements); + auto mesh = serac::mesh::refineAndDistribute(serac::buildCuboidMesh(elements_in_x, elements_in_y, elements_in_z, x_length, y_length, z_length), serial_refinements, parallel_refinements); auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); // create boundary domains for boundary conditions - auto fix_x = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(1)); + auto fix_x = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(5)); auto fix_y = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(2)); - auto fix_z = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(4)); + auto fix_z = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(1)); auto apply_displacement = serac::Domain::ofBoundaryElements(pmesh, serac::by_attr(3)); serac::Domain whole_mesh = serac::EntireDomain(pmesh); - serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 1}; -#ifndef MFEM_USE_STRUMPACK - SLIC_INFO_ROOT("Uniaxial app requires MFEM built with strumpack."); - return 1; -#endif + serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 0}; serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton, .relative_tol = 1.0e-10, @@ -117,6 +109,22 @@ int main(int argc, char* argv[]) double dt = max_time/(time_steps - 1); + // get nodes and dofs to compute total force + mfem::Array dof_list = apply_displacement.dof_list(&solid_solver.displacement().space()); + solid_solver.displacement().space().DofsToVDofs(0, dof_list); + + auto compute_tensile_force = [&dof_list](const serac::FiniteElementDual& reaction) -> double { + double R{}; + for (int i = 0; i < dof_list.Size(); i++) { + R += reaction(dof_list[i]); + } + return R; + }; + + std::ofstream file("uniaxial.txt"); + file << "# time displacement force" << std::endl; + file << 0 << " " << 0 << " " << 0 << std::endl; + for (int i = 0; i < time_steps; ++i) { SLIC_INFO_ROOT("------------------------------------------"); SLIC_INFO_ROOT(axom::fmt::format("TIME STEP {}", i)); @@ -127,8 +135,13 @@ int main(int argc, char* argv[]) // Output the sidre-based plot files solid_solver.outputStateToDisk(paraview_tag); + + double u = applied_displacement(serac::vec3{}, solid_solver.time())[0]; + double f = compute_tensile_force(solid_solver.dual("reactions")); + file << solid_solver.time() << " " << u << " " << f << std::endl; } + file.close(); serac::exitGracefully(); return 0; From beabdc70796d67eff04cb7fc8db900da1eef53e1 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Mon, 6 Jan 2025 09:14:57 -0800 Subject: [PATCH 07/32] Fix warning message that command-line options were not parsed --- examples/uniaxial/uniaxial_tension.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/examples/uniaxial/uniaxial_tension.cpp b/examples/uniaxial/uniaxial_tension.cpp index 035207d934..5583c8a7bb 100644 --- a/examples/uniaxial/uniaxial_tension.cpp +++ b/examples/uniaxial/uniaxial_tension.cpp @@ -52,9 +52,13 @@ int main(int argc, char* argv[]) // Mesh options app.add_option("--serial-refinements", serial_refinements, "Serial refinement steps", true); app.add_option("--parallel-refinements", parallel_refinements, "Parallel refinement steps", true); - app.add_option("--time-steps", time_steps, - "Number of time steps to divide simulation", true); + app.add_option("--time-steps", time_steps, "Number of time steps to divide simulation", true); app.add_option("--strain-rate", strain_rate, "Nominal strain rate", true); + app.set_help_flag("--help"); + + CLI11_PARSE(app, argc, argv); + + SLIC_INFO_ROOT(axom::fmt::format("strain rate: {}", strain_rate)); double max_time = max_strain/strain_rate; @@ -64,7 +68,10 @@ int main(int argc, char* argv[]) axom::sidre::DataStore datastore; serac::StateManager::initialize(datastore, simulation_tag + "_data"); - auto mesh = serac::mesh::refineAndDistribute(serac::buildCuboidMesh(elements_in_x, elements_in_y, elements_in_z, x_length, y_length, z_length), serial_refinements, parallel_refinements); + auto mesh = serac::mesh::refineAndDistribute( + serac::buildCuboidMesh(elements_in_x, elements_in_y, elements_in_z, + x_length, y_length, z_length), + serial_refinements, parallel_refinements); auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); // create boundary domains for boundary conditions @@ -113,7 +120,7 @@ int main(int argc, char* argv[]) mfem::Array dof_list = apply_displacement.dof_list(&solid_solver.displacement().space()); solid_solver.displacement().space().DofsToVDofs(0, dof_list); - auto compute_tensile_force = [&dof_list](const serac::FiniteElementDual& reaction) -> double { + auto compute_net_force = [&dof_list](const serac::FiniteElementDual& reaction) -> double { double R{}; for (int i = 0; i < dof_list.Size(); i++) { R += reaction(dof_list[i]); @@ -137,7 +144,7 @@ int main(int argc, char* argv[]) solid_solver.outputStateToDisk(paraview_tag); double u = applied_displacement(serac::vec3{}, solid_solver.time())[0]; - double f = compute_tensile_force(solid_solver.dual("reactions")); + double f = compute_net_force(solid_solver.dual("reactions")); file << solid_solver.time() << " " << u << " " << f << std::endl; } From 13d2503d136ad0b9f81b3473245e7ba41bc6b236 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Mon, 6 Jan 2025 15:31:43 -0800 Subject: [PATCH 08/32] Fix some of the NaNs generated in J2 caused by differentiating at singular points --- examples/uniaxial/uniaxial_tension.cpp | 3 +++ .../physics/materials/solid_material.hpp | 20 +++++++++++++------ 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/examples/uniaxial/uniaxial_tension.cpp b/examples/uniaxial/uniaxial_tension.cpp index 5583c8a7bb..15de518ddb 100644 --- a/examples/uniaxial/uniaxial_tension.cpp +++ b/examples/uniaxial/uniaxial_tension.cpp @@ -19,8 +19,11 @@ #include "serac/physics/state/state_manager.hpp" #include "serac/serac_config.hpp" +#include + int main(int argc, char* argv[]) { + feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); serac::initialize(argc, argv); constexpr int p = 2; diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index 067ca15af7..ac201460e6 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -281,7 +281,8 @@ struct J2SmallStrain { // (ii) admissibility const double eqps_old = state.accumulated_plastic_strain; auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_q) { - return trial_q - (3.0 * G + Hk) * delta_eqps - this->hardening(eqps_old + delta_eqps, delta_eqps/dt); + auto eqps_dot = dt > 0.0? delta_eqps/dt : 0.0*delta_eqps; + return trial_q - (3.0 * G + Hk) * delta_eqps - this->hardening(eqps_old + delta_eqps, eqps_dot); }; if (residual(0.0, get_value(q)) > tol * hardening.sigma_y) { // (iii) return mapping @@ -339,14 +340,15 @@ struct J2 { // small strain one. auto p = K * tr(Ee); auto s = 2.0 * G * dev(Ee); - auto q = sqrt(1.5) * norm(s); + double q_val = sqrt(1.5) * norm(get_value(s)); // (ii) admissibility const double eqps_old = state.accumulated_plastic_strain; auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_mises) { - return trial_mises - 3.0 * G * delta_eqps - this->hardening(eqps_old + delta_eqps, delta_eqps/dt); + auto eqps_dot = dt > 0.0? delta_eqps/dt : 0.0*delta_eqps; + return trial_mises - 3.0 * G * delta_eqps - this->hardening(eqps_old + delta_eqps, eqps_dot); }; - if (residual(0.0, get_value(q)) > tol * hardening.sigma_y) { + if (residual(0.0, q_val) > tol * hardening.sigma_y) { // (iii) return mapping // Note the tolerance for convergence is the same as the tolerance for entering the return map. @@ -354,8 +356,14 @@ struct J2 { // variables, the return map won't be repeated. ScalarSolverOptions opts{.xtol = 0, .rtol = tol * hardening.sigma_y, .max_iter = 25}; double lower_bound = 0.0; - double upper_bound = (get_value(q) - hardening(eqps_old, 0.0)) / (3.0 * G); - auto [delta_eqps, status] = solve_scalar_equation(residual, 0.0, lower_bound, upper_bound, opts, q); + double upper_bound = (q_val - hardening(eqps_old, 0.0)) / (3.0 * G); + // safe to compute derivative of mises stress now that we're in yielding branch + auto q = sqrt(1.5) * norm(s); + auto [delta_eqps, status] = solve_scalar_equation(residual, 0.0, lower_bound, upper_bound, opts, q); + + if (!status.converged) { + SLIC_WARNING("J2 solve failed to converge."); + } auto Np = 1.5 * s / q; From f84ed453cd56068960144c96a53da67e971c87e7 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Mon, 13 Jan 2025 11:49:37 -0800 Subject: [PATCH 09/32] Pass time incement only to material functor, instead of entire solid mechancis module --- src/serac/physics/base_physics.cpp | 1 + src/serac/physics/base_physics.hpp | 5 ++++- src/serac/physics/solid_mechanics.hpp | 17 ++++++++--------- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/serac/physics/base_physics.cpp b/src/serac/physics/base_physics.cpp index aca90b8dc7..1f0867ec2b 100644 --- a/src/serac/physics/base_physics.cpp +++ b/src/serac/physics/base_physics.cpp @@ -64,6 +64,7 @@ void BasePhysics::initializeBasePhysicsStates(int cycle, double time) timesteps_.clear(); time_ = time; + dt_ = 0.0; max_time_ = time; min_time_ = time; cycle_ = cycle; diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index bfee79de35..a97e9b9d23 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -618,7 +618,10 @@ class BasePhysics { */ double time_; - double time_prev_; + /** + * @brief Current time step + */ + double dt_; /** * @brief The maximum time reached for the forward solver diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 97265ddd4b..2082fae1da 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -740,15 +740,13 @@ class SolidMechanics, std::integer_se template struct MaterialStressFunctor { /// @brief Constructor for the functor - MaterialStressFunctor(Material material, SolidMechanics& owner) : material_(material), owner_(&owner) {} + MaterialStressFunctor(Material material, const double* dt) : material_(material), time_increment_(dt) {} /// @brief Material model Material material_; - /// @brief Pointer to the `SolidMechanics` instance that owns the current `MaterialStressFunctor` - /// - /// This is needed to compute the time increment from the enclosing - const SolidMechanics* owner_; + /// @brief Current time step + const double* time_increment_; /** * @brief Material stress response call @@ -771,9 +769,8 @@ class SolidMechanics, std::integer_se { auto du_dX = get(displacement); auto d2u_dt2 = get(acceleration); - const double dt = owner_->time_ - owner_->time_prev_; - auto stress = material_(state, dt, du_dX, params...); + auto stress = material_(state, *time_increment_, du_dX, params...); return serac::tuple{material_.density * d2u_dt2, stress}; } @@ -810,7 +807,7 @@ class SolidMechanics, std::integer_se { static_assert(std::is_same_v || std::is_same_v, "invalid quadrature data provided in setMaterial()"); - MaterialStressFunctor material_functor(material, *this); + MaterialStressFunctor material_functor(material, &dt_); residual_->AddDomainIntegral( Dimension{}, DependsOn<0, 1, @@ -1173,10 +1170,12 @@ class SolidMechanics, std::integer_se } } + // store the current time step for use in rate-dependent materials + dt_ = dt; + if (is_quasistatic_) { quasiStaticSolve(dt); } else { - time_prev_ = time_; // The current ode interface tracks 2 times, one internally which we have a handle to via time_, // and one here via the step interface. // We are ignoring this one, and just using the internal version for now. From 7e9a253eca9376a4cb2d583f432db60af15e805b Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Mon, 13 Jan 2025 15:39:09 -0800 Subject: [PATCH 10/32] Avoid a catastrophic cancellation in the tensor eigensolve --- src/serac/numerics/functional/tuple_tensor_dual_functions.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/serac/numerics/functional/tuple_tensor_dual_functions.hpp b/src/serac/numerics/functional/tuple_tensor_dual_functions.hpp index e44a315974..69b73c4ddf 100644 --- a/src/serac/numerics/functional/tuple_tensor_dual_functions.hpp +++ b/src/serac/numerics/functional/tuple_tensor_dual_functions.hpp @@ -861,10 +861,10 @@ inline SERAC_HOST_DEVICE tuple eig_symm(const mat3& A) double A11 = dot(s0, A_dev_s0); double A12 = dot(s0, A_dev_s1); - double A21 = dot(s1, A_dev_s0); + double A21 = A12; double A22 = dot(s1, A_dev_s1); - double delta = 0.5 * sgn(A11 - A22) * std::sqrt((A11 - A22) * (A11 - A22) + 4 * A12 * A21); + double delta = 0.5 * std::sqrt((A11 - A22) * (A11 - A22) + 4 * A12 * A21); eta(1) = 0.5 * (A11 + A22) - delta; eta(2) = 0.5 * (A11 + A22) + delta; From 6eee57fc6a0da4cf6bc76a73896527a7560e1444 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Mon, 13 Jan 2025 15:40:27 -0800 Subject: [PATCH 11/32] Change where time increment is set for material to avoid incorrect stiffness matrix in warm start --- src/serac/physics/solid_mechanics.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 2082fae1da..fd33b9975b 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -1170,9 +1170,6 @@ class SolidMechanics, std::integer_se } } - // store the current time step for use in rate-dependent materials - dt_ = dt; - if (is_quasistatic_) { quasiStaticSolve(dt); } else { @@ -1183,6 +1180,7 @@ class SolidMechanics, std::integer_se // but at the moment, the double times creates a lot of confusion, so // we short circuit the extra time here by passing a dummy time and ignoring it. double time_tmp = time_; + dt_ = dt; ode2_.Step(displacement_, velocity_, time_tmp, dt); } @@ -1521,7 +1519,6 @@ class SolidMechanics, std::integer_se // warm start must be called prior to the time update so that the previous Jacobians can be used consistently // throughout. warmStartDisplacement(dt); - time_prev_ = time_; time_ += dt; // this method is essentially equivalent to the 1-liner @@ -1641,7 +1638,7 @@ class SolidMechanics, std::integer_se } if (use_warm_start_) { - // Update the linearized Jacobian matrix + // Update external forcing auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, *parameters_[parameter_indices].state...); @@ -1669,6 +1666,7 @@ class SolidMechanics, std::integer_se } displacement_ += du_; + dt_ = dt; } }; From 26022f2b707ee36b0c2372f36216de5a5ade689a Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Mon, 13 Jan 2025 16:44:59 -0800 Subject: [PATCH 12/32] Add more output options to uiaxial tension --- examples/uniaxial/uniaxial_tension.cpp | 32 +++++++++++++++++--------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/examples/uniaxial/uniaxial_tension.cpp b/examples/uniaxial/uniaxial_tension.cpp index 15de518ddb..3054876e66 100644 --- a/examples/uniaxial/uniaxial_tension.cpp +++ b/examples/uniaxial/uniaxial_tension.cpp @@ -21,6 +21,14 @@ #include +template +void output(double u, double f, const Physics& solid, const std::string& paraview_tag, std::ofstream& file) +{ + solid.outputStateToDisk(paraview_tag); + file << solid.time() << " " << u << " " << f << std::endl; +} + + int main(int argc, char* argv[]) { feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); @@ -49,6 +57,7 @@ int main(int argc, char* argv[]) constexpr double eta = 1e-2; constexpr double max_strain = 3*strain_constant; + std::string output_filename = "uniaxial_fd.txt"; // Handle command line arguments axom::CLI::App app{"Plane strain uniaxial extension of a bar."}; @@ -57,11 +66,13 @@ int main(int argc, char* argv[]) app.add_option("--parallel-refinements", parallel_refinements, "Parallel refinement steps", true); app.add_option("--time-steps", time_steps, "Number of time steps to divide simulation", true); app.add_option("--strain-rate", strain_rate, "Nominal strain rate", true); + app.add_option("--output-file", output_filename, "Name for force-displacement output file", true); app.set_help_flag("--help"); CLI11_PARSE(app, argc, argv); SLIC_INFO_ROOT(axom::fmt::format("strain rate: {}", strain_rate)); + SLIC_INFO_ROOT(axom::fmt::format("time_steps: {}", time_steps)); double max_time = max_strain/strain_rate; @@ -114,9 +125,6 @@ int main(int argc, char* argv[]) solid_solver.completeSetup(); - std::string paraview_tag = simulation_tag + "_paraview"; - solid_solver.outputStateToDisk(paraview_tag); - double dt = max_time/(time_steps - 1); // get nodes and dofs to compute total force @@ -131,24 +139,26 @@ int main(int argc, char* argv[]) return R; }; - std::ofstream file("uniaxial.txt"); + std::string paraview_tag = simulation_tag + "_paraview"; + std::ofstream file(output_filename); file << "# time displacement force" << std::endl; - file << 0 << " " << 0 << " " << 0 << std::endl; + { + double u = applied_displacement(serac::vec3{}, solid_solver.time())[0]; + double f = compute_net_force(solid_solver.dual("reactions")); + output(u, f, solid_solver, paraview_tag, file); + } - for (int i = 0; i < time_steps; ++i) { + for (int i = 1; i < time_steps; ++i) { SLIC_INFO_ROOT("------------------------------------------"); SLIC_INFO_ROOT(axom::fmt::format("TIME STEP {}", i)); - SLIC_INFO_ROOT(axom::fmt::format("time = {} (out of {})", solid_solver.time(), max_time)); + SLIC_INFO_ROOT(axom::fmt::format("time = {} (out of {})", solid_solver.time() + dt, max_time)); serac::logger::flush(); solid_solver.advanceTimestep(dt); - // Output the sidre-based plot files - solid_solver.outputStateToDisk(paraview_tag); - double u = applied_displacement(serac::vec3{}, solid_solver.time())[0]; double f = compute_net_force(solid_solver.dual("reactions")); - file << solid_solver.time() << " " << u << " " << f << std::endl; + output(u, f, solid_solver, paraview_tag, file); } file.close(); From 0f386068c1b792ce51f6768e63c72b450538b513 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Wed, 15 Jan 2025 17:18:41 -0800 Subject: [PATCH 13/32] Change direction: make a new material setter for rate dependent materials so rate independent material signatures are preserved --- examples/uniaxial/uniaxial_tension.cpp | 5 +- .../physics/materials/solid_material.hpp | 8 +- src/serac/physics/solid_mechanics.hpp | 80 +++++++++++++++++-- src/serac/physics/tests/solid.cpp | 4 +- .../physics/tests/solid_multi_material.cpp | 2 +- .../physics/tests/solid_statics_patch.cpp | 8 +- 6 files changed, 88 insertions(+), 19 deletions(-) diff --git a/examples/uniaxial/uniaxial_tension.cpp b/examples/uniaxial/uniaxial_tension.cpp index 3054876e66..ca0f3adc66 100644 --- a/examples/uniaxial/uniaxial_tension.cpp +++ b/examples/uniaxial/uniaxial_tension.cpp @@ -50,6 +50,7 @@ int main(int argc, char* argv[]) constexpr double E = 1.0; constexpr double nu = 0.25; + constexpr double density = 1.0; constexpr double sigma_y = 0.001; constexpr double sigma_sat = 3.0*sigma_y; @@ -111,11 +112,11 @@ int main(int argc, char* argv[]) using Material = serac::solid_mechanics::J2; Hardening hardening{sigma_y, sigma_sat, strain_constant, eta}; - Material material{E, nu, hardening, .density = 1.0}; + Material material{E, nu, hardening, density}; auto internal_states = solid_solver.createQuadratureDataBuffer(Material::State{}, whole_mesh); - solid_solver.setMaterial(material, whole_mesh, internal_states); + solid_solver.setRateDependentMaterial(material, whole_mesh, internal_states); solid_solver.setFixedBCs(fix_x, serac::Component::X); solid_solver.setFixedBCs(fix_y, serac::Component::Y); diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index ac201460e6..ba39b95189 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -36,7 +36,7 @@ struct LinearIsotropic { * @return The Cauchy stress */ template - SERAC_HOST_DEVICE auto operator()(State& /* state */, double /* dt */, const tensor& du_dX) const + SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor& du_dX) const { auto I = Identity(); auto lambda = K - (2.0 / 3.0) * G; @@ -72,7 +72,7 @@ struct StVenantKirchhoff { * @return The Cauchy stress */ template - auto operator()(State&, double /* dt */, const tensor& grad_u) const + auto operator()(State&, const tensor& grad_u) const { static constexpr auto I = Identity(); auto F = grad_u + I; @@ -107,7 +107,7 @@ struct NeoHookean { * @return The Cauchy stress */ template - SERAC_HOST_DEVICE auto operator()(State& /* state */, double /* dt */, const tensor& du_dX) const + SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor& du_dX) const { using std::log1p; constexpr auto I = Identity(); @@ -140,7 +140,7 @@ struct NeoHookeanAdditiveSplit { * @return The first Piola stress */ template - SERAC_HOST_DEVICE auto operator()(State& /* state */, double /* dt */, const tensor& du_dX) const + SERAC_HOST_DEVICE auto operator()(State& /* state */, const tensor& du_dX) const { using std::pow; constexpr auto I = Identity(); diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index fd33b9975b..a1212bb20d 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -740,14 +740,11 @@ class SolidMechanics, std::integer_se template struct MaterialStressFunctor { /// @brief Constructor for the functor - MaterialStressFunctor(Material material, const double* dt) : material_(material), time_increment_(dt) {} + MaterialStressFunctor(Material material) : material_(material) {} /// @brief Material model Material material_; - /// @brief Current time step - const double* time_increment_; - /** * @brief Material stress response call * @@ -770,7 +767,7 @@ class SolidMechanics, std::integer_se auto du_dX = get(displacement); auto d2u_dt2 = get(acceleration); - auto stress = material_(state, *time_increment_, du_dX, params...); + auto stress = material_(state, du_dX, params...); return serac::tuple{material_.density * d2u_dt2, stress}; } @@ -807,7 +804,7 @@ class SolidMechanics, std::integer_se { static_assert(std::is_same_v || std::is_same_v, "invalid quadrature data provided in setMaterial()"); - MaterialStressFunctor material_functor(material, &dt_); + MaterialStressFunctor material_functor(material); residual_->AddDomainIntegral( Dimension{}, DependsOn<0, 1, @@ -826,6 +823,77 @@ class SolidMechanics, std::integer_se setMaterial(DependsOn<>{}, material, domain, qdata); } + /** + * Functor for materials that get dt as an argument + */ + template + struct RateDependentMaterialStressFunctor { + /// @brief Constructor for the functor + RateDependentMaterialStressFunctor(Material material, const double* dt) : material_(material), time_increment_(dt) {} + + /// @brief Material model + Material material_; + + /// @brief Current time step + const double* time_increment_; + + /** + * @brief Material stress response call + * + * @tparam X Spatial position type + * @tparam State state + * @tparam Displacement displacement + * @tparam Acceleration acceleration + * @tparam Params variadic parameters for call + * @param[in] state state + * @param[in] displacement displacement + * @param[in] acceleration acceleration + * @param[in] params parameter pack + * @return The calculated material response (tuple of volumetric heat capacity and thermal flux) for a linear + * isotropic material + */ + template + auto SERAC_HOST_DEVICE operator()(double, X, State& state, Displacement displacement, Acceleration acceleration, + Params... params) const + { + auto du_dX = get(displacement); + auto d2u_dt2 = get(acceleration); + + auto stress = material_(state, *time_increment_, du_dX, params...); + + return serac::tuple{material_.density * d2u_dt2, stress}; + } + }; + + + /** + * Set a material that gets dt as an argument + */ + template + void setRateDependentMaterial(DependsOn, const RateDependentMaterialType& material, Domain& domain, + qdata_type qdata = EmptyQData) + { + static_assert(std::is_same_v || std::is_same_v, + "invalid quadrature data provided in setMaterial()"); + RateDependentMaterialStressFunctor material_functor(material, &dt_); + residual_->AddDomainIntegral( + Dimension{}, + DependsOn<0, 1, + active_parameters + NUM_STATE_VARS...>{}, // the magic number "+ NUM_STATE_VARS" accounts for the + // fact that the displacement, acceleration, and shape + // fields are always-on and come first, so the `n`th + // parameter will actually be argument `n + NUM_STATE_VARS` + std::move(material_functor), domain, qdata); + } + + /// @overload + template + void setRateDependentMaterial(const RateDependentMaterialType& material, Domain& domain, + std::shared_ptr> qdata = EmptyQData) + { + setRateDependentMaterial(DependsOn<>{}, material, domain, qdata); + } + /** * @brief Set the underlying finite element state to a prescribed displacement * diff --git a/src/serac/physics/tests/solid.cpp b/src/serac/physics/tests/solid.cpp index 98ab94d9c3..48f1d667bc 100644 --- a/src/serac/physics/tests/solid.cpp +++ b/src/serac/physics/tests/solid.cpp @@ -66,7 +66,7 @@ void functional_solid_test_static_J2() using Hardening = solid_mechanics::LinearHardening; using Material = solid_mechanics::J2SmallStrain; - Hardening hardening{.sigma_y = 50.0, .Hi = 50.0}; + Hardening hardening{.sigma_y = 50.0, .Hi = 50.0, .eta = 0.0}; Material mat{ .E = 10000, // Young's modulus .nu = 0.25, // Poisson's ratio @@ -80,7 +80,7 @@ void functional_solid_test_static_J2() auto qdata = solid_solver.createQuadratureDataBuffer(initial_state); Domain whole_domain = EntireDomain(pmesh); - solid_solver.setMaterial(mat, whole_domain, qdata); + solid_solver.setRateDependentMaterial(mat, whole_domain, qdata); // prescribe zero displacement at the supported end of the beam, auto support = Domain::ofBoundaryElements(pmesh, by_attr(1)); diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index 6b56fbbff1..202174cb30 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -217,7 +217,7 @@ TEST(Solid, MultiMaterialWithState) auto qdata = solid.createQuadratureDataBuffer(initial_state, right); solid.setMaterial(mat_left, left); - solid.setMaterial(mat_right, right, qdata); + solid.setRateDependentMaterial(mat_right, right, qdata); Domain end_face = Domain::ofBoundaryElements(pmesh, by_attr(3)); solid.setTraction( diff --git a/src/serac/physics/tests/solid_statics_patch.cpp b/src/serac/physics/tests/solid_statics_patch.cpp index 7b6906c4db..99a42944dd 100644 --- a/src/serac/physics/tests/solid_statics_patch.cpp +++ b/src/serac/physics/tests/solid_statics_patch.cpp @@ -77,7 +77,7 @@ class AffineSolution { * @param essential_boundaries Boundary attributes on which essential boundary conditions are desired */ template - void applyLoads(const Material& material, SolidMechanics& sf, Domain essential_boundary, double dt) const + void applyLoads(const Material& material, SolidMechanics& sf, Domain essential_boundary) const { // essential BCs auto ebc_func = [*this](tensor X, double) { return this->eval(X); }; @@ -86,7 +86,7 @@ class AffineSolution { // natural BCs typename Material::State state; - tensor P = material(state, dt, A); + tensor P = material(state, A); auto traction = [P](auto, auto n0, auto) { return dot(P, n0); }; Domain entire_boundary = EntireBoundary(sf.mesh()); sf.setTraction(traction, entire_boundary); @@ -219,7 +219,7 @@ double solution_error(PatchBoundaryCondition bc) constexpr double dt = 1.0; Domain essential_boundary = Domain::ofBoundaryElements(pmesh, by_attr(essentialBoundaryAttributes(bc))); - exact_displacement.applyLoads(mat, solid, essential_boundary, dt); + exact_displacement.applyLoads(mat, solid, essential_boundary); // Finalize the data structures solid.completeSetup(); @@ -342,7 +342,7 @@ double pressure_error() constexpr double dt = 1.0; - tensor P = mat(state, dt, H); + tensor P = mat(state, H); auto F = H + Identity(); auto sigma = dot(P, transpose(F)) / det(F); double pressure = -sigma[0][0]; From 7f76133c9aac03664ee9ca380e6e7e6c0e413593 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 07:35:34 -0800 Subject: [PATCH 14/32] Fix build error by making material point simulator for rate dependent materials --- .../materials/material_verification_tools.hpp | 21 +++++++++++++++++-- .../physics/materials/tests/J2_material.cpp | 8 +++---- .../physics/tests/solid_statics_patch.cpp | 1 + 3 files changed, 24 insertions(+), 6 deletions(-) diff --git a/src/serac/physics/materials/material_verification_tools.hpp b/src/serac/physics/materials/material_verification_tools.hpp index 45fe47e415..30a81ff928 100644 --- a/src/serac/physics/materials/material_verification_tools.hpp +++ b/src/serac/physics/materials/material_verification_tools.hpp @@ -64,7 +64,6 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat output_history.reserve(num_steps); tensor dudx{}; - const double dt = t_max / double(num_steps - 1); for (size_t i = 0; i < num_steps; i++) { auto initial_guess = tensor{dudx[1][1], dudx[2][2]}; auto epsilon_yy_and_zz = find_root(sigma_yy_and_zz, initial_guess); @@ -72,7 +71,7 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat dudx[1][1] = epsilon_yy_and_zz[0]; dudx[2][2] = epsilon_yy_and_zz[1]; - auto stress = material(state, dt, dudx, parameter_functions(t)...); + auto stress = material(state, dudx, parameter_functions(t)...); output_history.push_back(tuple{t, dudx, stress, state}); t += dt; @@ -81,6 +80,24 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat return output_history; } +/** + * @brief Drive a rate-dependent material model thorugh a uniaxial tension experiment + * + * Drives material model through specified axial displacement gradient history. + * The time elaspses from 0 up to t_max. + * Currently only implemented for isotropic materials (or orthotropic materials with the + * principal axes aligned with the coordinate directions). + */ +template +auto uniaxial_stress_test_rate_dependent(double t_max, size_t num_steps, const MaterialType material, const StateType initial_state, + std::function epsilon_xx, const parameter_types... parameter_functions) +{ + const double dt = t_max / double(num_steps - 1); + auto mat_with_dt = [&material, dt](auto& state, auto du_dx, parameter_types... parameters) { return material(state, dt, du_dx, parameters...); }; + return uniaxial_stress_test(t_max, num_steps, mat_with_dt, initial_state, epsilon_xx, parameter_functions...); + +} + /** * @brief This function takes a material model (and associate state variables), * subjects it to a time history of stimuli, described by `functions ... f`, diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index 7d5c53acd3..24edb4e644 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -168,7 +168,7 @@ TEST(J2SmallStrain, Uniaxial) auto internal_state = Material::State{}; auto strain = [=](double t) { return sigma_y / E * t; }; - auto response_history = uniaxial_stress_test(2.0, 4, material, internal_state, strain); + auto response_history = uniaxial_stress_test_rate_dependent(2.0, 4, material, internal_state, strain); auto stress_exact = [=](double epsilon) { return epsilon < sigma_y / E ? E * epsilon : E / (E + Hi) * (sigma_y + Hi * epsilon); @@ -205,7 +205,7 @@ TEST(J2, Uniaxial) auto internal_state = Material::State{}; auto strain = [=](double t) { return sigma_y / E * t; }; - auto response_history = uniaxial_stress_test(2.0, 4, material, internal_state, strain); + auto response_history = uniaxial_stress_test_rate_dependent(2.0, 4, material, internal_state, strain); auto stress_exact = [=](double epsilon) { return epsilon < sigma_y / E ? E * epsilon : E / (E + Hi) * (sigma_y + Hi * epsilon); @@ -368,7 +368,7 @@ TEST(J2, PlotOutput) auto internal_state = Material::State{}; auto strain = [=](double t) { return max_strain*t/t_max; }; - auto response_history = uniaxial_stress_test(t_max, 100, material, internal_state, strain); + auto response_history = uniaxial_stress_test_rate_dependent(t_max, 100, material, internal_state, strain); std::ofstream file("J2_" + tag + ".txt"); @@ -378,7 +378,7 @@ TEST(J2, PlotOutput) double e = std::log1p(dudx[0][0]); // log strain double s = TK[0][0]; // Kirchhoff stress double pe = -std::log(get<3>(r).Fpinv[0][0]); // plastic strain - file << t << " " << e << " " << s << std::endl; + file << t << " " << e << " " << s << " " << pe << std::endl; } file.close(); }; diff --git a/src/serac/physics/tests/solid_statics_patch.cpp b/src/serac/physics/tests/solid_statics_patch.cpp index 9f45b90a5e..b7d7f7d1b4 100644 --- a/src/serac/physics/tests/solid_statics_patch.cpp +++ b/src/serac/physics/tests/solid_statics_patch.cpp @@ -360,6 +360,7 @@ double pressure_error() solid.completeSetup(); // Perform the quasi-static solve + constexpr double dt = 1.0; solid.advanceTimestep(dt); solid.outputStateToDisk(); From 235657949506953a0ee0f03f15fa44f6db9d4887 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 07:46:31 -0800 Subject: [PATCH 15/32] Clarify docstring --- src/serac/physics/solid_mechanics.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index d8761d0071..e295bcdfdb 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -820,7 +820,10 @@ class SolidMechanics, std::integer_se } /** - * Functor for materials that get dt as an argument + * Functor for materials that get dt as an argument + * + * Wraps materials that have the signature + * stress = material(state, dt, du_dX, params...) */ template struct RateDependentMaterialStressFunctor { From 8592dc96f522bbb5baa8b54a1ec24722b5e00937 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 13:24:33 -0800 Subject: [PATCH 16/32] Fix build warning for uninitialized variable - no effect on execution --- src/serac/physics/tests/solid_multi_material.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index 17220eee73..9aa914c805 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -198,7 +198,7 @@ TEST(Solid, MultiMaterialWithState) constexpr double Hi = E_right / 3.6; constexpr double sigma_y = 0.75 * applied_stress; - Hardening hardening{.sigma_y = sigma_y, .Hi = Hi}; + Hardening hardening{.sigma_y = sigma_y, .Hi = Hi, .eta = 0}; MaterialRight mat_right{ .E = E_right, // Young's modulus .nu = nu_right, // Poisson's ratio From b2ba2adc97366464be9c2e203590b03f53ff5f00 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 13:27:15 -0800 Subject: [PATCH 17/32] Remove test case that wasn't a real test, only executed to produce visual output --- .../physics/materials/tests/J2_material.cpp | 43 ------------------- 1 file changed, 43 deletions(-) diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index 24edb4e644..6d68347f5f 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -340,49 +340,6 @@ TEST(J2, FrameIndifference) ASSERT_LT(norm(error), 1e-13*norm(internal_state.Fpinv)); } -TEST(J2, PlotOutput) -{ - /* Log strain J2 plasticity has the nice feature that the exact uniaxial stress solution from - small strain plasticity are applicable, if you replace the lineasr strain with log strain - and use the Kirchhoff stress as the output. - */ - - using Hardening = solid_mechanics::VoceHardening; - using Material = solid_mechanics::J2; - - constexpr double E = 1.0; - constexpr double sigma_y = 0.001; - constexpr double sigma_sat = 3*sigma_y; - constexpr double strain_constant = 10*sigma_y/E; - constexpr double eta = 1e-2; - - constexpr double max_strain = 3*strain_constant; - constexpr double strain_rate = 1e-1; - const std::string tag = "m1"; - constexpr double t_max = max_strain/strain_rate; - - - - Hardening hardening{sigma_y, sigma_sat, strain_constant, eta}; - Material material{.E = E, .nu = 0.25, .hardening = hardening, .density = 1.0}; - - auto internal_state = Material::State{}; - auto strain = [=](double t) { return max_strain*t/t_max; }; - auto response_history = uniaxial_stress_test_rate_dependent(t_max, 100, material, internal_state, strain); - - std::ofstream file("J2_" + tag + ".txt"); - - for (auto r : response_history) { - auto [t, dudx, P, state] = r; - auto TK = dot(P, transpose(dudx + Identity<3>())); - double e = std::log1p(dudx[0][0]); // log strain - double s = TK[0][0]; // Kirchhoff stress - double pe = -std::log(get<3>(r).Fpinv[0][0]); // plastic strain - file << t << " " << e << " " << s << " " << pe << std::endl; - } - file.close(); -}; - } // namespace serac int main(int argc, char* argv[]) From 883edba42cf889953eccaa204a5a6c878e195bce Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 13:29:03 -0800 Subject: [PATCH 18/32] Make style --- examples/uniaxial/uniaxial_tension.cpp | 39 +++++++++---------- .../materials/material_verification_tools.hpp | 10 +++-- .../physics/materials/solid_material.hpp | 11 +++--- src/serac/physics/solid_mechanics.hpp | 18 +++++---- 4 files changed, 41 insertions(+), 37 deletions(-) diff --git a/examples/uniaxial/uniaxial_tension.cpp b/examples/uniaxial/uniaxial_tension.cpp index ca0f3adc66..5f28bb558b 100644 --- a/examples/uniaxial/uniaxial_tension.cpp +++ b/examples/uniaxial/uniaxial_tension.cpp @@ -28,7 +28,6 @@ void output(double u, double f, const Physics& solid, const std::string& paravie file << solid.time() << " " << u << " " << f << std::endl; } - int main(int argc, char* argv[]) { feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); @@ -42,7 +41,7 @@ int main(int argc, char* argv[]) constexpr int elements_in_x = 10; constexpr int elements_in_y = 1; constexpr int elements_in_z = 1; - + int serial_refinements = 0; int parallel_refinements = 0; int time_steps = 100; @@ -53,11 +52,11 @@ int main(int argc, char* argv[]) constexpr double density = 1.0; constexpr double sigma_y = 0.001; - constexpr double sigma_sat = 3.0*sigma_y; - constexpr double strain_constant = 10*sigma_y/E; + constexpr double sigma_sat = 3.0 * sigma_y; + constexpr double strain_constant = 10 * sigma_y / E; constexpr double eta = 1e-2; - constexpr double max_strain = 3*strain_constant; + constexpr double max_strain = 3 * strain_constant; std::string output_filename = "uniaxial_fd.txt"; // Handle command line arguments @@ -69,24 +68,23 @@ int main(int argc, char* argv[]) app.add_option("--strain-rate", strain_rate, "Nominal strain rate", true); app.add_option("--output-file", output_filename, "Name for force-displacement output file", true); app.set_help_flag("--help"); - + CLI11_PARSE(app, argc, argv); SLIC_INFO_ROOT(axom::fmt::format("strain rate: {}", strain_rate)); SLIC_INFO_ROOT(axom::fmt::format("time_steps: {}", time_steps)); - double max_time = max_strain/strain_rate; - + double max_time = max_strain / strain_rate; + // Create DataStore const std::string simulation_tag = "uniaxial"; const std::string mesh_tag = simulation_tag + "mesh"; axom::sidre::DataStore datastore; serac::StateManager::initialize(datastore, simulation_tag + "_data"); - auto mesh = serac::mesh::refineAndDistribute( - serac::buildCuboidMesh(elements_in_x, elements_in_y, elements_in_z, - x_length, y_length, z_length), - serial_refinements, parallel_refinements); + auto mesh = serac::mesh::refineAndDistribute( + serac::buildCuboidMesh(elements_in_x, elements_in_y, elements_in_z, x_length, y_length, z_length), + serial_refinements, parallel_refinements); auto& pmesh = serac::StateManager::setMesh(std::move(mesh), mesh_tag); // create boundary domains for boundary conditions @@ -98,15 +96,14 @@ int main(int argc, char* argv[]) serac::LinearSolverOptions linear_options{.linear_solver = serac::LinearSolver::Strumpack, .print_level = 0}; - serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton, - .relative_tol = 1.0e-10, - .absolute_tol = 1.0e-12, + serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton, + .relative_tol = 1.0e-10, + .absolute_tol = 1.0e-12, .max_iterations = 200, - .print_level = 1}; + .print_level = 1}; serac::SolidMechanics solid_solver( - nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, simulation_tag, - mesh_tag); + nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, simulation_tag, mesh_tag); using Hardening = serac::solid_mechanics::VoceHardening; using Material = serac::solid_mechanics::J2; @@ -121,12 +118,14 @@ int main(int argc, char* argv[]) solid_solver.setFixedBCs(fix_x, serac::Component::X); solid_solver.setFixedBCs(fix_y, serac::Component::Y); solid_solver.setFixedBCs(fix_z, serac::Component::Z); - auto applied_displacement = [strain_rate](serac::vec3, double t){ return serac::vec3{strain_rate*x_length*t, 0., 0.}; }; + auto applied_displacement = [strain_rate](serac::vec3, double t) { + return serac::vec3{strain_rate * x_length * t, 0., 0.}; + }; solid_solver.setDisplacementBCs(applied_displacement, apply_displacement, serac::Component::X); solid_solver.completeSetup(); - double dt = max_time/(time_steps - 1); + double dt = max_time / (time_steps - 1); // get nodes and dofs to compute total force mfem::Array dof_list = apply_displacement.dof_list(&solid_solver.displacement().space()); diff --git a/src/serac/physics/materials/material_verification_tools.hpp b/src/serac/physics/materials/material_verification_tools.hpp index 30a81ff928..9aedcc8fa7 100644 --- a/src/serac/physics/materials/material_verification_tools.hpp +++ b/src/serac/physics/materials/material_verification_tools.hpp @@ -89,13 +89,15 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat * principal axes aligned with the coordinate directions). */ template -auto uniaxial_stress_test_rate_dependent(double t_max, size_t num_steps, const MaterialType material, const StateType initial_state, - std::function epsilon_xx, const parameter_types... parameter_functions) +auto uniaxial_stress_test_rate_dependent(double t_max, size_t num_steps, const MaterialType material, + const StateType initial_state, std::function epsilon_xx, + const parameter_types... parameter_functions) { const double dt = t_max / double(num_steps - 1); - auto mat_with_dt = [&material, dt](auto& state, auto du_dx, parameter_types... parameters) { return material(state, dt, du_dx, parameters...); }; + auto mat_with_dt = [&material, dt](auto& state, auto du_dx, parameter_types... parameters) { + return material(state, dt, du_dx, parameters...); + }; return uniaxial_stress_test(t_max, num_steps, mat_with_dt, initial_state, epsilon_xx, parameter_functions...); - } /** diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index 2d164e15aa..ce2dc31c32 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -166,7 +166,7 @@ struct NeoHookeanAdditiveSplit { template auto overstress(double viscosity, T accumulated_plastic_strain_rate) { - return viscosity*accumulated_plastic_strain_rate; + return viscosity * accumulated_plastic_strain_rate; } /** @@ -211,7 +211,8 @@ struct PowerLawHardening { auto operator()(T1 accumulated_plastic_strain, T2 accumulated_plastic_strain_rate) const { using std::pow; - return sigma_y * pow(1.0 + accumulated_plastic_strain / eps0, 1.0 / n) + overstress(eta, accumulated_plastic_strain_rate); + return sigma_y * pow(1.0 + accumulated_plastic_strain / eps0, 1.0 / n) + + overstress(eta, accumulated_plastic_strain_rate); }; }; @@ -224,7 +225,7 @@ struct VoceHardening { double sigma_y; ///< yield strength double sigma_sat; ///< saturation value of flow strength double strain_constant; ///< The constant dictating how fast the exponential decays - double eta = 0.0; ///< linear viscosity for rate dependence + double eta = 0.0; ///< linear viscosity for rate dependence /** * @brief Computes the flow stress @@ -281,7 +282,7 @@ struct J2SmallStrain { // (ii) admissibility const double eqps_old = state.accumulated_plastic_strain; auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_q) { - auto eqps_dot = dt > 0.0? delta_eqps/dt : 0.0*delta_eqps; + auto eqps_dot = dt > 0.0 ? delta_eqps / dt : 0.0 * delta_eqps; return trial_q - (3.0 * G + Hk) * delta_eqps - this->hardening(eqps_old + delta_eqps, eqps_dot); }; if (residual(0.0, get_value(q)) > tol * hardening.sigma_y) { @@ -345,7 +346,7 @@ struct J2 { // (ii) admissibility const double eqps_old = state.accumulated_plastic_strain; auto residual = [eqps_old, G, dt, *this](auto delta_eqps, auto trial_mises) { - auto eqps_dot = dt > 0.0? delta_eqps/dt : 0.0*delta_eqps; + auto eqps_dot = dt > 0.0 ? delta_eqps / dt : 0.0 * delta_eqps; return trial_mises - 3.0 * G * delta_eqps - this->hardening(eqps_old + delta_eqps, eqps_dot); }; if (residual(0.0, q_val) > tol * hardening.sigma_y) { diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index c416b79150..26c9589604 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -828,7 +828,9 @@ class SolidMechanics, std::integer_se template struct RateDependentMaterialStressFunctor { /// @brief Constructor for the functor - RateDependentMaterialStressFunctor(Material material, const double* dt) : material_(material), time_increment_(dt) {} + RateDependentMaterialStressFunctor(Material material, const double* dt) : material_(material), time_increment_(dt) + { + } /// @brief Material model Material material_; @@ -855,7 +857,7 @@ class SolidMechanics, std::integer_se auto SERAC_HOST_DEVICE operator()(double, X, State& state, Displacement displacement, Acceleration acceleration, Params... params) const { - auto du_dX = get(displacement); + auto du_dX = get(displacement); auto d2u_dt2 = get(acceleration); auto stress = material_(state, *time_increment_, du_dX, params...); @@ -864,16 +866,16 @@ class SolidMechanics, std::integer_se } }; - /** * Set a material that gets dt as an argument */ template - void setRateDependentMaterial(DependsOn, const RateDependentMaterialType& material, Domain& domain, - qdata_type qdata = EmptyQData) + void setRateDependentMaterial(DependsOn, const RateDependentMaterialType& material, + Domain& domain, qdata_type qdata = EmptyQData) { - static_assert(std::is_same_v || std::is_same_v, - "invalid quadrature data provided in setMaterial()"); + static_assert( + std::is_same_v || std::is_same_v, + "invalid quadrature data provided in setMaterial()"); RateDependentMaterialStressFunctor material_functor(material, &dt_); residual_->AddDomainIntegral( Dimension{}, @@ -888,7 +890,7 @@ class SolidMechanics, std::integer_se /// @overload template void setRateDependentMaterial(const RateDependentMaterialType& material, Domain& domain, - std::shared_ptr> qdata = EmptyQData) + std::shared_ptr> qdata = EmptyQData) { setRateDependentMaterial(DependsOn<>{}, material, domain, qdata); } From b254304315e652bb309e16cedd0c5c89019db6fd Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 13:38:22 -0800 Subject: [PATCH 19/32] Fix doxygen warnings --- src/serac/physics/materials/solid_material.hpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index ce2dc31c32..33112b0495 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -163,6 +163,9 @@ struct NeoHookeanAdditiveSplit { double G; ///< shear modulus }; +/** + * Overstress for linear viscosity + */ template auto overstress(double viscosity, T accumulated_plastic_strain_rate) { @@ -175,7 +178,7 @@ auto overstress(double viscosity, T accumulated_plastic_strain_rate) struct LinearHardening { double sigma_y; ///< yield strength double Hi; ///< Isotropic hardening modulus - double eta; + double eta; ///< viscosity for linear rate sensitivity /** * @brief Computes the flow stress @@ -198,7 +201,7 @@ struct PowerLawHardening { double sigma_y; ///< yield strength double n; ///< hardening index in reciprocal form double eps0; ///< reference value of accumulated plastic strain - double eta = 0.0; + double eta = 0.0; ///< viscosity for linear rate sensitivity /** * @brief Computes the flow stress @@ -225,13 +228,15 @@ struct VoceHardening { double sigma_y; ///< yield strength double sigma_sat; ///< saturation value of flow strength double strain_constant; ///< The constant dictating how fast the exponential decays - double eta = 0.0; ///< linear viscosity for rate dependence + double eta = 0.0; ///< viscosity for linear rate sensitivity /** * @brief Computes the flow stress * - * @tparam T Number-like type for the argument - * @param accumulated_plastic_strain The uniaxial equivalent accumulated plastic strain + * @tparam T1 Number-like type + * @tparam T2 Number-like type + * @param epsilon_p The uniaxial equivalent accumulated plastic strain + * @param epsilon_p_dot The rate of the uniaxial equivalent accumulated plastic strain * @return Flow stress value */ template From d9bfc95099528d208f744139e1801d80ff01c4e6 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 13:38:55 -0800 Subject: [PATCH 20/32] Make style --- src/serac/physics/materials/solid_material.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index 33112b0495..c14883d27b 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -198,10 +198,10 @@ struct LinearHardening { * @brief Power-law isotropic hardening law */ struct PowerLawHardening { - double sigma_y; ///< yield strength - double n; ///< hardening index in reciprocal form - double eps0; ///< reference value of accumulated plastic strain - double eta = 0.0; ///< viscosity for linear rate sensitivity + double sigma_y; ///< yield strength + double n; ///< hardening index in reciprocal form + double eps0; ///< reference value of accumulated plastic strain + double eta = 0.0; ///< viscosity for linear rate sensitivity /** * @brief Computes the flow stress From 59193e03d379833223e234f8f7faaa4806347451 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 18 Jan 2025 13:47:49 -0800 Subject: [PATCH 21/32] fix more doxygen warnings --- src/serac/physics/materials/solid_material.hpp | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index c14883d27b..d5c0b5c2b5 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -183,8 +183,10 @@ struct LinearHardening { /** * @brief Computes the flow stress * - * @tparam T Number-like type for the argument + * @tparam T1 Number-like type + * @tparam T2 Number-like type * @param accumulated_plastic_strain The uniaxial equivalent accumulated plastic strain + * @param accumulated_plastic_strain_rate The rate of the uniaxial equivalent accumulated plastic strain * @return Flow stress value */ template @@ -206,8 +208,10 @@ struct PowerLawHardening { /** * @brief Computes the flow stress * - * @tparam T Number-like type for the argument + * @tparam T1 Number-like type + * @tparam T2 Number-like type * @param accumulated_plastic_strain The uniaxial equivalent accumulated plastic strain + * @param accumulated_plastic_strain_rate The rate of the uniaxial equivalent accumulated plastic strain * @return Flow stress value */ template @@ -235,15 +239,17 @@ struct VoceHardening { * * @tparam T1 Number-like type * @tparam T2 Number-like type - * @param epsilon_p The uniaxial equivalent accumulated plastic strain - * @param epsilon_p_dot The rate of the uniaxial equivalent accumulated plastic strain + * @param accumulated_plastic_strain The uniaxial equivalent accumulated plastic strain + * @param accumulated_plastic_strain_rate The rate of the uniaxial equivalent accumulated plastic strain * @return Flow stress value */ template - auto operator()(T1 epsilon_p, T2 epsilon_p_dot) const + auto operator()(T1 accumulated_plastic_strain, T2 accumulated_plastic_strain_rate) const { using std::exp; + auto& epsilon_p = accumulated_plastic_strain; auto Y_eq = sigma_sat - (sigma_sat - sigma_y) * exp(-epsilon_p / strain_constant); + auto& epsilon_p_dot = accumulated_plastic_strain_rate; auto Y_vis = overstress(eta, epsilon_p_dot); return Y_eq + Y_vis; }; From 33d8b9d31c1be974ee8b92b9a061e833686d887b Mon Sep 17 00:00:00 2001 From: Brandon Talamini <30813018+btalamini@users.noreply.github.com> Date: Mon, 27 Jan 2025 15:40:45 -0800 Subject: [PATCH 22/32] Update examples/uniaxial/CMakeLists.txt Co-authored-by: Chris White --- examples/uniaxial/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/uniaxial/CMakeLists.txt b/examples/uniaxial/CMakeLists.txt index a0195ddae9..1fb6e34750 100644 --- a/examples/uniaxial/CMakeLists.txt +++ b/examples/uniaxial/CMakeLists.txt @@ -10,7 +10,6 @@ blt_add_executable( NAME uniaxial_tension DEPENDS_ON serac_physics serac_mesh ) -# add serac_functional to DEPENDS_ON install( FILES From 789122b076ba5d2de93552fe6b4b5511c9915f4e Mon Sep 17 00:00:00 2001 From: Brandon Talamini <30813018+btalamini@users.noreply.github.com> Date: Mon, 27 Jan 2025 15:41:13 -0800 Subject: [PATCH 23/32] Fix typ in comment Co-authored-by: Chris White --- src/serac/physics/materials/material_verification_tools.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serac/physics/materials/material_verification_tools.hpp b/src/serac/physics/materials/material_verification_tools.hpp index 9aedcc8fa7..ae0290e07d 100644 --- a/src/serac/physics/materials/material_verification_tools.hpp +++ b/src/serac/physics/materials/material_verification_tools.hpp @@ -84,7 +84,7 @@ auto uniaxial_stress_test(double t_max, size_t num_steps, const MaterialType mat * @brief Drive a rate-dependent material model thorugh a uniaxial tension experiment * * Drives material model through specified axial displacement gradient history. - * The time elaspses from 0 up to t_max. + * The time elapses from 0 up to t_max. * Currently only implemented for isotropic materials (or orthotropic materials with the * principal axes aligned with the coordinate directions). */ From 88b7181922f670e1143fa807113cca2c0664c8bf Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Wed, 5 Feb 2025 14:17:06 -0800 Subject: [PATCH 24/32] Remove default values on some hardening classes so that they're all consistent --- src/serac/physics/materials/solid_material.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index d5c0b5c2b5..b83cfe0816 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -203,7 +203,7 @@ struct PowerLawHardening { double sigma_y; ///< yield strength double n; ///< hardening index in reciprocal form double eps0; ///< reference value of accumulated plastic strain - double eta = 0.0; ///< viscosity for linear rate sensitivity + double eta; ///< viscosity for linear rate sensitivity /** * @brief Computes the flow stress @@ -232,7 +232,7 @@ struct VoceHardening { double sigma_y; ///< yield strength double sigma_sat; ///< saturation value of flow strength double strain_constant; ///< The constant dictating how fast the exponential decays - double eta = 0.0; ///< viscosity for linear rate sensitivity + double eta; ///< viscosity for linear rate sensitivity /** * @brief Computes the flow stress From cf700d53a34c36912c809d772002f8b08d83d4b5 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Thu, 6 Feb 2025 09:55:23 -0800 Subject: [PATCH 25/32] Explicitly initialize viscosity parameter in all tests --- src/serac/physics/materials/tests/J2_material.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index 6d68347f5f..9ffe7690ad 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -236,7 +236,7 @@ TEST(J2, DerivativeCorrectness) using Hardening = solid_mechanics::PowerLawHardening; using Material = solid_mechanics::J2; - Hardening hardening{.sigma_y = 350e6, .n = 3, .eps0 = 0.00175}; + Hardening hardening{.sigma_y = 350e6, .n = 3, .eps0 = 0.00175, .eta = 0.0}; Material material{.E = 200e9, .nu = 0.25, .hardening = hardening, .density = 1.0}; // initialize internal state variables @@ -287,7 +287,7 @@ TEST(J2, FrameIndifference) using Hardening = solid_mechanics::VoceHardening; using Material = solid_mechanics::J2; - Hardening hardening{.sigma_y = 350e6, .sigma_sat = 700e6, .strain_constant = 0.01}; + Hardening hardening{.sigma_y = 350e6, .sigma_sat = 700e6, .strain_constant = 0.01, .eta = 0.0}; Material material{.E = 200.0e9, .nu = 0.25, .hardening = hardening, .density = 1.0}; // clang-format off From c440bb00f0de74508b84b111ab750690597e896d Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 8 Feb 2025 07:15:39 -0800 Subject: [PATCH 26/32] Add rigorous test of rate dependent case for J2 material --- .../physics/materials/tests/J2_material.cpp | 123 +++++++++++++++++- 1 file changed, 122 insertions(+), 1 deletion(-) diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index 9ffe7690ad..5f4045fd4e 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -21,7 +21,42 @@ namespace serac { -using namespace serac; +template +tensor log(tensor v) +{ + using std::log; + return make_tensor([&v](int i) { return log(v[i]); }); +} + +template +double mean(tensor v) +{ + double sum = 0; + for (int i = 0; i < n; i++) sum += v[i]; + return sum/n; +} + +template +double best_fit_slope(tensor x, tensor y) +{ + double x_bar = mean(x); + double y_bar = mean(y); + double numer = 0; + double denom = 0; + for (int i = 0; i < n; i++) { + numer += (x[i] - x_bar)*(y[i] - y_bar); + denom += (x[i] - x_bar)*(x[i] - x_bar); + } + return numer/denom; +} + +template +double compute_convergence_rate(tensor h, tensor e) +{ + auto log_h = log(h); + auto log_e = log(e); + return best_fit_slope(log_h, log_e); +}; /** * @brief a verification problem for the J2 material model, taken from example 2 of @@ -340,6 +375,92 @@ TEST(J2, FrameIndifference) ASSERT_LT(norm(error), 1e-13*norm(internal_state.Fpinv)); } +TEST(J2, UniaxialRateDependentMaterial) +{ + /* + Uniaxial stress with a rate-dependent material (J2 finite deformation) + Linear hardening, linear rate dependence. + There exists an exact solution. This test verifies that the numerically + integrated plasticity solution converges at first order with timestep + refinement. + */ + + using Hardening = solid_mechanics::LinearHardening; + using Material = solid_mechanics::J2; + + double E = 1.0; + double sigma_y = 0.1; + double Hi = E/20.0; + double eta = 0.5; + double tau = eta/(Hi + E); + constexpr double strain_rate = 1.0; + + Hardening hardening{.sigma_y = sigma_y, .Hi = Hi, .eta = eta}; + Material material{.E = E, .nu = 0.25, .hardening = hardening, .density = 1.0}; + + auto internal_state = Material::State{}; + + auto strain = [=](double t) { + double true_strain = strain_rate*t; + // convert true strain to nominal strain (that is, to du_dX) + return std::expm1(true_strain); + }; + + auto plastic_strain_exact = [=](double t) { + double epsilon = strain_rate*t; // true strain + if (epsilon < (sigma_y / E)) { + return 0.0; + } else { + double rate_independent_part = (E * epsilon - sigma_y) / (E + Hi); + double t_y = sigma_y/(E * strain_rate); + double rate_dependent_part = E*eta*strain_rate/std::pow(E + Hi, 2)*std::expm1((t_y - t)/tau); + return rate_independent_part + rate_dependent_part; + } + }; + + // exact solution for Mandel/Kirchhoff stresses (they coincide in this problem) + auto stress_exact = [=](double t) { + double epsilon_p = plastic_strain_exact(t); + double epsilon = strain_rate*t; + return E*(epsilon - epsilon_p); + }; + + auto compute_solution_errors = [&](int timesteps) { + const double time_span = 5.0; + double dt = time_span/(timesteps - 1); + auto response_history = uniaxial_stress_test_rate_dependent(time_span, timesteps, material, internal_state, strain); + double max_plastic_strain_error{}, max_stress_error{}; + for (auto r : response_history) { + auto [t, dudx, P, state] = r; + + auto TK = dot(P, transpose(dudx + Identity<3>())); // Kirchhoff stress + double stress_error = std::abs(TK[0][0] - stress_exact(t)); + max_stress_error = std::max(max_stress_error, stress_error); + + auto ep = state.accumulated_plastic_strain; + double plastic_strain_error = std::abs(ep - plastic_strain_exact(t)); + max_plastic_strain_error = std::max(max_plastic_strain_error, plastic_strain_error); + } + return tuple{max_stress_error, max_plastic_strain_error, dt}; + }; + + constexpr int max_refinements = 4; + tensor stress_errors; + tensor dt_sizes; + int timesteps = 21; + for (int i = 0; i < max_refinements; i++) { + auto [stress_error, plastic_strain_error, dt] = compute_solution_errors(timesteps); + // std::cout << "dt = " << dt << " stress error = " << stress_error << std::endl; + stress_errors[i] = stress_error; + dt_sizes[i] = dt; + timesteps *= 2; + } + + double convergence_rate = compute_convergence_rate(dt_sizes, stress_errors); + std::cout << "convergence rate is " << convergence_rate << std::endl; + EXPECT_NEAR(convergence_rate, 1.0, 0.05); +}; + } // namespace serac int main(int argc, char* argv[]) From 5878f1cdb49fc4f86c8d8fd6bad306990660662e Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 8 Feb 2025 11:32:57 -0800 Subject: [PATCH 27/32] Doxygen comments --- .../physics/materials/tests/J2_material.cpp | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index 5f4045fd4e..d548a03440 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -21,6 +21,9 @@ namespace serac { +/** + * Element-wise logarithm + */ template tensor log(tensor v) { @@ -28,6 +31,9 @@ tensor log(tensor v) return make_tensor([&v](int i) { return log(v[i]); }); } +/** + * Arithmetic mean + */ template double mean(tensor v) { @@ -36,6 +42,9 @@ double mean(tensor v) return sum/n; } +/** + * Find slope for least-square fit of a linear function to data + */ template double best_fit_slope(tensor x, tensor y) { @@ -50,14 +59,6 @@ double best_fit_slope(tensor x, tensor y) return numer/denom; } -template -double compute_convergence_rate(tensor h, tensor e) -{ - auto log_h = log(h); - auto log_e = log(e); - return best_fit_slope(log_h, log_e); -}; - /** * @brief a verification problem for the J2 material model, taken from example 2 of * R. M. Brannon ยท S. Leelavanichkul (2009) @@ -378,11 +379,12 @@ TEST(J2, FrameIndifference) TEST(J2, UniaxialRateDependentMaterial) { /* - Uniaxial stress with a rate-dependent material (J2 finite deformation) - Linear hardening, linear rate dependence. - There exists an exact solution. This test verifies that the numerically - integrated plasticity solution converges at first order with timestep - refinement. + Verify that the constitutive model converges to exact solution at first order + with timestep refinement. + + Uses the trick that an exact solution for small strain plasticity can be used + for the log strain finite deformation model if the strains are interpreted as + log strain and the stress is interpreted as the Mandel stress. */ using Hardening = solid_mechanics::LinearHardening; @@ -402,7 +404,7 @@ TEST(J2, UniaxialRateDependentMaterial) auto strain = [=](double t) { double true_strain = strain_rate*t; - // convert true strain to nominal strain (that is, to du_dX) + // convert true strain to nominal strain (that is, to du_dX[0][0]) return std::expm1(true_strain); }; @@ -412,7 +414,7 @@ TEST(J2, UniaxialRateDependentMaterial) return 0.0; } else { double rate_independent_part = (E * epsilon - sigma_y) / (E + Hi); - double t_y = sigma_y/(E * strain_rate); + double t_y = sigma_y/(E * strain_rate); // time at first yield double rate_dependent_part = E*eta*strain_rate/std::pow(E + Hi, 2)*std::expm1((t_y - t)/tau); return rate_independent_part + rate_dependent_part; } @@ -444,20 +446,18 @@ TEST(J2, UniaxialRateDependentMaterial) return tuple{max_stress_error, max_plastic_strain_error, dt}; }; - constexpr int max_refinements = 4; - tensor stress_errors; - tensor dt_sizes; + constexpr int refinements = 4; + tensor stress_errors; + tensor dt_sizes; int timesteps = 21; - for (int i = 0; i < max_refinements; i++) { + for (int i = 0; i < refinements; i++) { auto [stress_error, plastic_strain_error, dt] = compute_solution_errors(timesteps); - // std::cout << "dt = " << dt << " stress error = " << stress_error << std::endl; stress_errors[i] = stress_error; dt_sizes[i] = dt; timesteps *= 2; } - double convergence_rate = compute_convergence_rate(dt_sizes, stress_errors); - std::cout << "convergence rate is " << convergence_rate << std::endl; + double convergence_rate = best_fit_slope(log(dt_sizes), log(stress_errors)); EXPECT_NEAR(convergence_rate, 1.0, 0.05); }; From 1391e48c560ab7eb704cc527a6094a265c7d5a6b Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 8 Feb 2025 11:34:49 -0800 Subject: [PATCH 28/32] Make style --- src/serac/physics/materials/solid_material.hpp | 8 ++++---- src/serac/physics/materials/tests/J2_material.cpp | 10 +++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/serac/physics/materials/solid_material.hpp b/src/serac/physics/materials/solid_material.hpp index b83cfe0816..4abe2f7032 100644 --- a/src/serac/physics/materials/solid_material.hpp +++ b/src/serac/physics/materials/solid_material.hpp @@ -200,10 +200,10 @@ struct LinearHardening { * @brief Power-law isotropic hardening law */ struct PowerLawHardening { - double sigma_y; ///< yield strength - double n; ///< hardening index in reciprocal form - double eps0; ///< reference value of accumulated plastic strain - double eta; ///< viscosity for linear rate sensitivity + double sigma_y; ///< yield strength + double n; ///< hardening index in reciprocal form + double eps0; ///< reference value of accumulated plastic strain + double eta; ///< viscosity for linear rate sensitivity /** * @brief Computes the flow stress diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index d548a03440..b6ed13c15f 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -27,7 +27,7 @@ namespace serac { template tensor log(tensor v) { - using std::log; + using std::log; return make_tensor([&v](int i) { return log(v[i]); }); } @@ -39,7 +39,7 @@ double mean(tensor v) { double sum = 0; for (int i = 0; i < n; i++) sum += v[i]; - return sum/n; + return sum / n; } /** @@ -53,10 +53,10 @@ double best_fit_slope(tensor x, tensor y) double numer = 0; double denom = 0; for (int i = 0; i < n; i++) { - numer += (x[i] - x_bar)*(y[i] - y_bar); - denom += (x[i] - x_bar)*(x[i] - x_bar); + numer += (x[i] - x_bar) * (y[i] - y_bar); + denom += (x[i] - x_bar) * (x[i] - x_bar); } - return numer/denom; + return numer / denom; } /** From 634a35f8af829177b3f56d5b19e9de45c690f3f9 Mon Sep 17 00:00:00 2001 From: Brandon Talamini <30813018+btalamini@users.noreply.github.com> Date: Sat, 8 Feb 2025 11:38:08 -0800 Subject: [PATCH 29/32] Make type of initializer explicit Co-authored-by: Eric B. Chin --- src/serac/physics/tests/solid_multi_material.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serac/physics/tests/solid_multi_material.cpp b/src/serac/physics/tests/solid_multi_material.cpp index 9aa914c805..471c402b13 100644 --- a/src/serac/physics/tests/solid_multi_material.cpp +++ b/src/serac/physics/tests/solid_multi_material.cpp @@ -198,7 +198,7 @@ TEST(Solid, MultiMaterialWithState) constexpr double Hi = E_right / 3.6; constexpr double sigma_y = 0.75 * applied_stress; - Hardening hardening{.sigma_y = sigma_y, .Hi = Hi, .eta = 0}; + Hardening hardening{.sigma_y = sigma_y, .Hi = Hi, .eta = 0.0}; MaterialRight mat_right{ .E = E_right, // Young's modulus .nu = nu_right, // Poisson's ratio From e930fa76b42a6e78a9935d6f1ae816e677b9214c Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 8 Feb 2025 11:59:06 -0800 Subject: [PATCH 30/32] Add rate sensitivity parameter to input parsing --- src/serac/physics/materials/hardening_input.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/serac/physics/materials/hardening_input.cpp b/src/serac/physics/materials/hardening_input.cpp index cc7a72c581..188550add3 100644 --- a/src/serac/physics/materials/hardening_input.cpp +++ b/src/serac/physics/materials/hardening_input.cpp @@ -13,6 +13,7 @@ void HardeningInputOptions::defineInputFileSchema(axom::inlet::Container& contai // Shared between both hardening laws container.addString("law", "Name of the hardening law (e.g. PowerLawHardening)").required(true); container.addDouble("sigma_y", "Yield strength"); + container.addDouble("eta", "Plastic viscosity"); // PowerLawHardening container.addDouble("n", "Hardening index in reciprocal form"); @@ -30,12 +31,13 @@ void HardeningInputOptions::defineInputFileSchema(axom::inlet::Container& contai bool eps0_present = c.contains("eps0") && (c["eps0"].type() == double_type); bool sigma_sat_present = c.contains("sigma_sat") && (c["sigma_sat"].type() == double_type); bool strain_constant_present = c.contains("strain_constant") && (c["strain_constant"].type() == double_type); + bool eta_present = c.contains("eta") && (c["eta"].type() == double_type); std::string law = c["law"]; if (law == "PowerLawHardening") { - return sigma_y_present && n_present && eps0_present && !sigma_sat_present && !sigma_sat_present; + return sigma_y_present && n_present && eps0_present && eta_present && !sigma_sat_present && !sigma_sat_present; } else if (law == "VoceHardening") { - return sigma_y_present && !n_present && !eps0_present && sigma_sat_present && strain_constant_present; + return sigma_y_present && eta_present && !n_present && !eps0_present && sigma_sat_present && strain_constant_present; } return false; @@ -50,10 +52,10 @@ serac::var_hardening_t FromInlet::operator()(const axom: std::string law = base["law"]; if (law == "PowerLawHardening") { result = - serac::solid_mechanics::PowerLawHardening{.sigma_y = base["sigma_y"], .n = base["n"], .eps0 = base["eps0"]}; + serac::solid_mechanics::PowerLawHardening{.sigma_y = base["sigma_y"], .n = base["n"], .eps0 = base["eps0"], .eta = base["eta"]}; } else if (law == "VoceHardening") { result = serac::solid_mechanics::VoceHardening{ - .sigma_y = base["sigma_y"], .sigma_sat = base["sigma_sat"], .strain_constant = base["strain_constant"]}; + .sigma_y = base["sigma_y"], .sigma_sat = base["sigma_sat"], .strain_constant = base["strain_constant"], .eta = base["eta"]}; } return result; } From d35eb2c90cedc2665e875856b849d684555883fc Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 8 Feb 2025 12:06:25 -0800 Subject: [PATCH 31/32] Eliminate annoying warnings for signed/unsigned integer conversions --- src/serac/physics/materials/tests/J2_material.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serac/physics/materials/tests/J2_material.cpp b/src/serac/physics/materials/tests/J2_material.cpp index b6ed13c15f..404b90ab34 100644 --- a/src/serac/physics/materials/tests/J2_material.cpp +++ b/src/serac/physics/materials/tests/J2_material.cpp @@ -430,7 +430,7 @@ TEST(J2, UniaxialRateDependentMaterial) auto compute_solution_errors = [&](int timesteps) { const double time_span = 5.0; double dt = time_span/(timesteps - 1); - auto response_history = uniaxial_stress_test_rate_dependent(time_span, timesteps, material, internal_state, strain); + auto response_history = uniaxial_stress_test_rate_dependent(time_span, size_t(timesteps), material, internal_state, strain); double max_plastic_strain_error{}, max_stress_error{}; for (auto r : response_history) { auto [t, dudx, P, state] = r; From 7e4a6f52b928e987d8bb83ce1da7be9f85e2384c Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sat, 8 Feb 2025 12:11:44 -0800 Subject: [PATCH 32/32] Make style --- src/serac/physics/materials/hardening_input.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/serac/physics/materials/hardening_input.cpp b/src/serac/physics/materials/hardening_input.cpp index 188550add3..9702558be2 100644 --- a/src/serac/physics/materials/hardening_input.cpp +++ b/src/serac/physics/materials/hardening_input.cpp @@ -37,7 +37,8 @@ void HardeningInputOptions::defineInputFileSchema(axom::inlet::Container& contai if (law == "PowerLawHardening") { return sigma_y_present && n_present && eps0_present && eta_present && !sigma_sat_present && !sigma_sat_present; } else if (law == "VoceHardening") { - return sigma_y_present && eta_present && !n_present && !eps0_present && sigma_sat_present && strain_constant_present; + return sigma_y_present && eta_present && !n_present && !eps0_present && sigma_sat_present && + strain_constant_present; } return false; @@ -51,11 +52,13 @@ serac::var_hardening_t FromInlet::operator()(const axom: serac::var_hardening_t result; std::string law = base["law"]; if (law == "PowerLawHardening") { - result = - serac::solid_mechanics::PowerLawHardening{.sigma_y = base["sigma_y"], .n = base["n"], .eps0 = base["eps0"], .eta = base["eta"]}; + result = serac::solid_mechanics::PowerLawHardening{ + .sigma_y = base["sigma_y"], .n = base["n"], .eps0 = base["eps0"], .eta = base["eta"]}; } else if (law == "VoceHardening") { - result = serac::solid_mechanics::VoceHardening{ - .sigma_y = base["sigma_y"], .sigma_sat = base["sigma_sat"], .strain_constant = base["strain_constant"], .eta = base["eta"]}; + result = serac::solid_mechanics::VoceHardening{.sigma_y = base["sigma_y"], + .sigma_sat = base["sigma_sat"], + .strain_constant = base["strain_constant"], + .eta = base["eta"]}; } return result; }