diff --git a/README.md b/README.md index 4bd8b2d..23f2d55 100755 --- a/README.md +++ b/README.md @@ -52,6 +52,8 @@ material_context (facade) - `.reads()` — the parameters alone; `.reads_all_inputs()` — every input the material has at `finalize()` - `.stateful()` — keeps state of its own or acts outside its outputs: runs on every update - `.drives(local inputs...)` — reruns their source callbacks (a solver): runs after what those read +- `converges(residual, tolerance, max_passes)` — a fixed point over local edges: `update()` reruns the callbacks + between the local edge and `residual` until it converges, and throws `cutback_error` after `max_passes` - Debug checks (CMake option `NUMSIM_MATERIALS_DEBUG_CHECKS`, on except in release configurations) rerun the constants and throw when one changes (`ctx.verify_constants(bool)`), and throw when a callback reads an input it did not declare (`ctx.check_reads(bool)`); `ctx.is_constant("mat", "prop")` diff --git a/docs/materials/mean_field_damage/report.json b/docs/materials/mean_field_damage/report.json index 16eb2b0..4ab0816 100644 --- a/docs/materials/mean_field_damage/report.json +++ b/docs/materials/mean_field_damage/report.json @@ -1,7 +1,7 @@ { "title": "mean_field_strain and mean_field_composite_local_nonlinear with a damaging matrix", "include": ["../standard_small_strain.json"], - "notes": "Mori-Tanaka with local strains (incremental: the local strains advance with the concentrations of the current tangents, so after a damage-loading reversal they do not retrace and the composite keeps a residual stress at zero strain, visible against the intact Mori-Tanaka): mean_field_strain produces the phase strains from the composite increment and the concentration tensors of the algorithmic tangents (the matrix phase reads its tangent over a local edge, the previous pass); the damage chain and the stiff spheres are evaluated at those strains; mean_field_composite_local_nonlinear averages them. The evaluator repeats the pass until mt::residual is below tolerance. The published tangent sum c_r C_r : A_r is first order in the strain increment (exact as it vanishes, see test_mean_field_local); across the damage-surface kinks a central difference does not approximate it, so the tangent check is off here.", + "notes": "Mori-Tanaka with local strains (incremental: the local strains advance with the concentrations of the current tangents, so after a damage-loading reversal they do not retrace and the composite keeps a residual stress at zero strain, visible against the intact Mori-Tanaka): mean_field_strain produces the phase strains from the composite increment and the concentration tensors of the algorithmic tangents (the matrix phase reads its tangent over a local edge, the previous pass); the damage chain and the stiff spheres are evaluated at those strains; mean_field_composite_local_nonlinear averages them. update() repeats the pass until mt::residual is below the composite's tolerance (converges()). The published tangent sum c_r C_r : A_r is first order in the strain increment (exact as it vanishes, see test_mean_field_local); across the damage-surface kinks a central difference does not approximate it, so the tangent check is off here.", "parameters": {"Km": 3333.33, "Gm": 1153.85, "Ki": 166666.67, "Gi": 76923.08, "c": 0.3, "cm": 0.7, "kappa_0": 0.002, "beta": 30.0, "amplitude": 0.008, "tangent_every": 0}, "model": {"materials": [ {"type": "external_strain_source", "name": "strain"}, @@ -23,14 +23,14 @@ {"type": "exponential_damage_law", "name": "law", "yield_source": "yield", "kappa_0": "$kappa_0", "beta": "$beta"}, {"type": "isotropic_damage", "name": "damage", "elastic_source": "elastic", "damage_source": "law", "state_source": "state", "yield_source": "yield"}, {"type": "linear_stress", "name": "stiff", "tangent_source": "C_i", "strain_source": "mt_strain", "strain_property": "strain_spheres"}, - {"type": "mean_field_composite_local_nonlinear", "name": "mt", "phases": [ + {"type": "mean_field_composite_local_nonlinear", "name": "mt", "tolerance": 1e-10, "max_passes": 200, "phases": [ {"concentration": "mt_strain::concentration_matrix", "stress": "damage::damaged_stress", "tangent": "damage::damaged_tangent", "fraction": "matrix::volume_fraction", "used_tangent": "matrix::tangent"}, {"concentration": "mt_strain::concentration_spheres", "stress": "stiff::stress", "tangent": "C_i::tangent", "fraction": "spheres::volume_fraction"}]}, {"type": "phase", "name": "matrix_intact", "tangent_source": "C_m", "fraction_source": "cm_"}, {"type": "mean_field_composite_global", "name": "mt_intact", "phases": ["matrix_intact", "spheres"], "kernel": {"type": "mori_tanaka"}}, {"type": "linear_stress", "name": "intact", "tangent_source": "mt_intact", "strain_source": "strain"}]}, - "point": {"strain_source": "strain", "stress_source": "mt", "residual_source": "mt::residual", "residual_tolerance": 1e-10, "max_passes": 200}, + "point": {"strain_source": "strain", "stress_source": "mt"}, "record": {"d": "law::damage", "eps_m": "mt_strain::strain_matrix", "eps_i": "mt_strain::strain_spheres", "sig_m": "damage::damaged_stress", "sig_i": "stiff::stress", "res": "mt::residual", "sig_intact": "intact::stress"}, "checks": [ {"name": "local strains partition the composite strain", "cases": "*", "value": "(1 - c) * eps_m_11 + c * eps_i_11", "expected": "strain_11", "tolerance": 1e-10, "relative": true}, diff --git a/docs/materials/mean_field_reuss_plastic/report.json b/docs/materials/mean_field_reuss_plastic/report.json index 3bad506..abfc9d3 100644 --- a/docs/materials/mean_field_reuss_plastic/report.json +++ b/docs/materials/mean_field_reuss_plastic/report.json @@ -3,7 +3,7 @@ "include": [ "../standard_small_strain.json" ], - "notes": "Equal stress in series through the local pair: the concentrations A_r = C_r^-1 : C_R are formed with the algorithmic tangents of the previous pass, the matrix phase reads its tangent over a local edge, and the point repeats the pass until the residual is small. With a hardening phase the phase stresses agree to first order in the step (checked loosely here at 50 steps; see test_mean_field_damage for the refinement). The published tangent is first order in the step too, so the finite-difference check is off.", + "notes": "Equal stress in series through the local pair: the concentrations A_r = C_r^-1 : C_R are formed with the algorithmic tangents of the previous pass, the matrix phase reads its tangent over a local edge, and update() repeats the pass until the residual is below the composite's tolerance. With a hardening phase the phase stresses agree to first order in the step (checked loosely here at 50 steps; see test_mean_field_damage for the refinement). The published tangent is first order in the step too, so the finite-difference check is off.", "parameters": { "Km": 3333.33, "Gm": 1153.85, @@ -121,6 +121,8 @@ { "type": "mean_field_composite_local_nonlinear", "name": "mt", + "tolerance": 1e-10, + "max_passes": 200, "phases": [ { "concentration": "mt_strain::concentration_matrix", @@ -141,10 +143,7 @@ }, "point": { "strain_source": "strain", - "stress_source": "mt", - "residual_source": "mt::residual", - "residual_tolerance": 1e-10, - "max_passes": 200 + "stress_source": "mt" }, "record": { "eps_m": "mt_strain::strain_matrix", diff --git a/docs/materials/short_fibre_damage/report.json b/docs/materials/short_fibre_damage/report.json index 34b48cb..050dd25 100644 --- a/docs/materials/short_fibre_damage/report.json +++ b/docs/materials/short_fibre_damage/report.json @@ -1,7 +1,7 @@ { "title": "orientation_phases: damaging short fibres at their local strains", "include": ["../standard_small_strain.json"], - "notes": "Short fibres in an elastic matrix, orientation tensor diag(0.7, 0.25, 0.05), Lebedev degree 11: orientation_phases expands into 25 fibre phases (one per direction of the half rule), each with its own damage chain at its own local strain, Mori-Tanaka with an Eshelby spheroid per direction. sf localises with the current (damaged) tangents of the fibres (localisation nonlinear: the phase nodes read them over a local edge, the evaluator repeats the pass until sf::residual is below tolerance); sfl localises with the initial tangents (localisation linear: one pass); intact is short_fibre_composite without damage. Fibre 0 lies along x1, fibres 1 and 2 along x2 and x3. A softening fibre attracts strain under Mori-Tanaka, so the nonlinear localisation damages the fibres along the load earlier; under uniaxial strain it is also the softer composite (the linear localisation averages damaged stresses over strains formed with the intact concentrations, so no ordering holds in general). The published tangent sum c_r C_r : A_r holds the concentrations fixed, first order in the strain increment (see test_orientation_phases), so the central-difference tangent check is off.", + "notes": "Short fibres in an elastic matrix, orientation tensor diag(0.7, 0.25, 0.05), Lebedev degree 11: orientation_phases expands into 25 fibre phases (one per direction of the half rule), each with its own damage chain at its own local strain, Mori-Tanaka with an Eshelby spheroid per direction. sf localises with the current (damaged) tangents of the fibres (localisation nonlinear: the phase nodes read them over a local edge, update() repeats the pass until sf::residual is below the composite's tolerance); sfl localises with the initial tangents (localisation linear: one pass); intact is short_fibre_composite without damage. Fibre 0 lies along x1, fibres 1 and 2 along x2 and x3. A softening fibre attracts strain under Mori-Tanaka, so the nonlinear localisation damages the fibres along the load earlier; under uniaxial strain it is also the softer composite (the linear localisation averages damaged stresses over strains formed with the intact concentrations, so no ordering holds in general). The published tangent sum c_r C_r : A_r holds the concentrations fixed, first order in the strain increment (see test_orientation_phases), so the central-difference tangent check is off.", "parameters": {"Km": 3333.33, "Gm": 1153.85, "Kf": 166666.67, "Gf": 76923.08, "c": 0.2, "cm": 0.8, "aspect_ratio": 15.0, "kappa_0": 0.002, "beta": 30.0, "amplitude": 0.008, "tangent_every": 0}, "model": {"materials": [ {"type": "external_strain_source", "name": "strain"}, @@ -14,7 +14,7 @@ {"type": "phase", "name": "matrix", "tangent_source": "C_m", "fraction_source": "cm_"}, {"type": "linear_elasticity", "name": "elastic_m", "strain_producer_name": "sf_strain", "strain_property": "strain_matrix", "K": "$Km", "G": "$Gm"}, - {"type": "orientation_phases", "name": "sf", "localisation": "nonlinear", "strain_source": "strain", "matrix": "matrix", + {"type": "orientation_phases", "name": "sf", "localisation": "nonlinear", "tolerance": 1e-10, "max_passes": 200, "strain_source": "strain", "matrix": "matrix", "matrix_stress": "elastic_m::stress", "fraction_source": "c_", "orientation": [0.7, 0.25, 0.05, 0.0, 0.0, 0.0], "precision": 11, "kernel": {"type": "mori_tanaka", "eshelby": {"type": "eshelby_spheroid", "aspect_ratio": "$aspect_ratio", "nu_source": "moduli", "nu_property": "poisson_ratio"}}, "phase": {"tangent_source": "fibre_{q}", "tangent_property": "damaged_tangent", "tangent_edge": "local", "initial_tangent_source": "C_f"}, @@ -44,7 +44,7 @@ {"type": "short_fibre_composite", "name": "intact_tangent", "matrix": "matrix", "fibre": "fibre", "orientation": [0.7, 0.25, 0.05, 0.0, 0.0, 0.0], "precision": 11}, {"type": "linear_stress", "name": "intact", "tangent_source": "intact_tangent", "strain_source": "strain"} ]}, - "point": {"strain_source": "strain", "stress_source": "sf", "residual_source": "sf::residual", "residual_tolerance": 1e-10, "max_passes": 200}, + "point": {"strain_source": "strain", "stress_source": "sf"}, "record": {"d_along": "fl_0::damage", "d_x2": "fl_1::damage", "d_x3": "fl_2::damage", "dl_along": "lfl_0::damage", "res": "sf::residual", "sig_lin": "sfl::stress", "sig_intact": "intact::stress"}, "checks": [ diff --git a/docs/mean_field.md b/docs/mean_field.md index 9102971..d2640c2 100644 --- a/docs/mean_field.md +++ b/docs/mean_field.md @@ -241,15 +241,23 @@ Under Mori–Tanaka a softening fibre then attracts strain, so it damages more t ## 5. Who runs the fixed point -Local edges make every pass acyclic; the fixed point is reached by repeating `update()`: - -- `material_point_evaluator`: `residual_property` (`"material::property"`), `residual_tolerance`, `max_passes` - — the graph is updated again per host call until the residual is at most the tolerance; the UMAT and the - drivers go through it. In a report document: `"point": {..., "residual": "mt::residual", "tolerance": 1e-10, - "max_passes": 200}`. -- A solver's Newton loop re-evaluates the material each iteration; the material's own loop above still runs - so that the returned stress is consistent for the given increment. -- In a test: `while (ctx.get("mt", "residual") > tol) ctx.update();`. +Local edges make every pass acyclic. The engine runs the fixed point inside `update()`: + +- `mean_field_composite_local_nonlinear` declares its `residual` with `converges()`. Its parameters are + `tolerance` (default 1e-10) and `max_passes` (default 200). +- The engine reruns the iteration group until the residual is at most the tolerance. The group is the + callbacks from the phase that reads its tangent over the local edge to the composite. +- The callbacks before the group (the strain source, the constant tangents) and after it (recorders, plots) + run once per update. +- After `max_passes` passes, `update()` throws `cutback_error`, and the host cuts the step back. +- Every host converges without further setup: the UMAT, the drivers, the report and a test calling + `ctx.update()`. A solver's Newton loop re-evaluates the material each iteration, so the returned stress is + consistent for the given increment. +- Deprecated: the material point's `residual_source` (`residual_property` of `material_point_evaluator`, + `"point": {"residual_source": ...}` in a report) with `residual_tolerance` and `max_passes` repeats the whole + update. With a converging composite it finds the residual converged after the first update. In a report, set + `"tolerance"` and `"max_passes"` on the composite instead (on the `orientation_phases` entry for a + generated one). It is a Picard iteration on the tangents — linear convergence, fine for damage with moderate softening; near snap-back it can stall, which the residual then shows. The self-consistent kernel needs none of this: its own diff --git a/docs/uvwmaterial-mapping.md b/docs/uvwmaterial-mapping.md index f18a7ba..2a88954 100644 --- a/docs/uvwmaterial-mapping.md +++ b/docs/uvwmaterial-mapping.md @@ -35,7 +35,7 @@ is a `*_source` parameter. `tools/convert_uvwmaterial_card.py` converts the old | `eigenstrain_curing_material_base`, `degree_of_cure_dependent_material` | `thermal_eigenstrain` with the cure as source; `inclusion.fraction_source`; `weighted_sum` | | #137 | | `small_strain_rule_of_mixture_history` (Voigt) | `weighted_sum` with `constant_scalar` fractions: every phase reads the composite strain, the sum is Σ cᵢσᵢ, Σ cᵢCᵢ (`tests/test_voigt_composite.cpp`) | `terms`, weights | #117; not from JSON yet (no reader for `terms`) | | `small_strain_rule_of_mixture_history` (Reuss) | — | | #117, after the vector Newton lands | -| `small_strain_mean_field_composite_history`, `_plastic`, `_damage` | `mean_field_strain` + the phase materials + `mean_field_composite_local_linear` / `_nonlinear` | the strain node produces `strain_` (history) and the concentrations; the materials consume them (`strain_source`, `strain_property`); the stress node averages what they produce; a nonlinear phase reads its tangent over a local edge (`tangent_edge`) and the point closes the fixed point with `residual` | #116 | +| `small_strain_mean_field_composite_history`, `_plastic`, `_damage` | `mean_field_strain` + the phase materials + `mean_field_composite_local_linear` / `_nonlinear` | the strain node produces `strain_` (history) and the concentrations; the materials consume them (`strain_source`, `strain_property`); the stress node averages what they produce; a nonlinear phase reads its tangent over a local edge (`tangent_edge`) and the `_nonlinear` composite closes the fixed point in `update()` (`converges()` on its `residual`) | #116 | | `small_strain_short_fibre_history_composite` | — | | #118 | | `*_nonlocal*`, `element_deletion` | — | | #121 | | `finite_strain_mean_field_composite*` | — | | #122 | diff --git a/docs/writing_a_material.md b/docs/writing_a_material.md index 1fab5a8..bf004a5 100644 --- a/docs/writing_a_material.md +++ b/docs/writing_a_material.md @@ -156,6 +156,24 @@ Without `.stateful()`, a stepper that reads nothing would run once. A recorder o The engine then runs the solver after everything the driven callbacks read globally (the time step, the temperature), whatever order the deck creates them in. Reads of the solver's own outputs (a function reading the current unknown) add no edge: the solver fills those itself. - **What a driven callback must be:** a repeatable trial evaluation. Every call computes from the same committed history and its current inputs, and writes outputs only. It must not commit, count or record, because the solver calls it once per iteration. - **Driving any output of a group:** an output filled by another output's callback reruns that callback through `update_source()`. +- **A fixed point over a local edge:** a node that reads another over a local edge sees the previous pass. A material that measures how far two passes disagree declares that output with `converges()`, with the tolerance and the pass limit as parameters: + + ```cpp + base::converges(m_residual, base::template get_parameter("tolerance"), base::template get_parameter("max_passes")); + ``` + + At `finalize()` the engine forms an iteration group: the callbacks downstream of the consumer of a local edge and upstream of the residual. `update()` runs: + - first, the callbacks the group reads, once; + - then the group, until |residual| ≤ tolerance (a NaN never converges); + - last, the rest, once. A recorder or a plot therefore writes one row per update. + + After `max_passes` passes, `update()` throws `cutback_error` and names the residual. `update_property()` on a member runs its upstream once, without iterating. These rules throw at `finalize()` and name the node: + - `converges()` with no local edge behind the residual; + - a `.stateful()` callback inside a group; + - a second residual over the callbacks of a group; + - a tolerance below zero or `max_passes` below one. + + `mean_field_composite_local_nonlinear` declares its `residual` this way. - **`material_ref`:** a material reached through `material_ref` (a `local_newton` solver object) is not an input. The engine: - treats what the target reads as reads of every callback of the referencing material; - keeps those callbacks varying unless every output of the target is constant (a target with history, state or a bare output never is); diff --git a/include/numsim-materials/core/derivatives.h b/include/numsim-materials/core/derivatives.h new file mode 100644 index 0000000..1dc3b85 --- /dev/null +++ b/include/numsim-materials/core/derivatives.h @@ -0,0 +1,333 @@ +#ifndef NUMSIM_MATERIALS_CORE_DERIVATIVES_H +#define NUMSIM_MATERIALS_CORE_DERIVATIVES_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/history_property.h" +#include "numsim-materials/core/input_types.h" +#include "numsim-materials/core/property.h" +#include "numsim-materials/core/property_traits.h" + +/// Partial derivatives of a callback's outputs with respect to its reads (#217), packed into Mandel blocks. +/// +/// Every index pair maps to one Mandel component with the weight sqrt(2) off the diagonal, the layout of +/// solvers/unknown_layout.h. The weights make A : B = mandel(A) . mandel(B), so the chain rule over any +/// middle quantity (scalar, tensor2, tensor4) is a plain matrix product of the packed blocks. +namespace numsim::materials::derivatives { + +/// Rank of a value the graph differentiates: 0 for a scalar. +template +constexpr std::size_t rank_of() { + if constexpr (std::is_arithmetic_v) + return 0; + else + return V::rank(); +} + +/// The value behind an output or a read: T of history_property, input_property and input_history. +template +struct value_of { + using type = A; +}; +template +struct value_of> { + using type = T; +}; +template +struct value_of> { + using type = T; +}; +template +struct value_of> { + using type = T; +}; +template +using value_of_t = typename value_of::type; + +/// A value with a derivative: a scalar, or a symmetric tensor of rank 2 or 4 (minor symmetric). +template +concept node_value = std::is_same_v || std::is_same_v> || + std::is_same_v>; + +/// The type a derivative is read as: node_value's scalar, or a tensor of rank 2 to 8 (d tensor4 / d tensor4), of the same T and Dim. +template +concept derivative_value = node_value || std::is_same_v> || + std::is_same_v>; + +template +struct partial_value { + using type = tmech::tensor; +}; +template +struct partial_value { + using type = T; +}; +/// d out / d in: rank(out) + rank(in), a scalar for two scalars. +template +using partial_t = typename partial_value() + rank_of()>::type; + +/// Mandel components of a symmetric index pair. +template +inline constexpr std::size_t pair_width{Dim * (Dim + 1) / 2}; + +/// Mandel components of a value of rank `rank`, one Mandel index per index pair. +template +constexpr std::size_t width(std::size_t rank) { + std::size_t w{1}; + for (std::size_t pair{0}; pair < rank / 2; ++pair) + w *= pair_width; + return w; +} + +namespace detail { +/// The index pairs in tmech's Mandel order: 11, 22, 33, 23, 13, 12 (2D: 11, 22, 12). +template +constexpr auto pairs() { + if constexpr (Dim == 2) + return std::array, 3>{{{0, 0}, {1, 1}, {0, 1}}}; + else + return std::array, 6>{{{0, 0}, {1, 1}, {2, 2}, {1, 2}, {0, 2}, {0, 1}}}; +} + +/// The Mandel component of the index pair (i, j). +template +constexpr std::size_t component(std::size_t i, std::size_t j) { + if (i == j) + return i; + if constexpr (Dim == 2) + return 2; + else + return 6 - i - j; +} + +template +constexpr decltype(auto) entry(V& v, std::array const& idx, std::index_sequence) { + return v(idx[I]...); +} + +constexpr std::size_t power(std::size_t base, std::size_t exponent) { + std::size_t p{1}; + while (exponent-- > 0) + p *= base; + return p; +} +} // namespace detail + +/// Packs P = d out / d in (rank(out) = OutRank) into a width(OutRank) x width(rank(P) - OutRank) column-major block; +/// each index pair is projected onto its symmetric part. Index loops: packing is the boundary to the matrices. +template +void pack(P const& p, T* block) { + if constexpr (std::is_arithmetic_v

) { + *block = p; + } else { + constexpr std::size_t rank{P::rank()}, n{rank / 2}, W{pair_width}; + constexpr std::size_t rows{width(OutRank)}, cols{width(rank - OutRank)}; + constexpr auto pair{detail::pairs()}; + std::array m{}; + std::array idx{}; + for (std::size_t flat{0}; flat < rows * cols; ++flat) { + for (std::size_t k{n}, rest{flat}; k-- > 0; rest /= W) + m[k] = rest % W; // the first pair varies slowest + T sum{0}; + for (std::size_t swaps{0}; swaps < (std::size_t{1} << n); ++swaps) { + for (std::size_t k{0}; k < n; ++k) { + auto const flip{(swaps >> k) & 1}; + idx[2 * k] = pair[m[k]][flip]; + idx[2 * k + 1] = pair[m[k]][1 - flip]; + } + sum += detail::entry(p, idx, std::make_index_sequence{}); + } + T weight{T{1} / T(std::size_t{1} << n)}; + for (auto const mk : m) + if (mk >= Dim) + weight *= std::numbers::sqrt2_v; + block[flat / cols + (flat % cols) * rows] = weight * sum; + } + } +} + +/// The inverse of pack(): a minor-symmetric tensor (or a scalar) from its column-major Mandel block. +template +void unpack(T const* block, P& p) { + if constexpr (std::is_arithmetic_v

) { + p = *block; + } else { + constexpr std::size_t rank{P::rank()}, W{pair_width}; + constexpr std::size_t rows{width(OutRank)}, cols{width(rank - OutRank)}; + std::array idx{}; + for (std::size_t flat{0}; flat < detail::power(Dim, rank); ++flat) { + for (std::size_t k{rank}, rest{flat}; k-- > 0; rest /= Dim) + idx[k] = rest % Dim; + std::size_t mandel{0}; + T weight{1}; + for (std::size_t k{0}; k < rank; k += 2) { + auto const c{detail::component(idx[k], idx[k + 1])}; + mandel = mandel * W + c; + if (c >= Dim) + weight *= std::numbers::sqrt2_v; + } + detail::entry(p, idx, std::make_index_sequence{}) = block[mandel / cols + (mandel % cols) * rows] / weight; + } + } +} + +/// The value of `source` as V: a property's value or a history property's new value. +template +V& value_in(property_base& source) { + if (auto* plain{dynamic_cast*>(&source)}) + return plain->get(); + if (auto* history{dynamic_cast*>(&source)}) + return history->new_value(); + throw std::logic_error("derivatives: '" + source.traits().id.owner + "::" + source.traits().id.name + + "' holds another type than its reader"); +} + +/// The rank of the value `prop` holds when it has a derivative (node_value), else nothing. +template +std::optional property_rank(property_base const& prop) { + auto const holds{[&](std::type_identity) { + return dynamic_cast const*>(&prop) || dynamic_cast const*>(&prop); + }}; + if (holds(std::type_identity{})) + return 0; + if (holds(std::type_identity>{})) + return 2; + if (holds(std::type_identity>{})) + return 4; + return std::nullopt; +} + +/// A read the check perturbs one Mandel component at a time: its value is held by a property, given as `source`. +template +class read_probe { +public: + virtual ~read_probe() = default; + /// Mandel components of the read. + virtual std::size_t width() const noexcept = 0; + /// Mandel component k of the read's value, held by `source`. + virtual T read_component(property_base& source, std::size_t k) const = 0; + /// Adds h to Mandel component k of the read's value; the returned function restores the value. + virtual std::function perturb(property_base& source, std::size_t k, T h) const = 0; +}; + +template +class typed_read : public read_probe { +public: + std::size_t width() const noexcept override { return derivatives::width(rank_of()); } + + T read_component(property_base& source, std::size_t k) const override { + std::array(rank_of())> packed{}; + derivatives::pack()>(value_in(source), packed.data()); + return packed[k]; + } + + std::function perturb(property_base& source, std::size_t k, T h) const override { + auto& v{value_in(source)}; + std::function restore{[&v, old = v] { v = old; }}; + if constexpr (std::is_arithmetic_v) { + v += h; + } else { + std::array(rank_of())> unit{}; + unit[k] = h; + V step; + derivatives::unpack()>(unit.data(), step); + v = V{v + step}; + } + return restore; + } +}; + +/// One declared partial d output / d read of a callback: the value the callback fills, packed on demand. +template +class partial_block : public read_probe { +public: + partial_block(std::size_t rows, std::size_t cols) noexcept : m_rows(rows), m_cols(cols) {} + partial_block(partial_block const&) = delete; + partial_block& operator=(partial_block const&) = delete; + + std::size_t rows() const noexcept { return m_rows; } + std::size_t cols() const noexcept { return m_cols; } + std::size_t width() const noexcept override { return m_cols; } + /// The partial as a rows() x cols() column-major Mandel block. + virtual void pack(T* block) const = 0; + /// The output's value, rows() Mandel components. + virtual void pack_output(T* value) const = 0; + +private: + std::size_t m_rows, m_cols; +}; + +template +class typed_partial final : public partial_block { +public: + using output_value = value_of_t; + using read_value = value_of_t; + static constexpr std::size_t out_rank{rank_of()}; + + explicit typed_partial(Out const& output) + : partial_block{derivatives::width(out_rank), derivatives::width(rank_of())}, m_output(output) {} + + partial_t value{}; + + void pack(T* block) const override { derivatives::pack(value, block); } + + void pack_output(T* out) const override { + if constexpr (requires(Out const& o) { o.new_value(); }) + derivatives::pack(m_output.new_value(), out); + else + derivatives::pack(m_output, out); + } + + T read_component(property_base& source, std::size_t k) const override { return m_read.read_component(source, k); } + + std::function perturb(property_base& source, std::size_t k, T h) const override { return m_read.perturb(source, k, h); } + +private: + Out const& m_output; + typed_read m_read; +}; + +/// A requested derivative as the host reads it, unpacked from its accumulated block. +template +class derivative_slot_base { +public: + virtual ~derivative_slot_base() = default; + virtual void unpack(T const* block) = 0; +}; + +template +class derivative_slot final : public derivative_slot_base { +public: + explicit derivative_slot(std::size_t out_rank) : m_out_rank(out_rank) {} + + R value{}; + + void unpack(T const* block) override { + switch (m_out_rank) { + case 0: derivatives::unpack(block, value); break; + case 2: + if constexpr (rank_of() >= 2) + derivatives::unpack(block, value); + break; + case 4: + if constexpr (rank_of() >= 4) + derivatives::unpack(block, value); + break; + default: break; + } + } + +private: + std::size_t m_out_rank; +}; + +} // namespace numsim::materials::derivatives + +#endif // NUMSIM_MATERIALS_CORE_DERIVATIVES_H diff --git a/include/numsim-materials/core/material_base.h b/include/numsim-materials/core/material_base.h index 2ba1839..0109bac 100755 --- a/include/numsim-materials/core/material_base.h +++ b/include/numsim-materials/core/material_base.h @@ -169,8 +169,36 @@ class material_base : public material_interface { return *this; } + /// The callback fills d output / d read for each of `reads` (as in reads()) through partial() while + /// want_partials() is true. After reads(). + template + callback_handle& partials(Output const& output, Reads const&... reads) { + (each_read(reads, [&](auto const& read) { m_owner.declare_partial(m_outputs.front(), output, read); }), ...); + return *this; + } + + /// The outputs of the callback do not depend on `reads`: their partials are zero. After reads(). + template + callback_handle& independent_of(Reads const&... reads) { + (each_read(reads, [&](auto const& read) { m_owner.declare_independent(m_outputs.front(), read); }), ...); + return *this; + } + private: friend class material_base; + template + static void each_read(A const& read, F const& f) { + if constexpr (std::is_pointer_v) { + if (read) + f(*read); + } else if constexpr (detail::pointer_range) { + for (auto const* p : read) + if (p) + f(*p); + } else { + f(read); + } + } template static void sort_read(std::vector& inputs, std::vector& outputs, A const& read) { if constexpr (std::derived_from, input_wire_base>) @@ -202,6 +230,16 @@ class material_base : public material_interface { return handle; } + /// update() reruns the callbacks downstream of a local edge and upstream of `residual` (an output) until + /// |residual| <= tolerance, and throws cutback_error after max_passes. Both are references to parameters. + void converges(value_type const& residual, value_type const& tolerance, int const& max_passes) { + base::declare_converges(residual, tolerance, max_passes); + } + // the declaration keeps pointers to tolerance and max_passes: a temporary would dangle + void converges(value_type const&, value_type&&, int const&) = delete; + void converges(value_type const&, value_type const&, int&&) = delete; + void converges(value_type const&, value_type&&, int&&) = delete; + /// `output` has no callback of its own: the callback of `writer` fills it. Lets the engine skip a constant /// callback while another callback of the material keeps running. template diff --git a/include/numsim-materials/core/material_context.h b/include/numsim-materials/core/material_context.h index 78b97f3..cbd0311 100644 --- a/include/numsim-materials/core/material_context.h +++ b/include/numsim-materials/core/material_context.h @@ -82,9 +82,39 @@ class material_context { // --- Execution (delegated to property_engine) --- + /// Runs the callbacks, then accumulates the requested derivatives (derivative()). void update() { check_finalized("update"); m_engine.update(); + m_engine.differentiate(); + } + + /// d material::property / d seed_material::seed_property as R (rank of both added, a scalar for two scalars), + /// filled by every update() after the values. A requested path through a callback without its partial throws. + /// R is the value type or a tensor of rank 2 to 8 of the context's dimension. + template + requires derivatives::derivative_value + const R& derivative(const std::string& material, const std::string& property, const std::string& seed_material, + const std::string& seed_property) { + check_finalized("derivative"); + auto const* target{find_property(material, property)}; + auto const* seed{find_property(seed_material, seed_property)}; + if (!target || !seed) + throw std::runtime_error("derivative(): property '" + (target ? seed_material + "::" + seed_property : material + "::" + property) + + "' not found"); + return m_engine.template derivative(target, seed); + } + + /// Whether the seed is upstream of the property: otherwise derivative() is zero. Throws on an unknown name. + bool has_path(const std::string& material, const std::string& property, const std::string& seed_material, + const std::string& seed_property) { + check_finalized("has_path"); + auto const* target{find_property(material, property)}; + auto const* seed{find_property(seed_material, seed_property)}; + if (!target || !seed) + throw std::runtime_error("has_path(): property '" + (target ? seed_material + "::" + seed_property : material + "::" + property) + + "' not found"); + return m_engine.has_path(target, seed); } void update_property(const std::string& material, const std::string& property, diff --git a/include/numsim-materials/core/material_interface.h b/include/numsim-materials/core/material_interface.h index d6aed65..5149344 100644 --- a/include/numsim-materials/core/material_interface.h +++ b/include/numsim-materials/core/material_interface.h @@ -20,6 +20,7 @@ #include #include +#include "numsim-materials/core/derivatives.h" #include "numsim-materials/core/property_traits.h" #include "numsim-materials/core/property.h" #include "numsim-materials/core/history_property.h" @@ -27,6 +28,7 @@ #include "numsim-materials/core/property_registry_interface.h" #include "numsim-materials/core/material_ref.h" #include "numsim-materials/core/traits.h" +#include "numsim-materials/core/unknown_spec.h" namespace numsim::materials { @@ -61,6 +63,22 @@ bool same_value(T const& a, T const& b) { } } +/// How far `b` is from `a` relative to the larger of the two: 0 when equal, infinite when they differ and one is not finite. +template +double relative_deviation(T const& a, T const& b) { + double distance, size; + if constexpr (std::is_arithmetic_v) { + distance = std::abs(static_cast(a) - static_cast(b)); + size = std::max(std::abs(static_cast(a)), std::abs(static_cast(b))); + } else if constexpr (requires { tmech::norm(a - b); }) { + distance = static_cast(tmech::norm(a - b)); + size = static_cast(std::max(tmech::norm(a), tmech::norm(b))); + } else { + return a == b ? 0 : 1; + } + return distance == 0 ? 0 : std::isfinite(distance) && std::isfinite(size) ? distance / size : std::numeric_limits::infinity(); +} + /// Whether a value is finite; types without a norm are. template bool is_finite(T const& a) { @@ -76,6 +94,13 @@ bool is_finite(T const& a) { enum class verdict { same, changed, not_finite }; } // namespace detail +/// What the engine reports back to a residual system (solvers/system.h). +struct system_report { + int iterations{0}; ///< updates of the unknowns in the last solve + std::vector> zero_blocks; ///< (row, column) without a path, from finalize() + bool planned{false}; ///< finalize() has filled zero_blocks +}; + /// Virtual base for all materials: property registration, typed input wiring and material references; /// work happens in the callbacks bound to properties, there is no virtual update() to override. template @@ -209,6 +234,7 @@ class material_interface { using property_base_type = numsim_core::property_base; /// Copies a value; the returned check compares it with the current one. Empty for types never constant. using snapshot = std::function()>; + using deviation_probe = std::function()>; /// An output declared a function of the parameters, of `reads` (all inputs when unset) and of the outputs /// `own_reads` of this material (material_base::constant_output, callback_handle::reads). @@ -226,9 +252,24 @@ class material_interface { std::vector inputs; std::vector outputs; std::vector drives; ///< local inputs whose source callbacks it reruns + std::vector own_reads; ///< outputs of this material it reads }; const std::vector& declared_reads() const noexcept { return m_declared_reads; } + /// A residual system update() solves (solvers/system.h): the unknowns, outputs of this material, until the + /// residuals reach |R|_inf <= *tolerance. converges() declares one without unknowns. + struct system_declaration { + system_method method; + std::vector unknowns; ///< in column order + std::vector seeds; ///< start values of the scalar unknowns + std::vector residuals; ///< row i is solved for unknown i + property_base_type const* own_residual; ///< converges(): the output iterated to, max_iter counts passes + value_type const* tolerance; + int const* max_iter; + system_report* report; ///< null for converges() + }; + const std::vector& systems() const noexcept { return m_systems; } + /// The names of the materials this one reaches through material_ref. auto material_ref_targets() const { return m_material_refs | std::views::transform([](auto const& ref) -> const std::string& { return ref->target_name(); }); @@ -267,7 +308,126 @@ class material_interface { /// The engine's revision counter, bumped by set_parameter(); bound at finalize(). void bind_revision(std::size_t* revision) noexcept { m_revision = revision; } + /// d output / d read declared by a callback (callback_handle::partials()); the read is an input or an output of + /// this material, the block holds what partial() returns. + struct partial_declaration { + property_base_type const* callback; + property_base_type const* output; + input_wire_base const* input; ///< the read when it is an input ... + property_base_type const* own; ///< ... else the output of this material it is + void const* output_value; ///< the addresses partial() is called with + void const* read_value; + std::unique_ptr> block; + }; + const std::vector& partials() const noexcept { return m_partials; } + + /// A read whose partials are declared zero (callback_handle::independent_of()), as in partial_declaration; + /// `probe` perturbs the read for check_partials(), empty for a read that is not a scalar or tensor. + struct independent_declaration { + property_base_type const* callback; + input_wire_base const* input; + property_base_type const* own; + std::unique_ptr> probe; + }; + const std::vector& independent_reads() const noexcept { return m_independent; } + + /// The worst relative deviation of an output of a callback from a copy taken earlier. + struct output_deviation { + property_base_type const* output{nullptr}; + double relative{0}; + }; + /// Copies the outputs of the callback filling `callback` (with declared reads); the returned function reports the + /// largest relative deviation of any from its copy. Outputs of types without a comparison are not watched. + std::function watch_outputs(property_base_type const* callback) const { + auto const it{std::ranges::find(m_declared_reads, callback, &read_declaration::callback)}; + if (it == m_declared_reads.end()) + throw std::logic_error("watch_outputs(): the callback has no declared reads"); + std::vector>> watched; + for (auto const* output : it->outputs) + for (auto const& [address, entry] : m_outputs) + if (entry.property == output && entry.deviation) + watched.emplace_back(output, entry.deviation()); + return [watched = std::move(watched)] { + output_deviation worst; + for (auto const& [output, deviation] : watched) + if (auto const relative{deviation()}; relative > worst.relative || !worst.output) + worst = {output, relative}; + return worst; + }; + } + + /// Whether `callback` fills its partials when it runs: the engine sets it for the callbacks on a requested derivative. + void request_partials(property_base_type const* callback, bool on) { + std::erase(m_partials_requested, callback); + if (on) + m_partials_requested.push_back(callback); + } + bool partials_requested(property_base_type const* callback) const noexcept { + return std::ranges::contains(m_partials_requested, callback); + } + + /// Runs `run`, the update of `callback`, with want_partials() as requested for it; the engine wraps those callbacks. + template + void run_callback(property_base_type const* callback, Run const& run) { + struct restore { + bool& flag; + bool outer; + ~restore() { flag = outer; } + } const guard{m_want_partials, std::exchange(m_want_partials, partials_requested(callback))}; + run(); + } + protected: + /// True while the running callback is to fill its declared partials: a requested derivative passes through it. + bool want_partials() const noexcept { return m_want_partials; } + + /// The partial d output / d read the callback fills, declared by partials(): rank(output) + rank(read). + template + auto& partial(Out const& output, In const& read) { + using block_type = derivatives::typed_partial; + auto const it{std::ranges::find_if(m_partials, [&](auto const& p) { + return p.output_value == static_cast(&output) && p.read_value == static_cast(&read); + })}; + auto* block{it != m_partials.end() ? dynamic_cast(it->block.get()) : nullptr}; + if (!block) + throw std::logic_error("partial(): material '" + m_name + "' fills a partial its callback did not declare in partials()"); + return block->value; + } + + /// The callback filling `callback` (after its reads()) fills d output / d read through partial(). + template + void declare_partial(void const* callback, Out const& output, In const& read) { + static_assert(derivatives::node_value, value_type, Traits::Dim> && + derivatives::node_value, value_type, Traits::Dim>, + "partials(): outputs and reads are scalars or symmetric tensors of rank 2 or 4"); + auto const& [d, entry]{declared_read_of("partials", callback)}; + auto const& out{output_entry_of("partials", &output)}; + if (!std::ranges::contains(d->outputs, out.property)) + throw std::logic_error("partials(): '" + key_of(out) + "' is not an output of the callback of '" + key_of(entry) + "'"); + auto const [input, own]{read_of("partials", *d, entry, read)}; + if (std::ranges::any_of(m_partials, [&](auto const& p) { return p.output == out.property && p.input == input && p.own == own; })) + throw std::logic_error("partials(): '" + key_of(out) + "' declares a partial twice"); + if (is_independent(d->callback, input, own)) + throw std::logic_error("partials(): '" + key_of(entry) + "' declares a read both in partials() and independent_of()"); + m_partials.push_back({d->callback, out.property, input, own, &output, &read, + std::make_unique>(output)}); + } + + /// The outputs of the callback filling `callback` do not depend on `read`: its partials are zero. It holds for every + /// output of the callback, and check_partials() perturbs the read to verify it. + template + void declare_independent(void const* callback, In const& read) { + auto const& [d, entry]{declared_read_of("independent_of", callback)}; + auto const [input, own]{read_of("independent_of", *d, entry, read)}; + if (is_independent(d->callback, input, own) || + std::ranges::any_of(m_partials, [&](auto const& p) { return p.callback == d->callback && p.input == input && p.own == own; })) + throw std::logic_error("independent_of(): '" + key_of(entry) + "' declares a read twice or also in partials()"); + std::unique_ptr> probe; + if constexpr (derivatives::node_value, value_type, Traits::Dim>) + probe = std::make_unique>>(); + m_independent.push_back({d->callback, input, own, std::move(probe)}); + } + template T& add_property(std::string const& property_name, T&& data = T()) { T& value{m_property_handler.template add_property(name(), property_name, std::forward(data))}; @@ -279,14 +439,22 @@ class material_interface { return !detail::is_finite(value) ? not_finite : detail::same_value(copy, value) ? same : changed; }}; }; - m_outputs[&value] = {*m_property_handler.find(name(), property_name), std::move(take)}; + deviation_probe deviation; + if constexpr (detail::verifiable) + deviation = [&value] { return std::function{[&value, copy = value] { return detail::relative_deviation(copy, value); }}; }; + m_outputs[&value] = {*m_property_handler.find(name(), property_name), std::move(take), std::move(deviation)}; return value; } template history_property& add_history(std::string const& property_name, T&& data = T()) { auto& value{m_property_handler.template add_history(name(), property_name, std::forward(data))}; - m_outputs[&value] = {*m_property_handler.find(name(), property_name), {}}; // never constant + deviation_probe deviation; + if constexpr (detail::verifiable) + deviation = [&value] { + return std::function{[&value, copy = value.new_value()] { return detail::relative_deviation(copy, value.new_value()); }}; + }; + m_outputs[&value] = {*m_property_handler.find(name(), property_name), {}, std::move(deviation)}; // never constant return value; } @@ -350,7 +518,8 @@ class material_interface { auto const entries{outputs | std::views::transform([this](void const* output) { return &output_entry_of("reads", output); }) | std::ranges::to()}; m_declared_reads.push_back({callback.property, inputs, - entries | std::views::transform(&output_entry::property) | std::ranges::to(), {}}); + entries | std::views::transform(&output_entry::property) | std::ranges::to(), {}, + own_reads}); if (verifiable && std::ranges::none_of(entries, [this](auto const* entry) { return entry->property->is_history() || std::ranges::contains(m_stateful, entry->property); })) @@ -385,6 +554,29 @@ class material_interface { } } + /// `residual` is the value add_property() returned; `tolerance` and `max_passes` are read on every update. + void declare_converges(value_type const& residual, value_type const& tolerance, int const& max_passes) { + auto const& entry{output_entry_of("converges", &residual)}; + if (std::ranges::contains(m_systems, entry.property, &system_declaration::own_residual)) + throw std::logic_error("converges(): '" + key_of(entry) + "' is declared twice"); + m_systems.push_back({system_method::picard, {}, {}, {}, entry.property, &tolerance, &max_passes, nullptr}); + } + + /// `unknowns` are values add_property() or add_history() returned, `residuals` inputs of this material. + void declare_system(system_method method, std::vector const& unknowns, std::vector seeds, + std::vector residuals, value_type const& tolerance, int const& max_iter, + system_report& report) { + if (!std::ranges::all_of(residuals, [this](auto const* r) { return std::ranges::contains(m_typed_inputs, r, &std::unique_ptr::get); })) + throw std::logic_error("system(): material '" + m_name + "' names a residual that is not one of its inputs"); + std::vector properties; + for (auto const* unknown : unknowns) + properties.push_back(output_entry_of("system", unknown).property); + m_systems.push_back({method, std::move(properties), std::move(seeds), std::move(residuals), nullptr, &tolerance, &max_iter, &report}); + } + void declare_converges(value_type const&, value_type&&, int const&) = delete; + void declare_converges(value_type const&, value_type const&, int&&) = delete; + void declare_converges(value_type const&, value_type&&, int&&) = delete; + void declare_written_by(void const* output, void const* writer) { auto const& out{output_entry_of("written_by", output)}; auto const& by{output_entry_of("written_by", writer)}; @@ -432,12 +624,14 @@ class material_interface { struct output_entry { property_base_type const* property; snapshot take; + deviation_probe deviation; ///< copies the value; the returned function gives its relative deviation from the copy }; std::unordered_map m_outputs; ///< value address -> property std::vector m_constant_outputs; std::vector m_read_declared; ///< callbacks with declared reads std::vector m_declared_reads; std::vector m_stateful; ///< outputs never constant + std::vector m_systems; std::vector> m_written_by, m_callback_order; /// The entry of an output of this material; declarations end at finalize(), when the engine reads them. @@ -466,6 +660,35 @@ class material_interface { bool has_writer(output_entry const& entry) const { return std::ranges::contains(m_written_by, entry.property, &decltype(m_written_by)::value_type::first); } + /// The read declaration of the callback filling `callback`, which partials() and independent_of() follow. + std::pair declared_read_of(char const* what, void const* callback) { + auto const& entry{output_entry_of(what, callback)}; + auto const it{std::ranges::find(m_declared_reads, entry.property, &read_declaration::callback)}; + if (it == m_declared_reads.end()) + throw std::logic_error(std::string{what} + "(): '" + key_of(entry) + "' declares partials before its reads"); + return {&*it, entry}; + } + /// A read of `d` as (input, own output); throws when the callback does not read it. + template + std::pair read_of(char const* what, read_declaration const& d, + output_entry const& entry, In const& read) const { + if constexpr (std::derived_from) { + if (std::ranges::contains(d.inputs, static_cast(&read))) + return {&read, nullptr}; + } else if (auto const it{m_outputs.find(&read)}; + it != m_outputs.end() && std::ranges::contains(d.own_reads, it->second.property)) { + return {nullptr, it->second.property}; + } + throw std::logic_error(std::string{what} + "(): '" + key_of(entry) + "' names a read its callback does not declare in reads()"); + } + bool is_independent(property_base_type const* callback, input_wire_base const* input, property_base_type const* own) const { + return std::ranges::any_of(m_independent, [&](auto const& i) { return i.callback == callback && i.input == input && i.own == own; }); + } + + std::vector m_partials; + std::vector m_independent; + std::vector m_partials_requested; ///< callbacks filling their partials + bool m_want_partials{false}; std::size_t* m_revision{nullptr}; std::vector m_reads_all; ///< callbacks reading every input (declare_reads_all) }; diff --git a/include/numsim-materials/core/partial_check.h b/include/numsim-materials/core/partial_check.h new file mode 100644 index 0000000..d90678a --- /dev/null +++ b/include/numsim-materials/core/partial_check.h @@ -0,0 +1,159 @@ +#ifndef NUMSIM_MATERIALS_CORE_PARTIAL_CHECK_H +#define NUMSIM_MATERIALS_CORE_PARTIAL_CHECK_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/material_context.h" + +namespace numsim::materials { + +namespace partial_check_detail { +/// Runs `f` when the scope ends, also while an exception unwinds it; what `f` throws is dropped. +template +struct on_exit { + F f; + ~on_exit() { + try { + f(); + } catch (...) { + } + } +}; +template +on_exit(F) -> on_exit; +} // namespace partial_check_detail + +/// Compares every partial `material` declares with central differences of its own callback, one Mandel component +/// of the read at a time, at the state the context holds (after an update()), and perturbs every read declared +/// independent_of() to see that no output of the callback moves. Returns one line per wrong declaration; the +/// outputs and reads are restored, also when a perturbed callback throws. +template +std::vector check_partials(material_context& ctx, std::string const& material, + typename Traits::value_type tolerance = 1e-6) { + using T = typename Traits::value_type; + using partial_check_detail::on_exit; + constexpr T eps{std::numeric_limits::epsilon()}; + constexpr T scale_floor{1e-3}; // the scale of a read that is zero, or nearly + constexpr T independence_step{1e-3}; + auto* mat{ctx.find(material)}; + if (!mat) + throw std::invalid_argument("check_partials(): no material '" + material + "'"); + + auto const run_of{[&](auto const* callback) { + auto* runs{ctx.find_property(material, callback->traits().id.name)}; + return [runs] { runs->traits().update(); }; + }}; + auto const source_of{[&](auto const& d) { + return d.input ? ctx.find_property(d.input->source_owner(), d.input->source_name()) + : ctx.find_property(material, d.own->traits().id.name); + }}; + auto const key_of{[](auto const* source) { return source->traits().id.owner + "::" + source->traits().id.name; }}; + // the largest component of the read: one step for all of them + auto const read_scale{[&](auto const& read, auto& source) { + T scale{scale_floor}; + for (std::size_t k{0}; k < read.width(); ++k) + scale = std::max(scale, std::abs(read.read_component(source, k))); + return scale; + }}; + + std::vector mismatches; + auto const& partials{mat->partials()}; + std::vector callbacks; + for (auto const& p : partials) + if (!std::ranges::contains(callbacks, p.callback)) + callbacks.push_back(p.callback); + + for (auto const* callback : callbacks) { + auto const run{run_of(callback)}; + on_exit rerun{run}; // the outputs at the state the check found + bool const requested{mat->partials_requested(callback)}; + mat->request_partials(callback, true); + { + on_exit reset{[&] { mat->request_partials(callback, requested); }}; + run(); + } + auto of_callback{partials | std::views::filter([callback](auto const& q) { return q.callback == callback; })}; + // packed before any perturbed rerun overwrites them + std::vector> packed; + for (auto const& p : of_callback) + p.block->pack(packed.emplace_back(p.block->rows() * p.block->cols()).data()); + for (auto const& [p, analytic] : std::views::zip(of_callback, packed)) { + auto const& block{*p.block}; + std::vector plus(block.rows()), minus(block.rows()); + auto* source{source_of(p)}; + auto const read{key_of(source)}; + T scale{0}; + for (auto const a : analytic) + scale = std::max(scale, std::abs(a)); + T const read_size{read_scale(block, *source)}; + T const h{std::sqrt(eps) * read_size}; + for (std::size_t k{0}; k < block.cols(); ++k) { + for (auto [step, values] : {std::pair{h, &plus}, std::pair{-h, &minus}}) { + on_exit restore{block.perturb(*source, k, step)}; + run(); + block.pack_output(values->data()); + } + T error{0}, size{scale}, output{0}; + for (std::size_t r{0}; r < block.rows(); ++r) { + T const fd{(plus[r] - minus[r]) / (2 * h)}; + error = std::max(error, std::abs(fd - analytic[r + k * block.rows()])); + size = std::max(size, std::abs(fd)); + output = std::max({output, std::abs(plus[r]), std::abs(minus[r])}); + } + // Rounding of the differences, for an output no larger than the analytic block can produce: a large + // output must not hide a small partial. + T const allowance{10 * eps * std::min(output, scale * (read_size + h)) / h}; + if (error > tolerance * size + allowance) { + mismatches.push_back(std::format("'{}::{}' wrt '{}': Mandel column {} differs from central differences by {:g} " + "(largest entry {:g})", + material, p.output->traits().id.name, read, k, error, size)); + break; + } + } + } + } + + std::vector independents; + for (auto const& i : mat->independent_reads()) + if (i.probe && !std::ranges::contains(independents, i.callback)) + independents.push_back(i.callback); + for (auto const* callback : independents) { + auto const run{run_of(callback)}; + on_exit rerun{run}; + for (auto const& i : mat->independent_reads()) { + if (i.callback != callback || !i.probe) + continue; + auto* source{source_of(i)}; + auto const read{key_of(source)}; + T const h{independence_step * read_scale(*i.probe, *source)}; + run(); // the outputs at the state the check found + auto const deviation{mat->watch_outputs(callback)}; + bool moved{false}; + for (std::size_t k{0}; k < i.probe->width() && !moved; ++k) + for (T const step : {h, -h}) { + on_exit restore{i.probe->perturb(*source, k, step)}; + run(); + auto const worst{deviation()}; + if (worst.relative > tolerance) { + mismatches.push_back(std::format("'{}::{}' is declared independent_of '{}', but Mandel component {} of the read " + "changes it by {:g} (relative)", + material, worst.output->traits().id.name, read, k, worst.relative)); + moved = true; + break; + } + } + } + } + return mismatches; +} + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_CORE_PARTIAL_CHECK_H diff --git a/include/numsim-materials/core/property_engine.h b/include/numsim-materials/core/property_engine.h index 5d38990..3764bef 100644 --- a/include/numsim-materials/core/property_engine.h +++ b/include/numsim-materials/core/property_engine.h @@ -2,9 +2,12 @@ #define PROPERTY_ENGINE_H #include +#include #include +#include #include #include +#include #include #include #include @@ -19,6 +22,9 @@ #include #include #include +#include +#include "numsim-materials/core/derivatives.h" +#include "numsim-materials/core/errors.h" #include "numsim-materials/core/property_traits.h" #include "numsim-materials/core/property.h" #include "numsim-materials/core/history_property.h" @@ -84,21 +90,26 @@ class property_engine { for (auto const* mat : interfaces) m_stateful.insert(mat->stateful_outputs().begin(), mat->stateful_outputs().end()); collect_history_properties(properties); + form_groups(interfaces); find_constants(interfaces, writers); if constexpr (NUMSIM_MATERIALS_DEBUG_CHECKS) track_reads(interfaces); route_to_writers(interfaces); + index_partials(interfaces); + plan_systems(); } - /// Run the callbacks still to run, in topological order: a constant one only until it has run. + /// Run the callbacks still to run, in topological order: a constant one only until it has run. A residual system + /// (solvers/system.h, material_base::converges) is solved between what it reads and the rest. void update() { guarded("update", [this] { refresh_constants(); if (m_verify) std::ranges::for_each(m_retired, [this](property_base* prop) { verify(prop); }); - for (auto* prop : m_active) - prop->traits().update(); + for (auto const& s : m_stages) + s.system ? solve(*s.system, s.callbacks) : run(s.callbacks); retire([](property_base*) { return true; }); + m_evaluated = true; }); } @@ -135,8 +146,15 @@ class property_engine { if (m_verify) for (auto* prop : m_retired | std::views::filter(runs)) verify(prop); - for (auto* prop : m_active | std::views::filter(runs)) - prop->traits().update(); + for (auto const& s : m_stages) { + if (!s.system || s.system->unknowns.empty()) { // a converges() group runs once, as a member asked for directly + for (auto* prop : s.callbacks | std::views::filter(runs)) + prop->traits().update(); + } else if (std::ranges::any_of(s.callbacks, runs) || std::ranges::any_of(s.system->unknowns, runs)) { + auto const members{s.callbacks | std::views::filter([&](auto* prop) { return !exclude.contains(prop); }) | std::ranges::to()}; + solve(*s.system, members); + } + } retire(runs); }); } @@ -234,6 +252,61 @@ class property_engine { return g; } + /// d target / d seed as R (rank(target) + rank(seed), a scalar for two scalars), filled by differentiate(). + /// Resolved once: a path from the seed to the target through a callback without the partial it needs throws. + /// On an evaluated context the callbacks on the path run once more and fill the result at once. + template + requires derivatives::derivative_value + const R& derivative(const property_base* target, const property_base* seed) { + using slot_type = derivatives::derivative_slot; + auto const what{[&] { return "derivative(): d '" + prop_key(target->traits().id) + "' / d '" + prop_key(seed->traits().id) + "'"; }}; + if (m_running) + throw std::logic_error(what() + ": called from a callback"); + if (auto const it{std::ranges::find_if(m_plans, [&](auto const& p) { return p.targets.front() == target && p.seed == seed; })}; + it != m_plans.end()) { + if (auto* slot{dynamic_cast(it->slot.get())}) + return slot->value; + throw std::logic_error(what() + " is requested as another type already"); + } + auto const target_rank{derivatives::property_rank(*target)}; + auto const seed_rank{derivatives::property_rank(*seed)}; + if (!target_rank || !seed_rank) + throw std::logic_error(what() + ": both must be scalars or symmetric tensors of rank 2 or 4"); + if (*target_rank + *seed_rank != derivatives::rank_of()) + throw std::logic_error(std::format("{} has rank {}, the requested type rank {}", what(), *target_rank + *seed_rank, + derivatives::rank_of())); + auto plan{plan_derivative(target, seed, *seed_rank, what())}; + auto slot{std::make_unique(*target_rank)}; + std::vector const zero(derivatives::width(*target_rank) * derivatives::width(*seed_rank)); + slot->unpack(zero.data()); + auto& value{slot->value}; + plan.slot = std::move(slot); + for (auto const& step : plan.steps) + step.material->request_partials(step.callback, true); + m_plans.push_back(std::move(plan)); + reactivate(); // a constant callback fills its partials once more + if (m_evaluated) + fill_now(m_plans.back()); + return value; + } + + /// Whether `seed` is upstream of `target` (or the target itself): a derivative of one by the other can be non-zero. + /// Without a path derivative() stays zero. + bool has_path(const property_base* target, const property_base* seed) const { + return upstream_of(target).contains(const_cast(seed)); + } + + /// Accumulates the requested derivatives from the partials the last update filled: D[out] += P * D[in], through a + /// residual system by the implicit function theorem. + void differentiate() { + if (m_plans.empty()) + return; + for (auto& [block, packed] : m_blocks) + block->pack(packed.data()); + for (auto& plan : m_plans) + finish(plan); + } + void dump() const { std::println(stderr, "=== Property execution order ({} properties) ===", m_property_execution_order.size()); @@ -577,6 +650,7 @@ class property_engine { } catch (...) { m_running = false; m_failed = true; + m_evaluated = false; throw; } m_running = false; @@ -586,6 +660,7 @@ class property_engine { void reactivate() { m_active = m_callbacks; m_retired.clear(); + distribute(); } void refresh_constants() { @@ -605,6 +680,7 @@ class property_engine { })}; m_retired.insert(m_retired.end(), done.begin(), done.end()); m_active.erase(done.begin(), done.end()); + distribute(); } /// Reruns a constant callback and throws when one of its material's constant outputs changes. @@ -671,6 +747,363 @@ class property_engine { } std::unordered_set m_routed; ///< outputs whose update reruns their writer + using value_type = typename Traits::value_type; + using matrix = Eigen::Matrix; + using column_vector = Eigen::Matrix; + using tensor2 = tmech::tensor; + using system_declaration = typename material_interface_type::system_declaration; + struct system_stage; + + /// Callbacks run in sequence (no system) or as a residual system's members until it converges. + struct stage { + system_stage* system{nullptr}; + std::vector callbacks; ///< the active ones, in topological order + }; + + static void run(std::vector const& callbacks) { + for (auto* prop : callbacks) + prop->traits().update(); + } + + /// Throws when tolerance or max_iter (max_passes) of a system is unusable (NaN included). + static void check_settings(system_stage const& s) { + auto const converges{s.declared->own_residual != nullptr}; + if (!(*s.declared->tolerance >= 0) || *s.declared->max_iter < (converges ? 1 : 0)) + throw std::invalid_argument(describe(s) + " needs a tolerance >= 0 and " + (converges ? "max_passes >= 1" : "max_iter >= 0")); + } + + /// Sets the unknowns, reruns the members and gathers the residuals until |R|_inf <= tolerance (a non-finite residual never + /// converges); max_iter counts the updates of the unknowns, converges() counts passes. + void solve(system_stage& s, std::vector const& callbacks) { + auto const& d{*s.declared}; + auto const limit{d.own_residual ? *d.max_iter - 1 : *d.max_iter}; + check_settings(s); + start(s); + for (int it{0};; ++it) { + if (d.report) + d.report->iterations = it; + scatter(s, s.x); + run(callbacks); + gather(s, s.R); + auto const norm{s.R.hasNaN() ? std::numeric_limits::quiet_NaN() : s.R.template lpNorm()}; + if (std::isfinite(norm) && norm <= *d.tolerance) + return; + if (it >= limit || !s.R.allFinite()) + throw cutback_error(d.own_residual ? std::format("'{}' did not reach {:g} in {} pass{} (last {:g})", s.name, *d.tolerance, it + 1, + it == 0 ? "" : "es", norm) + : std::format("system '{}' did not reach {:g} in {} iteration{} (last |R| {:g})", s.name, + *d.tolerance, it, it == 1 ? "" : "s", norm)); + step(s, callbacks, it); + } + } + + /// One update of the unknowns from the residual at x. + void step(system_stage& s, std::vector const& callbacks, int it) { + switch (s.declared->method) { + case system_method::picard: + if (s.width > 0) + s.x -= s.R; + return; + case system_method::newton: jacobian(s); break; + case system_method::newton_fd: jacobian_by_differences(s, callbacks); break; + } + s.lu.compute(s.J); + column_vector const dx{s.lu.solve(s.R)}; + if (!dx.allFinite()) + throw cutback_error(std::format("system '{}': the Jacobian is singular at iteration {}", s.name, it)); + s.x -= dx; + } + + /// J = dR/dx from the partials the members filled at x. + void jacobian(system_stage& s) { + auto& plan{*s.jacobian}; + for (auto const b : plan.blocks) + m_blocks[b].first->pack(m_blocks[b].second.data()); + accumulate(plan); + s.J.setZero(s.height, s.width); + for (std::size_t i{0}; i < s.residuals.size(); ++i) + if (auto const node{plan.target_nodes[i]}) + s.J.middleRows(s.rows[i], plan.D[*node].rows()) = plan.D[*node]; + } + + /// J by forward differences, one Mandel component of x at a time, with a step relative to the size of x (1 when + /// x is zero); the unknowns are back at x afterwards, also when a member throws. + void jacobian_by_differences(system_stage& s, std::vector const& callbacks) { + s.J.resize(s.height, s.width); + column_vector x{s.x}, R(s.height); + auto const norm{s.x.template lpNorm()}; + auto const scale{norm > 0 ? norm : value_type{1}}; + try { + for (Eigen::Index j{0}; j < s.width; ++j) { + auto const h{std::sqrt(std::numeric_limits::epsilon()) * std::max(std::abs(s.x[j]), scale)}; + x[j] = s.x[j] + h; + scatter(s, x); + run(callbacks); + gather(s, R); + s.J.col(j) = (R - s.R) / (x[j] - s.x[j]); + x[j] = s.x[j]; + } + } catch (...) { + scatter(s, s.x); + throw; + } + scatter(s, s.x); + } + + /// x from the seeds, a history unknown from its value at the end of the last step. + void start(system_stage& s) { + s.x.setZero(s.width); + for (std::size_t i{0}; i < s.unknowns.size(); ++i) + if (s.unknowns[i]->is_history()) + pack_value(*s.unknowns[i], s.ranks[i], s.x.data() + s.columns[i], true); + else if (s.ranks[i] == 0) + s.x[s.columns[i]] = s.declared->seeds[i]; + } + + void scatter(system_stage& s, column_vector const& x) { + for (std::size_t i{0}; i < s.unknowns.size(); ++i) + unpack_value(x.data() + s.columns[i], s.ranks[i], *s.unknowns[i]); + } + + void gather(system_stage& s, column_vector& R) { + for (std::size_t i{0}; i < s.residuals.size(); ++i) + pack_value(*s.residuals[i], s.residual_ranks[i], R.data() + s.rows[i], false); + } + + /// A scalar or a symmetric tensor into its Mandel components; `old`: a history property's old value. + static void pack_value(property_base& prop, std::size_t rank, value_type* out, bool old) { + auto const value{[&](std::type_identity) -> V const& { + if (old) + return dynamic_cast&>(prop).old_value(); + return derivatives::value_in(prop); + }}; + if (rank == 0) + *out = value(std::type_identity{}); + else + derivatives::pack(value(std::type_identity{}), out); + } + + static void unpack_value(value_type const* in, std::size_t rank, property_base& prop) { + if (rank == 0) + derivatives::value_in(prop) = *in; + else + derivatives::unpack(in, derivatives::value_in(prop)); + } + + /// A declaration with its properties: the unknowns and residuals in Mandel rows and columns. + std::unique_ptr make_system(const material_interface_type& mat, const system_declaration& d) { + auto s{std::make_unique()}; + s->declared = &d; + s->name = d.own_residual ? prop_key(d.own_residual->traits().id) : mat.name(); + auto const rank{[&](property_base* prop) { + auto const r{derivatives::property_rank(*prop)}; + if (!r || *r > 2) + throw std::logic_error(describe(*s) + ": '" + prop_key(prop->traits().id) + "' is neither a scalar nor a symmetric tensor"); + return *r; + }}; + for (auto const* unknown : d.unknowns) { + s->unknowns.push_back(find_in_graph(mat.name(), unknown->traits().id.name)); + s->ranks.push_back(rank(s->unknowns.back())); + s->columns.push_back(s->width); + s->width += static_cast(derivatives::width(s->ranks.back())); + } + if (d.own_residual) + s->residuals.push_back(find_in_graph(mat.name(), d.own_residual->traits().id.name)); + for (auto const* input : d.residuals) + s->residuals.push_back(find_in_graph(input->source_owner(), input->source_name())); + for (auto* residual : s->residuals) { + s->residual_ranks.push_back(rank(residual)); + s->rows.push_back(s->height); + s->height += static_cast(derivatives::width(s->residual_ranks.back())); + } + s->R.resize(s->height); + return s; + } + + static std::string describe(system_stage const& s) { + return s.declared->own_residual ? "converges(): '" + s.name + "'" : "system '" + s.name + "'"; + } + + /// Per system, by residual order. The members of a system with unknowns are the callbacks downstream of its + /// unknowns and upstream of its residuals. converges(): the callbacks on a path from the consumer Q of a closed + /// local edge P->Q (P downstream of Q, upstream of the residual) to the residual; its group must not share a + /// callback with an earlier one (nested loops). Each stage runs after the callbacks it reads, before the rest; + /// a system with unknowns shares no member. + void form_groups(const std::vector& interfaces) { + m_systems.clear(); + for (auto const* mat : interfaces) + for (auto const& d : mat->systems()) + m_systems.push_back(make_system(*mat, d)); + auto const position{[this](auto const& s) { + return std::ranges::min(s->residuals | std::views::transform([this](auto const* r) { + return std::ranges::find(m_property_execution_order, r) - m_property_execution_order.begin(); + })); + }}; + std::ranges::stable_sort(m_systems, {}, position); + + std::unordered_map> downstream_edges; + std::unordered_map> writers_of; + for (auto const& e : m_edges) { + if (e.kind != edge_kind::local) + downstream_edges[e.from].push_back(e.to); + if (e.kind == edge_kind::written_by) + writers_of[e.to].push_back(e.from); + } + // the callbacks that fill a property: itself, or its writers + auto const callbacks_filling{[&](const property_base* prop) { + return has_callback(prop) ? std::vector{prop} : writers_of[prop]; + }}; + std::unordered_map> reach_cache; + auto const reach{[&](const property_base* from) -> const std::unordered_set& { + if (auto const it{reach_cache.find(from)}; it != reach_cache.end()) + return it->second; + std::unordered_set seen{from}; + std::vector todo{from}; + while (!todo.empty()) { + auto const* prop{todo.back()}; + todo.pop_back(); + for (auto const* next : downstream_edges[prop]) + if (seen.insert(next).second) + todo.push_back(next); + } + return reach_cache.emplace(from, std::move(seen)).first->second; + }}; + auto const callbacks{m_property_execution_order | std::views::filter(has_callback) | std::ranges::to()}; + std::unordered_map member_of; + std::vector>> groups; + for (auto const& owned_system : m_systems) { + auto& s{*owned_system}; + auto const converges{s.declared->own_residual != nullptr}; + auto const loop{s.unknowns.empty()}; // converges(), or picard iterating a local loop + check_settings(s); + std::unordered_set upstream; + for (auto const* residual : s.residuals) + upstream.insert(upstream_of(residual).begin(), upstream_of(residual).end()); + std::unordered_set path; + std::unordered_set owners; // earlier groups whose loops this residual also sees + auto const converging_owner{[&](const property_base* prop) { + auto const it{member_of.find(prop)}; + return it != member_of.end() && it->second->unknowns.empty() ? it->second : nullptr; + }}; + if (!loop) { + std::vector todo{s.unknowns.begin(), s.unknowns.end()}; + while (!todo.empty()) { + auto const* prop{todo.back()}; + todo.pop_back(); + for (auto const* next : downstream_edges[prop]) + if (path.insert(next).second) + todo.push_back(next); + } + } else + for (auto const& e : m_edges) { + if (e.kind != edge_kind::local) + continue; + auto const producers{callbacks_filling(e.from)}; + for (auto const* q : callbacks_filling(e.to)) { + auto const closing{std::ranges::find_if(producers, [&](auto const* p) { return upstream.contains(p) && reach(q).contains(p); })}; + if (closing == producers.end()) + continue; + // the callbacks of the loop q ... p + auto const body{callbacks | std::views::filter([&](auto const* prop) { return reach(q).contains(prop) && reach(prop).contains(*closing); }) | + std::ranges::to>()}; + auto const owned{std::ranges::count_if(body, [&](auto const* prop) { return converging_owner(prop) != nullptr; })}; + if (owned == static_cast(body.size())) { // the loop of an earlier group + owners.insert(converging_owner(body.front())); + continue; + } + if (owned != 0) { + auto const* other{converging_owner(*std::ranges::find_if(body, [&](auto const* prop) { return converging_owner(prop) != nullptr; }))}; + throw std::logic_error((converges ? "converges(): the loop " : "system '" + s.name + "': the loop ") + prop_key((*closing)->traits().id) + + " -> " + prop_key(q->traits().id) + (converges ? " of '" + s.name + "'" : std::string{}) + + " runs through the iteration group of '" + other->name + "': nested loops are not supported"); + } + path.insert(reach(q).begin(), reach(q).end()); + } + } + std::vector members; + for (auto* prop : callbacks) { + if (!upstream.contains(prop) || !path.contains(prop)) + continue; + if (auto const it{member_of.find(prop)}; it != member_of.end()) { + if (!s.unknowns.empty() || !it->second->unknowns.empty()) + throw std::logic_error(describe(s) + ": '" + prop_key(prop->traits().id) + "' is a member of " + describe(*it->second) + + " too; systems are disjoint"); + continue; + } + if (m_stateful.contains(prop)) + throw std::logic_error(converges ? "converges(): '" + prop_key(prop->traits().id) + + "' is stateful() and inside the iteration group of '" + s.name + + "'; a pass must not count or record" + : describe(s) + ": '" + prop_key(prop->traits().id) + + "' is stateful() and between its unknowns and residuals; a pass must not count or record"); + members.push_back(prop); + } + if (members.empty()) + throw std::logic_error(!loop ? describe(s) + ": no callback lies between its unknowns and its residuals" + : describe(s) + " " + + (owners.empty() ? std::string{"has no closed local loop behind it: nothing to iterate (a local edge whose producer is " + "downstream of its consumer and upstream of the residual)"} + : "iterates no loop outside the group of '" + (*owners.begin())->name + + "': one residual per group")); + s.members.assign(members.begin(), members.end()); + for (auto* prop : members) + member_of.emplace(prop, &s); + for (auto const* unknown : s.unknowns) + member_of.emplace(unknown, &s); // what reads an unknown waits for its system + groups.emplace_back(&s, std::move(members)); + } + + // each system after the systems it reads; before it, the callbacks it reads that have not run yet + m_stages.clear(); + m_stage_of.clear(); + std::unordered_set placed; + auto const place{[&](auto&& props, system_stage* system) { + for (auto* prop : props) + m_stage_of[prop] = m_stages.size(); + m_stages.push_back({system, {}}); + if (system) + placed.insert(system); + }}; + while (!groups.empty()) { + auto const reads{[&](auto const& group) { + std::unordered_set all; + for (auto const* member : group.second) + all.insert(upstream_of(member).begin(), upstream_of(member).end()); + return all; + }}; + auto const next{std::ranges::find_if(groups, [&](auto const& group) { + return std::ranges::none_of(reads(group), [&](auto const* prop) { + auto const it{member_of.find(prop)}; + return it != member_of.end() && it->second != group.first && !placed.contains(it->second); + }); + })}; + if (next == groups.end()) + throw std::logic_error(groups.front().first->declared->own_residual && groups.back().first->declared->own_residual + ? "converges(): the iteration groups of '" + groups.front().first->name + "' and '" + + groups.back().first->name + "' read each other" + : describe(*groups.front().first) + " and " + describe(*groups.back().first) + " read each other"); + auto const before{reads(*next)}; + place(callbacks | std::views::filter([&](auto const* prop) { + return before.contains(prop) && !member_of.contains(prop) && !m_stage_of.contains(prop); + }), + nullptr); + place(next->second, next->first); + groups.erase(next); + } + place(callbacks | std::views::filter([&](auto const* prop) { return !m_stage_of.contains(prop); }), nullptr); + } + + /// The active callbacks into their stages. + void distribute() { + for (auto& s : m_stages) + s.callbacks.clear(); + for (auto* prop : m_active) + m_stages[m_stage_of.at(prop)].callbacks.push_back(prop); + } + + std::vector m_stages; ///< in execution order + std::unordered_map m_stage_of; ///< callback -> its stage + /// The upstream set of a property, computed once per finalize(). const std::unordered_set& upstream_of(const property_base* target) const { if (auto const it{m_upstream.find(target)}; it != m_upstream.end()) @@ -702,6 +1135,323 @@ class property_engine { return result; } + // --- derivatives (#217) --- + using read_declaration = typename material_interface_type::read_declaration; + + struct derivative_step { + std::size_t out, in; ///< nodes of the plan + std::size_t block; ///< in m_blocks + material_interface_type* material; + const property_base* callback; + }; + /// A residual system a derivative passes: its unknowns are seeds in `column`..., its residuals the targets from + /// `first_target` on; fold() replaces their columns by the implicit function theorem. + struct system_fold { + system_stage* system; + Eigen::Index column; + std::size_t first_target; + matrix J, B, dx_dp; ///< reused by every fold of the plan + Eigen::PartialPivLU lu; + }; + /// D = d node / d seeds over the nodes between the seeds and the targets that carry one, in topological order. + struct derivative_plan { + const property_base* seed{nullptr}; ///< of a requested derivative, whose target is targets.front() + std::vector targets; + std::vector> target_nodes; + std::vector> seeds; ///< (first column, width) of node k + Eigen::Index seed_width{0}; ///< the columns of `seed`, first + std::vector D; + std::vector steps; ///< in topological order of their outputs + std::vector blocks; ///< the partials the steps read, in m_blocks + std::vector folds; ///< in the order of the systems' stages + std::unique_ptr> slot; + }; + + /// A declared system as update() solves it, in Mandel components: unknown i in the columns from columns[i], + /// its residual in the rows from rows[i]. + struct system_stage { + const system_declaration* declared{nullptr}; + std::string name; ///< the system's, or the residual's for converges() + std::vector unknowns, residuals; + std::vector members; ///< the callbacks it reruns, from finalize() + std::vector ranks, residual_ranks; + std::vector columns, rows; + Eigen::Index width{0}, height{0}; + column_vector x, R; + matrix J; + Eigen::PartialPivLU lu; + std::optional jacobian; ///< newton: J from the partials, the unknowns as seeds + }; + std::vector> m_systems; + + /// D from the packed blocks, folded through the systems it passes, into the slot of a requested derivative. + void finish(derivative_plan& plan) { + accumulate(plan); + fold(plan); + if (auto const target{plan.target_nodes.front()}) + plan.slot->unpack(plan.D[*target].data()); + } + + /// Reruns the callbacks of a new plan at the values held and accumulates it. A stateful callback is not rerun: its + /// partials come with the next update(). + void fill_now(derivative_plan& plan) { + auto const on_plan{[&](auto const* prop) { + return std::ranges::any_of(plan.steps, [&](auto const& step) { return step.callback == prop; }); + }}; + if (std::ranges::any_of(m_callbacks, [&](auto const* prop) { return on_plan(prop) && m_stateful.contains(prop); })) + return; + guarded("derivative", [&] { + for (auto* prop : m_callbacks | std::views::filter(on_plan)) + prop->traits().update(); + }); + for (auto& [block, packed] : m_blocks) + block->pack(packed.data()); + finish(plan); + } + + /// Indexes the declared reads by output and runs every callback declaring partials with want_partials() as requested. + void index_partials(const std::vector& interfaces) { + m_plans.clear(); + m_blocks.clear(); + m_read_declarations.clear(); + for (auto* mat : interfaces) { + for (auto const& d : mat->declared_reads()) + for (auto const* out : d.outputs) + m_read_declarations.emplace(out, std::pair{mat, &d}); + for (auto const& p : mat->partials()) + if (auto* prop{find_in_graph(mat->name(), p.callback->traits().id.name)}; m_partials_wrapped.insert(prop).second) { + auto& update{prop->traits().update}; + update = [inner = std::move(update), mat, callback = p.callback] { mat->run_callback(callback, inner); }; + } + } + } + + /// Per system, its structural zero blocks (no path from the unknown to the residual) and, for newton, the plan + /// of J = dR/dx; the members fill their partials on every pass. + void plan_systems() { + for (auto& owned : m_systems) { + auto& s{*owned}; + if (auto* report{s.declared->report}) { + report->zero_blocks.clear(); + report->planned = true; + for (std::size_t i{0}; i < s.residuals.size(); ++i) + for (std::size_t j{0}; j < s.unknowns.size(); ++j) + if (!upstream_of(s.residuals[i]).contains(s.unknowns[j])) + report->zero_blocks.emplace_back(i, j); + } + if (s.declared->method != system_method::newton) + continue; + std::unordered_set up; + std::vector> seeds; + for (auto const* residual : s.residuals) + up.insert(upstream_of(residual).begin(), upstream_of(residual).end()); + for (std::size_t j{0}; j < s.unknowns.size(); ++j) + seeds.emplace_back(s.unknowns[j], static_cast(derivatives::width(s.ranks[j]))); + s.jacobian = plan_columns(std::vector(s.residuals.begin(), s.residuals.end()), seeds, up, describe(s) + ": its Jacobian"); + for (auto const& step : s.jacobian->steps) + step.material->request_partials(step.callback, true); + } + } + + /// d target / d seed: the systems whose unknowns the target reads (also through another system) add their + /// unknowns as seeds and their residuals as targets, so fold() can apply the implicit function theorem. + derivative_plan plan_derivative(const property_base* target, const property_base* seed, std::size_t seed_rank, std::string const& what) { + auto const seed_width{static_cast(derivatives::width(seed_rank))}; + auto up{upstream_of(target)}; + std::vector passed; + for (bool grown{true}; grown;) { + grown = false; + for (auto const& st : m_stages) { + auto* s{st.system}; + if (s && !std::ranges::contains(passed, s) && !std::ranges::contains(s->unknowns, seed) && + std::ranges::any_of(s->unknowns, [&](auto* u) { return up.contains(u); })) { + passed.push_back(s); + grown = true; + for (auto const* residual : s->residuals) + up.insert(upstream_of(residual).begin(), upstream_of(residual).end()); + } + } + } + if (!up.contains(const_cast(seed))) + passed.clear(); + auto const stage_index{[this](system_stage* s) { return std::ranges::find(m_stages, s, &stage::system) - m_stages.begin(); }}; + std::ranges::sort(passed, {}, stage_index); + for (auto const* s : passed) + if (std::ranges::contains(s->members, seed) || + std::ranges::any_of(m_edges, [&](auto const& e) { return e.kind == edge_kind::written_by && e.to == seed && std::ranges::contains(s->members, e.from); })) + throw std::logic_error(what + ": '" + prop_key(seed->traits().id) + "' is rerun by " + describe(*s) + + " and solved again; seed an input of the system or its unknown instead"); + + std::vector targets{target}; + std::vector> seeds{{seed, seed_width}}; + std::vector folds; + auto column{seed_width}; + for (auto* s : passed) { + folds.push_back({s, column, targets.size(), matrix::Zero(s->height, s->width), matrix::Zero(s->height, seed_width), {}, {}}); + for (std::size_t j{0}; j < s->unknowns.size(); ++j) + seeds.emplace_back(s->unknowns[j], static_cast(derivatives::width(s->ranks[j]))); + column += s->width; + targets.insert(targets.end(), s->residuals.begin(), s->residuals.end()); + } + auto plan{plan_columns(std::move(targets), seeds, up, what)}; + plan.seed = seed; + plan.seed_width = seed_width; + plan.folds = std::move(folds); + return plan; + } + + /// D = d node / d seeds over `up`: seed k is the identity in its own columns. Every node a seed reaches needs + /// the partial of each read that carries one, declared or declared zero (independent_of()). + derivative_plan plan_columns(std::vector targets, std::vector> const& seeds, + std::unordered_set const& up, std::string const& what) { + derivative_plan plan; + plan.targets = std::move(targets); + std::unordered_map node; + std::unordered_map via; ///< the read each node is first reached over + Eigen::Index columns{0}; + for (auto const& [seed, width] : seeds) { + node.emplace(seed, plan.seeds.size()); + plan.seeds.emplace_back(columns, width); + columns += width; + } + for (auto const& [first, width] : plan.seeds) + plan.D.emplace_back(width, columns); + for (auto const* prop : m_property_execution_order | std::views::filter([&](auto* p) { return up.contains(p) && !node.contains(p); })) { + std::vector reached; + if (auto const edges{m_reverse_edges.find(const_cast(prop))}; edges != m_reverse_edges.end()) + std::ranges::copy_if(edges->second, std::back_inserter(reached), [&](auto const* src) { return node.contains(src); }); + if (reached.empty()) + continue; + auto const declared{m_read_declarations.find(prop)}; + if (declared == m_read_declarations.end()) + throw missing_partial(what, plan.targets, reached.front(), prop, via, "its callback declares no reads"); + auto* mat{declared->second.first}; + auto const& d{*declared->second.second}; + std::vector steps; + std::vector sources; + Eigen::Index width{0}; + auto const take{[&](input_wire_base const* input, const property_base* own, const property_base* src) { + sources.push_back(src); + if (!node.contains(src)) + return; + auto const same_read{[&](auto const& p) { return p.input == input && p.own == own; }}; + auto const& partials{mat->partials()}; + if (auto const p{std::ranges::find_if(partials, [&](auto const& q) { return q.output == prop && same_read(q); })}; + p != partials.end()) { + steps.push_back({0, node.at(src), block_of(*p->block), mat, p->callback}); + width = static_cast(p->block->rows()); + via.try_emplace(prop, src); + } else if (!std::ranges::any_of(mat->independent_reads(), [&](auto const& i) { return i.callback == d.callback && same_read(i); })) { + throw missing_partial(what, plan.targets, src, prop, via, "its callback declares neither partials() nor independent_of() for that read"); + } + }}; + for (auto const* input : d.inputs) + if (input->edge_kind() == EdgeKind::Global) + if (auto const* src{find_in_graph(input->source_owner(), input->source_name())}) + take(input, nullptr, src); + for (auto const* own : d.own_reads) + take(nullptr, own, own); + for (auto const* src : reached) + if (!std::ranges::contains(sources, src) && !std::ranges::contains(d.outputs, src)) + throw missing_partial(what, plan.targets, src, prop, via, "it is reached over a '" + std::string{to_string(kind_of(src, prop))} + + "' edge, not a declared read, which carries no partials"); + if (steps.empty()) + continue; + auto const index{plan.D.size()}; + plan.D.emplace_back(width, columns); + for (auto const& step : steps) { + plan.steps.push_back({index, step.in, step.block, step.material, step.callback}); + if (!std::ranges::contains(plan.blocks, step.block)) + plan.blocks.push_back(step.block); + } + node.emplace(prop, index); + } + for (auto const* target : plan.targets) { + auto const it{node.find(target)}; + plan.target_nodes.push_back(it != node.end() ? std::optional{it->second} : std::nullopt); + } + return plan; + } + + /// D[out] += P * D[in] from the packed partials, the seeds the identity in their columns. + void accumulate(derivative_plan& plan) { + for (std::size_t k{0}; k < plan.D.size(); ++k) { + plan.D[k].setZero(); + if (k < plan.seeds.size()) + plan.D[k].middleCols(plan.seeds[k].first, plan.seeds[k].second).setIdentity(); + } + for (auto const& step : plan.steps) + plan.D[step.out].noalias() += m_blocks[step.block].second * plan.D[step.in]; + } + + /// Per system passed, at its converged point R(x, p) = 0: dx/dp = -J^-1 B with J = dR/dx and B = dR/dp taken + /// with the unknowns as seeds, then every node's D_p += D_x dx/dp. A singular J throws a cutback_error. + static void fold(derivative_plan& plan) { + auto const p{plan.seed_width}; + for (auto& f : plan.folds) { + auto const& s{*f.system}; + f.J.setZero(); + f.B.setZero(); + for (std::size_t i{0}; i < s.residuals.size(); ++i) + if (auto const node{plan.target_nodes[f.first_target + i]}) { + auto const& d{plan.D[*node]}; + f.J.middleRows(s.rows[i], d.rows()) = d.middleCols(f.column, s.width); + f.B.middleRows(s.rows[i], d.rows()) = d.leftCols(p); + } + f.lu.compute(f.J); + f.dx_dp = -f.lu.solve(f.B); + if (!f.dx_dp.allFinite()) + throw cutback_error(std::format("derivative(): system '{}' has a singular Jacobian at its converged point", s.name)); + for (auto& d : plan.D) + d.leftCols(p).noalias() += d.middleCols(f.column, s.width) * f.dx_dp; + } + } + + /// The packed block of a partial, shared by every plan through it. + std::size_t block_of(const derivatives::partial_block& block) { + if (auto const it{std::ranges::find(m_blocks, &block, &std::pair*, matrix>::first)}; + it != m_blocks.end()) + return static_cast(it - m_blocks.begin()); + m_blocks.emplace_back(&block, matrix(static_cast(block.rows()), static_cast(block.cols()))); + return m_blocks.size() - 1; + } + + edge_kind kind_of(const property_base* from, const property_base* to) const { + auto const it{std::ranges::find_if(m_edges, [&](auto const& e) { return e.from == from && e.to == to; })}; + return it != m_edges.end() ? it->kind : edge_kind::dependency; + } + + /// Names the path seed -> ... -> from -> at -> ... -> target in the error, toward the first target behind `at`. + std::logic_error missing_partial(std::string const& what, std::vector const& targets, const property_base* from, + const property_base* at, const std::unordered_map& via, + const std::string& why) const { + std::vector path{from, at}; + for (auto it{via.find(from)}; it != via.end(); it = via.find(it->second)) + path.insert(path.begin(), it->second); + auto const target{std::ranges::find_if(targets, [&](auto const* t) { return upstream_of(t).contains(const_cast(at)); })}; + if (target != targets.end()) { + auto const& up{upstream_of(*target)}; + for (auto const* p{at}; p != *target;) { // any edge into the target's upstream leads on to the target + auto const next{std::ranges::find_if(m_edges, [&](auto const& e) { + return e.from == p && e.kind != edge_kind::local && up.contains(const_cast(e.to)); + })}; + if (next == m_edges.end()) + break; + path.push_back(p = next->to); + } + } + std::string text; + for (auto const* p : path) + text += (text.empty() ? "" : " -> ") + prop_key(p->traits().id); + return std::logic_error(what + " passes '" + prop_key(at->traits().id) + "' from '" + prop_key(from->traits().id) + "', but " + why + + "; path " + text); + } + + std::vector m_plans; + std::vector*, matrix>> m_blocks; ///< packed once per differentiate() + std::unordered_map> m_read_declarations; + std::unordered_set m_partials_wrapped; ///< callbacks wrapped by index_partials() + std::vector m_execution_order; std::vector m_property_execution_order; std::unordered_map m_key_to_prop; // for string-based find @@ -723,6 +1473,7 @@ class property_engine { bool m_verify{NUMSIM_MATERIALS_DEBUG_CHECKS != 0}; bool m_fold{true}; bool m_running{false}; ///< an update is in progress + bool m_evaluated{false}; ///< an update() ran to its end since the last failure bool m_failed{false}; ///< the last update threw std::unique_ptr m_check_reads{std::make_unique(true)}; ///< stable address for the wrapped callbacks }; diff --git a/include/numsim-materials/core/unknown_spec.h b/include/numsim-materials/core/unknown_spec.h index 949101f..8f65129 100644 --- a/include/numsim-materials/core/unknown_spec.h +++ b/include/numsim-materials/core/unknown_spec.h @@ -29,6 +29,24 @@ struct unknown_spec { /// Used to declare structurally-zero blocks, which are then never wired. using block_ref = std::pair; +/// How a residual system (solvers/system.h) updates its unknowns. +enum class system_method { newton, newton_fd, picard }; + +/// An unknown of a residual system: it starts from `seed` (a tensor from zero) or, with `history`, from its +/// value at the end of the last step, kept per point. +struct system_unknown { + std::string name; + unknown_kind kind{unknown_kind::scalar}; + double seed{0}; + bool history{false}; +}; + +/// One equation of a residual system: the residual "material::property" and the unknown it is solved for. +struct system_equation { + std::string residual; + std::string unknown; +}; + } // namespace numsim::materials #endif // UNKNOWN_SPEC_H diff --git a/include/numsim-materials/default_materials.h b/include/numsim-materials/default_materials.h index fa0dce2..e5cc17b 100644 --- a/include/numsim-materials/default_materials.h +++ b/include/numsim-materials/default_materials.h @@ -6,6 +6,7 @@ #include "numsim-materials/solvers/backward_euler.h" #include "numsim-materials/solvers/local_newton.h" #include "numsim-materials/solvers/rk_integrator.h" +#include "numsim-materials/solvers/system.h" #include "numsim-materials/solvers/vector_newton.h" #include "numsim-materials/materials/scalar_stepper.h" #include "numsim-materials/materials/constant_scalar.h" @@ -194,6 +195,8 @@ void register_default_materials() { // A single registered type covers every unknown combination — the layout is // dispatched from the "unknowns" parameter, not baked into the template. factory.template register_type>("vector_newton"); + // unknowns and equations from any nodes, solved by the engine (#218) + factory.template register_type>("system"); // Unregistered until now, so the whole Runge-Kutta family -- every explicit // and implicit tableau -- was unreachable from a document. Nothing caught it // because the conformance scan looked only at materials/, not solvers/. diff --git a/include/numsim-materials/driver/material_point.h b/include/numsim-materials/driver/material_point.h index 4eb4b29..480c1f1 100644 --- a/include/numsim-materials/driver/material_point.h +++ b/include/numsim-materials/driver/material_point.h @@ -44,7 +44,8 @@ class material_point { optional>(para, "inputs", {}); // further "material::property" history kept out of load_history() / store_history() optional>(para, "host_owned", {}); - // "material::property": evaluate() updates again until its magnitude is at most the tolerance + // deprecated, removed later: "material::property" evaluate() updates again until its magnitude is at most the + // tolerance; a material declaring converges() is converged after the one update optional(para, "residual_source", std::string{}); optional(para, "residual_tolerance", value_type{1e-10}); optional(para, "max_passes", 100); @@ -131,6 +132,28 @@ class material_point { return output(src.material, src.property); } + /// d output / d input, filled by evaluate(): output "material::property", input an input of this point (its history + /// output: the strain, the temperature) or "material::property". T is a scalar for two scalars, else a tensor of + /// the summed rank (stress by strain: rank 4). + template + requires derivatives::derivative_value + T const &derivative(std::string const &output, std::string const &input) const { + auto const out{connection_source::parse(output)}; + find(out.material, out.property); + auto const seed{seed_of(input)}; + find(seed.material, seed.property); + return m_ctx.template derivative(out.material, out.property, seed.material, seed.property); + } + + /// Whether the input (as in derivative()) is upstream of the output: otherwise derivative() is zero. + bool has_path(std::string const &output, std::string const &input) const { + auto const out{connection_source::parse(output)}; + find(out.material, out.property); + auto const seed{seed_of(input)}; + find(seed.material, seed.property); + return m_ctx.has_path(out.material, out.property, seed.material, seed.property); + } + /// A history output, for both ends of the step. template history_property const &history(std::string const &material, std::string const &property) const { @@ -142,7 +165,8 @@ class material_point { // --- Evaluation --- - /// One graph update, repeated to the fixed point when a residual is configured. + /// One graph update; with the deprecated residual_source, repeated until that residual is converged (a graph + /// whose fixed point is declared with converges() is converged after the first). void evaluate() { if (!m_props_readers.empty() && !m_props_bound) throw fatal_error("material_point: this model reads its material constants per call — call bind_props() " @@ -198,6 +222,22 @@ class material_point { return parameters; } + /// An input's history output, or "material::property". + connection_source seed_of(std::string const &input) const { + if (input.contains("::")) + return connection_source::parse(input); + if (!std::ranges::contains(m_inputs, input)) + throw umat::fatal_error("material_point: '" + input + "' is not an input of this material point"); + std::vector history; + for (auto const *prop : m_ctx.property_execution_order()) + if (prop->is_history() && prop->traits().id.owner == input) + history.push_back(prop->traits().id.name); + if (history.size() != 1) + throw umat::fatal_error("material_point: input '" + input + "' has " + std::to_string(history.size()) + + " history outputs; name the one to differentiate by as 'material::property'"); + return {input, history.front()}; + } + property_base *find(std::string const &material, std::string const &property) const { if (auto *prop = m_ctx.find_property(material, property)) return prop; diff --git a/include/numsim-materials/io/json_parameter_converter.h b/include/numsim-materials/io/json_parameter_converter.h index 146569a..0e03276 100644 --- a/include/numsim-materials/io/json_parameter_converter.h +++ b/include/numsim-materials/io/json_parameter_converter.h @@ -201,6 +201,53 @@ json_reader_registry make_default_json_registry() { return result; }); + // a list of objects whose keys are among `allowed`: the message names the parameter and the key + const auto require_list = [](const JsonType& j, const std::string& key, std::initializer_list allowed) { + if (!adapter::is_array(j)) + throw std::runtime_error(key + " is a list of objects"); + for (const auto& elem : j) { + if (!adapter::is_object(elem)) + throw std::runtime_error(key + " is a list of objects"); + adapter::for_each_key(elem, [&](const std::string& name) { + if (std::ranges::none_of(allowed, [&](const char* a) { return name == a; })) + throw std::runtime_error(key + ": unknown key '" + name + "'"); + }); + } + }; + // the unknowns of a residual system (solvers/system.h): "seed" is a number or "history" + reg.template add>([require_list](const JsonType& j, const std::string& key) -> std::any { + std::vector result; + require_list(j, key, {"name", "kind", "seed"}); + for (const auto& elem : j) { + auto& u = result.emplace_back(); + u.name = adapter::template get(adapter::at(elem, "name")); + const auto kind = adapter::template get(adapter::at(elem, "kind")); + if (kind == "sym_tensor") + u.kind = unknown_kind::sym_tensor; + else if (kind != "scalar") + throw std::runtime_error(key + ": unknown '" + u.name + "' has the kind '" + kind + "' (scalar or sym_tensor)"); + if (!adapter::contains(elem, "seed")) + continue; + const JsonType seed = adapter::at(elem, "seed"); + if (adapter::is_number(seed)) + u.seed = adapter::template get(seed); + else if (adapter::is_string(seed) && adapter::template get(seed) == "history") + u.history = true; + else + throw std::runtime_error(key + ": the seed of '" + u.name + "' is a number or \"history\""); + } + return result; + }); + // its equations: {"residual": "material::property", "for": unknown}, without "for" when it has no unknowns + reg.template add>([require_list](const JsonType& j, const std::string& key) -> std::any { + std::vector result; + require_list(j, key, {"residual", "for"}); + for (const auto& elem : j) + result.push_back({adapter::template get(adapter::at(elem, "residual")), + adapter::contains(elem, "for") ? adapter::template get(adapter::at(elem, "for")) : std::string{}}); + return result; + }); + // a vector of string pairs, shared by "zero_blocks" and "terms": the message names the key it was given // and says what an entry must be in that parameter's words const auto string_pairs = [](const JsonType& j, const std::string& key, const char* entry) -> std::any { diff --git a/include/numsim-materials/io/orientation_phases.h b/include/numsim-materials/io/orientation_phases.h index 4ed65d7..3aec05f 100644 --- a/include/numsim-materials/io/orientation_phases.h +++ b/include/numsim-materials/io/orientation_phases.h @@ -29,7 +29,8 @@ namespace numsim::materials { /// "fraction_source", "orientation" [a11, a22, a33, a12, a13, a23], "per_direction", "stress" ("{q}" /// template); optional "precision" (11), "phase", "fraction_property", "matrix_tangent" and "tangent" /// (default the phase nodes' tangents), "localisation" (linear | nonlinear: the fibres localise with their -/// current tangents, read by the "phase" template over a local edge) and "matrix_used_tangent" (nonlinear). +/// current tangents, read by the "phase" template over a local edge), "matrix_used_tangent", "tolerance" and +/// "max_passes" (nonlinear: the composite's fixed point). /// Sources are "material::property". template std::vector expand_orientation_phases(JsonType const &e) { @@ -64,8 +65,9 @@ std::vector expand_orientation_phases(JsonType const &e) { if (nonlinear && phase_template.value("tangent_edge", std::string{}) != "local") throw std::invalid_argument("orientation_phases: a nonlinear localisation reads the fibre tangents over a local edge: the \"phase\" " "template needs \"tangent_edge\": \"local\" and an \"initial_tangent_source\""); - if (!nonlinear && e.contains("matrix_used_tangent")) - throw std::invalid_argument("orientation_phases: matrix_used_tangent goes with \"localisation\": \"nonlinear\""); + for (char const *key : {"matrix_used_tangent", "tolerance", "max_passes"}) + if (!nonlinear && e.contains(key)) + throw std::invalid_argument(std::string{"orientation_phases: "} + key + " goes with \"localisation\": \"nonlinear\""); JsonType const &deck{require("per_direction")}; if (!deck.is_array() || deck.empty()) throw std::invalid_argument("orientation_phases: per_direction lists the entries of one direction"); @@ -125,6 +127,9 @@ std::vector expand_orientation_phases(JsonType const &e) { out.front()["property"] = e.at("fraction_property"); out.push_back({{"type", "mean_field_strain"}, {"name", strain_name}, {"strain_source", require("strain_source")}, {"phases", phases}, {"kernel", kernel}}); out.push_back({{"type", nonlinear ? "mean_field_composite_local_nonlinear" : "mean_field_composite_local_linear"}, {"name", name}, {"phases", locals}}); + for (char const *key : {"tolerance", "max_passes"}) + if (e.contains(key)) + out.back()[key] = e.at(key); return out; } diff --git a/include/numsim-materials/materials/linear_elasticity.h b/include/numsim-materials/materials/linear_elasticity.h index 8111a6f..82acd96 100755 --- a/include/numsim-materials/materials/linear_elasticity.h +++ b/include/numsim-materials/materials/linear_elasticity.h @@ -42,7 +42,9 @@ class linear_elasticity { update_tangent(); // readable before the first update base::add_callback(&linear_elasticity::update_tangent, m_C).reads(); - base::add_callback(&linear_elasticity::update_stress, m_sig).reads(m_C, m_eps, m_eigenstrain); + base::add_callback(&linear_elasticity::update_stress, m_sig) + .reads(m_C, m_eps, m_eigenstrain) + .partials(m_sig, m_C, m_eps, m_eigenstrain); } /** @@ -69,6 +71,19 @@ class linear_elasticity m_sig = tmech::dcontract(m_C, tensor2{m_eps.get() - m_eigenstrain->get()}); else m_sig = tmech::dcontract(m_C, m_eps.get()); + if (base::want_partials()) + update_partials(); + } + + /// d sigma / d eps = C, d sigma / d eps_eig = -C, d sigma / d C = Isym (x) (eps - eps_eig). + void update_partials() { + base::partial(m_sig, m_eps) = m_C; + if (m_eigenstrain) { + base::partial(m_sig, *m_eigenstrain) = value_type{-1} * m_C; + base::partial(m_sig, m_C) = tmech::otimes(sym_identity(), tensor2{m_eps.get() - m_eigenstrain->get()}); + } else { + base::partial(m_sig, m_C) = tmech::otimes(sym_identity(), m_eps.get()); + } } void update_tangent(){ diff --git a/include/numsim-materials/materials/mean_field_composite_local_nonlinear.h b/include/numsim-materials/materials/mean_field_composite_local_nonlinear.h index 5d1bc0f..9b80f2f 100644 --- a/include/numsim-materials/materials/mean_field_composite_local_nonlinear.h +++ b/include/numsim-materials/materials/mean_field_composite_local_nonlinear.h @@ -2,6 +2,7 @@ #define NUMSIM_MATERIALS_MEAN_FIELD_COMPOSITE_LOCAL_NONLINEAR_H #include +#include #include #include #include @@ -14,10 +15,11 @@ namespace numsim::materials { /// A composite of nonlinear phases at their local strains: the localisation used the tangents of the /// previous pass ("used_tangent", the phase node's on a local edge), so "residual" = max_r ||C_r - C_used,r|| -/// / ||C_r|| says when the passes agree. Outputs "stress", "tangent", "residual". +/// / ||C_r|| says when the passes agree. update() repeats the passes until it is at most "tolerance" +/// (converges()). Outputs "stress", "tangent", "residual". /// /// Parameters: "phases" [{"concentration", "stress", "tangent", "fraction", "used_tangent" (the phases whose -/// tangent changes)}: "material::property"]. +/// tangent changes)}: "material::property"], "tolerance" (1e-10), "max_passes" (200: then a cutback_error). template class mean_field_composite_local_nonlinear final : public detail::mean_field_composite_local_base, Traits> { @@ -36,6 +38,14 @@ class mean_field_composite_local_nonlinear final if (m_used.empty()) throw std::invalid_argument(base::where() + "no phase names the 'used_tangent' its strain was localised with"); base::register_outputs(m_residual).reads(base::m_A, base::m_sigma, base::m_C, base::m_c, m_used | std::views::values); + base::converges(m_residual, base::template get_parameter("tolerance"), base::template get_parameter("max_passes")); + } + + static typename base::input_parameter_controller parameters() { + auto para{base::parameters()}; + optional(para, "tolerance", value_type{1e-10}); + optional(para, "max_passes", 200); + return para; } static constexpr char const *type_name() { return "mean_field_composite_local_nonlinear"; } @@ -49,7 +59,9 @@ class mean_field_composite_local_nonlinear final m_residual = 0; for (auto const &[r, used] : m_used) { auto const &C{base::m_C[r]->get()}; - m_residual = std::max(m_residual, tmech::norm(tensor4{C - used->get()}) / tmech::norm(C)); + auto const change{tmech::norm(tensor4{C - used->get()}) / tmech::norm(C)}; + if (std::isnan(change) || change > m_residual) // std::max drops a NaN + m_residual = change; } } diff --git a/include/numsim-materials/materials/thermal_eigenstrain.h b/include/numsim-materials/materials/thermal_eigenstrain.h index edae4f6..060c1f5 100644 --- a/include/numsim-materials/materials/thermal_eigenstrain.h +++ b/include/numsim-materials/materials/thermal_eigenstrain.h @@ -30,7 +30,7 @@ class thermal_eigenstrain final : public material_base("reference_temperature")), m_T(base::template add_input_history( connection_source{base::template get_parameter("temperature_source"), "state"}, EdgeKind::Global)) { - base::add_callback(&thermal_eigenstrain::compute, m_eigenstrain).reads(m_alpha, m_T); + base::add_callback(&thermal_eigenstrain::compute, m_eigenstrain).reads(m_alpha, m_T).partials(m_eigenstrain, m_alpha, m_T); } static input_parameter_controller parameters() { @@ -41,7 +41,13 @@ class thermal_eigenstrain final : public material_base(); } + void compute() { + m_eigenstrain = m_alpha.get() * (m_T.new_value() - m_T_ref) * tmech::eye(); + if (base::want_partials()) { + base::partial(m_eigenstrain, m_alpha) = (m_T.new_value() - m_T_ref) * tmech::eye(); + base::partial(m_eigenstrain, m_T) = m_alpha.get() * tmech::eye(); + } + } private: tensor2 &m_eigenstrain; diff --git a/include/numsim-materials/report/report_json.h b/include/numsim-materials/report/report_json.h index 15cc526..2f45263 100644 --- a/include/numsim-materials/report/report_json.h +++ b/include/numsim-materials/report/report_json.h @@ -186,7 +186,8 @@ inline std::vector output_stems(nlohmann::json const &doc) { /// "point": the parameters of the material-point driver (strain_source or deformation_gradient_source, /// stress_source, optional tangent_source, stress_property, tangent_property ("damaged_stress") and time_source, -/// and optionally tolerance, max_iter, fd_step, and the fixed point residual_source, residual_tolerance, max_passes). +/// and optionally tolerance, max_iter, fd_step, and the deprecated fixed point residual_source, residual_tolerance, +/// max_passes: a composite declaring converges() takes its tolerance and max_passes as parameters). template typename Traits::ParameterHandler parse_point(nlohmann::json const &doc) { using T = typename Traits::value_type; diff --git a/include/numsim-materials/solvers/system.h b/include/numsim-materials/solvers/system.h new file mode 100644 index 0000000..bfdc31b --- /dev/null +++ b/include/numsim-materials/solvers/system.h @@ -0,0 +1,144 @@ +#ifndef NUMSIM_MATERIALS_SOLVERS_SYSTEM_H +#define NUMSIM_MATERIALS_SOLVERS_SYSTEM_H + +#include +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/material_base.h" +#include "numsim-materials/core/parameter_helpers.h" +#include "numsim-materials/core/unknown_spec.h" + +namespace numsim::materials { + +/// A residual system set up at run time (#218), "type": "system": its unknowns are outputs of this node, its +/// residuals ("equations": [{"residual": "material::property", "for": unknown}]) outputs of any nodes reading them. +/// update() solves it by "method" newton (J from the partials), newton_fd or picard (x -= R; without unknowns +/// the group of converges()); "seed": "history" starts an unknown from its value at the end of the last step. +template +class residual_system 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; + using tensor2 = tmech::tensor; + + template + explicit residual_system(Args&&... args) + : base(std::forward(args)...), m_tolerance(base::template get_parameter("tolerance")), + m_max_iter(base::template get_parameter("max_iter")) { + auto const& unknowns{base::template get_parameter>("unknowns")}; + auto const& equations{base::template get_parameter>("equations")}; + auto const method{check(unknowns, equations)}; + std::vector outputs; + std::vector seeds; + std::vector residuals; + for (auto const& u : unknowns) { + auto const& residual{std::ranges::find(equations, u.name, &system_equation::unknown)->residual}; + if (u.kind == unknown_kind::scalar) + outputs.push_back(add(u, residual, residuals)); + else + outputs.push_back(add(u, residual, residuals)); + seeds.push_back(static_cast(u.seed)); + m_names.push_back(u.name); + } + if (unknowns.empty()) + for (auto const& e : equations) + residuals.push_back(&input(e.residual)); + base::declare_system(method, outputs, std::move(seeds), std::move(residuals), m_tolerance, m_max_iter, m_report); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + optional(para, "method", std::string{"newton"}); + optional>(para, "unknowns", {}); + required>(para, "equations"); + optional(para, "tolerance", value_type{1e-10}); + optional(para, "max_iter", 50); + return para; + } + + /// The updates of x the last solve took. + int iterations() const noexcept { return m_report.iterations; } + + /// The Jacobian blocks (row unknown, column unknown) without a path from the column unknown to the residual + /// of the row: structurally zero (after finalize()). + std::vector zero_blocks() const { + if (!m_report.planned) + throw std::logic_error("system '" + base::name() + "': zero_blocks() is known after finalize()"); + return m_report.zero_blocks | + std::views::transform([this](auto const& b) { return block_ref{m_names[b.first], m_names[b.second]}; }) | + std::ranges::to(); + } + +private: + /// The unknown `u` (history or not) and the input of its residual; returns the output. + template + void const* add(system_unknown const& u, std::string const& residual, std::vector& residuals) { + residuals.push_back(&input(residual)); + if (u.history) + return &base::template add_history_output(u.name); + return &base::template add_output(u.name); + } + + template + input_property const& input(std::string const& residual) { + auto const split{residual.find("::")}; + return base::template add_input(residual.substr(0, split), residual.substr(split + 2), EdgeKind::Local); + } + + [[noreturn]] void fail(std::string const& what) const { + throw std::invalid_argument("system '" + base::name() + "': " + what); + } + + /// Validates the parameters before any property is created; returns the method. + system_method check(std::vector const& unknowns, std::vector const& equations) const { + auto const& name{base::template get_parameter("method")}; + auto const method{name == "newton" ? system_method::newton : name == "newton_fd" ? system_method::newton_fd + : name == "picard" ? system_method::picard + : throw std::invalid_argument("system '" + base::name() + "': unknown method '" + + name + "' (newton, newton_fd or picard)")}; + if (unknowns.empty() && method != system_method::picard) + fail("'unknowns' is empty; only picard iterates without unknowns"); + if (equations.empty()) + fail("'equations' is empty"); + if (!(m_tolerance >= 0) || m_max_iter < 0) + fail("needs a tolerance >= 0 and max_iter >= 0"); + for (auto it{equations.begin()}; it != equations.end(); ++it) { + auto const& e{*it}; + if (std::ranges::contains(equations.begin(), it, e.residual, &system_equation::residual)) + fail("residual '" + e.residual + "' is used by more than one equation: one residual per unknown"); + auto const split{e.residual.find("::")}; + if (split == std::string::npos || split == 0 || split + 2 >= e.residual.size()) + fail("residual '" + e.residual + "' is not \"material::property\""); + if (unknowns.empty() && !e.unknown.empty()) + fail("the equation for '" + e.unknown + "': a system without unknowns has no \"for\""); + if (!unknowns.empty() && !std::ranges::contains(unknowns, e.unknown, &system_unknown::name)) + fail("residual '" + e.residual + "' is for '" + e.unknown + "', which is not an unknown"); + } + for (auto it{unknowns.begin()}; it != unknowns.end(); ++it) { + if (it->name.empty()) + fail("an unknown has an empty name"); + if (std::ranges::contains(unknowns.begin(), it, it->name, &system_unknown::name)) + fail("unknown '" + it->name + "' is declared twice"); + if (it->kind == unknown_kind::sym_tensor && it->seed != 0) + fail("unknown '" + it->name + "': the seed of a sym_tensor is 0 or \"history\""); + auto const count{std::ranges::count(equations, it->name, &system_equation::unknown)}; + if (count != 1) + fail("unknown '" + it->name + "' has " + std::to_string(count) + " equations, needs one"); + } + return method; + } + + value_type const& m_tolerance; + int const& m_max_iter; + std::vector m_names; + system_report m_report; +}; + +} // namespace numsim::materials + +#endif // NUMSIM_MATERIALS_SOLVERS_SYSTEM_H diff --git a/include/numsim-materials/umat/material_point_evaluator.h b/include/numsim-materials/umat/material_point_evaluator.h index b8d0383..2be3c09 100644 --- a/include/numsim-materials/umat/material_point_evaluator.h +++ b/include/numsim-materials/umat/material_point_evaluator.h @@ -71,9 +71,9 @@ class material_point_evaluator { std::optional tangent_source{}; /// external_deformation_gradient_source (finite strain). APPENDED, see tangent_source. std::string deformation_gradient_source{}; - /// Optional "material::property" of a residual: the graph is updated again until it is at most - /// residual_tolerance (a fixed point closed over local edges, e.g. nonlinear mean-field phases). - /// APPENDED, see tangent_source. + /// Deprecated, removed later: optional "material::property" of a residual, the graph is updated again until + /// it is at most residual_tolerance. A material declaring converges() (mean_field_composite_local_nonlinear) + /// closes its fixed point in the update itself. APPENDED, see tangent_source. std::string residual_property{}; value_type residual_tolerance{1e-10}; int max_passes{100}; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1844164..149783a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -81,6 +81,7 @@ add_numsim_test(test_orientation test_orientation.cpp) add_numsim_test(test_short_fibre_composite test_short_fibre_composite.cpp) add_numsim_test(test_incremental_elasticity test_incremental_elasticity.cpp) add_numsim_test(test_thermal_eigenstrain test_thermal_eigenstrain.cpp) +add_numsim_test(test_partial_derivatives test_partial_derivatives.cpp) add_numsim_test(test_kamal_sourour_rate test_kamal_sourour_rate.cpp) add_numsim_test(test_curing_eigenstrain test_curing_eigenstrain.cpp) add_numsim_test(test_damage_laws test_damage_laws.cpp) @@ -88,6 +89,9 @@ add_numsim_test(test_uvwmaterial_cards test_uvwmaterial_cards.cpp) add_numsim_test(test_voigt_composite test_voigt_composite.cpp) add_numsim_test(test_hencky_elasticity test_hencky_elasticity.cpp) add_numsim_test(test_material_point test_material_point.cpp) +add_numsim_test(test_iteration_groups test_iteration_groups.cpp) +add_numsim_test(test_residual_system test_residual_system.cpp) +add_numsim_test(test_system_plasticity test_system_plasticity.cpp) # Data dumper for plotting (not a test — standalone executable) add_executable(plot_data plot_data.cpp) diff --git a/tests/test_iteration_groups.cpp b/tests/test_iteration_groups.cpp new file mode 100644 index 0000000..c0a9a4c --- /dev/null +++ b/tests/test_iteration_groups.cpp @@ -0,0 +1,652 @@ +// Iteration groups: a fixed point closed over local edges runs inside update() (#214). +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/errors.h" +#include "numsim-materials/core/material_base.h" +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/core/parameter_helpers.h" +#include "numsim-materials/driver/material_point.h" +#include "numsim-materials/materials/constant_scalar.h" +#include "numsim-materials/postprocessing/property_recorder.h" +#include "numsim-materials/solvers/system.h" + +namespace { + +namespace nm = numsim::materials; +using policy = nm::material_policy_default; +using T = policy::value_type; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; + +/// The names of the callbacks in the order they ran. +std::vector g_log; + +/// "target" = goal::value; "calls" counts its runs (stateful: never folded, never in a group). +template +class probe final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit probe(Args &&...args) + : base(std::forward(args)...), m_target(base::template add_output("target")), + m_calls(base::template add_output("calls")), m_goal(base::template add_input("goal", "value")) { + base::add_callback(&probe::compute, m_target, m_calls).reads(m_goal).stateful(); + } + + static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; } + + void compute() { + g_log.push_back(base::name()); + m_target = m_goal.get(); + m_calls += 1; + } + +private: + T &m_target; + T &m_calls; + nm::input_property const &m_goal; +}; + +/// "used" = "::x" of the previous pass, over a local edge (as a phase reads its tangent). +template +class guess final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit guess(Args &&...args) + : base(std::forward(args)...), m_used(base::template add_output("used")), + m_x(base::template add_input(base::template get_parameter("source"), + base::template get_parameter("property"), nm::EdgeKind::Local)) { + base::add_callback(&guess::copy, m_used).reads(m_x); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::required(para, "source"); + nm::optional(para, "property", std::string{"x"}); + return para; + } + + void copy() { + g_log.push_back(base::name()); + m_used = m_x.get(); + } + +private: + T &m_used; + nm::input_property const &m_x; +}; + +/// x = (used + target) / 2 and residual = |x - used|: every pass halves the distance to the target. It +/// declares converges(residual), and with "twice" a second residual "check" over the same callbacks. +template +class relaxer final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit relaxer(Args &&...args) + : base(std::forward(args)...), m_x(base::template add_output("x")), + m_residual(base::template add_output("residual")), m_check(base::template add_output("check")), + m_used(base::template add_input(base::template get_parameter("guess"), "used")), + m_target(target(base::template get_parameter("target_source"))) { + auto handle{base::add_callback(&relaxer::compute, m_x, m_residual, m_check).reads(m_used, m_target)}; + if (base::template get_parameter("stateful")) + handle.stateful(); + auto const &tolerance{base::template get_parameter("tolerance")}; + auto const &max_passes{base::template get_parameter("max_passes")}; + if (base::template get_parameter("converges")) + base::converges(m_residual, tolerance, max_passes); + if (base::template get_parameter("twice")) + base::converges(m_check, tolerance, max_passes); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "guess", std::string{"guess"}); + nm::optional(para, "target_source", std::string{"probe::target"}); + nm::optional(para, "tolerance", T{1e-8}); + nm::optional(para, "max_passes", 100); + nm::optional(para, "stateful", false); + nm::optional(para, "twice", false); + nm::optional(para, "converges", true); + return para; + } + + void compute() { + g_log.push_back(base::name()); + m_x = (m_used.get() + m_target.get()) / 2; + m_residual = std::abs(m_x - m_used.get()); + m_check = m_residual; + } + +private: + T &m_x; + T &m_residual; + T &m_check; + nm::input_property const &m_used; + nm::input_property const &m_target; + + nm::input_property const &target(std::string const &source) { + auto const src{nm::connection_source::parse(source)}; + return base::template add_input(src.material, src.property); + } +}; + +/// "residual" = 1, declared converging with no local edge behind it. +template +class lonely final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit lonely(Args &&...args) + : base(std::forward(args)...), m_residual(base::template add_output("residual")) { + base::add_callback(&lonely::compute, m_residual).reads(); + base::converges(m_residual, base::template get_parameter("tolerance"), base::template get_parameter("max_passes")); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "tolerance", T{1e-8}); + nm::optional(para, "max_passes", 10); + return para; + } + + void compute() { m_residual = 1; } + +private: + T &m_residual; +}; + +/// "x" = "" (its output), a callback that only passes a value on. +template +class follow final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit follow(Args &&...args) + : base(std::forward(args)...), m_x(base::template add_output("x")), m_in(source(base::template get_parameter("source"))) { + base::add_callback(&follow::copy, m_x).reads(m_in); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::required(para, "source"); + return para; + } + + void copy() { m_x = m_in.get(); } + +private: + T &m_x; + nm::input_property const &m_in; + + nm::input_property const &source(std::string const &s) { + auto const src{nm::connection_source::parse(s)}; + return base::template add_input(src.material, src.property); + } +}; + +/// Two halves of a solver (as j2_plasticity and its hardening law): each reads the other over a local edge, and +/// each counts its runs (stateful: never folded, and an error inside a group). +template +class half final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit half(Args &&...args) + : base(std::forward(args)...), m_out(base::template add_output("out")), m_calls(base::template add_output("calls")), + m_in(base::template add_input(base::template get_parameter("other"), "out", nm::EdgeKind::Local)) { + base::add_callback(&half::compute, m_out, m_calls).reads(m_in).stateful(); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::required(para, "other"); + return para; + } + + void compute() { + m_out = m_in.get() + 1; + m_calls += 1; + } + + // a deleted overload keeps the pointers the declaration stores from dangling + template + static constexpr bool accepts_references = requires(Self &s, T &r, T const &tol, int const &n) { s.converges(r, tol, n); }; + template + static constexpr bool accepts_literals = requires(Self &s, T &r) { s.converges(r, 1e-10, 100); }; + template + static constexpr bool accepts_literal_tolerance = requires(Self &s, T &r, int const &n) { s.converges(r, 1e-10, n); }; + template + static constexpr bool accepts_literal_passes = requires(Self &s, T &r, T const &tol) { s.converges(r, tol, 100); }; + template + static constexpr bool declares_literals = requires(Self &s, T &r) { s.declare_converges(r, 1e-10, 100); }; + template + static constexpr bool declares_references = requires(Self &s, T &r, T const &tol, int const &n) { s.declare_converges(r, tol, n); }; + +private: + T &m_out; + T &m_calls; + nm::input_property const &m_in; +}; + +static_assert(half::accepts_references<> && half::declares_references<>); +static_assert(!half::accepts_literals<> && !half::accepts_literal_tolerance<> && !half::accepts_literal_passes<>, + "converges() keeps pointers to its arguments: a temporary would dangle"); +static_assert(!half::declares_literals<>); + +void constant(ctx_type &ctx, std::string const &name, T value) { + param_type p; + p.insert("name", name); + p.insert("value", value); + ctx.create>(p); +} + +/// goal -> probe -> [guess -> relax] (the group, closed by guess reading relax::x over a local edge), and with +/// `recorded` a property_recorder of relax::x after it. `relax` adds parameters of the relaxer. +std::unique_ptr build(param_type relax = {}, bool recorded = false, bool finalize = true) { + auto owner{std::make_unique()}; + auto &ctx{*owner}; + constant(ctx, "goal", 2); + param_type p; + p.insert("name", "probe"); + ctx.create>(p); + p.clear(); + p.insert("name", "guess"); + p.insert("source", "relax"); + ctx.create>(p); + relax.insert("name", "relax"); + ctx.create>(relax); + if (recorded) { + p.clear(); + p.insert("name", "rec"); + p.insert>("scalar_sources", {"relax::x"}); + ctx.create>(p); + } + if (finalize) + ctx.finalize(); + return owner; +} + +/// The message of the exception of type E that `f` throws; empty when it throws none. +template +std::string thrown(F &&f) { + try { + f(); + } catch (E const &e) { + return e.what(); + } + return {}; +} + +TEST(IterationGroups, UpdateIteratesToTheResidual) { + auto ctx{build()}; + ctx->update(); + EXPECT_LE(ctx->get("relax", "residual"), 1e-8); + EXPECT_NEAR(ctx->get("relax", "x"), 2, 1e-7); +} + +TEST(IterationGroups, ACallbackBeforeTheGroupRunsOncePerUpdate) { + auto ctx{build()}; + for (int n{1}; n <= 3; ++n) { + ctx->update(); + EXPECT_EQ(ctx->get("probe", "calls"), n) << "the group took about 27 passes, the probe ran once"; + } +} + +TEST(IterationGroups, ARecorderAfterTheGroupWritesOneRowPerUpdate) { + auto ctx{build({}, true)}; + for (int n{0}; n < 3; ++n) + ctx->update(); + auto const &rec{dynamic_cast &>(*ctx->find("rec"))}; + ASSERT_EQ(rec.buffer().rows(), 3u); + EXPECT_NEAR(rec.buffer().at(0, 0), 2, 1e-7) << "the row holds the converged value"; +} + +// residual_source is deprecated: the material point finds the residual converged after the one update +TEST(IterationGroups, TheMaterialPointResidualSourceStillWorks) { + auto ctx{build({}, true)}; + param_type p; + p.insert("residual_source", "relax::residual"); + p.insert("residual_tolerance", 1e-8); + nm::material_point point{*ctx, p}; + for (int n{0}; n < 3; ++n) + point.evaluate(); + EXPECT_EQ(dynamic_cast &>(*ctx->find("rec")).buffer().rows(), 3u); + EXPECT_EQ(ctx->get("probe", "calls"), 3); +} + +TEST(IterationGroups, TooFewPassesThrowACutbackNamingTheResidual) { + param_type relax; + relax.insert("max_passes", 3); + auto ctx{build(relax)}; + auto const message{thrown([&] { ctx->update(); })}; + EXPECT_NE(message.find("relax::residual"), std::string::npos) << message; + EXPECT_NE(message.find("3 passes"), std::string::npos) << message; + EXPECT_THROW(ctx->update(), std::logic_error) << "the context is failed until revert()"; + ctx->revert(); + dynamic_cast &>(*ctx->find("relax")).set_parameter("max_passes", 100); + EXPECT_NO_THROW(ctx->update()) << "the parameter is read on every update"; +} + +TEST(IterationGroups, ANaNResidualNeverConverges) { + auto ctx{build()}; + dynamic_cast &>(*ctx->find("goal")).set_parameter("value", std::numeric_limits::quiet_NaN()); + auto const message{thrown([&] { ctx->update(); })}; + EXPECT_NE(message.find("in 1 pass "), std::string::npos) << "stops at the first pass: " << message; +} + +TEST(IterationGroups, UpdatePropertyOnAMemberRunsItsUpstreamOnce) { + auto ctx{build()}; + ctx->update_property("relax", "residual"); + EXPECT_EQ(ctx->get("relax", "x"), 1) << "one pass from x = 0: no iteration"; + EXPECT_EQ(ctx->get("relax", "residual"), 1); +} + +TEST(IterationGroups, ConstantsUpstreamStayFoldedAndFollowSetParameter) { + auto ctx{build()}; + ctx->update(); + EXPECT_TRUE(ctx->is_constant("goal", "value")); + dynamic_cast &>(*ctx->find("goal")).set_parameter("value", 4); + ctx->update(); + EXPECT_NEAR(ctx->get("relax", "x"), 4, 1e-7); + EXPECT_LE(ctx->get("relax", "residual"), 1e-8); +} + +TEST(IterationGroups, AGroupDownstreamOfAnotherRunsAfterIt) { + auto ctx{build({}, false, false)}; + param_type p; + p.insert("name", "guess2"); + p.insert("source", "relax2"); + ctx->create>(p); + p.clear(); + p.insert("name", "relax2"); + p.insert("guess", "guess2"); + p.insert("target_source", "relax::x"); // the converged x of the first group + ctx->create>(p); + ctx->finalize(); + ctx->update(); + EXPECT_LE(ctx->get("relax", "residual"), 1e-8); + EXPECT_LE(ctx->get("relax2", "residual"), 1e-8); + EXPECT_NEAR(ctx->get("relax2", "x"), 2, 1e-7); +} + +TEST(IterationGroups, AGroupDownstreamOfAnotherRunsAfterItInTheCallLog) { + auto ctx{build({}, false, false)}; + param_type p; + p.insert("name", "guess2"); + p.insert("source", "relax2"); + ctx->create>(p); + p.clear(); + p.insert("name", "relax2"); + p.insert("guess", "guess2"); + p.insert("target_source", "relax::x"); + ctx->create>(p); + ctx->finalize(); + g_log.clear(); + ctx->update(); + auto const first{[&](std::string const &name) { return std::ranges::find(g_log, name) - g_log.begin(); }}; + auto const last{[&](std::string const &name) { return std::ranges::find_last(g_log, name).begin() - g_log.begin(); }}; + EXPECT_EQ(first("probe"), 0); + EXPECT_EQ(std::ranges::count(g_log, "probe"), 1); + EXPECT_LT(last("relax"), first("guess2")) << "the second group starts after the first has converged"; + EXPECT_LT(last("guess"), first("guess2")); + EXPECT_GT(std::ranges::count(g_log, "relax2"), 1) << "and iterates"; + EXPECT_EQ(g_log.back(), "relax2"); +} + +TEST(IterationGroups, GroupsThatReadEachOtherThrowAtFinalize) { + param_type relax; + relax.insert("target_source", "guess2::used"); // each group reads the front of the other + auto ctx{build(relax, false, false)}; + param_type p; + p.insert("name", "guess2"); + p.insert("source", "relax2"); + ctx->create>(p); + p.clear(); + p.insert("name", "relax2"); + p.insert("guess", "guess2"); + p.insert("target_source", "guess::used"); + ctx->create>(p); + auto const message{thrown([&] { ctx->finalize(); })}; + EXPECT_NE(message.find("relax::residual"), std::string::npos) << message; + EXPECT_NE(message.find("relax2::residual"), std::string::npos) << message; + EXPECT_NE(message.find("read each other"), std::string::npos) << message; +} + +TEST(IterationGroups, ALocalEdgeWhoseProducerIsOutsideTheGroupThrowsAtFinalize) { + ctx_type ctx; // guess reads probe::target of the previous pass, but probe is not downstream of guess + constant(ctx, "goal", 2); + param_type p; + p.insert("name", "probe"); + ctx.create>(p); + p.clear(); + p.insert("name", "guess"); + p.insert("source", "probe"); + p.insert("property", "target"); + ctx.create>(p); + p.clear(); + p.insert("name", "relax"); + ctx.create>(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_NE(message.find("relax::residual"), std::string::npos) << message; + EXPECT_NE(message.find("closed local loop"), std::string::npos) << message; +} + +TEST(IterationGroups, AResidualInTheMiddleOfTheLoopThrowsAtFinalize) { + ctx_type ctx; // guess -> relax -> follow -> (local) guess: the producer of the edge is behind the residual + constant(ctx, "goal", 2); + param_type p; + p.insert("name", "probe"); + ctx.create>(p); + p.clear(); + p.insert("name", "guess"); + p.insert("source", "follow"); + ctx.create>(p); + p.clear(); + p.insert("name", "relax"); + ctx.create>(p); + p.clear(); + p.insert("name", "follow"); + p.insert("source", "relax::x"); + ctx.create>(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_NE(message.find("relax::residual"), std::string::npos) << message; + EXPECT_NE(message.find("closed local loop"), std::string::npos) << message; +} + +TEST(IterationGroups, ANestedLoopThrowsAtFinalize) { + ctx_type ctx; // the inner loop (gi, ri) feeds the outer one (go, ro), which feeds the inner one back through ro::x + constant(ctx, "goal", 2); + param_type p; + p.insert("name", "probe"); + ctx.create>(p); + p.clear(); + p.insert("name", "go"); + p.insert("source", "ro"); + ctx.create>(p); + p.clear(); + p.insert("name", "gi"); + p.insert("source", "ri"); + ctx.create>(p); + p.clear(); + p.insert("name", "ri"); + p.insert("guess", "gi"); + p.insert("target_source", "go::used"); + ctx.create>(p); + p.clear(); + p.insert("name", "ro"); + p.insert("guess", "go"); + p.insert("target_source", "ri::x"); + ctx.create>(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_NE(message.find("ri::residual"), std::string::npos) << message; + EXPECT_NE(message.find("ro::residual"), std::string::npos) << message; + EXPECT_NE(message.find("nested"), std::string::npos) << message; +} + +TEST(IterationGroups, ALoopThroughTheGroupOfAnotherResidualThrowsAtFinalize) { + param_type relax; // relax2 -> (local) guess3 -> relax (target_source) -> relax2: the loop runs through relax's group + relax.insert("target_source", "guess3::used"); + auto ctx{build(relax, false, false)}; + param_type p; + p.insert("name", "guess3"); + p.insert("source", "relax2"); + ctx->create>(p); + p.clear(); + p.insert("name", "relax2"); + p.insert("guess", "guess3"); + p.insert("target_source", "relax::x"); + ctx->create>(p); + auto const message{thrown([&] { ctx->finalize(); })}; + EXPECT_NE(message.find("relax2::residual"), std::string::npos) << message; + EXPECT_NE(message.find("relax::residual"), std::string::npos) << message; + EXPECT_NE(message.find("nested"), std::string::npos) << message; +} + +TEST(IterationGroups, ASolverWithLocalEdgesUpstreamOfAResidualRunsOncePerUpdate) { + ctx_type ctx; // a and b read each other over local edges, as j2_plasticity and its hardening law do + constant(ctx, "goal", 2); + param_type p; + p.insert("name", "a"); + p.insert("other", "b"); + ctx.create>(p); + p.clear(); + p.insert("name", "b"); + p.insert("other", "a"); + ctx.create>(p); + p.clear(); + p.insert("name", "probe"); + ctx.create>(p); + p.clear(); + p.insert("name", "guess"); + p.insert("source", "relax"); + ctx.create>(p); + p.clear(); + p.insert("name", "relax"); + p.insert("target_source", "b::out"); + ctx.create>(p); + ASSERT_NO_THROW(ctx.finalize()) << "neither is a member of the group: both are stateful"; + for (int n{1}; n <= 3; ++n) { + ctx.update(); + EXPECT_EQ(ctx.get("a", "calls"), n); + EXPECT_EQ(ctx.get("b", "calls"), n); + EXPECT_LE(ctx.get("relax", "residual"), 1e-8); + } +} + +TEST(IterationGroups, AnInfiniteResidualIsNotConvergedByAnInfiniteTolerance) { + param_type relax; + relax.insert("tolerance", std::numeric_limits::infinity()); + auto ctx{build(relax)}; + dynamic_cast &>(*ctx->find("goal")).set_parameter("value", std::numeric_limits::infinity()); + auto const message{thrown([&] { ctx->update(); })}; + EXPECT_NE(message.find("in 1 pass "), std::string::npos) << "stops at the first pass: " << message; +} + +TEST(IterationGroups, SettingsChangedAfterFinalizeAreCheckedAtEveryUpdate) { + auto ctx{build()}; + auto &relax{dynamic_cast &>(*ctx->find("relax"))}; + relax.set_parameter("tolerance", -1); + EXPECT_NE(thrown([&] { ctx->update(); }).find("relax::residual"), std::string::npos); + ctx->revert(); + relax.set_parameter("tolerance", 1e-8); + relax.set_parameter("max_passes", 0); + EXPECT_NE(thrown([&] { ctx->update(); }).find("relax::residual"), std::string::npos); + ctx->revert(); + relax.set_parameter("max_passes", 100); + EXPECT_NO_THROW(ctx->update()); +} + +TEST(IterationGroups, ConvergesWithoutALocalEdgeBehindItThrowsAtFinalize) { + ctx_type ctx; + param_type p; + p.insert("name", "alone"); + ctx.create>(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_NE(message.find("alone::residual"), std::string::npos) << message; + EXPECT_NE(message.find("closed local loop"), std::string::npos) << message; +} + +TEST(IterationGroups, AStatefulCallbackInsideAGroupThrowsAtFinalize) { + param_type relax; + relax.insert("stateful", true); + auto ctx{build(relax, false, false)}; + auto const message{thrown([&] { ctx->finalize(); })}; + EXPECT_NE(message.find("relax::x"), std::string::npos) << message; + EXPECT_NE(message.find("stateful"), std::string::npos) << message; +} + +TEST(IterationGroups, OverlappingGroupsThrowAtFinalize) { + param_type relax; + relax.insert("twice", true); + auto ctx{build(relax, false, false)}; + auto const message{thrown([&] { ctx->finalize(); })}; + EXPECT_NE(message.find("relax::residual"), std::string::npos) << message; + EXPECT_NE(message.find("relax::check"), std::string::npos) << message; +} + +TEST(IterationGroups, InvalidSettingsThrowAtFinalize) { + for (auto const &[key, value] : {std::pair{"tolerance", -1.0}, std::pair{"tolerance", std::numeric_limits::quiet_NaN()}}) { + param_type relax; + relax.insert(key, value); + auto ctx{build(relax, false, false)}; + EXPECT_NE(thrown([&] { ctx->finalize(); }).find("relax::residual"), std::string::npos) << value; + } + param_type relax; + relax.insert("max_passes", 0); + auto ctx{build(relax, false, false)}; + EXPECT_NE(thrown([&] { ctx->finalize(); }).find("relax::residual"), std::string::npos); +} + +// converges() is shorthand for a picard system without unknowns (#218) +TEST(IterationGroups, APicardSystemWithoutUnknownsIsTheGroup) { + auto const with_system{[](int max_iter) { + param_type relax; + relax.insert("converges", false); + auto ctx{build(relax, false, false)}; + param_type p; + p.insert("name", "fix"); + p.insert("method", "picard"); + p.insert>("equations", {{"relax::residual", ""}}); + p.insert("tolerance", 1e-8); + p.insert("max_iter", max_iter); + ctx->create>(p); + ctx->finalize(); + return ctx; + }}; + auto ctx{with_system(99)}; + auto reference{build()}; + for (int n{1}; n <= 3; ++n) { + ctx->update(); + reference->update(); + EXPECT_EQ(ctx->get("relax", "x"), reference->get("relax", "x")); + EXPECT_EQ(ctx->get("probe", "calls"), n); + if (n == 1) { + EXPECT_GT(dynamic_cast &>(*ctx->find("fix")).iterations(), 10) << "the passes of the first update"; + } + } + auto short_of_it{with_system(2)}; + EXPECT_NE(thrown([&] { short_of_it->update(); }).find("'fix'"), std::string::npos); +} + +} // namespace diff --git a/tests/test_material_conformance.cpp b/tests/test_material_conformance.cpp index 5d94152..fe2b539 100644 --- a/tests/test_material_conformance.cpp +++ b/tests/test_material_conformance.cpp @@ -5,8 +5,10 @@ #include #include #include +#include "numsim-materials/core/partial_check.h" #include "numsim-materials/default_materials.h" #include "numsim-materials/io/json_material_factory.h" +#include "numsim-materials/umat/external_state_source.h" /// Generic conformance checks over every material, so a new one is covered /// without a bespoke test file. @@ -153,4 +155,62 @@ TEST(MaterialConformance, MissingRequiredParametersThrow) { } } +/// Materials whose callbacks declare partials or independence (.partials(), .independent_of()), checked by finite differences below. +const std::set kPartialsChecked = {"linear_elasticity", "thermal_eigenstrain"}; + +/// Header stems under materials/ and solvers/ that declare partials or independent reads. +std::set headers_declaring_partials() { + namespace fs = std::filesystem; + std::set out; + for (const char* directory : {NUMSIM_MATERIALS_INCLUDE_DIR, NUMSIM_SOLVERS_INCLUDE_DIR}) + for (const auto& e : fs::directory_iterator(fs::path{directory})) { + if (e.path().extension() != ".h") continue; + std::ifstream in(e.path()); + const std::string src{std::istreambuf_iterator(in), std::istreambuf_iterator()}; + if (src.find(".partials(") != std::string::npos || src.find(".independent_of(") != std::string::npos) + out.insert(e.path().stem().string()); + } + return out; +} + +TEST(MaterialConformance, EveryMaterialDeclaringPartialsIsChecked) { + EXPECT_EQ(headers_declaring_partials(), kPartialsChecked) + << "a material declaring partials or independent reads joins kPartialsChecked and the model of DeclaredPartialsMatchFiniteDifferences"; +} + +/// Every declared partial against central differences of its own callback, at a state off every symmetry. +TEST(MaterialConformance, DeclaredPartialsMatchFiniteDifferences) { + nm::material_context ctx; + for (const auto* source : {"strain", "temperature"}) { + policy::ParameterHandler p; + p.insert("name", source); + if (std::string{source} == "strain") + ctx.create>(p); + else + ctx.create>(p); + } + for (const auto* document : { + R"({"type": "constant_scalar", "name": "alpha", "value": 1.2e-5})", + R"({"type": "thermal_eigenstrain", "name": "thermal_eigenstrain", "alpha_source": "alpha::value", "reference_temperature": 20.0})", + R"({"type": "linear_elasticity", "name": "linear_elasticity", "strain_producer_name": "strain", "eigenstrain_source": "thermal_eigenstrain", "K": 150.0, "G": 70.0})"}) { + const nlohmann::json j = nlohmann::json::parse(document); + nm::create_from_json(ctx, j); + } + ctx.finalize(); + tmech::tensor strain; + strain(0, 0) = 1e-3, strain(1, 1) = -4e-4, strain(2, 2) = 2e-4; + strain(0, 1) = strain(1, 0) = 3e-4, strain(0, 2) = strain(2, 0) = -1e-4, strain(1, 2) = strain(2, 1) = 5e-4; + dynamic_cast*>(ctx.find("strain"))->bind({}, strain); + dynamic_cast*>(ctx.find("temperature"))->bind(20, 75); + ctx.update(); + std::set checked; + for (const auto* material : ctx.materials()) + if (!material->partials().empty()) { + checked.insert(material->name()); + for (const auto& mismatch : nm::check_partials(ctx, material->name())) + ADD_FAILURE() << mismatch; + } + EXPECT_EQ(checked, kPartialsChecked); +} + } // namespace diff --git a/tests/test_mean_field_damage.cpp b/tests/test_mean_field_damage.cpp index a52f519..0ce4534 100644 --- a/tests/test_mean_field_damage.cpp +++ b/tests/test_mean_field_damage.cpp @@ -1,9 +1,12 @@ #include #include +#include #include #include +#include #include #include "numsim-materials/core/elastic_constants.h" +#include "numsim-materials/core/errors.h" #include "numsim-materials/core/material_context.h" #include "numsim-materials/materials/constant_scalar.h" #include "numsim-materials/materials/eshelby_sphere.h" @@ -25,6 +28,7 @@ #include "numsim-materials/driver/material_point_driver.h" #include "numsim-materials/driver/external_state_source.h" #include "numsim-materials/umat/json_model.h" +#include "numsim-materials/umat/material_point_evaluator.h" namespace { @@ -137,13 +141,19 @@ std::unique_ptr build(bool nonlinear, std::string const &kernel = "mor ctx.create>(p); p.clear(); p.insert("name", "mt"); - p.insert("phases", {{{"concentration", "mt_strain::concentration_matrix"}, - {"stress", nonlinear ? "damage::damaged_stress" : "elastic::stress"}, - {"tangent", nonlinear ? "damage::damaged_tangent" : "elastic::tangent"}, - {"fraction", "matrix::volume_fraction"}, {"used_tangent", "matrix::tangent"}}, - {{"concentration", "mt_strain::concentration_spheres"}, {"stress", "stiff::stress"}, - {"tangent", "C_i::tangent"}, {"fraction", "spheres::volume_fraction"}}}); - ctx.create>(p); + nm::phase_sources phases{{{"concentration", "mt_strain::concentration_matrix"}, + {"stress", nonlinear ? "damage::damaged_stress" : "elastic::stress"}, + {"tangent", nonlinear ? "damage::damaged_tangent" : "elastic::tangent"}, + {"fraction", "matrix::volume_fraction"}}, + {{"concentration", "mt_strain::concentration_spheres"}, {"stress", "stiff::stress"}, + {"tangent", "C_i::tangent"}, {"fraction", "spheres::volume_fraction"}}}; + if (nonlinear) + phases[0]["used_tangent"] = "matrix::tangent"; + p.insert("phases", phases); + if (nonlinear) + ctx.create>(p); + else // a converging composite needs the local edge + ctx.create>(p); ctx.finalize(); return owner; } @@ -154,23 +164,13 @@ tensor2 uniaxial(T e) { return eps; } -/// One step: passes until the phase tangents the strains were formed with agree with the materials'. -int converge(ctx_type &ctx, T tol = 1e-10) { - int passes{1}; - ctx.update(); - while (ctx.get("mt", "residual") > tol && passes < 200) { - ctx.update(); - ++passes; - } - return passes; -} - -/// Loads the composite in equal steps to the strain `eps`, committing every step. +/// Loads the composite in equal steps to the strain `eps`, committing every step: each update() repeats the +/// passes until the phase tangents the strains were formed with agree with the materials' (converges()). void load(ctx_type &ctx, tensor2 const &eps, int steps) { auto &strain{*dynamic_cast *>(ctx.find("strain"))}; for (int n{1}; n <= steps; ++n) { strain.bind(tensor2{eps * T(n - 1) / T(steps)}, tensor2{eps * T(n) / T(steps)}); - converge(ctx); + ctx.update(); ctx.commit(); } } @@ -223,6 +223,7 @@ TEST(MeanFieldDamage, TangentIsTheConcentrationWeightedPhaseTangent) { TEST(MeanFieldDamage, ResidualDecreasesOverThePassesOfAStep) { auto ctx{build(true)}; load(*ctx, uniaxial(0.004), 4); + ctx->find("mt")->set_parameter("tolerance", std::numeric_limits::infinity()); // one pass per update auto &strain{*dynamic_cast *>(ctx->find("strain"))}; strain.bind(uniaxial(0.004), uniaxial(0.006)); // a large damaging step: several passes ctx->update(); @@ -263,6 +264,53 @@ TEST(MeanFieldDamage, TangentMatchesFiniteDifferencesOnASmallStep) { EXPECT_LT(d.tangent_error(res), 1e-2); } +// the fixed point inside update() takes bitwise the passes of the material point's loop it replaces: there, the +// composite converges after every pass and residual_source repeats the whole update +TEST(MeanFieldDamage, TheEngineIteratesBitwiseAsTheMaterialPointDid) { + for (auto const &control : {nm::load_control::strain(), nm::load_control::uniaxial_stress()}) { + auto engine{build(true)}, host{build(true)}; + host->find("mt")->set_parameter("tolerance", std::numeric_limits::infinity()); + auto in_engine{mt_point()}; + in_engine.insert("residual_source", std::string{}); + nm::material_point_driver e(*engine, in_engine), h(*host, mt_point()); + for (int n{1}; n <= 24; ++n) { + T const eps{(n <= 16 ? n : 32 - n) * 5e-4}; // loading, then unloading + auto const a{e.step(uniaxial(eps), tensor2{}, control, n - 1, 1)}, b{h.step(uniaxial(eps), tensor2{}, control, n - 1, 1)}; + ASSERT_TRUE(a.converged) << n; + ASSERT_TRUE(b.converged) << n; + EXPECT_EQ(a.residuals, b.residuals) << n; + EXPECT_EQ(a.stress, b.stress) << n; + EXPECT_EQ(engine->get("mt", "tangent"), host->get("mt", "tangent")) << n; + EXPECT_EQ(engine->get("law", "damage"), host->get("law", "damage")) << n; + EXPECT_EQ(engine->get("mt", "residual"), host->get("mt", "residual")) << n; + e.accept(); + h.accept(); + } + EXPECT_GT(engine->get("law", "damage"), 0); + } +} + +// the UMAT closes the fixed point without residual_property: the composite converges in update() +TEST(MeanFieldDamage, TheUmatEvaluatorConvergesWithoutAResidualProperty) { + auto ctx{build(true)}, reference{build(true)}; + nm::umat::material_point_evaluator::config cfg; + cfg.strain_source = "strain"; + cfg.stress_source = "mt"; + nm::umat::material_point_evaluator eval(*ctx, cfg); + std::vector statev(eval.nstatv(), 0); + T stran[6]{}, stress[6]{}, ddsdde[36]{}; + T const dstran[6]{1e-3, 0, 0, 0, 0, 0}; + for (int n{1}; n <= 6; ++n) { + eval.evaluate({.stran = stran, .dstran = dstran, .stress = stress, .ddsdde = ddsdde, .statev = statev}); + EXPECT_LE(ctx->get("mt", "residual"), 1e-10) << n; + stran[0] += dstran[0]; + } + load(*reference, uniaxial(0.006), 6); + EXPECT_GT(reference->get("law", "damage"), 0); + T const expected{reference->get("mt", "stress")(0, 0)}; + EXPECT_NEAR(stress[0], expected, 1e-10 * std::abs(expected)); +} + TEST(MeanFieldDamage, StressControlReEvaluatesFromTheCommittedLocalStrains) { auto ctx{build(true)}; nm::material_point_driver d(*ctx, mt_point()); @@ -367,6 +415,24 @@ TEST(MeanFieldDamage, ReussEqualisesThePhaseStressesAsTheStepVanishes) { EXPECT_LT(fine, 1e-5); } +TEST(MeanFieldDamage, ANaNTangentIsNotConverged) { + auto ctx{build(true)}; + dynamic_cast &>(*ctx->find("elastic")).set_parameter("K", std::numeric_limits::quiet_NaN()); + auto &strain{*dynamic_cast *>(ctx->find("strain"))}; + strain.bind(tensor2{}, uniaxial(1e-3)); + EXPECT_THROW(ctx->update(), nm::cutback_error) << "the residual of a NaN tangent is NaN, not 0"; +} + +TEST(MeanFieldDamage, AZeroTangentIsNotConverged) { + auto ctx{build(true)}; + auto &elastic{dynamic_cast &>(*ctx->find("elastic"))}; + elastic.set_parameter("K", 0); + elastic.set_parameter("G", 0); + auto &strain{*dynamic_cast *>(ctx->find("strain"))}; + strain.bind(tensor2{}, uniaxial(1e-3)); + EXPECT_THROW(ctx->update(), nm::cutback_error) << "norm(C) = 0 makes the residual 0 / 0"; +} + TEST(MeanFieldDamage, RejectsAnInconsistentWiring) { { ctx_type ctx; // a local tangent edge without an initial tangent diff --git a/tests/test_orientation_phases.cpp b/tests/test_orientation_phases.cpp index 825d5aa..df54075 100644 --- a/tests/test_orientation_phases.cpp +++ b/tests/test_orientation_phases.cpp @@ -321,6 +321,12 @@ TEST(NonlinearFibres, ExpandIntoTheNonlinearComposite) { EXPECT_EQ(named(expand(e), "sf_3").at("tangent_source"), "fibre_3"); e["matrix_used_tangent"] = "matrix::tangent"; EXPECT_EQ(named(expand(e), "sf").at("phases")[0].at("used_tangent"), "matrix::tangent"); + EXPECT_FALSE(named(expand(e), "sf").contains("tolerance")) << "the composite's defaults"; + e["tolerance"] = 1e-12; + e["max_passes"] = 50; + json const converging = named(expand(e), "sf"); + EXPECT_EQ(converging.at("tolerance"), 1e-12); + EXPECT_EQ(converging.at("max_passes"), 50); } TEST(NonlinearFibres, NeedALocalTangentEdge) { @@ -338,6 +344,9 @@ TEST(NonlinearFibres, NeedALocalTangentEdge) { e["localisation"] = "linear"; e["matrix_used_tangent"] = "matrix::tangent"; // only a nonlinear composite reads it EXPECT_THROW((void)expand(e), std::invalid_argument); + e.erase("matrix_used_tangent"); + e["tolerance"] = 1e-12; // only a nonlinear composite iterates + EXPECT_THROW((void)expand(e), std::invalid_argument); } /// The driver on "strain" -> "sf", closing the local fixed point on "sf::residual". diff --git a/tests/test_partial_derivatives.cpp b/tests/test_partial_derivatives.cpp new file mode 100644 index 0000000..ba05c0a --- /dev/null +++ b/tests/test_partial_derivatives.cpp @@ -0,0 +1,780 @@ +#include +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/material_base.h" +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/core/parameter_helpers.h" +#include "numsim-materials/core/partial_check.h" +#include "numsim-materials/core/tensor_identities.h" +#include "numsim-materials/driver/material_point.h" +#include "numsim-materials/materials/constant_scalar.h" +#include "numsim-materials/materials/linear_elasticity.h" +#include "numsim-materials/materials/thermal_eigenstrain.h" +#include "numsim-materials/umat/external_state_source.h" + +// Partial derivatives per callback and their forward accumulation (#217). + +namespace { + +namespace nm = numsim::materials; +using policy = nm::material_policy_default; +using T = policy::value_type; +using tensor2 = tmech::tensor; +using tensor4 = tmech::tensor; +using tensor6 = tmech::tensor; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; + +tensor2 const I{tmech::eye()}; +tensor4 const Is{nm::sym_identity()}; + +tensor2 sym(T a, T b, T c, T d, T e, T f) { + tensor2 x; + x(0, 0) = a, x(1, 1) = b, x(2, 2) = c; + x(1, 2) = x(2, 1) = d, x(0, 2) = x(2, 0) = e, x(0, 1) = x(1, 0) = f; + return x; +} +tensor2 const B2{sym(1, -2, 0.5, 0.3, -0.7, 1.1)}; +tensor2 const C2{sym(0.4, 0.2, -1, 0.6, 0.1, -0.3)}; +tensor2 const E2{sym(-0.5, 1.5, 0.25, -0.2, 0.9, 0.35)}; +tensor4 const A4{tensor4{2 * Is + 0.5 * tmech::otimes(I, I) + tmech::otimes(B2, B2)}}; +tensor4 const B4{tensor4{3 * Is - tmech::otimes(I, I) + 0.5 * tmech::otimes(C2, C2)}}; + +// One law per edge of the contraction table: value and its partial in tmech. +struct norm_trace { // scalar of tensor2 + using in = tensor2; + using out = T; + static T value(tensor2 const &x) { return tmech::dcontract(x, x) + tmech::trace(x); } + static tensor2 partial(tensor2 const &x) { return 2 * x + I; } +}; +struct square { // scalar of scalar + using in = T; + using out = T; + static T value(T s) { return s * s; } + static T partial(T s) { return 2 * s; } +}; +struct scale_b { // tensor2 of scalar + using in = T; + using out = tensor2; + static tensor2 value(T s) { return s * B2; } + static tensor2 partial(T) { return B2; } +}; +struct map_a { // tensor2 of tensor2 + using in = tensor2; + using out = tensor2; + static tensor2 value(tensor2 const &x) { return tmech::dcontract(A4, x); } + static tensor4 partial(tensor2 const &) { return A4; } +}; +struct map_b { + using in = tensor2; + using out = tensor2; + static tensor2 value(tensor2 const &x) { return tmech::dcontract(B4, x); } + static tensor4 partial(tensor2 const &) { return B4; } +}; +struct project_c { // scalar of tensor2, linear + using in = tensor2; + using out = T; + static T value(tensor2 const &x) { return tmech::dcontract(C2, x); } + static tensor2 partial(tensor2 const &) { return C2; } +}; +struct dyad { // tensor4 of tensor2: Y = tr(x) C (x) x + using in = tensor2; + using out = tensor4; + static tensor4 value(tensor2 const &x) { return tmech::trace(x) * tmech::otimes(C2, x); } + static tensor6 partial(tensor2 const &x) { + return tensor6{tmech::trace(x) * tmech::otimes(C2, Is) + tmech::otimes(tensor4{tmech::otimes(C2, x)}, I)}; + } +}; +struct contract_e { // tensor2 of tensor4: z = Y : E + using in = tensor4; + using out = tensor2; + static tensor2 value(tensor4 const &Y) { return tmech::dcontract(Y, E2); } + static tensor6 partial(tensor4 const &) { return tmech::otimes(Is, E2); } +}; +struct constant_law { // tensor2 of tensor2, independent of its read + using in = tensor2; + using out = tensor2; + static tensor2 value(tensor2 const &) { return B2; } + static tensor4 partial(tensor2 const &) { return tensor4{}; } +}; +// y = c tr(x)^3, about 1e8 at a strain of 1e-3; the partial is off by `error` (relative). +template +struct cubic_strain { + using in = tensor2; + using out = T; + static constexpr T c{1e19}; + static T value(tensor2 const &x) { return c * std::pow(tmech::trace(x), 3); } + static tensor2 partial(tensor2 const &x) { + return tensor2{(1 + Error * T{1e-6}) * 3 * c * std::pow(tmech::trace(x), 2) * I}; + } +}; +struct throwing_map { // throws once the read leaves the state it was armed at + using in = tensor2; + using out = tensor2; + static inline tensor2 armed{}; + static inline bool active{false}; + static tensor2 value(tensor2 const &x) { + if (active && tmech::norm(tensor2{x - armed}) > 0) + throw std::runtime_error("perturbed state"); + return tmech::dcontract(A4, x); + } + static tensor4 partial(tensor2 const &) { return A4; } +}; +struct wrong_map { // declares A4 for B4 : x + using in = tensor2; + using out = tensor2; + static tensor2 value(tensor2 const &x) { return tmech::dcontract(B4, x); } + static tensor4 partial(tensor2 const &) { return A4; } +}; + +/// y = Law::value(x) for x read from "x_source"; "partials" is "declared" (the default), "independent" or "none". +template +class law_node final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit law_node(Args &&...args) + : base(std::forward(args)...), m_y(base::template add_output("y")), + m_source(nm::connection_source::parse(base::template get_parameter("x_source"))), + m_x(base::template add_input(m_source.material, m_source.property, nm::EdgeKind::Global)), + m_mode(base::template get_parameter("partials")) { + auto handle{base::add_callback(&law_node::compute, m_y)}; + handle.reads(m_x); + if (m_mode == "declared") + handle.partials(m_y, m_x); + else if (m_mode == "independent") + handle.independent_of(m_x); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::required(para, "x_source"); + nm::optional(para, "partials", std::string{"declared"}); + return para; + } + + static inline std::function hook; ///< called by every run, when set + + void compute() { + ++m_runs; + if (hook) + hook(); + m_y = Law::value(m_x.get()); + if (base::want_partials()) { + base::partial(m_y, m_x) = Law::partial(m_x.get()); + ++m_fills; + } + } + + int runs() const noexcept { return m_runs; } + int fills() const noexcept { return m_fills; } + +private: + typename Law::out &m_y; + nm::connection_source const m_source; + nm::input_property const &m_x; + std::string const &m_mode; + int m_runs{0}; + int m_fills{0}; +}; + +/// Declares its partials wrongly, as "mistake" says. +template +class misdeclared final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit misdeclared(Args &&...args) + : base(std::forward(args)...), m_y(base::template add_output("y")), + m_z(base::template add_output("z")), + m_x(base::template add_input("strain", "strain", nm::EdgeKind::Global)), + m_w(base::template add_input("temperature", "state", nm::EdgeKind::Global)) { + auto const &mistake{base::template get_parameter("mistake")}; + base::add_callback(&misdeclared::compute_z, m_z).reads(m_x); + auto handle{base::add_callback(&misdeclared::compute, m_y)}; + if (mistake == "before_reads") + handle.partials(m_y, m_x); + handle.reads(m_x); + if (mistake == "unread") + handle.partials(m_y, m_w); + else if (mistake == "twice") + handle.partials(m_y, m_x).partials(m_y, m_x); + else if (mistake == "both") + handle.partials(m_y, m_x).independent_of(m_x); + else if (mistake == "other_output") + handle.partials(m_z, m_x); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::required(para, "mistake"); + return para; + } + + void compute() { m_y = m_x.get(); } + void compute_z() { m_z = m_x.get(); } + +private: + tensor2 &m_y; + tensor2 &m_z; + nm::input_property const &m_x; + nm::input_property const &m_w; +}; + +/// y = a::y + b::y, both read and differentiated. +template +class sum_node final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit sum_node(Args &&...args) + : base(std::forward(args)...), m_y(base::template add_output("y")), + m_a(base::template add_input("a", "y", nm::EdgeKind::Global)), + m_b(base::template add_input("b", "y", nm::EdgeKind::Global)) { + base::add_callback(&sum_node::compute, m_y).reads(m_a, m_b).partials(m_y, m_a, m_b); + } + + static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; } + + void compute() { + m_y = m_a.get() + m_b.get(); + if (base::want_partials()) { + base::partial(m_y, m_a) = Is; + base::partial(m_y, m_b) = Is; + } + } + +private: + tensor2 &m_y; + nm::input_property const &m_a; + nm::input_property const &m_b; +}; + +/// y = (a::y : a::y) b::y: the partial with respect to b depends on a. +template +class product_node final : public nm::material_base, Traits> { +public: + using base = nm::material_base, Traits>; + using input_parameter_controller = typename base::input_parameter_controller; + + template + explicit product_node(Args &&...args) + : base(std::forward(args)...), m_y(base::template add_output("y")), + m_a(base::template add_input("a", "y", nm::EdgeKind::Global)), + m_b(base::template add_input("b", "y", nm::EdgeKind::Global)) { + base::add_callback(&product_node::compute, m_y).reads(m_a, m_b).partials(m_y, m_a, m_b); + } + + static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; } + + void compute() { + auto const aa{tmech::dcontract(m_a.get(), m_a.get())}; + m_y = aa * m_b.get(); + if (base::want_partials()) { + base::partial(m_y, m_a) = tensor4{T{2} * tmech::otimes(m_b.get(), m_a.get())}; + base::partial(m_y, m_b) = tensor4{aa * Is}; + } + } + +private: + tensor2 &m_y; + nm::input_property const &m_a; + nm::input_property const &m_b; +}; + +/// A context with the sources "strain" and "temperature" and law nodes added by node(). +struct chain { + std::unique_ptr owner{std::make_unique()}; + ctx_type &ctx{*owner}; + + chain() { + param_type p; + p.insert("name", "strain"); + ctx.create>(p); + p.clear(); + p.insert("name", "temperature"); + ctx.create>(p); + } + + template + law_node &node(std::string const &name, std::string const &x_source, std::string const &partials = "declared") { + param_type p; + p.insert("name", name); + p.insert("x_source", x_source); + p.insert("partials", partials); + return ctx.create>(p); + } + + void load(tensor2 const &strain, T temperature) { + dynamic_cast *>(ctx.find("strain"))->bind(tensor2{}, strain); + dynamic_cast *>(ctx.find("temperature"))->bind(T{0}, temperature); + ctx.update(); + } +}; + +tensor2 const x0{sym(0.3, -0.1, 0.2, 0.05, -0.15, 0.25)}; +constexpr T s0{1.7}; + +template +T distance(A const &a, B const &b) { + if constexpr (std::is_arithmetic_v) + return std::abs(a - b); + else + return tmech::norm(a - b); +} + +// --- the contraction table, one hand-made chain per row --- + +TEST(PartialDerivatives, ScalarOfScalarOfTensor) { + chain c; + c.node("y", "strain::strain"); + c.node("z", "y::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "strain", "strain")}; + c.load(x0, s0); + tensor2 const expected{2 * norm_trace::value(x0) * norm_trace::partial(x0)}; + EXPECT_LT(distance(d, expected), 1e-12); +} + +TEST(PartialDerivatives, TensorOfScalarOfTensor) { + chain c; + c.node("y", "strain::strain"); + c.node("z", "y::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "strain", "strain")}; + c.load(x0, s0); + EXPECT_LT(distance(d, tensor4{tmech::otimes(B2, norm_trace::partial(x0))}), 1e-12); +} + +TEST(PartialDerivatives, TensorOfTensorOfTensor) { + chain c; + c.node("y", "strain::strain"); + c.node("z", "y::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "strain", "strain")}; + c.load(x0, s0); + EXPECT_LT(distance(d, tensor4{tmech::dcontract(B4, A4)}), 1e-12); +} + +TEST(PartialDerivatives, ScalarOfTensorOfTensor) { + chain c; + c.node("y", "strain::strain"); + c.node("z", "y::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "strain", "strain")}; + c.load(x0, s0); + EXPECT_LT(distance(d, tensor2{tmech::dcontract(C2, A4)}), 1e-12); +} + +TEST(PartialDerivatives, TensorOfTensorOfScalar) { + chain c; + c.node("y", "temperature::state"); + c.node("z", "y::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "temperature", "state")}; + c.load(x0, s0); + EXPECT_LT(distance(d, tensor2{tmech::dcontract(A4, B2)}), 1e-12); +} + +TEST(PartialDerivatives, TensorOfFourthOrderTensorOfTensor) { + chain c; + c.node("y", "strain::strain"); + c.node("z", "y::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "strain", "strain")}; + c.load(x0, s0); + // z = tr(x) (x : E) C + tensor4 const expected{tmech::trace(x0) * tmech::otimes(C2, E2) + tmech::dcontract(x0, E2) * tmech::otimes(C2, I)}; + EXPECT_LT(distance(d, expected), 1e-12); +} + +TEST(PartialDerivatives, ScalarOfScalarIsAScalar) { + chain c; + c.node("y", "temperature::state"); + c.node("z", "y::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "temperature", "state")}; + c.load(x0, s0); + EXPECT_NEAR(d, 4 * s0 * s0 * s0, 1e-12); // z = s^4 +} + +TEST(PartialDerivatives, ParallelPathsAreSummed) { + chain c; + c.node("a", "strain::strain"); + c.node("b", "strain::strain"); + param_type p; + p.insert("name", "s"); + c.ctx.create>(p); + c.node("z", "s::y"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("z", "y", "strain", "strain")}; + c.load(x0, s0); + EXPECT_LT(distance(d, tensor2{tmech::dcontract(C2, tensor4{A4 + B4})}), 1e-12); +} + +TEST(PartialDerivatives, NoPathGivesZero) { + chain c; + c.node("y", "strain::strain"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("y", "y", "temperature", "state")}; + c.load(x0, s0); + EXPECT_EQ(tmech::norm(d), 0); +} + +TEST(PartialDerivatives, TheResultFollowsTheState) { + chain c; + c.node("y", "strain::strain"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("y", "y", "strain", "strain")}; + for (T const f : {1.0, -2.0, 0.5}) { + c.load(tensor2{f * x0}, s0); + EXPECT_LT(distance(d, norm_trace::partial(tensor2{f * x0})), 1e-12); + } +} + +// --- declarations and their errors --- + +TEST(PartialDerivatives, AMissingPartialOnARequestedPathThrowsAndNamesThePath) { + chain c; + c.node("a", "strain::strain"); + c.node("b", "a::y", "none"); + c.node("c", "b::y"); + c.ctx.finalize(); + try { + c.ctx.derivative("c", "y", "strain", "strain"); + FAIL() << "a requested path without a partial was accepted"; + } catch (std::logic_error const &e) { + std::string const what{e.what()}; + EXPECT_TRUE(what.contains("'b::y'")) << what; + EXPECT_TRUE(what.contains("strain::strain -> a::y -> b::y -> c::y")) << what; + } +} + +TEST(PartialDerivatives, AMissingPartialOffTheRequestedPathsDoesNotThrow) { + chain c; + c.node("a", "strain::strain"); + c.node("b", "a::y", "none"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("a", "y", "strain", "strain")}; + EXPECT_NO_THROW(c.load(x0, s0)); + EXPECT_LT(distance(d, A4), 1e-12); +} + +TEST(PartialDerivatives, IndependentOfDeclaresAZeroPartial) { + chain c; + c.node("a", "strain::strain"); + c.node("b", "a::y", "independent"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("b", "y", "strain", "strain")}; + c.load(x0, s0); + EXPECT_EQ(tmech::norm(d), 0); +} + +TEST(PartialDerivatives, TheResultTypeMustMatchTheRanks) { + chain c; + c.node("a", "strain::strain"); + c.ctx.finalize(); + EXPECT_THROW(c.ctx.derivative("a", "y", "strain", "strain"), std::logic_error); + EXPECT_THROW(c.ctx.derivative("a", "y", "strain", "nothing"), std::runtime_error); +} + +TEST(PartialDerivatives, PartialsAreDeclaredForReadsAndOutputsOfTheCallback) { + for (auto const *mistake : {"unread", "twice", "both", "other_output", "before_reads"}) { + chain c; + param_type p; + p.insert("name", "m"); + p.insert("mistake", mistake); + EXPECT_THROW(c.ctx.create>(p), std::logic_error) << mistake; + } +} + +// --- independent_of() is verified, not trusted (review of #227) --- + +TEST(PartialDerivatives, AWrongIndependentOfFailsTheFiniteDifferenceCheck) { + chain c; + c.node("a", "strain::strain"); + c.node("liar", "a::y", "independent"); // B4 : x does depend on x + c.ctx.finalize(); + c.load(x0, s0); + auto const mismatches{nm::check_partials(c.ctx, "liar")}; + ASSERT_FALSE(mismatches.empty()); + EXPECT_TRUE(mismatches.front().contains("liar::y")) << mismatches.front(); + EXPECT_TRUE(mismatches.front().contains("a::y")) << mismatches.front(); + EXPECT_TRUE(mismatches.front().contains("independent_of")) << mismatches.front(); + EXPECT_LT(distance(c.ctx.get("liar", "y"), map_b::value(c.ctx.get("a", "y"))), 1e-15) << "the state is restored"; +} + +TEST(PartialDerivatives, ATrueIndependentOfPassesTheFiniteDifferenceCheck) { + chain c; + c.node("a", "strain::strain"); + c.node("k", "a::y", "independent"); + c.ctx.finalize(); + c.load(x0, s0); + EXPECT_TRUE(nm::check_partials(c.ctx, "k").empty()); +} + +// --- the step follows the read, the allowance the analytic block (review of #227) --- + +TEST(PartialDerivatives, ACorrectPartialPassesAtTheScaleOfAStrain) { + chain c; + c.node>("n", "strain::strain"); + c.ctx.finalize(); + c.load(tensor2{1e-3 * x0}, s0); + ASSERT_GT(std::abs(c.ctx.get("n", "y")), 1e7); + auto const mismatches{nm::check_partials(c.ctx, "n")}; + EXPECT_TRUE(mismatches.empty()) << mismatches.front(); +} + +TEST(PartialDerivatives, ASmallRelativeErrorIsDetectedAtOutputsOfOneHundredMillion) { + chain c; + c.node>("n", "strain::strain"); // 2e-5 relative + c.ctx.finalize(); + c.load(tensor2{1e-3 * x0}, s0); + ASSERT_GT(std::abs(c.ctx.get("n", "y")), 1e7); + EXPECT_FALSE(nm::check_partials(c.ctx, "n").empty()); +} + +// --- a throwing rerun leaves the context as it was (review of #227) --- + +TEST(PartialDerivatives, AThrowingPerturbedRerunRestoresReadsAndOutputs) { + chain c; + c.node("t", "strain::strain"); + c.ctx.finalize(); + c.load(x0, s0); + throwing_map::armed = x0; + throwing_map::active = true; + EXPECT_THROW(nm::check_partials(c.ctx, "t"), std::runtime_error); + throwing_map::active = false; + EXPECT_EQ(tmech::norm(tensor2{dynamic_cast *>(c.ctx.find("strain"))->strain().new_value() - x0}), 0) << "the read"; + EXPECT_LT(distance(c.ctx.get("t", "y"), throwing_map::value(x0)), 1e-15) << "the output"; +} + +// --- the result type is the one the context can fill (review of #227) --- + +template +concept derivable = requires(ctx_type &c, std::string const &name) { c.template derivative(name, name, name, name); }; + +TEST(PartialDerivatives, OnlyTheScalarAndTensorsOfTheContextAreDerivativeTypes) { + static_assert(derivable); + static_assert(derivable); + static_assert(derivable); + static_assert(derivable); + static_assert(!derivable>); // another dimension would write out of bounds + static_assert(!derivable>); + static_assert(!derivable); + static_assert(!derivable); + static_assert(!derivable); +} + + +TEST(PartialDerivatives, AFourthOrderResultOfATensorIsRankSix) { + chain c; + c.node("d", "strain::strain"); + c.ctx.finalize(); + auto const &d{c.ctx.derivative("d", "y", "strain", "strain")}; + c.load(x0, s0); + EXPECT_LT(distance(d, dyad::partial(x0)), 1e-12); +} + +// --- derivative() on an evaluated context, and the path query (review of #227) --- + +TEST(PartialDerivatives, ADerivativeRequestedAfterAnUpdateIsFilledAtOnce) { + chain c; + c.node("a", "strain::strain"); + c.node("b", "a::y"); + c.ctx.finalize(); + c.load(x0, s0); + auto const &d{c.ctx.derivative("b", "y", "strain", "strain")}; + EXPECT_LT(distance(d, tensor4{tmech::dcontract(B4, A4)}), 1e-12); + c.load(tensor2{2 * x0}, s0); + EXPECT_LT(distance(d, tensor4{tmech::dcontract(B4, A4)}), 1e-12) << "and by the next update"; +} + +TEST(PartialDerivatives, ADerivativeFromACallbackIsRefused) { + chain c; + c.node("a", "strain::strain"); + c.ctx.finalize(); + bool refused{false}; + law_node::hook = [&] { + try { + c.ctx.derivative("a", "y", "strain", "strain"); + } catch (std::logic_error const &) { + refused = true; + } + }; + c.load(x0, s0); + law_node::hook = nullptr; + EXPECT_TRUE(refused); +} + +TEST(PartialDerivatives, HasPathTellsIndependenceFromATypo) { + chain c; + c.node("y", "strain::strain"); + c.ctx.finalize(); + EXPECT_TRUE(c.ctx.has_path("y", "y", "strain", "strain")); + EXPECT_TRUE(c.ctx.has_path("y", "y", "y", "y")); + EXPECT_FALSE(c.ctx.has_path("y", "y", "temperature", "state")); + EXPECT_FALSE(c.ctx.has_path("strain", "strain", "y", "y")); + EXPECT_THROW(c.ctx.has_path("y", "y", "temperature", "nothing"), std::runtime_error); +} + +// --- zero cost when nothing is requested --- + +TEST(PartialDerivatives, NoPartialsAreComputedWithoutARequest) { + chain c; + auto &a{c.node("a", "strain::strain")}; + auto &b{c.node("b", "a::y")}; + auto &q{c.node("q", "strain::strain")}; + c.ctx.finalize(); + for (int i{0}; i < 3; ++i) + c.load(x0, s0); + EXPECT_EQ(a.runs(), 3); + EXPECT_EQ(a.fills() + b.fills() + q.fills(), 0); + + c.ctx.derivative("b", "y", "strain", "strain"); // fills once at the held state + EXPECT_EQ(a.fills(), 1); + EXPECT_EQ(b.fills(), 1); + for (int i{0}; i < 2; ++i) + c.load(x0, s0); + EXPECT_EQ(a.fills(), 3); + EXPECT_EQ(b.fills(), 3); + EXPECT_EQ(q.fills(), 0) << "q is off the requested path"; +} + +// --- the per-node finite-difference check --- + +TEST(PartialDerivatives, TheFiniteDifferenceCheckPassesCorrectPartials) { + chain c; + c.node("n", "strain::strain"); + c.node("s", "n::y"); + c.node("sb", "s::y"); + c.node("a", "sb::y"); + c.node("p", "a::y"); + c.node("d", "a::y"); + c.node("e", "d::y"); + c.ctx.finalize(); + c.load(x0, s0); + for (auto const *name : {"n", "s", "sb", "a", "p", "d", "e"}) + EXPECT_TRUE(nm::check_partials(c.ctx, name).empty()) << name; +} + +TEST(PartialDerivatives, TheFiniteDifferenceCheckComparesWithThePartialsAtTheState) { + chain c; + c.node("a", "strain::strain"); + c.node("b", "temperature::state"); + param_type p; + p.insert("name", "q"); + c.ctx.create>(p); + c.ctx.finalize(); + c.ctx.derivative("q", "y", "strain", "strain"); // keeps the partials of q requested + c.load(x0, s0); + auto const mismatches{nm::check_partials(c.ctx, "q")}; + EXPECT_TRUE(mismatches.empty()) << mismatches.front(); +} + +TEST(PartialDerivatives, TheFiniteDifferenceCheckNamesAWrongPartial) { + chain c; + c.node("w", "strain::strain"); + c.ctx.finalize(); + c.load(x0, s0); + auto const mismatches{nm::check_partials(c.ctx, "w")}; + ASSERT_FALSE(mismatches.empty()); + EXPECT_TRUE(mismatches.front().contains("w::y")) << mismatches.front(); + EXPECT_TRUE(mismatches.front().contains("strain::strain")) << mismatches.front(); + EXPECT_LT(distance(c.ctx.get("w", "y"), wrong_map::value(x0)), 1e-15) << "the check restores the state"; +} + +// --- thermal eigenstrain with linear elasticity: the first DDSDDT --- + +constexpr T K{100}, G{40}, alpha{1e-5}, T_ref{20}; + +struct thermoelastic { + std::unique_ptr owner{std::make_unique()}; + ctx_type &ctx{*owner}; + std::unique_ptr> point; + + thermoelastic() { + param_type p; + p.insert("name", "strain"); + ctx.create>(p); + p.clear(); + p.insert("name", "temperature"); + ctx.create>(p); + p.clear(); + p.insert("name", "alpha"); + p.insert("value", alpha); + ctx.create>(p); + p.clear(); + p.insert("name", "thermal"); + p.insert("alpha_source", "alpha::value"); + p.insert("reference_temperature", T_ref); + ctx.create>(p); + p.clear(); + p.insert("name", "elastic"); + p.insert("strain_producer_name", "strain"); + p.insert("eigenstrain_source", "thermal"); + p.insert("K", K); + p.insert("G", G); + ctx.create>(p); + ctx.finalize(); + param_type q; + q.insert>("inputs", {"strain", "temperature"}); + point = std::make_unique>(ctx, q); + } + + void evaluate(tensor2 const &strain, T temperature) { + point->reset(); + point->input>("strain").bind(tensor2{}, strain); + point->input>("temperature").bind(T_ref, temperature); + point->evaluate(); + } +}; + +TEST(PartialDerivatives, StressWithRespectToTemperatureMatchesFiniteDifferences) { + thermoelastic m; + auto const &ddsddt{m.point->derivative("elastic::stress", "temperature")}; + tensor2 const eps{sym(1e-3, -2e-4, 5e-4, 1e-4, -3e-4, 2e-4)}; + T const theta{80}, h{1e-3}; + m.evaluate(eps, theta + h); + tensor2 const plus{m.point->output("elastic::stress")}; + m.evaluate(eps, theta - h); + tensor2 const minus{m.point->output("elastic::stress")}; + m.evaluate(eps, theta); + tensor2 const fd{(plus - minus) / (2 * h)}; + EXPECT_LT(tmech::norm(tensor2{ddsddt - fd}), 1e-8 * tmech::norm(fd)); + EXPECT_LT(tmech::norm(tensor2{ddsddt + 3 * K * alpha * I}), 1e-12) << "-C : alpha I"; +} + +TEST(PartialDerivatives, StressWithRespectToStrainEqualsTheTangent) { + thermoelastic m; + auto const &ddsdde{m.point->derivative("elastic::stress", "strain")}; + m.evaluate(sym(1e-3, -2e-4, 5e-4, 1e-4, -3e-4, 2e-4), 60); + auto const &tangent{m.point->output("elastic::tangent")}; + EXPECT_LT(tmech::norm(tensor4{ddsdde - tangent}), 1e-12 * tmech::norm(tangent)); +} + +TEST(PartialDerivatives, TheMaterialPointNamesUnknownInputs) { + thermoelastic m; + EXPECT_THROW(m.point->derivative("elastic::stress", "pressure"), nm::umat::fatal_error); + EXPECT_THROW(m.point->derivative("elastic::nothing", "temperature"), nm::umat::fatal_error); +} + +TEST(PartialDerivatives, ThermoelasticPartialsPassTheFiniteDifferenceCheck) { + thermoelastic m; + m.evaluate(sym(1e-3, -2e-4, 5e-4, 1e-4, -3e-4, 2e-4), 60); + for (auto const *name : {"thermal", "elastic"}) + for (auto const &mismatch : nm::check_partials(m.ctx, name)) + ADD_FAILURE() << mismatch; +} + +} // namespace diff --git a/tests/test_residual_system.cpp b/tests/test_residual_system.cpp new file mode 100644 index 0000000..29bb302 --- /dev/null +++ b/tests/test_residual_system.cpp @@ -0,0 +1,1100 @@ +// Residual systems set up at run time: unknowns and equations from any nodes (#218). +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/errors.h" +#include "numsim-materials/core/material_base.h" +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/core/parameter_helpers.h" +#include "numsim-materials/core/tensor_identities.h" +#include "numsim-materials/default_materials.h" +#include "numsim-materials/driver/material_point.h" +#include "numsim-materials/solvers/system.h" +#include "numsim-materials/solvers/vector_newton.h" +#include "numsim-materials/umat/external_state_source.h" + +#if __has_include() +#define NUMSIM_HAVE_JSON 1 +#include +#include "numsim-materials/io/json_material_factory.h" +#endif + +namespace { + +namespace nm = numsim::materials; +using policy = nm::material_policy_default; +using T = policy::value_type; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; +using tensor2 = tmech::tensor; +using tensor4 = tmech::tensor; +using system_type = nm::residual_system; +using solver_type = nm::vector_newton; + +enum class system_mode { linear, nonlinear, singular, decoupled }; + +/// The rows "x" ("rx"), "y" ("ry") or both of the two-scalar systems of test_vector_newton, with their partials. +/// The x row of the decoupled system does not read y at all. +class scalar_rows final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit scalar_rows(Args &&...args) + : base(std::forward(args)...), m_mode(static_cast(base::get_parameter("mode"))), + m_a(base::get_parameter("a")), m_b(base::get_parameter("b")), m_rows(base::get_parameter("rows")), + m_x(base::add_input(base::get_parameter("solver"), "x")), + m_y(m_rows == "x" && m_mode == system_mode::decoupled ? nullptr + : &base::add_input(base::get_parameter("solver"), "y")) { + if (m_rows != "y") + m_rx = &base::add_output("rx"); + if (m_rows != "x") + m_ry = &base::add_output("ry"); + auto handle{base::add_callback(&scalar_rows::compute, m_rx, m_ry)}; + handle.reads(m_x, m_y); + if (m_rx) + handle.partials(*m_rx, m_x, m_y); + if (m_ry) + handle.partials(*m_ry, m_x, m_y); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "solver", std::string{"solver"}); + nm::optional(para, "rows", std::string{"xy"}); + nm::optional(para, "mode", 0); + nm::optional(para, "a", T{0}); + nm::optional(para, "b", T{0}); + return para; + } + + void compute() { + T const x{m_x.get()}, y{m_y ? m_y->get() : T{0}}; + T rx{}, ry{}, jxx{}, jxy{}, jyx{}, jyy{}; + switch (m_mode) { + case system_mode::linear: + rx = 2 * x + y - m_a, jxx = 2, jxy = 1; + ry = x + 3 * y - m_b, jyx = 1, jyy = 3; + break; + case system_mode::nonlinear: + rx = x * x + y - m_a, jxx = 2 * x, jxy = 1; + ry = x + y * y - m_b, jyx = 1, jyy = 2 * y; + break; + case system_mode::singular: + rx = x + y - m_a, jxx = 1, jxy = 1; + ry = x + y - m_b, jyx = 1, jyy = 1; + break; + case system_mode::decoupled: + rx = 2 * x - m_a, jxx = 2, jxy = 0; + ry = x + 3 * y - m_b, jyx = 1, jyy = 3; + break; + } + if (m_rx) + *m_rx = rx; + if (m_ry) + *m_ry = ry; + if (!base::want_partials()) + return; + if (m_rx) { + base::partial(*m_rx, m_x) = jxx; + if (m_y) + base::partial(*m_rx, *m_y) = jxy; + } + if (m_ry) { + base::partial(*m_ry, m_x) = jyx; + base::partial(*m_ry, *m_y) = jyy; + } + } + +private: + system_mode m_mode; + T const &m_a; + T const &m_b; + std::string const &m_rows; + nm::input_property const &m_x; + nm::input_property const *m_y; + T *m_rx{nullptr}; + T *m_ry{nullptr}; +}; + +/// R_g = g + tr(B) - a, R_B = B - g P + c g^2 I, as mixed_system_t of test_vector_newton, with its partials. +template +class mixed_rows final : public nm::material_base, Policy> { +public: + using base = nm::material_base, Policy>; + using input_parameter_controller = typename base::input_parameter_controller; + static constexpr auto D = base::Dim; + using t2 = tmech::tensor; + using t4 = tmech::tensor; + + template + explicit mixed_rows(Args &&...args) + : base(std::forward(args)...), m_rg(base::template add_output("rg")), m_rB(base::template add_output("rB")), + m_a(base::template get_parameter("a")), m_c(base::template get_parameter("c")), + m_g(base::template add_input("solver", "g")), m_B(base::template add_input("solver", "B")) { + if constexpr (D == 2) + m_P = t2{0.5, 0.1, 0.1, 0.7}; + else + m_P = t2{0.5, 0.1, 0.2, 0.1, 0.7, 0.3, 0.2, 0.3, 0.9}; + base::add_callback(&mixed_rows::compute, m_rg, m_rB).reads(m_g, m_B).partials(m_rg, m_g, m_B).partials(m_rB, m_g, m_B); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "a", T{1}); + nm::optional(para, "c", T{0}); + return para; + } + + t2 const &P() const noexcept { return m_P; } + + void compute() { + auto const g{m_g.get()}; + auto const &B{m_B.get()}; + auto const I{tmech::eye()}; + m_rg = g + tmech::trace(B) - m_a; + m_rB = B - g * m_P + (m_c * g * g) * I; + if (!base::want_partials()) + return; + base::partial(m_rg, m_g) = T{1}; + base::partial(m_rg, m_B) = t2{I}; + base::partial(m_rB, m_g) = t2{-m_P + (2 * m_c * g) * I}; + base::partial(m_rB, m_B) = nm::sym_identity(); + } + +private: + T &m_rg; + t2 &m_rB; + T const &m_a; + T const &m_c; + nm::input_property const &m_g; + nm::input_property const &m_B; + t2 m_P{}; +}; + +/// R_g = g + tr(B) - a - tr(eps), R_B = B - g P, with eps from strain::strain: strain_coupled_system of +/// test_vector_newton with the strain as a seed instead of a member. +class strain_rows final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit strain_rows(Args &&...args) + : base(std::forward(args)...), m_rg(base::add_output("rg")), m_rB(base::add_output("rB")), + m_a(base::get_parameter("a")), m_g(base::add_input("solver", "g")), m_B(base::add_input("solver", "B")), + m_eps(base::add_input("strain", "strain")) { + base::add_callback(&strain_rows::compute, m_rg, m_rB) + .reads(m_g, m_B, m_eps) + .partials(m_rg, m_g, m_B, m_eps) + .partials(m_rB, m_g, m_B, m_eps); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "a", T{2}); + return para; + } + + static tensor2 P() { return tensor2{0.5, 0.2, 0.1, 0.2, 0.7, 0.3, 0.1, 0.3, 0.9}; } + + void compute() { + auto const I{tmech::eye()}; + m_rg = m_g.get() + tmech::trace(m_B.get()) - m_a - tmech::trace(m_eps.get()); + m_rB = m_B.get() - m_g.get() * P(); + if (!base::want_partials()) + return; + base::partial(m_rg, m_g) = T{1}; + base::partial(m_rg, m_B) = tensor2{I}; + base::partial(m_rg, m_eps) = tensor2{-I}; + base::partial(m_rB, m_g) = tensor2{-P()}; + base::partial(m_rB, m_B) = nm::sym_identity(); + base::partial(m_rB, m_eps) = tensor4{}; + } + +private: + T &m_rg; + tensor2 &m_rB; + T const &m_a; + nm::input_property const &m_g; + nm::input_property const &m_B; + nm::input_property const &m_eps; +}; + +/// r = x^3 + x - load::state for the unknown x of "solver"; with "picard", r = x - cos(load + x) / 2. +class cubic_row final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit cubic_row(Args &&...args) + : base(std::forward(args)...), m_r(base::add_output("r")), m_picard(base::get_parameter("picard")), + m_x(base::add_input("solver", "x")), m_load(base::add_input("load", "state")) { + auto handle{base::add_callback(&cubic_row::compute, m_r).reads(m_x, m_load)}; + if (base::get_parameter("declares_partials")) + handle.partials(m_r, m_x, m_load); + if (base::get_parameter("counts")) + handle.stateful(); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "picard", false); + nm::optional(para, "declares_partials", true); + nm::optional(para, "counts", false); + return para; + } + + void compute() { + auto const x{m_x.get()}, load{m_load.get()}; + m_r = m_picard ? x - std::cos(load + x) / 2 : x * x * x + x - load; + if (base::want_partials()) { + base::partial(m_r, m_x) = m_picard ? 1 + std::sin(load + x) / 2 : 3 * x * x + 1; + base::partial(m_r, m_load) = m_picard ? std::sin(load + x) / 2 : T{-1}; + } + } + +private: + T &m_r; + bool const &m_picard; + nm::input_property const &m_x; + nm::input_property const &m_load; +}; + +/// r = x^2 - load::state for the unknown x of "solver": a double root, and a singular Jacobian, at load 0. +class square_row final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit square_row(Args &&...args) + : base(std::forward(args)...), m_r(base::add_output("r")), m_x(base::add_input("solver", "x")), + m_load(base::add_input("load", "state")) { + base::add_callback(&square_row::compute, m_r).reads(m_x, m_load).partials(m_r, m_x, m_load); + } + + static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; } + + void compute() { + m_r = m_x.get() * m_x.get() - m_load.get(); + if (base::want_partials()) { + base::partial(m_r, m_x) = 2 * m_x.get(); + base::partial(m_r, m_load) = T{-1}; + } + } + +private: + T &m_r; + nm::input_property const &m_x; + nm::input_property const &m_load; +}; + +/// r = t^3 + t - 2 with t = x / scale for the unknown x of "solver": the root x = scale, whatever the scale. +class scaled_row final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit scaled_row(Args &&...args) + : base(std::forward(args)...), m_r(base::add_output("r")), m_scale(base::get_parameter("scale")), + m_x(base::add_input("solver", "x")) { + base::add_callback(&scaled_row::compute, m_r).reads(m_x).partials(m_r, m_x); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "scale", T{1}); + return para; + } + + void compute() { + T const t{m_x.get() / m_scale}; + m_r = t * t * t + t - 2; + if (base::want_partials()) + base::partial(m_r, m_x) = (3 * t * t + 1) / m_scale; + } + +private: + T &m_r; + T const &m_scale; + nm::input_property const &m_x; +}; + +/// r = "value" whatever x is: a residual that is NaN or infinite. +class poison_row final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit poison_row(Args &&...args) + : base(std::forward(args)...), m_r(base::add_output("r")), m_value(base::get_parameter("value")), + m_x(base::add_input("solver", "x")) { + base::add_callback(&poison_row::compute, m_r).reads(m_x); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "value", T{0}); + return para; + } + + void compute() { m_r = m_value + (m_x.get() - m_x.get()); } + +private: + T &m_r; + T const &m_value; + nm::input_property const &m_x; +}; + +/// r = y^3 + y - solver::x for the unknown y of "second": a system reading the unknown of another. +class chain_row final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit chain_row(Args &&...args) + : base(std::forward(args)...), m_r(base::add_output("r")), m_y(base::add_input("second", "y")), + m_x(base::add_input("solver", "x")) { + base::add_callback(&chain_row::compute, m_r).reads(m_y, m_x).partials(m_r, m_y, m_x); + } + + static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; } + + void compute() { + T const y{m_y.get()}; + m_r = y * y * y + y - m_x.get(); + if (base::want_partials()) { + base::partial(m_r, m_y) = 3 * y * y + 1; + base::partial(m_r, m_x) = T{-1}; + } + } + +private: + T &m_r; + nm::input_property const &m_y; + nm::input_property const &m_x; +}; + +/// A residual system "solver" over `unknowns` and `equations`. +system_type &make_system(ctx_type &ctx, std::vector unknowns, std::vector equations, + std::string const &method = "newton", int max_iter = 50, T tolerance = 1e-10, + std::string const &name = "solver") { + param_type p; + p.insert("name", name); + p.insert("method", method); + p.insert>("unknowns", std::move(unknowns)); + p.insert>("equations", std::move(equations)); + p.insert("max_iter", max_iter); + p.insert("tolerance", tolerance); + return ctx.create(p); +} + +std::vector xy() { return {{"x", nm::unknown_kind::scalar}, {"y", nm::unknown_kind::scalar}}; } + +/// x, y of the scalar system in one node "sys" (or split into "fx" and "fy"), solved by a system. +struct scalar_rig { + std::unique_ptr owner{std::make_unique()}; + ctx_type &ctx{*owner}; + system_type *system{nullptr}; + + scalar_rig(system_mode mode, T a, T b, int max_iter = 50, bool split = false, std::string const &method = "newton") { + std::vector equations{split ? std::vector{{"fx::rx", "x"}, {"fy::ry", "y"}} + : std::vector{{"sys::ry", "y"}, {"sys::rx", "x"}}}; + system = &make_system(ctx, xy(), equations, method, max_iter); + for (auto const &[name, rows] : split ? std::vector>{{"fx", "x"}, {"fy", "y"}} + : std::vector>{{"sys", "xy"}}) { + param_type p; + p.insert("name", name); + p.insert("rows", rows); + p.insert("mode", static_cast(mode)); + p.insert("a", a); + p.insert("b", b); + ctx.create(p); + } + ctx.finalize(); + } + + T x() { return ctx.get("solver", "x"); } + T y() { return ctx.get("solver", "y"); } +}; + +/// vector_newton on the same scalar system: the root and the fewest max_iter it converges with. +struct newton_reference { + T x{}, y{}; + int iterations{-1}; +}; + +class scalar_system_2 final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit scalar_system_2(Args &&...args) + : base(std::forward(args)...), m_rx(base::add_output("residual_x")), m_ry(base::add_output("residual_y")), + m_jxx(base::add_output("jacobian_x_x")), m_jxy(base::add_output("jacobian_x_y")), + m_jyx(base::add_output("jacobian_y_x")), m_jyy(base::add_output("jacobian_y_y")), + m_rows(static_cast(base::get_parameter("mode")), base::get_parameter("a"), base::get_parameter("b")), + m_x(base::add_input("solver", "x", nm::EdgeKind::Local)), m_y(base::add_input("solver", "y", nm::EdgeKind::Local)) { + base::add_callback(&scalar_system_2::compute, m_rx, m_ry, m_jxx, m_jxy, m_jyx, m_jyy).reads(m_x, m_y); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "mode", 0); + nm::optional(para, "a", T{0}); + nm::optional(para, "b", T{0}); + return para; + } + + void compute() { + auto const [mode, a, b]{m_rows}; + T const x{m_x.get()}, y{m_y.get()}; + switch (mode) { + case system_mode::linear: + m_rx = 2 * x + y - a, m_jxx = 2, m_jxy = 1; + m_ry = x + 3 * y - b, m_jyx = 1, m_jyy = 3; + break; + case system_mode::nonlinear: + m_rx = x * x + y - a, m_jxx = 2 * x, m_jxy = 1; + m_ry = x + y * y - b, m_jyx = 1, m_jyy = 2 * y; + break; + case system_mode::singular: + m_rx = x + y - a, m_jxx = 1, m_jxy = 1; + m_ry = x + y - b, m_jyx = 1, m_jyy = 1; + break; + case system_mode::decoupled: + m_rx = 2 * x - a, m_jxx = 2, m_jxy = 0; + m_ry = x + 3 * y - b, m_jyx = 1, m_jyy = 3; + break; + } + } + +private: + T &m_rx, &m_ry, &m_jxx, &m_jxy, &m_jyx, &m_jyy; + std::tuple m_rows; + nm::input_property const &m_x; + nm::input_property const &m_y; +}; + +newton_reference vector_newton_on(system_mode mode, T a, T b) { + newton_reference ref; + for (int max_iter{0}; max_iter <= 50 && ref.iterations < 0; ++max_iter) { + ctx_type ctx; + param_type p; + p.insert("name", "solver"); + p.insert("function", "sys"); + p.insert("max_iter", max_iter); + p.insert("graph_driven", false); + p.insert>("unknowns", {{"x", nm::unknown_kind::scalar}, {"y", nm::unknown_kind::scalar}}); + auto &solver{ctx.create(p)}; + p.clear(); + p.insert("name", "sys"); + p.insert("mode", static_cast(mode)); + p.insert("a", a); + p.insert("b", b); + ctx.create(p); + ctx.finalize(); + solver.solve(); + if (solver.converged()) + ref = {ctx.get("solver", "x"), ctx.get("solver", "y"), max_iter}; + } + return ref; +} + +/// The message of the exception of type E that `f` throws; empty when it throws none. +template +std::string thrown(F &&f) { + try { + f(); + } catch (E const &e) { + return e.what(); + } + return {}; +} + +// --- the cases of test_vector_newton, as systems --- + +TEST(ResidualSystem, TakesTheIterationsOfVectorNewtonAndLandsOnItsRoot) { + struct problem { + system_mode mode; + T a, b; + }; + for (auto const &[mode, a, b] : {problem{system_mode::linear, 5, 10}, problem{system_mode::nonlinear, 3, 5}, + problem{system_mode::nonlinear, 5.248, 27.4}, problem{system_mode::decoupled, 5, 10}}) { + auto const ref{vector_newton_on(mode, a, b)}; + ASSERT_GE(ref.iterations, 0); + scalar_rig rig{mode, a, b}; + rig.ctx.update(); + EXPECT_EQ(rig.system->iterations(), ref.iterations) << static_cast(mode); + EXPECT_NEAR(rig.x(), ref.x, 1e-14); + EXPECT_NEAR(rig.y(), ref.y, 1e-14); + if (ref.iterations > 0) { + scalar_rig short_of_it{mode, a, b, ref.iterations - 1}; + EXPECT_THROW(short_of_it.ctx.update(), nm::cutback_error) << "max_iter counts Newton updates, as in vector_newton"; + } + } +} + +TEST(ResidualSystem, LinearSystemIsExactAndConvergesOnExactlyMaxIter) { + scalar_rig rig{system_mode::linear, 5, 10, 1}; + EXPECT_NO_THROW(rig.ctx.update()); + EXPECT_EQ(rig.system->iterations(), 1); + EXPECT_NEAR(rig.x(), 1, 1e-12); + EXPECT_NEAR(rig.y(), 3, 1e-12); +} + +TEST(ResidualSystem, NonConvergenceThrowsACutbackNamingTheSystem) { + scalar_rig rig{system_mode::nonlinear, 3, 5, 0}; + auto const message{thrown([&] { rig.ctx.update(); })}; + EXPECT_TRUE(message.contains("'solver'")) << message; +} + +TEST(ResidualSystem, ASingularJacobianThrowsACutbackAndLeavesTheUnknownsFinite) { + scalar_rig rig{system_mode::singular, 4, 7}; + EXPECT_THROW(rig.ctx.update(), nm::cutback_error); + EXPECT_TRUE(std::isfinite(rig.x()) && std::isfinite(rig.y())); +} + +TEST(ResidualSystem, SolvesArePointLocal) { + scalar_rig fresh{system_mode::nonlinear, 3, 5}; + fresh.ctx.update(); + scalar_rig shared{system_mode::nonlinear, 5.248, 27.4}; + shared.ctx.update(); + auto *sys{shared.ctx.find("sys")}; + sys->set_parameter("a", 3); + sys->set_parameter("b", 5); + shared.ctx.update(); + EXPECT_DOUBLE_EQ(shared.x(), fresh.x()) << "the seed is 0, not the previous root"; + EXPECT_DOUBLE_EQ(shared.y(), fresh.y()); +} + +TEST(ResidualSystem, AFailedSolveDoesNotMoveTheNext) { + scalar_rig reference{system_mode::nonlinear, 3, 5}; + reference.ctx.update(); + scalar_rig rig{system_mode::nonlinear, -10, -10}; + EXPECT_THROW(rig.ctx.update(), nm::cutback_error); + rig.ctx.revert(); + rig.ctx.find("sys")->set_parameter("a", 3); + rig.ctx.find("sys")->set_parameter("b", 5); + rig.ctx.update(); + EXPECT_DOUBLE_EQ(rig.x(), reference.x()); +} + +template +void mixed_case(T c, T tolerance) { + using t2 = tmech::tensor; + nm::material_context ctx; + param_type p; + p.insert("name", "solver"); + p.insert>("unknowns", {{"g", nm::unknown_kind::scalar}, {"B", nm::unknown_kind::sym_tensor}}); + p.insert>("equations", {{"sys::rB", "B"}, {"sys::rg", "g"}}); + auto &system{ctx.template create>(p)}; + p.clear(); + p.insert("name", "sys"); + p.insert("a", 2); + p.insert("c", c); + auto &sys{ctx.template create>(p)}; + ctx.finalize(); + ctx.update(); + T const trP{tmech::trace(sys.P())}, s{1 + trP}, D{T(Policy::Dim)}; + T const g{c == 0 ? 2 / s : (s - std::sqrt(s * s - 4 * D * c * 2)) / (2 * D * c)}; + EXPECT_NEAR(ctx.template get("solver", "g"), g, tolerance); + t2 const B{g * sys.P() - (c * g * g) * tmech::eye()}; + EXPECT_TRUE(tmech::almost_equal(ctx.template get("solver", "B"), B, tolerance)); + if (c == 0) + EXPECT_EQ(system.iterations(), 1) << "linear: one update"; + else + EXPECT_EQ(system.iterations(), 4) << "2.0 -> 4.2e-2 -> 2.1e-4 -> 5.6e-9 -> 2.2e-16, as vector_newton measured"; +} + +TEST(ResidualSystem, MixedScalarAndSymmetricTensorSystem) { mixed_case(0, 1e-10); } +TEST(ResidualSystem, NonlinearMixedSystemIteratesQuadratically) { mixed_case(0.1, 1e-12); } +TEST(ResidualSystem, MixedSystemIn2D) { mixed_case(0, 1e-10); } + +// --- equations from any nodes --- + +TEST(ResidualSystem, SplitOverTwoResidualNodesConvergesIdentically) { + for (auto const mode : {system_mode::linear, system_mode::nonlinear, system_mode::decoupled}) { + T const a{mode == system_mode::nonlinear ? T{3} : T{5}}, b{mode == system_mode::nonlinear ? T{5} : T{10}}; + scalar_rig one{mode, a, b}; + scalar_rig split{mode, a, b, 50, true}; + one.ctx.update(); + split.ctx.update(); + EXPECT_EQ(split.system->iterations(), one.system->iterations()); + EXPECT_EQ(split.x(), one.x()); + EXPECT_EQ(split.y(), one.y()); + } +} + +TEST(ResidualSystem, StructuralZeroBlocksAreTheBlocksWithoutAPath) { + scalar_rig split{system_mode::decoupled, 5, 10, 50, true}; + EXPECT_EQ(split.system->zero_blocks(), (std::vector{{"x", "y"}})) + << "what vector_newton declared as zero_blocks: [[\"x\", \"y\"]]"; + EXPECT_TRUE(scalar_rig(system_mode::linear, 5, 10, 50, true).system->zero_blocks().empty()); + EXPECT_TRUE(scalar_rig(system_mode::decoupled, 5, 10).system->zero_blocks().empty()) << "one node reads both"; + split.ctx.update(); + EXPECT_NEAR(split.x(), 2.5, 1e-12); + EXPECT_NEAR(split.y(), 2.5, 1e-12); +} + +TEST(ResidualSystem, NewtonFdFindsTheSameRoot) { + scalar_rig fd{system_mode::nonlinear, 3, 5, 50, false, "newton_fd"}; + fd.ctx.update(); + EXPECT_NEAR(fd.x(), 1, 1e-10); + EXPECT_NEAR(fd.y(), 2, 1e-10); +} + +// --- the implicit function theorem through the system --- + +struct strain_rig { + std::unique_ptr owner{std::make_unique()}; + ctx_type &ctx{*owner}; + std::unique_ptr> point; + + explicit strain_rig(std::string const &method = "newton") { + param_type p; + p.insert("name", "strain"); + ctx.create>(p); + make_system(ctx, {{"g", nm::unknown_kind::scalar}, {"B", nm::unknown_kind::sym_tensor}}, {{"sys::rg", "g"}, {"sys::rB", "B"}}, + method); + p.clear(); + p.insert("name", "sys"); + ctx.create(p); + ctx.finalize(); + param_type q; + q.insert>("inputs", {"strain"}); + point = std::make_unique>(ctx, q); + } + + void evaluate(tensor2 const &eps) { + point->reset(); + point->input>("strain").bind(tensor2{}, eps); + point->evaluate(); + } +}; + +tensor2 const eps0{0.03, 0.01, 0.00, 0.01, 0.05, 0.02, 0.00, 0.02, 0.07}; + +TEST(ResidualSystemTangent, MatchesTheAnalyticStrainDerivative) { + strain_rig rig; + auto const &dg{rig.point->derivative("solver::g", "strain")}; + auto const &dB{rig.point->derivative("solver::B", "strain")}; + rig.evaluate(eps0); + T const denom{1 + tmech::trace(strain_rows::P())}; + EXPECT_NEAR(rig.ctx.get("solver", "g"), (2 + tmech::trace(eps0)) / denom, 1e-12); + EXPECT_TRUE(tmech::almost_equal(dg, tensor2{tmech::eye() / denom}, 1e-12)); + tensor4 const dB_ref{tmech::otimes(strain_rows::P(), tmech::eye()) / denom}; + EXPECT_TRUE(tmech::almost_equal(dB, dB_ref, 1e-12)) << "minor-symmetric, major-asymmetric"; + ASSERT_GT(std::abs(dB_ref(0, 0, 1, 1) - dB_ref(1, 1, 0, 0)), 1e-3); +} + +TEST(ResidualSystemTangent, AMemberReadOutsideGetsTheTotalDerivative) { + strain_rig rig; + auto const &drg{rig.point->derivative("sys::rg", "strain")}; + rig.evaluate(eps0); + EXPECT_LT(tmech::norm(drg), 1e-12) << "R stays zero along the solution: dR/deps|x + dR/dx dx/deps = 0"; +} + +TEST(ResidualSystemTangent, AMissingPartialOnTheNewtonPathThrowsAtFinalize) { + ctx_type ctx; + param_type p; + p.insert("name", "load"); + ctx.create>(p); + make_system(ctx, {{"x", nm::unknown_kind::scalar}}, {{"f::r", "x"}}); + p.clear(); + p.insert("name", "f"); + p.insert("declares_partials", false); + ctx.create(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_TRUE(message.contains("'f::r'")) << message; + EXPECT_TRUE(message.contains("solver::x")) << message; +} + +TEST(ResidualSystemTangent, NewtonFdNeedsNoPartials) { + ctx_type ctx; + param_type p; + p.insert("name", "load"); + auto &load{ctx.create>(p)}; + make_system(ctx, {{"x", nm::unknown_kind::scalar}}, {{"f::r", "x"}}, "newton_fd"); + p.clear(); + p.insert("name", "f"); + p.insert("declares_partials", false); + ctx.create(p); + ctx.finalize(); + load.bind(0, 2); + ctx.update(); + EXPECT_NEAR(ctx.get("solver", "x"), 1, 1e-10); +} + +// --- picard --- + +TEST(ResidualSystem, PicardWithUnknownsIteratesXMinusR) { + ctx_type ctx; + param_type p; + p.insert("name", "load"); + auto &load{ctx.create>(p)}; + auto &picard{make_system(ctx, {{"x", nm::unknown_kind::scalar}}, {{"f::r", "x"}}, "picard", 200, 1e-12)}; + p.clear(); + p.insert("name", "f"); + p.insert("picard", true); + ctx.create(p); + ctx.finalize(); + load.bind(0, 0.3); + ctx.update(); + T const x{ctx.get("solver", "x")}; + EXPECT_NEAR(x, std::cos(0.3 + x) / 2, 1e-12); + EXPECT_GT(picard.iterations(), 5) << "linear convergence: a contraction, not a Newton step"; +} + +// --- validation --- + +TEST(ResidualSystemValidation, RejectsInconsistentUnknownsAndEquations) { + using unknowns = std::vector; + using equations = std::vector; + struct bad { + unknowns u; + equations e; + std::string method; + char const *says; + }; + for (auto const &[u, e, method, says] : { + bad{{{"x", nm::unknown_kind::scalar}, {"x", nm::unknown_kind::scalar}}, {{"a::r", "x"}, {"a::s", "x"}}, "newton", "'x'"}, + bad{{{"", nm::unknown_kind::scalar}}, {{"a::r", ""}}, "newton", "empty"}, + bad{{}, {}, "newton", "unknowns"}, + bad{xy(), {{"a::r", "x"}}, "newton", "'y'"}, + bad{xy(), {{"a::r", "x"}, {"a::s", "z"}}, "newton", "'z'"}, + bad{xy(), {{"a::r", "x"}, {"a::s", "x"}}, "newton", "'x'"}, + bad{xy(), {{"a::r", "x"}, {"a::r", "y"}}, "newton", "'a::r'"}, + bad{xy(), {{"a::r", "x"}, {"a_s", "y"}}, "newton", "a_s"}, + bad{xy(), {{"a::r", "x"}, {"a::s", "y"}}, "gauss_seidel", "gauss_seidel"}, + bad{{}, {{"a::r", "x"}}, "picard", "'x'"}, + bad{{{"B", nm::unknown_kind::sym_tensor, 1.0}}, {{"a::r", "B"}}, "newton", "seed"}, + }) { + ctx_type ctx; + auto const message{thrown([&] { make_system(ctx, u, e, method); })}; + EXPECT_TRUE(message.contains(says)) << says << ": " << message; + EXPECT_TRUE(message.contains("'solver'")) << message; + } +} + +TEST(ResidualSystemValidation, AResidualOfAnotherKindFailsToWire) { + ctx_type ctx; + make_system(ctx, {{"B", nm::unknown_kind::sym_tensor}}, {{"load::state", "B"}}); + param_type p; + p.insert("name", "load"); + ctx.create>(p); + EXPECT_THROW(ctx.finalize(), std::runtime_error) << "a scalar residual for a tensor unknown fails to wire"; +} + +TEST(ResidualSystemValidation, AStatefulMemberThrowsAtFinalize) { + ctx_type ctx; + param_type p; + p.insert("name", "load"); + ctx.create>(p); + make_system(ctx, {{"x", nm::unknown_kind::scalar}}, {{"f::r", "x"}}); + p.clear(); + p.insert("name", "f"); + p.insert("counts", true); + ctx.create(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_TRUE(message.contains("f::r")) << message; + EXPECT_TRUE(message.contains("stateful")) << message; +} + +TEST(ResidualSystemValidation, SystemsAreDisjoint) { + ctx_type ctx; // f::r reads the unknowns of both + make_system(ctx, {{"x", nm::unknown_kind::scalar}}, {{"f::r", "x"}}); + make_system(ctx, {{"state", nm::unknown_kind::scalar}}, {{"f::r", "state"}}, "newton", 50, 1e-10, "load"); + param_type p; + p.insert("name", "f"); + ctx.create(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_TRUE(message.contains("'solver'") && message.contains("'load'")) << message; +} + +TEST(ResidualSystemValidation, NoPathFromTheUnknownsToTheResidualsThrowsAtFinalize) { + ctx_type ctx; + make_system(ctx, {{"x", nm::unknown_kind::scalar}}, {{"load::state", "x"}}); + param_type p; + p.insert("name", "load"); + ctx.create>(p); + auto const message{thrown([&] { ctx.finalize(); })}; + EXPECT_TRUE(message.contains("'solver'")) << message; +} + +// --- "seed": "history" --- + +/// x of every step of `loads` at a host-owned point, its history in `buffer`, and the iterations each took. +struct history_rig { + std::unique_ptr owner{std::make_unique()}; + ctx_type &ctx{*owner}; + system_type *system{nullptr}; + std::unique_ptr> point; + + explicit history_rig(bool history) { + param_type p; + p.insert("name", "load"); + ctx.create>(p); + system = &make_system(ctx, {{"x", nm::unknown_kind::scalar, 0, history}}, {{"f::r", "x"}}, "newton", 50, 1e-12); + p.clear(); + p.insert("name", "f"); + ctx.create(p); + ctx.finalize(); + param_type q; + q.insert>("inputs", {"load"}); + point = std::make_unique>(ctx, q); + } + + std::pair step(std::vector &buffer, T load) { + buffer.resize(point->history_width()); + point->reset(); + point->load_history(buffer); + point->input>("load").bind(0, load); + point->evaluate(); + point->store_history(buffer); + return {ctx.get("solver", "x"), system->iterations()}; + } +}; + +TEST(ResidualSystemHistorySeed, TwoPointsAlternatelyTakeTheIterationsOfEachAlone) { + std::vector const loads_a{1, 2, 3, 3.5}, loads_b{-2, -1, 0.5, 4}; + auto const alone{[](std::vector const &loads, bool history) { + history_rig rig{history}; + std::vector buffer; + std::vector> out; + for (auto const load : loads) + out.push_back(rig.step(buffer, load)); + return out; + }}; + auto const a{alone(loads_a, true)}, b{alone(loads_b, true)}; + + history_rig shared{true}; + EXPECT_EQ(shared.point->history_width(), 1u) << "the unknown is kept as history"; + std::vector buffer_a, buffer_b; + for (std::size_t i{0}; i < loads_a.size(); ++i) { + EXPECT_EQ(shared.step(buffer_a, loads_a[i]), a[i]) << "point a, step " << i; + EXPECT_EQ(shared.step(buffer_b, loads_b[i]), b[i]) << "point b, step " << i; + } + + int warm{0}, cold{0}; + for (auto const &[x, iterations] : alone(loads_a, true)) + warm += iterations; + for (auto const &[x, iterations] : alone(loads_a, false)) + cold += iterations; + EXPECT_LT(warm, cold) << "each step starts from the last converged x"; +} + +// --- review of the residual systems --- + +/// A "solver" over the scalar x (seed `seed`) whose residual "f::r" is the row `Row`, and a "load" source. +template +struct row_rig { + std::unique_ptr owner{std::make_unique()}; + ctx_type &ctx{*owner}; + system_type *system{nullptr}; + nm::external_scalar_source *load{nullptr}; + + row_rig(std::string const &method, int max_iter, T tolerance, T seed, param_type row = {}) { + param_type p; + p.insert("name", "load"); + load = &ctx.create>(p); + system = &make_system(ctx, {{"x", nm::unknown_kind::scalar, static_cast(seed)}}, {{"f::r", "x"}}, method, max_iter, tolerance); + row.insert("name", "f"); + ctx.create(row); + ctx.finalize(); + } + + T x() { return ctx.get("solver", "x"); } +}; + +TEST(ResidualSystemReview, TheIterationsAreReportedWhenTheSolveThrows) { + auto const ref{vector_newton_on(system_mode::nonlinear, 3, 5)}; + ASSERT_GT(ref.iterations, 1); + scalar_rig rig{system_mode::nonlinear, 3, 5, ref.iterations - 1}; + EXPECT_THROW(rig.ctx.update(), nm::cutback_error); + EXPECT_EQ(rig.system->iterations(), ref.iterations - 1) << "the updates taken before the throw"; +} + +TEST(ResidualSystemReview, ZeroBlocksBeforeFinalizeThrows) { + ctx_type ctx; + auto &system{make_system(ctx, xy(), {{"fx::rx", "x"}, {"fy::ry", "y"}})}; + auto const message{thrown([&] { (void)system.zero_blocks(); })}; + EXPECT_TRUE(message.contains("'solver'") && message.contains("finalize")) << message; +} + +TEST(ResidualSystemReview, NewtonFdLeavesTheUnknownsAtTheIterateWhenTheJacobianIsSingular) { + scalar_rig rig{system_mode::singular, 4, 7, 50, false, "newton_fd"}; + EXPECT_THROW(rig.ctx.update(), nm::cutback_error); + EXPECT_EQ(rig.x(), 0) << "not x + h"; + EXPECT_EQ(rig.y(), 0); +} + +TEST(ResidualSystemReview, ANonFiniteResidualNeverConverges) { + T const inf{std::numeric_limits::infinity()}; + for (T const value : {std::numeric_limits::quiet_NaN(), inf, -inf}) + for (T const tolerance : {T{1e-10}, inf}) + for (std::string const method : {"picard", "newton_fd"}) { + param_type row; + row.insert("value", value); + row_rig rig{method, 5, tolerance, 0, row}; + auto const message{thrown([&] { rig.ctx.update(); })}; + EXPECT_TRUE(message.contains("'solver'")) << value << ", tolerance " << tolerance << ", " << method << ": " << message; + } +} + +TEST(ResidualSystemReview, AConvergedSeedNeedsNoUpdateWithToleranceAndMaxIterZero) { + for (std::string const method : {"newton", "newton_fd", "picard"}) { + row_rig rig{method, 0, 0, 1}; + rig.load->bind(1, 1); // x = 1 solves x^2 = 1 exactly + EXPECT_NO_THROW(rig.ctx.update()) << method; + EXPECT_EQ(rig.system->iterations(), 0) << method; + EXPECT_EQ(rig.x(), 1); + row_rig off{method, 0, 0, 1}; + off.load->bind(1, 2); + EXPECT_THROW(off.ctx.update(), nm::cutback_error) << method << ": not converged, and no update allowed"; + } +} + +TEST(ResidualSystemReview, NewtonFdConvergesAsFastAsNewtonOnSmallUnknowns) { + for (T const scale : {T{1e-9}, T{1}, T{1e6}}) { + param_type row; + row.insert("scale", scale); + row_rig newton{"newton", 50, 1e-10, 0.5 * scale, row}; + newton.ctx.update(); + EXPECT_NEAR(newton.x() / scale, 1, 1e-9); + row_rig fd{"newton_fd", 50, 1e-10, 0.5 * scale, row}; + EXPECT_NO_THROW(fd.ctx.update()) << "scale " << scale; + EXPECT_NEAR(fd.x() / scale, 1, 1e-9) << "scale " << scale; + EXPECT_LE(fd.system->iterations(), newton.system->iterations() + 1) << "scale " << scale; + } +} + +TEST(ResidualSystemTangent, NewtonFdAndPicardGiveTheImplicitFunctionTangent) { + row_rig fd{"newton_fd", 50, 1e-12, 0}; + auto const &dx_fd{fd.ctx.derivative("solver", "x", "load", "state")}; + fd.load->bind(0, 2); + fd.ctx.update(); + EXPECT_NEAR(fd.x(), 1, 1e-10); + EXPECT_NEAR(dx_fd, 1 / (3 * fd.x() * fd.x() + 1), 1e-10); + + param_type row; + row.insert("picard", true); + row_rig picard{"picard", 500, 1e-13, 0, row}; + auto const &dx{picard.ctx.derivative("solver", "x", "load", "state")}; + auto const at{[&](T load) { + picard.load->bind(0, load); + picard.ctx.update(); + return picard.x(); + }}; + T const load{0.3}, h{1e-5}; + T const plus{at(load + h)}, minus{at(load - h)}; + T const x{at(load)}; + T const s{std::sin(load + x)}; + EXPECT_NEAR(dx, -(s / 2) / (1 + s / 2), 1e-10); + EXPECT_NEAR(dx, (plus - minus) / (2 * h), 1e-7); +} + +TEST(ResidualSystemTangent, TwoCoupledSystemsMatchCentralDifferences) { + ctx_type ctx; + param_type p; + p.insert("name", "load"); + auto &load{ctx.create>(p)}; + make_system(ctx, {{"x", nm::unknown_kind::scalar}}, {{"f::r", "x"}}, "newton", 50, 1e-13); + make_system(ctx, {{"y", nm::unknown_kind::scalar}}, {{"g::r", "y"}}, "newton", 50, 1e-13, "second"); + p.clear(); + p.insert("name", "f"); + ctx.create(p); + p.clear(); + p.insert("name", "g"); + ctx.create(p); + ctx.finalize(); + auto const &dx{ctx.derivative("solver", "x", "load", "state")}; + auto const &dy{ctx.derivative("second", "y", "load", "state")}; + auto const at{[&](T value) { + load.bind(0, value); + ctx.update(); + return std::pair{ctx.get("solver", "x"), ctx.get("second", "y")}; + }}; + T const value{3}, h{1e-5}; + auto const [xp, yp]{at(value + h)}; + auto const [xm, ym]{at(value - h)}; + auto const [x, y]{at(value)}; + EXPECT_NEAR(dx, (xp - xm) / (2 * h), 1e-7); + EXPECT_NEAR(dy, (yp - ym) / (2 * h), 1e-7) << "the second system reads the unknown of the first"; + EXPECT_NEAR(dy, 1 / (3 * y * y + 1) / (3 * x * x + 1), 1e-10); +} + +TEST(ResidualSystemTangent, ASingularJacobianAtTheConvergedPointThrowsACutbackNamingTheSystem) { + row_rig rig{"newton", 50, 1e-10, 0}; + rig.ctx.derivative("solver", "x", "load", "state"); + rig.load->bind(0, 0); // x = 0 is a double root, and the seed already solves it + auto const message{thrown([&] { rig.ctx.update(); })}; + EXPECT_TRUE(message.contains("'solver'")) << message; +} + +TEST(ResidualSystemTangent, ASeedThatIsAMemberOutputOfAPassedSystemThrows) { + strain_rig rig; + auto const message{thrown([&] { rig.ctx.derivative("solver", "g", "sys", "rg"); })}; + EXPECT_TRUE(message.contains("sys::rg") && message.contains("'solver'")) << message; +} + +#ifdef NUMSIM_HAVE_JSON +TEST(ResidualSystemJson, ADocumentBuildsTheSystem) { + nm::register_default_materials(); + nm::material_factory::instance().template register_type("scalar_rows"); + auto const doc = nlohmann::json::parse(R"({"materials": [ + {"type": "system", "name": "solver", "method": "newton", "tolerance": 1e-12, "max_iter": 25, + "unknowns": [{"name": "x", "kind": "scalar", "seed": 0}, {"name": "y", "kind": "scalar", "seed": "history"}], + "equations": [{"residual": "sys::rx", "for": "x"}, {"residual": "sys::ry", "for": "y"}]}, + {"type": "scalar_rows", "name": "sys", "mode": 0, "a": 5.0, "b": 10.0}]})"); + ctx_type ctx; + for (auto const &m : doc["materials"]) + nm::create_from_json(ctx, m); + ctx.finalize(); + ctx.update(); + EXPECT_NEAR(ctx.get("solver", "x"), 1, 1e-12); + EXPECT_NEAR(ctx.get("solver", "y"), 3, 1e-12); + EXPECT_TRUE(ctx.find_property("solver", "y")->is_history()); + EXPECT_FALSE(ctx.find_property("solver", "x")->is_history()); +} + +TEST(ResidualSystemJson, RejectsAnUnknownSeedOrKind) { + nm::register_default_materials(); + for (auto const *unknown : {R"({"name": "x", "kind": "scalar", "seed": "warm"})", R"({"name": "x", "kind": "vector"})"}) { + auto const doc = nlohmann::json::parse(std::string{R"({"type": "system", "name": "solver", "unknowns": [)"} + unknown + + R"(], "equations": [{"residual": "sys::rx", "for": "x"}]})"); + ctx_type ctx; + auto const message{thrown([&] { nm::create_from_json(ctx, doc); })}; + EXPECT_TRUE(message.contains("unknowns") && message.contains("'x'")) << unknown << ": " << message; + } +} + +TEST(ResidualSystemJson, TheReadersRequireListsAndRejectUnknownKeys) { + nm::register_default_materials(); + struct bad { + char const *unknowns, *equations, *says; + }; + for (auto const &[unknowns, equations, says] : { + bad{R"({"name": "x", "kind": "scalar"})", R"([{"residual": "sys::rx", "for": "x"}])", "unknowns"}, + bad{R"(["x"])", R"([{"residual": "sys::rx", "for": "x"}])", "unknowns"}, + bad{R"([{"name": "x", "kind": "scalar", "seeed": 1}])", R"([{"residual": "sys::rx", "for": "x"}])", "seeed"}, + bad{R"([{"name": "x", "kind": "scalar"}])", R"({"residual": "sys::rx", "for": "x"})", "equations"}, + bad{R"([{"name": "x", "kind": "scalar"}])", R"([{"residual": "sys::rx", "forr": "x"}])", "forr"}, + }) { + auto const doc = nlohmann::json::parse(std::string{R"({"type": "system", "name": "solver", "unknowns": )"} + unknowns + + R"(, "equations": )" + equations + "}"); + ctx_type ctx; + auto const message{thrown([&] { nm::create_from_json(ctx, doc); })}; + EXPECT_TRUE(message.contains(says)) << unknowns << " / " << equations << ": " << message; + } +} +#endif + +} // namespace diff --git a/tests/test_system_plasticity.cpp b/tests/test_system_plasticity.cpp new file mode 100644 index 0000000..fce1dea --- /dev/null +++ b/tests/test_system_plasticity.cpp @@ -0,0 +1,317 @@ +// j2_kinematic_plasticity (#166) as a residual system (#218): a flow node holds the residuals and the history, +// a stress node reads the plastic strain, and the system solves {dlambda, dbeta}. The tangent is the +// derivative the engine accumulates through the system, compared with the shipped material. +#include +#include +#include +#include +#include +#include +#include "numsim-materials/core/material_base.h" +#include "numsim-materials/core/material_context.h" +#include "numsim-materials/core/parameter_helpers.h" +#include "numsim-materials/core/partial_check.h" +#include "numsim-materials/core/tensor_identities.h" +#include "numsim-materials/core/unknown_spec.h" +#include "numsim-materials/driver/external_state_source.h" +#include "numsim-materials/materials/j2_kinematic_plasticity.h" +#include "numsim-materials/materials/linear_isotropic_hardening.h" +#include "numsim-materials/materials/plasticity_utils.h" +#include "numsim-materials/solvers/system.h" +#include "numsim-materials/solvers/vector_newton.h" + +namespace { + +namespace nm = numsim::materials; +using policy = nm::material_policy_default; +using T = policy::value_type; +using ctx_type = nm::material_context; +using param_type = policy::ParameterHandler; +using t2 = tmech::tensor; +using t4 = tmech::tensor; + +constexpr T K{166.67}, G{76.92}, sigma_0{1.0}, H_iso{5.0}, C{30.0}; + +T equivalent(t2 const &x) { return std::sqrt(T{1.5} * tmech::dcontract(x, x)); } + +/// Yield and Armstrong-Frederick residuals of j2_kinematic_plasticity at the unknowns of "plastic", and the +/// history they imply. An elastic trial state turns the yield row into dlambda = 0. +class kinematic_flow final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit kinematic_flow(Args &&...args) + : base(std::forward(args)...), m_r_l(base::add_output("r_l")), m_r_b(base::add_output("r_b")), + m_eps_p(base::add_history_output("plastic_strain")), m_kappa(base::add_history_output("equivalent_plastic_strain")), + m_beta(base::add_history_output("back_stress")), m_gamma(base::get_parameter("gamma")), + m_strain(base::add_input("stepper", "strain")), m_dl(base::add_input("plastic", "dlambda")), + m_db(base::add_input("plastic", "dbeta")) { + base::add_callback(&kinematic_flow::compute, m_r_l, m_r_b, m_eps_p, m_kappa, m_beta) + .reads(m_strain, m_dl, m_db) + .partials(m_r_l, m_strain, m_dl, m_db) + .partials(m_r_b, m_strain, m_dl, m_db) + .partials(m_eps_p, m_strain, m_dl, m_db); + } + + static input_parameter_controller parameters() { + input_parameter_controller para{base::parameters()}; + nm::optional(para, "gamma", T{0}); + return para; + } + + void compute() { + auto const &IIdev{nm::dev_projector()}; + auto const dl{m_dl.get()}; + auto const kappa_n{m_kappa.old_value()}; + t2 const beta_n{m_beta.old_value()}; + t2 const s_tr{2 * G * tmech::dev(t2{m_strain.get() - m_eps_p.old_value()})}; + t2 const beta{beta_n + m_db.get()}; + t2 const eta{s_tr - tmech::dev(beta)}; + auto const q{equivalent(eta)}; + bool const has_direction{q > 1e-10 * sigma_0}; + t2 const N{has_direction ? t2{T{1.5} * eta / q} : t2{}}; + t4 const A{has_direction ? t4{(T{1.5} * IIdev - tmech::otimes(N, N)) / q} : t4{}}; + bool const plastic{equivalent(t2{s_tr - beta_n}) - sigma_0 - H_iso * kappa_n > 0}; + auto const c{T{2} / 3 * C}; + + m_r_l = plastic ? q - 3 * G * dl - sigma_0 - H_iso * (kappa_n + dl) : dl; + m_r_b = t2{(1 + m_gamma * dl) * beta - beta_n - (c * dl) * N}; + m_eps_p.new_value() = t2{m_eps_p.old_value() + dl * N}; + m_kappa.new_value() = kappa_n + dl; + m_beta.new_value() = beta; + if (!base::want_partials()) + return; + base::partial(m_r_l, m_strain) = plastic ? t2{2 * G * N} : t2{}; + base::partial(m_r_l, m_dl) = plastic ? -3 * G - H_iso : T{1}; + base::partial(m_r_l, m_db) = plastic ? t2{-N} : t2{}; + base::partial(m_r_b, m_strain) = t4{-(2 * G * c * dl) * A}; + base::partial(m_r_b, m_dl) = t2{m_gamma * beta - c * N}; + base::partial(m_r_b, m_db) = t4{(1 + m_gamma * dl) * nm::sym_identity() + (c * dl) * A}; + base::partial(m_eps_p, m_strain) = t4{(2 * G * dl) * A}; + base::partial(m_eps_p, m_dl) = N; + base::partial(m_eps_p, m_db) = t4{-dl * A}; + } + +private: + T &m_r_l; + t2 &m_r_b; + nm::history_property &m_eps_p; + nm::history_property &m_kappa; + nm::history_property &m_beta; + T const &m_gamma; + nm::input_property const &m_strain; + nm::input_property const &m_dl; + nm::input_property const &m_db; +}; + +/// sigma = C_e : (eps - flow::plastic_strain). +class elastic_stress final : public nm::material_base { +public: + using base = nm::material_base; + + template + explicit elastic_stress(Args &&...args) + : base(std::forward(args)...), m_stress(base::add_output("stress")), m_strain(base::add_input("stepper", "strain")), + m_eps_p(base::add_input("flow", "plastic_strain")) { + base::add_callback(&elastic_stress::compute, m_stress).reads(m_strain, m_eps_p).partials(m_stress, m_strain, m_eps_p); + } + + static input_parameter_controller parameters() { return input_parameter_controller{base::parameters()}; } + + void compute() { + m_stress = tmech::dcontract(m_C, t2{m_strain.get() - m_eps_p.get()}); + if (base::want_partials()) { + base::partial(m_stress, m_strain) = m_C; + base::partial(m_stress, m_eps_p) = t4{-m_C}; + } + } + +private: + t2 &m_stress; + nm::input_property const &m_strain; + nm::input_property const &m_eps_p; + t4 const m_C{nm::plasticity_detail::make_isotropic_tangent(K, G)}; +}; + +/// The shipped material on vector_newton, as test_j2_kinematic_plasticity builds it. +struct shipped { + ctx_type ctx; + nm::external_strain_source *src{}; + + explicit shipped(T gamma) { + param_type p; + p.insert("name", "stepper"); + src = &ctx.create>(p); + p.clear(); + p.insert("name", "solver"); + p.insert("function", "kin"); + p.insert("tolerance", T{1e-12}); + p.insert("strain_derivative", true); + p.insert("graph_driven", false); + p.insert>("unknowns", {{"dlambda", nm::unknown_kind::scalar}, {"dbeta", nm::unknown_kind::sym_tensor}}); + ctx.create>(p); + p.clear(); + p.insert("name", "hardening"); + p.insert("source", "kin"); + p.insert("K", H_iso); + ctx.create>(p); + p.clear(); + p.insert("name", "kin"); + p.insert("hardening_source", "hardening"); + p.insert("strain_source", "stepper"); + p.insert("solver_source", "solver"); + p.insert("K", K); + p.insert("G", G); + p.insert("sigma_0", sigma_0); + p.insert("C", C); + p.insert("gamma", gamma); + ctx.create>(p); + ctx.finalize(); + } + + void step(t2 const &eps) { + src->bind(eps, eps); + ctx.update(); + } + t2 const &stress() { return ctx.get("kin", "stress"); } + t4 const &tangent() { return ctx.get("kin", "tangent"); } + t2 const &back_stress() { return ctx.get("kin", "back_stress"); } +}; + +/// The same model as a system over a flow node and a stress node. +struct as_system { + ctx_type ctx; + nm::external_strain_source *src{}; + nm::residual_system *system{}; + t4 const *tangent{}; + + explicit as_system(T gamma) { + param_type p; + p.insert("name", "stepper"); + src = &ctx.create>(p); + p.clear(); + p.insert("name", "plastic"); + p.insert("tolerance", T{1e-12}); + p.insert>("unknowns", {{"dlambda", nm::unknown_kind::scalar}, {"dbeta", nm::unknown_kind::sym_tensor}}); + p.insert>("equations", {{"flow::r_l", "dlambda"}, {"flow::r_b", "dbeta"}}); + system = &ctx.create>(p); + p.clear(); + p.insert("name", "flow"); + p.insert("gamma", gamma); + ctx.create(p); + p.clear(); + p.insert("name", "kin"); + ctx.create(p); + ctx.finalize(); + tangent = &ctx.derivative("kin", "stress", "stepper", "strain"); + } + + void step(t2 const &eps) { + src->bind(eps, eps); + ctx.update(); + } + t2 const &stress() { return ctx.get("kin", "stress"); } + t2 const &back_stress() { return ctx.get("flow", "back_stress"); } +}; + +t2 strain_a() { + t2 e; + e.fill(T{0}); + e(0, 0) = 0.010; + e(1, 1) = -0.003; + e(0, 1) = e(1, 0) = 0.004; + return e; +} +t2 strain_b() { + t2 e; + e.fill(T{0}); + e(0, 0) = -0.002; + e(2, 2) = 0.006; + e(0, 1) = e(1, 0) = 0.012; + return e; +} + +template +T relative(A const &a, A const &b) { + return tmech::norm(A{a - b}) / tmech::norm(b); +} + +TEST(SystemPlasticity, MatchesTheShippedStressAndTangent) { + for (T const gamma : {0.0, 40.0}) { + shipped ref{gamma}; + as_system sys{gamma}; + for (auto const &eps : {t2{0.1 * strain_a()}, strain_a(), strain_b(), t2{0.9 * strain_b()}}) { + ref.step(eps); + sys.step(eps); + EXPECT_LT(relative(sys.stress(), ref.stress()), 1e-12) << "gamma " << gamma; + EXPECT_LT(relative(*sys.tangent, ref.tangent()), 1e-12) << "gamma " << gamma; + if (tmech::norm(ref.back_stress()) > 0) { + EXPECT_LT(relative(sys.back_stress(), ref.back_stress()), 1e-12); + } + ref.ctx.commit(); + sys.ctx.commit(); + } + EXPECT_GT(tmech::norm(sys.back_stress()), 0) << "the path yields"; + } +} + +TEST(SystemPlasticity, OnlyRecoveryNeedsMoreThanOneNewtonStep) { + for (T const gamma : {0.0, 40.0}) { + as_system sys{gamma}; + sys.step(strain_a()); + EXPECT_EQ(sys.system->iterations() == 1, gamma == 0) << "Prager is a radial return: one step from the trial state"; + } +} + +TEST(SystemPlasticity, NoZeroBlocksAndThePartialsMatchFiniteDifferences) { + as_system sys{40.0}; + EXPECT_TRUE(sys.system->zero_blocks().empty()); + sys.step(strain_a()); + sys.ctx.commit(); + sys.step(strain_b()); + for (auto const *name : {"flow", "kin"}) + for (auto const &mismatch : nm::check_partials(sys.ctx, name)) + ADD_FAILURE() << mismatch; +} + +TEST(SystemPlasticity, UpdatePropertySolvesTheSystemLikeUpdate) { + for (T const gamma : {0.0, 40.0}) { + as_system full{gamma}, single{gamma}; + full.step(strain_a()); + single.src->bind(strain_a(), strain_a()); + single.ctx.update_property("kin", "stress"); + EXPECT_GT(tmech::norm(full.back_stress()), 0) << "the path yields"; + EXPECT_LT(relative(single.stress(), full.stress()), 1e-14) << "gamma " << gamma; + EXPECT_LT(relative(single.back_stress(), full.back_stress()), 1e-14) << "gamma " << gamma; + } +} + +TEST(SystemPlasticity, TheTangentMatchesCentralDifferencesOfTheStress) { + for (T const gamma : {0.0, 40.0}) { + as_system sys{gamma}; + sys.step(t2{0.1 * strain_a()}); + sys.ctx.commit(); + sys.step(strain_a()); + ASSERT_GT(tmech::norm(sys.back_stress()), 0) << "the step yields"; + t4 const tangent{*sys.tangent}; + constexpr T h{1e-6}; + for (int i{0}; i < 3; ++i) + for (int j{i}; j < 3; ++j) { + t2 E; + E.fill(T{0}); + E(i, j) = E(j, i) = i == j ? T{1} : T{0.5}; + auto const stress_at{[&](T step) { + sys.ctx.revert(); + sys.step(t2{strain_a() + step * E}); + return t2{sys.stress()}; + }}; + t2 const plus{stress_at(h)}, minus{stress_at(-h)}; + t2 const fd{(T{1} / (2 * h)) * (plus - minus)}; + t2 const analytic{tmech::dcontract(tangent, E)}; + EXPECT_LT(relative(fd, analytic), 1e-5) << "gamma " << gamma << ", component " << i << j; + } + } +} + +} // namespace