diff --git a/FESA.sln b/FESA.sln index 64295f3..03a1a75 100644 --- a/FESA.sln +++ b/FESA.sln @@ -14,6 +14,8 @@ Project("{BC8A1FFA-BEE3-4634-8014-F334798102B3}") = "fesa_reference_tests", "msv EndProject Project("{BC8A1FFA-BEE3-4634-8014-F334798102B3}") = "gtest", "msvc\gtest.vcxproj", "{03194BA4-3F3B-5399-BCAF-FB4FA47300CE}" EndProject +Project("{BC8A1FFA-BEE3-4634-8014-F334798102B3}") = "fesa_jet_baseline_tests", "msvc\fesa_jet_baseline_tests.vcxproj", "{5FD398DA-6058-4CE3-9440-C9E0FDD180C3}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|x64 = Debug|x64 @@ -44,5 +46,9 @@ Global {03194BA4-3F3B-5399-BCAF-FB4FA47300CE}.Debug|x64.Build.0 = Debug|x64 {03194BA4-3F3B-5399-BCAF-FB4FA47300CE}.Release|x64.ActiveCfg = Release|x64 {03194BA4-3F3B-5399-BCAF-FB4FA47300CE}.Release|x64.Build.0 = Release|x64 + {5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Debug|x64.ActiveCfg = Debug|x64 + {5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Debug|x64.Build.0 = Debug|x64 + {5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Release|x64.ActiveCfg = Release|x64 + {5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Release|x64.Build.0 = Release|x64 EndGlobalSection EndGlobal diff --git a/docs/calculation-performance/README.md b/docs/calculation-performance/README.md index 746ba4e..3ed2955 100644 --- a/docs/calculation-performance/README.md +++ b/docs/calculation-performance/README.md @@ -10,10 +10,114 @@ Debug·Release 네 빌드는 모두 종료0이며 CLI와 세 테스트 실행 Debug는 `mkl_tbb_threadd`·`tbb12_debug`, Release는 `mkl_tbb_thread`·`tbb12`를 사용하고 HDF5의 `libmmd.dll`을 유지한다. 근거는 `.build/mkl-tbb-20261020-01/`에 보관했다. -사용자 지시로 이번 설정 변경에서는 테스트·모델·수치 비교·성능 측정·VTune을 실행하지 -않았다. 실제 DLL loading, 수치 정확성, 속도 향상과 OpenMP 대기 제거는 미검증이다. +사용자 지시로 당시 설정 변경에서는 테스트·모델·수치 비교·성능 측정·VTune을 실행하지 +않았다. 실제 DLL loading, 수치 정확성, 속도 향상과 OpenMP 대기 제거는 당시 미검증이었다. 기존 OpenMP 전용 local-thread/TLS helper와 HHT cache 진단 코드는 유지했으며 TBB에서의 -동등성도 확인하지 않았다. 아래 과거 시간·검증 기록을 새 TBB 구성의 결과로 승계하지 않는다. +동등성도 당시 확인하지 않았다. 이후 Jet 최적화의 새 TBB 검증·측정은 다음 절과 구분한다. +아래 과거 OpenMP 시간·검증 기록을 새 TBB 구성의 결과로 승계하지 않는다. + +## Jet Product 최적화와 현재 검증 + +2026-10-21 TBB 기준 `9ae7e6d414bfe36c20433c9dfa9fa4095aa9f7bd`와 최종 제품 source +`685f5e6a826c7bb66b9cef90e086ef4b2685252d`를 새로 비교했다. 후속 `3ede11b`는 +검증 근거만 기록한 빈 commit이며 제품 tree는 같다. 24자유도 Shell Jet과 private 누적 +helper만 수정했고 물리·공개 API·schema·수렴 기준은 유지했다. 기본 직접 작은 drilling +스프링과 기존 상대 비틀림을 모두 보존했으며 공식 시간 모델은 직접 방식을 사용했다. + +| 검토 후보 | 채택 결과와 이유 | +| --- | --- | +| 외적·Hessian 결합 | 모든 MXCSR 예외가 masked이고 첫 MKL 외적이 성공한 경우에만 그 전치를 재사용한다. 두 번째 외적 호출을 줄이되 각 성분의 덧셈 순서를 유지한다. Unmasked 또는 외적 거부 시 기존 순서로 계산한다. | +| 누적의 중간 저장 | 소비되는 임시 Jet이 할당된 Hessian을 단독 소유할 때만 재사용한다. 공유·동일 객체 alias·unmasked 예외는 기존 경로를 유지한다. | +| 상수 연산 특화 | 단순 scalar 대체는 원본 192개 경계 그룹 중 184개에서 signed zero 등의 결과가 달라 비채택했다. 기존 constant provenance와 implicit-zero는 유지한다. | +| 요청 없는 방향 채널 생략 | 비채택했다. 사용하지 않는 고차 채널도 overflow·예외·finite 실패를 결정하므로 단순 생략은 실패를 성공으로 바꿀 수 있다. 필요한 방향·혼합 도함수는 모두 유지한다. | +| 메모리 관리 | 서로 배타적인 빠른 경로와 기존 경로가 scratch 배열 하나를 공유한다. 단순 monotonic arena는 저장량 증가로 비채택했으며 범용 recycler는 구현·측정하지 않았다. | + +대표 local 평가의 Hessian 할당은 직접 방식에서 방향 없음 5416→4047회, 방향 있음 +6452→4911회로 줄었다. 기존 방식은 7456→5733회, 7636→5897회다. MSVC Release 곱셈의 +예약 stack은 통합 후보의 7696→5488 bytes로 2208 bytes 감소했지만 원래 5472 bytes보다 +16 bytes 크고 `__chkstk`도 유지한다. 단순 arena는 직접 평가의 payload를 331776에서 +18648576 bytes로 늘려 부적합했다. 이 할당·stack 수치를 전체 시간 단축률로 환산하지 않는다. +부동소수점 예외의 첫 trap, 수명·공유·동시·중첩 평가와 기존 private helper의 관성 consumer를 +독립 검토했고 FMA·재결합·전역 compiler 옵션은 추가하지 않았다. + +### 정식 모델 시간과 증분 기록 + +원본 `reference/shellNL-dyn-meshfine/shellNL-dyn-meshfine.inp`의 hash가 같은 새 복사본에서 +MSVC x64 Release·동일 MKL/TBB runtime DLL·기본 스레드 환경·전원 조건으로 각각 한 번 +실행했다. Ryzen7 5800H의 8 cores/16 logical processors 환경이며 별도 워밍업·VTune·동시 +빌드·시험·모델 실행은 없다. 타이머는 `Popen` 요청 직전부터 `wait` 반환 직후까지로, +읽기·계산·복구·로그·HDF5 close/finalize를 포함하고 준비·빌드·비교는 제외한다. + +| 실행 | 전체 초 | 분·초 | +| --- | ---: | ---: | +| 새 TBB 기준 | 4134.657782 | 68분 54.66초 | +| 최종 개선본 | 3747.788261 | 62분 27.79초 | + +이 단일 쌍에서 386.869521초(6분 26.87초), 9.356748% 감소·speedup 1.103226배를 확인해 +`improved`로 판정했다. 변경하지 않은 Native Hex8 선형 사례에서도 약4.9% 시간 감소가 +있었고, 긴 비선형·동일 CMake 경로의 Hex8/열 사례는 기준과 훨씬 가까웠다. 이는 실행 변동의 +문맥이며 공식 시간을 보정하거나 분산·통계 신뢰도를 추정하는 데 사용하지 않았다. 현재 +구체적인 환경·수치·반복 불일치가 없어 추가 쌍이나 profile은 실행하지 않았다. 한 쌍의 +실측이며 반복 통계·다른 기계에서의 보장·특정 함수만의 기여를 입증한 것은 아니다. +기존 FEM-REQ-PERF-001..009와 Shell 100초 목표는 변경하지 않았고 여전히 미달이다. + +두 실행 모두 증분0의 초기화와 양의 증분1..100을 구분했다. 모든 증분의 횟수가 같으며, +1..71은 Newton 보정 2회·로그 잔차 평가 3회, 72..100은 1회·2회다. 양의 증분 합은 +171회·271회이고 초기화의 잔차 기록 1개는 별도다. 이는 고유 Newton 로그 기록 수이며 +숨은 root·trial·assembly의 전체 호출 수가 아니다. 개선본의 초기화 9.044633초 + +스텝 합 3737.305479초 + 마무리 1.438149초는 전체 시간과 일치한다. + +스텝 시간은 10ms 간격의 외부 accepted 로그 관측으로 얻은 경과시간이며 순수 요소 계산 +시간이 아니다. 개선본의 최대 accepted 경계 관측 불확실성은 0.0329873초이고 누락·중복· +observer 오류·구간 합계 오차는 없다. 이 경계 불확실성은 스케줄링이나 실행 변동의 상한이 +아니다. [기준 CSV](../../.build/jet-product-20261021-01/baseline/time_steps.csv)와 +[개선본 CSV](../../.build/jet-product-20261021-01/candidate/time_steps.csv)에 101개 행의 +시간·보정·관측 평가 수와 경계 간격을 보존했다. + +### 정확성·회귀 결과와 남은 한계 + +정식 모델은 729절점·101프레임·schema2를 유지했다. HDF5의 2235개 dataset·27524개 +attribute를 읽고 12358개 수치 leaf의 15066400개 성분을 대조해 finite·signed zero· +정수 identity를 포함한 수치 bytes가 모두 같음을 확인했다. 문자열 차이는 열거한 입력 +경로 6개뿐이다. 전체729절점 내부 재현성과 외부9절점 참조의 U101 blocking family 통과를 +구분하며, UR100 failed warning family·경고200개는 그대로다. + +| 새 최종 검사 | 실제 실행 결과 | +| --- | --- | +| Native MSVC x64 Debug 전체 | 1035개 = 1018 PASS / 기존16 FAIL / 기존1 SKIP | +| 새 CMake Debug 전체 CTest | 1036개 = 1019 PASS / 기존16 FAIL / 기존1 SKIP | +| Python | 164 PASS | +| 세 구성의 독립·원본 Jet 대조 | 구성별 5종, 총15개 명령 PASS | + +Native 단위915·통합74·참조46개와 CTest의 실제 등록·실행 identity를 대조했다. CTest에는 +같은1035개 외 독립 Jet 검사1개가 포함된다. Native 단위 exit1·CTest exit8은 기준에서 +재현된 실패 때문이며 전체 green이 아니다. 모든 기존 실패 assertion·수치·skip 진단은 +명시한 source 경로와 GTest 시간 줄 외에 동일했고 신규 회귀는 없다. 원본 Hex8 동해석 +101프레임 두 사례, 열 2차원 격자 수렴과 슬래브12개를 포함해 실행했다. 각 구성의 독립 +항등식·원본42 private/64 public 사례·곱셈3072·누적192·alias8 대조도 통과했으며 기존 +snapshot은 재생성하지 않았다. 결과 비교 도구는 동일 복사본/허용 경로 차이 2개와 값· +signed zero·ID·예상 밖 문자열·NaN 검출 5개 대조 사례를 통과했다. + +참조46개는 물리 모델46개의 정확성 통과를 뜻하지 않는다. 실제 분류는 positive Shell +검사3개·기타2개·policy/error 계약40개·기존 열 NT11 불일치 진단1개다. NT11의 87 failed +rows·정규화 RMS0.049674084615673193은 기준과 동일한 수치 실패이며 진단 검사의 PASS를 +열 정확성 PASS로 바꾸지 않았다. 정식 `reference/shellNL-meshfine`와 +`reference/shellNL-meshfine_02`가 없고 후자의 전용 비교 계약도 없다. `reference/temp`는 +사용하지 않았으며 정식 원본36개 파일의 전후 hash는 같다. 기존16 FAIL·1 SKIP·자료 누락, +HHT 물리·시간 정확도와 열 필수 참조의 미해결 범위가 남아 전체 물리 검증은 미완료다. + +기존 `.build/cmake-debug` 캐시는 다른 source root를 가리켜 그대로 보존했다. 이번에는 +동일 preset에 `-B .build/cmake-jet-debug`를 지정하고 그 경로를 build·CTest `--test-dir`에 +사용했다. 실제231개 compile command와 전체 검사 등록이 명시한 source/build 경로 외에는 +기준과 같음을 확인했다. Native Release 명령의 include 경로 따옴표 차이는 Windows argv로 +대조했고 옵션·순서·값은 같다. 과거 compiler binary SHA는 확보되지 않았으므로 동일 hash를 +주장하지 않으며 보존된 compiler/pass 버전·경로·실제 옵션·SDK/runtime와 현재 hash를 기록했다. + +전체 근거는 [최종 비교](../../.build/jet-product-20261021-01/comparison.json), +[검사·실행 요약](../../.build/jet-product-20261021-01/t6/T6-summary.json), +[독립 검토](../../.build/jet-product-20261021-01/review/t6-final.json)에 있다. 각 실행의 source· +실행물·입력·결과 hash와 raw·즉시 exit를 보존했으며, 자료는 로컬 `.build`에 있고 Git에 +포함하지 않는다. 문서 갱신 뒤 제품 source가 바뀌지 않았으므로 build/test를 반복하지 않았다. ## 채택한 계산 경로와 수명 diff --git a/docs/hht-implicit-dynamics/README.md b/docs/hht-implicit-dynamics/README.md index d6209b4..edae672 100644 --- a/docs/hht-implicit-dynamics/README.md +++ b/docs/hht-implicit-dynamics/README.md @@ -9,13 +9,22 @@ polar/swing/twist를 생략하고 physical 회전·관성·고정 null 기저· `sum_e kd*gamma_ei*(n_i·d_i0)-n_i·Fext_rot`다. 그 같은 식의 gradient·Hessian· 필요한 3차 방향미분을 사용한다. 초기 director를 현재 null로 고정하지 않는다. -현재 source `6bbbb2b`를 전체 실행한 Debug Native/CTest는 각각1035개 중 +2026-10-21 현재 Jet 최적화 제품 source는 `685f5e6`이다. 새 TBB 기준 `9ae7e6d`와의 +전체 회귀에서 Native1035개는 1018 PASS·기존16 FAIL·1 SKIP, CTest1036개는 +1019 PASS·같은16 FAIL·1 SKIP, Python164개는 PASS이며 신규 회귀는 없다. +직접 작은 스프링의 정식729절점·101프레임 수치 bytes와 모든 증분의 Newton 기록 수는 +기준과 같고 외부9절점 U101 families는 PASS, UR 경고200개는 유지했다. 동일 Release 조건의 +공식 한 쌍은 4134.657782→3747.788261초로 9.36% 짧게 측정됐다. 반복 통계나 전체 물리· +시간 정확도 통과를 뜻하지 않으며 기존16 FAIL·1 SKIP·정식 참조 누락은 남는다. +후보별 채택·비채택과 측정 범위는 [계산 성능 결과](../calculation-performance/README.md)에 있다. + +2026-10-20 source `6bbbb2b`를 당시 전체 실행한 Debug Native/CTest는 각각1035개 중 1024 PASS·기존10 FAIL·1 INCONCLUSIVE SKIP이며 별도 Python164개는 PASS다. 정식729node·101frame Direct Release 해석과 외부9node U의101 families 비교는 CLI0/비교0, 최대 정규화 RMS0.2282326248391787%로 통과했고 UR 경고200개를 유지했다. 기존 방식은 별도 Debug의 관련 검사와 실제 소형3frame HHT 출력을 확인했다. 729node 기존 방식의 외부 U 재해석은 이번에 하지 않았다. 과거 결과를 현재 PASS로 승계하지 않는다. -기존10개 물리 실패·N7 미판정·정식 자료 누락이 남으므로 전체 검증 완료는 아니다. +당시10개 물리 실패·N7 미판정·정식 자료 누락으로 전체 검증 완료는 아니었다. 아래 2026-10-19 수용 기록은 이전 수렴 실험의 결정이며 현재 전체 PASS를 뜻하지 않는다. 2026-10-19 후속 사용자 승인으로 이번 수렴 규칙 실험은 **Abaqus 변위 검증 통과·실패10개 보류로 조건부 완료**한다. diff --git a/msvc/fesa_jet_baseline_tests.vcxproj b/msvc/fesa_jet_baseline_tests.vcxproj new file mode 100644 index 0000000..0ab92c0 --- /dev/null +++ b/msvc/fesa_jet_baseline_tests.vcxproj @@ -0,0 +1,21 @@ + + + Debugx64Releasex64 + {5FD398DA-6058-4CE3-9440-C9E0FDD180C3}Win32Proj10.0 + + Applicationv145true + + + + + $(FesaRoot)src\fesa;$(FESA_MKL_ROOT)\include;$(FESA_TBB_ROOT)\include;$(FESA_HDF5_ROOT)\include;%(AdditionalIncludeDirectories) + MKL_ILP64;__TBB_NO_IMPLICIT_LINKAGE=1;H5_BUILT_AS_DYNAMIC_LIB;%(PreprocessorDefinitions) + + + + + {3AEC39D9-4641-5161-9107-265EB50CA882}true + + + + diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp index e8d41e9..158d59e 100644 --- a/src/fesa/elements/nonlinear_mitc4_shell.cpp +++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp @@ -1,5 +1,7 @@ #include "fesa/elements/nonlinear_mitc4_shell.h" +#include + #include #include #include @@ -384,6 +386,9 @@ class Hessian { /// @brief Identifies only the immutable zero Hessian, independently of seeds. bool IsImplicitZero() const { return !values_; } + /// @brief Allows a consumed temporary to reuse only its own materialized data. + bool HasUniqueStorage() const { return values_ && values_.use_count() == 1; } + /// @brief Gets fresh output storage without redundant zero initialization. /// @note Only fresh operator results use this pointer; all 576 entries must /// be written before any read or copy of the result. @@ -412,11 +417,14 @@ class Hessian { /// @brief Forms an ordinary derivative outer product without a Matrix copy. /// @note Rejected nonfinite operands propagate to the existing Jet finite gate. -void Outer(const std::array& left, +bool Outer(const std::array& left, const std::array& right, double* result) { if (!dense_blas_internal::OuterProduct(left.data(), left.size(), right.data(), - right.size(), result)) + right.size(), result)) { std::fill_n(result, kHessianSize, std::numeric_limits::quiet_NaN()); + return false; + } + return true; } /// @brief Carries exact first and second trial-coordinate variations. @@ -483,6 +491,30 @@ Jet operator+(const Jet& left, const Jet& right) { return result; } +/// @brief Adds to a consumed temporary without allocating another Hessian. +/// @note Shared/live aliases and unmasked exceptions retain the original path. +Jet operator+(Jet&& left, const Jet& right) { + if ((_mm_getcsr() & _MM_MASK_MASK) != _MM_MASK_MASK || &left == &right || + !left.second.HasUniqueStorage()) + return static_cast(left) + right; + left.value = left.value + right.value; + left.directional_first = left.directional_first + right.directional_first; + left.directional_second = left.directional_second + right.directional_second; + for (std::size_t index = 0U; index < kGlobalDofCount; ++index) { + left.first[index] = left.first[index] + right.first[index]; + left.mixed_first[index] = left.mixed_first[index] + right.mixed_first[index]; + left.mixed_second[index] = left.mixed_second[index] + right.mixed_second[index]; + } + left.constant_derivatives = left.constant_derivatives && + right.constant_derivatives && + std::isfinite(left.value); + const double* right_second = right.second.ReadData(); + double* left_second = left.second.Data(); + for (std::size_t index = 0U; index < kHessianSize; ++index) + left_second[index] = left_second[index] + right_second[index]; + return std::move(left); +} + /// @brief Subtracts two second-order trial-coordinate quantities. Jet operator-(const Jet& left, const Jet& right) { Jet result{}; @@ -559,14 +591,45 @@ Jet operator*(const Jet& left, const Jet& right) { left.value * right_second[index]; return result; } - Outer(left.first, right.first, result_second); - std::array reverse; - Outer(right.first, left.first, reverse.data()); + const bool masked_exceptions = + (_mm_getcsr() & _MM_MASK_MASK) == _MM_MASK_MASK; + const bool outer_accepted = Outer(left.first, right.first, result_second); + std::array scratch; + if (masked_exceptions && outer_accepted) { + // Retain earlier upper entries for later lower entries before overwriting + // them. Reuse the MKL k=1 transpose while preserving the original row-major + // arithmetic order as well as each entry's four-term association. + std::size_t saved = 0U; + for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { + for (std::size_t column = 0U; column < row; ++column) { + const std::size_t index = row * kGlobalDofCount + column; + result_second[index] = left_second[index] * right.value + + result_second[index] + scratch[saved++] + + left.value * right_second[index]; + } + const std::size_t diagonal = row * kGlobalDofCount + row; + const double product = result_second[diagonal]; + result_second[diagonal] = left_second[diagonal] * right.value + + product + product + left.value * right_second[diagonal]; + for (std::size_t column = row + 1U; column < kGlobalDofCount; ++column) { + const std::size_t index = row * kGlobalDofCount + column; + const std::size_t transpose = column * kGlobalDofCount + row; + const double forward = result_second[index]; + const double reverse = result_second[transpose]; + scratch[column * (column - 1U) / 2U + row] = forward; + result_second[index] = left_second[index] * right.value + + forward + reverse + left.value * right_second[index]; + } + } + return result; + } + // Preserve the reverse input-validation order after a rejected first outer. + Outer(right.first, left.first, scratch.data()); for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { const std::size_t index = row * kGlobalDofCount + column; result_second[index] = left_second[index] * right.value + - result_second[index] + reverse[index] + + result_second[index] + scratch[index] + left.value * right_second[index]; } } @@ -1068,7 +1131,7 @@ Jet FrameComponent(const JetVector3& left, const JetTensor& tensor, Jet result{}; for (std::size_t row = 0U; row < 3U; ++row) { for (std::size_t column = 0U; column < 3U; ++column) { - result = result + left[row] * tensor[row][column] * right[column]; + result = std::move(result) + left[row] * tensor[row][column] * right[column]; } } return result; @@ -1146,9 +1209,9 @@ Jet PhysicalEnergy(const std::array& coordinates, shear * strain[4U]}; Jet density{}; for (std::size_t component = 0U; component < strain.size(); ++component) { - density = density + strain[component] * stress[component]; + density = std::move(density) + strain[component] * stress[component]; } - energy = energy + 0.5 * reference.reference_jacobian * density; + energy = std::move(energy) + 0.5 * reference.reference_jacobian * density; } } if (!IsFinite(energy)) { diff --git a/src/fesa/math/jet_arithmetic.h b/src/fesa/math/jet_arithmetic.h index 439e439..986e74f 100644 --- a/src/fesa/math/jet_arithmetic.h +++ b/src/fesa/math/jet_arithmetic.h @@ -3,6 +3,7 @@ #include #include +#include namespace fesa::math_internal { /// @brief Multiplies fixed Jet matrices in row/column/index order. @@ -17,7 +18,7 @@ std::array, 3> MultiplyJetMatrix( for (std::size_t column = 0; column < 3; ++column) { result[row][column] = constant(0.0); for (std::size_t index = 0; index < 3; ++index) - result[row][column] = result[row][column] + + result[row][column] = std::move(result[row][column]) + left[row][index] * right[index][column]; } return result; @@ -31,7 +32,7 @@ std::array MultiplyJetMatrixVector( std::array result{constant(0.0),constant(0.0),constant(0.0)}; for (std::size_t row = 0; row < 3; ++row) for (std::size_t column = 0; column < 3; ++column) - result[row] = result[row] + matrix[row][column] * vector[column]; + result[row] = std::move(result[row]) + matrix[row][column] * vector[column]; return result; } } // namespace fesa::math_internal diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index af52699..94358c3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -179,6 +179,17 @@ target_compile_definitions( FESA_TEST_BINARY_DIR="${PROJECT_BINARY_DIR}" ) +add_executable(fesa_jet_baseline_tests + unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp) +target_include_directories(fesa_jet_baseline_tests PRIVATE + "${PROJECT_SOURCE_DIR}/src/fesa") +target_compile_definitions(fesa_jet_baseline_tests PRIVATE + MKL_ILP64 __TBB_NO_IMPLICIT_LINKAGE=1 H5_BUILT_AS_DYNAMIC_LIB) +target_compile_options(fesa_jet_baseline_tests PRIVATE + $<$:/fp:precise /EHsc>) +target_link_libraries(fesa_jet_baseline_tests PRIVATE + fesa_solver Fesa::MKL Fesa::TBB Fesa::HDF5) + if(WIN32) if(NOT TARGET ${FESA_MKL_THREAD_TARGET} OR NOT TARGET MKL::mkl_core) message(FATAL_ERROR "The oneMKL runtime files required by tests were not resolved") @@ -226,6 +237,7 @@ if(WIN32) fesa_stage_test_runtime(fesa_unit_tests) fesa_stage_test_runtime(fesa_integration_tests) fesa_stage_test_runtime(fesa_reference_tests) + fesa_stage_test_runtime(fesa_jet_baseline_tests) unset(_fesa_mkl_dispatch_dll) unset(_fesa_mkl_dispatch_dlls) @@ -258,4 +270,11 @@ add_custom_target( fesa_unit_tests fesa_integration_tests fesa_reference_tests + fesa_jet_baseline_tests ) + +# Verification never creates or replaces an immutable baseline snapshot. +add_test(NAME NonlinearMitc4Jet.IndependentOracle + COMMAND fesa_jet_baseline_tests --verify) +set_tests_properties(NonlinearMitc4Jet.IndependentOracle PROPERTIES + LABELS "unit;jet-baseline") diff --git a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp new file mode 100644 index 0000000..a05264d --- /dev/null +++ b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp @@ -0,0 +1,935 @@ +// Standalone private baseline fixture: compile this translation unit instead of +// nonlinear_mitc4_shell.cpp, then link the ordinary solver archive. No +// production source, public interface, or official CLI instrumentation is +// needed. Build target: fesa_jet_baseline_tests (MSBuild/CMake). --verify runs +// independent checks without writing evidence. Compare immutable evidence with +// --compare PATH, Debug to Debug or Release to Release; never regenerate a +// snapshot during regression. See nonlinear_mitc4_shell_jet_baseline.md. +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../../../src/fesa/elements/nonlinear_mitc4_shell.cpp" + +namespace fesa { +namespace { +using Coefficients = std::array; +std::fenv_t initial_environment; +unsigned initial_mxcsr; + +std::string FpState() { + unsigned control = 0; + _controlfp_s(&control, 0, 0); + std::ostringstream out; + out << std::fetestexcept(FE_ALL_EXCEPT) << ':' << std::fegetround() << ':' + << _statusfp() << ':' << control << ':' << _mm_getcsr(); + return out.str(); +} + +std::string ResetFp() { + std::fesetenv(&initial_environment); + std::feclearexcept(FE_ALL_EXCEPT); + _clearfp(); + _mm_setcsr(initial_mxcsr & ~0x3fu); + return FpState(); +} + +struct Monomial { + double coefficient; + std::array powers{}; +}; + +// Independent oracle: differentiate monomial exponent lists, never Jet rules. +double Differentiate(const std::vector& polynomial, + const Coefficients& point, + const std::vector& derivative) { + long double sum = 0; + for (auto term : polynomial) { + long double value = term.coefficient; + for (auto index : derivative) { + value *= term.powers[index]; + if (term.powers[index] == 0) { + value = 0; + break; + } + --term.powers[index]; + } + for (std::size_t i = 0; i < 24; ++i) + for (unsigned p = 0; p < term.powers[i]; ++p) value *= point[i]; + sum += value; + } + return static_cast(sum); +} + +std::uint64_t Bits(double value) { + std::uint64_t result; + std::memcpy(&result, &value, sizeof(result)); + return result; +} +void Check(double actual, double expected, const char* channel) { + if (!std::isfinite(actual) || + std::abs(actual - expected) > 2e-12 * std::max(1.0, std::abs(expected))) + throw std::runtime_error(std::string("independent oracle: ") + channel); +} +Jet Polynomial(const std::vector& terms, + const std::array& variables) { + Jet result = Constant(0); + for (const auto& term : terms) { + Jet value = Constant(term.coefficient); + for (std::size_t i = 0; i < 24; ++i) + for (unsigned p = 0; p < term.powers[i]; ++p) + value = value * variables[i]; + result = result + value; + } + return result; +} +std::string Snapshot(const char* name, const Jet& jet, int exceptions, + unsigned status, unsigned mxcsr) { + std::ostringstream out; + out << name << ' ' << jet.constant_derivatives << ' ' + << jet.second.IsImplicitZero() << ' ' << IsFinite(jet) << ' ' + << exceptions << ' ' << status << ' ' << mxcsr << std::hex; + auto put = [&](double x) { + if (std::isnan(x)) + out << " nan"; + else + out << ' ' << Bits(x); + }; + put(jet.value); + for (auto x : jet.first) put(x); + for (auto x : jet.second) put(x); + put(jet.directional_first); + put(jet.directional_second); + for (auto x : jet.mixed_first) put(x); + for (auto x : jet.mixed_second) put(x); + return out.str() + "\n"; +} +template +std::string Observe(const char* name, F expression, Jet* returned = nullptr) { + const auto before = ResetFp(); + Jet value = expression(); + const int exceptions = std::fetestexcept(FE_ALL_EXCEPT); + const unsigned status = _statusfp(); + const unsigned mxcsr = _mm_getcsr(); + const auto after = FpState(); + if (returned != nullptr) *returned = value; + return std::string("initial-fp=") + before + " final-fp=" + after + " " + + Snapshot(name, value, exceptions, status, mxcsr); +} + +std::string Verify(bool corrupt_oracle) { + Coefficients x{}, direction{}; + std::array variables{}; + std::vector terms; + for (std::size_t i = 0; i < 24; ++i) { + x[i] = (static_cast(i % 7) - 3) / 8; + direction[i] = (static_cast(i % 5) - 2) / 16; + variables[i] = Variable(x[i], i, direction[i]); + Monomial term{(i % 2 ? -1.0 : 1.0) / 8}; + term.powers[i] = 3; + term.powers[(i + 1) % 24] += 2; + terms.push_back(term); + Monomial cross{0.5}; + cross.powers[i] = 1; + cross.powers[(i + 11) % 24] += 1; + terms.push_back(cross); + } + Jet polynomial = Polynomial(terms, variables); + if (corrupt_oracle) polynomial.value += 1; // Deliberate fixture-only mutant. + Check(polynomial.value, Differentiate(terms, x, {}), "value"); + double d1 = 0, d2 = 0; + for (std::size_t i = 0; i < 24; ++i) { + const double gradient = Differentiate(terms, x, {i}); + Check(polynomial.first[i], gradient, "gradient"); + d1 += gradient * direction[i]; + double mixed1 = 0, mixed2 = 0; + for (std::size_t j = 0; j < 24; ++j) { + const double hessian = Differentiate(terms, x, {i, j}); + Check(polynomial.second.ReadData()[24 * i + j], hessian, "Hessian"); + d2 += hessian * direction[i] * direction[j]; + mixed1 += hessian * direction[j]; + for (std::size_t k = 0; k < 24; ++k) + mixed2 += + Differentiate(terms, x, {i, j, k}) * direction[j] * direction[k]; + } + Check(polynomial.mixed_first[i], mixed1, "mixed_first"); + Check(polynomial.mixed_second[i], mixed2, "mixed_second"); + } + Check(polynomial.directional_first, d1, "directional_first"); + Check(polynomial.directional_second, d2, "directional_second"); + // Affine composition has a rank-one Hessian; its third derivative contracts + // independently to the mixed-second channel. + Jet affine = Constant(2); + Coefficients weights{}; + double slope = 0; + for (std::size_t i = 0; i < 24; ++i) { + weights[i] = (static_cast(i % 3) - 1) / 64; + affine = affine + weights[i] * variables[i]; + slope += weights[i] * direction[i]; + } + auto unary = [&](const Jet& jet, double value, double a, double b, double c) { + Check(jet.value, value, "unary value"); + Check(jet.directional_first, a * slope, "unary direction1"); + Check(jet.directional_second, b * slope * slope, "unary direction2"); + for (std::size_t i = 0; i < 24; ++i) { + Check(jet.first[i], a * weights[i], "unary gradient"); + Check(jet.mixed_first[i], b * weights[i] * slope, "unary mixed1"); + Check(jet.mixed_second[i], c * weights[i] * slope * slope, + "unary mixed2"); + for (std::size_t j = 0; j < 24; ++j) + Check(jet.second.ReadData()[24 * i + j], b * weights[i] * weights[j], + "unary Hessian"); + } + }; + const double t = affine.value, r = std::sqrt(t); + unary(Sin(affine), std::sin(t), std::cos(t), -std::sin(t), -std::cos(t)); + unary(Cos(affine), std::cos(t), -std::sin(t), -std::cos(t), std::sin(t)); + unary(Reciprocal(affine), 1 / t, -1 / (t * t), 2 / (t * t * t), + -6 / (t * t * t * t)); + unary(Sqrt(affine), r, 0.5 / r, -0.25 / (r * r * r), + 0.375 / (r * r * r * r * r)); + const double den = 1 + t * t; + unary(Atan2(affine, Constant(1)), std::atan(t), 1 / den, -2 * t / (den * den), + (6 * t * t - 2) / (den * den * den)); + const Jet original = polynomial; + const Jet alias = original; + const auto before = Snapshot("alias", alias, 0, 0, 0); + const Jet transformed = alias * alias + original; + if (before != Snapshot("alias", original, 0, 0, 0)) + throw std::runtime_error("alias modified"); + std::string out; + std::vector>> futures; + for (int i = 0; i < 4; ++i) + futures.push_back(std::async(std::launch::async, [&] { + Jet value; + const auto evidence = Observe( + "shared-full", [&] { return original * original + original; }, + &value); + return std::make_pair(std::move(value), evidence); + })); + for (std::size_t i = 0; i < futures.size(); ++i) { + const auto worker = futures[i].get(); + if (Snapshot("shared", worker.first, 0, 0, 0) != + Snapshot("shared", transformed, 0, 0, 0)) + throw std::runtime_error("shared complete-channel snapshot"); + out += "worker-" + std::to_string(i) + " " + worker.second; + } + auto add = [&](const char* name, auto f) { out += Observe(name, f); }; + add("constant", [] { return Constant(-0.0); }); + add("variable-zero-direction", [] { return Variable(-0.0, 23); }); + add("variable-direction", [] { return Variable(0.25, 0, -0.5); }); + add("polynomial", [&] { return Polynomial(terms, variables); }); + add("polynomial-zero-direction", [&] { + auto zero_direction = variables; + for (auto& variable : zero_direction) variable.directional_first = 0; + return Polynomial(terms, zero_direction); + }); + add("actual-zero-hessian", [&] { return original - original; }); + add("implicit-zero-nonconstant", + [] { return Variable(1, 0) - Variable(1, 0); }); + add("implicit-add", [] { return Variable(1, 0) + Constant(-1); }); + add("implicit-subtract", [] { return Variable(1, 0) - Constant(1); }); + for (int kind = 0; kind < 4; ++kind) { + const auto name = "surviving-copy-" + std::to_string(kind); + add(name.c_str(), [&] { + Jet survivor; + std::string expected; + { + Jet source = kind == 0 ? Constant(-0.0) + : kind == 1 ? Variable(0.5, 7, -0.25) + : Polynomial(terms, variables); + if (kind == 3) source = source - source; + const Jet another_alias = source; + expected = Snapshot("survivor", source, 0, 0, 0); + survivor = another_alias; + } + if (Snapshot("survivor", survivor, 0, 0, 0) != expected) + throw std::runtime_error("surviving copy complete-channel snapshot"); + return survivor; + }); + } + add("constant-multiply", [] { return Constant(-0.0) * Constant(3); }); + add("constant-compose", [] { return Sin(Constant(0)); }); + add("scale", [&] { return -0.5 * original; }); + add("multiply", [&] { return original * original; }); + add("compose-sin", [&] { return Sin(original); }); + add("compose-cos", [&] { return Cos(original); }); + add("divide", [&] { return original / affine; }); + add("sqrt", [&] { return Sqrt(affine); }); + add("atan2", [&] { return Atan2(original, affine); }); + add("dot", [&] { + return Dot(JetVector3{original, affine, variables[0]}, + JetVector3{affine, variables[1], original}); + }); + add("cancellation", + [] { return (Variable(1e16, 0, 1) + Constant(1)) - Constant(1e16); }); + add("overflow", [] { + return Variable(std::numeric_limits::max(), 0, 1) * Constant(2); + }); + add("subnormal", [] { + return std::numeric_limits::denorm_min() * Variable(0.5, 0, 0.5); + }); + add("nan", [] { + return Constant(std::numeric_limits::quiet_NaN()) * Variable(1, 0); + }); + add("inf", [] { + return Constant(std::numeric_limits::infinity()) * Constant(0); + }); + add("singular", [] { return Reciprocal(Variable(0, 0, 1)); }); + add("singular-zero-direction", [] { return Reciprocal(Variable(0, 0)); }); + add("constant-singular-reciprocal", [] { return Reciprocal(Constant(0)); }); + add("constant-singular-sqrt", [] { return Sqrt(Constant(0)); }); + add("negative-sqrt", [] { return Sqrt(Variable(-1, 0, 1)); }); + add("after-failure", [&] { return Polynomial(terms, variables); }); + // Exercise the production MKL k=1 outer path directly, including entries + // whose rounded bits agree while their IEEE exception state can differ. + auto outer = [](double a, double b) { + Jet jet{}; + Coefficients left{}, right{}; + for (std::size_t i = 0; i < 24; ++i) { + left[i] = i % 2 ? -a : a; + right[i] = i % 3 ? -b : b; + } + Outer(left, right, jet.second.OverwriteData()); + return jet; + }; + add("outer-signed-zero", [&] { return outer(-0.0, 1); }); + add("outer-underflow", + [&] { return outer(std::numeric_limits::denorm_min(), 0.5); }); + add("outer-overflow", + [&] { return outer(std::numeric_limits::max(), 2); }); + add("outer-invalid", + [&] { return outer(std::numeric_limits::infinity(), 0); }); + if (!Constant(1).constant_derivatives || + !Constant(1).second.IsImplicitZero() || + Variable(1, 0).constant_derivatives || + !Variable(1, 0).second.IsImplicitZero() || + polynomial.second.IsImplicitZero()) + throw std::runtime_error("provenance"); + return out; +} + +// Separate supplemental evidence: never changes the original T1 snapshot. +// A changed Hessian association, a lost transpose term, or a changed MKL zero +// convention must fail the bit/classification and complete FP-state comparison. +std::string ProductSnapshots(bool omit_transpose) { + const std::array boundary{ + 0.0, -0.0, 1.0, -1.0, 0.1, -0.3, + std::numeric_limits::denorm_min(), + -std::numeric_limits::denorm_min(), + std::numeric_limits::min(), + std::nextafter(std::numeric_limits::min(), 0.0), + std::numeric_limits::max(), + -std::numeric_limits::max(), + std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::signaling_NaN()}; + std::string out; + for (const int rounding : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, + FE_TOWARDZERO}) { + for (const unsigned flush : {0u, 0x40u, 0x8000u, 0x8040u}) { + for (const bool sticky : {false, true}) { + auto reset = [&] { + ResetFp(); + std::fesetround(rounding); + _mm_setcsr((_mm_getcsr() & ~0x8040u) | flush); + if (sticky) std::feraiseexcept(FE_ALL_EXCEPT); + return FpState(); + }; + for (std::size_t kind = 0; kind < 96; ++kind) { + ResetFp(); // Construct inputs outside the observed environment. + Jet left{}, right{}; + left.value = 0.75; + right.value = -1.25; + left.directional_first = 0.25; + right.directional_first = -0.5; + left.directional_second = -0.125; + right.directional_second = 0.375; + std::uint64_t random = 0x719b245af51d6723ULL + kind; + auto coefficient = [&] { + random ^= random << 13; + random ^= random >> 7; + random ^= random << 17; + return (static_cast(random % 2000001) - 1000000) / + 1048576.0; + }; + for (std::size_t i = 0; i < 24; ++i) { + left.first[i] = kind < 32 ? boundary[(kind + i) % 16] + : coefficient(); + right.first[i] = kind < 32 ? boundary[(kind / 2 + 3 * i) % 16] + : coefficient(); + left.mixed_first[i] = coefficient(); + right.mixed_first[i] = coefficient(); + left.mixed_second[i] = coefficient(); + right.mixed_second[i] = coefficient(); + } + // Include dense, asymmetric inputs, implicit zeros and aliased input. + if (kind % 3 != 0) { + double* l = left.second.OverwriteData(); + double* r = right.second.OverwriteData(); + for (std::size_t i = 0; i < 576; ++i) { + l[i] = coefficient(); + r[i] = coefficient(); + } + } + if (kind == 32) left = Constant(-0.0); + if (kind == 33) right = Constant(0.0); + if (kind == 34) left.value = boundary[10]; + if (kind == 35) right.value = boundary[12]; + if (kind == 36) left.second.Data()[17] = boundary[10]; + if (kind == 37) right.second.Data()[24] = boundary[10]; + if (kind == 38) left.second.Data()[17] = boundary[14]; + if (kind == 39) right = left; // Shared immutable Hessian storage. + if (kind == 40) { + left.first[0] = boundary[12]; + right.first[0] = boundary[15]; + } + if (kind == 41) { + left.first[0] = boundary[15]; + right.first[0] = boundary[12]; + } + if (kind >= 42 && kind < 64) { + for (std::size_t i = 0; i < 24; ++i) { + left.first[i] = boundary[(kind + i) % 12]; + right.first[i] = boundary[(kind / 2 + 3 * i) % 12]; + if (kind == 42) { + left.first[i] = boundary[6]; + right.first[i] = i % 2 ? -0.5 : 0.5; + } + if (kind == 43) { + left.first[i] = boundary[10]; + right.first[i] = i % 2 ? -2.0 : 2.0; + } + if (kind == 44) { + left.first[i] = i % 2 ? -0.0 : 0.0; + right.first[i] = i % 3 ? -1.0 : 1.0; + } + } + } + const auto left_before = Snapshot("input", left, 0, 0, 0); + const auto right_before = Snapshot("input", right, 0, 0, 0); + const auto before = reset(); + Jet product = left * right; + if (omit_transpose && kind >= 64) { + // Deliberately omit the reverse outer term in one finite entry. + // This mutant tests arithmetic detection, not just file integrity. + product.second.Data()[1] = + left.second.ReadData()[1] * right.value + + left.first[0] * right.first[1] + + left.value * right.second.ReadData()[1]; + } + const int exceptions = std::fetestexcept(FE_ALL_EXCEPT); + const unsigned status = _statusfp(), mxcsr = _mm_getcsr(); + const auto after = FpState(); + out += "product " + std::to_string(kind) + " initial-fp=" + before + + " final-fp=" + after + " " + + Snapshot("result", product, exceptions, status, mxcsr); + if (left_before != Snapshot("input", left, 0, 0, 0) || + right_before != Snapshot("input", right, 0, 0, 0)) + throw std::runtime_error("supplemental product modified input"); + + // Independently pin the actual MKL operand-swap/transpose contract. + // Capture FP state before inspecting values (sNaN classification can + // itself raise INVALID). Each call starts from the same environment. + Jet forward{}, reverse{}; + reset(); + Outer(left.first, right.first, forward.second.OverwriteData()); + const auto forward_fp = FpState(); + reset(); + Outer(right.first, left.first, reverse.second.OverwriteData()); + const auto reverse_fp = FpState(); + for (std::size_t row = 0; row < 24; ++row) { + for (std::size_t column = 0; column < 24; ++column) { + const double a = forward.second.ReadData()[row * 24 + column]; + const double b = reverse.second.ReadData()[column * 24 + row]; + if (!(std::isnan(a) && std::isnan(b)) && Bits(a) != Bits(b)) + throw std::runtime_error("outer transpose coefficient mismatch"); + } + } + out += "forward-outer-fp=" + forward_fp + + " reverse-outer-fp=" + reverse_fp + "\n"; + } + } + } + } + ResetFp(); + return out; +} + +// Public request-mode evidence for later direction-channel separation. Existing +// tangent tests check nonzero-direction derivatives but do not freeze the full +// omitted-direction response inventory and invalid-direction diagnostics. +std::string PublicSnapshots() { + std::array nodes; + const std::array, 4> coordinates{ + {{{-1, -1, 0}}, {{1, -1, 0}}, {{1, 1, 0}}, {{-1, 1, 0}}}}; + std::array pointers{}; + std::array frames{}; + NonlinearMitc4TrialDofs trial{}; + Vector zero(24), direction(24), wrong(23), nan(24), inf(24); + for (std::size_t i = 0; i < 4; ++i) { + nodes[i] = {{"Jet-Baseline", static_cast(i + 1), + std::to_string(i + 1)}, + coordinates[i], + {"jet-baseline.inp", i + 1}}; + pointers[i] = &nodes[i]; + frames[i] = {static_cast(i), {0, 0, 1}, {1, 0, 0}, {0, 1, 0}}; + trial.displacement[i] = {0.01 * i, -0.02 * i, 0.015 * i}; + trial.rotation_vector[i] = {0.02 * i, -0.01 * i, 0.005 * i}; + } + for (std::size_t i = 0; i < 24; ++i) + direction[i] = (static_cast(i % 5) - 2) / 8; + nan[7] = std::numeric_limits::quiet_NaN(); + inf[7] = std::numeric_limits::infinity(); + const auto created = NonlinearMitc4Shell::Create( + pointers, frames, {"Section", 0.2, 0, {"jet-baseline.inp", 20}}, + {"Material", 120, 0.25, {"jet-baseline.inp", 21}}); + if (!created.HasValue()) throw std::runtime_error("public fixture creation"); + std::ostringstream out; + for (const auto mode : {Mitc4DrillingMode::kDirectSpring, + Mitc4DrillingMode::kLegacyCorotatedTwist}) { + const auto shell = + mitc4_drilling_internal::TestAccess::WithMode(created.Value(), mode); + auto observe = [&](const char* name, const NonlinearMitc4TrialDofs& input, + const Vector* request, + std::optional expected_success) { + const auto before = ResetFp(); + const auto result = + EvaluateInternalResidualAndTangent(shell, input, request); + const int exceptions = std::fetestexcept(FE_ALL_EXCEPT); + const unsigned status = _statusfp(), mxcsr = _mm_getcsr(); + if (expected_success && result.HasValue() != *expected_success) + throw std::runtime_error(std::string("public outcome: ") + name); + out << "public " << static_cast(mode) << ' ' << name << ' ' + << result.HasValue() << ' ' << exceptions << ' ' << status << ' ' + << mxcsr << " initial-fp=" << before << " final-fp=" << FpState(); + auto put = [&](double value) { + if (!std::isfinite(value)) + throw std::runtime_error("public nonfinite field"); + out << ' ' << std::hex << Bits(value) << std::dec; + }; + auto vector = [&](const Vector& value) { + out << " v" << value.Size(); + for (std::size_t i = 0; i < value.Size(); ++i) put(value[i]); + }; + auto matrix = [&](const Matrix& value) { + out << " m" << value.Rows() << 'x' << value.Columns(); + for (std::size_t i = 0; i < value.Rows() * value.Columns(); ++i) + put(value.Data()[i]); + }; + if (result.HasValue()) { + const auto& r = result.Value(); + vector(r.internal_residual); + matrix(r.internal_tangent); + vector(r.internal_force); + put(r.physical_strain_energy); + vector(r.physical_force); + vector(r.drilling_force); + matrix(r.drilling_tangent); + matrix(r.drilling_metric); + put(r.drilling_energy); + vector(r.force_directional_second); + vector(r.drilling_force_directional_second); + vector(r.drilling_constraint); + matrix(r.drilling_constraint_gradient); + for (const auto& hessian : r.drilling_constraint_hessian) + matrix(hessian); + for (const auto& third : + r.drilling_constraint_third_directional_gradient) + vector(third); + put(r.recovery.physical_strain_energy); + for (const auto& point : r.recovery.gauss_points) { + out << " location=" << static_cast(point.location); + for (double value : point.natural_coordinates) put(value); + for (const auto& axis : {point.local_frame.e1, point.local_frame.e2, + point.local_frame.e3}) + for (double value : axis) put(value); + for (double value : point.generalized_strain) put(value); + for (double value : point.section_resultant) put(value); + for (const auto& stress : point.stress) { + out << " section=" << static_cast(stress.position); + put(stress.zeta); + for (double value : stress.components) put(value); + } + } + } else { + const auto& failure = result.GetStatus(); + out << " category=" + << (failure.Category() ? static_cast(*failure.Category()) + : -1); + for (const auto& diagnostic : failure.Diagnostics()) + out << ' ' << static_cast(diagnostic.severity) << ' ' + << std::quoted(diagnostic.code) << ' ' + << std::quoted(diagnostic.location.file.generic_string()) << ' ' + << diagnostic.location.line << ' ' + << std::quoted(diagnostic.keyword) << ' ' + << std::quoted(diagnostic.entity_identity) << ' ' + << std::quoted(diagnostic.message); + } + out << '\n'; + }; + observe("none", trial, nullptr, true); + observe("zero", trial, &zero, true); + observe("nonzero", trial, &direction, true); + observe("invalid-dimension", trial, &wrong, false); + observe("nan-direction", trial, &nan, false); + observe("inf-direction", trial, &inf, false); + auto invalid = trial; + invalid.displacement[1] = {std::numeric_limits::quiet_NaN(), 0, 0}; + observe("nan-trial-none", invalid, nullptr, false); + observe("none-after-failure", trial, nullptr, true); + auto requests = [&](const std::string& name, + const NonlinearMitc4TrialDofs& input) { + observe((name + "-none").c_str(), input, nullptr, std::nullopt); + observe((name + "-zero").c_str(), input, &zero, std::nullopt); + observe((name + "-nonzero").c_str(), input, &direction, std::nullopt); + }; + requests("zero-rotation", NonlinearMitc4TrialDofs{}); + for (const double threshold : {1e-8, 0.25}) { + const double angle = std::sqrt(threshold); + const std::array angles{ + std::nextafter(angle, 0.0), angle, + std::nextafter(angle, std::numeric_limits::infinity())}; + for (std::size_t i = 0; i < angles.size(); ++i) { + NonlinearMitc4TrialDofs boundary{}; + for (auto& rotation : boundary.rotation_vector) + rotation = {angles[i], 0, 0}; + requests((threshold == 0.25 ? "rotation-squared-quarter-" + : "rotation-squared-small-") + + std::to_string(i), + boundary); + } + } + auto translated = trial; + for (auto& displacement : translated.displacement) + displacement = displacement + Vector3{1e9, -2e9, 3e9}; + requests("large-translation", translated); + } + return out.str(); +} + +__declspec(noinline) double Consume(const Jet& value) { + double result = + value.value + value.directional_first + value.directional_second; + for (double entry : value.first) result += entry; + for (double entry : value.second) result += entry; + for (double entry : value.mixed_first) result += entry; + for (double entry : value.mixed_second) result += entry; + return result + value.constant_derivatives + value.second.IsImplicitZero(); +} + +// T3 evidence is separate from both immutable earlier snapshot formats. These +// calls also exercise the original overload when built against the T1 source. +std::string AccumulationSnapshots(bool reassociate) { + const std::array boundary{ + 0.0, -0.0, 1.0, -1.0, 0.1, -0.3, 1e16, + std::numeric_limits::denorm_min(), + std::numeric_limits::max(), + std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN(), + std::numeric_limits::signaling_NaN()}; + std::string out; + const std::array rounding{FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, + FE_TOWARDZERO}; + for (const int mode : rounding) { + for (const unsigned flush : {0u, 0x8040u}) { + for (const unsigned sticky : {0u, 0x21u}) { + for (std::size_t kind = 0; kind < boundary.size(); ++kind) { + ResetFp(); + Jet left = Variable(boundary[kind], kind, -0.25); + Jet right = Variable(0.5, (kind + 1) % 24, 0.125); + if (kind % 3 == 0) left = Constant(boundary[kind]); + if (kind % 3 == 1) { + left.second.Data()[25] = boundary[kind]; + right.second.Data()[48] = -0.5; + } + left.mixed_first[3] = boundary[kind]; + left.mixed_second[4] = -0.0; + const Jet saved_left = left; + const Jet saved_right = right; + const auto before_left = Snapshot("input", left, 0, 0, 0); + const auto before_right = Snapshot("input", right, 0, 0, 0); + const auto environment = [&] { + ResetFp(); + std::fesetround(mode); + _mm_setcsr((_mm_getcsr() & ~0x8040u) | flush | sticky); + }; + const auto observe = [&](const char* name, auto expression) { + environment(); + const auto before = FpState(); + const Jet result = expression(); + const int exceptions = std::fetestexcept(FE_ALL_EXCEPT); + const unsigned status = _statusfp(), mxcsr = _mm_getcsr(); + out += std::to_string(mode) + ":" + std::to_string(flush) + ":" + + std::to_string(sticky) + ":" + std::to_string(kind) + " " + + before + " " + FpState() + " " + + Snapshot(name, result, exceptions, status, mxcsr); + }; + const JetVector3 a{left, right, Constant(-1)}; + const JetVector3 b{right, Constant(1), left}; + observe("dot", [&] { + if (reassociate) + return a[0] * b[0] + (a[1] * b[1] + a[2] * b[2]); + return Dot(a, b); + }); + for (std::size_t component = 0; component < 3; ++component) + observe("cross", [&] { return Cross(a, b)[component]; }); + const JetMatrix3 matrix{{a, b, a}}; + const JetMatrix3 other{{b, a, b}}; + // Every matrix/vector result is consumed over the whole fixture; + // boundary groups rotate the selected entry to keep evidence bounded. + observe("matrix", [&] { + return Multiply(matrix, other)[kind % 3][(kind / 3) % 3]; + }); + observe("matrix-vector", [&] { + return Multiply(matrix, b)[kind % 3]; + }); + observe("frame", [&] { return FrameComponent(a, matrix, b); }); + observe("fixed-frame", [&] { + return FrameComponent(ConstantVector(Vector3{0.5, -1, 0}), matrix, + ConstantVector(Vector3{-0.0, 0.25, 1})); + }); + observe("energy-accumulation", [&] { + Jet energy{}; + for (std::size_t i = 0; i < 3; ++i) + energy = std::move(energy) + a[i] * b[i]; + return energy; + }); + observe("shared-rvalue", [&] { + Jet value = left; + return std::move(value) + right; + }); + observe("same-object-rvalue", [&] { + Jet value = right; + return std::move(value) + value; + }); + observe("unique-rvalue", [&] { + Jet value = left * right; + return std::move(value) + right; + }); + observe("constant-product", [&] { return Constant(-0.5) * left; }); + observe("scalar-product", [&] { return -0.5 * left; }); + if (before_left != Snapshot("input", saved_left, 0, 0, 0) || + before_left != Snapshot("input", left, 0, 0, 0) || + before_right != Snapshot("input", saved_right, 0, 0, 0) || + before_right != Snapshot("input", right, 0, 0, 0)) + throw std::runtime_error("accumulation input alias modified"); + } + } + } + } + ResetFp(); + return out; +} + +std::string AccumulationAliasSnapshots(bool overwrite_input) { + std::string out; + for (int kind = 0; kind < 4; ++kind) { + Jet value = Variable(kind == 0 ? -0.0 : 0.5, 0, -0.25) * + Variable(0.75, 1, 0.5); + if (kind == 1) value = value - value; + if (kind == 2) value.second.Data()[0] = + std::numeric_limits::max(); + const auto before = Snapshot("live-input", value, 0, 0, 0); + out += Observe("unique-same-object", [&] { + if (overwrite_input) value.first[0] = 17; // Fixture-only negative control. + return std::move(value) + value; + }); + if (before != Snapshot("live-input", value, 0, 0, 0)) + throw std::runtime_error("same-object live operand modified"); + Jet survivor = value; + out += Observe("shared-move", [&] { + return std::move(value) + Constant(-1); + }); + if (before != Snapshot("live-input", survivor, 0, 0, 0)) + throw std::runtime_error("surviving shared operand modified"); + } + return out; +} + +void AccumulationBench() { + const Jet left = Variable(0.5, 0, 0.25) * Variable(0.75, 1, -0.5); + const Jet right = Variable(0.25, 2, 0.125) * Variable(-0.5, 3, 0.25); + const JetVector3 a{left, right, left}, b{right, left, right}; + const JetMatrix3 matrix{{a, b, a}}, other{{b, a, b}}; + const auto e1 = ConstantVector(Vector3{0.5, -1, 0.25}); + const auto e2 = ConstantVector(Vector3{-0.25, 0.5, 1}); + volatile double consumed = 0; + for (int kind = 0; kind < 6; ++kind) { + const auto start = std::chrono::steady_clock::now(); + for (int i = 0; i < 3000; ++i) { + if (kind == 0) consumed = consumed + Consume(Dot(a, b)); + if (kind == 1) { + const auto result = Cross(a, b); + for (const auto& value : result) consumed = consumed + Consume(value); + } + if (kind == 2) { + const auto result = Multiply(matrix, other); + for (const auto& row : result) + for (const auto& value : row) consumed = consumed + Consume(value); + } + if (kind == 3) { + const auto result = Multiply(matrix, b); + for (const auto& value : result) consumed = consumed + Consume(value); + } + if (kind == 4) + consumed = consumed + Consume(FrameComponent(e1, matrix, e2)); + if (kind == 5) { + Jet density{}; + for (std::size_t component = 0; component < 3; ++component) + density = std::move(density) + a[component] * b[component]; + consumed = consumed + Consume(0.5 * density); + } + } + std::cout << "accumulation-bench," << kind << ",3000," + << std::chrono::duration( + std::chrono::steady_clock::now() - start).count() + << ",allocation-included,setup-excluded," << consumed << '\n'; + } +} + +void Bench() { + volatile double consumed = 0; + for (int kind = 0; kind < 8; ++kind) { + Jet left = + kind == 0 ? Constant(0.5) : Variable(0.5, 0, kind % 2 ? 0.0 : 0.25); + Jet right = Variable(0.75, 1, kind == 0 || kind % 2 ? 0.0 : -0.5); + if (kind >= 3) { + left = left * left; + right = right * right; + } + if (kind == 5) left = left - left; + if (kind == 7) { + left = Constant(0.5); + right = Constant(0.75); + } + const auto start = std::chrono::steady_clock::now(); + for (int i = 0; i < 20000; ++i) { + left.value += (i % 2 ? -1.0 : 1.0) / 1048576; + const Jet value = kind == 6 ? Dot(JetVector3{left, right, left}, + JetVector3{right, left, right}) + : left * right; + consumed = consumed + Consume(value); + } + std::cout << "bench," << kind << ",20000," + << std::chrono::duration( + std::chrono::steady_clock::now() - start) + .count() + << ",allocation-included," << consumed << '\n'; + } +} +} // namespace +} // namespace fesa + +int main(int argc, char** argv) { + try { + std::fegetenv(&fesa::initial_environment); + fesa::initial_mxcsr = _mm_getcsr(); + if (argc == 2 && std::string(argv[1]) == "--bench") { + fesa::Bench(); + return 0; + } + if (argc == 2 && std::string(argv[1]) == "--accumulation-bench") { + fesa::AccumulationBench(); + return 0; + } + if (argc >= 3 && (std::string(argv[1]) == "--accumulation-write" || + std::string(argv[1]) == "--accumulation-compare")) { + const bool alias = argc > 3 && std::string(argv[3]) == "--unique-alias"; + auto actual = alias ? fesa::AccumulationAliasSnapshots( + argc > 4 && std::string(argv[4]) == + "--overwrite-input") + : fesa::AccumulationSnapshots( + argc > 3 && std::string(argv[3]) == + "--reassociate"); + if (std::string(argv[1]) == "--accumulation-write") { + if (std::ifstream(argv[2], std::ios::binary)) + throw std::runtime_error("accumulation snapshot already exists"); + std::ofstream output(argv[2], std::ios::binary); + output << actual; + output.close(); + if (!output) throw std::runtime_error("accumulation snapshot write failure"); + } else { + std::ifstream input(argv[2], std::ios::binary); + const std::string expected((std::istreambuf_iterator(input)), {}); + if (expected != actual) + throw std::runtime_error("accumulation snapshot mismatch"); + } + std::cout << argv[1] + << (alias ? " passed; 8 unique/same-object/shared alias cases" + : " passed; 192 boundary/environment groups, Dot/" + "Cross/matrix/vector/frame/accumulation/constants") + << ", 651 channels and FP state\n"; + return 0; + } + if (argc >= 3 && (std::string(argv[1]) == "--product-write" || + std::string(argv[1]) == "--product-compare")) { + auto actual = fesa::ProductSnapshots( + argc > 3 && std::string(argv[3]) == "--omit-transpose"); + if (argc > 3 && std::string(argv[3]) == "--mutate") + actual[actual.find('\n') - 1] ^= 1; + if (std::string(argv[1]) == "--product-write") { + if (std::ifstream(argv[2], std::ios::binary)) + throw std::runtime_error("supplemental snapshot already exists"); + std::ofstream output(argv[2], std::ios::binary); + output << actual; + output.close(); + if (!output) throw std::runtime_error("supplemental snapshot write failure"); + } else { + std::ifstream input(argv[2], std::ios::binary); + const std::string expected((std::istreambuf_iterator(input)), {}); + if (expected != actual) + throw std::runtime_error("supplemental product snapshot mismatch"); + } + std::cout << argv[1] << " passed; 3072 product cases and paired MKL outer " + "checks, 651 channels, input aliases, rounding/FTZ/DAZ/sticky\n"; + return 0; + } + const bool verify_only = argc == 2 && std::string(argv[1]) == "--verify"; + if (argc < 3 && !verify_only) + throw std::runtime_error( + "--verify OR --write|--compare PATH [--mutate|--corrupt-oracle]"); + std::string actual = + fesa::Verify(argc > 3 && std::string(argv[3]) == "--corrupt-oracle"); + actual += fesa::PublicSnapshots(); + if (argc > 3 && std::string(argv[3]) == "--mutate") + actual[actual.find('\n') - 1] ^= 1; + if (verify_only) { + // CTest and ordinary regression do not create baseline files. + } else if (std::string(argv[1]) == "--write") { + if (std::ifstream(argv[2], std::ios::binary)) + throw std::runtime_error("snapshot already exists"); + std::ofstream output(argv[2], std::ios::binary); + output << actual; + output.close(); + if (!output) throw std::runtime_error("snapshot write failure"); + } else if (std::string(argv[1]) == "--compare") { + std::ifstream input(argv[2], std::ios::binary); + const std::string expected((std::istreambuf_iterator(input)), {}); + if (expected != actual) + throw std::runtime_error("snapshot mismatch (bits/provenance/fenv)"); + } else { + throw std::runtime_error("unknown verification mode"); + } + std::cout + << argv[1] + << " passed; independent oracle and 42 private cases " + "(651 floating channels each), 64 public request/error cases\n"; + return 0; + } catch (const std::exception& error) { + std::cerr << error.what() << '\n'; + return 1; + } +} diff --git a/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md new file mode 100644 index 0000000..5dfa0ba --- /dev/null +++ b/tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md @@ -0,0 +1,104 @@ +# Nonlinear MITC4 Jet 기준 검사 + +`fesa_jet_baseline_tests`는 제품의 private 구현을 별도 번역 단위에 포함하고 +일반 solver archive를 연결한다. 제품 source와 공개 API를 계측하지 않는다. +독립 다항식 지수 미분 및 미분 항등식과 원래 경로의 비트 snapshot은 별도 근거다. + +VS x64 환경과 저장소의 승인된 dependency 환경변수에서 실행한다. + +```bat +MSBuild msvc\fesa_jet_baseline_tests.vcxproj /m:2 /p:Configuration=Debug /p:Platform=x64 +.build\native\Debug\tests\fesa_jet_baseline_tests.exe --verify +cmake --preset windows-msvc-ninja-debug +cmake --build --preset windows-msvc-ninja-debug --target fesa_jet_baseline_tests --parallel 2 +ctest --preset windows-msvc-ninja-debug -R "^NonlinearMitc4Jet\.IndependentOracle$" --output-on-failure +``` + +Release는 Native의 `Configuration=Release` 또는 해당 CMake Release preset을 쓴다. +기존 unit/integration/reference 세 GTest 실행물의 발견 identity와 이 전용 검사의 +선택 수는 구분한다. CTest의 `--verify`는 baseline 파일을 생성하지 않는다. + +승인된 기준 revision에서만 새 경로에 `--write `을 한 번 실행한다. +기존 파일은 덮어쓰지 않는다. 이후에는 원래 파일과 hash를 보존하고 같은 빌드 경로와 +구성끼리 비교한다. 일반 회귀나 후보 코드에서 `--write`로 기대값을 만들지 않는다. + +```bat +.build\native\Debug\tests\fesa_jet_baseline_tests.exe --compare "" +.build\native\Release\tests\fesa_jet_baseline_tests.exe --compare "" +.build\cmake-debug\tests\fesa_jet_baseline_tests.exe --compare "" +``` + +`--compare --mutate`와 `--corrupt-oracle`는 각각 snapshot 비교와 독립 +oracle의 fixture 전용 오류 주입이며 exit 1이 기대값이다. `--bench`는 별도 micro +측정이다. 연산의 할당·해제와 교대 입력값 갱신을 포함하고 fixture 생성은 제외한다. +전체 채널을 읽는 noinline checksum과 volatile consume을 포함해 결과 제거를 막는다. +공식 모델 CLI의 시간과 혼합하지 않는다. + +snapshot은 finite 값·signed zero·Inf 비트, NaN 분류, provenance와 implicit Hessian, +복원한 초기 FP 환경 및 연산 뒤 FE/rounding/control word/full MXCSR를 기록한다. +공개 반환의 방향 요청 유무·실패 진단·회복량 identity도 포함한다. 수식의 정확성은 +독립 oracle로 검사하고, snapshot은 기준의 계산 순서와 관측 의미 보존을 검사한다. + +T2 곱셈 경계는 기존 T1 snapshot과 별도 CLI로 보존한다. 원본 제품 source의 hash를 +확인하고 보존한 원본 archive를 연결한 실행물에서만 새 경로에 한 번 생성한다. +후보는 `--product-compare`만 실행한다. 기존 파일을 덮어쓰거나 T1 기대값을 확장하지 않는다. + +```bat + --product-write + --product-compare + --product-compare --mutate + --product-compare --omit-transpose +``` + +마지막 두 명령은 각각 출력 변조와 reverse 외적 항 누락의 fixture 전용 음성 대조로 +exit 1이 기대값이다. 이 그룹은 96개 입력을 4 rounding × 4 FTZ/DAZ × 2 초기 sticky +환경에서 실행한 3072개 곱셈이다. 전체 651 채널, 비대칭·공유·implicit Hessian과 +입력 불변성, provenance, 전체 FP 상태를 비교한다. finite 전용 ±0·subnormal·underflow· +overflow gradient와 비유한 입력 거부를 구분하며, 실제 MKL 두 외적의 operand 교환과 +transpose 성분도 비교한다. 두 호출의 FP 상태는 각각 기록하므로 입력 검사 조기 종료로 +차이가 있는 경우를 수치 동등성으로 숨기지 않는다. NaN payload는 비교하지 않는다. +이 그룹은 예외를 mask한 상태의 관측이다. Unmasked 환경의 최초 trap 순서는 별도 +SEH/process 관측이나 생성 코드 근거가 필요하며, source의 row-major 순서만으로 보장하지 않는다. + +T2의 무가드 단일 외적 후보는 unmasked 최초 INVALID/OVERFLOW 순서가 달라 비채택했다. +채택한 guarded 경로는 MXCSR의 여섯 예외가 모두 masked이고 첫 외적이 성공한 경우에만 +전치 성분과 276개 상삼각 scratch를 재사용한다. Unmasked 환경과 외적 거부는 기존 두 호출, +576개 reverse 외적과 row-major 결합을 유지한다. Masked snapshot과 별도 원본/후보 SEH +관측을 함께 확인하며, 정상 경로의 임시 쓰기 감소를 전체 stack 예약 감소로 보고하지 않는다. + +T5는 서로 배타적인 276개 상삼각 보관과 576개 reverse 외적에 같은 local +`std::array` scratch를 사용한다. 입력·공유 Hessian의 저장소와 수명은 +바꾸지 않는다. 두 분기에서 쓰기 전 읽기가 없는지 검토하고 위 3구성 비교 및 별도 +6개 unmasked 최초 trap을 확인한다. Stack 예약량은 실제 Release production object의 +frame과 `__chkstk`를 확인하며, heap 할당 감소나 전체 모델 속도 향상으로 해석하지 않는다. +별도 복사 source의 Values 생성·소멸 계측은 생존 시작량을 포함한 peak payload와 +증가량을 구분한다. 동시 호출 peak는 스케줄에 따라 달라지고, 계측 시간은 정상 Release +측정에서 제외한다. Hessian 복사 연산 횟수는 implicit/null 복사도 포함하며 strong-ref +증가 횟수나 CPU 비용을 뜻하지 않는다. Global pool이나 trial 간 cache는 도입하지 않는다. + +T3 누적 검사는 기존 T1/T2 snapshot을 확장하지 않고 별도 파일을 쓴다. 원본 source와 +보존 archive의 출처를 확인한 실행물에서만 새 경로에 생성한다. 실제 CMake 경로는 +원본 source staging을 CMake로 configure/build한 실행물과 구분해 연결한다. + +```bat + --accumulation-write + --accumulation-compare + --accumulation-compare --reassociate + --accumulation-write --unique-alias + --accumulation-compare --unique-alias + --accumulation-compare --unique-alias --overwrite-input + --accumulation-bench +``` + +일반 그룹은 12개 경계 × 4 rounding × 2 FTZ/DAZ × 2 초기 sticky의 192개 조합으로 +Dot·Cross·3×3 matrix/vector·일반/고정 frame·energy 형태 누적·상수 Jet/scalar를 +검사한다. 별도 8개 alias 그룹은 유일한 materialized Hessian을 가진 동일 객체, +surviving copy와 공유 rvalue를 검사한다. `--reassociate`와 `--overwrite-input`는 +fixture 전용 음성 대조이며 exit 1이 기대값이다. 모든 비교는 651개 채널·provenance· +implicit Hessian·FP 상태를 포함한다. 상수 Jet와 scalar를 서로 같은 연산으로 가정하지 않는다. + +`--accumulation-bench`는 여섯 연산 종류마다 3,000회, 전체 결과 채널 checksum을 +소비한다. Fixture 준비는 제외하고 연산의 할당·해제와 checksum은 포함한다. 이 일반 +실행물의 시간과 별도 복사 source의 Hessian materialization 계측을 구분한다. +Unmasked 최초 trap은 별도 SEH/process 관측으로, 곱셈 단계와 실제 덧셈 단계의 +검증을 구분한다. 전체 모델 시간·stack 예약 감소를 뜻하지 않는다.