diff --git a/include/numsim-materials/default_materials.h b/include/numsim-materials/default_materials.h index 1ba4672..d433970 100644 --- a/include/numsim-materials/default_materials.h +++ b/include/numsim-materials/default_materials.h @@ -7,6 +7,7 @@ #include "numsim-materials/solvers/vector_newton.h" #include "numsim-materials/materials/scalar_stepper.h" #include "numsim-materials/materials/constant_scalar.h" +#include "numsim-materials/materials/props_scalar.h" #include "numsim-materials/materials/isotropic_tangent.h" #include "numsim-materials/materials/linear_elasticity.h" #include "numsim-materials/materials/linear_stress.h" @@ -73,6 +74,7 @@ void register_default_materials() { factory.template register_type>("scalar_stepper"); factory.template register_type>("linear_elasticity"); factory.template register_type>("constant_scalar"); + factory.template register_type>("props_scalar"); factory.template register_type>("isotropic_tangent"); factory.template register_type>("linear_stress"); factory.template register_type>("autocatalytic_reaction"); diff --git a/include/numsim-materials/io/json_parameter_converter.h b/include/numsim-materials/io/json_parameter_converter.h index a6aff2f..f6d60f5 100644 --- a/include/numsim-materials/io/json_parameter_converter.h +++ b/include/numsim-materials/io/json_parameter_converter.h @@ -93,6 +93,7 @@ json_reader_registry make_default_json_registry() { reg.template add(); reg.template add(); reg.template add(); + reg.template add(); reg.template add(); reg.template add(); diff --git a/include/numsim-materials/materials/props_scalar.h b/include/numsim-materials/materials/props_scalar.h new file mode 100644 index 0000000..2dd68bf --- /dev/null +++ b/include/numsim-materials/materials/props_scalar.h @@ -0,0 +1,76 @@ +#ifndef NUMSIM_MATERIALS_PROPS_SCALAR_H +#define NUMSIM_MATERIALS_PROPS_SCALAR_H + +#include +#include +#include +#include + +#include "numsim-materials/core/material_base.h" + +namespace numsim::materials { + +/// A scalar taken from the host's material-constants array on every call. +/// +/// constant_scalar bakes its number in at construction, which is right when the +/// constants are fixed for a material name (Abaqus PROPS). This one records only +/// WHICH slot it owns, for hosts whose constants vary — CalculiX interpolates +/// them by temperature. +/// +/// Rebuilding the graph per call would also be correct, since nothing is +/// retained between calls, but it measures 22x an evaluation (6.7 us against +/// 302 ns) where reading a slot is free. +/// +/// bind() DEREFERENCES and keeps no pointer: a host array may be a per-call +/// temporary, and a kept pointer would read a dead stack slot next call — +/// usually returning the right number, because the slot is commonly reused. +/// +/// Plain property, no update callback, as constant_scalar: nothing reaches +/// statev_map, and the value is in place before ctx.update() so ordering +/// cannot matter. +/// +/// Parameters: +/// "name": material name +/// "index": which host constant this publishes, 0-based +template +class props_scalar final : public material_base, Traits> { +public: + using base = material_base, Traits>; + using value_type = typename base::value_type; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit props_scalar(Args&&... args) + : base(std::forward(args)...), + m_value(base::template add_output("value")), + m_index(base::template get_parameter("index")) { + // Until the first bind(), rather than whatever the storage held. + m_value = value_type{}; + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + para.template insert("index").template add(); + return para; + } + + /// Copy this material's constant out of the host's array. + void bind(std::span props) { + if (m_index >= props.size()) + throw std::out_of_range( + "props_scalar '" + base::name() + "': wants constant " + + std::to_string(m_index) + " but only " + std::to_string(props.size()) + + " were supplied"); + m_value = props[m_index]; + } + + [[nodiscard]] std::size_t index() const noexcept { return m_index; } + +private: + value_type& m_value; + const std::size_t m_index; +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_PROPS_SCALAR_H diff --git a/include/numsim-materials/umat/json_model.h b/include/numsim-materials/umat/json_model.h index e8be7d4..21118c2 100644 --- a/include/numsim-materials/umat/json_model.h +++ b/include/numsim-materials/umat/json_model.h @@ -43,6 +43,10 @@ /// /// Named "constants", not "props": a PROPERTY here is a graph node, and it is /// what the deck calls them (*USER MATERIAL, CONSTANTS=). +/// +/// An entry naming a props_scalar binds the SLOT instead: that material reads +/// the number from the host each call rather than having it baked in. Swap the +/// type and nothing else in the document changes. namespace numsim::materials::umat { /// Host-driven source materials. Kept out of register_default_materials() so @@ -69,14 +73,31 @@ void ensure_materials_registered() { }); } -/// The parameter a binding will actually write. -/// -/// Its own function because it is the single place that has to stay in step -/// with the substitution loop below — validating one name and writing another -/// is how a target ends up half-checked. -inline std::string bound_parameter(const nlohmann::json& /*material*/, - const connection_source& binding) { - return binding.property; +/// What a binding writes into the document. One function for both the +/// validation and the substitution below — checking one parameter name and +/// writing another is how a target ends up half-verified. +struct constant_binding_target { + std::string parameter; + /// props_scalar is told its SLOT and reads the number itself each call; + /// everything else has the number written in now. + bool writes_slot{false}; +}; + +inline constant_binding_target bound_parameter( + const nlohmann::json& material, const connection_source& binding) { + if (material.value("type", std::string{}) == "props_scalar") { + // The target names the property the constant arrives on — "value", as for + // constant_scalar — so swapping the type does not force "constants" to be + // rewritten. The parameter written ("index") is therefore not what the + // document says. + if (binding.property != "value") + throw fatal_error( + "json_model: a props_scalar target names the property it publishes, " + "which is \"value\" — got \"" + binding.property + + "\"; the slot comes from the entry's position in \"constants\""); + return {"index", true}; + } + return {binding.property, false}; } /// Reject a target naming a parameter the material does not declare. @@ -102,7 +123,7 @@ void require_declared_parameter(const nlohmann::json& material, // Declared AND numeric. Every material declares "name", and most declare // *_source strings, so checking mere existence accepts targets that can only // fail later — with a JSON type error rather than anything about decks. - const auto wanted = bound_parameter(material, binding); + const auto wanted = bound_parameter(material, binding).parameter; const auto schema = factory.schema(type); std::vector numeric; bool declared = false, wanted_is_numeric = false; @@ -231,10 +252,16 @@ typename umat_registry::builder make_json_builder( // Into a copy, so the registered document stays a template. nlohmann::json doc = parsed; for (std::size_t i = 0; i < bindings.size(); ++i) - for (auto& material : doc["materials"]) - if (material.contains("name") && - material["name"].get() == bindings[i].material) - material[bindings[i].property] = props[i]; + for (auto& material : doc["materials"]) { + if (!material.contains("name") || + material["name"].get() != bindings[i].material) + continue; + // Resolved by the function that validated the target, so the two + // cannot drift. + const auto target = bound_parameter(material, bindings[i]); + material[target.parameter] = + target.writes_slot ? nlohmann::json(i) : nlohmann::json(props[i]); + } for (const auto& material : doc["materials"]) create_from_json(ctx, material); diff --git a/include/numsim-materials/umat/material_point_evaluator.h b/include/numsim-materials/umat/material_point_evaluator.h index 7f2bb4c..04ce250 100644 --- a/include/numsim-materials/umat/material_point_evaluator.h +++ b/include/numsim-materials/umat/material_point_evaluator.h @@ -1,6 +1,7 @@ #ifndef NUMSIM_MATERIALS_UMAT_MATERIAL_POINT_EVALUATOR_H #define NUMSIM_MATERIALS_UMAT_MATERIAL_POINT_EVALUATOR_H +#include #include #include #include @@ -11,6 +12,7 @@ #include #include "numsim-materials/core/material_context.h" +#include "numsim-materials/materials/props_scalar.h" #include "numsim-materials/umat/errors.h" #include "numsim-materials/umat/external_state_source.h" #include "numsim-materials/umat/statev_map.h" @@ -127,6 +129,16 @@ class material_point_evaluator { } m_statev = std::make_unique>(m_ctx, exclusions); + + // Collected rather than configured: a hand-listed set is a second thing to + // keep in step with the graph, and a missed entry is a stale modulus. + for (auto* material : m_ctx.materials()) + if (auto* reader = dynamic_cast*>(material)) { + m_props_readers.push_back(reader); + m_props_needed = std::max(m_props_needed, reader->index() + 1); + } + m_live_slots.assign(m_props_needed, false); + for (const auto* reader : m_props_readers) m_live_slots[reader->index()] = true; } /// Doubles this material needs in STATEV — what *DEPVAR must be at least. @@ -187,6 +199,16 @@ class material_point_evaluator { value_type dtime, value_type* stress6, value_type* tangent36, std::span drot = {}) { + // Every reader publishes zero until bind_props() runs: moduli of zero, an + // all-zero DDSDDE, and a host that fails to converge with nothing naming + // the cause. Catches never-bound, not stale — the plane-stress solve + // re-runs this per iterate without re-binding, so the two are + // indistinguishable from here. The registry binds on every call. + if (!m_props_readers.empty() && !m_props_bound) + throw fatal_error( + "material_point_evaluator: this model reads its material constants " + "per call — call bind_props() before evaluating"); + // STATEV is the only state store: reload it every call, so a repeated call // on an unconverged iterate starts from t_n, not from the previous trial. m_statev->unpack(statev); @@ -207,6 +229,37 @@ class material_point_evaluator { tangent_to_buffer(*m_tangent, tangent36); } + /// Copy the host's material constants into the graph, once per host call and + /// before evaluating. Separate from `call` because the plane-stress solve + /// runs the graph repeatedly for one call and the constants do not change + /// between iterates. A no-op without props_scalar materials. + void bind_props(std::span props) { + if (m_props_readers.empty()) return; + if (props.size() < m_props_needed) + throw fatal_error( + "material_point_evaluator: the model reads " + + std::to_string(m_props_needed) + + " host constants but this call supplied " + + std::to_string(props.size()) + " — check *USER MATERIAL, CONSTANTS="); + for (auto* reader : m_props_readers) reader->bind(props); + m_props_bound = true; + } + + /// True when any constant is read per call rather than baked into the graph. + [[nodiscard]] bool has_live_props() const noexcept { + return !m_props_readers.empty(); + } + + /// True when host constant @p slot is read per call, so a change there is + /// intended rather than a contradiction. + /// + /// Per slot, not per model: a document may mix the two binding times, and + /// answering for the whole model would wave the BAKED constants through as + /// well, leaving them silently at their first value. + [[nodiscard]] bool is_live_prop(std::size_t slot) const noexcept { + return slot < m_live_slots.size() && m_live_slots[slot]; + } + /// Write the updated history back. No commit(): the host owns the timestep. void store_statev(value_type* statev) const { m_statev->pack(statev); } @@ -344,6 +397,10 @@ class material_point_evaluator { const numsim_core::history_property* m_plastic_strain{nullptr}; std::unique_ptr> m_statev; + std::vector*> m_props_readers; + std::vector m_live_slots; + std::size_t m_props_needed{0}; + bool m_props_bound{false}; }; } // namespace numsim::materials::umat diff --git a/include/numsim-materials/umat/plane_stress_evaluator.h b/include/numsim-materials/umat/plane_stress_evaluator.h index a39b665..26b9465 100644 --- a/include/numsim-materials/umat/plane_stress_evaluator.h +++ b/include/numsim-materials/umat/plane_stress_evaluator.h @@ -66,6 +66,16 @@ class plane_stress_evaluator { /// Iterations the last evaluate() needed, for diagnostics. [[nodiscard]] int last_iterations() const noexcept { return m_last_iters; } + /// Once for the whole out-of-plane solve: the constants do not depend on the + /// iterate. + void bind_props(std::span props) { + m_inner->bind_props(props); + } + + [[nodiscard]] bool has_live_props() const noexcept { + return m_inner->has_live_props(); + } + void evaluate(const call& c) { if (c.statev.size() < nstatv()) throw fatal_error( diff --git a/include/numsim-materials/umat/umat_interface.h b/include/numsim-materials/umat/umat_interface.h index b10ce03..449113c 100644 --- a/include/numsim-materials/umat/umat_interface.h +++ b/include/numsim-materials/umat/umat_interface.h @@ -179,10 +179,13 @@ class umat_registry { void evaluate(std::string_view cmname, std::span props, const call& c) { auto& ts = thread_state_for(cmname, props); - if (c.ec == element_case::plane_stress) + if (c.ec == element_case::plane_stress) { + ts.ps->bind_props(props); ts.ps->evaluate(c); - else + } else { + ts.solid->bind_props(props); ts.solid->evaluate(c); + } } /// Drop this thread's cached contexts. Only needed if models are @@ -207,10 +210,9 @@ class umat_registry { std::unique_ptr ctx; std::unique_ptr solid; std::unique_ptr ps; - /// How many constants the context was built from. The graph is built once - /// and reused, so a later call arriving with a different count would mean - /// the cached parameters no longer describe this material. - std::size_t nprops{0}; + /// What the context was built from; the graph is reused, so different + /// constants on a later call would no longer describe this material. + std::vector props; }; static std::unordered_mapsecond.nprops != props.size()) + // PROPS cannot vary for one material name, so anything different + // contradicts the graph the constants were baked into. Values, not just + // the count: same length with different numbers is the case that reaches + // a material. Two faults, two messages — the check only ever explains. + if (it->second.props.size() != props.size()) throw fatal_error( - "numsim UMAT: material '" + std::string(key) + - "' was built from " + std::to_string(it->second.nprops) + - " constants but this call supplies " + - std::to_string(props.size()) + - " — PROPS must be constant for a given material name"); + "numsim UMAT: material '" + std::string(key) + "' was built from " + + std::to_string(it->second.props.size()) + + " constants but this call supplies " + std::to_string(props.size()) + + " — NPROPS cannot vary for a given material name"); + + // Then per SLOT, skipping the ones the model reads live. Asking + // has_live_props() for the whole model instead would wave a mixed + // document's BAKED constants through as well: they would keep the first + // call's values while the live ones tracked the deck, silently, which is + // the defect this check exists to catch. + for (std::size_t i = 0; i < props.size(); ++i) + if (!it->second.solid->is_live_prop(i) && + it->second.props[i] != props[i]) + throw fatal_error( + "numsim UMAT: material '" + std::string(key) + "' constant " + + std::to_string(i + 1) + " was baked into the graph as " + + std::to_string(it->second.props[i]) + " but this call supplies " + + std::to_string(props[i]) + + " — PROPS must be constant for a given material name; use " + "distinct *MATERIAL names for distinct constants, or a " + "props_scalar for a constant that genuinely varies per call"); return it->second; } @@ -269,7 +287,7 @@ class umat_registry { thread_state ts; ts.ctx = std::make_unique(); m.build(*ts.ctx, props); - ts.nprops = props.size(); + ts.props.assign(props.begin(), props.end()); if (!ts.ctx->is_finalized()) throw fatal_error( "the builder returned without calling finalize() on the context"); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f14a377..4425a79 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -25,6 +25,7 @@ add_numsim_test(test_plane_stress_evaluator test_plane_stress_evaluator.cpp) add_numsim_test(test_umat_interface test_umat_interface.cpp) add_numsim_test(test_tangent_generator test_tangent_generator.cpp) add_numsim_test(test_json_model test_json_model.cpp) +add_numsim_test(test_props_scalar test_props_scalar.cpp) target_link_libraries(test_umat_interface PRIVATE Threads::Threads) # Data dumper for plotting (not a test — standalone executable) diff --git a/tests/test_props_scalar.cpp b/tests/test_props_scalar.cpp new file mode 100644 index 0000000..e8026c7 --- /dev/null +++ b/tests/test_props_scalar.cpp @@ -0,0 +1,508 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/materials/constant_scalar.h" +#include "numsim-materials/materials/isotropic_tangent.h" +#include "numsim-materials/materials/linear_stress.h" +#include "numsim-materials/materials/props_scalar.h" +#include "numsim-materials/umat/json_model.h" + +// The Fortran-callable symbol, so live constants are exercised through the real +// ABI and not only through the C++ evaluator. +NUMSIM_MATERIALS_DEFINE_UMAT(numsim::materials::material_policy_default) + +namespace { + +namespace nm = numsim::materials; +namespace u = numsim::materials::umat; + +using policy = nm::material_policy_default; +using T = policy::value_type; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; +using registry = u::umat_registry; +using tensor2 = tmech::tensor; +using tensor4 = tmech::tensor; + +/// C1111 for an isotropic tangent, the quantity every test here reads back. +constexpr T c1111(T K, T G) { return K + 4.0 * G / 3.0; } + +// --------------------------------------------------------------------------- +// The material on its own +// --------------------------------------------------------------------------- + +TEST(PropsScalar, PublishesTheConstantAtItsIndex) { + ctx_type ctx; + param_type p; + p.insert("name", "K"); + p.insert("index", 1); + auto& k = ctx.create>(p); + ctx.finalize(); + + const T props[3] = {11.0, 22.0, 33.0}; + k.bind(std::span(props, 3)); + EXPECT_DOUBLE_EQ(ctx.get("K", "value"), 22.0); +} + +/// A definite value before the first bind, not whatever the storage held. +TEST(PropsScalar, IsZeroBeforeTheFirstBind) { + ctx_type ctx; + param_type p; + p.insert("name", "K"); + p.insert("index", 0); + ctx.create>(p); + ctx.finalize(); + EXPECT_DOUBLE_EQ(ctx.get("K", "value"), 0.0); +} + +TEST(PropsScalar, RejectsAnIndexPastTheSuppliedConstants) { + ctx_type ctx; + param_type p; + p.insert("name", "K"); + p.insert("index", 5); + auto& k = ctx.create>(p); + ctx.finalize(); + + const T props[2] = {1.0, 2.0}; + EXPECT_THROW(k.bind(std::span(props, 2)), std::out_of_range); +} + +/// Plain, not history: no STATEV slot. external_scalar_source would cost one +/// per constant. +TEST(PropsScalar, CostsNoStatevSlot) { + ctx_type ctx; + param_type p; + p.insert("name", "strain_in"); + ctx.create>(p); + p.clear(); + p.insert("name", "K"); + p.insert("index", 0); + ctx.create>(p); + p.clear(); + p.insert("name", "G"); + p.insert("index", 1); + ctx.create>(p); + p.clear(); + p.insert("name", "stiffness"); + p.insert("K_source", "K"); + p.insert("G_source", "G"); + ctx.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("tangent_source", "stiffness"); + p.insert("strain_source", "strain_in"); + ctx.create>(p); + ctx.finalize(); + + u::material_point_evaluator::config cfg; + cfg.strain_source = "strain_in"; + cfg.stress_source = "elastic"; + cfg.tangent_source = "stiffness"; + u::material_point_evaluator eval(ctx, cfg); + + EXPECT_EQ(eval.nstatv(), 0u); + EXPECT_TRUE(eval.has_live_props()); +} + +// --------------------------------------------------------------------------- +// Through the evaluator +// --------------------------------------------------------------------------- + +/// A model whose two moduli are live deck constants. +void build_live(ctx_type& ctx) { + param_type p; + p.insert("name", "strain_in"); + ctx.create>(p); + p.clear(); + p.insert("name", "K"); + p.insert("index", 0); + ctx.create>(p); + p.clear(); + p.insert("name", "G"); + p.insert("index", 1); + ctx.create>(p); + p.clear(); + p.insert("name", "stiffness"); + p.insert("K_source", "K"); + p.insert("G_source", "G"); + ctx.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("tangent_source", "stiffness"); + p.insert("strain_source", "strain_in"); + ctx.create>(p); + ctx.finalize(); +} + +u::material_point_evaluator::config live_config() { + u::material_point_evaluator::config cfg; + cfg.strain_source = "strain_in"; + cfg.stress_source = "elastic"; + cfg.tangent_source = "stiffness"; + return cfg; +} + +/// One uniaxial call, returning DDSDDE(1,1). +T call_once(u::material_point_evaluator& eval, + std::span props) { + const T stran[6] = {0, 0, 0, 0, 0, 0}; + const T dstran[6] = {0.001, 0, 0, 0, 0, 0}; + std::vector stress(6, 0.0), ddsdde(36, 0.0), statev; + u::material_point_evaluator::call c; + c.stran = stran; + c.dstran = dstran; + c.stress = stress; + c.ddsdde = ddsdde; + c.statev = statev; + eval.bind_props(props); + eval.evaluate(c); + return ddsdde[0]; +} + +/// The point: one graph, two constant sets, two stiffnesses, no rebuild. +TEST(PropsScalar, ChangedConstantsChangeTheTangentWithoutARebuild) { + ctx_type ctx; + build_live(ctx); + u::material_point_evaluator eval(ctx, live_config()); + + const T soft[2] = {100.0, 40.0}; + const T stiff[2] = {300.0, 140.0}; + + EXPECT_NEAR(call_once(eval, soft), c1111(100.0, 40.0), 1e-9); + EXPECT_NEAR(call_once(eval, stiff), c1111(300.0, 140.0), 1e-9); + // and back, so the second answer is not simply "the last one wins forever" + EXPECT_NEAR(call_once(eval, soft), c1111(100.0, 40.0), 1e-9); +} + +/// The stored-pointer trap. The buffer is overwritten AFTER bind_props and +/// before the graph runs: dereferencing at bind time makes that irrelevant, +/// keeping the pointer makes the tangent follow the clobbered values. +TEST(PropsScalar, ReadsTheConstantsAtBindTimeNotAtUpdateTime) { + ctx_type ctx; + build_live(ctx); + u::material_point_evaluator eval(ctx, live_config()); + + std::vector buffer{100.0, 40.0}; + eval.bind_props(buffer); + + // Everything the host promised is gone by the time the graph runs. + buffer[0] = -1.0e9; + buffer[1] = -1.0e9; + + const T stran[6] = {0, 0, 0, 0, 0, 0}; + const T dstran[6] = {0.001, 0, 0, 0, 0, 0}; + std::vector stress(6, 0.0), ddsdde(36, 0.0), statev; + u::material_point_evaluator::call c; + c.stran = stran; + c.dstran = dstran; + c.stress = stress; + c.ddsdde = ddsdde; + c.statev = statev; + eval.evaluate(c); + + EXPECT_NEAR(ddsdde[0], c1111(100.0, 40.0), 1e-9) + << "the constants must be copied at bind(), not read through a kept " + "pointer during update()"; +} + +/// Every reader publishes zero until bind_props runs, so an unbound model would +/// evaluate with moduli of zero and an all-zero DDSDDE. +TEST(PropsScalar, EvaluatingBeforeBindingIsFatal) { + ctx_type ctx; + build_live(ctx); + u::material_point_evaluator eval(ctx, live_config()); + + const T stran[6] = {0, 0, 0, 0, 0, 0}; + const T dstran[6] = {0.001, 0, 0, 0, 0, 0}; + std::vector stress(6, 0.0), ddsdde(36, 0.0), statev; + u::material_point_evaluator::call c; + c.stran = stran; + c.dstran = dstran; + c.stress = stress; + c.ddsdde = ddsdde; + c.statev = statev; + + EXPECT_THROW(eval.evaluate(c), u::fatal_error); + + // Bound, it evaluates normally — the guard must not be a one-way latch. + const T props[2] = {100.0, 40.0}; + eval.bind_props(props); + EXPECT_NO_THROW(eval.evaluate(c)); + EXPECT_NEAR(ddsdde[0], c1111(100.0, 40.0), 1e-9); +} + +/// The one place binding meets an ITERATIVE evaluator. If a change moved the +/// bind inside the loop, or dropped the forward, only this notices. +TEST(PropsScalar, PlaneStressUsesTheLiveConstants) { + ctx_type ctx; + build_live(ctx); + u::plane_stress_evaluator ps(ctx, live_config(), {}); + + constexpr T K = 100.0, G = 40.0; + const T props[2] = {K, G}; + const T stran[3] = {0, 0, 0}; + const T dstran[3] = {0.001, 0, 0}; + std::vector stress(3, 0.0), ddsdde(9, 0.0), statev(ps.nstatv(), 0.0); + u::material_point_evaluator::call c; + c.stran = stran; + c.dstran = dstran; + c.stress = stress; + c.ddsdde = ddsdde; + c.statev = statev; + c.ec = u::element_case::plane_stress; + + ps.bind_props(props); + ps.evaluate(c); + + // Condensed modulus, derived from the constants the host supplied. + const T E = 9 * K * G / (3 * K + G); + const T nu = (3 * K - 2 * G) / (2 * (3 * K + G)); + EXPECT_NEAR(ddsdde[0], E / (1 - nu * nu), 1e-9); + EXPECT_NEAR(stress[2], 0.0, 1e-10) << "sigma_33 must be driven to zero"; +} + +/// Too few constants is a setup error: fatal, not a cutback. +TEST(PropsScalar, TooFewConstantsIsFatalAtTheEvaluator) { + ctx_type ctx; + build_live(ctx); + u::material_point_evaluator eval(ctx, live_config()); + + const T only_one[1] = {100.0}; + EXPECT_THROW(eval.bind_props(std::span(only_one, 1)), u::fatal_error); +} + +/// A baked-constant model reports itself as such and pays for nothing. +TEST(PropsScalar, BakedConstantModelsHaveNoLiveProps) { + ctx_type ctx; + param_type p; + p.insert("name", "strain_in"); + ctx.create>(p); + p.clear(); + p.insert("name", "K"); + p.insert("value", 100.0); + ctx.create>(p); + p.clear(); + p.insert("name", "G"); + p.insert("value", 40.0); + ctx.create>(p); + p.clear(); + p.insert("name", "stiffness"); + p.insert("K_source", "K"); + p.insert("G_source", "G"); + ctx.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("tangent_source", "stiffness"); + p.insert("strain_source", "strain_in"); + ctx.create>(p); + ctx.finalize(); + + u::material_point_evaluator eval(ctx, live_config()); + EXPECT_FALSE(eval.has_live_props()); + // bind_props is then a no-op, including for an empty array. + EXPECT_NO_THROW(eval.bind_props({})); +} + +// --------------------------------------------------------------------------- +// Through JSON and the real umat_ entry point +// --------------------------------------------------------------------------- + +/// Identical to the baked document but for the material type — the claim the +/// JSON layer makes, so it is asserted rather than described. +const char* kLive = R"({ + "materials": [ + {"type": "external_strain_source", "name": "strain_in"}, + {"type": "props_scalar", "name": "K"}, + {"type": "props_scalar", "name": "G"}, + {"type": "isotropic_tangent", "name": "stiffness", + "K_source": "K", "G_source": "G"}, + {"type": "linear_stress", "name": "elastic", + "tangent_source": "stiffness", "strain_source": "strain_in"} + ], + "constants": ["K::value", "G::value"] +})"; + +registry::config json_config() { + registry::config cfg; + cfg.strain_source = "strain_in"; + cfg.stress_source = "elastic"; + cfg.tangent_source = "stiffness"; + return cfg; +} + +/// K baked, G live. A document may mix the two binding times, so the +/// consistency check has to be answered per slot rather than per model. +const char* kMixed = R"({ + "materials": [ + {"type": "external_strain_source", "name": "strain_in"}, + {"type": "constant_scalar", "name": "K", "value": 0}, + {"type": "props_scalar", "name": "G"}, + {"type": "isotropic_tangent", "name": "stiffness", + "K_source": "K", "G_source": "G"}, + {"type": "linear_stress", "name": "elastic", + "tangent_source": "stiffness", "strain_source": "strain_in"} + ], + "constants": ["K::value", "G::value"] +})"; + +struct Registration { + Registration() { + u::register_json_model("LIVEELASTIC", kLive, json_config()); + // One name per test: each warms the cache itself, so sharing would let + // test order decide the outcome. + u::register_json_model("MIXEDBAKED", kMixed, json_config()); + u::register_json_model("MIXEDLIVE", kMixed, json_config()); + } +}; +const Registration registration_{}; + +struct fortran_name { + char buf[80]; + explicit fortran_name(const std::string& s) { + for (auto& c : buf) c = ' '; + for (std::size_t i = 0; i < s.size() && i < 80; ++i) buf[i] = s[i]; + } +}; + +T uniaxial_tangent(const std::string& name, const T* props, int nprops) { + const fortran_name cm(name); + T statev[1] = {0}; + const T stran[6] = {0, 0, 0, 0, 0, 0}; + const T dstran[6] = {0.001, 0, 0, 0, 0, 0}; + T stress[6] = {0}, ddsdde[36] = {0}, pnewdt = 1.0; + T sse = 0, spd = 0, scd = 0, rpl = 0, ddsddt[6] = {0}, drplde[6] = {0}; + T drpldt = 0; + const T time[2] = {0, 0}; + T dtime = 0.1; + const T temp = 0, dtemp = 0, predef = 0, dpred = 0, celent = 1; + const T coords[3] = {0}, drot[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + const T dfg[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + int noel = 1, npt = 1, layer = 1, kspt = 1, jstep = 1, kinc = 1; + int ndi = 3, nshr = 3, ntens = 6, nstatv = 0; + + umat_(stress, statev, ddsdde, &sse, &spd, &scd, &rpl, ddsddt, drplde, &drpldt, + stran, dstran, time, &dtime, &temp, &dtemp, &predef, &dpred, cm.buf, + &ndi, &nshr, &ntens, &nstatv, props, &nprops, coords, drot, &pnewdt, + &celent, dfg, dfg, &noel, &npt, &layer, &kspt, &jstep, &kinc, 80); + + EXPECT_DOUBLE_EQ(pnewdt, 1.0); + return ddsdde[0]; +} + +/// The registry's plane-stress dispatch, which binds through a different call +/// than the solid path. PropsScalar.PlaneStressUsesTheLiveConstants covers the +/// evaluator's forward; nothing covered the registry line that calls it, so +/// deleting it passed the whole suite. +TEST(PropsScalarJson, PlaneStressBindsThroughTheRegistry) { + constexpr T K = 100.0, G = 40.0; + const T props[2] = {K, G}; + const fortran_name cm("LIVEELASTIC"); + + T statev[2] = {0, 0}; + const T stran[3] = {0, 0, 0}; + const T dstran[3] = {0.001, 0, 0}; + T stress[3] = {0}, ddsdde[9] = {0}, pnewdt = 1.0; + T sse = 0, spd = 0, scd = 0, rpl = 0, ddsddt[3] = {0}, drplde[3] = {0}; + T drpldt = 0; + const T time[2] = {0, 0}; + T dtime = 0.1; + const T temp = 0, dtemp = 0, predef = 0, dpred = 0, celent = 1; + const T coords[3] = {0}, drot[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + const T dfg[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1}; + int noel = 1, npt = 1, layer = 1, kspt = 1, jstep = 1, kinc = 1; + int ndi = 2, nshr = 1, ntens = 3, nstatv = 2, nprops = 2; + + umat_(stress, statev, ddsdde, &sse, &spd, &scd, &rpl, ddsddt, drplde, &drpldt, + stran, dstran, time, &dtime, &temp, &dtemp, &predef, &dpred, cm.buf, + &ndi, &nshr, &ntens, &nstatv, props, &nprops, coords, drot, &pnewdt, + &celent, dfg, dfg, &noel, &npt, &layer, &kspt, &jstep, &kinc, 80); + + EXPECT_DOUBLE_EQ(pnewdt, 1.0); + const T E = 9 * K * G / (3 * K + G); + const T nu = (3 * K - 2 * G) / (2 * (3 * K + G)); + EXPECT_NEAR(ddsdde[0], E / (1 - nu * nu), 1e-9) + << "the constants must reach the material on the plane-stress path too"; +} + +/// Position is the slot: no "index" in the document, yet K takes PROPS[0]. +TEST(PropsScalarJson, ConstantsPositionBindsTheSlot) { + const T props[2] = {250.0, 90.0}; + EXPECT_NEAR(uniaxial_tangent("LIVEELASTIC", props, 2), c1111(250.0, 90.0), + 1e-9); +} + +/// The target names the property the constant arrives on — "value" — not the +/// "index" the builder writes. Identical spelling to constant_scalar is what +/// lets the type be swapped, so the wrong one is rejected. +TEST(PropsScalarJson, RejectsATargetNamingIndexRatherThanTheProperty) { + const char* names_index = R"({ + "materials": [{"type": "props_scalar", "name": "K"}], + "constants": ["K::index"] + })"; + EXPECT_THROW(u::make_json_builder(names_index), u::fatal_error); + + const char* names_property = R"({ + "materials": [{"type": "props_scalar", "name": "K"}], + "constants": ["K::value"] + })"; + EXPECT_NO_THROW(u::make_json_builder(names_property)); +} + +/// The registry's consistency check stands down for live constants: nothing +/// was baked in for a changed array to contradict. +TEST(PropsScalarJson, ChangedDeckConstantsAreHonouredNotRejected) { + const T first[2] = {100.0, 40.0}; + const T second[2] = {300.0, 140.0}; + + EXPECT_NEAR(uniaxial_tangent("LIVEELASTIC", first, 2), c1111(100.0, 40.0), + 1e-9); + EXPECT_NEAR(uniaxial_tangent("LIVEELASTIC", second, 2), c1111(300.0, 140.0), + 1e-9); +} + +// --------------------------------------------------------------------------- +// Mixed documents: baked and live constants side by side +// --------------------------------------------------------------------------- + +/// A changed BAKED constant is still fatal beside a live one. Answering +/// has_live_props() per model let this through: G tracked the deck while K +/// kept its first value, with no diagnostic. +TEST(PropsScalarJson, ChangingABakedConstantIsFatalEvenBesideALiveOne) { + const T first[2] = {100.0, 40.0}; + const T changed_both[2] = {300.0, 140.0}; // K baked, G live + + ASSERT_NEAR(uniaxial_tangent("MIXEDBAKED", first, 2), c1111(100.0, 40.0), + 1e-9); + + int fatal_count = 0; + static int* counter = &fatal_count; + u::set_fatal_handler([](const char*) { ++*counter; }); + const T got = uniaxial_tangent("MIXEDBAKED", changed_both, 2); + u::set_fatal_handler(nullptr); + + EXPECT_EQ(fatal_count, 1) + << "a changed baked constant must be reported, not served from the " + "cached graph because some OTHER constant is live"; + // The wrong-but-plausible answer this used to return. + EXPECT_FALSE(std::abs(got - c1111(100.0, 140.0)) < 1e-9) + << "K kept its first value while G followed the deck"; +} + +/// The other half, so the fix cannot pass by rejecting every mixed document. +TEST(PropsScalarJson, ChangingOnlyTheLiveConstantIsHonouredInAMixedDocument) { + const T first[2] = {100.0, 40.0}; + const T live_changed[2] = {100.0, 140.0}; // K unchanged, G changed + + EXPECT_NEAR(uniaxial_tangent("MIXEDLIVE", first, 2), c1111(100.0, 40.0), + 1e-9); + EXPECT_NEAR(uniaxial_tangent("MIXEDLIVE", live_changed, 2), + c1111(100.0, 140.0), 1e-9); +} + +} // namespace diff --git a/tests/test_umat_interface.cpp b/tests/test_umat_interface.cpp index 72dafc2..5f7078c 100644 --- a/tests/test_umat_interface.cpp +++ b/tests/test_umat_interface.cpp @@ -166,6 +166,9 @@ struct Registration { // against an already-built name the NPROPS-consistency check fires first // and require_props is never reached. registry::instance().register_model("COLDNAME", build_deck_elastic, de); + // One test only: it warms the cache itself, so a shared name would let + // test order decide the outcome. + registry::instance().register_model("VALUEPROBE", build_deck_elastic, de); // Deliberately lower-case, to prove the registry folds case on both sides. registry::config lc; @@ -765,7 +768,43 @@ TEST(UmatInterface, ChangingNpropsForTheSameNameIsFatal) { call_umat("STIFF", stress, statev.data(), ddsdde, stran, dstran, 0.0, 0.1, 3, 3, 6, 0, &pnewdt, nullptr, nullptr, nullptr, nullptr, three, 3); EXPECT_EQ(FatalProbe::count, 1); - EXPECT_NE(FatalProbe::last.find("constant"), std::string::npos) + // Both counts, so NPROPS=3 is not left to be matched against one number. + EXPECT_NE(FatalProbe::last.find("2 constants"), std::string::npos) + << FatalProbe::last; + EXPECT_NE(FatalProbe::last.find("supplies 3"), std::string::npos) + << FatalProbe::last; +} + +/// Same count, different numbers — the case that reaches a material. A count +/// check accepts it and serves the first call's stiffness forever: a converged +/// analysis with the wrong moduli. +TEST(UmatInterface, ChangingPropsValuesForTheSameNameIsFatal) { + FatalProbe probe; + std::vector statev(1, 0.0); + const T stran[6] = {0, 0, 0, 0, 0, 0}; + const T dstran[6] = {0.001, 0, 0, 0, 0, 0}; + T stress[6] = {0}, ddsdde[36] = {0}, pnewdt = 1.0; + const T soft[2] = {100.0, 40.0}; + const T stiff[2] = {300.0, 140.0}; + + // Builds the context, and shows which constants it was built from. + call_umat("VALUEPROBE", stress, statev.data(), ddsdde, stran, dstran, 0.0, 0.1, + 3, 3, 6, 0, &pnewdt, nullptr, nullptr, nullptr, nullptr, soft, 2); + ASSERT_EQ(FatalProbe::count, 0); + ASSERT_NEAR(ddsdde[0], 100.0 + 4.0 * 40.0 / 3.0, 1e-9); + + // Same name, same count, different values. + call_umat("VALUEPROBE", stress, statev.data(), ddsdde, stran, dstran, 0.0, 0.1, + 3, 3, 6, 0, &pnewdt, nullptr, nullptr, nullptr, nullptr, stiff, 2); + EXPECT_EQ(FatalProbe::count, 1) + << "a same-length PROPS with different numbers must be reported, not " + "silently served from the cached graph"; + // Must name which constant disagrees, and both values. + EXPECT_NE(FatalProbe::last.find("constant 1"), std::string::npos) + << FatalProbe::last; + EXPECT_NE(FatalProbe::last.find("100.0"), std::string::npos) + << FatalProbe::last; + EXPECT_NE(FatalProbe::last.find("300.0"), std::string::npos) << FatalProbe::last; }