From 6a6c8c33c8d5814f3b22fcedfeef8cf85f8ab57b Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 17:38:27 +0900 Subject: [PATCH] perf(hex8): defer optional paired tangent until a valid linear trial --- include/fesa/elements/hex8_element.h | 15 ++- src/fesa/elements/hex8_element.cpp | 12 +- src/fesa/elements/hex8_high_low.cpp | 32 +++-- src/fesa/elements/hex8_linear_tangent_cache.h | 21 ++++ tests/unit/elements/hex8_high_low_test.cpp | 109 ++++++++++++++++++ 5 files changed, 167 insertions(+), 22 deletions(-) create mode 100644 src/fesa/elements/hex8_linear_tangent_cache.h diff --git a/include/fesa/elements/hex8_element.h b/include/fesa/elements/hex8_element.h index 2cab0e2..a9062b7 100644 --- a/include/fesa/elements/hex8_element.h +++ b/include/fesa/elements/hex8_element.h @@ -2,6 +2,7 @@ #define FESA_ELEMENTS_HEX8_ELEMENT_H_ #include +#include #include #include "fesa/elements/element.h" @@ -9,6 +10,9 @@ #include "fesa/results/solid_result_records.h" namespace fesa { +namespace hex8_internal { +struct LinearTangentCache; +} /// @brief Owns all local quantities from one validated linear or TL trial. struct Hex8Response { @@ -90,7 +94,8 @@ class Hex8Element final : public Element { bool nonlinear) const; /// @brief Uses one paired gradient for energy, force and analytic Hessian. Result EvaluateHighLow( - const Vector& high, const Vector& low, bool nonlinear) const; + const Vector& high, const Vector& low, bool nonlinear, + bool prepare_linear_cache = true) const; /// @brief Attaches stable source context to a model failure. Status Failure(const char* message) const; @@ -103,11 +108,9 @@ class Hex8Element final : public Element { SourceLocation location_; std::array points_{}; Matrix stiffness_{24U, 24U}; - // Optional immutable Hessian from the same paired quadrature graph. Failure - // to build it does not reject models supported by the original F64 paths. - bool paired_linear_tangent_ready_{false}; - Matrix paired_linear_tangent_high_{24U, 24U}; - Matrix paired_linear_tangent_low_{24U, 24U}; + // Optional immutable paired Hessian, prepared once on the first valid paired + // linear trial. Copies share only their identical owned reference definition. + std::shared_ptr paired_linear_tangent_; }; } // namespace fesa diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index 26a9b78..f2bb5d2 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -7,6 +7,7 @@ #include "../math/dense_blas_internal.h" #include "hex8_geometry.h" +#include "hex8_linear_tangent_cache.h" namespace fesa { namespace { @@ -116,7 +117,9 @@ Hex8Element::Hex8Element(EntityIndex index, ElementDofLayout layout, kappa_{kappa}, mu_{mu}, density_{density}, - location_{std::move(location)} {} + location_{std::move(location)}, + paired_linear_tangent_{ + std::make_shared()} {} Status Hex8Element::Failure(const char* message) const { return Status::Failure(FailureCategory::kModel, @@ -273,13 +276,6 @@ Result Hex8Element::Create(EntityIndex index, if (backend_inputs_failed) return Result::Failure(candidate.Failure( "Hex8 dense MKL operands failed finite validation.")); - const Vector zero{24U}; - const auto paired = candidate.EvaluateHighLow(zero, zero, false); - if (paired.HasValue()) { - candidate.paired_linear_tangent_high_ = paired.Value().high.tangent; - candidate.paired_linear_tangent_low_ = paired.Value().tangent_low; - candidate.paired_linear_tangent_ready_ = true; - } return Result::Success(std::move(candidate)); } diff --git a/src/fesa/elements/hex8_high_low.cpp b/src/fesa/elements/hex8_high_low.cpp index 26c6b8b..4d377e7 100644 --- a/src/fesa/elements/hex8_high_low.cpp +++ b/src/fesa/elements/hex8_high_low.cpp @@ -2,10 +2,12 @@ #include #include +#include #include #include "../math/high_low_arithmetic.h" #include "hex8_geometry.h" +#include "hex8_linear_tangent_cache.h" #include "hex8_pair_arithmetic.h" namespace fesa { @@ -53,7 +55,8 @@ std::array Pack(const Tensor& a) { } // namespace Result Hex8Element::EvaluateHighLow( - const Vector& high, const Vector& low, bool nonlinear) const { + const Vector& high, const Vector& low, bool nonlinear, + bool prepare_linear_cache) const { using math_internal::IsNormalized; if (!math_internal::ArithmeticEnvironmentSupported()) return Result::Failure(Failure( @@ -85,6 +88,19 @@ Result Hex8Element::EvaluateHighLow( message = "Hex8 full high+low Jacobian interval arithmetic is nonfinite."; return Result::Failure(Failure(message)); } + if (!nonlinear && prepare_linear_cache) { + // Trial guards above run before optional reference preparation. The nested + // zero trial bypasses this once entry and uses the original paired graph. + std::call_once(paired_linear_tangent_->preparation, [this] { + const Vector zero{24U}; + const auto paired = EvaluateHighLow(zero, zero, false, false); + if (paired.HasValue()) { + paired_linear_tangent_->high = paired.Value().high.tangent; + paired_linear_tangent_->low = paired.Value().tangent_low; + paired_linear_tangent_->ready = true; + } + }); + } Hex8HighLowResponse response; std::array force{}; @@ -196,7 +212,7 @@ Result Hex8Element::EvaluateHighLow( for (std::size_t i = 0; i < 3; ++i) { force[3U * a + i] += Dot(p[i], gradients[a]) * HighLow{point.weight}; } - if (!nonlinear && paired_linear_tangent_ready_) continue; + if (!nonlinear && paired_linear_tangent_->ready) continue; std::array fg{}; for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) @@ -243,7 +259,7 @@ Result Hex8Element::EvaluateHighLow( } } } - if (nonlinear || !paired_linear_tangent_ready_) { + if (nonlinear || !paired_linear_tangent_->ready) { std::array bulk_fg{}; for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) @@ -260,10 +276,10 @@ Result Hex8Element::EvaluateHighLow( Failure("Hex8 high+low energy is unsupported or nonfinite.")); response.high.strain_energy = energy.high; response.strain_energy_low = energy.low; - if (!nonlinear && paired_linear_tangent_ready_) { - response.high.material_tangent = paired_linear_tangent_high_; - response.high.tangent = paired_linear_tangent_high_; - response.tangent_low = paired_linear_tangent_low_; + if (!nonlinear && paired_linear_tangent_->ready) { + response.high.material_tangent = paired_linear_tangent_->high; + response.high.tangent = paired_linear_tangent_->high; + response.tangent_low = paired_linear_tangent_->low; } for (std::size_t i = 0; i < 24; ++i) { if (!IsNormalized(force[i])) @@ -271,7 +287,7 @@ Result Hex8Element::EvaluateHighLow( Failure("Hex8 high+low force is unsupported or nonfinite.")); response.high.internal_force[i] = force[i].high; response.internal_force_low[i] = force[i].low; - if (!nonlinear && paired_linear_tangent_ready_) continue; + if (!nonlinear && paired_linear_tangent_->ready) continue; for (std::size_t j = 0; j < 24; ++j) { const auto tangent = material[i][j] + geometric[i][j]; if (!IsNormalized(material[i][j]) || !IsNormalized(geometric[i][j]) || diff --git a/src/fesa/elements/hex8_linear_tangent_cache.h b/src/fesa/elements/hex8_linear_tangent_cache.h new file mode 100644 index 0000000..6e5087c --- /dev/null +++ b/src/fesa/elements/hex8_linear_tangent_cache.h @@ -0,0 +1,21 @@ +#ifndef FESA_SRC_FESA_ELEMENTS_HEX8_LINEAR_TANGENT_CACHE_H_ +#define FESA_SRC_FESA_ELEMENTS_HEX8_LINEAR_TANGENT_CACHE_H_ + +#include + +#include "fesa/math/matrix.h" + +namespace fesa::hex8_internal { +/// @brief Shares one immutable paired Hessian among identical element copies. +/// @note Only call_once preparation writes. Linear readers synchronize through +/// the flag; nonlinear trials do not access this optional state. At most two +/// 24-by-24 matrices (9216 value bytes) are retained after successful preparation. +struct LinearTangentCache { + std::once_flag preparation; + bool ready{false}; + Matrix high{0U, 0U}; + Matrix low{0U, 0U}; +}; +} // namespace fesa::hex8_internal + +#endif // FESA_SRC_FESA_ELEMENTS_HEX8_LINEAR_TANGENT_CACHE_H_ diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index 6ae57cc..68deaa1 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -3,10 +3,13 @@ #include #include +#include #include #include +#include #include #include +#include #include #include "fesa/materials/isotropic_linear_elastic_material.h" @@ -93,6 +96,112 @@ fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) { return fesa::Hex8Element::Create(0, layout, coordinates, material, 2).Value(); } +void ExpectSamePairResponse(const fesa::Hex8HighLowResponse& actual, + const fesa::Hex8HighLowResponse& expected) { + EXPECT_EQ(PairBits(actual.high.strain_energy), + PairBits(expected.high.strain_energy)); + EXPECT_EQ(PairBits(actual.strain_energy_low), + PairBits(expected.strain_energy_low)); + ASSERT_EQ(actual.high.tangent.Rows(), 24U); + ASSERT_EQ(actual.tangent_low.Rows(), 24U); + for (std::size_t i = 0; i < 24; ++i) { + EXPECT_EQ(PairBits(actual.high.internal_force[i]), + PairBits(expected.high.internal_force[i])); + EXPECT_EQ(PairBits(actual.internal_force_low[i]), + PairBits(expected.internal_force_low[i])); + for (std::size_t j = 0; j < 24; ++j) { + EXPECT_EQ(PairBits(actual.high.tangent(i, j)), + PairBits(expected.high.tangent(i, j))); + EXPECT_EQ(PairBits(actual.tangent_low(i, j)), + PairBits(expected.tangent_low(i, j))); + } + } + for (std::size_t point = 0; point < 8; ++point) { + EXPECT_EQ(actual.high.points[point].integration_point, + expected.high.points[point].integration_point); + for (std::size_t i = 0; i < 6; ++i) { + EXPECT_EQ(PairBits(actual.high.points[point].cauchy_stress[i]), + PairBits(expected.high.points[point].cauchy_stress[i])); + EXPECT_EQ(PairBits(actual.high.points[point].strain[i]), + PairBits(expected.high.points[point].strain[i])); + } + } +} + +TEST(Hex8HighLow, ConcurrentFirstLinearTrialsAcrossCopiesAndMoveKeepBits) { + const auto element = UnitBrick(); + auto copied = element; + const auto moved = std::move(copied); + fesa::Vector high(24), low_a(24), low_b(24); + for (std::size_t i = 0; i < 24; ++i) high[i] = 1; + for (auto a : {1, 2, 5, 6}) { + low_a[3U * a] = std::ldexp(1.0, -60); + low_b[3U * a] = -std::ldexp(1.0, -59); + } + const auto oracle = UnitBrick(); + const auto expected_a = oracle.EvaluateLinearHighLow(high, low_a); + const auto expected_b = oracle.EvaluateLinearHighLow(high, low_b); + ASSERT_TRUE(expected_a.HasValue()); + ASSERT_TRUE(expected_b.HasValue()); + std::promise start; + const auto signal = start.get_future().share(); + std::array>, 4> calls; + for (std::size_t i = 0; i < calls.size(); ++i) + calls[i] = std::async(std::launch::async, [&, i] { + signal.wait(); + return (i % 2U == 0U ? element : moved) + .EvaluateLinearHighLow(high, i < 2U ? low_a : low_b); + }); + start.set_value(); + for (std::size_t i = 0; i < calls.size(); ++i) { + const auto response = calls[i].get(); + ASSERT_TRUE(response.HasValue()); + ExpectSamePairResponse(response.Value(), + i < 2U ? expected_a.Value() : expected_b.Value()); + } +} + +TEST(Hex8HighLow, FirstLinearPreparationRecoversAfterTrialAndEnvironmentFailure) { + ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported()); + const auto element = UnitBrick(); + fesa::Vector high(24), low(24); + EXPECT_FALSE(element.EvaluateLinearHighLow(fesa::Vector(23), low).HasValue()); + high[0] = std::numeric_limits::infinity(); + EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue()); + high[0] = 0; + struct RestoreEnvironment { + int rounding = std::fegetround(); +#if defined(_M_X64) || defined(__SSE2__) + unsigned csr = _mm_getcsr(); +#endif + ~RestoreEnvironment() { + std::fesetround(rounding); +#if defined(_M_X64) || defined(__SSE2__) + _mm_setcsr(csr); +#endif + } + } restore; + ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0); + const auto bad_environment = element.EvaluateLinearHighLow(high, low); + ASSERT_FALSE(bad_environment.HasValue()); + EXPECT_EQ(bad_environment.GetStatus().Diagnostics()[0].message, + "Hex8 high+low requires nearest rounding and gradual underflow."); + ASSERT_EQ(std::fesetround(FE_TONEAREST), 0); +#if defined(_M_X64) || defined(__SSE2__) + _mm_setcsr(restore.csr | 0x8040U); + EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue()); + _mm_setcsr(restore.csr); +#endif + const auto recovered = element.EvaluateLinearHighLow(high, low); + ASSERT_TRUE(recovered.HasValue()); + high[3] = -2; + EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue()); + high[3] = 0; + const auto after_failure = element.EvaluateLinearHighLow(high, low); + ASSERT_TRUE(after_failure.HasValue()); + ExpectSamePairResponse(after_failure.Value(), recovered.Value()); +} + TEST(Hex8HighLow, LowOnlyAffineStrainSurvivesUniformHighTranslation) { auto brick = UnitBrick(); fesa::Vector high(24), low(24);