From e878379c30c0e991e5ef531b886ecd87b7994ce0 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 9 Oct 2026 23:55:01 +0900 Subject: [PATCH] perf(hex8): cache immutable paired linear Hessian --- docs/hex8-structural-analysis/lessons.md | 7 ++- .../numerical-review.md | 3 ++ include/fesa/elements/hex8_element.h | 5 +++ src/fesa/elements/hex8_element.cpp | 7 +++ src/fesa/elements/hex8_high_low.cpp | 13 +++++- tests/unit/elements/hex8_high_low_test.cpp | 44 +++++++++++++++++++ 6 files changed, 76 insertions(+), 3 deletions(-) diff --git a/docs/hex8-structural-analysis/lessons.md b/docs/hex8-structural-analysis/lessons.md index 7fde4ba..10b402e 100644 --- a/docs/hex8-structural-analysis/lessons.md +++ b/docs/hex8-structural-analysis/lessons.md @@ -241,7 +241,7 @@ Production의 principal refinement·history·ULP 탐색은 이 비교에 포함 2026-10-09 승인 A 후, low를 누락하는 대조 kernel에서 실제 요소 assertion3개가 실패했다. 전체 pair로 H/Green strain/응력/힘/에너지/접선을 계산한 kernel은 native build0/CTest0이다. -원래 요소27개와 새 검사15개의 최종 묶음42/42가 통과했고, quarter-turn 작은 전단의 힘 +원래 요소27개와 새 검사17개의 최종 묶음44/44가 통과했고, quarter-turn 작은 전단의 힘 covariance와 독립 exact-dyadic primitive oracle을 포함한다. Native 산술6/6은 RN/FTZ/DAZ 실행 전제와 normal high에 최소 subnormal FMA error가 붙는 경계까지 확인했다. @@ -251,6 +251,11 @@ covariance와 독립 exact-dyadic primitive oracle을 포함한다. Native 산 중앙 에너지 차분의 StVK cubic 오차는 step²로 감소한다. 이를 exact 미분값으로 간주해 과도한 test bound를 적용했던 실패를 구분하고 분석적 truncation bound를 적용했다. +큰 반복 비용을 줄이기 위해 선형 Hessian은 같은 paired zero-state graph의 결과를 +immutable cache로 보관한다. Cache 생성 실패는 기존 Float64 Create 지원을 거부하지 않는다. +변형·실패 trial 후에도 high/low Hessian은 같고 힘·에너지·복구와 비선형 접선은 같은 전체 +trial에서 계산한다. 독립 reviewer가 이 캐시와 native44/44의 원시 로그도 확인했다. + 독립 numerical reviewer는 실제 shared arithmetic·kernel·geometry의 한정된 범위를 PASS로 판정했다. 재현 로그는 `.harness/hex8-evidence/precision/`에 보존한다. 이 결과는 동적 전체 상태·시간수렴·원본 endpoint·공식 출력의 완료를 뜻하지 않는다. 현재 checklist는 미완료다. diff --git a/docs/hex8-structural-analysis/numerical-review.md b/docs/hex8-structural-analysis/numerical-review.md index bc8d9a7..e874602 100644 --- a/docs/hex8-structural-analysis/numerical-review.md +++ b/docs/hex8-structural-analysis/numerical-review.md @@ -16,6 +16,9 @@ rule, isotropic material Hessian·geometric 항·Cauchy 복구와 full X+high+lo 산술 검사 로그를 확인했다. 해당 범위 bounded PASS이며 새 actionable finding은 없다. 일반 correctly-rounded106-bit 보장이나 T6 전체 state/시간수렴/원본 PASS가 아니다. T6 consumer·backend·worker 경계와 rollback은 해당 구현이 고정된 후 별도 검토한다. +이후 같은 reviewer가 optional immutable paired-linear Hessian cache를 검토했다. +동일 graph·선형 force/energy/recovery·비선형 trial Hessian·기존 F64 Create 지원을 유지하며 +native44/44 원시 로그를 확인했다. 이 캐시 변경도 bounded PASS다. ## HEX8-NR-001 — 영 회전 경계와 필수 원본 충돌 diff --git a/include/fesa/elements/hex8_element.h b/include/fesa/elements/hex8_element.h index 47f05c7..a2fd838 100644 --- a/include/fesa/elements/hex8_element.h +++ b/include/fesa/elements/hex8_element.h @@ -101,6 +101,11 @@ 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}; }; } // namespace fesa diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index 9c3d4ff..e7f711f 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -180,6 +180,13 @@ Result Hex8Element::Create(EntityIndex index, return Result::Failure( candidate.Failure("Hex8 reference stiffness has an unrepresentable " "positive diagonal.")); + 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 efba073..c6577be 100644 --- a/src/fesa/elements/hex8_high_low.cpp +++ b/src/fesa/elements/hex8_high_low.cpp @@ -151,12 +151,15 @@ Result Hex8Element::EvaluateHighLow( response.high.points[g] = {element_index_, static_cast(g + 1U), point.natural, point.position, Pack(stress), Pack(strain)}; - std::array fg{}; for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) { - fg[a][i] = Dot(f[i], gradients[a]); force[3U * a + i] += Dot(p[i], gradients[a]) * HighLow{point.weight}; } + 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) + fg[a][i] = Dot(f[i], gradients[a]); for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { const auto dot = Dot(gradients[a], gradients[b]); @@ -183,12 +186,18 @@ 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_; + } for (std::size_t i = 0; i < 24; ++i) { if (!IsNormalized(force[i])) return Result::Failure( 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; 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/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index 9029006..e837e0c 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -65,6 +65,26 @@ TEST(Hex8HighLow, NonlinearLowGradientProducesSecondOrderForce) { 1e-12 * strain * strain); } +TEST(Hex8HighLow, VolumetricAndTransverseTermsRetainLowStrain) { + auto brick = UnitBrick(0, 0.25); + fesa::Vector high(24), low(24); + const double strain = std::ldexp(1.0, -60); + for (std::size_t a = 0; a < 8; ++a) high[3*a] = 1; + for (auto a : {1, 2, 5, 6}) low[3*a] = strain; + const auto linear = brick.EvaluateLinearHighLow(high, low); + const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low); + ASSERT_TRUE(linear.HasValue()); ASSERT_TRUE(nonlinear.HasValue()); + // Uniform axial strain: lambda=mu=4.8, each nodal face integral is 1/4. + EXPECT_NEAR(linear.Value().high.internal_force[3], 3.6*strain, 1e-13*strain); + EXPECT_NEAR(linear.Value().high.internal_force[7], 1.2*strain, 1e-13*strain); + EXPECT_NEAR(linear.Value().high.strain_energy, 7.2*strain*strain, + 1e-13*strain*strain); + const auto quadratic = Subtract( + {nonlinear.Value().high.internal_force[3], nonlinear.Value().internal_force_low[3]}, + {linear.Value().high.internal_force[3], linear.Value().internal_force_low[3]}); + EXPECT_NEAR(quadratic.high, 5.4*strain*strain, 1e-12*strain*strain); +} + TEST(Hex8HighLow, RejectsInvalidLowShapeAndOverlappingComponents) { auto brick = UnitBrick(); fesa::Vector high(24), low(24); @@ -90,6 +110,30 @@ TEST(Hex8HighLow, LinearAndNonlinearZeroStateHaveSameFullHessian) { } } +TEST(Hex8HighLow, CachedLinearHessianIsIndependentOfTrialAndPriorFailure) { + auto brick = UnitBrick(0, 0.25); + fesa::Vector zero(24), high(24), low(24); + high[3] = 0.1; + low[3] = std::ldexp(1.0, -60); + const auto reference = brick.EvaluateLinearHighLow(zero, zero); + ASSERT_TRUE(reference.HasValue()); + EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue()); + const auto strained = brick.EvaluateLinearHighLow(high, low); + ASSERT_TRUE(strained.HasValue()); + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) { + EXPECT_EQ(reference.Value().high.tangent(i,j), strained.Value().high.tangent(i,j)); + EXPECT_EQ(reference.Value().tangent_low(i,j), strained.Value().tangent_low(i,j)); + } + const auto nonlinear_zero = brick.EvaluateNonlinearHighLow(zero, zero); + ASSERT_TRUE(nonlinear_zero.HasValue()); + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) { + EXPECT_EQ(reference.Value().high.tangent(i,j), nonlinear_zero.Value().high.tangent(i,j)); + EXPECT_EQ(reference.Value().tangent_low(i,j), nonlinear_zero.Value().tangent_low(i,j)); + } +} + TEST(Hex8HighLow, CoupledBricksTransmitLowStrainAcrossSharedFace) { const auto left = UnitBrick(), right = UnitBrick(1); fesa::Vector high(24), left_low(24), right_low(24);