From d5e08f868bbc8b288da6ed377839fa153aebb8ec Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Thu, 10 Sep 2026 02:24:04 +0900 Subject: [PATCH] feat: complete first MITC4 FEAST eigen implementation validation --- docs/ADR.md | 18 + docs/TOLERANCE.md | 20 + docs/mitc4-shell-feast-eigen2/ADR.md | 10 + .../FEATURED_AGENTS.md | 5 +- docs/mitc4-shell-feast-eigen2/PRD.md | 42 +- .../implementation-first-version-review.md | 67 +++ .../implementation-plan.md | 8 +- .../implementation-verification.md | 48 ++- docs/mitc4-shell-feast-eigen2/io-mapping.md | 22 +- docs/mitc4-shell-feast-eigen2/io.md | 11 +- .../numerical-review-cases.md | 10 +- .../reference-model.md | 22 +- docs/mitc4-shell-feast-eigen2/requirements.md | 4 +- .../solvers/eigen/mkl_feast_eigen_solver.cpp | 120 +++++- .../eigen/mkl_feast_eigen_solver_test_seam.h | 70 ++++ tests/CMakeLists.txt | 1 + .../reference/eigen_reference_comparison.cpp | 6 +- tests/reference/eigen_reference_comparison.h | 2 + .../eigen_reference_comparison_test.cpp | 28 ++ .../reference/reference_tolerance_policy.cpp | 39 +- tests/reference/reference_tolerance_policy.h | 7 + .../reference_tolerance_policy_test.cpp | 55 +++ .../analysis/eigen_mass_assembly_test.cpp | 180 +++++++++ .../mitc4_shell_mass_contract_test.cpp | 341 +++++++++++++++- .../hdf5/eigen_hdf5_results_writer_test.cpp | 18 + .../eigen/mkl_feast_eigen_solver_test.cpp | 380 +++++++++++++++++- 26 files changed, 1452 insertions(+), 82 deletions(-) create mode 100644 docs/mitc4-shell-feast-eigen2/implementation-first-version-review.md create mode 100644 src/fesa/solvers/eigen/mkl_feast_eigen_solver_test_seam.h create mode 100644 tests/unit/analysis/eigen_mass_assembly_test.cpp diff --git a/docs/ADR.md b/docs/ADR.md index 9d88805..a79b125 100644 --- a/docs/ADR.md +++ b/docs/ADR.md @@ -402,6 +402,24 @@ owner가 수정하며 새 근거나 변경 없이 같은 명령·agent 호출을 [Feature 문서 관리](FEATURE_DOCUMENTATION.md), CLI·설정과 실패 진단은 [Harness](HARNESS.md) 및 [동작 과정](HARNESS_WORKFLOW.md)을 따른다. + +### ADR-029: MITC4 eigen 첫 구현의 외부 Value 비교에 모드별 5%를 적용한다 + +- 상태: 사용자 명시 승인, 2026-09-10. Owner: 메인 assistant. +- 승인 근거: “5% 이내면 첫 구현이므로 Abaqus와의 고유치 정확도 향상은 추후에 진행”, “구현을 마무리”. +- 결정: [TOLERANCE](TOLERANCE.md)의 `mitc4-eigen-first-implementation` 예외를 + `mitc4-shell-feast-eigen2` 외부 scalar Value 비교에만 적용한다. 기존 모드별 blocking + singleton의 상대오차와 RMS를 모두 inclusive5%로 판정한다. +- 이유: 첫 구현의 외부 정확도 목표를 사용자 선택에 맞춘다. 모든 mode를 한 family로 + 합치지 않고 각 모드에서 독립적으로5%를 검사한다. +- 영향: 일반 `ReferenceTolerancePolicy::Evaluate`의 RMS1%는 유지하며, eigen 전용 named + entry point가 같은 metric/사전검사와 RMS5%를 사용한다. Comparator report에 적용 정책과 + RMS 한계를 기록한다. HDF5 schema·reference bytes·내부 잔차·기존 static 비교는 유지한다. +- 한계: 외부 고유치 정확도 향상과 S3는 후속 과제다. 나머지 필수 correctness/실패 보존 + 검사를 면제하거나 release를 승인하지 않는다. 1% 기준의 과거 실패를 소급 PASS로 바꾸지 않는다. +- 검증: [Feature ADR-012](mitc4-shell-feast-eigen2/ADR.md#mfe2-adr-012--첫-구현의-외부-고유치-5-수용)과 + [구현 결과](mitc4-shell-feast-eigen2/implementation-verification.md)에 경계·음성·회귀와 fresh 원본 비교를 기록한다. + ## 후속 ADR 작성 기준 프로젝트 전체의 dependency 방향, 공통 identity·수명·결과 의미, 기술 baseline, diff --git a/docs/TOLERANCE.md b/docs/TOLERANCE.md index 1733aec..139c984 100644 --- a/docs/TOLERANCE.md +++ b/docs/TOLERANCE.md @@ -6,6 +6,7 @@ - status: `approved-and-implemented` - effective_date: `2026-08-18` - decision_record: `docs/ADR.md`의 ADR-022 +- scoped_exception: `mitc4-eigen-first-implementation`, 2026-09-10 사용자 승인, [ADR-029](ADR.md#adr-029). ## 목적 @@ -34,6 +35,25 @@ component 구성과 최종 판정 영향은 해당 기능의 `requirements.md`, 독립적인 absolute-error tolerance는 사용하지 않는다. Absolute error는 near-zero 행을 판정하고 결과를 진단하기 위해서만 사용한다. +## 승인된 MITC4 eigen 첫 구현 예외 + +2026-09-10 사용자는 `mitc4-shell-feast-eigen2`의 첫 구현에서 Abaqus 고유치와의 차이가 +5% 이내이면 수용하고 정확도 향상은 후속 진행하도록 승인했다. 해당 기능의 외부 scalar +`Value` 비교에만 named policy **`mitc4-eigen-first-implementation`**을 적용한다. + +- 기존 source step/modal frame별 **blocking singleton Value family**를 유지한다. +- 각 모드의 상대오차는 `<=0.05`, family relative RMS도 `<=0.05`로 판정한다. + 비영 scalar singleton에서는 두 metric이 같으며 어느 모드도 다른 모드의 scale로 판정하지 않는다. +- 아래 scale·near-zero·zero-scale exact-zero·identity/finite 사전검사 및 결정적인 보고를 유지한다. +- 비교 보고서는 `policy_id=mitc4-eigen-first-implementation`과 + `relative_rms_tolerance=0.05`를 기록하여 기존 공통 RMS1% 판정과 구별한다. +- 다른 기능과 MITC4 정적 U/UR/reaction 비교는 기존 공통 상대오차5%·RMS1%를 그대로 사용한다. + 내부 eigen 잔차·정규화·직교성 `1e-10`도 변경하지 않는다. + +이 예외는 고유치의 외부 acceptance 변경이며 reference 보정·물리식 변경·S3 해소 또는 +Abaqus 내부 알고리즘 동등성의 근거가 아니다. 기존 1% 기준의 실패 기록은 당시 evidence로 +보존한다. 예외의 축소·폐지와 새 기준 적용도 아래 변경 관리 절차를 따른다. + ## Comparison family 수치 scale은 개별 행이나 component마다 만들지 않고 comparison family마다 계산한다. diff --git a/docs/mitc4-shell-feast-eigen2/ADR.md b/docs/mitc4-shell-feast-eigen2/ADR.md index e26c04f..0b5d770 100644 --- a/docs/mitc4-shell-feast-eigen2/ADR.md +++ b/docs/mitc4-shell-feast-eigen2/ADR.md @@ -97,3 +97,13 @@ - 영향: [정식화 E5](formulation-eigen.md), [수치 controls C5](numerical-review-controls.md), backend와 HDF5 writer의 control 검증을 함께 갱신한다. 기존 schema의 실제 `base_fpm[5]` 기록값이 바뀌며 field·shape·단위와 reference mapping은 유지한다. - 검증: [실행 결과](implementation-verification.md) 및 [독립 수치 검토](implementation-feast-residual-review.md)에 새 실행을 기록한다. 기존 trace 실패 evidence는 당시 결과로 보존한다. - 범위: S3 정식화·부분 회전 BC·reaction·평형 기준은 변경하지 않고 알려진 충돌을 후속 과제로 유지한다. S3를 해결하거나 unsupported로 재분류한 결정이 아니다. 외부 reference와 공통 tolerance도 유지한다. + +## MFE2-ADR-012 — 첫 구현의 외부 고유치 5% 수용 + +- 상태: 사용자 명시 승인, 2026-09-10. Owner: 메인 assistant. +- 승인 근거: “5% 이내면 첫 구현이므로 Abaqus와의 고유치 정확도 향상은 추후에 진행”, “구현을 마무리”. +- 결정: [공통 ADR-029](../ADR.md#adr-029)와 [TOLERANCE 예외](../TOLERANCE.md#승인된-mitc4-eigen-첫-구현-예외)의 `mitc4-eigen-first-implementation` 정책을 적용한다. Mode별 Value relative/RMS limit은 inclusive5%다. +- 영향: Requirement022, reference model, I/O mapping과 comparator/report를 갱신한다. 같은 모드 identity·singleton family·원본 artifact와 내부 all-c1e-10 기준을 유지한다. 다른 기능의 공통 RMS1%는 변경하지 않는다. +- 범위: 첫 구현의 외부 고유치 정확도 향상 및 ADR-011의 S3는 후속 과제로 유지한다. 나머지 필수 구현/실패 보존 검사는 실제 evidence로 마무리하며 미실행을 통과로 표시하지 않는다. +- 결과: [구현 결과](implementation-verification.md)와 [첫 구현 독립 검토](implementation-first-version-review.md)에 새 정책의 fresh 검증 및 제한을 기록한다. 과거1% 비교 실패는 당시 계약의 이력으로 보존한다. +- 통합 승인: 사용자는 구현·검증 완료 후 `dev` 병합과 push를 명시 요청했다. 메인은 통과한 변경을 커밋하고 `dev`에 반영하며 관련 없는 사용자 worktree 변경은 보존한다. 외부 release/tag 승인으로 확대하지 않는다. diff --git a/docs/mitc4-shell-feast-eigen2/FEATURED_AGENTS.md b/docs/mitc4-shell-feast-eigen2/FEATURED_AGENTS.md index bf2f74c..de26d68 100644 --- a/docs/mitc4-shell-feast-eigen2/FEATURED_AGENTS.md +++ b/docs/mitc4-shell-feast-eigen2/FEATURED_AGENTS.md @@ -2,8 +2,8 @@ - Feature ID: `mitc4-shell-feast-eigen2` - Owner: 메인 assistant -- Requirement baseline revision: `r1` -- 현재 작업: requirement/formulation/I/O r1과 [ADR-011](ADR.md#mfe2-adr-011--feast-residual-stopping과-s3-후속-분리)의 FEAST residual stopping 변경을 적용한다. C5 controls는 r2이며 S3는 후속 과제로 유지한다. +- Requirement baseline revision: `r2` (외부 고유치5% 수용) +- 현재 상태: [ADR-011](ADR.md#mfe2-adr-011--feast-residual-stopping과-s3-후속-분리)의 residual stopping과 [ADR-012](ADR.md#mfe2-adr-012--첫-구현의-외부-고유치-5-수용)의 외부5% 기준으로 첫 구현 검증을 완료했다. 전체417개 검사와 독립 검토는 PASS다. C5 controls·I/O mapping은 r2이며 S3와 정확도 향상은 후속 과제로 유지한다. - 목적: MITC4 질량행렬과 FEAST 고유진동 해석의 요구사항을 이론·수학·수치 검토·입출력 계약으로 구체화한다. - 현재 checklist: [PRD](PRD.md#현재-checklist). @@ -41,5 +41,6 @@ | [구현 수치 검토](implementation-numerical-review.md), [검토 이력](implementation-numerical-review-history.md), [S3 정적 검토](implementation-numerical-review-s3-static.md) | numerical-review-agent | 구현의 독립 수치 대조와 계약 충돌 | | [구현 물리 검토](implementation-physics-review.md) | physics-evaluation-agent | 실제 retained 결과의 물리 evidence·범위와 기존 문제 | | [FEAST 종료 제어 검토](implementation-feast-residual-review.md) | numerical-review-agent | 사용자 승인 residual stopping의 독립 검토와 새 실행 판정 | +| [첫 구현 검토](implementation-first-version-review.md) | numerical-review-agent | 승인된5% 범위·필수 coverage 보강·실행 및 제한의 독립 확인 | 반환은 결과·변경 파일·정적 확인·미해결 항목으로 충분하다. Git 경로 선택·커밋과 최종 결과 확인은 메인이 수행한다. diff --git a/docs/mitc4-shell-feast-eigen2/PRD.md b/docs/mitc4-shell-feast-eigen2/PRD.md index 29d6b18..bb92be0 100644 --- a/docs/mitc4-shell-feast-eigen2/PRD.md +++ b/docs/mitc4-shell-feast-eigen2/PRD.md @@ -2,8 +2,8 @@ - Feature ID: `mitc4-shell-feast-eigen2` - Owner: 메인 assistant -- Requirement baseline revision: `r1` -- 상태: requirement/formulation/NR/I/O r1 문서 완료. 후속 사용자 요청에 따라 implementation·실행 검증 진행 중이다. +- Requirement baseline revision: `r2` (022의 외부 accuracy 기준 변경) +- 상태: **사용자 승인 첫 구현 완료**. 전체 417개 검사와 모드별 외부5% 비교, 독립 수치 검토를 통과했다. S3와 후속 정확도 개선은 보류한다. 기준: [프로젝트 지침](../../AGENTS.md), [프로젝트 PRD](../PRD.md), [공통 ADR](../ADR.md#adr-028). 지역 지침은 [FEATURED_AGENTS](FEATURED_AGENTS.md), @@ -16,7 +16,7 @@ 변경사항을 Git에 커밋하도록 요청했다. 이어 **모든 단계를 진행하지 말고 requirement 단계만 진행**하도록 범위를 명시했다. -이후 사용자가 **research부터 io-definition까지** 진행하도록 요청했다. 현재 범위는 research, +이후 사용자가 **research부터 io-definition까지** 진행하도록 요청했다. 당시 범위는 research, formulation, 독립 numerical review, I/O 계약 문서와 단계별 커밋이다. 구현 계획·코드, solver build/test·reference 비교와 release 감사는 포함하지 않는다. 별도 coordinator나 phase 상태기계를 만들지 않는다. @@ -27,12 +27,19 @@ solver build/test·reference 비교와 release 감사는 포함하지 않는다. 2026-09-09 KST 사용자는 **implementation 단계를 진행**하도록 요청하고 **implementation 단계의 30분 제한은 없다**고 명시했다. 현재 범위는 [구현 계획](implementation-plan.md)에 따른 코드·테스트·실제 FEAST 및 필수 reference/회귀 검증과 작업별 커밋이다. 앞의 30분 제한은 -완료된 문서 작업의 이력이다. 원본 reference 생성·변경과 push/release는 포함하지 않는다. +완료된 문서 작업의 이력이다. 원본 reference 생성·변경과 release는 포함하지 않는다. Push는 아래 후속 요청으로 승인되었다. 2026-09-10 사용자는 FEAST 종료 조건을 trace에서 residual로 변경하고 S3는 추후 고려하도록 명시했다. [ADR-011](ADR.md#mfe2-adr-011--feast-residual-stopping과-s3-후속-분리)에 따라 backend·저장 control 검증·관련 문서를 갱신하며, 기존 독립 잔차와 외부 tolerance, 질량·drilling 수식은 유지한다. +이어 사용자는 “5% 이내면 첫 구현이므로 Abaqus와의 고유치 정확도 향상은 추후에 진행”하고 +구현을 마무리하도록 요청했다. [ADR-012](ADR.md#mfe2-adr-012--첫-구현의-외부-고유치-5-수용)의 +모드별 외부5% 정책을 적용하고 나머지 필수 검증을 보강한다. S3 보류는 유지한다. + +사용자는 구현과 검증 완료 후 `dev` 병합 및 push도 명시 승인했다. 관련 변경을 커밋하고 +검증한 tree를 `dev`에 병합한 뒤 원격 `dev`에 반영한다. 기존 사용자 변경과 reference 원본은 보존한다. + ## 솔버 기능의 요구 범위 요약 기존 사용자 결정에 기록된 다음 의미를 새 요구사항의 출발점으로 사용한다. @@ -79,24 +86,29 @@ solver build/test·reference 비교와 release 감사는 포함하지 않는다. 이 checklist는 문서 작업 상태다. Solver 구현·실제 FEAST 실행·회귀 및 외부 reference 비교의 통과를 뜻하지 않는다. -- [ ] 질량·stabilization과 독립 analytical/invariant 테스트 구현·검증. -- [ ] 실제 FEAST backend·count/cluster/zero/실패 controls와 runtime 검증. +- [x] 질량·stabilization과 독립 analytical/invariant 테스트 구현·검증. S3 public static 충돌은 아래 후속 항목으로 분리한다. +- [x] 실제 FEAST backend·count/cluster/zero/실패 controls와 runtime 검증. 실행한 대표 실패 경로와 미주입 한계는 결과 보고서에 구분한다. - [x] Density/Frequency semantic·source identity·진단 및 입력 테스트. 기존 정적 PREPRINT 회귀를 보존했다. -- [ ] Constraint-first 축소·modal recovery·별도 eigen state/procedure 구현. -- [ ] HDF5 additive schema·candidate 보존·Value comparator 및 통합 테스트. -- [ ] 제외 예제 사용 없이 남는 static/내부 회귀·전체 build/test·필수 reference 실행. -- [ ] 관련 독립 numerical/physics 확인, 결과·한계 기록과 커밋. +- [x] Constraint-first 축소·modal recovery·별도 eigen state/procedure 구현. +- [x] HDF5 additive schema·candidate 보존·Value comparator 및 통합 테스트. +- [x] 제외 예제 사용 없이 남는 static/내부 회귀·전체 build/test·필수 reference 실행: 417/417 PASS. +- [x] 관련 독립 numerical/physics 확인과 결과·한계 기록. 첫 구현 독립 판정은 승인된 범위에서 PASS다. - [x] 승인된 FEAST residual stopping과 HDF5 control 검증을 함께 적용하고 변경 전 실패·변경 후 검사를 확인한다. - [x] 원본 eigen의 전체 c개 독립 잔차·복구·HDF5 생성을 확인하고 종료 제어의 독립 수치 검토를 받는다. -- [ ] 원본 eigen의 외부 Value 불일치 원인을 조사한다. 현재 Mode2–10의 모드별 RMS가 1%를 초과한다. +- [x] 원본 eigen 10개 Value가 사용자 승인 모드별5% 기준을 충족한다. 최대 상대오차는 약4.261331%다. + +후속 보류 항목: + +- [ ] Abaqus와의 외부 Value 정확도를 개선한다. 이전 RMS1% 기준의 Mode2–10 불일치는 이력으로 보존한다. - [ ] S3 public static 정적 평형 계약 충돌은 사용자 요청에 따라 후속 검토한다. 현재 구현·실행은 [implementation-verification](implementation-verification.md)에 기록했다. -실제 mass/reduction/FEAST/procedure/HDF5 경로와 focused 검사가 있으며, 체크되지 않은 항목이 -전부 미구현이라는 뜻은 아니다. 원본 eigen의 전체 c residual 실패는 residual stopping으로 -해소했다. 현재 외부 Value 불일치와 S3 public static 계약 충돌, 세부 case의 미완료 evidence가 -남으므로 기능 전체 acceptance는 완료하지 않았다. +실제 mass/reduction/FEAST/procedure/HDF5 경로와 전체 실행을 확인했고 원본 eigen의 전체 c +residual 실패는 residual stopping으로 해소했다. [첫 구현 독립 검토](implementation-first-version-review.md)는 +승인된5%·S3 보류 범위에서 PASS다. 실제512회 vendor 호출과 HDF5 close-error 주입은 수행하지 +않았으며 대표 실패 경로·정적 대조의 범위를 보고서에 기록했다. 기존 nonlinear recovery 문제와 +제외 pressure 예제의 coverage 부재를 해결한 것으로 해석하지 않는다. ## 이전 문서 단계의 결과와 구현 입력 diff --git a/docs/mitc4-shell-feast-eigen2/implementation-first-version-review.md b/docs/mitc4-shell-feast-eigen2/implementation-first-version-review.md new file mode 100644 index 0000000..f99c103 --- /dev/null +++ b/docs/mitc4-shell-feast-eigen2/implementation-first-version-review.md @@ -0,0 +1,67 @@ +# 첫 구현 acceptance 독립 검토 + +- Feature: `mitc4-shell-feast-eigen2`; Owner: `numerical-review-agent`; 2026-09-10. +- 판정: **`PASS` — 사용자 승인 첫 구현 범위**. S3와 후속 정확도 개선은 보류 상태다. +- 적용: 사용자 승인 외부 Value 5%, [PRD](PRD.md), [ADR](ADR.md), + [TOLERANCE](../TOLERANCE.md), [reference-model](reference-model.md), + [필수 cases](numerical-review-cases.md), [내부 controls](numerical-review-controls.md). + +## 승인 경계와 수치 판단 + +원본 Mode1–10은 source `Step-2`의 독립 modal-frame singleton이다. 첫 구현 profile +`mitc4-eigen-first-implementation`은 각 `Value`의 reference-only relative error와 RMS를 +각각 .05로 검사한다. 경계 포함이며 family를 합치거나 mode를 평균하지 않는다. +다른 공통 comparison의 RMS .01, missing/extra/duplicate/nonfinite 거부, zero-scale exact-zero, +내부 전체 c개의 residual/Gram `1e-10`은 유지한다. + +[종료 제어 검토](implementation-feast-residual-review.md)의 원본 최대 Value 오차 +`4.2613310845%`는 승인된 새 5% 범위 안이다. 새 comparator의 실제 focused 실행에서 +원본 10개 family가 통과했다. 기존 1% FAIL 기록은 당시 판정으로 보존한다. 수식·DOF·단위·적분· +mass·whitening·복구·state update는 이번 외부 정책으로 바뀌지 않는다. + +## 필수 보강과 추가 검증 한계 + +| 분류 | 요구·현재 evidence·필요한 보강 | +| --- | --- | +| 확인 — backend 실패 | REQ012/019, C21/C22. Fatal count·optional info1 fallback, final count 불일치·변경 없는 solve 반복 금지, info2의 더 큰 구간 retry, 요청 N 밖 고유벡터의 실제 residual 초과 및 nonfinite/Gram 거부 통과. Production이 사용하는 동일 budget helper의 511/512·363/364, size product·PTRDIFF/container 경계, 실제 count/solve workspace bad_alloc·thread 복원 및 초기 count3회 한계를 검사 | +| 확인 — raw·순위·반복군 | C18/C19/C21. 증가 upper의 감소 count 거부, raw -1e-12와 -1e-9 보존 및 -1e-9 바로 아래 실패 통과. 작은 isolated repeated군을 세 번 실제 FEAST로 풀어 순서·history와 승인된 Frobenius projector/sqrt(2r) metric을 비교 | +| 확인 — mass oracle·단위 | REQ005/014/026, C02/C03/C05. 48개 rational Ck/각 row·first spatial moment, theta×d 부호의 C1 cross-block·C2 rotary block과 두24×24 질량, pairwise Jensen의1/2 및2*C1 항을 대조. Scalar는 m_e(t/2)^k, matrix는 D Frobenius norm합, Jensen 차이는 qMLq+qMCq로 정규화. Proper Q+shift를 geometry/director부터 재계산하여 두 mass의 공변성 통과 | +| 확인 — 공유 node assembly | REQ003/004/007, C06. rho2/t.5와 rho5/t1의 공유 node URX만 free인 patch에서 mass7/64, local K합/m와 modal amplitude²*mass1, element source order 변화의 결과를 확인. Model-level BC로 fixture를 바로잡은 뒤 통과; parser 지원 범위를 넓히지 않음 | +| 확인 — final 보존 | REQ019, C24. Eigen candidate 변조·empty pencil 보존 및 실제 valid HDF5에 Win32 delete-sharing lock을 걸어 close/self-check 뒤 replace 실패·원본 bytes 보존 통과. 공통 AtomicFinalizeValidatedCandidate의 실패 경로와 연결 | +| 확인 — 5% 한정 정책 | C25. ±5% 경계와 바로 밖, zero/nonfinite·identity 거부, common RMS1% 보존, actual profile metadata와 원본 10 family fresh 비교 통과. 아래 focused evidence와 source 대조 | +| 사용자 명시 보류 | S3의 public static moment 충돌. Direct pencil/covariance 검사의 성공과 구별하고 지원 정적 사례의 성공으로 표시하지 않는다 | +| 추가 검증의 한계 | 모든 조건수·geometry injectivity·mesh spectrum 수렴차수·vendor/OS 오류 조합·임의 반복군 방향은 보장하지 않는다. H5Fclose 오류 자체와512회 actual vendor 호출은 미실행이다. CloseChecked 실패는 RequireHdf5→catch/cleanup→replace전 종료를 정적으로 대조했고 budget은 실제 사용 helper의 경계 검사다. 대표 분기 검증을 모든 오류 조합 전수 검사로 확대하지 않는다 | + +검토 대상 tests: [mass contract](../../tests/unit/elements/mitc4_shell_mass_contract_test.cpp), +[mass/kernel](../../tests/unit/elements/mitc4_shell_test.cpp), +[공유 node assembly](../../tests/unit/analysis/eigen_mass_assembly_test.cpp), +[reduction](../../tests/unit/analysis/eigen_reduction_test.cpp), +[backend](../../tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp), +[procedure](../../tests/integration/analysis/linear_frequency_procedure_test.cpp), +[writer](../../tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp), +[common finalize](../../tests/unit/io/hdf5/hdf5_writer_components_test.cpp). +Source의 guard 존재와 실제 실패 주입 통과를 구분한다. 검토자는 source/tests·상위 계약·reference를 +수정하지 않고 메인의 fresh 결과를 대조한다. + +## 실행 evidence + +- `.harness/eigen-first-implementation-focused.log`: comparator/policy 10/10 PASS, + 원본 실제 해석→HDF5→10개 family 비교 포함. 정책 코드의 공통 Evaluate는 RMS .01을 + 유지하고 named entry는 blocking singleton `Value`만 허용한다. +- 독립 재계산: [기존 rational coefficient 적분](numerical-review-evidence.md)의 Python을 + 표준 라이브러리로 새 실행하고 mass test의 48개 rational 상수를 Fraction으로 대조하여 + 전부 exact 일치를 확인했다. 이종 두 요소 공유 node의 rotary mass도 + `2*(1/2)^3/48+5/48=7/64 [M L²]`로 독립 확인했다. 이는 C++ 실행 대체가 아니다. +- `.harness/eigen-first-implementation-coverage.log`: 31개 중30 PASS. Backend24개·writer + 통과, C06 fixture의 BC 위치만 실패했다. `.harness/eigen-first-implementation-mass.log`에서 + fixture와 정규화 보강 뒤 mass6/6 PASS를 직접 확인했다. +- `.harness/eigen-first-implementation-boundaries.log`: raw negative와 count monotonicity 경계2/2 PASS. +- `.harness/eigen-first-implementation-full-run.log`: MSVC x64 Debug의 configure/build/ + discovery/CTest 모두 exit0, **417/417 PASS**, 166.86초. Unit339/integration42/reference36. + 같은 build의 `Testing/Temporary/LastTest.log`에서도 Passed417/Failed0와 원본 reference + 실행을 직접 확인했다. 과거 C06 실패를 담은 `LastTestsFailed.log`는 최신 판정 근거로 사용하지 않는다. + +이번 검토에서 추가 numerical correction은 필요하지 않다. 독립 내부 기준과 새 외부 +정책의 실행 연결을 확인했으며, 모든 vendor·OS 실패 조합이나 S3 해결·Abaqus 내부 절차 +동등성을 승인한 판정은 아니다. [이전 물리 검토](implementation-physics-review.md)의 기존 +nonlinear recovery 제한도 이 변경으로 해소됐다고 주장하지 않는다. diff --git a/docs/mitc4-shell-feast-eigen2/implementation-plan.md b/docs/mitc4-shell-feast-eigen2/implementation-plan.md index 246a0bf..ec44567 100644 --- a/docs/mitc4-shell-feast-eigen2/implementation-plan.md +++ b/docs/mitc4-shell-feast-eigen2/implementation-plan.md @@ -32,7 +32,8 @@ 4. 영향받는 shell/shellNL/linear-limit, nonlinear tangent/objectivity/평형·state, pressure synthetic·기타 기존 회귀를 보존한다. 제외된 복합 테스트의 retained 관찰은 분리한다. 5. 기능 완료 시 구성된 전체 build/test와 필수 reference를 한 번 실행하고 독립 물리 검토와 - 실제 결과·실패·미실행을 기록한다. 실패 수치를 expected로 복사하거나 tolerance를 완화하지 않는다. + 실제 결과·실패·미실행을 기록한다. 실패 수치를 expected로 복사하지 않는다. 외부 고유치 + 기준은 [ADR-012](ADR.md#mfe2-adr-012--첫-구현의-외부-고유치-5-수용)의 사용자 승인5%를 적용한다. 현재 `.harness/config.json`은 `mfe-workspace-msvc-ninja-debug` configure/build/test preset과 `.harness/build-ninja-msvc-debug`를 선택한다. 개인 preset의 dependency 경로는 commit하지 않는다. @@ -40,5 +41,6 @@ Focused 실행은 해당 CMake target/CTest filter, 전체 검증은 `python scr 메인이 빌드 실행을 직렬화한다. C++에는 MSVC x64 Debug `/W4 /WX` 기준을 적용한다. 실제 runtime·known spectra와 원본 모델을 실행했다. 현재 결과·수정 범위와 미완료 검증은 -[implementation-verification](implementation-verification.md)에 기록한다. 원본 eigen all-c -residual과 S3 정적 계약 충돌을 사용자 결정 경계로 반환하며 기능 완료로 표시하지 않는다. +[implementation-verification](implementation-verification.md)에 기록한다. All-c residual은 +ADR-011의 residual stopping으로 해소했다. ADR-012의 첫 구현 범위에서 외부5% 비교와 +필수 mass/실패 보존 evidence를 확인하며 S3는 명시된 후속 과제로 유지한다. diff --git a/docs/mitc4-shell-feast-eigen2/implementation-verification.md b/docs/mitc4-shell-feast-eigen2/implementation-verification.md index ab2635b..8ace32b 100644 --- a/docs/mitc4-shell-feast-eigen2/implementation-verification.md +++ b/docs/mitc4-shell-feast-eigen2/implementation-verification.md @@ -1,13 +1,57 @@ # MITC4 eigen 구현 및 실행 결과 - Feature: `mitc4-shell-feast-eigen2`; Owner: 메인 assistant; 갱신 2026-09-10 KST. -- 현재 판정: **FEAST 잔차 초과 해소, 외부 Value 비교와 기능 acceptance 미완료**. +- 현재 판정: **사용자 승인 첫 구현 범위 PASS**. FEAST 잔차 초과를 해소했고 외부10개 Value가 모드별5% 기준을 충족한다. S3와 후속 정확도 개선은 보류한다. 아래 2026-09-10 결과가 최신이며 2026-09-09 trace 실행은 당시 evidence로 보존한다. - 기준: [프로젝트 지침](../../AGENTS.md), [PRD](PRD.md), [ADR](ADR.md), [정식화](formulation.md), [수치 controls](numerical-review-controls.md), [I/O](io.md). - 독립 검토: [수치](implementation-numerical-review.md), [S3 정적 계약](implementation-numerical-review-s3-static.md), - [물리](implementation-physics-review.md), [FEAST 종료 변경](implementation-feast-residual-review.md). + [물리](implementation-physics-review.md), [FEAST 종료 변경](implementation-feast-residual-review.md), [첫 구현 최종 검토](implementation-first-version-review.md). + +## 2026-09-10 첫 구현의 5% 수용과 검증 보강 + +사용자는 첫 구현의 Abaqus Value 차이를 모드별5%까지 수용하고 정확도 향상과 S3는 후속으로 +진행하도록 명시했다. [공통 ADR-029](../ADR.md#adr-029), [Feature ADR-012](ADR.md#mfe2-adr-012--첫-구현의-외부-고유치-5-수용)에 따라 +`mitc4-eigen-first-implementation`의 상대오차/RMS를 각각5%로 적용한다. 기존 singleton +identity·공통 사전검사와 내부 all-c1e-10은 유지하며 다른 기능은 기존 RMS1%다. +비교 보고서에는 policy_id와 relative_rms_tolerance를 기록한다. 아래 과거1% 판정은 소급 수정하지 않는다. + +| 보강 | 실제 관찰 | +| --- | --- | +| C25 외부 policy | 합성 output의 ±5% 경계와 바로 밖, zero/nonfinite·nonmodal family 거부, 다른 공통 policy RMS1% 보존. 변경 전3%/5% 비교 실패를 확인한 뒤 새 정책으로 통과 | +| C02/03 mass | Warped C0/C1/C2 48개 rational 성분·row·first spatial moment, 24×24 consistent/lumped cross block 및 독립 pairwise Jensen energy. 차원별 scale과 D-congruence Frobenius 사용 | +| C05/06 covariance·assembly | Proper Q+shift 재계산, 서로 다른 rho/t/material 두 요소의 공유 절점 mass7/64·modal normalization·source order 변경 | +| C19/21 backend | Raw tiny-negative/정확한 -tau_zero 보존과 경계 바로 밖 거부, fatal/count 모순·monotonicity·returned mismatch·nonfinite·요청 N 밖의 residual/Gram 실패, certified fallback·변경된 구간 retry | +| C22 budget·storage | Production-used private helper의 511/512·363/364·SIZE_MAX 및 곱/저장 한계, 실제 allocation try 안의 bad_alloc 주입, initial info1의 정확히3회 종료·MKL thread 복원 | +| C24 final 보존 | Eigen writer의 valid final을 Win32 delete-sharing lock으로 보호하여 실제 replacement 실패를 만들고 byte 보존·self-check 확인. 공통 finalizer의 기존 실패 테스트도 유지 | + +오류 주입은 `src/fesa/solvers/eigen/`의 private thread-local seam으로만 수행한다. +일반 실행에서는 hook이 없으며 vendor 호출·수치식·control·공개 API를 유지한다. Bounds/storage +helper는 기존 검사를 옮겨 production과 테스트가 같은 조건을 사용한다. 모든 vendor/OS 오류 +조합을 실행한 것은 아니며 actual512회 vendor 호출이나 실제 HDF5 close-error 주입은 하지 않았다. +Close 실패는 기존 checked-close→exception→candidate 정리/기존 final 보존 제어 흐름을 정적으로 대조했다. + +| 실행 | 결과 | +| --- | --- | +| policy 변경 전 focused | 기존 RMS1% 때문에3% 및 ±5% 사례 실패. `.harness/eigen-first-implementation-red.log` | +| policy/원본 focused | 10/10 PASS, 원본10개 Value 포함. `.harness/eigen-first-implementation-focused.log` | +| coverage 최초 focused | 30/31 PASS. C06 합성 BC 위치만 input 계약 위반으로 실패했고 model-level로 수정 | +| 수정 mass focused | 6/6 PASS. `.harness/eigen-first-implementation-mass.log` | +| 최종 monotonicity/raw 경계 focused | 2/2 PASS. `.harness/eigen-first-implementation-boundaries.log` | +| 구성된 전체 build/test | `python scripts/execute.py` exit0, configure/build/discovery/test 모두 exit0. **417/417 PASS**: unit339·integration42·reference36, 166.86초. `.harness/eigen-first-implementation-full-run.log`, 상세 `.harness/eigen-first-implementation-full-ctest.log` | + +최종 전체 실행 후 공식 HDF5를 read-only로 다시 추출했다 +(`.harness/eigen-first-implementation-metrics.log`). 원본10개 Value의 최대 상대오차는 +`4.261331%`, 전체 c=22 최대 독립 잔차는 `1.605467212112365e-13`으로 아래 residual +실행과 같다. 새 profile의 외부 비교와 내부 `1e-10` 기준을 모두 통과했다. Compiler warning을 +추가하지 않았으며 Python Harness suite는 변경 범위가 아니므로 실행하지 않았다. + +추가 테스트 작성 중 nested initializer와 lambda capture의 컴파일 오류를 수정한 뒤 MSVC x64 +Debug build가 성공했다. 실패한 빌드 뒤 등록되지 않은 테스트의 `No tests were found`를 PASS로 +계산하지 않았다. [첫 구현 독립 검토](implementation-first-version-review.md)는 새 oracle와 +failure/acceptance evidence·전체 실행을 대조하고 **사용자 승인 첫 구현 범위 PASS**로 판정했다. +임의 모델 정확도나 외부 release를 보장하는 판정은 아니다. ## 2026-09-10 FEAST residual stopping 변경 결과 diff --git a/docs/mitc4-shell-feast-eigen2/io-mapping.md b/docs/mitc4-shell-feast-eigen2/io-mapping.md index 24c6290..8104d0e 100644 --- a/docs/mitc4-shell-feast-eigen2/io-mapping.md +++ b/docs/mitc4-shell-feast-eigen2/io-mapping.md @@ -1,6 +1,6 @@ # MITC4 eigen 내부 mapping·reference projection -- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `io-definition-agent`; Contract revision: `r1`. +- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `io-definition-agent`; Contract revision: `r2` (외부 Value acceptance만 변경). - 기본 보고서: [io.md](io.md). 수학 원본: [formulation-eigen](formulation-eigen.md). ## 내부 소유권과 identity @@ -112,14 +112,24 @@ Tolerance 전에 다음을 순서대로 완료한다. [NR reference-model](reference-model.md)이 정의하는 family는 `(case,source step,modal-frame=mode_id,quantity=eigenvalue, dimension=time^-2,coordinates=scalar,impact=blocking)`다. Mode별 frame은 10개 singleton -family이며 하나의 storage table이라고 서로 합치지 않는다. [공통 TOLERANCE](../TOLERANCE.md)의 -reference-only S, inclusive near-zero/relative row branch, relative RMS, zero-scale 처리를 -그대로 적용한다. 현재 양수 singleton의 실효 기준은 mode별 RMS 1%이며 5% row만으로 통과하지 -않는다. S=0은 exact-zero FESA만 성공한다. I/O는 새 tolerance를 선택하지 않는다. +family이며 하나의 storage table이라고 서로 합치지 않는다. 2026-09-10 사용자 명시 5% 승인과 +[MFE2-ADR-012](ADR.md)에 따라 [공통 TOLERANCE](../TOLERANCE.md)의 승인된 한정 정책 +`mitc4-eigen-first-implementation`을 이 Feature의 외부 Value 비교에 적용한다. +Reference-only S와 near-zero 분류·row branch는 유지하고 상대오차와 relative RMS의 한계를 +각각 `0.05`로 적용한다. 같은 identity의 양수 singleton에서는 두 metric이 모두 +`abs(FESA - reference) / abs(reference)`이므로 mode별 5% 이하(inclusive)가 통과 조건이다. +Mode 사이 평균·collapse·재정렬로 실패를 숨기지 않는다. S=0은 exact-zero FESA만 성공하며, +missing/extra/duplicate/nonfinite/schema/identity 오류는 위 prechecks에서 먼저 거부한다. +다른 공통 family의 RMS `0.01`과 전체 c개 내부 독립 검증 `1e-10`은 유지한다. +이는 승인된 첫 구현 acceptance이며 후속 정확도 개선이나 S3 해결의 evidence가 아니다. Report는 input/reference/result identity, freshness evidence, 전체 prechecks, source key/ component/family/count/S, raw reference/FESA/error·branch·row verdict, RMS·family verdict, -blocking failure와 overall verdict를 mode 순서로 남긴다. Zero-scale 미정 metric은 named state로 +blocking failure와 overall verdict를 mode 순서로 남긴다. Comparator report는 +`policy_id=mitc4-eigen-first-implementation`과 `relative_rms_tolerance=0.05`를 명시해 +공통 RMS 1% 정책과 구별한다. 임의 tolerance를 지정하는 generic policy argument는 노출하지 않는다. +이 report 계약 변경은 HDF5 field·shape·단위·schema version을 변경하지 않는다. +Zero-scale 미정 metric은 named state로 표시하고 NaN/Inf/floor를 넣지 않는다. 필요 시 HDF5 projection CSV는 보조 debug view일 뿐이다. 외부 eigen gate는 Value뿐이다. N개의 global vectors는 내부 constraint/residual/normalization diff --git a/docs/mitc4-shell-feast-eigen2/io.md b/docs/mitc4-shell-feast-eigen2/io.md index d0e267a..127eae9 100644 --- a/docs/mitc4-shell-feast-eigen2/io.md +++ b/docs/mitc4-shell-feast-eigen2/io.md @@ -1,6 +1,6 @@ # MITC4 eigen I/O 계약 -- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `io-definition-agent`; Contract revision: `r1`. +- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `io-definition-agent`; Contract revision: `r2` (외부 Value acceptance만 변경). - 입력 baseline: [requirements r1](requirements.md), [formulation r1](formulation.md). - 기준: [프로젝트 지침](../../AGENTS.md), [PRD](../PRD.md), [ADR-028](../ADR.md#adr-028), [ARCHITECTURE](../ARCHITECTURE.md), [Feature ADR](ADR.md). @@ -14,6 +14,13 @@ lifecycle는 유지한다. Eigen discriminator·density·modal datasets는 addit N개의 global mass-normalized vectors는 내부 검증을 위해 필수이며 외부 shape gate는 없다. 전체 c개의 내부 validation evidence는 공개 mode N개와 분리한다. +2026-09-10 사용자 명시 5% 승인에 따라 첫 구현 외부 Value 비교에는 +`mitc4-eigen-first-implementation` 정책을 적용한다. 같은 identity의 mode별 singleton에서 +상대오차와 relative RMS를 모두 5% 이하(inclusive)로 판정한다. 근거는 +[공통 TOLERANCE](../TOLERANCE.md)와 [MFE2-ADR-012](ADR.md)이며 I/O가 새 기준을 선택한 것이 아니다. +다른 공통 comparison family의 RMS 1%, 전체 c개 내부 독립 검증 `1e-10`, HDF5 schema는 유지한다. +정확도 개선과 이미 분리한 S3는 후속 과제이며 이 acceptance 변경으로 해결을 주장하지 않는다. + | 상세 문서 | 권위와 내용 | | --- | --- | | [io-input](io-input.md) | Supported syntax/defaults, source scope, rho/BC mapping, no-op allowlist, fail-closed diagnostics | @@ -54,7 +61,7 @@ Requirement prefix는 `FEM-REQ-MITC4-EIGEN2-`다. 아래 검증을 문서 작업 | IO-T06 | 실제 FEAST known zero/repeated/cluster/n1/N=n·full capacity, optional abandoned count와 safe fallback, count/reset/budget/index/allocation/Gram/residual 실패; 최종 N과 c 구별 | 008–013,019 | | IO-T07 | 기존 accepted state/final 준비 후 solve/validation/log/lock/write/flush/close/reopen/self-check/replace failure injection; 기존 final 보존과 최초 실패 final 미생성 | 012,019 | | IO-T08 | TXT reorder·같은 Value의 다른 mode 허용; missing/extra/duplicate/nonfinite/malformed/step/quantity/unit/mass/backend/schema 오류를 tolerance 전에 거부; Freq 변화는 Value 판정 불변 | 020–021 | -| IO-T09 | Common inclusive row/RMS 경계·양옆, near-zero와 zero-scale; singleton은 row5% pass/RMS1% fail을 구별 | 022 | +| IO-T09 | Named first-implementation 정책의 singleton row/RMS 5% inclusive 경계·양옆, near-zero와 exact zero-scale; report policy_id/relative_rms_tolerance와 다른 공통 family RMS1% 유지 확인 | 022 | | IO-T10 | Density 없는 static 입력과 schema0/1 유지, retained shell/shellNL/linear-limit source identity·U/UR/reaction 분류와 fresh 회귀, 제외 예제0건 | 023–025,027 | | IO-T11 | Consistent→lumped 원본 연결과 같은 mass의 실제 FEAST/검증; mass oracle·inertia/null·coordinate/NR-MFE-001 및 deterministic assembly/row/diagnostic | 005–008,011,013–014,023,026–027 | diff --git a/docs/mitc4-shell-feast-eigen2/numerical-review-cases.md b/docs/mitc4-shell-feast-eigen2/numerical-review-cases.md index bf3967c..6c72254 100644 --- a/docs/mitc4-shell-feast-eigen2/numerical-review-cases.md +++ b/docs/mitc4-shell-feast-eigen2/numerical-review-cases.md @@ -1,8 +1,9 @@ # 필수 검증 사례와 요구사항 연결 -- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `numerical-review-agent`; revision `r1`. +- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `numerical-review-agent`; revision `r2` (2026-09-10, E2-C25의 승인된 첫 구현 profile 반영). - 기본 보고서: [numerical-review](numerical-review.md); 판정 수치: [controls](numerical-review-controls.md). -- 아래는 후속 구현의 필수 관찰이다. 이번에 실제 실행한 수학 oracle는 [evidence](numerical-review-evidence.md)뿐이다. +- 아래는 구현의 필수 관찰이다. 문서 r1 시점의 실제 수학 oracle는 [evidence](numerical-review-evidence.md)이며, + 최신 실행·미실행 분류는 [첫 구현 검토](implementation-first-version-review.md)를 따른다. [이전 case C01–C24](../mitc4-shell-feast-eigen/numerical-review-cases.md)의 synthetic fixture와 독립 oracle를 읽고 새 formulation에 대조해 아래에 채택한다. 숫자·fixture 상세는 그 문서를 @@ -36,7 +37,7 @@ Old C25–C27의 실행·portfolio 지시는 아래 새 정의로 대체한다. | E2-C22 | 실제 endpoint/cluster probes 및 mock count512/solve3/초기count3 한계, budget 직전 fallback·예약 불충분 실패. n²/nm0/nnz+1/index/byte overflow·allocation 실패. a0 below spectrum은 독립 bound로 검증; empty count 성공을 필수로 기대하지 않음 | | E2-C23 | 원본 single Step-2/N10/density/blank bounds/normalization, source Lanczos/SIM/acoustic-free 의미; unsupported preload/acoustic/mixed/damping/비영 BC 등의 silent no-op 금지. I/O allowlist로 구체화 | | E2-C24 | Source renumbering, HDF5 self-check와 final normalized modes/used mass/backend/count evidence; solve/validation/write/close/replace failure의 기존 state/final 보존. 정확한 layout은 I/O | -| E2-C25 | [Reference model](reference-model.md)의 원본 Value 10개, per-modal-frame singleton families·identity/finite 전검사·inclusive row/RMS/nearzero/zero scale; reordered rows/동일값 distinct IDs 허용, missing/extra/duplicate 거부 | +| E2-C25 | [Reference model](reference-model.md)의 원본 Value 10개, per-modal-frame singleton families·identity/finite 전검사. 승인된 `mitc4-eigen-first-implementation`에서 inclusive row/RMS 각각5%, 바로 밖은 실패; nearzero/zero scale 유지. 다른 공통 family RMS1% 보존; reordered rows/동일값 distinct IDs 허용, missing/extra/duplicate 거부 | | E2-C26 | Retained shell/shellNL/linear-limit fresh solve·기존 U blocking/UR warning·reaction review-only, density-free static, partial-BC response, Newton tangent/objectivity/force/moment/free residual/energy/rollback·pressure synthetic·기타 기존 회귀 | | E2-C27 | 향후 구현 시 구성 MSVC x64 Debug 전체 CMake/CTest와 actual FEAST/필수 reference를 한 번 연결. Source identity로 등록/선택/실행/제외를 구분하고 관련 없는 회귀 보존. 고정 과거 test 수 또는 중복 manual/automatic 실행 요구 없음 | @@ -44,7 +45,8 @@ E2-C17/18/22는 mock만으로 충족하지 않는다. E2-C21의 실패 주입은 바꾸는 행위가 아니라 adapter 경계 테스트다. Expected eigenvalues로 production bounds를 설정하거나 동일 production 함수를 oracle로 호출하지 않는다. Actual FEAST의 G01/G02 installed probes에는 occupied/empty/endpoint count 반환·flag reset·full m0=n/zero matrix· -workspace/ABI/thread 복원이 포함된다. 아직 **NOT RUN**이다. +workspace/ABI/thread 복원이 포함된다. r1 문서 시점에는 **NOT RUN**이었으며, +이후 실행 여부는 [첫 구현 검토](implementation-first-version-review.md)의 evidence와 구분한다. E2-C26은 excluded original/copied/cached examples를 절대로 사용하지 않는다. 제외된 혼합 시험에 있던 retained shellNL force/moment balance·free force/moment residual diff --git a/docs/mitc4-shell-feast-eigen2/reference-model.md b/docs/mitc4-shell-feast-eigen2/reference-model.md index 25cff20..86e5b6e 100644 --- a/docs/mitc4-shell-feast-eigen2/reference-model.md +++ b/docs/mitc4-shell-feast-eigen2/reference-model.md @@ -1,8 +1,8 @@ # 외부 reference projection과 tolerance 적용 -- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `numerical-review-agent`; revision `r1`. +- Feature ID: `mitc4-shell-feast-eigen2`; Owner: `numerical-review-agent`; revision `r2` (2026-09-10 첫 구현의 외부 Value 5% 승인). - 기본 보고서: [numerical-review](numerical-review.md); 기준: [requirements](requirements.md), [TOLERANCE](../TOLERANCE.md). -- 판정: 비교 계약 정의 완료. 실제 FESA 결과 comparison은 **NOT RUN**. +- 판정: 승인된 비교 계약 갱신. 실제 최신 판정은 [첫 구현 검토](implementation-first-version-review.md)에 기록한다. ## Eigen의 지정 원본과 logical identity @@ -45,10 +45,20 @@ component `Value`, `blocking`이다. 모든 mode를 한 frame으로 합치는 의미를 삭제하며 낮은 mode의 scale을 높은 mode로 바꾸므로 채택하지 않는다. 이는 과거 singleton 결론을 자동 승계하거나 결과 오차를 본 뒤 고른 family가 아니다. -각 family의 S=abs(reference Value), B=.01S, 행별 relative limit .05 및 relative RMS .01을 -그대로 적용한다. 이번 양의 singleton에서 RMS=relative error이므로 실효 각 Value limit은 -1%다. 그렇더라도 행 판정과 RMS를 각각 기록한다. Inclusive boundary, near-zero, -zero-reference-scale의 exact-zero 및 named failure는 공통 문서 그대로다. +2026-09-10 사용자는 “5% 이내면 첫 구현이므로 Abaqus와의 고유치 정확도 향상은 추후에 +진행하겠습니다. 구현을 마무리해주세요.”라고 승인했다. [Feature ADR](ADR.md)와 +[공통 TOLERANCE의 한정 정책](../TOLERANCE.md)에 따라 `mitc4-eigen-first-implementation` +profile을 이 Feature의 외부 `Value`에만 적용한다. + +각 family의 S=abs(reference Value), B=.01S와 행별 relative limit .05는 유지하며 +relative RMS limit을 .05로 적용한다. 이번 양의 singleton에서 RMS=relative error이므로 +각 Value는 **상대오차 5% 이하(inclusive)**일 때 통과한다. 행 판정과 RMS를 각각 기록하고 +profile ID·실제 RMS 한계 .05를 보고한다. Near-zero, zero-reference-scale의 exact-zero, +named failure와 identity/finite 전검사는 공통 문서 그대로다. 다른 Feature·static family의 +RMS .01, 이 Feature의 전체 c개 내부 검증 1e-10은 변경하지 않는다. + +기존 r1의 singleton RMS1% 실패 기록은 당시 판정으로 보존한다. 이번 승인은 첫 구현의 +외부 정확도 기준 변경이며 S3 해결·다른 필수 검증 면제·Abaqus 내부 알고리즘 동등성 승인이 아니다. Frequency 변환, Freq fitting, mode shapes/변위/회전/stress/reaction은 새 eigen 외부 gate가 아니다. File/dataset/component/schema 누락, wrong step/frame/quantity/backend/mass variant, diff --git a/docs/mitc4-shell-feast-eigen2/requirements.md b/docs/mitc4-shell-feast-eigen2/requirements.md index b27a8ee..e4b248b 100644 --- a/docs/mitc4-shell-feast-eigen2/requirements.md +++ b/docs/mitc4-shell-feast-eigen2/requirements.md @@ -1,7 +1,7 @@ # MITC4 질량행렬·FEAST eigen 요구사항 - Feature ID: `mitc4-shell-feast-eigen2` -- Requirement baseline revision: `r1` +- Requirement baseline revision: `r2` (2026-09-10, 사용자 승인에 따른 022 외부 accuracy 기준 변경) - Owner: `requirement-agent` - 상태: 요구사항 기준선 작성. 구현·수치 검증 또는 readiness 통과 기록이 아니다. @@ -83,7 +83,7 @@ I/O layout은 후속 전문 owner가 필요한 시점에 정의한다. 외부 to | FEM-REQ-MITC4-EIGEN2-019 | State와 final HDF5는 검증된 candidate만 반영하고 제어된 실패에서 기존 유효 결과를 보존해야 한다. | Solve·validation·writer/finalization 실패 후 이전 state/final hash 보존; live log는 기존 `_log.txt` 독립 lifecycle 유지 | P, N; Formulation, I/O, Implementation | | FEM-REQ-MITC4-EIGEN2-020 | 새 eigen 외부 검증은 지정 원본 TXT의 `Value`만 derived lumped mass로 계산한 FESA HDF5 고유치와 비교해야 한다. | 같은 지정 `.inp`·TXT, Mode 1~10 포함; mass variant 사후 선택 금지. Eigen `Freq`·변위·회전·stress·reaction·mode shape 외부 gate를 추가하지 않음. 승인된 기존 static 회귀는 별도 Requirement 024의 기존 계약 사용 | U, R; NR, I/O | | FEM-REQ-MITC4-EIGEN2-021 | 비교기는 required identity의 missing/extra/duplicate/nonfinite와 quantity/schema 오류를 tolerance 전에 거부해야 한다. | 행 재정렬은 같은 mapping, mode label 중복·누락·추가·비유한 `Value`는 실패; 같은 값의 다른 mode는 허용; comparator 음성 테스트 | P, R; NR, I/O, Implementation | -| FEM-REQ-MITC4-EIGEN2-022 | 외부 고유치 acceptance는 공통 `docs/TOLERANCE.md`의 scale·row·RMS·zero-scale 정책을 적용해야 한다. | NR/I/O가 공통 정책의 같은 step/frame 규칙 안에서 mode identity·family를 정의하고 모든 요청 Value의 blocking 비교 통과; frequency로 대체하거나 tolerance를 조정하지 않음 | U, P; NR, I/O | +| FEM-REQ-MITC4-EIGEN2-022 | 첫 구현의 외부 고유치 acceptance는 `docs/TOLERANCE.md`의 승인된 `mitc4-eigen-first-implementation` 예외를 적용해야 한다. | 같은 step/modal frame별 singleton Value의 relative/RMS가 각각5% 이하일 때 통과. Scale·zero-scale·identity/finite 사전검사와 보고는 공통 정책 유지. Policy ID와 RMS5%를 명시하며 frequency 대체·reference 보정 금지. 승인 근거는 Feature ADR-012 | U, P; NR, I/O | | FEM-REQ-MITC4-EIGEN2-023 | 변경 동작은 의미 있는 테스트와 구성된 MSVC x64 Debug CMake/CTest 검증으로 확인해야 한다. | 가능한 버그 재현, 변경 중 관련 검사, 기능 완료 시 전체 build/test와 필수 reference 비교 1회; 새 warning 없음과 실제 FEAST runtime evidence. 아래 제외 예제 사용 0건 및 남은 필수 coverage와 실행 identity 대조; 형식적 RED·수동/자동 중복 실행을 요구하지 않음 | U, P; Implementation | | FEM-REQ-MITC4-EIGEN2-024 | 기존 선형·비선형 정적해석은 승인된 stabilization 개선 외의 입력·출력·물리 의미를 보존하고 영향받는 correctness·reference 회귀를 통과해야 한다. | 아래 portfolio의 `shell`·`shellNL` 및 linear-limit 역할을 fresh solve로 검증하고 기존 comparison identity·판정 정책 유지; tangent/선형극한/평형·reaction/energy/state 검증. 제외 예제는 복사본·cache로도 재사용하지 않음. 개선 전 bitwise equality를 요구하지 않고 density 없는 기존 정적 모델도 유지 | U, P, L, N; Formulation, NR, Implementation, Physics | | FEM-REQ-MITC4-EIGEN2-025 | 솔버 기능 완료 판단은 적용 계약과 모든 필수 검증 evidence를 추적하여 내려야 한다. | 각 요구와 독립 numerical review·analytical/invariant·필요한 physics/reference evidence 연결; 필수 실패·미실행·stale를 통과로 집계하지 않음. 사용자 제외·과거 결과·현재 실행 구분; 기존 명시 감사 요청은 아래 범위 경계대로 보존 | U, P; 메인, NR, Physics | diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp index c2ece52..883a327 100644 --- a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp +++ b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp @@ -17,23 +17,23 @@ #include #include "fesa/math/sparse_matrix.h" +#include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h" namespace fesa { + +namespace internal { namespace { - -constexpr std::uint64_t kSeed = 0x4D495443344D4645ULL; constexpr std::size_t kMaximumCountCalls = 512U; - -Status Failure(const std::string& code, const std::string& identity, - const std::string& message) { - return Status::Failure( - FailureCategory::kSolver, - {{Severity::kError, code, {{}, 0U}, "FEAST", identity, message}}); +thread_local FeastBackendTestHook* g_feast_backend_test_hook = nullptr; +thread_local FeastWorkspaceTestHook* g_feast_workspace_test_hook = nullptr; } -bool FitsMkl(const std::size_t value) { - return value <= - static_cast((std::numeric_limits::max)()); +bool CanCallFeast(const std::size_t used) { + return used < kMaximumCountCalls; +} + +bool CanStartUpperCandidate(const std::size_t used) { + return used <= 384U - 21U; } bool ProductFits(const std::size_t left, const std::size_t right) { @@ -48,6 +48,78 @@ bool FitsDoubleStorage(const std::size_t count) { sizeof(double); } +ScopedFeastBackendTestHook::ScopedFeastBackendTestHook( + FeastBackendTestHook hook) + : hook_(std::move(hook)), previous_(g_feast_backend_test_hook) { + g_feast_backend_test_hook = &hook_; +} + +ScopedFeastBackendTestHook::~ScopedFeastBackendTestHook() { + g_feast_backend_test_hook = previous_; +} + +ScopedFeastWorkspaceTestHook::ScopedFeastWorkspaceTestHook( + FeastWorkspaceTestHook hook) + : hook_(std::move(hook)), previous_(g_feast_workspace_test_hook) { + g_feast_workspace_test_hook = &hook_; +} + +ScopedFeastWorkspaceTestHook::~ScopedFeastWorkspaceTestHook() { + g_feast_workspace_test_hook = previous_; +} + +void ApplyFeastWorkspaceTestHook(const bool count_workspace) { + if (g_feast_workspace_test_hook != nullptr) + (*g_feast_workspace_test_hook)(count_workspace); +} + +void ApplyFeastBackendTestHook(const FeastBackendTestCall& call, MKL_INT* m0, + MKL_INT* returned, MKL_INT* loop, + double* epsout, MKL_INT* info, + double* eigenvalues, + const std::size_t eigenvalue_count, + double* vectors, + const std::size_t vector_count, + double* residuals, + const std::size_t residual_count) { + if (g_feast_backend_test_hook == nullptr) return; + FeastBackendTestOutputs outputs{ + static_cast(*info), static_cast(*m0), + static_cast(*returned), static_cast(*loop), *epsout, + std::vector(eigenvalues, eigenvalues + eigenvalue_count), + std::vector(vectors, vectors + vector_count), + std::vector(residuals, residuals + residual_count)}; + (*g_feast_backend_test_hook)(call, outputs); + *info = static_cast(outputs.info); + *m0 = static_cast(outputs.m0); + *returned = static_cast(outputs.returned); + *loop = static_cast(outputs.loop); + *epsout = outputs.epsout; + if (outputs.eigenvalues.size() == eigenvalue_count) + std::copy(outputs.eigenvalues.begin(), outputs.eigenvalues.end(), eigenvalues); + if (outputs.vectors.size() == vector_count) + std::copy(outputs.vectors.begin(), outputs.vectors.end(), vectors); + if (outputs.residuals.size() == residual_count) + std::copy(outputs.residuals.begin(), outputs.residuals.end(), residuals); +} +} // namespace internal + +namespace { + +constexpr std::uint64_t kSeed = 0x4D495443344D4645ULL; + +Status Failure(const std::string& code, const std::string& identity, + const std::string& message) { + return Status::Failure( + FailureCategory::kSolver, + {{Severity::kError, code, {{}, 0U}, "FEAST", identity, message}}); +} + +bool FitsMkl(const std::size_t value) { + return value <= + static_cast((std::numeric_limits::max)()); +} + double NextUp(const double value) { return std::nextafter(value, (std::numeric_limits::infinity)()); } @@ -307,12 +379,12 @@ Result MklFeastEigenSolver::Solve( "feast-invalid-dimension", "matrix-or-request", "FEAST requires 1 <= requested modes <= square matrix dimension.")); } - if (!FitsMkl(n) || !ProductFits(n, n)) { + if (!FitsMkl(n) || !internal::ProductFits(n, n)) { return Result::Failure(Failure( "feast-integer-overflow", "matrix-shape", "FEAST workspace dimensions exceed the private backend range.")); } - if (!FitsDoubleStorage(n) || !FitsDoubleStorage(n * n)) { + if (!internal::FitsDoubleStorage(n) || !internal::FitsDoubleStorage(n * n)) { return Result::Failure(Failure( "eigen-allocation-failed", "count-workspace", "The required count workspace exceeds addressable double storage.")); @@ -405,6 +477,7 @@ Result MklFeastEigenSolver::Solve( std::vector count_res; std::vector count_x; try { + internal::ApplyFeastWorkspaceTestHook(true); count_e.assign(n, 0.0); count_res.assign(n, 0.0); count_x.assign(n * n, 0.0); @@ -416,7 +489,7 @@ Result MklFeastEigenSolver::Solve( auto Count = [&](const double call_lower, const double call_upper, const std::string& purpose, std::size_t* count) -> Status { - if (calls.size() >= kMaximumCountCalls) + if (!internal::CanCallFeast(calls.size())) return Failure("eigen-search-budget-exhausted", "count-calls", "The bounded FEAST count-call budget was exhausted."); std::array fpm{}; @@ -438,6 +511,10 @@ Result MklFeastEigenSolver::Solve( b.Value().columns.data(), fpm.data(), &epsout, &loop, &call_lower, &call_upper, &m0, count_e.data(), count_x.data(), &returned, count_res.data(), &info); + internal::ApplyFeastBackendTestHook( + {true, calls.size(), call_lower, call_upper, n}, &m0, &returned, + &loop, &epsout, &info, count_e.data(), count_e.size(), count_x.data(), + count_x.size(), count_res.data(), count_res.size()); const bool count_available = info == 0 && m0 >= 0 && static_cast(m0) <= n; FeastBackendCall event{calls.size(), @@ -538,7 +615,7 @@ Result MklFeastEigenSolver::Solve( for (std::size_t midpoint_attempt = 0U; upper_search_termination != "count-target" && midpoint_attempt < 32U; ++midpoint_attempt) { - if (calls.size() + 21U > 384U) { + if (!internal::CanStartUpperCandidate(calls.size())) { upper_search_termination = "optimization-budget"; break; } @@ -646,11 +723,12 @@ Result MklFeastEigenSolver::Solve( const std::size_t padded_count = (std::max)(count + 8U, (3U * count + 1U) / 2U); std::size_t capacity = count == n ? n : (std::min)(n, padded_count); - if (!ProductFits(n, capacity)) + if (!internal::ProductFits(n, capacity)) return Result::Failure( Failure("feast-workspace-overflow", "solve-workspace", "The FEAST solve workspace size overflows.")); - if (!FitsDoubleStorage(n * capacity) || !FitsDoubleStorage(capacity)) { + if (!internal::FitsDoubleStorage(n * capacity) || + !internal::FitsDoubleStorage(capacity)) { return Result::Failure(Failure( "eigen-allocation-failed", "solve-workspace", "The required solve workspace exceeds addressable double storage.")); @@ -663,11 +741,12 @@ Result MklFeastEigenSolver::Solve( bool solve_accepted = false; std::size_t solve_capacity = capacity; for (std::size_t attempt = 0U; attempt < 3U; ++attempt) { - if (calls.size() >= kMaximumCountCalls) + if (!internal::CanCallFeast(calls.size())) return Result::Failure( Failure("eigen-search-budget-exhausted", "solve-call", "The total backend call budget is exhausted before solve.")); try { + internal::ApplyFeastWorkspaceTestHook(false); eigenvalues.assign(solve_capacity, 0.0); residuals.assign(solve_capacity, 0.0); vectors.assign(n * solve_capacity, 0.0); @@ -695,6 +774,11 @@ Result MklFeastEigenSolver::Solve( b.Value().columns.data(), fpm.data(), &epsout, &loop, &fixed_lower, &solve_upper, &m0, eigenvalues.data(), vectors.data(), &returned, residuals.data(), &info); + internal::ApplyFeastBackendTestHook( + {false, calls.size(), fixed_lower, solve_upper, solve_capacity}, &m0, + &returned, &loop, &epsout, &info, eigenvalues.data(), + eigenvalues.size(), vectors.data(), vectors.size(), residuals.data(), + residuals.size()); const bool solve_outputs_available = info == 0 && loop >= 0 && std::isfinite(epsout) && returned >= 0 && static_cast(returned) <= solve_capacity; diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver_test_seam.h b/src/fesa/solvers/eigen/mkl_feast_eigen_solver_test_seam.h new file mode 100644 index 0000000..dcfbd09 --- /dev/null +++ b/src/fesa/solvers/eigen/mkl_feast_eigen_solver_test_seam.h @@ -0,0 +1,70 @@ +#ifndef FESA_SOLVERS_EIGEN_MKL_FEAST_EIGEN_SOLVER_TEST_SEAM_H_ +#define FESA_SOLVERS_EIGEN_MKL_FEAST_EIGEN_SOLVER_TEST_SEAM_H_ + +#include +#include +#include + +namespace fesa::internal { + +bool CanCallFeast(std::size_t used); +bool CanStartUpperCandidate(std::size_t used); +bool ProductFits(std::size_t left, std::size_t right); +bool FitsDoubleStorage(std::size_t count); + +// Private test seam. It is applied only after the production oneMKL call has +// returned and is thread-local so production solver state remains unchanged. +struct FeastBackendTestCall { + bool count; + std::size_t ordinal; + double lower; + double upper; + std::size_t capacity; +}; + +struct FeastBackendTestOutputs { + int info; + std::ptrdiff_t m0; + std::ptrdiff_t returned; + int loop; + double epsout; + std::vector eigenvalues; + std::vector vectors; + std::vector residuals; +}; + +using FeastBackendTestHook = + std::function; +using FeastWorkspaceTestHook = std::function; + +class ScopedFeastBackendTestHook final { + public: + explicit ScopedFeastBackendTestHook(FeastBackendTestHook hook); + ~ScopedFeastBackendTestHook(); + + ScopedFeastBackendTestHook(const ScopedFeastBackendTestHook&) = delete; + ScopedFeastBackendTestHook& operator=(const ScopedFeastBackendTestHook&) = + delete; + + private: + FeastBackendTestHook hook_; + FeastBackendTestHook* previous_; +}; + +class ScopedFeastWorkspaceTestHook final { + public: + explicit ScopedFeastWorkspaceTestHook(FeastWorkspaceTestHook hook); + ~ScopedFeastWorkspaceTestHook(); + + ScopedFeastWorkspaceTestHook(const ScopedFeastWorkspaceTestHook&) = delete; + ScopedFeastWorkspaceTestHook& operator=( + const ScopedFeastWorkspaceTestHook&) = delete; + + private: + FeastWorkspaceTestHook hook_; + FeastWorkspaceTestHook* previous_; +}; + +} // namespace fesa::internal + +#endif // FESA_SOLVERS_EIGEN_MKL_FEAST_EIGEN_SOLVER_TEST_SEAM_H_ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index dde1cc8..4a08197 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -7,6 +7,7 @@ add_executable( unit/analysis/analysis_model_test.cpp unit/analysis/analysis_state_test.cpp unit/analysis/eigen_reduction_test.cpp + unit/analysis/eigen_mass_assembly_test.cpp unit/analysis/nonlinear_analysis_state_test.cpp unit/analysis/nonlinear_static_control_test.cpp unit/analysis/newton_iteration_log_test.cpp diff --git a/tests/reference/eigen_reference_comparison.cpp b/tests/reference/eigen_reference_comparison.cpp index d1e6209..6a80dd5 100644 --- a/tests/reference/eigen_reference_comparison.cpp +++ b/tests/reference/eigen_reference_comparison.cpp @@ -388,6 +388,9 @@ Result EigenReferenceComparison::Compare( report.case_id = reference_case.case_id; report.source_step_name = frequency.source_step_name_raw; report.requested_mode_count = requested; + report.policy_id = + ReferenceTolerancePolicy::kMitc4EigenFirstImplementationPolicyId; + report.relative_rms_tolerance = ReferenceTolerancePolicy::kRelativeTolerance; std::vector families; families.reserve(modes.size()); for (const auto& mode : modes) { @@ -404,7 +407,8 @@ Result EigenReferenceComparison::Compare( true, {{mode.eigenvalue, row->second}}}); } - const auto evaluation = ReferenceTolerancePolicy::Evaluate(families); + const auto evaluation = + ReferenceTolerancePolicy::EvaluateMitc4EigenFirstImplementation(families); if (!evaluation.valid || evaluation.families.size() != families.size()) { Fail("nonfinite-comparison-value", "Tolerance evaluation rejected modal values."); diff --git a/tests/reference/eigen_reference_comparison.h b/tests/reference/eigen_reference_comparison.h index be4be6e..cb6c002 100644 --- a/tests/reference/eigen_reference_comparison.h +++ b/tests/reference/eigen_reference_comparison.h @@ -34,6 +34,8 @@ struct EigenReferenceComparisonReport { std::vector rows; std::vector families; bool passed; + std::string policy_id; + double relative_rms_tolerance; }; /// @brief Compares modal eigenvalues without changing the approved input or diff --git a/tests/reference/eigen_reference_comparison_test.cpp b/tests/reference/eigen_reference_comparison_test.cpp index e61c479..8b3dc6b 100644 --- a/tests/reference/eigen_reference_comparison_test.cpp +++ b/tests/reference/eigen_reference_comparison_test.cpp @@ -5,9 +5,12 @@ #include #include +#include #include #include +#include #include +#include #include #include "fesa/io/abaqus/input_reader.h" @@ -117,7 +120,32 @@ TEST(EigenReferenceComparison, auto reference_case = fesa::test::EigenReferenceCase{"x", inp, txt, h5}; auto ok = fesa::test::EigenReferenceComparison::Compare(reference_case); ASSERT_TRUE(ok.HasValue()); + EXPECT_EQ(ok.Value().policy_id, "mitc4-eigen-first-implementation"); + EXPECT_DOUBLE_EQ(ok.Value().relative_rms_tolerance, 0.05); EXPECT_EQ(ok.Value().families.size(), 2U); + // Only synthetic FESA output is changed; each TXT mode remains an independent + // scalar family, including at both inclusive five-percent boundaries. + for (const auto& sample : + std::vector>{ + {103.0, true}, {105.0, true}, {95.0, true}, + {std::nextafter(105.0, std::numeric_limits::infinity()), false}, + {std::nextafter(95.0, 0.0), false}}) { + SCOPED_TRACE(sample.first); + f = H5Fopen(h5.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + ds = H5Dopen2(f, "/steps/Step-1/eigen/modes", H5P_DEFAULT); + t = H5Dget_type(ds); + v[0].eigenvalue = sample.first; + ASSERT_GE(H5Dwrite(ds, t, H5S_ALL, H5S_ALL, H5P_DEFAULT, v), 0); + H5Tclose(t); + H5Dclose(ds); + H5Fclose(f); + const auto comparison = + fesa::test::EigenReferenceComparison::Compare(reference_case); + ASSERT_TRUE(comparison.HasValue()); + EXPECT_EQ(comparison.Value().passed, sample.second); + EXPECT_EQ(comparison.Value().families[0].rms_passed, sample.second); + EXPECT_TRUE(comparison.Value().families[1].passed); + } struct TextMutation { std::string text; std::string expected_code; diff --git a/tests/reference/reference_tolerance_policy.cpp b/tests/reference/reference_tolerance_policy.cpp index e6b750d..cecdc51 100644 --- a/tests/reference/reference_tolerance_policy.cpp +++ b/tests/reference/reference_tolerance_policy.cpp @@ -19,10 +19,9 @@ bool IsFiniteFamily(const ReferenceToleranceFamily& family) { }); } -} - -ReferenceToleranceEvaluation ReferenceTolerancePolicy::Evaluate( - const std::vector& families) { +ReferenceToleranceEvaluation EvaluateFamilies( + const std::vector& families, + const double relative_rms_tolerance) { ReferenceToleranceEvaluation evaluation{true, {}, true, {}}; evaluation.families.reserve(families.size()); @@ -41,7 +40,8 @@ ReferenceToleranceEvaluation ReferenceTolerancePolicy::Evaluate( decision.reference_scale = (std::max)(decision.reference_scale, std::abs(row.reference_value)); } - decision.near_zero_band = kNearZeroRatio * decision.reference_scale; + decision.near_zero_band = + ReferenceTolerancePolicy::kNearZeroRatio * decision.reference_scale; double error_norm = 0.0; for (std::size_t index = 0U; index < family.rows.size(); ++index) { @@ -74,12 +74,14 @@ ReferenceToleranceEvaluation ReferenceTolerancePolicy::Evaluate( row_decision.passed = absolute_error <= row_decision.threshold; } else { row_decision.threshold = - kRelativeTolerance * std::abs(row.reference_value); + ReferenceTolerancePolicy::kRelativeTolerance * + std::abs(row.reference_value); row_decision.relative_error = absolute_error / std::abs(row.reference_value); row_decision.relative_error_applicable = true; row_decision.branch = ReferenceToleranceBranch::kRelative; - row_decision.passed = row_decision.relative_error <= kRelativeTolerance; + row_decision.passed = row_decision.relative_error <= + ReferenceTolerancePolicy::kRelativeTolerance; } decision.rows.push_back(row_decision); } @@ -100,7 +102,7 @@ ReferenceToleranceEvaluation ReferenceTolerancePolicy::Evaluate( if (!std::isfinite(decision.relative_rms)) { return {false, "nonfinite-comparison-value", false, {}}; } - decision.rms_passed = decision.relative_rms <= kRelativeRmsTolerance; + decision.rms_passed = decision.relative_rms <= relative_rms_tolerance; } decision.passed = rows_passed && decision.rms_passed; if (decision.blocking && !decision.passed) { @@ -112,4 +114,25 @@ ReferenceToleranceEvaluation ReferenceTolerancePolicy::Evaluate( return evaluation; } +} // namespace + +ReferenceToleranceEvaluation ReferenceTolerancePolicy::Evaluate( + const std::vector& families) { + return EvaluateFamilies(families, kRelativeRmsTolerance); +} + +ReferenceToleranceEvaluation +ReferenceTolerancePolicy::EvaluateMitc4EigenFirstImplementation( + const std::vector& families) { + if (families.empty() || + !std::all_of(families.begin(), families.end(), [](const auto& family) { + return family.blocking && family.rows.size() == 1U && + family.components.size() == 1U && + family.components.front() == "Value"; + })) { + return {false, "invalid-eigen-tolerance-family", false, {}}; + } + return EvaluateFamilies(families, kRelativeTolerance); +} + } diff --git a/tests/reference/reference_tolerance_policy.h b/tests/reference/reference_tolerance_policy.h index 5cda054..5d41d19 100644 --- a/tests/reference/reference_tolerance_policy.h +++ b/tests/reference/reference_tolerance_policy.h @@ -63,9 +63,16 @@ class ReferenceTolerancePolicy { static constexpr double kNearZeroRatio = 0.01; static constexpr double kRelativeTolerance = 0.05; static constexpr double kRelativeRmsTolerance = 0.01; + static constexpr const char* kMitc4EigenFirstImplementationPolicyId = + "mitc4-eigen-first-implementation"; static ReferenceToleranceEvaluation Evaluate( const std::vector& families); + + /// @brief Applies the approved five-percent RMS limit to blocking singleton + /// Value families from the MITC4 eigen first implementation only. + static ReferenceToleranceEvaluation EvaluateMitc4EigenFirstImplementation( + const std::vector& families); }; } diff --git a/tests/reference/reference_tolerance_policy_test.cpp b/tests/reference/reference_tolerance_policy_test.cpp index 90335d0..bb56f05 100644 --- a/tests/reference/reference_tolerance_policy_test.cpp +++ b/tests/reference/reference_tolerance_policy_test.cpp @@ -115,4 +115,59 @@ TEST(ReferenceTolerancePolicy, EXPECT_TRUE(evaluation.families.empty()); } +TEST(ReferenceTolerancePolicy, + EigenFirstImplementationDoesNotRelaxCommonPolicy) { + const std::vector families{ + MakeFamily("mode-1", {"Value"}, true, {{105.0, 100.0}})}; + const auto common = ReferenceTolerancePolicy::Evaluate(families); + const auto eigen = + ReferenceTolerancePolicy::EvaluateMitc4EigenFirstImplementation(families); + ASSERT_TRUE(common.valid); + ASSERT_TRUE(eigen.valid); + EXPECT_FALSE(common.passed); + EXPECT_TRUE(eigen.passed); + EXPECT_DOUBLE_EQ(eigen.families[0].relative_rms, 0.05); + EXPECT_DOUBLE_EQ(common.families[0].relative_rms, + eigen.families[0].relative_rms); +} + +TEST(ReferenceTolerancePolicy, + EigenFirstImplementationPreservesZeroAndFiniteChecks) { + for (const auto& row : std::vector{ + {0.0, 0.0}, {1.0, 0.0}, + {std::numeric_limits::quiet_NaN(), 1.0}}) { + const auto result = + ReferenceTolerancePolicy::EvaluateMitc4EigenFirstImplementation( + {MakeFamily("mode-1", {"Value"}, true, {row})}); + if (!std::isfinite(row.fesa_value)) { + EXPECT_FALSE(result.valid); + EXPECT_FALSE(result.passed); + } else { + ASSERT_TRUE(result.valid); + EXPECT_EQ(result.passed, row.fesa_value == 0.0); + if (row.fesa_value != 0.0) { + EXPECT_EQ(result.families[0].diagnostic_code, + "zero-reference-scale-nonzero-error"); + } + } + } +} + +TEST(ReferenceTolerancePolicy, + EigenFirstImplementationRejectsNonmodalFamilies) { + for (const auto& families : + std::vector>{ + {}, + {MakeFamily("mode-1", {"U1"}, true, {{100.0, 100.0}})}, + {MakeFamily("mode-1", {"Value"}, false, {{100.0, 100.0}})}, + {MakeFamily("modes", {"Value"}, true, + {{100.0, 100.0}, {100.0, 100.0}})}}) { + const auto result = + ReferenceTolerancePolicy::EvaluateMitc4EigenFirstImplementation(families); + EXPECT_FALSE(result.valid); + EXPECT_FALSE(result.passed); + EXPECT_EQ(result.diagnostic_code, "invalid-eigen-tolerance-family"); + } +} + } diff --git a/tests/unit/analysis/eigen_mass_assembly_test.cpp b/tests/unit/analysis/eigen_mass_assembly_test.cpp new file mode 100644 index 0000000..d8f580f --- /dev/null +++ b/tests/unit/analysis/eigen_mass_assembly_test.cpp @@ -0,0 +1,180 @@ +#include "fesa/analysis/linear_frequency_procedure.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/prepared_analysis.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/mitc4_shell.h" +#include "fesa/model/domain.h" + +namespace { + +std::string StatusText(const fesa::Status& status) { + std::string result; + for (const auto& diagnostic : status.Diagnostics()) + result += diagnostic.code + " " + diagnostic.message + "\n"; + return result; +} + +fesa::Result SharedNodeDomain(const bool reverse_source_order) { + std::ostringstream deck; + deck << R"inp(*Part,name=Panel +*Node +1,0,0,0 +2,1,0,0 +3,2,0,0 +4,0,1,0 +5,1,1,0 +6,2,1,0 +*Element,type=S4 +)inp"; + if (reverse_source_order) + deck << "20,2,3,6,5\n10,1,2,5,4\n"; + else + deck << "10,1,2,5,4\n20,2,3,6,5\n"; + deck << R"inp(*Elset,elset=Thin +10 +*Elset,elset=Thick +20 +*Shell Section,elset=Thin,material=Light +0.5 +*Shell Section,elset=Thick,material=Heavy +1.0 +*End Part +*Assembly,name=Assembly +*Instance,name=Panel-1,part=Panel +*End Instance +*Nset,nset=Fixed,instance=Panel-1 +1,2,3,4,6 +*Nset,nset=Shared,instance=Panel-1 +5 +*End Assembly +*Material,name=Light +*Elastic +12,0 +*Density +2 +*Material,name=Heavy +*Elastic +12,0 +*Density +5 +*Boundary +Fixed,1,6 +Shared,1,3 +Shared,5,6 +*Step,name=Step-2,nlgeom=NO,perturbation +*Frequency +1 +*End Step +)inp"; + static std::atomic sequence{0U}; + const auto path = std::filesystem::temp_directory_path() / + ("fesa-eigen-mass-" + + std::to_string(std::chrono::steady_clock::now() + .time_since_epoch() + .count()) + + "-" + std::to_string(sequence.fetch_add(1U)) + ".inp"); + { + std::ofstream stream{path}; + stream << deck.str(); + } + auto prepared = fesa::PrepareAnalysisInput(path); + std::error_code ignored; + std::filesystem::remove(path, ignored); + if (!prepared.HasValue()) + return fesa::Result::Failure(prepared.GetStatus()); + return fesa::Result::Success( + std::move(prepared.Value().domain)); +} + +struct SharedDofPencil { + double stiffness{0.0}; + double mass{0.0}; +}; + +fesa::Result DirectSharedDofPencil( + const fesa::Domain& domain) { + const fesa::ElementFactory factory; + SharedDofPencil result; + for (std::size_t index = 0U; index < domain.ShellElements().Size(); ++index) { + const auto& definition = domain.ShellElements()[index]; + const auto runtime = factory.Create(definition, domain); + if (!runtime.HasValue()) + return fesa::Result::Failure(runtime.GetStatus()); + const auto* shell = + dynamic_cast(runtime.Value().get()); + if (shell == nullptr) + return fesa::Result::Failure(fesa::Status::Failure( + fesa::FailureCategory::kModel, + {{fesa::Severity::kError, "invalid-test-shell", {}, "ELEMENT", {}, + "Frequency mass fixture did not create an MITC4 shell."}})); + const auto local = std::find(definition.NodeIndices().begin(), + definition.NodeIndices().end(), 4U); + if (local == definition.NodeIndices().end()) + return fesa::Result::Failure(fesa::Status::Failure( + fesa::FailureCategory::kModel, + {{fesa::Severity::kError, "invalid-test-shell", {}, "ELEMENT", {}, + "Frequency mass fixture lost its shared node."}})); + const auto mass = + shell->Mass(definition.MaterialIndex() == 0U ? 2.0 : 5.0); + const auto stiffness = shell->Stiffness(); + if (!mass.HasValue()) + return fesa::Result::Failure(mass.GetStatus()); + if (!stiffness.HasValue()) + return fesa::Result::Failure(stiffness.GetStatus()); + const std::size_t dof = + 6U * static_cast( + std::distance(definition.NodeIndices().begin(), local)) + + 3U; + result.mass += mass.Value().lumped_global24(dof, dof); + result.stiffness += stiffness.Value().stabilized_global24(dof, dof); + } + return fesa::Result::Success(result); +} + +// E2-C06: sections and densities differ across the shared source node. The +// one free physical rotation makes the executor's reduced pencil exactly the +// independently scattered local diagonal sum. +TEST(EigenMassAssembly, ScattersHeterogeneousSharedNodeMassInSourceOrder) { + constexpr double analytic_shared_rotary_mass = 7.0 / 64.0; + std::array eigenvalues{}; + for (std::size_t order = 0U; order < eigenvalues.size(); ++order) { + const auto domain = SharedNodeDomain(order != 0U); + ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus()); + const auto expected = DirectSharedDofPencil(domain.Value()); + ASSERT_TRUE(expected.HasValue()) << StatusText(expected.GetStatus()); + EXPECT_LE(std::abs(expected.Value().mass - analytic_shared_rotary_mass) / + analytic_shared_rotary_mass, + 1e-12); + const auto state = + fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(state.HasValue()) << StatusText(state.GetStatus()); + ASSERT_EQ(state.Value().free_dof_count, 1U); + ASSERT_EQ(state.Value().modes.size(), 1U); + const double expected_eigenvalue = + expected.Value().stiffness / analytic_shared_rotary_mass; + EXPECT_LE(std::abs(state.Value().modes[0U].eigenvalue - expected_eigenvalue) / + (std::abs(state.Value().modes[0U].eigenvalue) + + std::abs(expected_eigenvalue)), + 1e-12); + const double amplitude = state.Value().vectors(0U, 6U * 4U + 3U); + EXPECT_LE(std::abs(amplitude * amplitude * analytic_shared_rotary_mass - + 1.0), + 1e-12); + eigenvalues[order] = state.Value().modes[0U].eigenvalue; + } + EXPECT_NEAR(eigenvalues[0U], eigenvalues[1U], 1e-12); +} + +} // namespace diff --git a/tests/unit/elements/mitc4_shell_mass_contract_test.cpp b/tests/unit/elements/mitc4_shell_mass_contract_test.cpp index c847150..d0697f9 100644 --- a/tests/unit/elements/mitc4_shell_mass_contract_test.cpp +++ b/tests/unit/elements/mitc4_shell_mass_contract_test.cpp @@ -1,5 +1,6 @@ #include +#include #include #include #include @@ -12,6 +13,90 @@ namespace { using Coordinates = std::array; +double FrobeniusNorm(const fesa::Matrix& matrix, const double length) { + double norm = 0.0; + for (std::size_t row = 0U; row < matrix.Rows(); ++row) { + for (std::size_t column = 0U; column < matrix.Columns(); ++column) { + const double row_scale = row % 6U < 3U ? length : 1.0; + const double column_scale = column % 6U < 3U ? length : 1.0; + norm = std::hypot(norm, row_scale * matrix(row, column) * column_scale); + } + } + return norm; +} + +double MaximumEdgeLength(const std::array& nodes) { + double result = 0.0; + for (std::size_t node = 0U; node < nodes.size(); ++node) { + const auto& first = nodes[node].coordinates; + const auto& second = nodes[(node + 1U) % nodes.size()].coordinates; + result = (std::max)(result, std::hypot( + std::hypot(first[0U] - second[0U], first[1U] - second[1U]), + first[2U] - second[2U])); + } + return result; +} + +Coordinates Rotate(const Coordinates& value) { + return {-value[1U], + 0.6 * value[0U] - 0.8 * value[2U], + 0.8 * value[0U] + 0.6 * value[2U]}; +} + +Coordinates DirectorCross(const Coordinates& director, + const Coordinates& rotation) { + return {rotation[1U] * director[2U] - rotation[2U] * director[1U], + rotation[2U] * director[0U] - rotation[0U] * director[2U], + rotation[0U] * director[1U] - rotation[1U] * director[0U]}; +} + +double Dot(const Coordinates& left, const Coordinates& right) { + return left[0U] * right[0U] + left[1U] * right[1U] + + left[2U] * right[2U]; +} + +fesa::Matrix RotateCongruently(const fesa::Matrix& original) { + fesa::Matrix result{24U, 24U}; + for (std::size_t row = 0U; row < 24U; ++row) { + for (std::size_t column = 0U; column < 24U; ++column) { + const std::size_t row_base = 6U * (row / 6U) + + (row % 6U < 3U ? 0U : 3U); + const std::size_t column_base = 6U * (column / 6U) + + (column % 6U < 3U ? 0U : 3U); + const std::size_t row_component = row % 3U; + const std::size_t column_component = column % 3U; + for (std::size_t left = 0U; left < 3U; ++left) { + for (std::size_t right = 0U; right < 3U; ++right) { + const Coordinates basis_left{ + left == 0U ? 1.0 : 0.0, left == 1U ? 1.0 : 0.0, + left == 2U ? 1.0 : 0.0}; + const Coordinates basis_right{ + right == 0U ? 1.0 : 0.0, right == 1U ? 1.0 : 0.0, + right == 2U ? 1.0 : 0.0}; + result(row, column) += Rotate(basis_left)[row_component] * + original(row_base + left, column_base + right) * + Rotate(basis_right)[column_component]; + } + } + } + } + return result; +} + +void ExpectScaledCovariance(const fesa::Matrix& original, + const fesa::Matrix& rotated, + const double length) { + const auto expected = RotateCongruently(original); + fesa::Matrix difference{24U, 24U}; + for (std::size_t row = 0U; row < 24U; ++row) + for (std::size_t column = 0U; column < 24U; ++column) + difference(row, column) = rotated(row, column) - expected(row, column); + const double ratio = FrobeniusNorm(difference, length) / + (FrobeniusNorm(rotated, length) + + FrobeniusNorm(expected, length)); + EXPECT_LE(ratio, 1e-12); +} + fesa::Node Node(std::int64_t label, Coordinates coordinates) { return {{"I", label, std::to_string(label)}, coordinates, {"mass.inp", 1U}}; } @@ -42,23 +127,119 @@ double Quadratic(const fesa::Matrix& matrix, const fesa::Vector& vector) { return vector.Dot(matrix.Multiply(vector)); } -// E2-C02: values derive independently in numerical-review-evidence.md. +// E2-C02: all values derive independently in numerical-review-evidence.md. TEST(Mitc4ShellMassContract, PreservesWarpedMomentsRowsAndFirstSpatialMoment) { const auto mass = WarpedShell().Mass(3.0); ASSERT_TRUE(mass.HasValue()); const auto& value = mass.Value(); - EXPECT_NEAR(value.consistent_moments[0U](0U, 0U), 17609.0 / 28800.0, 1e-12); - EXPECT_NEAR(value.consistent_moments[1U](0U, 0U), 11003.0 / 1440000.0, 1e-12); - EXPECT_NEAR(value.consistent_moments[2U](0U, 0U), 88177.0 / 6912000.0, 1e-12); - double total_m = 0.0, total_h = 0.0, total_j = 0.0; - for (const auto& row : value.nodal_moments) { - total_m += row[0U]; - total_h += row[1U]; - total_j += row[2U]; + const double expected[3][4][4] = { + {{17609.0 / 28800.0, 3337.0 / 11520.0, 1547.0 / 11520.0, + 37.0 / 128.0}, + {3337.0 / 11520.0, 15793.0 / 28800.0, 7177.0 / 28800.0, + 1547.0 / 11520.0}, + {1547.0 / 11520.0, 7177.0 / 28800.0, 4297.0 / 9600.0, + 14239.0 / 57600.0}, + {37.0 / 128.0, 1547.0 / 11520.0, 14239.0 / 57600.0, + 15683.0 / 28800.0}}, + {{11003.0 / 1440000.0, 5053.0 / 1440000.0, 9839.0 / 5760000.0, + 21667.0 / 5760000.0}, + {5053.0 / 1440000.0, 27769.0 / 4320000.0, 53411.0 / 17280000.0, + 9839.0 / 5760000.0}, + {9839.0 / 5760000.0, 53411.0 / 17280000.0, 5687.0 / 960000.0, + 28673.0 / 8640000.0}, + {21667.0 / 5760000.0, 9839.0 / 5760000.0, 28673.0 / 8640000.0, + 63913.0 / 8640000.0}}, + {{88177.0 / 6912000.0, 27847.0 / 4608000.0, 7747.0 / 2764800.0, + 41693.0 / 6912000.0}, + {27847.0 / 4608000.0, 1757.0 / 153600.0, 1331.0 / 256000.0, + 7747.0 / 2764800.0}, + {7747.0 / 2764800.0, 1331.0 / 256000.0, 21521.0 / 2304000.0, + 23773.0 / 4608000.0}, + {41693.0 / 6912000.0, 7747.0 / 2764800.0, 23773.0 / 4608000.0, + 5237.0 / 460800.0}}}; + constexpr double element_mass = 7743.0 / 1600.0; + constexpr double half_thickness = 0.25; + for (std::size_t moment = 0U; moment < 3U; ++moment) { + const double scale = element_mass * std::pow(half_thickness, moment); + for (std::size_t row = 0U; row < 4U; ++row) { + double expected_row = 0.0; + for (std::size_t column = 0U; column < 4U; ++column) { + EXPECT_LE(std::abs(value.consistent_moments[moment](row, column) - + expected[moment][row][column]) / + scale, + 1e-12) + << "moment=" << moment << " row=" << row + << " column=" << column; + expected_row += expected[moment][row][column]; + } + EXPECT_LE(std::abs(value.nodal_moments[row][moment] - expected_row) / + scale, + 1e-12) + << "moment=" << moment << " row=" << row; + } } - EXPECT_NEAR(total_m, 7743.0 / 1600.0, 1e-12); - EXPECT_NEAR(total_h, 739.0 / 12000.0, 1e-12); - EXPECT_NEAR(total_j, 517.0 / 5120.0, 1e-12); + const std::array coordinates{{ + {{-1.0, -1.0, 0.1}}, {{1.0, -1.0, -0.1}}, + {{1.0, 1.0, 0.2}}, {{-1.0, 1.0, -0.1}}}}; + const std::array directors{{ + {{0.0, 0.0, 1.0}}, {{0.6, 0.0, 0.8}}, + {{2.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0}}, {{0.0, 0.6, 0.8}}}}; + const Coordinates reference = coordinates.front(); + const Coordinates expected_first{{-79747.0 / 360000.0, + -685951.0 / 2880000.0, + 49651.0 / 320000.0}}; + Coordinates actual_first{}; + for (std::size_t node = 0U; node < 4U; ++node) { + for (std::size_t component = 0U; component < 3U; ++component) { + actual_first[component] += + value.nodal_moments[node][0U] * + (coordinates[node][component] - reference[component]) + + value.nodal_moments[node][1U] * directors[node][component]; + } + } + double difference = 0.0; + for (std::size_t component = 0U; component < 3U; ++component) + difference = std::hypot(difference, actual_first[component] - + (expected_first[component] - + element_mass * reference[component])); + EXPECT_LE(difference / (element_mass * std::hypot(2.0, 0.3)), 1e-12); +} + +// E2-C05: rebuild both masses after a proper rotation and translation. +TEST(Mitc4ShellMassContract, RecomputesWarpedMassWithScaledCovariance) { + const std::array nodes{ + Node(1, {-1.0, -1.0, 0.1}), Node(2, {1.0, -1.0, -0.1}), + Node(3, {1.0, 1.0, 0.2}), Node(4, {-1.0, 1.0, -0.1})}; + const std::array directors{{ + {{0.0, 0.0, 1.0}}, {{0.6, 0.0, 0.8}}, + {{2.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0}}, {{0.0, 0.6, 0.8}}}}; + std::array transformed_nodes = nodes; + std::array transformed_directors{}; + for (std::size_t node = 0U; node < nodes.size(); ++node) { + transformed_nodes[node].coordinates = Rotate(nodes[node].coordinates); + transformed_nodes[node].coordinates[0U] += 7.0; + transformed_nodes[node].coordinates[1U] -= 3.0; + transformed_nodes[node].coordinates[2U] += 5.0; + transformed_directors[node] = Rotate(directors[node]); + } + const fesa::ShellSection section{"S", 0.5, 0U, {"mass.inp", 2U}}; + const fesa::LinearElasticMaterial material{"M", 12.0, 0.0, + {"mass.inp", 3U}}; + const auto original = fesa::Mitc4Shell::Create(Pointers(nodes), directors, + section, material); + const auto transformed = fesa::Mitc4Shell::Create( + Pointers(transformed_nodes), transformed_directors, section, material); + ASSERT_TRUE(original.HasValue()); + ASSERT_TRUE(transformed.HasValue()); + const auto original_mass = original.Value().Mass(3.0); + const auto transformed_mass = transformed.Value().Mass(3.0); + ASSERT_TRUE(original_mass.HasValue()); + ASSERT_TRUE(transformed_mass.HasValue()); + const double length = MaximumEdgeLength(nodes); + ExpectScaledCovariance(original_mass.Value().consistent_global24, + transformed_mass.Value().consistent_global24, length); + ExpectScaledCovariance(original_mass.Value().lumped_global24, + transformed_mass.Value().lumped_global24, length); } // E2-C03: scalar row-sum lumping retains uniform translation energy and is @@ -80,6 +261,142 @@ TEST(Mitc4ShellMassContract, SatisfiesUniformAndJensenEnergyContracts) { Quadratic(mass.Value().consistent_global24, field)); } +// E2-C03: construct both 24-by-24 mass variants from the independently +// checked scalar moments, including the C1 translation/rotation coupling. +TEST(Mitc4ShellMassContract, MapsScalarMomentsIntoBothMassBlocks) { + const auto mass = WarpedShell().Mass(3.0); + ASSERT_TRUE(mass.HasValue()); + const std::array directors{{ + {{0.0, 0.0, 1.0}}, {{0.6, 0.0, 0.8}}, + {{2.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0}}, {{0.0, 0.6, 0.8}}}}; + fesa::Matrix expected_consistent{24U, 24U}; + fesa::Matrix expected_lumped{24U, 24U}; + for (std::size_t first = 0U; first < 4U; ++first) { + for (std::size_t second = 0U; second < 4U; ++second) { + for (std::size_t row = 0U; row < 3U; ++row) { + for (std::size_t column = 0U; column < 3U; ++column) { + const double identity = row == column ? 1.0 : 0.0; + const Coordinates unit_column{{column == 0U ? 1.0 : 0.0, + column == 1U ? 1.0 : 0.0, + column == 2U ? 1.0 : 0.0}}; + const Coordinates unit_row{{row == 0U ? 1.0 : 0.0, + row == 1U ? 1.0 : 0.0, + row == 2U ? 1.0 : 0.0}}; + const double sj = DirectorCross(directors[second], unit_column)[row]; + const double si = DirectorCross(directors[first], unit_row)[column]; + expected_consistent(6U * first + row, 6U * second + column) = + mass.Value().consistent_moments[0U](first, second) * identity; + expected_consistent(6U * first + row, 6U * second + 3U + column) = + mass.Value().consistent_moments[1U](first, second) * sj; + expected_consistent(6U * first + 3U + row, 6U * second + column) = + mass.Value().consistent_moments[1U](first, second) * si; + expected_consistent(6U * first + 3U + row, + 6U * second + 3U + column) = + mass.Value().consistent_moments[2U](first, second) * + Dot(DirectorCross(directors[first], unit_row), + DirectorCross(directors[second], unit_column)); + } + } + } + for (std::size_t row = 0U; row < 3U; ++row) { + for (std::size_t column = 0U; column < 3U; ++column) { + const double identity = row == column ? 1.0 : 0.0; + const Coordinates unit_column{{column == 0U ? 1.0 : 0.0, + column == 1U ? 1.0 : 0.0, + column == 2U ? 1.0 : 0.0}}; + const Coordinates unit_row{{row == 0U ? 1.0 : 0.0, + row == 1U ? 1.0 : 0.0, + row == 2U ? 1.0 : 0.0}}; + const double sj = DirectorCross(directors[first], unit_column)[row]; + const double si = DirectorCross(directors[first], unit_row)[column]; + expected_lumped(6U * first + row, 6U * first + column) = + mass.Value().nodal_moments[first][0U] * identity; + expected_lumped(6U * first + row, 6U * first + 3U + column) = + mass.Value().nodal_moments[first][1U] * sj; + expected_lumped(6U * first + 3U + row, 6U * first + column) = + mass.Value().nodal_moments[first][1U] * si; + expected_lumped(6U * first + 3U + row, 6U * first + 3U + column) = + mass.Value().nodal_moments[first][2U] * + Dot(DirectorCross(directors[first], unit_row), + DirectorCross(directors[first], unit_column)); + } + } + } + fesa::Vector velocity{24U}; + for (std::size_t node = 0U; node < 4U; ++node) { + velocity[6U * node] = static_cast(node + 1U); + velocity[6U * node + 1U] = -0.5 * static_cast(node); + velocity[6U * node + 3U] = 0.25; + } + for (std::size_t row = 0U; row < 24U; ++row) { + for (std::size_t column = 0U; column < 24U; ++column) { + EXPECT_NEAR(mass.Value().consistent_global24(row, column), + expected_consistent(row, column), 1e-12); + EXPECT_NEAR(mass.Value().lumped_global24(row, column), + expected_lumped(row, column), 1e-12); + } + } + fesa::Matrix consistent_difference{24U, 24U}; + fesa::Matrix lumped_difference{24U, 24U}; + for (std::size_t row = 0U; row < 24U; ++row) { + for (std::size_t column = 0U; column < 24U; ++column) { + consistent_difference(row, column) = + mass.Value().consistent_global24(row, column) - + expected_consistent(row, column); + lumped_difference(row, column) = mass.Value().lumped_global24(row, column) - + expected_lumped(row, column); + } + } + const double length = std::hypot(2.0, 0.3); + EXPECT_LE(FrobeniusNorm(consistent_difference, length) / + (FrobeniusNorm(mass.Value().consistent_global24, length) + + FrobeniusNorm(expected_consistent, length)), + 1e-12); + EXPECT_LE(FrobeniusNorm(lumped_difference, length) / + (FrobeniusNorm(mass.Value().lumped_global24, length) + + FrobeniusNorm(expected_lumped, length)), + 1e-12); + double pairwise_work = 0.0; + for (std::size_t first = 0U; first < 4U; ++first) { + for (std::size_t second = 0U; second < 4U; ++second) { + Coordinates displacement_difference{}; + Coordinates director_velocity_difference{}; + for (std::size_t component = 0U; component < 3U; ++component) { + displacement_difference[component] = + velocity[6U * first + component] - + velocity[6U * second + component]; + } + Coordinates first_rotation{{velocity[6U * first + 3U], + velocity[6U * first + 4U], + velocity[6U * first + 5U]}}; + Coordinates second_rotation{{velocity[6U * second + 3U], + velocity[6U * second + 4U], + velocity[6U * second + 5U]}}; + const auto first_director_velocity = + DirectorCross(directors[first], first_rotation); + const auto second_director_velocity = + DirectorCross(directors[second], second_rotation); + for (std::size_t component = 0U; component < 3U; ++component) + director_velocity_difference[component] = + first_director_velocity[component] - + second_director_velocity[component]; + pairwise_work += 0.5 * ( + mass.Value().consistent_moments[0U](first, second) * + Dot(displacement_difference, displacement_difference) + + 2.0 * mass.Value().consistent_moments[1U](first, second) * + Dot(displacement_difference, director_velocity_difference) + + mass.Value().consistent_moments[2U](first, second) * + Dot(director_velocity_difference, director_velocity_difference)); + } + } + const double consistent_work = + Quadratic(mass.Value().consistent_global24, velocity); + const double lumped_work = Quadratic(mass.Value().lumped_global24, velocity); + EXPECT_LE(std::abs((lumped_work - consistent_work) - pairwise_work) / + (lumped_work + consistent_work), + 1e-12); +} + // E2-C04: assemble the physical velocity work over successively finer meshes. TEST(Mitc4ShellMassContract, ReproducesAffineSquareVelocityWork) { const std::array directors{{{{0.0, 0.0, 1.0}}, diff --git a/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp index b48152c..926d667 100644 --- a/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp +++ b/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp @@ -1,5 +1,7 @@ +#define NOMINMAX #include "fesa/io/hdf5/eigen_hdf5_results_writer.h" +#include #include #include @@ -113,6 +115,22 @@ TEST(EigenHdf5ResultsWriter, WritesModalSchemaAndRejectsTamperedState) { const std::string retained{std::istreambuf_iterator{retained_file}, {}}; EXPECT_TRUE(retained == original); retained_file.close(); + // Denying delete sharing forces the actual atomic replacement to fail after + // writing and validating a fresh candidate, while retaining a valid final. + const HANDLE lock = CreateFileW(path.c_str(), GENERIC_READ, + FILE_SHARE_READ | FILE_SHARE_WRITE, nullptr, + OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, nullptr); + ASSERT_NE(lock, INVALID_HANDLE_VALUE); + const auto locked_write = writer.Write(path, domain, state, {}); + EXPECT_NE(CloseHandle(lock), 0); + ASSERT_FALSE(locked_write.IsOk()); + ASSERT_FALSE(locked_write.Diagnostics().empty()); + EXPECT_EQ(locked_write.Diagnostics()[0].code, "hdf5-finalization-failure"); + std::ifstream locked_final{path, std::ios::binary}; + const std::string after_failure{std::istreambuf_iterator{locked_final}, + {}}; + EXPECT_EQ(after_failure, original); + locked_final.close(); ASSERT_TRUE( fesa::hdf5_internal::SelfCheckEigenFile(path, domain, state, {}).IsOk()); const hid_t corrupt_file = diff --git a/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp b/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp index 5581709..b3ef36d 100644 --- a/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp +++ b/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp @@ -5,11 +5,15 @@ #include #include +#include +#include +#include #include #include #include "fesa/fem/dof_manager.h" #include "fesa/math/sparse_matrix.h" +#include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h" namespace { @@ -54,6 +58,59 @@ void RequireSuccess(const fesa::Result& result) { << diagnostics.front().message; } +void ExpectFailureCode(const fesa::Result& result, + const char* code) { + ASSERT_FALSE(result.HasValue()); + ASSERT_FALSE(result.GetStatus().Diagnostics().empty()); + EXPECT_EQ(result.GetStatus().Diagnostics().front().code, code); +} + +void ExpectSameBackendCalls(const std::vector& first, + const std::vector& second) { + ASSERT_EQ(first.size(), second.size()); + for (std::size_t index = 0U; index < first.size(); ++index) { + EXPECT_EQ(first[index].ordinal, second[index].ordinal); + EXPECT_EQ(first[index].kind, second[index].kind); + EXPECT_EQ(first[index].purpose, second[index].purpose); + EXPECT_DOUBLE_EQ(first[index].lower, second[index].lower); + EXPECT_DOUBLE_EQ(first[index].upper, second[index].upper); + EXPECT_EQ(first[index].allocated_m0, second[index].allocated_m0); + EXPECT_EQ(first[index].info, second[index].info); + EXPECT_EQ(first[index].count_available, second[index].count_available); + EXPECT_EQ(first[index].count, second[index].count); + EXPECT_EQ(first[index].outcome, second[index].outcome); + EXPECT_EQ(first[index].decision_reason, second[index].decision_reason); + } +} + +void ExpectSameProjector(const fesa::Matrix& first, const fesa::Matrix& second, + const std::size_t columns) { + ASSERT_EQ(first.Rows(), second.Rows()); + double squared_difference = 0.0; + for (std::size_t row = 0U; row < first.Rows(); ++row) { + for (std::size_t column = 0U; column < first.Rows(); ++column) { + double first_value = 0.0; + double second_value = 0.0; + for (std::size_t mode = 0U; mode < columns; ++mode) { + first_value += first(row, mode) * first(column, mode); + second_value += second(row, mode) * second(column, mode); + } + const double difference = first_value - second_value; + squared_difference += difference * difference; + } + } + EXPECT_LE(std::sqrt(squared_difference / + (2.0 * static_cast(columns))), + 1.0e-6); +} + +void ExpectThreadSettingRestored( + const fesa::Result& result, const char* code) { + ExpectFailureCode(result, code); + const int restored = mkl_set_num_threads_local(2); + EXPECT_EQ(restored, 2); +} + } // namespace TEST(MklFeastEigenSolver, SolvesN1AndRetainsBackendEvidence) { @@ -109,15 +166,28 @@ TEST(MklFeastEigenSolver, PreservesRepeatedMultiplicityAndStableOrder) { const auto first = solver.Solve(matrix, 3U, 1.0); const auto second = solver.Solve(matrix, 3U, 1.0); + const auto third = solver.Solve(matrix, 3U, 1.0); RequireSuccess(first); RequireSuccess(second); - if (!first.HasValue() || !second.HasValue()) return; + RequireSuccess(third); + if (!first.HasValue() || !second.HasValue() || !third.HasValue()) return; ASSERT_EQ(first.Value().eigenvalues.size(), 3U); EXPECT_NEAR(first.Value().eigenvalues[0U], 2.0, 1.0e-10); EXPECT_NEAR(first.Value().eigenvalues[1U], 2.0, 1.0e-10); EXPECT_NEAR(first.Value().eigenvalues[2U], 5.0, 1.0e-10); EXPECT_EQ(first.Value().column_ids, second.Value().column_ids); + EXPECT_EQ(first.Value().column_ids, third.Value().column_ids); + EXPECT_EQ(first.Value().eigenvalues, second.Value().eigenvalues); + EXPECT_EQ(first.Value().eigenvalues, third.Value().eigenvalues); + ExpectSameBackendCalls(first.Value().backend_calls, + second.Value().backend_calls); + ExpectSameBackendCalls(first.Value().backend_calls, + third.Value().backend_calls); + ExpectSameProjector(first.Value().eigenvectors, second.Value().eigenvectors, + 2U); + ExpectSameProjector(first.Value().eigenvectors, third.Value().eigenvectors, + 2U); ExpectUnitGram(first.Value().eigenvectors); } @@ -251,3 +321,311 @@ TEST(MklFeastEigenSolver, ResolvesCloseDistinctSpectrumAfterGivensCongruence) { 1e-9); ExpectUnitGram(result.Value().eigenvectors); } + +TEST(MklFeastEigenSolver, RejectsFatalCountStatusWithoutRetry) { + std::size_t hook_calls = 0U; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + ++hook_calls; + if (call.count && call.ordinal == 0U) outputs.info = -1; + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectFailureCode(result, "eigen-feast-failed"); + EXPECT_EQ(hook_calls, 1U); +} + +TEST(MklFeastEigenSolver, BoundsCountCallsAndUpperSearchReservations) { + constexpr std::size_t kMaximum = (std::numeric_limits::max)(); + EXPECT_TRUE(fesa::internal::CanCallFeast(511U)); + EXPECT_FALSE(fesa::internal::CanCallFeast(512U)); + EXPECT_FALSE(fesa::internal::CanCallFeast(kMaximum)); + EXPECT_TRUE(fesa::internal::CanStartUpperCandidate(363U)); + EXPECT_FALSE(fesa::internal::CanStartUpperCandidate(364U)); + EXPECT_FALSE(fesa::internal::CanStartUpperCandidate(kMaximum)); +} + +TEST(MklFeastEigenSolver, BoundsPrivateWorkspaceProductsBeforeAllocation) { + constexpr std::size_t kMaximum = (std::numeric_limits::max)(); + EXPECT_TRUE(fesa::internal::ProductFits(0U, kMaximum)); + EXPECT_TRUE(fesa::internal::ProductFits(1U, kMaximum)); + EXPECT_FALSE(fesa::internal::ProductFits(kMaximum, 2U)); + + const std::size_t byte_limit = static_cast( + (std::numeric_limits::max)() / sizeof(double)); + const std::size_t storage_limit = + (std::min)(byte_limit, std::vector{}.max_size()); + EXPECT_TRUE(fesa::internal::FitsDoubleStorage(0U)); + EXPECT_TRUE(fesa::internal::FitsDoubleStorage(storage_limit)); + EXPECT_FALSE(fesa::internal::FitsDoubleStorage(kMaximum)); + if (storage_limit < kMaximum) + EXPECT_FALSE(fesa::internal::FitsDoubleStorage(storage_limit + 1U)); +} + +TEST(MklFeastEigenSolver, TranslatesCountWorkspaceAllocationFailure) { + const int prior = mkl_set_num_threads_local(2); + fesa::internal::ScopedFeastWorkspaceTestHook hook([](const bool count) { + if (count) throw std::bad_alloc(); + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectThreadSettingRestored(result, "eigen-allocation-failed"); + static_cast(mkl_set_num_threads_local(prior)); +} + +TEST(MklFeastEigenSolver, TranslatesSolveWorkspaceAllocationFailure) { + const int prior = mkl_set_num_threads_local(2); + fesa::internal::ScopedFeastWorkspaceTestHook hook([](const bool count) { + if (!count) throw std::bad_alloc(); + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectThreadSettingRestored(result, "eigen-allocation-failed"); + static_cast(mkl_set_num_threads_local(prior)); +} + +TEST(MklFeastEigenSolver, StopsInitialInfoOneAfterThreeCountCalls) { + std::size_t count_calls = 0U; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (call.count) { + ++count_calls; + outputs.info = 1; + } + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectFailureCode(result, "eigen-feast-failed"); + EXPECT_EQ(count_calls, 3U); +} + +TEST(MklFeastEigenSolver, RejectsFinalCountHistoryMismatch) { + fesa::internal::ScopedFeastBackendTestHook hook( + [](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + // n=1 has total enclosure at ordinal 0 and final count at ordinal 1. + if (call.count && call.ordinal == 1U) outputs.m0 = 0; + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectFailureCode(result, "eigen-count-inconsistent"); +} + +TEST(MklFeastEigenSolver, RejectsDecreasingCountAcrossDistinctUpperGuards) { + constexpr std::size_t kSize = 11U; + std::vector values(kSize * kSize, 0.0); + for (std::size_t index = 0U; index < kSize; ++index) + values[index * kSize + index] = static_cast(index + 1U); + fesa::internal::ScopedFeastBackendTestHook hook( + [kSize](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + // Ordinals 1 and 2 are the first candidate's increasing upper guards. + if (call.count && call.ordinal == 1U) outputs.m0 = kSize; + if (call.count && call.ordinal == 2U) outputs.m0 = kSize - 1U; + }); + + const auto result = + fesa::MklFeastEigenSolver{}.Solve(MakeDenseCsr(kSize, values), 1U, 1.0); + + ExpectFailureCode(result, "eigen-count-inconsistent"); +} + +TEST(MklFeastEigenSolver, PreservesRawNegativeBoundaryAndRejectsBelowIt) { + const double threshold = -1.0e-9; + const std::vector accepted_negatives{ + -1.0e-12, threshold}; + for (const double negative : accepted_negatives) { + fesa::internal::ScopedFeastBackendTestHook hook( + [negative](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) { + ASSERT_EQ(outputs.eigenvalues.size(), 3U); + outputs.eigenvalues = {negative, 0.0, 1.0}; + } + }); + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(3U, {negative, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 1.0}), + 3U, 1.0); + RequireSuccess(result); + if (!result.HasValue()) continue; + ASSERT_EQ(result.Value().eigenvalues.size(), 3U); + EXPECT_DOUBLE_EQ(result.Value().eigenvalues[0U], negative); + EXPECT_DOUBLE_EQ(result.Value().eigenvalues[1U], 0.0); + EXPECT_DOUBLE_EQ(result.Value().eigenvalues[2U], 1.0); + } + + const double below_threshold = + std::nextafter(threshold, -(std::numeric_limits::infinity)()); + fesa::internal::ScopedFeastBackendTestHook hook( + [below_threshold](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) outputs.eigenvalues = {below_threshold, 0.0, 1.0}; + }); + const auto rejected = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(3U, {below_threshold, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 1.0}), + 3U, 1.0); + ExpectFailureCode(rejected, "negative-eigenvalue-spectrum"); +} + +TEST(MklFeastEigenSolver, DoesNotRepeatUnchangedSolveAfterReturnedMismatch) { + std::size_t solve_calls = 0U; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) { + ++solve_calls; + outputs.returned = 0; + } + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectFailureCode(result, "eigen-feast-failed"); + EXPECT_EQ(solve_calls, 1U); +} + +TEST(MklFeastEigenSolver, RetriesSolveInfoTwoOnlyWithLargerInterval) { + constexpr std::size_t kSize = 11U; + std::vector values(kSize * kSize, 0.0); + for (std::size_t index = 0U; index < kSize; ++index) + values[index * kSize + index] = static_cast(index + 1U); + std::vector solve_uppers; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) { + solve_uppers.push_back(call.upper); + if (solve_uppers.size() == 1U) outputs.info = 2; + } + }); + + const auto result = + fesa::MklFeastEigenSolver{}.Solve(MakeDenseCsr(kSize, values), 1U, 1.0); + + RequireSuccess(result); + if (!result.HasValue()) return; + ASSERT_EQ(solve_uppers.size(), 2U); + EXPECT_LT(solve_uppers[0U], solve_uppers[1U]); +} + +TEST(MklFeastEigenSolver, RejectsNonfiniteSolveOutputWithoutRetry) { + std::size_t solve_calls = 0U; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) { + ++solve_calls; + outputs.eigenvalues[0U] = + (std::numeric_limits::quiet_NaN)(); + } + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectFailureCode(result, "invalid-eigenvalue-result"); + EXPECT_EQ(solve_calls, 1U); +} + +TEST(MklFeastEigenSolver, RejectsResidualFailureWithoutSolveRetry) { + std::size_t solve_calls = 0U; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) { + ++solve_calls; + // N=1 still requires validation of the returned all-c=2 modes. + outputs.vectors[2U] = std::sqrt(0.5); + outputs.vectors[3U] = std::sqrt(0.5); + } + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(2U, {1.0, 0.0, 0.0, 2.0}), 1U, 1.0); + + ExpectFailureCode(result, "eigenpair-validation-failed"); + EXPECT_NE(result.GetStatus().Diagnostics().front().message.find( + "Independent whitened residual exceeds"), + std::string::npos); + EXPECT_EQ(solve_calls, 1U); +} + +TEST(MklFeastEigenSolver, RejectsNonfiniteVectorAndEpsoutWithoutRetry) { + for (const bool nonfinite_vector : {true, false}) { + std::size_t solve_calls = 0U; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) { + ++solve_calls; + if (nonfinite_vector) + outputs.vectors[0U] = + (std::numeric_limits::quiet_NaN)(); + else + outputs.epsout = (std::numeric_limits::quiet_NaN)(); + } + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(1U, {2.0}), 1U, 1.0); + + ExpectFailureCode(result, "invalid-eigenvalue-result"); + EXPECT_EQ(solve_calls, 1U); + } +} + +TEST(MklFeastEigenSolver, RejectsGramFailureWithoutSolveRetry) { + std::size_t solve_calls = 0U; + fesa::internal::ScopedFeastBackendTestHook hook( + [&](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + if (!call.count) { + ++solve_calls; + ASSERT_EQ(outputs.vectors.size(), 4U); + outputs.vectors = {1.0, 0.0, 1.0, 0.0}; + } + }); + + const auto result = fesa::MklFeastEigenSolver{}.Solve( + MakeDenseCsr(2U, {0.0, 0.0, 0.0, 0.0}), 2U, 1.0); + + ExpectFailureCode(result, "eigenpair-validation-failed"); + EXPECT_EQ(solve_calls, 1U); +} + +TEST(MklFeastEigenSolver, UsesCertifiedEnclosureAfterOptionalInfoOne) { + constexpr std::size_t kSize = 11U; + std::vector values(kSize * kSize, 0.0); + for (std::size_t index = 0U; index < kSize; ++index) + values[index * kSize + index] = static_cast(index + 1U); + fesa::internal::ScopedFeastBackendTestHook hook( + [](const fesa::internal::FeastBackendTestCall& call, + fesa::internal::FeastBackendTestOutputs& outputs) { + // The first upper guard follows the total enclosing count. + if (call.count && call.ordinal == 1U) outputs.info = 1; + }); + + const auto result = + fesa::MklFeastEigenSolver{}.Solve(MakeDenseCsr(kSize, values), 1U, 1.0); + + RequireSuccess(result); + if (!result.HasValue()) return; + EXPECT_EQ(result.Value().upper_search_termination, "optional-status"); + ASSERT_GE(result.Value().backend_calls.size(), 2U); + EXPECT_EQ(result.Value().backend_calls[1U].outcome, "optional-abandoned"); +}