diff --git a/docs/calculation-performance/hex8-performance.md b/docs/calculation-performance/hex8-performance.md new file mode 100644 index 0000000..52eaf39 --- /dev/null +++ b/docs/calculation-performance/hex8-performance.md @@ -0,0 +1,255 @@ +# Hex8 계산 경로 개선 기록 + +이 기록은 성능 요구사항의 Hex8 구현 범위와 개발 검증을 설명한다. +기준 작업 소스는 `8c7c41729eb9eaa904f6739b54d34c37b1bd555f`다. +MKL 경계 수치 소스는 `71e08169ecfc3eee8372cd831b799b87aa47c78b`에서 +검토·검사했으며 이후 `43c747f`는 오류 우선순위 검사 두 개만 추가했다. +`6d709fd`는 trial-local 반복 연산 제거로 별도 검토·검사했다. +후속 `774bfbc`는 geometric 첫 곱을 같은 point 안에서 재사용하며, `6bc5cd3`은 +아래 선택적 paired Hessian의 준비 수명을 바꾼 소스다. `a08df0f`는 동시 호출 +검사의 시작 대기 조건만 추가했고 같은 production을 새 Debug/Release로 검증했다. +후속 병합 `45aac19`에는 root의 공통 finite classifier `47d2a81`이 포함된다. +이 병합 뒤의 소비자 재빌드와 동적 개발 실행은 문서 commit `decd959`에서 수행했다. +이전 검사 결과를 새 소스의 통과로 대입하지 않는다. +공식 기준 Release 소스와 이 개발용 Debug 대조 소스를 혼용하지 않는다. + +## 실제 병목에 대응하는 변경 + +기준 VTune의 `TwoSum → Add → Multiply → hex8_high_low::operator* → +EvaluateHighLow` 호출 경로에서 확인한 pair 산술 비용을 줄였다. +이 단계에서는 공통 high+low 산술을 바꾸지 않고 Hex8 private `hex8_pair_arithmetic.h`에 +두 helper를 추가했다. + +- `AddPair`는 양쪽 입력의 finite·normalized 검사를 먼저 수행한다. Canonical + zero operand만 항등원으로 처리하고 결과 low의 signed zero를 원래 Add와 + 같은 +0으로 만든다. 두 항이 비zero면 원래 Add를 그대로 호출한다. +- `MultiplyPair`는 같은 입력 검사를 유지한다. Low가 zero인 경우에만 정확히 + zero인 cross-product와 반복되는 zero Add를 생략한다. Nonzero component의 + 곱 순서는 원래와 같으며 `TwoProduct`의 least-bit underflow 판정, FMA tail, + overflow 실패를 유지한다. 두 low가 모두 비zero면 원래 Multiply를 호출한다. + +Normalized pair에 zero를 더하는 원래 연산은 값과 비zero low를 보존하며 low의 +zero 부호만 canonicalize한다. 한 번 수행한 뒤 반복해도 같은 결과이므로 이 +항등 연산만 줄일 수 있다. 임의의 작은 값이나 거의 zero인 값을 제거하지 않는다. +Signed zero·비정규화 pair·비유한 값은 shortcut 전에 원래 guard에서 거부한다. +[최종 독립 수치 검토](../../.harness/performance/t3/source-numerical-review-71e0816.md)는 +이 증명과 아래 MKL·오류 수명 경계를 확인했다. 검토 PASS는 실제 실행이나 공식 +전체 성능 합격을 대신하지 않는다. + +동일 trial point에서 `deviatoric_coefficient*fg`, `mu*fg`, `mu*fft`의 첫 곱을 +각각 24·24·9개 준비한다. 원래 576개 tangent entry마다 같은 첫 곱을 반복했던 +경로에서만 재사용하고 두 번째 곱과 덧셈 association은 유지한다. Bulk의 같은 +첫 곱도 row별 24개 준비한다. 다음 trial에서는 새로 계산하며 reference나 다른 +trial의 constitutive 값을 가져오지 않는다. +[이 변경의 독립 검토](../../.harness/performance/t3/source-numerical-review-6d709fd.md)는 +같은 입력·연산으로 생성한 pair의 bit 및 invalid propagation이 유지됨을 확인했다. + +Geometric의 `gradient[a][j]*s[l][j]`도 b에 무관하므로 point별 72개 첫 곱으로 +준비한다. 원래 뒤의 `gradient[b][l]` 곱과 j/l 누적 순서는 그대로다. +이는 nonlinear trial point의 automatic storage 1,152 bytes이며 global 또는 +다른 trial의 cache가 아니다. 새로운 조기 finite 반환이나 clamp를 넣지 않는다. + +## 선택적 paired Hessian의 준비 수명 + +현재 정적 모델은 Float64 `Evaluate`를 사용한다. 개발 static VTune에서 paired +`EvaluateHighLow`는 `Create`가 선택적 zero-Hessian을 미리 준비하는 경로에 +나타났다. 초기 보고에서 static에 `EvaluateSolids`를 연결했던 설명은 잘못이었다. +`EvaluateSolids`는 원래 dynamic profile의 경로이며 이번 static에는 없다. +[소유 경로별 집계](../../.harness/performance/t3/vtune-owner-attribution.json)는 +Create의 descendant sampled CPU 5.15735초, writer 2.750676초, static state +2.626941초를 보였다. 중첩 그룹이며 각 메서드의 header-leaf 자체 비용은 제외되어 +완전한 inclusive 총합이 아니다. 또한 `71e0816` binary 자료를 후속 후보의 비용으로 +소급하지 않는다. + +기존 Create는 선택적 paired 준비가 실패해도 F64 모델을 거부하지 않았다. 따라서 +현재 후보는 모든 필수 Create geometry/material/reference stiffness/finite/positive +diagonal 검사를 유지하면서 이 선택적 준비만 첫 유효 `EvaluateLinearHighLow` +호출로 옮긴다. Float64 static·state·writer는 사용하지 않는 pair 준비를 하지 않는다. +Nonlinear paired 계산은 이 cache의 ready/data를 읽거나 준비하지 않는다. + +준비는 실제 trial의 환경·shape·normalized relative state·전체 current geometry +검사를 통과한 뒤 실행한다. `call_once`의 zero-trial callback은 같은 private +paired quadrature graph를 쓰며 재귀적인 once 진입을 명시적으로 건너뛴다. +Unsupported rounding/FTZ/DAZ 호출은 once에 들어가지 않으므로 환경 복원 뒤 정상 +호출이 준비할 수 있다. 성공한 high와 low Hessian을 모두 저장한 뒤에만 ready를 +설정한다. 정상 trial 입력을 cache에 저장하지 않는다. + +Private prepared-data에는 once flag·ready와 두 matrix만 있다. 같은 owned +reference coordinates/points/material을 복사한 Element끼리 공유한다. Caller나 +Domain 주소에 의존하지 않고 callback의 `this`를 이후에 저장하지 않는다. +성공한 저장량은 24×24 high/low 각 하나, logical value 9,216 bytes와 metadata다. +다른 trial의 cache/history가 누적되지 않는다. Const 첫 호출의 동시성은 once로 +동기화한다. 실제 서로 다른 trial과 copy·move의 첫 호출 검사는 모든 worker가 +시작 대기에 도달한 뒤 함께 해제하며, 원래 전체 pair snapshot도 별도로 대조했다. +Once 횟수와 최대 저장량의 source 구조 증명은 native 결과와 구별한다. +[수명과 동기화의 독립 검토](../../.harness/performance/t3/source-numerical-review-6bc5cd3.md)는 +`6bc5cd3` production 소스를 확인한 자료다. + +## MKL 적용과 연산 순서 + +| 경로 | 일반 연산 / MKL 경계 | 보존한 의미 | +| --- | --- | --- | +| `Hex8Element::Create` | Deviatoric·shear·bulk의 24×24 rank-one outer → `Matrix::Multiply` / DGEMM k=1 | Coefficient를 왼쪽 값에 먼저 곱하고 기존 point·entry 순서로 누적 | +| `Hex8Element::Evaluate` | 같은 24×24 outer와 node/component permutation → DGEMM k=1 | 기존 material/geometric 분리와 최종 tangent 합 | +| `GradientProduct` | Node Gram의 192 entry products → private `MultiplyDiagonal` / DSBMV k=0 | 원래 `(p0+p1)+p2` association | +| `GradientProduct`의 geometric 경로 | 64개 entry의 9 terms, 두 단계 entry products → DSBMV k=0 | 원래 `ga[j]*s[l][j]` 뒤 `gb[l]` 곱, 원래 j/l 순서의 누적 | +| `Hex8Element::ComputeMass` | 8×8 `rho*N_a*N_b` outer → DGEMM k=1 | Weight를 뒤에 곱하고 전체 consistent matrix의 행합을 BC 축소 전에 계산 | + +Gram과 geometric 연산을 일반 GEMM으로 재결합했던 후보는 채택하지 않았다. +중간 `G*Sᵀ`의 큰 합이 원래 term별 연산에는 없는 overflow를 만들 수 있기 +때문이다. 현재 entry products는 MKL을 사용하면서 원래 association과 ordered +sum을 유지한다. Workspace는 단발 192 또는 576개의 연속 값이며 지속 cache나 +global dense matrix를 만들지 않는다. + +공통 Matrix와 private diagonal adapter의 finite 입력 검사를 유지한다. 전역 +thread 설정이나 vendor dispatch, PARDISO 제어, controller·수렴 기준은 바꾸지 +않았다. 작은 MKL 호출의 비용을 아직 공식 전체 실행의 개선 효과로 단정하지 않는다. + +## 실패한 계산과 진단 + +MKL 입력 검증이 이미 실패한 point 또는 bulk block에서만 원래 scalar 식을 +진단용으로 재현한다. 이는 원래 실패 위치와 뒤따르는 gate의 우선순위를 유지하기 +위한 경로다. `backend_inputs_failed`가 남으므로 기존 검사들이 모두 통과해도 +정상 element/response를 반환할 수 없다. 이 경우 명시적인 MKL finite-validation +실패로 종료하며 직접 계산한 성공 결과로 대체하지 않는다. + +Create의 volume·stiffness finite·positive diagonal 순서, Evaluate의 point tensor· +최종 energy·각 row force/tangent 검사를 유지한다. 뒤 point의 geometry나 tensor +실패가 먼저 선택되는 의미도 보존한다. 진단용 candidate는 caller state, accepted +history, 최종 HDF5로 발행하지 않는다. 비유한 operand를 zero로 바꿔 backend를 +실행하거나 tolerance를 완화하지 않는다. + +Density는 Create에서 finite positive로 검증되고 Gauss shape는 불변이며 (0,1)에 +있다. 따라서 mass의 `rho*N_a` 입력은 finite다. Product와 weight 뒤 finite 검사, +consistent 누적 검사, 전체 행합의 positive/finite 검사는 각각 원래대로 남아 있다. + +## 특수 정밀도와 소유권 + +| 파일 / 연산 | 크기·정밀도 | 예외 근거 | +| --- | --- | --- | +| `hex8_high_low.cpp`의 tensor·force·energy·24×24 tangent | Normalized binary64 high+low | Low 상태와 underflow-risk tail 정책을 Float64 BLAS로 대체할 수 없음. Shortcut은 원래 bit oracle와 독립 증명으로 제한 | +| `hex8_geometry.cpp`의 whole-cell geometry certificate | Fixed3 좌표, outward interval | 전체 cell 양성 판정·depth 10·budget 65536을 유지. 변경 없음 | +| `hex8_element.cpp`의 F/S/P/FFT/stress 좌표 tensor | Fixed3×3 좌표·구성 관계 | 좌표 contraction이며 일반 node outer와 구별. 변경 없음 | +| 같은 파일의 displacement 차이와 force contraction | FMA residual 및 보상 누적 | 큰 translation에서 tiny strain을 보존하는 오차항과 association 유지 | +| `solid_static_procedure.cpp`의 uncancelled arithmetic scale | 원래 CSR products의 절대값·오차 인증 누적 | 물리 matrix action이 아니라 기존 linear residual의 roundoff certificate. 변경 없음 | +| MKL 검증 실패 block의 원래 scalar 식 | 실패 진단 전용 | 정상 성공 반환이 금지된 진단 재현. 성공 matrix action의 대체 경로가 아님 | + +Element는 reference coordinates·gradients·weights를 값으로 소유한다. 새 helper나 +MKL workspace는 호출 수명이며 trial 상태를 주소로 재사용하지 않는다. Relative +high+low gradient0은 다른 gradient들의 negative pair sum을 그대로 사용한다. +Mean-volume energy, full reference row-sum mass, point identity, RN64 출력 및 +geometry/trial 검사를 변경하지 않았다. 조립의 deterministic tuple reduction과 +parallel buffer 소유권도 변경하지 않았다. + +## 실제 개발 검증 + +| 범위 | 실제 결과 | 근거 | +| --- | --- | --- | +| 새 MultiplyPair oracle, production은 원래 경로 | Debug configure/build 0, Hex8 46/46 PASS | [Native](../../.harness/performance/t3/before-native.json), [log](../../.harness/performance/t3/before-checks.txt) | +| Production `71e0816` | Debug build 0, Hex8·static·frequency·paired state/kinematics 84/84 PASS, native 0 | [Native](../../.harness/performance/t3/after-native.json), [log](../../.harness/performance/t3/after-checks.txt) | +| `43c747f` tests-only delta | Build 0, overflow 진단 우선순위 2/2 PASS, native 0 | [Native](../../.harness/performance/t3/priority-native.json), [log](../../.harness/performance/t3/priority-checks.txt) | +| `6d709fd` trial-local 반복 연산 제거 | Debug build 0, Hex8 48/48 PASS, native 0 | [Native](../../.harness/performance/t3/cache-native.json), [log](../../.harness/performance/t3/cache-checks.txt) | +| 동일 TU의 원본/현재 Debug 진단 대조 | 60 finite-input 조합, 출력 306행 동일, 각각 native 0 | [Native](../../.harness/performance/t3/diagnostic-native.json), [comparison](../../.harness/performance/t3/diagnostic-comparison.json), [original manifest](../../.harness/performance/t3/original-debug/manifest.json) | +| `a08df0f` 선택적 준비 수명 | Fresh Debug configure/build 0, Hex8·paired state/kinematics·serial/TBB 59/59 PASS, native 0 | [Native](../../.harness/performance/t3/lazy-native.json), [검사](../../.harness/performance/t3/lazy-checks.txt), [dependencies](../../.harness/performance/t3/lazy-dependencies.txt) | +| 같은 TU의 원본/current pair snapshot | Flat/warped, linear/nonlinear 전체 response의 raw bits 9,804행 동일, 양쪽 compile/run native 0 | [Native](../../.harness/performance/t3/pair-snapshot-native.json), [대조](../../.harness/performance/t3/pair-snapshot-comparison.json) | +| `decd959` 공통 finite classifier 통합 | Correct-prefix Debug configure/build/deps 0, HighLowArithmetic·Hex8·paired·serial/TBB 68/68 PASS, native 0 | [Native](../../.harness/performance/t3/postmerge-debug-native.json), [검사](../../.harness/performance/t3/postmerge-debug-checks.txt) | + +Pair oracle는 고정된 원래 산술 header를 사용하여 raw high/low bits를 비교한다. +Signed zero·half-ulp·invalid·최소 subnormal·최대값·overflow·underflow 경계, +모든 power-of-two exponent와 양 부호, 12,000개의 deterministic raw random +입력을 포함한다. 기존 analytical·energy/force/tangent·rigid·warped patch·mass· +geometry·cached failure·state 보존 검사와 함께 실행했다. + +별도 Release 개발 probe는 원본 정적 입력을 각각 새 scratch로 복사하여 워밍업 +없이 소스별 1회 실행했다. 최초 `43c747f`는 9.4669573초, 반복 연산을 줄인 +`6d709fd`는 5.2753221초, 선택적 준비를 지연한 `a08df0f`는 1.7391958초였다. +각각 native 0, identity·finite·완전한 결과 검사와 +기존 전체 U/constrained RF/point S 비교를 통과했다. +[최초 실행](../../.harness/performance/t3/development-static.json), +[후속 실행](../../.harness/performance/t3/development-static-cache.json), +[후속 비교](../../.harness/performance/t3/development-static-cache/run-1/comparison.json), +[선택적 준비 지연 실행](../../.harness/performance/t3/development-static-lazy.json), +[해당 비교](../../.harness/performance/t3/development-static-lazy/run-1/comparison.json)를 +모두 보존한다. 서로 다른 개발 소스의 값이며 빠른 값을 골라 공식 표본으로 쓰지 +않는다. 마지막 개발 값은 2초 이내지만 최종 통합 소스의 공식 목표 합격은 아직 아니다. + +최초 개발 Release binary의 [VTune 함수 자료](../../.harness/performance/t3/vtune-functions.csv)와 +[호출 경로](../../.harness/performance/t3/vtune-callstacks.csv), +[실제 명령과 native 기록](../../.harness/performance/t3/vtune-native.json)을 보존한다. +Software sampling 10ms·stack 수집·hardware characterization 비활성 설정이다. +Collector와 두 export는 native 0이며 target native는 독립 포착하지 않아 null이다. +CPU 17.227초 중 effective 8.093초, spin 8.171초였고 iomp worker 6.605811초, +LowestBitExponent 0.904216초, TwoSum 0.458413초, MultiplyPair 0.364517초였다. +Fpclassify 29곳의 exclusive sampled CPU 합은 1.61345초다. Inclusive 호출 경로의 +시간과 합산하지 않고 프로파일 경과 시간을 공식 목표 판정에 대입하지 않는다. + +원본 Debug library와 candidate를 같은 TU·compiler flags·runtime DLL로 대조했다. +원본 library의 Hex8 production은 `8c7c417`이며 새 standalone helper 검사만 +있던 build다. 이 진단의 컴파일/실행을 CTest 종료 확인 전에 시작하여 약 10초 +CPU 작업이 겹쳤다. [겹침 기록](../../.harness/performance/t3/diagnostic-overlap.json)을 +보존하며 해당 elapsed는 성능 evidence에서 제외한다. 이미 완료된 동일 진단을 +반복하지 않았다. Finite 조합의 동일 결과를 모든 극단 입력의 증명으로 확대하지 않는다. + +기존 Debug/Release Ninja의 MSVC include prefix가 깨졌고 조사한 Hex8 object의 +header dependency 수가 0이었다. [실제 조사](../../.harness/performance/t3/dependency-audit.txt)를 +보존한다. 산술 변경의 cpp 직접 재컴파일 로그는 있지만 header-only rebuild의 +정확성은 보장하지 않는다. `1d77e6b`는 private include guard만 고친 별도 commit으로 +산술 본문은 같다. 최종 검증과 공식 Release는 정상 dependency를 갖춘 fresh build가 +필요하며 기존 개발 build를 최종 검증으로 확대하지 않는다. +VSLANG=1033만으로는 현재 설치의 Korean MSVC resource 문제를 해결하지 못했다. +Root가 확인한 실제 prefix hook으로 `a08df0f`의 fresh Debug/Release build를 +구성했다. 각각 Hex8 high+low object의 Ninja dependencies 15개와 private/public +header를 확인했다. [Release 기록](../../.harness/performance/t3/lazy-release-native.json)은 +configure/build와 dependency 검사의 실제 명령을 보존한다. 이번 private class layout의 +변경에는 factory·state·writer 등 모든 소비자를 새로 컴파일해야 하며 옛 object를 +수동 복사하거나 선택적으로 섞지 않는다. + +Pair snapshot의 원본 executable은 `8c7c417` 전체 include tree를 동결하여 원본 +library와 함께 컴파일했다. Current executable은 current header와 fresh current +library를 사용한다. 같은 probe TU와 Debug flags를 쓰되 두 executable의 response만 +대조했다. Current header를 원본 library에 연결하면 `Hex8Element`의 private layout +변경으로 ABI가 맞지 않으므로 이 조합은 사용하지 않았다. + +Root의 별도 `47d2a81` finite bit-classification 소스는 `45aac19`에서 통합했다. +Root의 독립 검토와 62개 검사는 별도 공통 산술 자료다. 이 병합 뒤 올바른 include +prefix metadata를 가진 기존 fresh tree에서 Debug 15개, Release 8개 build edge로 +실제 header 소비자를 재컴파일했다. Hex8 high+low dependencies 15개를 다시 확인했다. +[Release 재빌드](../../.harness/performance/t3/postmerge-release-native.json)의 +configure/build/deps 모두 native 0이다. 이전 개발 시간이나 snapshot을 이 병합 +소스의 결과로 소급하지 않는다. + +다음 동적 개발 실행은 재컴파일한 고정 소스의 원본 `hexa_dynnonlinear.inp`를 +워밍업 없이 한 번 실행하고 전체 경과 시간을 잰다. 제한은 150초이며 완료하면 +기존 전체 101-frame U 비교를 수행한다. 100초 초과 또는 timeout이면 같은 소스로 +30초 제한 VTune software sampling을 한 번 수집하는 범위가 승인되어 있다. +실제 `decd959` 실행은 150.0362372초에서 timeout으로 종료했다. 외부 측정기와 +종료시킨 target native는 각각 1이며 결과 HDF5가 없어 101-frame 비교는 미실행이다. +입력 원본/scratch 바이트와 실행 전후 source/header/binary/build/DLL identity는 +같았다. [실행 자료](../../.harness/performance/t3/development-dynamic-lazy.json)와 +[실제 target 기록](../../.harness/performance/t3/development-dynamic-lazy/run-1/native.json)을 +보존한다. 100초 목표 미달이며 timeout을 완료된 표본으로 처리하지 않는다. + +후속 software sampling은 같은 binary로 30초 한 번만 수집했다. PowerShell이 +VTune의 정상 stderr를 `NativeCommandError`로 처리하여 wrapper native 1로 +끝났고 collector/target의 actual exit는 포착하지 못했다. 해당 값은 unknown이다. +[실패 기록](../../.harness/performance/t3/vtune-dynamic-lazy-wrapper-failure.json)을 +보존하고 수집을 반복하지 않았다. 기존 raw 자료만 finalize/export했으며 두 export와 +그 wrapper는 native 0이다. 수집과 export 후 관련 프로세스가 모두 끝났음을 확인했다. +[함수 자료](../../.harness/performance/t3/vtune-dynamic-lazy-functions.csv), +[호출 경로](../../.harness/performance/t3/vtune-dynamic-lazy-callstacks.csv), +[export native](../../.harness/performance/t3/vtune-dynamic-lazy-export-native.json)는 +부분 구간의 개발 profile이며 전체 실행 시간 자료가 아니다. + +최신 sampled CPU는 LowestBitExponent 51.037592초, TwoSum 27.917520초, +TwoProduct 13.335204초, Multiply 10.732785초, MultiplyPair 7.678654초다. +LowestBitExponent에서 TwoProduct·Multiply·MultiplyPair를 거쳐 nonlinear +EvaluateHighLow와 EvaluateSolids로 이어지는 경로가 확인된다. 여러 inline 함수가 +같은 주소의 서로 다른 attribution으로 표시되어 IsFiniteBinary64의 행들을 단순 +합산하지 않는다. 이전 source/profile의 sampled CPU와 이 부분 자료를 직접 +성능 비율로 비교하지 않는다. 공통 최하위 bit 검색의 정확한 같은 결과를 유지하면서 +비용을 줄이는 후보를 root 소유 범위의 후속 검토로 전달했다. + +84개 실행과 추가 두 검사 실행은 서로 다른 검사 범위다. 이를 하나의 최종 전체 +검사 실행으로 표시하지 않는다. 최종 전체 Debug/CTest, Hex8 필수 원본 5개 비교, +동적 101-frame 비교, 최종 VTune과 공식 단일 시간 측정은 아직 이 작업의 개발 +검증으로 완료되지 않았다. diff --git a/include/fesa/elements/hex8_element.h b/include/fesa/elements/hex8_element.h index 2cab0e2..ebc6a7b 100644 --- a/include/fesa/elements/hex8_element.h +++ b/include/fesa/elements/hex8_element.h @@ -2,6 +2,7 @@ #define FESA_ELEMENTS_HEX8_ELEMENT_H_ #include +#include #include #include "fesa/elements/element.h" @@ -9,6 +10,9 @@ #include "fesa/results/solid_result_records.h" namespace fesa { +namespace hex8_internal { +struct LinearTangentCache; +} // namespace hex8_internal /// @brief Owns all local quantities from one validated linear or TL trial. struct Hex8Response { @@ -90,7 +94,8 @@ class Hex8Element final : public Element { bool nonlinear) const; /// @brief Uses one paired gradient for energy, force and analytic Hessian. Result EvaluateHighLow( - const Vector& high, const Vector& low, bool nonlinear) const; + const Vector& high, const Vector& low, bool nonlinear, + bool prepare_linear_cache = true) const; /// @brief Attaches stable source context to a model failure. Status Failure(const char* message) const; @@ -103,11 +108,9 @@ class Hex8Element final : public Element { SourceLocation location_; std::array points_{}; Matrix stiffness_{24U, 24U}; - // Optional immutable Hessian from the same paired quadrature graph. Failure - // to build it does not reject models supported by the original F64 paths. - bool paired_linear_tangent_ready_{false}; - Matrix paired_linear_tangent_high_{24U, 24U}; - Matrix paired_linear_tangent_low_{24U, 24U}; + // Optional immutable paired Hessian, prepared once on the first valid paired + // linear trial. Copies share only their identical owned reference definition. + std::shared_ptr paired_linear_tangent_; }; } // namespace fesa diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index b715f5d..71de886 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -1,10 +1,14 @@ #include "fesa/elements/hex8_element.h" #include +#include #include +#include #include +#include "../math/dense_blas_internal.h" #include "hex8_geometry.h" +#include "hex8_linear_tangent_cache.h" namespace fesa { namespace { @@ -26,6 +30,51 @@ std::array Pack(const Tensor& value) { double Dot(const Triple& a, const Triple& b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } +/// @brief Forms a scaled 24-component outer product through dense MKL. +Matrix Outer(const std::array& values, double coefficient) { + Matrix left{24U, 1U}; + Matrix right{1U, 24U}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) { + left(3U * a + i, 0) = coefficient * values[a][i]; + right(0, 3U * a + i) = values[a][i]; + } + return left.Multiply(right); +} +/// @brief Uses MKL entry products and the original ordered Gram/stress sums. +Matrix GradientProduct(const std::array& gradients, + const Tensor* stress = nullptr) { + const std::size_t terms = stress == nullptr ? 3U : 9U; + Vector left{64U * terms}; + Vector right{64U * terms}; + Vector final_factor{64U * terms}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t b = 0; b < 8; ++b) + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < (stress == nullptr ? 1U : 3U); ++l) { + const auto entry = (8U * a + b) * terms + + (stress == nullptr ? j : 3U * j + l); + left[entry] = gradients[a][j]; + right[entry] = stress == nullptr ? gradients[b][j] : (*stress)[l][j]; + final_factor[entry] = gradients[b][l]; + } + auto products = dense_blas_internal::MultiplyDiagonal(left, right); + if (stress != nullptr) + products = dense_blas_internal::MultiplyDiagonal(products, final_factor); + Matrix result{8U, 8U}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t b = 0; b < 8; ++b) { + const auto first = (8U * a + b) * terms; + if (stress == nullptr) { + result(a, b) = (products[first] + products[first + 1U]) + + products[first + 2U]; + } else { + for (std::size_t term = 0; term < terms; ++term) + result(a, b) += products[first + term]; + } + } + return result; +} /// @brief Retains product and cancellation errors in fixed-order contractions. void AccumulateProduct(double left, double right, double& sum, double& correction) { @@ -69,7 +118,9 @@ Hex8Element::Hex8Element(EntityIndex index, ElementDofLayout layout, kappa_{kappa}, mu_{mu}, density_{density}, - location_{std::move(location)} {} + location_{std::move(location)}, + paired_linear_tangent_{ + std::make_shared()} {} Status Hex8Element::Failure(const char* message) const { return Status::Failure(FailureCategory::kModel, @@ -149,6 +200,7 @@ Result Hex8Element::Create(EntityIndex index, // Zero Hessian of the same reference-relative mean-volume energy. std::array integrated_gradients{}; double volume = 0.0; + bool backend_inputs_failed = false; for (const auto& point : candidate.points_) { auto gradients = point.gradients; gradients[0] = {}; @@ -158,30 +210,60 @@ Result Hex8Element::Create(EntityIndex index, for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) integrated_gradients[a][i] += point.weight * gradients[a][i]; + Matrix deviatoric{0U, 0U}; + Matrix shear{0U, 0U}; + Matrix gram{0U, 0U}; + bool point_products_failed = false; + try { + deviatoric = Outer(gradients, -(2.0 / 3.0) * candidate.mu_); + shear = Outer(gradients, candidate.mu_); + gram = GradientProduct(gradients); + } catch (const std::invalid_argument&) { + backend_inputs_failed = true; + point_products_failed = true; + } for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { - const auto& ga = gradients[a]; - const auto& gb = gradients[b]; - const double dot = Dot(ga, gb); for (std::size_t i = 0; i < 3; ++i) for (std::size_t k = 0; k < 3; ++k) { - candidate.stiffness_(3U * a + i, 3U * b + k) += - point.weight * (-(2.0 / 3.0) * candidate.mu_ * ga[i] * gb[k] + - candidate.mu_ * ga[k] * gb[i] + - (i == k ? candidate.mu_ * dot : 0.0)); + if (point_products_failed) { + // Diagnostic replay only: this candidate can never succeed. + const auto& ga = gradients[a]; + const auto& gb = gradients[b]; + candidate.stiffness_(3U * a + i, 3U * b + k) += + point.weight * + (-(2.0 / 3.0) * candidate.mu_ * ga[i] * gb[k] + + candidate.mu_ * ga[k] * gb[i] + + (i == k ? candidate.mu_ * Dot(ga, gb) : 0.0)); + } else { + candidate.stiffness_(3U * a + i, 3U * b + k) += + point.weight * (deviatoric(3U * a + i, 3U * b + k) + + shear(3U * a + k, 3U * b + i) + + (i == k ? candidate.mu_ * gram(a, b) : 0.0)); + } } } } if (!std::isfinite(volume) || volume <= 0.0) return Result::Failure( candidate.Failure("Hex8 reference volume is nonfinite.")); - for (std::size_t a = 0; a < 8; ++a) - for (std::size_t b = 0; b < 8; ++b) - for (std::size_t i = 0; i < 3; ++i) - for (std::size_t k = 0; k < 3; ++k) - candidate.stiffness_(3U * a + i, 3U * b + k) += - (candidate.kappa_ / volume) * integrated_gradients[a][i] * - integrated_gradients[b][k]; + Matrix bulk{0U, 0U}; + try { + bulk = Outer(integrated_gradients, candidate.kappa_ / volume); + } catch (const std::invalid_argument&) { + backend_inputs_failed = true; + bulk = Matrix{24U, 24U}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t b = 0; b < 8; ++b) + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t k = 0; k < 3; ++k) + bulk(3U * a + i, 3U * b + k) = + (candidate.kappa_ / volume) * integrated_gradients[a][i] * + integrated_gradients[b][k]; + } + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) + candidate.stiffness_(i, j) += bulk(i, j); for (std::size_t i = 0; i < 24; ++i) for (std::size_t j = 0; j < 24; ++j) if (!std::isfinite(candidate.stiffness_(i, j))) @@ -192,13 +274,9 @@ Result Hex8Element::Create(EntityIndex index, return Result::Failure( candidate.Failure("Hex8 reference stiffness has an unrepresentable " "positive diagonal.")); - const Vector zero{24U}; - const auto paired = candidate.EvaluateHighLow(zero, zero, false); - if (paired.HasValue()) { - candidate.paired_linear_tangent_high_ = paired.Value().high.tangent; - candidate.paired_linear_tangent_low_ = paired.Value().tangent_low; - candidate.paired_linear_tangent_ready_ = true; - } + if (backend_inputs_failed) + return Result::Failure(candidate.Failure( + "Hex8 dense MKL operands failed finite validation.")); return Result::Success(std::move(candidate)); } @@ -220,17 +298,30 @@ Result Hex8Element::ComputeMass() const { return Result::Failure( Failure("Hex8 mass requires reference density.")); Matrix consistent{24U, 24U}; - for (const auto& point : points_) + for (const auto& point : points_) { + Matrix left{8U, 1U}; + Matrix right{1U, 8U}; + for (std::size_t a = 0; a < 8; ++a) { + left(a, 0) = *density_ * point.shape[a]; + right(0, a) = point.shape[a]; + } + Matrix products{0U, 0U}; + try { + products = left.Multiply(right); + } catch (const std::invalid_argument&) { + return Result::Failure( + Failure("Hex8 reference mass is nonfinite.")); + } for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { - const double value = - *density_ * point.shape[a] * point.shape[b] * point.weight; + const double value = products(a, b) * point.weight; if (!std::isfinite(value)) return Result::Failure( Failure("Hex8 reference mass is nonfinite.")); for (std::size_t i = 0; i < 3; ++i) consistent(3U * a + i, 3U * b + i) += value; } + } Matrix mass{24U, 24U}; for (std::size_t i = 0; i < 24; ++i) for (std::size_t j = 0; j < 24; ++j) { @@ -345,6 +436,7 @@ Result Hex8Element::Evaluate(const Vector& u, Failure("Hex8 mean volume is nonfinite.")); response.strain_energy = (0.5 * kappa_ * mean_trace) * (volume * mean_trace); std::array integrated_fg{}; + bool backend_inputs_failed = false; for (std::size_t g = 0; g < 8; ++g) { const auto& point = points_[g]; const auto& f = deformations[g]; @@ -382,6 +474,30 @@ Result Hex8Element::Evaluate(const Vector& u, for (std::size_t i = 0; i < 3; ++i) integrated_fg[a][i] += point.weight * Dot(f[i], gradients[g][a]); } + Matrix deviatoric{0U, 0U}; + Matrix shear{0U, 0U}; + Matrix gram{0U, 0U}; + Matrix geometric{0U, 0U}; + bool point_products_failed = false; + if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(p) || + !Finite(fft) || !Finite(stress)) + return Result::Failure( + Failure("Hex8 trial tensor is nonfinite.")); + if (nonlinear) { + std::array fg{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) + fg[a][i] = Dot(f[i], gradients[g][a]); + try { + deviatoric = Outer(fg, -(2.0 / 3.0) * mu_); + shear = Outer(fg, mu_); + gram = GradientProduct(gradients[g]); + geometric = GradientProduct(gradients[g], &s); + } catch (const std::invalid_argument&) { + backend_inputs_failed = true; + point_products_failed = true; + } + } for (std::size_t a = 0; a < 8; ++a) { const auto& ga = gradients[g][a]; for (std::size_t i = 0; i < 3; ++i) { @@ -400,42 +516,60 @@ Result Hex8Element::Evaluate(const Vector& u, if (nonlinear) for (std::size_t b = 0; b < 8; ++b) { const auto& gb = gradients[g][b]; - double geometric = 0.0; - for (std::size_t j = 0; j < 3; ++j) - for (std::size_t l = 0; l < 3; ++l) - geometric += ga[j] * s[l][j] * gb[l]; + double diagnostic_geometric = 0.0; + if (point_products_failed) + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < 3; ++l) + diagnostic_geometric += ga[j] * s[l][j] * gb[l]; for (std::size_t i = 0; i < 3; ++i) for (std::size_t k = 0; k < 3; ++k) { const auto row = 3U * a + i; const auto column = 3U * b + k; - response.material_tangent(row, column) += - point.weight * - (-(2.0 / 3.0) * mu_ * Dot(f[i], ga) * Dot(f[k], gb) + - mu_ * Dot(f[i], gb) * Dot(f[k], ga) + - mu_ * fft[i][k] * Dot(ga, gb)); + if (point_products_failed) { + // Reproduce original diagnostic priority, never acceptance. + response.material_tangent(row, column) += + point.weight * + (-(2.0 / 3.0) * mu_ * Dot(f[i], ga) * Dot(f[k], gb) + + mu_ * Dot(f[i], gb) * Dot(f[k], ga) + + mu_ * fft[i][k] * Dot(ga, gb)); + } else { + response.material_tangent(row, column) += + point.weight * + (deviatoric(row, column) + + shear(3U * b + i, 3U * a + k) + + mu_ * fft[i][k] * gram(a, b)); + } if (i == k) response.geometric_tangent(row, column) += - geometric * point.weight; + (point_products_failed ? diagnostic_geometric + : geometric(a, b)) * point.weight; } } } - if (!Finite(f) || !Finite(strain) || !Finite(s) || !Finite(p) || - !Finite(fft) || !Finite(stress)) - return Result::Failure( - Failure("Hex8 trial tensor is nonfinite.")); response.points[g] = {element_index_, static_cast(g + 1U), point.natural, point.position, Pack(stress), Pack(strain)}; } for (std::size_t i = 0; i < 24; ++i) response.internal_force[i] += force_correction[i]; - if (nonlinear) - for (std::size_t a = 0; a < 8; ++a) - for (std::size_t b = 0; b < 8; ++b) - for (std::size_t i = 0; i < 3; ++i) - for (std::size_t k = 0; k < 3; ++k) - response.material_tangent(3 * a + i, 3 * b + k) += - (kappa_ / volume) * integrated_fg[a][i] * integrated_fg[b][k]; + if (nonlinear) { + Matrix bulk{0U, 0U}; + try { + bulk = Outer(integrated_fg, kappa_ / volume); + } catch (const std::invalid_argument&) { + backend_inputs_failed = true; + bulk = Matrix{24U, 24U}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t b = 0; b < 8; ++b) + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t k = 0; k < 3; ++k) + bulk(3U * a + i, 3U * b + k) = + (kappa_ / volume) * integrated_fg[a][i] * integrated_fg[b][k]; + } + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) + response.material_tangent(i, j) += bulk(i, j); + } if (!std::isfinite(response.strain_energy)) return Result::Failure( Failure("Hex8 trial energy is nonfinite.")); @@ -454,6 +588,9 @@ Result Hex8Element::Evaluate(const Vector& u, Failure("Hex8 trial tangent is nonfinite.")); } } + if (backend_inputs_failed) + return Result::Failure( + Failure("Hex8 dense MKL operands failed finite validation.")); return Result::Success(std::move(response)); } } // namespace fesa diff --git a/src/fesa/elements/hex8_high_low.cpp b/src/fesa/elements/hex8_high_low.cpp index e18fe8c..4d377e7 100644 --- a/src/fesa/elements/hex8_high_low.cpp +++ b/src/fesa/elements/hex8_high_low.cpp @@ -2,10 +2,13 @@ #include #include +#include #include #include "../math/high_low_arithmetic.h" #include "hex8_geometry.h" +#include "hex8_linear_tangent_cache.h" +#include "hex8_pair_arithmetic.h" namespace fesa { namespace { @@ -13,12 +16,14 @@ using math_internal::HighLow; using Triple = std::array; using Tensor = std::array; -HighLow operator+(HighLow a, HighLow b) { return math_internal::Add(a, b); } +HighLow operator+(HighLow a, HighLow b) { + return hex8_internal::AddPair(a, b); +} HighLow operator-(HighLow a, HighLow b) { return math_internal::Subtract(a, b); } HighLow operator*(HighLow a, HighLow b) { - return math_internal::Multiply(a, b); + return hex8_internal::MultiplyPair(a, b); } HighLow operator/(HighLow a, HighLow b) { return math_internal::Divide(a, b); @@ -50,7 +55,8 @@ std::array Pack(const Tensor& a) { } // namespace Result Hex8Element::EvaluateHighLow( - const Vector& high, const Vector& low, bool nonlinear) const { + const Vector& high, const Vector& low, bool nonlinear, + bool prepare_linear_cache) const { using math_internal::IsNormalized; if (!math_internal::ArithmeticEnvironmentSupported()) return Result::Failure(Failure( @@ -82,6 +88,19 @@ Result Hex8Element::EvaluateHighLow( message = "Hex8 full high+low Jacobian interval arithmetic is nonfinite."; return Result::Failure(Failure(message)); } + if (!nonlinear && prepare_linear_cache) { + // Trial guards above run before optional reference preparation. The nested + // zero trial bypasses this once entry and uses the original paired graph. + std::call_once(paired_linear_tangent_->preparation, [this] { + const Vector zero{24U}; + const auto paired = EvaluateHighLow(zero, zero, false, false); + if (paired.HasValue()) { + paired_linear_tangent_->high = paired.Value().high.tangent; + paired_linear_tangent_->low = paired.Value().tangent_low; + paired_linear_tangent_->ready = true; + } + }); + } Hex8HighLowResponse response; std::array force{}; @@ -193,7 +212,7 @@ Result Hex8Element::EvaluateHighLow( for (std::size_t i = 0; i < 3; ++i) { force[3U * a + i] += Dot(p[i], gradients[a]) * HighLow{point.weight}; } - if (!nonlinear && paired_linear_tangent_ready_) continue; + if (!nonlinear && paired_linear_tangent_->ready) continue; std::array fg{}; for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) @@ -201,6 +220,24 @@ Result Hex8Element::EvaluateHighLow( for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) integrated_fg[a][i] += fg[a][i] * HighLow{point.weight}; + // Reuse identical first products within this trial point only. The second + // multiplication and the original left association remain unchanged. + std::array deviatoric_fg{}, shear_fg{}; + Tensor shear_fft{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) { + deviatoric_fg[a][i] = deviatoric_coefficient * fg[a][i]; + shear_fg[a][i] = HighLow{mu_} * fg[a][i]; + } + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t k = 0; k < 3; ++k) + shear_fft[i][k] = HighLow{mu_} * fft[i][k]; + std::array gradient_stress{}; + if (nonlinear) + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t j = 0; j < 3; ++j) + for (std::size_t l = 0; l < 3; ++l) + gradient_stress[a][j][l] = gradients[a][j] * s[l][j]; for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) { const auto dot = Dot(gradients[a], gradients[b]); @@ -208,36 +245,41 @@ Result Hex8Element::EvaluateHighLow( if (nonlinear) for (std::size_t j = 0; j < 3; ++j) for (std::size_t l = 0; l < 3; ++l) - geo += gradients[a][j] * s[l][j] * gradients[b][l]; + geo += gradient_stress[a][j][l] * gradients[b][l]; for (std::size_t i = 0; i < 3; ++i) for (std::size_t k = 0; k < 3; ++k) { const auto row = 3U * a + i; const auto column = 3U * b + k; - const auto mat = deviatoric_coefficient * fg[a][i] * fg[b][k] + - HighLow{mu_} * fg[b][i] * fg[a][k] + - HighLow{mu_} * fft[i][k] * dot; + const auto mat = deviatoric_fg[a][i] * fg[b][k] + + shear_fg[b][i] * fg[a][k] + + shear_fft[i][k] * dot; material[row][column] += mat * HighLow{point.weight}; if (i == k) geometric[row][column] += geo * HighLow{point.weight}; } } } - if (nonlinear || !paired_linear_tangent_ready_) + if (nonlinear || !paired_linear_tangent_->ready) { + std::array bulk_fg{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) + bulk_fg[a][i] = bulk_coefficient * integrated_fg[a][i]; for (std::size_t a = 0; a < 8; ++a) for (std::size_t b = 0; b < 8; ++b) for (std::size_t i = 0; i < 3; ++i) for (std::size_t k = 0; k < 3; ++k) material[3 * a + i][3 * b + k] += - bulk_coefficient * integrated_fg[a][i] * integrated_fg[b][k]; + bulk_fg[a][i] * integrated_fg[b][k]; + } if (!IsNormalized(energy)) return Result::Failure( Failure("Hex8 high+low energy is unsupported or nonfinite.")); response.high.strain_energy = energy.high; response.strain_energy_low = energy.low; - if (!nonlinear && paired_linear_tangent_ready_) { - response.high.material_tangent = paired_linear_tangent_high_; - response.high.tangent = paired_linear_tangent_high_; - response.tangent_low = paired_linear_tangent_low_; + if (!nonlinear && paired_linear_tangent_->ready) { + response.high.material_tangent = paired_linear_tangent_->high; + response.high.tangent = paired_linear_tangent_->high; + response.tangent_low = paired_linear_tangent_->low; } for (std::size_t i = 0; i < 24; ++i) { if (!IsNormalized(force[i])) @@ -245,7 +287,7 @@ Result Hex8Element::EvaluateHighLow( Failure("Hex8 high+low force is unsupported or nonfinite.")); response.high.internal_force[i] = force[i].high; response.internal_force_low[i] = force[i].low; - if (!nonlinear && paired_linear_tangent_ready_) continue; + if (!nonlinear && paired_linear_tangent_->ready) continue; for (std::size_t j = 0; j < 24; ++j) { const auto tangent = material[i][j] + geometric[i][j]; if (!IsNormalized(material[i][j]) || !IsNormalized(geometric[i][j]) || diff --git a/src/fesa/elements/hex8_linear_tangent_cache.h b/src/fesa/elements/hex8_linear_tangent_cache.h new file mode 100644 index 0000000..6e5087c --- /dev/null +++ b/src/fesa/elements/hex8_linear_tangent_cache.h @@ -0,0 +1,21 @@ +#ifndef FESA_SRC_FESA_ELEMENTS_HEX8_LINEAR_TANGENT_CACHE_H_ +#define FESA_SRC_FESA_ELEMENTS_HEX8_LINEAR_TANGENT_CACHE_H_ + +#include + +#include "fesa/math/matrix.h" + +namespace fesa::hex8_internal { +/// @brief Shares one immutable paired Hessian among identical element copies. +/// @note Only call_once preparation writes. Linear readers synchronize through +/// the flag; nonlinear trials do not access this optional state. At most two +/// 24-by-24 matrices (9216 value bytes) are retained after successful preparation. +struct LinearTangentCache { + std::once_flag preparation; + bool ready{false}; + Matrix high{0U, 0U}; + Matrix low{0U, 0U}; +}; +} // namespace fesa::hex8_internal + +#endif // FESA_SRC_FESA_ELEMENTS_HEX8_LINEAR_TANGENT_CACHE_H_ diff --git a/src/fesa/elements/hex8_pair_arithmetic.h b/src/fesa/elements/hex8_pair_arithmetic.h new file mode 100644 index 0000000..d2fdd5f --- /dev/null +++ b/src/fesa/elements/hex8_pair_arithmetic.h @@ -0,0 +1,40 @@ +#ifndef FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_ARITHMETIC_H_ +#define FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_ARITHMETIC_H_ + +#include "../math/high_low_arithmetic.h" + +namespace fesa::hex8_internal { +/// @brief Preserves original zero-add normalization without redundant TwoSums. +/// @note Valid nonzero terms retain the original addition association. +inline math_internal::HighLow AddPair(math_internal::HighLow left, + math_internal::HighLow right) { + using namespace math_internal; + if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow(); + if (left.high == 0 && left.low == 0) + return {right.high, right.low == 0 ? 0.0 : right.low}; + if (right.high == 0 && right.low == 0) + return {left.high, left.low == 0 ? 0.0 : left.low}; + return Add(left, right); +} +/// @brief Omits exact-zero cross products without changing nonzero association. +/// @note A retained zero addition performs the original canonicalization and +/// invalid-result check; its repeated application is idempotent. +inline math_internal::HighLow MultiplyPair(math_internal::HighLow left, + math_internal::HighLow right) { + using namespace math_internal; + if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow(); + if (left.low != 0 && right.low != 0) return Multiply(left, right); + auto product = TwoProduct(left.high, right.high); + if (right.low != 0) { + product = AddPair(product, TwoProduct(left.high, right.low)); + } else if (left.low != 0) { + // The original zero addition precedes the nonzero low-high term. + product = AddPair(product, HighLow{}); + } + if (left.low != 0) + product = AddPair(product, TwoProduct(left.low, right.high)); + return AddPair(product, HighLow{}); +} +} // namespace fesa::hex8_internal + +#endif // FESA_SRC_FESA_ELEMENTS_HEX8_PAIR_ARITHMETIC_H_ diff --git a/tests/unit/elements/hex8_element_test.cpp b/tests/unit/elements/hex8_element_test.cpp index 6ce64ae..7c055af 100644 --- a/tests/unit/elements/hex8_element_test.cpp +++ b/tests/unit/elements/hex8_element_test.cpp @@ -811,6 +811,36 @@ TEST(Hex8Element, FiniteInputsThatOverflowDerivedArithmeticFailClosed) { std::string::npos); } +TEST(Hex8Element, BackendOperandRejectionRetainsOriginalStiffnessDiagnostic) { + Coordinates coordinates{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) + coordinates[a][i] = 0.01 * kSigns[a][i]; + // Original native probe: finite moduli/volume, but mu*gradient overflows. + // The failed diagnostic candidate must not become a scalar success fallback. + const auto element = Make(coordinates, 1e308, -0.5); + ASSERT_FALSE(element.HasValue()); + EXPECT_EQ(element.GetStatus().Diagnostics()[0].code, "invalid-hex8"); + EXPECT_EQ(element.GetStatus().Diagnostics()[0].message, + "Hex8 reference stiffness is nonfinite."); +} + +TEST(Hex8Element, TensorFailureRetainsPriorityOverEnergyAndTangentOverflow) { + Coordinates coordinates{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) + coordinates[a][i] = 0.01 * kSigns[a][i]; + const auto element = Make(coordinates, 1e200, -0.5); + ASSERT_TRUE(element.HasValue()); + const auto displacement = Affine(coordinates, + Tensor{{{1e100, 0, 0}, {0, 1e100, 0}, {0, 0, 1e100}}}); + const auto response = element.Value().EvaluateNonlinear(displacement); + ASSERT_FALSE(response.HasValue()); + EXPECT_EQ(response.GetStatus().Diagnostics()[0].code, "invalid-hex8"); + EXPECT_EQ(response.GetStatus().Diagnostics()[0].message, + "Hex8 trial tensor is nonfinite."); +} + // Independent Gauss-Legendre abscissae from polynomial root iteration. This // test does not call the kernel's shapes, Jacobian, inverse or quadrature. std::vector> LegendreRule(unsigned count) { diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp index add510b..91e008d 100644 --- a/tests/unit/elements/hex8_high_low_test.cpp +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -2,16 +2,88 @@ #include +#include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include #include "fesa/materials/isotropic_linear_elastic_material.h" #include "../../../src/fesa/math/high_low_arithmetic.h" #include "../../../src/fesa/elements/hex8_geometry.h" +#include "../../../src/fesa/elements/hex8_pair_arithmetic.h" +#include "../math/high_low_arithmetic_original.h" namespace { using fesa::math_internal::HighLow; using fesa::math_internal::Subtract; +std::uint64_t PairBits(double value) { + std::uint64_t bits = 0; + std::memcpy(&bits, &value, sizeof(bits)); + return bits; +} + +void ExpectOriginalProduct(HighLow left, HighLow right) { + const auto actual = fesa::hex8_internal::MultiplyPair(left, right); + const auto expected = fesa::high_low_original::Multiply( + {left.high, left.low}, {right.high, right.low}); + EXPECT_EQ(PairBits(actual.high), PairBits(expected.high)); + EXPECT_EQ(PairBits(actual.low), PairBits(expected.low)); + const auto actual_sum = fesa::hex8_internal::AddPair(left, right); + const auto expected_sum = fesa::high_low_original::Add( + {left.high, left.low}, {right.high, right.low}); + EXPECT_EQ(PairBits(actual_sum.high), PairBits(expected_sum.high)); + EXPECT_EQ(PairBits(actual_sum.low), PairBits(expected_sum.low)); +} + +TEST(Hex8HighLow, ScalarCrossProductShortcutMatchesOriginalBits) { + ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported()); + const double minimum = std::numeric_limits::denorm_min(); + const double maximum = std::numeric_limits::max(); + const double infinity = std::numeric_limits::infinity(); + const double nan = std::numeric_limits::quiet_NaN(); + std::vector pairs{{}, {-0.0, 0.0}, {0.0, -0.0}, {-0.0, -0.0}, + {1, -0.0}, {-1, -0.0}, {1, 1}, {0, minimum}, {infinity, 0}, + {0, infinity}, {nan, 0}, {0, nan}, {maximum}, {-maximum}, + {minimum}, {-minimum}, {std::numeric_limits::min()}, + {0.5}, {std::nextafter(0.5, 1.0)}, + {1, std::ldexp(1.0, -53)}, + {std::nextafter(1.0, 2.0), -std::ldexp(1.0, -53)}}; + for (const auto left : pairs) + for (const auto right : pairs) ExpectOriginalProduct(left, right); + for (int exponent = -1074; exponent <= 1023; ++exponent) { + const double value = std::ldexp(1.0, exponent); + for (double sign : {-1.0, 1.0}) + for (const auto right : pairs) { + ExpectOriginalProduct(HighLow{sign * value}, right); + ExpectOriginalProduct(right, HighLow{sign * value}); + } + } + std::mt19937_64 random(0x48e81234U); + for (int sample = 0; sample < 12000; ++sample) { + double left = 0; + double right = 0; + const auto left_bits = random(); + const auto right_bits = random(); + std::memcpy(&left, &left_bits, sizeof(left)); + std::memcpy(&right, &right_bits, sizeof(right)); + const auto pair = fesa::math_internal::TwoSum(left, right); + ExpectOriginalProduct(pair, HighLow{}); + ExpectOriginalProduct(HighLow{}, pair); + ExpectOriginalProduct(pair, HighLow{right}); + ExpectOriginalProduct(HighLow{right}, pair); + ExpectOriginalProduct(pair, pair); + } +} + fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) { fesa::ElementDofLayout layout; layout.node_indices = {0, 1, 2, 3, 4, 5, 6, 7}; @@ -27,6 +99,116 @@ fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) { return fesa::Hex8Element::Create(0, layout, coordinates, material, 2).Value(); } +void ExpectSamePairResponse(const fesa::Hex8HighLowResponse& actual, + const fesa::Hex8HighLowResponse& expected) { + EXPECT_EQ(PairBits(actual.high.strain_energy), + PairBits(expected.high.strain_energy)); + EXPECT_EQ(PairBits(actual.strain_energy_low), + PairBits(expected.strain_energy_low)); + ASSERT_EQ(actual.high.tangent.Rows(), 24U); + ASSERT_EQ(actual.tangent_low.Rows(), 24U); + for (std::size_t i = 0; i < 24; ++i) { + EXPECT_EQ(PairBits(actual.high.internal_force[i]), + PairBits(expected.high.internal_force[i])); + EXPECT_EQ(PairBits(actual.internal_force_low[i]), + PairBits(expected.internal_force_low[i])); + for (std::size_t j = 0; j < 24; ++j) { + EXPECT_EQ(PairBits(actual.high.tangent(i, j)), + PairBits(expected.high.tangent(i, j))); + EXPECT_EQ(PairBits(actual.tangent_low(i, j)), + PairBits(expected.tangent_low(i, j))); + } + } + for (std::size_t point = 0; point < 8; ++point) { + EXPECT_EQ(actual.high.points[point].integration_point, + expected.high.points[point].integration_point); + for (std::size_t i = 0; i < 6; ++i) { + EXPECT_EQ(PairBits(actual.high.points[point].cauchy_stress[i]), + PairBits(expected.high.points[point].cauchy_stress[i])); + EXPECT_EQ(PairBits(actual.high.points[point].strain[i]), + PairBits(expected.high.points[point].strain[i])); + } + } +} + +TEST(Hex8HighLow, ConcurrentFirstLinearTrialsAcrossCopiesAndMoveKeepBits) { + const auto element = UnitBrick(); + auto copied = element; + const auto moved = std::move(copied); + fesa::Vector high(24), low_a(24), low_b(24); + for (std::size_t i = 0; i < 24; ++i) high[i] = 1; + for (auto a : {1, 2, 5, 6}) { + low_a[3U * a] = std::ldexp(1.0, -60); + low_b[3U * a] = -std::ldexp(1.0, -59); + } + const auto oracle = UnitBrick(); + const auto expected_a = oracle.EvaluateLinearHighLow(high, low_a); + const auto expected_b = oracle.EvaluateLinearHighLow(high, low_b); + ASSERT_TRUE(expected_a.HasValue()); + ASSERT_TRUE(expected_b.HasValue()); + std::promise start; + const auto signal = start.get_future().share(); + std::atomic arrived{0U}; + std::array>, 4> calls; + for (std::size_t i = 0; i < calls.size(); ++i) + calls[i] = std::async(std::launch::async, [&, i] { + arrived.fetch_add(1U, std::memory_order_release); + signal.wait(); + return (i % 2U == 0U ? element : moved) + .EvaluateLinearHighLow(high, i < 2U ? low_a : low_b); + }); + while (arrived.load(std::memory_order_acquire) != calls.size()) + std::this_thread::yield(); + start.set_value(); + for (std::size_t i = 0; i < calls.size(); ++i) { + const auto response = calls[i].get(); + ASSERT_TRUE(response.HasValue()); + ExpectSamePairResponse(response.Value(), + i < 2U ? expected_a.Value() : expected_b.Value()); + } +} + +TEST(Hex8HighLow, FirstLinearPreparationRecoversAfterTrialAndEnvironmentFailure) { + ASSERT_TRUE(fesa::math_internal::ArithmeticEnvironmentSupported()); + const auto element = UnitBrick(); + fesa::Vector high(24), low(24); + EXPECT_FALSE(element.EvaluateLinearHighLow(fesa::Vector(23), low).HasValue()); + high[0] = std::numeric_limits::infinity(); + EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue()); + high[0] = 0; + struct RestoreEnvironment { + int rounding = std::fegetround(); +#if defined(_M_X64) || defined(__SSE2__) + unsigned csr = _mm_getcsr(); +#endif + ~RestoreEnvironment() { + std::fesetround(rounding); +#if defined(_M_X64) || defined(__SSE2__) + _mm_setcsr(csr); +#endif + } + } restore; + ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0); + const auto bad_environment = element.EvaluateLinearHighLow(high, low); + ASSERT_FALSE(bad_environment.HasValue()); + EXPECT_EQ(bad_environment.GetStatus().Diagnostics()[0].message, + "Hex8 high+low requires nearest rounding and gradual underflow."); + ASSERT_EQ(std::fesetround(FE_TONEAREST), 0); +#if defined(_M_X64) || defined(__SSE2__) + _mm_setcsr(restore.csr | 0x8040U); + EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue()); + _mm_setcsr(restore.csr); +#endif + const auto recovered = element.EvaluateLinearHighLow(high, low); + ASSERT_TRUE(recovered.HasValue()); + high[3] = -2; + EXPECT_FALSE(element.EvaluateLinearHighLow(high, low).HasValue()); + high[3] = 0; + const auto after_failure = element.EvaluateLinearHighLow(high, low); + ASSERT_TRUE(after_failure.HasValue()); + ExpectSamePairResponse(after_failure.Value(), recovered.Value()); +} + TEST(Hex8HighLow, LowOnlyAffineStrainSurvivesUniformHighTranslation) { auto brick = UnitBrick(); fesa::Vector high(24), low(24);