baram2584/FESADev
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

merge: integrate validated Hex8 performance changes

Browse Source
This commit is contained in:
KOKO\Mimi
2026-10-14 18:49:26 +09:00
parent 85ae2eb7be e85b96e651
commit 4365f133e1
8 changed files with 777 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/calculation-performance/hex8-performance.md
+255
View File
@@ -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과 공식 단일 시간 측정은 아직 이 작업의 개발
검증으로 완료되지 않았다.
include/fesa/elements/hex8_element.h
+9 -6
View File
@@ -2,6 +2,7 @@
#define FESA_ELEMENTS_HEX8_ELEMENT_H_
#include <array>
#include <memory>
#include <optional>
#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<Hex8HighLowResponse> 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<ReferencePoint, 8> 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<hex8_internal::LinearTangentCache> paired_linear_tangent_;
};
} // namespace fesa
src/fesa/elements/hex8_element.cpp
+183 -46
View File
@@ -1,10 +1,14 @@
#include "fesa/elements/hex8_element.h"
#include <cmath>
#include <memory>
#include <set>
#include <stdexcept>
#include <utility>
#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<double, 6> 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<Triple, 8>& 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<Triple, 8>& 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<hex8_internal::LinearTangentCache>()} {}
Status Hex8Element::Failure(const char* message) const {
return Status::Failure(FailureCategory::kModel,
@@ -149,6 +200,7 @@ Result<Hex8Element> Hex8Element::Create(EntityIndex index,
// Zero Hessian of the same reference-relative mean-volume energy.
std::array<Triple, 8> 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> 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<Hex8Element>::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> Hex8Element::Create(EntityIndex index,
return Result<Hex8Element>::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<Hex8Element>::Failure(candidate.Failure(
"Hex8 dense MKL operands failed finite validation."));
return Result<Hex8Element>::Success(std::move(candidate));
}
@@ -220,17 +298,30 @@ Result<ElementStiffnessContribution> Hex8Element::ComputeMass() const {
return Result<ElementStiffnessContribution>::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<ElementStiffnessContribution>::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<ElementStiffnessContribution>::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<Hex8Response> Hex8Element::Evaluate(const Vector& u,
Failure("Hex8 mean volume is nonfinite."));
response.strain_energy = (0.5 * kappa_ * mean_trace) * (volume * mean_trace);
std::array<Triple, 8> 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<Hex8Response> 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<Hex8Response>::Failure(
Failure("Hex8 trial tensor is nonfinite."));
if (nonlinear) {
std::array<Triple, 8> 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<Hex8Response> 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<Hex8Response>::Failure(
Failure("Hex8 trial tensor is nonfinite."));
response.points[g] = {element_index_, static_cast<int>(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<Hex8Response>::Failure(
Failure("Hex8 trial energy is nonfinite."));
@@ -454,6 +588,9 @@ Result<Hex8Response> Hex8Element::Evaluate(const Vector& u,
Failure("Hex8 trial tangent is nonfinite."));
}
}
if (backend_inputs_failed)
return Result<Hex8Response>::Failure(
Failure("Hex8 dense MKL operands failed finite validation."));
return Result<Hex8Response>::Success(std::move(response));
}
} // namespace fesa
src/fesa/elements/hex8_high_low.cpp
+57 -15
View File
@@ -2,10 +2,13 @@
#include <array>
#include <cmath>
#include <mutex>
#include <utility>
#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<HighLow, 3>;
using Tensor = std::array<Triple, 3>;
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<double, 6> Pack(const Tensor& a) {
} // namespace
Result<Hex8HighLowResponse> 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<Hex8HighLowResponse>::Failure(Failure(
@@ -82,6 +88,19 @@ Result<Hex8HighLowResponse> Hex8Element::EvaluateHighLow(
message = "Hex8 full high+low Jacobian interval arithmetic is nonfinite.";
return Result<Hex8HighLowResponse>::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<HighLow, 24> force{};
@@ -193,7 +212,7 @@ Result<Hex8HighLowResponse> 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<Triple, 8> fg{};
for (std::size_t a = 0; a < 8; ++a)
for (std::size_t i = 0; i < 3; ++i)
@@ -201,6 +220,24 @@ Result<Hex8HighLowResponse> 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<Triple, 8> 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<Tensor, 8> 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<Hex8HighLowResponse> 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<Triple, 8> 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<Hex8HighLowResponse>::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<Hex8HighLowResponse> 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]) ||
src/fesa/elements/hex8_linear_tangent_cache.h
+21
View File
@@ -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 <mutex>
#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_
src/fesa/elements/hex8_pair_arithmetic.h
+40
View File
@@ -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_
tests/unit/elements/hex8_element_test.cpp
+30
View File
@@ -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<std::array<double, 2>> LegendreRule(unsigned count) {
tests/unit/elements/hex8_high_low_test.cpp
+182
View File
@@ -2,16 +2,88 @@
#include <gtest/gtest.h>
#include <array>
#include <atomic>
#include <cfenv>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <future>
#include <limits>
#include <random>
#include <thread>
#include <utility>
#include <vector>
#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<double>::denorm_min();
const double maximum = std::numeric_limits<double>::max();
const double infinity = std::numeric_limits<double>::infinity();
const double nan = std::numeric_limits<double>::quiet_NaN();
std::vector<HighLow> 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<double>::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<void> start;
const auto signal = start.get_future().share();
std::atomic<std::size_t> arrived{0U};
std::array<std::future<fesa::Result<fesa::Hex8HighLowResponse>>, 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<double>::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);