diff --git a/docs/hex8-structural-analysis/ADR.md b/docs/hex8-structural-analysis/ADR.md index 1c66764..23a21c8 100644 --- a/docs/hex8-structural-analysis/ADR.md +++ b/docs/hex8-structural-analysis/ADR.md @@ -57,3 +57,15 @@ warning/no-op으로 수용한다. Mixed1–6 range의1–3은 적용하고4–6 **영향:** 실제 DOF·constraint·state/output은3성분이다. 비영 회전 BC·모든 회전 CLOAD, DOF범위 밖이나 의미 있는 미지원 BC 옵션은 계속 거부한다. 원본을 수정·preprocess하지 않는다. Parser에서 경고·translational 적용·회전 equation 부재·거부 경계를 행동 테스트로 확인한다. + +## HEX8-ADR-008 — High+low 내부 물리 상태와 Float64 observable + +**승인:** 2026-10-09 사용자는 제한된 단일·두 좌표 및 보상 가산 관찰 실패 이후 +A([내부 정밀도 변경안](precision-proposal.md))를 선택했다. +**선택:** Solid HHT 내부 상태·predictor·힘·rates·수지·gate에 normalized high+low를 +일관되게 사용하고 accepted history/rollback도 전체 상태를 보존한다. 기존 Float64 +backend와 HDF5 field/shape/dtype은 유지하되 공식 값은 전체 물리량의 반올림 observable이다. +**영향:** 기존 Si/exact-zero·strict gate·HHT 계수·기준 M·고정 h/schedule·100회와 외부 +tolerance는 유지한다. Kernel/geometry에 low를 전달하지 않는 update carry와 구분한다. +작은 독립 산술/Hex8·시간수렴·상태 보존 검증과 긍정적 실제 실행 이후 필수 원본101개 +endpoint·HDF5/외부 비교·전체 회귀를 진행한다. 승인과 수렴 완료를 혼동하지 않는다. diff --git a/docs/hex8-structural-analysis/FEATURED_AGENTS.md b/docs/hex8-structural-analysis/FEATURED_AGENTS.md index 2bcaeeb..42794b9 100644 --- a/docs/hex8-structural-analysis/FEATURED_AGENTS.md +++ b/docs/hex8-structural-analysis/FEATURED_AGENTS.md @@ -16,7 +16,8 @@ static/frequency/HHT 절차를 추가한다. 코드 존재와 기능 검증 완 | [reference model](reference-model.md) | read-only 원본 inventory·quantity/key/point/time mapping | | [numerical review](numerical-review.md) | 작성자와 독립된 정식화/controls 검토의 결과 | | [lessons](lessons.md) | 정밀도 손실·수렴 실패와 수정의 근거 및 재사용할 주의점 | -| [precision proposal](precision-proposal.md) | 미채택 내부 정밀도 변경안·결과 의미·승인 전 검증 범위 | +| [precision proposal](precision-proposal.md) | 승인된 내부 정밀도 변경안·결과 의미·구현 전 검증 범위 | +| [precision formulation](precision-formulation.md) | HHT 전체 high+low 상태·동일 gradient 미분·산술 범위·반올림 출력 계약 | 메인이 문서·공유 연결·완료 판단을 소유한다. Kernel owner는 Hex8 local 계산/복구 payload와 직접 consumer를, semantic owner는 Domain/parser/3성분 DOF/load/assembly를 소유한다. diff --git a/docs/hex8-structural-analysis/PRD.md b/docs/hex8-structural-analysis/PRD.md index c3b84d8..f2b2df9 100644 --- a/docs/hex8-structural-analysis/PRD.md +++ b/docs/hex8-structural-analysis/PRD.md @@ -424,6 +424,11 @@ Generalized sparse FEAST의 공개 진입점은 ## 현재 checklist +2026-10-09 사용자 A 승인은 [내부 정밀도 변경안](precision-proposal.md)의 전체 high+low +물리 상태와 Float64 반올림 observable을 FEM-REQ-HEX8-015·017–020에 적용한다. +HHT 식·기준 consistent M·입력·고정 schedule·100회·Si/exact-zero·원래 strict gate 및 +외부 tolerance는 유지한다. 아래 구현·검증 checklist는 이 승인만으로 완료되지 않는다. + - [x] 수식·입출력·numeric controls·원본 mapping 정리 및 독립 numerical review 완료 - [x] Hex8 K0/consistent M/비선형 f/Kmat/Kgeo/Kt·형상·복구 kernel 검증 - [x] Solid Domain/parser/3N DOF/ordered assembly/amplitude 검증 diff --git a/docs/hex8-structural-analysis/formulation.md b/docs/hex8-structural-analysis/formulation.md index 0320733..79da771 100644 --- a/docs/hex8-structural-analysis/formulation.md +++ b/docs/hex8-structural-analysis/formulation.md @@ -4,6 +4,9 @@ 물리 선택은 [ADR](ADR.md)에 있다. 이 문서는 승인된 식의 구현용 유도이며 외부 결과에 맞추어 physics를 바꾸는 권한이 아니다. 수치 controls는 [별도 원본](numeric-controls.md), 출력 tensor/identity는 [I/O](io.md)를 따른다. +2026-10-09 승인된 Solid HHT 내부 정밀도 변경의 구현 계약은 +[precision formulation](precision-formulation.md)이다. 동적 paired 경로의 상대변위 gradient, +전체 상태·계수 연산과 반올림 출력 의미는 해당 문서를 함께 적용한다. ## Shape·기준 geometry·적분 diff --git a/docs/hex8-structural-analysis/io.md b/docs/hex8-structural-analysis/io.md index a95da5d..73a921d 100644 --- a/docs/hex8-structural-analysis/io.md +++ b/docs/hex8-structural-analysis/io.md @@ -7,6 +7,13 @@ FEM-REQ-HEX8-008–011,016–020,022–023다. Public CLI·path/lease/exit는 ## Source 입력과 실패 경계 +Solid HHT의 [승인된 high+low 계산](precision-proposal.md)은 내부 전체 물리 상태를 +accepted history와 다음 trial에 유지한다. 공식 HDF5의 기존 field·shape·Float64 dtype은 +변경하지 않으며 값은 전체 물리량의 round-to-Float64 observable이다. Self-check는 +직렬화가 그 반올림과 일치하는지 확인한다. 내부 gate PASS를 rounded 배열 자체의 +동일 gate PASS로 해석하지 않는다. 로그에 내부 산술을 명시하고 직렬화 비유한 값· +overflow는 실패로 처리해 기존 final을 보존한다. Restart는 지원하지 않는다. + 기존 single Part/identity Instance/set/target subset을 재사용하고 source spelling/location/ instance/label과 internal index를 구분한다. C3D8의8 distinct nodes와 SOLID SECTION의 elset/ material 연결, constant E/nu 및 procedure별 density capability를 검증한다. Blank section row는 diff --git a/docs/hex8-structural-analysis/lessons.md b/docs/hex8-structural-analysis/lessons.md index 63c4bbf..7fde4ba 100644 --- a/docs/hex8-structural-analysis/lessons.md +++ b/docs/hex8-structural-analysis/lessons.md @@ -228,10 +228,29 @@ Production의 principal refinement·history·ULP 탐색은 이 비교에 포함 게이트를 검증한 것이다. TEMP source는 실행 전 복원하고 실제 object를 재컴파일했다. 이 결과는 production carry 도입이나 반복 횟수 확대를 뒷받침하지 않는다. -Low를 물리 변위로 사용하는 [내부 정밀도 변경안](precision-proposal.md)은 미채택이다. -단순 carry와 구분하여 전체 계산 상태 및 Float64 출력의 반올림 의미를 결정해야 한다. -현재 물리 기준이나 원본을 조용히 변경하지 않는다. +이 진단 종료 시점에는 low를 물리 변위로 사용하는 변경안을 채택하지 않았다. +이후 사용자가 2026-10-09 선택 A로 [내부 정밀도 변경안](precision-proposal.md)을 승인했다. +단순 carry와 구분해 전체 계산 상태를 high+low로 유지하고 Float64는 반올림 observable로 +출력한다. 기존 실패 evidence는 새 방식의 통과 증거가 아니며 물리 기준·원본은 유지한다. 원시 capture·정확 산술·backsolve·이웃 예측·native replay의 재현 자료와 독립 raw 검산은 로컬 `.harness/hex8-evidence/t6-independent-frozen/`에 보존한다. 기능 checklist는 그대로 미완료다. + +## 승인된 전체 high+low 요소의 첫 검증 + +2026-10-09 승인 A 후, low를 누락하는 대조 kernel에서 실제 요소 assertion3개가 실패했다. +전체 pair로 H/Green strain/응력/힘/에너지/접선을 계산한 kernel은 native build0/CTest0이다. +원래 요소27개와 새 검사15개의 최종 묶음42/42가 통과했고, quarter-turn 작은 전단의 힘 +covariance와 독립 exact-dyadic primitive oracle을 포함한다. Native 산술6/6은 RN/FTZ/DAZ +실행 전제와 normal high에 최소 subnormal FMA error가 붙는 경계까지 확인했다. + +동적 선형과 비선형은 같은 상대변위 gradient의 미분을 사용한다. 기존 rounded K0를 새 힘에 +혼합하지 않는다. 기준 절점0의 signed 방향 차분은 접선·에너지 gradient와 일치했고, +두 brick의 공유 면 힘, linear/NL zero-state Hessian 및 low만 남은 현재 셀의 ±폭을 검사했다. +중앙 에너지 차분의 StVK cubic 오차는 step²로 감소한다. 이를 exact 미분값으로 간주해 +과도한 test bound를 적용했던 실패를 구분하고 분석적 truncation bound를 적용했다. + +독립 numerical reviewer는 실제 shared arithmetic·kernel·geometry의 한정된 범위를 PASS로 +판정했다. 재현 로그는 `.harness/hex8-evidence/precision/`에 보존한다. 이 결과는 동적 전체 +상태·시간수렴·원본 endpoint·공식 출력의 완료를 뜻하지 않는다. 현재 checklist는 미완료다. diff --git a/docs/hex8-structural-analysis/numeric-controls.md b/docs/hex8-structural-analysis/numeric-controls.md index ae05468..d866e2c 100644 --- a/docs/hex8-structural-analysis/numeric-controls.md +++ b/docs/hex8-structural-analysis/numeric-controls.md @@ -7,7 +7,11 @@ ## 계산·형상 인증 모든 input/파생 값/결과를 finite 검증하고 positive Jacobian·밀도·material capability를 확인한다. -Production 계산 precision은 Float64다. Whole-cell 인증은 outward interval coefficient bounds와 +정적·고유치와 기존 계산은 Float64다. 사용자 승인 [정밀도 변경안](precision-proposal.md)에 따라 +Solid HHT의 내부 물리 상태와 힘·rates·수지·gate 계산은 normalized high+low를 사용한다. +Float64 sparse backend와 반올림 출력은 유지하며, 전체 상태 gate를 rounded 출력만으로 +재계산한 gate PASS로 표현하지 않는다. 구현 검증은 아직 진행 중이다. +Whole-cell 인증은 outward interval coefficient bounds와 subdivision이며 estimate-only positive 검사와 abs(det J)를 허용하지 않는다. Depth default10, 처리한 cell budget default65536이며 uncertain 시 fail-closed다. 이 값은 resource/precision tuning이고 accepted geometry의 positive 조건을 완화하는 tolerance가 아니다. @@ -40,11 +44,14 @@ h²·da/L, h²·dA/L를 모두 확인한다. Residual/correction limit1e-10, upd max iteration100이 default다. Update scale은 동일 단위의 uncancelled formula term 합이며 scale0은 exact-zero다. 빈 free vector는 vacuous convergence이며 required prescribed/recovery validation은 유지한다. -### 표현 가능한 변위에서의 내부 정제 +### 이전 Float64 구현의 내부 정제 기록 검토된 별도 동적 구현 `f661864`는 위 판정을 유지하며 다음 정제를 사용한다. 관련76개 검사와 독립 코드 검토는 통과했지만 두 동적 원본은 세 번째 증분에서 실패했다. 이전 구현 `b394ae7`의 원본 실패와 변경 근거는 [lessons](lessons.md)에 기록한다. +아래는 이전 경로의 검토 기록이며 새 paired 경로의 실행 지침이 아니다. +새 경로는 [precision formulation](precision-formulation.md)의 전체 상태 Newton을 적용하고, +모든 원래 free 성분을 푼다. 이전 holding/history/ULP 탐색을 자동으로 이식하지 않는다. 물리 잔차가 통과한 좌표에서 sub-ULP 갱신 소실이 발생하면, 물리 잔차와 각 좌표의 fresh correction4개가 모두 통과한 좌표를 hold한다. 다음 반복부터 원래 effective tangent의 diff --git a/docs/hex8-structural-analysis/numerical-review.md b/docs/hex8-structural-analysis/numerical-review.md index 48fdb89..bc8d9a7 100644 --- a/docs/hex8-structural-analysis/numerical-review.md +++ b/docs/hex8-structural-analysis/numerical-review.md @@ -6,6 +6,17 @@ 필수 원본/입력 계약 충돌은2026-10-09 사용자 선택 A로 해소했다. 구현 전 문서 blocker0이며 실제 solver 수치 PASS가 아니다. +## 2026-10-09 — 승인된 high+low의 한정된 검토 + +`hex8_t5_numerical_spec` reviewer가 [precision formulation](precision-formulation.md)과 +실제 shared arithmetic·paired kernel·geometry를 독립 검토했다. Canonical zero와 scalar +−0 lifting의 두 결함은 수정했고 native fixture로 검증했다. 동일 상대변위 gradient의 chain +rule, isotropic material Hessian·geometric 항·Cauchy 복구와 full X+high+low interval은 +일관된다. 산술의 실제 graph bound와 exact-dyadic 기대값을 별도로 검토했으며 native6/6 +산술 검사 로그를 확인했다. 해당 범위 bounded PASS이며 새 actionable finding은 없다. +일반 correctly-rounded106-bit 보장이나 T6 전체 state/시간수렴/원본 PASS가 아니다. +T6 consumer·backend·worker 경계와 rollback은 해당 구현이 고정된 후 별도 검토한다. + ## HEX8-NR-001 — 영 회전 경계와 필수 원본 충돌 검토 당시 PRD008/I/O는 BC/CLOAD DOF4–6을 값0이어도 거부했다. 그러나 모든 필수 INP에 diff --git a/docs/hex8-structural-analysis/precision-formulation.md b/docs/hex8-structural-analysis/precision-formulation.md new file mode 100644 index 0000000..9a799fd --- /dev/null +++ b/docs/hex8-structural-analysis/precision-formulation.md @@ -0,0 +1,280 @@ +# Hex8 동해석 high+low 정식화 + +적용: 사용자가 2026-10-09 선택 A로 승인한 [정밀도 제안](precision-proposal.md). +기본 물리식·DOF·부호·적분점 순서는 [formulation](formulation.md), 판정 원본은 +[numeric-controls](numeric-controls.md), 영속 출력은 [I/O](io.md)다. 본 문서는 +승인된 정밀도 변경의 구현 계약이며 구현 완료나 독립 수치 승인이 아니다. +기존 선형/비선형 정적·고유치의 Float64 경로를 변경하는 계약이 아니다. + +## 산술과 고정 계수 + +내부 scalar \(z=(z_h,z_l)\)의 물리값은 실수 합 \(\widehat z=z_h+z_l\)다. +두 성분은 finite IEEE binary64이며, RN은 round-to-nearest, ties-to-even이다. +정규화는 error-free `TwoSum` 계열의 재분해로 같은 합을 보존하고 +\(z_h=\operatorname{RN}(\widehat z)\)를 만족시킨다. 영은 유일한 \((+0,+0)\)으로 +정규화한다. Nonoverlap 조건은 이 RN 식으로 정의하며 binade 경계에서 임의의 +대칭 ULP 간격을 가정하지 않는다. 입력 binary64 \(c\)는 \((c,0)\)으로 올린다. + +아래 식의 덧셈·곱셈·나눗셈·축약은 별도 표시가 없으면 전체 pair 연산이다. +연산의 목표는 전체 피연산자의 실수 연산값이며 실제 결과에는 pair 산술 반올림이 있다. +덧셈은 high/low 네 항과 취소 오차, 곱셈은 hh/hl/lh/ll 네 곱과 product error를 +포함하고 정규화한다. 나눗셈은 nonzero 전체 분모로 계산하며 high quotient만을 +반환하지 않는다. Quotient refinement의 잔차 \(x-yq\)에도 양쪽 low와 product +error가 참여한다. `TwoSum`/FMA product-error 및 고정 순서의 residual refinement는 +가능한 구현이며 구체적인 storage/backend/vendor 타입은 지정하지 않는다. + +이 표현은 IEEE binary128도, 모든 연산의 correctly-rounded 106-bit 보장도 아니다. +EFT 전제가 성립하는 normal 범위에서 연산 graph의 오차는 same-unit uncancelled operand +scale에 대해 binary64 단일 연산의 1차 오차를 회수한 2차 수준이어야 한다. 취소 후 작은 +결과값만을 relative error 분모로 삼지 않는다. 구현별 graph의 roundoff bound를 독립 exact +dyadic/rational oracle와 대조한다. 최종 잔차가 작다는 사실로 산술 오차 bound를 대신하지 +않는다. Subnormal에는 gradual underflow를 적용하고 EFT의 exact-error 보장이 깨지는 +항은 실제 binary64 반올림 오차로 따로 포함한다. Exponent 범위를 확장하지 않으며 최소 +양의 quantum은 \(2^{-1074}\)다. Underflow로 이 quantum보다 작은 결과가 반올림되는 +것과 임의의 low/잔차 clamp를 구별한다. FTZ/DAZ 또는 재결합이 전제를 깨는 실행은 지원 +산술의 증거가 아니다. 입력·중간·결과의 NaN/Inf, overflow, zero divisor와 정상화 실패는 +candidate 실패다. 유한한 최종 실수값만으로 overflow한 중간 연산을 통과시키지 않는다. + +Underflow로 소실된 비영 항을 exact-zero 평형이나 양의 형상 인증의 근거로 쓰지 않는다. +특히 high product가 정상/비영이라는 검사만으로 hl/lh/ll product와 FMA residual의 +underflow가 배제되지는 않는다. 구현한 graph가 해당 소실의 오차를 보수적으로 추적해 +zero/부호/gate를 보증하지 못하는 입력은 지원 산술 범위 실패다. 모든 EFT intermediate를 +지원한다는 가정 없이 실제 지원 범위와 실패 경계를 독립 검사에 남긴다. + +Exact-zero, 부호, abs, max, 비교와 동일 상태 판별은 정규화한 전체 pair를 사용한다. +비영 low를 지닌 값을 high만 읽어 zero로 판정하지 않는다. Pair 합산 순서는 기존 source, +element, local component 및 CSR 순서로 결정적으로 고정한다. + +### 실제 primitive graph의 오차와 지원 경계 + +구현은 `src/fesa/math/high_low_arithmetic.h`이며 원시 두 성분 입력을 조용히 정상화하지 +않는다. F64 scalar lifting의 ±0만 canonical(+0,+0)로 만든다. Native 실행에서 RN 및 +MXCSR의 FTZ/DAZ 비활성을 검사하고, 지원하지 않는 환경을 자동 변경하지 않는다. + +Add의 head/tail/middle/join `TwoSum`이 EFT인 지원 범위에서 최종 exact 합과의 차이는 +`j=joined.low, m=middle.low, t=tail.low`, `z1=RN(j+m)`, `z2=RN(z1+t)`의 두 반올림뿐이다. +따라서 absolute bound는 `Badd=|z1-(j+m)|+|z2-(z1+t)|`다. Normal 합에 대해서는 +`u(|j+m|+|z1+t|)`, `u=2^-53`으로 상계한다. Subnormal 덧셈은 quantum-grid에서 exact인 +지원 사례를 별도로 검사한다. 취소된 최종 결과만으로 bound를 나누지 않는다. + +TwoProduct는 각 비영 입력의 실제 least dyadic bit exponent 합이 −1074 이상인 경우만 +지원한다. FMA residual까지 표현될 수 있는 보수적 범위이며 high product만 검사하지 않는다. +지원된 hh/hl/lh/ll 네 product는 exact pair이고 Multiply의 absolute bound는 이어지는 +세 Add의 `Badd` 합이다. Overflow 또는 지원되지 않는 exact tail은 invalid로 전파한다. + +Divide는 첫 quotient 후 두 번 full residual refinement를 사용한다. 각 step에서 현재 stored +pair의 실수값을 n,d,q라 하고 `rtilde=Subtract(n,Multiply(d,q))=(rh,rl)`, +`Br=Bmultiply(d,q)+Bsubtract(n,p)`, `c=RN(rh/dh)`로 둔다. `Bq`는 `Add(q,c)`의 bound, +`Bc`는 scalar quotient RN bound이면 새 quotient의 absolute error는 +`Bq+Bc+(|rh*dl|+|rl*dh|)/(|dh|*|d|)+Br/|d|` 이하이다. +분모는 exact pair의 양의 절댓값 또는 `|dh|-|dl|`의 보수적 양의 lower bound를 사용한다. +Normal quotient의 `Bc=u*|rh/dh|`, 지원 subnormal의 exact rational rounding cell bound를 +구별한다. 두 번째 step도 새 full residual을 다시 계산한다. Quotient가 비영 residual을0으로 +소실하거나 중간 연산이 지원 범위를 벗어나면 실패한다. 이 bound는 universal correctly-rounded +106-bit 보장이 아니다. 독립 exact rational fixtures와 binary scaling으로 실제 graph를 검증한다. + +최소 계수 정밀도 선택은 다음과 같다. 기존 reference cache의 F64 \(N_a,G_a,w_g\), +Lamé 상수와 좌표, 기존 조립 F64 consistent \(M\), 기존 load assembler가 endpoint에서 +평가한 F64 \(F_n\)를 exact dyadic 입력으로 올린다. 형상/재료/질량 적분을 새로 정의하거나 +높은 차수 quadrature로 바꾸지 않는다. 단, 아래 상대변위의 미분에 필요한 \(\widehat G_0\)와 +dynamic paired linear operator는 동일 pair 산술로 구성한다. 기존 rounded K0를 paired +힘에 혼합하지 않는다. 일정 F64 \(\alpha,\beta,\gamma,h\)는 기존 생성값을 올리며 +\(c=1+\alpha\), \(b=\beta h^2\)와 모든 파생 계수는 공통 pair 산술로 계산한다. +\(b>0,c>0\) 및 finite를 요구한다. 기존 coefficient 생성·검증과 +`time=double(n)*h`, `N*h==T` 및 endpoint 하중 평가 의미는 유지한다. + +## 같은 상대변위에서의 energy·force·Hessian + +절점 index \(a=0,\ldots,7\), 성분 \(i=1,2,3\), node-major 24 DOF를 사용한다. +\(u,v\)의 단위는 length, length/time이고 \(a,A\)는 length/time²다. +기존 8개 reference point 각각에서 + +\[ + d_a=u_a-u_0,\quad H=\sum_{a=1}^7d_a\otimes G_a, + \qquad \widehat G_a=G_a\ (a>0),\quad + \widehat G_0=-\sum_{a=1}^7G_a. +\] + +모든 차·곱·합은 pair다. \(\delta H=\sum_{a=0}^7\delta u_a\otimes\widehat G_a\)가 +**실제로 평가한 상대변위 식의 미분**이다. Force/tangent에 기존 \(G_0\)를 따로 사용하면 +이 미분이 아니다. F64 cache를 exact dyadic로 보면 기존 \(\sum_aG_a=0\)은 반올림 때문에 +exact일 필요가 없다. \(\widehat G_0\) 선택은 기존 상대변위 weak form의 일관된 수치 평가이며 +잔차 projection이나 사후 힘 보정이 아니다. Pair 자체의 유한 반올림 오차는 남는다. +\(\widehat G_0\) 생성은 local node1부터7까지 고정 순서의 paired sum과 negate다. +수식의 실수 합과 그 pair 근사 오차를 구별하고, 임의 reduction 순서에서 저장된 +\(\sum_a\widehat G_a\)가 exact0이라고 주장하지 않는다. Node0의 미분 오차도 이 graph의 +산술 bound에 포함한다. + +선형 동해석에서는 \(\epsilon=\operatorname{sym}H\), +\(\sigma=\lambda\operatorname{tr}\epsilon I+2\mu\epsilon\)이고 + +\[ + U_{lin}=\sum_g[\tfrac12\lambda(\operatorname{tr}\epsilon)^2+ + \mu\epsilon:\epsilon]w_g, + \quad f_{ai}=\sum_g\sigma_{iJ}\widehat G_{aJ}w_g, + \quad K_{lin}=\sum_g\widehat B_g^TC\widehat B_gw_g. +\] + +\(\widehat B\)는 위 \(\widehat G\)로 만든 기존 engineering-shear B다. +비선형 동해석은 \(F=I+H\)이며 tiny strain의 불필요한 I 취소를 피하는 동등식 + +\[ + E_G=\operatorname{sym}H+\tfrac12H^TH,\quad + S=\lambda\operatorname{tr}(E_G)I+2\mu E_G,\quad P=FS, +\] +\[ + U=\sum_g[\tfrac12\lambda(\operatorname{tr}E_G)^2+\mu E_G:E_G]w_g, + \qquad f_{ai}=\sum_gP_{iJ}\widehat G_{aJ}w_g. +\] + +\(H,E_G,S,P,F,U,f\)와 8점 recovery는 low를 보존한다. H를 만든 직후 F64로 축소하는 +경로는 이 계약을 충족하지 않는다. \(\mathbb C\)는 기본 정식화의 같은 StVK tensor이며 + +\[ + (K_{geo})_{ai,bk}=\sum_g\delta_{ik}\widehat G_{aJ}S_{LJ} + \widehat G_{bL}w_g, +\] +\[ + (K_{mat})_{ai,bk}=\sum_g\widehat G_{aJ}F_{iI}\mathbb C_{IJLN} + F_{kN}\widehat G_{bL}w_g, + \qquad K_t=K_{mat}+K_{geo}. +\] + +각 식은 동일 trial·동일 계수·동일 2×2×2 적분을 사용한다. 에너지의 실수 polynomial +미분을 pair 산술로 평가한 tangent이며 유한 반올림 함수의 미분을 주장하지 않는다. +Zero state의 paired \(K_t=K_{lin}\)이고 geometric 항은0이다. Paired \(K_{lin}\)은 기존 +cached F64 K0와 coefficient-roundoff 수준에서 다를 수 있다. Byte identity나 별도 +rounded K0와의 exact equality를 요구하지 않는다. 차이는 실제 측정하고 기본 linear +weak form·energy gradient·tangent 불변량으로 독립 검토한다. 대칭 강제 평균이나 +잔차 floor로 이 검토를 대신하지 않는다. Strained tangent가 SPD라고 가정하지 않는다. + +## 형상·질량 action·복구 + +현재 절점은 실수값 \(x_a=X_a+u_{ah}+u_{al}\)이다. Whole-cell 인증의 outward-rounded +interval은 세 항 전체를 포함한다. F64 `X+RN(u)`를 먼저 만든 뒤 singleton interval로 +바꾸지 않는다. Interval \(I(X)+I(u_h)+I(u_l)\) 또는 exact-sum expansion의 보수적 +enclosure를 사용하며 이후 determinant, Bernstein 변환·subdivision도 outward rounding한다. +양의 lower bound를 인증하지 못하면 기존 depth10/cell budget65536에서 uncertain failure다. +Low를 누락한 형상을 인증하거나 small low를 무조건 안전한 것으로 간주하지 않는다. +Point \(\det F>0\)도 full tensor에서 확인하며 유한성/부호가 불확실하면 accept하지 않는다. + +비선형 Cauchy \(\sigma=FSF^T/\det F\), Green strain, 선형 Cauchy/infinitesimal strain은 +pair로 복구한다. 8점 global tensor의 실제 shear `[11,22,33,12,13,23]`를 유지한다. +\(M\)은 기준 형상·밀도에서 한 번 조립한 **기존 F64 matrix 그 자체**이고 변형/시간에 따라 +변하지 않는다. 모든 action \((Mz)_i=\sum_j(M_{ij},0)z_j\)는 full z와 pair product/reduction이다. +F64 \(M\)으로 먼저 rounded high action을 계산하고 low를 별도 residual 보정으로 넣지 않는다. + +\[ + T=\tfrac12v^TMv,\quad p=\sum_a(Mv)_a,\quad + L_O=\sum_ax_a\times(Mv)_a,\quad \dot L_O=\sum_ax_a\times(Ma)_a. +\] + +Current position·momentum·reaction·work·energy balance까지 pair로 계산한다. +\(W_{ext,n+1}=W_{ext,n}+\tfrac12(F_n+F_{n+1})^T(u_{n+1}-u_n)\), +balance error는 기존 \(T+U-T_0-U_0-W_{ext}\)다. HHT numerical dissipation을 지우지 않는다. +선형 current moment의 \(\sum u_a\times f_a=O(u^2)\) 해석과 TL 객관성/기준-arm identity는 +기본 정식화를 유지한다. F64 coefficient 오차와 pair 산술 오차를 독립 검사의 scale에서 +구별하고, finite-rotation 결과를 exact0으로 보정하지 않는다. + +## HHT 전체 state와 Newton + +Accepted \(u_n,v_n,a_n,A_n,f_n,F_n,r_n,W_{ext,n}\), predictor 및 candidate는 +low를 포함한다. 초기 \(u_0=v_0=0\), \(A_0=a_0\), +\(M_{ff}a_{0f}=F_f(0+)-f_f(0)\)다. 초기 backend 해도 lift한 뒤 full endpoint로 검증한다. +Constrained \(u,v,a,A\)는 항상 canonical exact0이다. + +\[ + u_p=u_n+hv_n+h^2(\tfrac12-\beta)A_n,\qquad + v_p=v_n+h(1-\gamma)A_n, +\] +\[ + A_{n+1}=(u_{n+1}-u_p)/b,\quad + v_{n+1}=v_p+\gamma hA_{n+1},\quad + a_{n+1}=(A_{n+1}+\alpha a_n)/c. +\] + +Candidate의 rates·element response·mass actions·recovery는 이 동일 전체 상태에서 평가한다. +Physical residual은 \(r=Ma+f-F\), reaction은 constrained projection, free reaction은 exact0이다. +다음 두 weighted 표현을 별도로 계산한다. + +\[ + R=MA_{n+1}+c(f_{n+1}-F_{n+1})-\alpha(f_n-F_n), + \qquad R_{eq}=c r_{n+1}-\alpha r_n, +\] +\[ + R-R_{eq}=M d_A,\quad d_A=A_{n+1}-c a_{n+1}+\alpha a_n, + \qquad K_{eff}=M/b+cK_t. +\] + +마지막 식은 rates의 roundoff defect를 드러낸 항등식이다. Accepted old residual을0으로 +대체하지 않는다. \(R_{eq}\)를 RHS로 쓰는 경우에도 원래 weighted 식과의 일치 gate를 +통과해야 하고 tangent는 같은 weighted 표현의 \(K_{eff}\)다. +Endpoint residual만 사용하면 그 Jacobian은 \(K_{eff}/c\)이며 서로 혼용하지 않는다. +Backend에는 최종 조립한 paired \(K_{eff}\)와 RHS를 RN64로 전달하고 F64 보정량을 pair로 +올려 \(u\leftarrow u+\delta u\)에 누적한다. High-only 상태에서 tangent를 재계산하지 않는다. +이 backend solve는 원래 tangent의 rounded approximation이며 full residual/physical gate를 +대체하지 않는다. Directional tangent 검사는 paired analytic derivative와 rounded backend +matrix의 coefficient-rounding 차이를 따로 평가한다. Backend/nonfinite 실패는 rollback한다. + +Physical criterion은 **모든 원래 free 성분**의 +\(S_i=|(Ma)_i|+|f_i|+|F_i|\), \(\max_i|r_i|/S_i\le10^{-10}\)다. +\(S_i=0\)이면 \(r_i=0\)을 exact로 요구한다. Norm·분모·gate도 full pair다. +Weighted identity scale은 기존 uncancelled 정의를 그대로 pair로 평가한다. +즉 각 row의 \(\sum_j[|M_{ij}A_j|+c|M_{ij}a_j|+|\alpha M_{ij}a_{n,j}|+ +|M_{ij}d_{A,j}|]\), \(c(|f_i|+|F_i|)+|\alpha|(|f_{n,i}|+|F_{n,i}|)\), +\(c|r_i|+|\alpha r_{n,i}|\)의 합이다. \(|R-R_{eq}|/scale\le10^{-12}\), +zero scale은 exact-zero branch이며 finite 검사를 먼저 수행한다. + +Fresh correction의 네 channel은 +\(\delta A=\delta u/b\), \(\delta v=\gamma h\delta A\), +\(\delta a=\delta A/c\)에서 각각 +\(\|\delta u\|_\infty/L,h\|\delta v\|_\infty/L, +h^2\|\delta a\|_\infty/L,h^2\|\delta A\|_\infty/L\le10^{-10}\)다. +동일 full \(\delta u\)와 공통 계수를 쓰며 새 dummy channel은 없다. +Update gate는 위 u/v/A 등식의 좌변과 전개된 각 우변 항의 same-unit abs 합을 scale로 +사용해 defect/scale \(\le10^{-12}\), scale0은 exact-zero다. Increment당100회 제한, +고정 h와 기존 schedule, 빈 free의 vacuous convergence를 유지한다. +이미 존재하는 hold/revisit 제어를 사용한다면 final gate는 모든 원래 free 성분을 검사하고 +동일 상태 판별에 전체 u pair와 mask를 포함한다. Carry·holding으로 physics state의 low를 +대체하지 않는다. + +## Commit·observable·검증의 의미 + +Accepted/trial 복사는 모든 high/low와 이전 residual·누적 work를 보존한다. Geometry, +recovery, finite, update/weighted identity, physical/correction gate가 모두 통과한 candidate만 +commit한다. 실패는 accepted state·predictor 원본·reference cache·기존 최종 파일을 보존한다. +다음 증분은 published F64 배열 대신 accepted 전체 pair를 사용한다. Restart는 지원하지 않는다. + +각 공개 F64 값은 \(\operatorname{RN}(z_h+z_l)\)이며 high만 저장한다는 가정 대신 올바른 +최종 반올림을 검증한다. Schema·shape·source identity는 기존 계약을 유지한다. +Recovery를 rounded u에서 재계산해 published full-recovery observable을 바꾸지 않는다. +내부 gate PASS는 published F64 u/v/a를 재입력해 같은 gate를 통과했다는 주장과 다르다. +Self-check는 schema/identity/finite 및 내부 물리량의 올바른 RN64 serialization을 확인한다. +로그는 내부 high+low 산술과 적용 gate를 명시한다. 출력 overflow/nonfinite는 실패다. + +검증은 다음 반증 가능한 사례를 포함한다. 실행/통과 기록은 이 문서의 범위가 아니다. + +- Exact dyadic/rational oracle: \(x=(1,2^{-60})\), \(y=(1,-2^{-60})\)에서 + \(x-y=2^{-59}\), \(x\,2^{60}-2^{60}=1\), binary-power division·normalization, + cancellation·binade/tie·subnormal·overflow/zero-divisor를 검사한다. High collapse 대조는 + 앞 두 항등식에서 실패해야 한다. 일반 곱/비이진 나눗셈은 실제 graph 오차 bound로 검사한다. +- 작은 brick 및 coupled 두 Hex8: 공통 high 병진과 \(2^{-60}\) 규모의 비균일 low field로 + nonzero H/E/f/U와 이웃 element/global row 전달을 확인한다. 에너지 차분·signed directional + tangent의 step sweep(기준 node0의 독립 perturbation 포함), paired linear–NL zero-state 일치, + rigid translation/rotation·covariance와 + 기존 잠김/warped quadrature 한계를 검사한다. High-only 대조가 같은 물리장을 만들지 못하는 + 유의한 force/energy scale을 선택한다. +- Coupled consistent M에 low velocity/acceleration을 주어 off-diagonal inertia·T·momentum 및 + current moment를 검사한다. 비영 old \(r_n\)의 weighted 항등식·네 correction·세 update 식과 + full-state oracle를 대조한다. \(-\alpha r_n\)를 생략하는 대조는 실패해야 한다. +- Exact oracle로 \(X+u_h+u_l\)의 interval 포함을 확인하고 low가 현재 면 간격을 만드는 + 작은 형상 사례를 포함한다. Positive/invalid/uncertain 구분 및 low 누락의 false certificate를 + 검사한다. 불확실한 형상을 interval 대신 estimate로 승인하지 않는다. +- Accepted low가 다음 predictor/rates에 참여하는 copy/commit 사례, 실패 trial의 rollback, + rounded output/readback 및 geometry/backend/output 실패 보존을 검사한다. 기존 시간수렴, + 에너지·운동량·반력·외부 원본 기준도 유지한다. + +미해결 evidence는 구현한 primitive graph의 normal/subnormal 오차 bound와 지원 범위의 +독립 검토, paired coefficient 차이의 실제 측정, 위 불변량 및 원본 endpoint 검증이다. +추가 physics/acceptance 선택은 이 문서에서 만들지 않는다. 산술 bound나 작은 coupled 사례가 +충족되지 않으면 tolerance를 바꾸거나 더 높은 정밀도를 자동 채택하지 않고 근거를 반환한다. diff --git a/docs/hex8-structural-analysis/precision-proposal.md b/docs/hex8-structural-analysis/precision-proposal.md index 5bf706b..6c279b3 100644 --- a/docs/hex8-structural-analysis/precision-proposal.md +++ b/docs/hex8-structural-analysis/precision-proposal.md @@ -1,7 +1,8 @@ -# Hex8 동해석 내부 정밀도 변경 제안 — 미채택 +# Hex8 동해석 내부 정밀도 변경안 — 승인됨 -이 문서는 사용자 검토용 제안이다. 현재 [PRD](PRD.md)·[정식화](formulation.md)· -[numeric controls](numeric-controls.md)의 승인 계약을 변경하거나 구현 완료를 뜻하지 않는다. +2026-10-09 사용자가 A를 선택해 내부 high+low 물리 상태와 Float64 출력의 반올림 의미를 +승인했다. 이 문서는 [PRD](PRD.md)·[정식화](formulation.md)· +[numeric controls](numeric-controls.md)의 변경 근거이며, 구현 완료를 뜻하지 않는다. 실패 근거와 검증 범위는 [lessons](lessons.md)에 기록한다. ## 변경이 필요한 이유 @@ -13,7 +14,7 @@ 갱신 중 사라진 덧셈을 다음 반복에 누적하는 carry와, 실제 변위를 high+low로 표현하는 물리 상태는 다르다. Carry 실험의 kernel·rates·형상은 rounded high만 사용한다. -아래 제안은 low도 물리 계산에 사용하는 별도 정밀도 계약이며, 원본 수렴을 보장하지 않는다. +아래 승인안은 low도 물리 계산에 사용하는 별도 정밀도 계약이며, 원본 수렴을 보장하지 않는다. ## 제안하는 계산과 결과의 의미 @@ -37,8 +38,8 @@ 검증하고, 로그에 사용한 내부 산술을 명시한다. Restart는 계속 지원하지 않는다. Float64 직렬화의 overflow·비유한 값은 실패로 처리하며 기존 최종 파일을 보존한다. -현재 accepted state와 published Float64 state를 동일 계산 입력으로 취급하는 의미가 -달라지므로, 위 결과 의미와 내부 정밀도 변경을 함께 승인해야 한다. +승인한 변경에서는 accepted state가 high+low 전체 물리 상태이며, published Float64 state는 +그 반올림 observable이다. 같은 계산 입력이라는 이전 의미를 적용하지 않는다. ## 구현 전에 고정할 검증 diff --git a/docs/tracking/findings.md b/docs/tracking/findings.md index 93c0fdf..206cfc7 100644 --- a/docs/tracking/findings.md +++ b/docs/tracking/findings.md @@ -12,8 +12,9 @@ Float64 해의 부재를 증명하지 않으며, 수렴 기준·원본을 변경할 근거로 쓰지 않는다. - 후속 방향: [실패 근거](../hex8-structural-analysis/lessons.md)를 보존하고, 전체 물리 상태의 high+low 계산과 Float64 출력의 반올림 의미를 정하는 - [미채택 변경안](../hex8-structural-analysis/precision-proposal.md)을 사용자와 결정한다. - 변경을 승인하기 전에는 기존 precision·acceptance 계약을 유지한다. + [승인된 변경안](../hex8-structural-analysis/precision-proposal.md)을 구현한다. + 사용자가 2026-10-09 선택 A로 승인했다. 산술·요소의 첫 native 검사는 통과했으나 + 전체 동적 state·시간수렴·원본 endpoint 검증은 남아 있다. Acceptance 기준은 유지한다. ## HHT shell의 시간수렴 diff --git a/include/fesa/elements/hex8_element.h b/include/fesa/elements/hex8_element.h index a90476f..47f05c7 100644 --- a/include/fesa/elements/hex8_element.h +++ b/include/fesa/elements/hex8_element.h @@ -20,6 +20,16 @@ struct Hex8Response { std::array points{}; }; +/// @brief Normalized full-state response for Solid HHT high+low arithmetic. +/// @note Point rows are rounded observables; forces, Hessian and energy retain +/// their low components until global assembly and physical validation. +struct Hex8HighLowResponse { + Hex8Response high; + Vector internal_force_low{24U}; + Matrix tangent_low{24U, 24U}; + double strain_energy_low{0.0}; +}; + /// @brief Implements full-integration C3D8 linear and TL StVK mechanics. /// @note Reference geometry is owned and immutable; no trial is committed. class Hex8Element final : public Element { @@ -45,6 +55,15 @@ class Hex8Element final : public Element { Result EvaluateLinear(const Vector& displacement) const; /// @brief Returns one certified StVK trial and its analytic energy Hessian. Result EvaluateNonlinear(const Vector& displacement) const; + /// @brief Evaluates linear energy/force/Hessian on the complete paired state. + /// @note Uses the same reference-relative gradient operator as the paired TL + /// path. Inputs must have 24 finite normalized high+low translations. + Result EvaluateLinearHighLow( + const Vector& high, const Vector& low) const; + /// @brief Evaluates certified TL StVK mechanics without collapsing low state. + /// @note Failure preserves immutable reference data and all caller state. + Result EvaluateNonlinearHighLow( + const Vector& high, const Vector& low) const; /// @brief Returns the constant reference consistent mass; density required. Result ComputeMass() const; /// @brief Returns immutable positive reference Gauss volume weights. @@ -67,6 +86,9 @@ class Hex8Element final : public Element { /// @brief Evaluates without mutating the cached reference or previous values. Result Evaluate(const Vector& displacement, 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; /// @brief Attaches stable source context to a model failure. Status Failure(const char* message) const; diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index ad74a2e..bd85b6e 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -39,6 +39,7 @@ add_library( core/status.cpp elements/element_factory.cpp elements/hex8_element.cpp + elements/hex8_high_low.cpp elements/hex8_geometry.cpp elements/euler_beam_3d.cpp elements/finite_rotation.cpp diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index d930850..9c3d4ff 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -241,6 +241,14 @@ Result Hex8Element::EvaluateLinear(const Vector& u) const { Result Hex8Element::EvaluateNonlinear(const Vector& u) const { return Evaluate(u, true); } +Result Hex8Element::EvaluateLinearHighLow( + const Vector& high, const Vector& low) const { + return EvaluateHighLow(high, low, false); +} +Result Hex8Element::EvaluateNonlinearHighLow( + const Vector& high, const Vector& low) const { + return EvaluateHighLow(high, low, true); +} Result Hex8Element::Evaluate(const Vector& u, bool nonlinear) const { if (u.Size() != 24U) diff --git a/src/fesa/elements/hex8_geometry.cpp b/src/fesa/elements/hex8_geometry.cpp index e758daf..e7eff40 100644 --- a/src/fesa/elements/hex8_geometry.cpp +++ b/src/fesa/elements/hex8_geometry.cpp @@ -133,7 +133,7 @@ std::pair Split(const Coefficients& original, } // namespace Certificate Certify(const Hex8Element::Coordinates& coordinates, - const Vector* displacement) { + const Vector* displacement, const Vector* displacement_low) { std::array, 8> x{}; for (std::size_t a = 0; a < 8; ++a) for (std::size_t i = 0; i < 3; ++i) { @@ -145,6 +145,11 @@ Certificate Certify(const Hex8Element::Coordinates& coordinates, if (!std::isfinite(u)) return Certificate::kNonfinite; x[a][i] = Add(x[a][i], {u, u}); } + if (displacement_low) { + const auto low = (*displacement_low)[3U * a + i]; + if (!std::isfinite(low)) return Certificate::kNonfinite; + x[a][i] = Add(x[a][i], {low, low}); + } if (!Finite(x[a][i])) return Certificate::kNonfinite; } Coefficients coefficients{}; diff --git a/src/fesa/elements/hex8_geometry.h b/src/fesa/elements/hex8_geometry.h index 54a8a44..6e670c7 100644 --- a/src/fesa/elements/hex8_geometry.h +++ b/src/fesa/elements/hex8_geometry.h @@ -16,7 +16,8 @@ enum class Certificate { kPositive, kNonpositive, kUncertain, kNonfinite }; /// @brief Certifies the closed natural cube with interval Bernstein bounds. /// @note Current addition rounding is enclosed when displacement is supplied. Certificate Certify(const Hex8Element::Coordinates& coordinates, - const Vector* displacement = nullptr); + const Vector* displacement = nullptr, + const Vector* displacement_low = nullptr); /// @brief Evaluates trilinear shapes and natural derivatives in source order. void Shape(const Triple& natural, std::array& values, std::array& derivatives); diff --git a/src/fesa/elements/hex8_high_low.cpp b/src/fesa/elements/hex8_high_low.cpp new file mode 100644 index 0000000..efba073 --- /dev/null +++ b/src/fesa/elements/hex8_high_low.cpp @@ -0,0 +1,206 @@ +#include "fesa/elements/hex8_element.h" + +#include +#include +#include + +#include "../math/high_low_arithmetic.h" +#include "hex8_geometry.h" + +namespace fesa { +namespace { +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 math_internal::Subtract(a, b); +} +HighLow operator*(HighLow a, HighLow b) { + return math_internal::Multiply(a, b); +} +HighLow operator/(HighLow a, HighLow b) { + return math_internal::Divide(a, b); +} +HighLow& operator+=(HighLow& a, HighLow b) { return a = a + b; } + +/// @brief Contracts complete pairs in fixed component order. +HighLow Dot(const Triple& a, const Triple& b) { + HighLow result; + for (std::size_t i = 0; i < 3; ++i) result += a[i] * b[i]; + return result; +} +HighLow Determinant(const Tensor& a) { + return a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) - + a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) + + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]); +} +bool Finite(const Tensor& a) { + for (const auto& row : a) + for (auto value : row) + if (!math_internal::IsNormalized(value)) return false; + return true; +} +std::array Pack(const Tensor& a) { + using math_internal::Rounded; + return {Rounded(a[0][0]), Rounded(a[1][1]), Rounded(a[2][2]), + Rounded(a[0][1]), Rounded(a[0][2]), Rounded(a[1][2])}; +} +} // namespace + +Result Hex8Element::EvaluateHighLow( + const Vector& high, const Vector& low, bool nonlinear) const { + using math_internal::IsNormalized; + if (!math_internal::ArithmeticEnvironmentSupported()) + return Result::Failure(Failure( + "Hex8 high+low requires nearest rounding and gradual underflow.")); + if (high.Size() != 24U || low.Size() != 24U) + return Result::Failure( + Failure("Hex8 high+low trial requires two 24-component vectors.")); + std::array relative{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) { + const auto row = 3U * a + i; + const HighLow value{high[row], low[row]}; + const HighLow origin{high[i], low[i]}; + if (!IsNormalized(value) || !IsNormalized(origin)) + return Result::Failure( + Failure("Hex8 high+low translation is not finite normalized.")); + relative[row] = value - origin; + if (!IsNormalized(relative[row])) + return Result::Failure( + Failure("Hex8 high+low relative translation is unsupported.")); + } + const auto certificate = hex8_internal::Certify(coordinates_, &high, &low); + if (certificate != hex8_internal::Certificate::kPositive) { + const char* message = "Hex8 full high+low current geometry is nonpositive."; + if (certificate == hex8_internal::Certificate::kUncertain) + message = "Hex8 full high+low Jacobian positivity is unresolved at " + "depth 10 or cell budget 65536."; + if (certificate == hex8_internal::Certificate::kNonfinite) + message = "Hex8 full high+low Jacobian interval arithmetic is nonfinite."; + return Result::Failure(Failure(message)); + } + + Hex8HighLowResponse response; + std::array force{}; + std::array, 24> material{}, geometric{}; + HighLow energy; + for (std::size_t g = 0; g < 8; ++g) { + const auto& point = points_[g]; + std::array gradients{}; + Triple gradient_sum{}; + Tensor h{}; + for (std::size_t a = 1; a < 8; ++a) + for (std::size_t j = 0; j < 3; ++j) { + gradients[a][j] = HighLow{point.gradients[a][j]}; + gradient_sum[j] += gradients[a][j]; + for (std::size_t i = 0; i < 3; ++i) + h[i][j] += relative[3U * a + i] * gradients[a][j]; + } + // This is the derivative of the relative displacement expression above. + // Existing rounded G0 and cached K0 are not substituted into this operator. + for (std::size_t j = 0; j < 3; ++j) + gradients[0][j] = math_internal::Negate(gradient_sum[j]); + Tensor f{}; + Tensor strain{}, s{}, p{}, fft{}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) { + f[i][j] = nonlinear ? h[i][j] : HighLow{}; + if (i == j) f[i][j] += HighLow{1}; + strain[i][j] = HighLow{0.5} * (h[i][j] + h[j][i]); + if (nonlinear) + for (std::size_t k = 0; k < 3; ++k) + strain[i][j] += HighLow{0.5} * h[k][i] * h[k][j]; + } + const HighLow trace = strain[0][0] + strain[1][1] + strain[2][2]; + HighLow density_energy = HighLow{0.5} * HighLow{lambda_} * trace * trace; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) { + s[i][j] = HighLow{2} * HighLow{mu_} * strain[i][j]; + if (i == j) s[i][j] += HighLow{lambda_} * trace; + density_energy += HighLow{mu_} * strain[i][j] * strain[i][j]; + } + energy += density_energy * HighLow{point.weight}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) { + for (std::size_t k = 0; k < 3; ++k) { + p[i][j] += f[i][k] * s[k][j]; + fft[i][j] += f[i][k] * f[j][k]; + } + } + const HighLow det = Determinant(f); + if (!IsNormalized(det) || math_internal::Compare(det, HighLow{}) != 1) + return Result::Failure(Failure( + "Hex8 full deformation determinant is not finite positive.")); + Tensor stress = s; + if (nonlinear) { + stress = {}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) { + for (std::size_t k = 0; k < 3; ++k) + stress[i][j] += p[i][k] * f[j][k]; + stress[i][j] = stress[i][j] / det; + } + } + if (!Finite(h) || !Finite(f) || !Finite(strain) || !Finite(s) || + !Finite(p) || !Finite(fft) || !Finite(stress)) + return Result::Failure(Failure( + "Hex8 high+low tensor arithmetic is outside its supported range.")); + response.high.points[g] = {element_index_, static_cast(g + 1U), + point.natural, point.position, + Pack(stress), Pack(strain)}; + std::array fg{}; + for (std::size_t a = 0; a < 8; ++a) + for (std::size_t i = 0; i < 3; ++i) { + fg[a][i] = Dot(f[i], gradients[a]); + force[3U * a + i] += Dot(p[i], gradients[a]) * HighLow{point.weight}; + } + 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]); + HighLow geo; + 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]; + 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 = HighLow{lambda_} * fg[a][i] * fg[b][k] + + HighLow{mu_} * fg[b][i] * fg[a][k] + + HighLow{mu_} * fft[i][k] * dot; + material[row][column] += mat * HighLow{point.weight}; + if (i == k) + geometric[row][column] += geo * HighLow{point.weight}; + } + } + } + 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; + for (std::size_t i = 0; i < 24; ++i) { + if (!IsNormalized(force[i])) + return Result::Failure( + Failure("Hex8 high+low force is unsupported or nonfinite.")); + response.high.internal_force[i] = force[i].high; + response.internal_force_low[i] = force[i].low; + 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]) || + !IsNormalized(tangent)) + return Result::Failure( + Failure("Hex8 high+low tangent is unsupported or nonfinite.")); + response.high.material_tangent(i, j) = material[i][j].high; + response.high.geometric_tangent(i, j) = geometric[i][j].high; + response.high.tangent(i, j) = tangent.high; + response.tangent_low(i, j) = tangent.low; + } + } + return Result::Success(std::move(response)); +} +} // namespace fesa diff --git a/src/fesa/math/high_low_arithmetic.h b/src/fesa/math/high_low_arithmetic.h new file mode 100644 index 0000000..c724d1e --- /dev/null +++ b/src/fesa/math/high_low_arithmetic.h @@ -0,0 +1,162 @@ +#ifndef FESA_MATH_HIGH_LOW_ARITHMETIC_H_ +#define FESA_MATH_HIGH_LOW_ARITHMETIC_H_ + +#include +#include +#include +#include +#if defined(_M_X64) || defined(__SSE2__) +#include +#endif + +namespace fesa::math_internal { +/// @brief Checks the rounding and gradual-underflow premises without mutation. +inline bool ArithmeticEnvironmentSupported() { + if (std::fegetround() != FE_TONEAREST) return false; +#if defined(_M_X64) || defined(__SSE2__) + // MXCSR FTZ, DAZ and rounding-control bits must all be clear. + if ((_mm_getcsr() & 0xE040U) != 0U) return false; +#endif + return true; +} +/// @brief Two binary64 components of one physical value; raw parts are kept. +/// @note Arithmetic requires normalized finite parts and nearest rounding. +/// This is an expansion with bounded rounding, not an IEEE quad type. +struct HighLow { + double high{0}; + double low{0}; + HighLow() = default; + // Scalar lifting canonicalizes exact signed zero; raw pair inputs do not. + HighLow(double value) : high(value == 0 ? 0.0 : value) {} + HighLow(double high_part, double low_part) + : high(high_part), low(low_part) {} +}; + +/// @brief Gives an invalid arithmetic result for a checked caller to reject. +inline HighLow InvalidHighLow() { + return {std::numeric_limits::quiet_NaN(), 0}; +} +/// @brief Checks both physical components, before comparisons or publication. +inline bool IsFinite(HighLow value) { + return std::isfinite(value.high) && std::isfinite(value.low); +} +/// @brief Recovers the exact error of supported finite binary64 addition. +inline HighLow TwoSum(double left, double right) { + const double sum = left + right; + const double bridge = sum - left; + const double restored_left = sum - bridge; + const double error = (left - restored_left) + (right - bridge); + if (!std::isfinite(sum) || !std::isfinite(bridge) || + !std::isfinite(restored_left) || !std::isfinite(error)) + return InvalidHighLow(); + return sum == 0 && error == 0 ? HighLow{} : HighLow{sum, error}; +} +/// @brief Requires the raw pair to equal its nearest-rounding normalization. +inline bool IsNormalized(HighLow value) { + if (!IsFinite(value)) return false; + if (value.high == 0 && value.low == 0) + return !std::signbit(value.high) && !std::signbit(value.low); + const auto normalized = TwoSum(value.high, value.low); + return IsFinite(normalized) && normalized.high == value.high && + normalized.low == value.low; +} +/// @brief Rounds the full physical value to one finite output value. +inline double Rounded(HighLow value) { + return IsNormalized(value) ? value.high + value.low + : std::numeric_limits::quiet_NaN(); +} +/// @brief Adds normalized expansions; rounding of discarded tail is bounded. +inline HighLow Add(HighLow left, HighLow right) { + if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow(); + const auto head = TwoSum(left.high, right.high); + const auto tail = TwoSum(left.low, right.low); + const auto middle = TwoSum(head.low, tail.high); + const auto joined = TwoSum(head.high, middle.high); + const double remainder = (joined.low + middle.low) + tail.low; + if (!IsFinite(head) || !IsFinite(tail) || !IsFinite(middle) || + !IsFinite(joined) || !std::isfinite(remainder)) + return InvalidHighLow(); + return TwoSum(joined.high, remainder); +} +/// @brief Negates both components without collapsing the physical value. +inline HighLow Negate(HighLow value) { + if (!IsNormalized(value)) return InvalidHighLow(); + return value.high == 0 && value.low == 0 + ? HighLow{} + : HighLow{-value.high, -value.low}; +} +/// @brief Subtracts full expansions, including cancellations between highs. +inline HighLow Subtract(HighLow left, HighLow right) { + return Add(left, Negate(right)); +} +/// @brief Returns the exponent of the least nonzero binary64 significand bit. +/// @pre value is finite and nonzero; gradual underflow is enabled. +inline int LowestBitExponent(double value) { + int exponent = 0; + const double fraction = std::frexp(std::abs(value), &exponent); + auto significand = static_cast(std::ldexp(fraction, 53)); + int trailing = 0; + while ((significand & 1U) == 0U) { + significand >>= 1U; + ++trailing; + } + return exponent - 53 + trailing; +} +/// @brief Gives a product and its FMA residual when the dyadic tail fits. +/// @note Rejects underflow-risk exact tails, even if the rounded high is finite. +inline HighLow TwoProduct(double left, double right) { + if (!std::isfinite(left) || !std::isfinite(right)) return InvalidHighLow(); + if (left == 0 || right == 0) return {}; + if (LowestBitExponent(left) + LowestBitExponent(right) < -1074) + return InvalidHighLow(); + const double product = left * right; + if (!std::isfinite(product)) return InvalidHighLow(); + const double error = std::fma(left, right, -product); + if (!std::isfinite(error)) return InvalidHighLow(); + return TwoSum(product, error); +} +/// @brief Retains all four component products in deterministic order. +inline HighLow Multiply(HighLow left, HighLow right) { + if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow(); + auto product = TwoProduct(left.high, right.high); + product = Add(product, TwoProduct(left.high, right.low)); + product = Add(product, TwoProduct(left.low, right.high)); + return Add(product, TwoProduct(left.low, right.low)); +} +/// @brief Orders canonical values; 2 means unordered/unsupported input. +inline int Compare(HighLow left, HighLow right) { + if (!IsNormalized(left) || !IsNormalized(right)) return 2; + if (left.high < right.high) return -1; + if (left.high > right.high) return 1; + if (left.low < right.low) return -1; + if (left.low > right.low) return 1; + return 0; +} +/// @brief Takes the absolute value using the full physical sign. +inline HighLow Abs(HighLow value) { + const int sign = Compare(value, HighLow{}); + return sign == 2 ? InvalidHighLow() : (sign < 0 ? Negate(value) : value); +} +/// @brief Refines a binary64 quotient twice using the full residual. +/// @note Unsupported products/quotient underflow propagate an invalid result. +inline HighLow Divide(HighLow numerator, HighLow denominator) { + if (!IsNormalized(numerator) || !IsNormalized(denominator) || + Compare(denominator, HighLow{}) == 0) + return InvalidHighLow(); + if (Compare(numerator, HighLow{}) == 0) return {}; + const double first = numerator.high / denominator.high; + if (!std::isfinite(first) || first == 0) return InvalidHighLow(); + auto quotient = HighLow{first}; + for (int iteration = 0; iteration < 2; ++iteration) { + const auto residual = Subtract(numerator, Multiply(denominator, quotient)); + if (!IsNormalized(residual)) return InvalidHighLow(); + const double correction = residual.high / denominator.high; + if (!std::isfinite(correction) || (residual.high != 0 && correction == 0)) + return InvalidHighLow(); + quotient = Add(quotient, HighLow{correction}); + } + return quotient; +} +} // namespace fesa::math_internal + +#endif // FESA_MATH_HIGH_LOW_ARITHMETIC_H_ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d277846..3da1acd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -38,6 +38,7 @@ add_executable( unit/core/status_test.cpp unit/elements/element_factory_test.cpp unit/elements/hex8_element_test.cpp + unit/elements/hex8_high_low_test.cpp unit/elements/euler_beam_3d_test.cpp unit/elements/finite_rotation_test.cpp unit/elements/mitc4_shell_test.cpp @@ -48,6 +49,7 @@ add_executable( unit/elements/thermal_shell_test.cpp unit/fem/dof_manager_test.cpp unit/math/dense_blas_internal_test.cpp + unit/math/high_low_arithmetic_test.cpp unit/math/matrix_test.cpp unit/math/sparse_matrix_test.cpp unit/math/vector3_test.cpp diff --git a/tests/unit/elements/hex8_high_low_test.cpp b/tests/unit/elements/hex8_high_low_test.cpp new file mode 100644 index 0000000..9029006 --- /dev/null +++ b/tests/unit/elements/hex8_high_low_test.cpp @@ -0,0 +1,254 @@ +#include "fesa/elements/hex8_element.h" + +#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" + +namespace { +using fesa::math_internal::HighLow; +using fesa::math_internal::Subtract; + +fesa::Hex8Element UnitBrick(double offset = 0, double nu = 0) { + fesa::ElementDofLayout layout; + layout.node_indices = {0, 1, 2, 3, 4, 5, 6, 7}; + layout.components_per_node = {fesa::DofComponent::kUx, + fesa::DofComponent::kUy, + fesa::DofComponent::kUz}; + fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0}, + {1, 1, 0}, {0, 1, 0}, + {0, 0, 1}, {1, 0, 1}, + {1, 1, 1}, {0, 1, 1}}}; + for (auto& x : coordinates) x[0] += offset; + fesa::IsotropicLinearElasticMaterial material("test", 12, nu, {}); + return fesa::Hex8Element::Create(0, layout, coordinates, material, 2).Value(); +} + +TEST(Hex8HighLow, LowOnlyAffineStrainSurvivesUniformHighTranslation) { + auto brick = UnitBrick(); + fesa::Vector high(24), low(24); + const double strain = std::ldexp(1.0, -60); + for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1; + for (auto a : {1, 2, 5, 6}) low[3 * a] = strain; + const auto response = brick.EvaluateLinearHighLow(high, low); + ASSERT_TRUE(response.HasValue()); + for (std::size_t a = 0; a < 8; ++a) { + const bool positive = a == 1 || a == 2 || a == 5 || a == 6; + EXPECT_NEAR(response.Value().high.internal_force[3 * a], + (positive ? 3 : -3) * strain, 1e-13 * strain); + } + EXPECT_NEAR(response.Value().high.strain_energy, 6 * strain * strain, + 1e-13 * strain * strain); + for (const auto& row : response.Value().high.points) + EXPECT_NEAR(row.strain[0], strain, 1e-13 * strain); +} + +TEST(Hex8HighLow, NonlinearLowGradientProducesSecondOrderForce) { + auto brick = UnitBrick(); + fesa::Vector high(24), low(24); + const double strain = std::ldexp(1.0, -60); + for (std::size_t a = 0; a < 8; ++a) high[3 * a] = 1; + for (auto a : {1, 2, 5, 6}) low[3 * a] = strain; + const auto linear = brick.EvaluateLinearHighLow(high, low); + const auto nonlinear = brick.EvaluateNonlinearHighLow(high, low); + ASSERT_TRUE(linear.HasValue()); + ASSERT_TRUE(nonlinear.HasValue()); + const auto& l = linear.Value(); + const auto& n = nonlinear.Value(); + const auto increment = Subtract( + HighLow{n.high.internal_force[3], n.internal_force_low[3]}, + HighLow{l.high.internal_force[3], l.internal_force_low[3]}); + EXPECT_NEAR(increment.high, 4.5 * strain * strain, + 1e-12 * strain * strain); +} + +TEST(Hex8HighLow, RejectsInvalidLowShapeAndOverlappingComponents) { + auto brick = UnitBrick(); + fesa::Vector high(24), low(24); + EXPECT_FALSE(brick.EvaluateLinearHighLow(high, fesa::Vector(23)).HasValue()); + high[0] = 1; + low[0] = 1; + EXPECT_FALSE(brick.EvaluateNonlinearHighLow(high, low).HasValue()); +} + +TEST(Hex8HighLow, LinearAndNonlinearZeroStateHaveSameFullHessian) { + auto brick = UnitBrick(); + fesa::Vector zero(24); + const auto linear = brick.EvaluateLinearHighLow(zero, zero); + const auto nonlinear = brick.EvaluateNonlinearHighLow(zero, zero); + ASSERT_TRUE(linear.HasValue()); + ASSERT_TRUE(nonlinear.HasValue()); + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) { + EXPECT_EQ(linear.Value().high.tangent(i, j), + nonlinear.Value().high.tangent(i, j)); + EXPECT_EQ(linear.Value().tangent_low(i, j), + nonlinear.Value().tangent_low(i, j)); + } +} + +TEST(Hex8HighLow, CoupledBricksTransmitLowStrainAcrossSharedFace) { + const auto left = UnitBrick(), right = UnitBrick(1); + fesa::Vector high(24), left_low(24), right_low(24); + const double strain = std::ldexp(1.0, -60); + for (std::size_t a = 0; a < 8; ++a) { + high[3 * a] = 1; + const bool positive = a == 1 || a == 2 || a == 5 || a == 6; + left_low[3 * a] = positive ? strain : 0; + right_low[3 * a] = positive ? 2 * strain : strain; + } + const auto l = left.EvaluateLinearHighLow(high, left_low); + const auto r = right.EvaluateLinearHighLow(high, right_low); + ASSERT_TRUE(l.HasValue()); + ASSERT_TRUE(r.HasValue()); + const std::array face_left{1, 2, 5, 6}; + const std::array face_right{0, 3, 4, 7}; + for (std::size_t a = 0; a < 4; ++a) { + const auto i = 3 * face_left[a], j = 3 * face_right[a]; + const auto assembled = fesa::math_internal::Add( + {l.Value().high.internal_force[i], l.Value().internal_force_low[i]}, + {r.Value().high.internal_force[j], r.Value().internal_force_low[j]}); + EXPECT_NEAR(assembled.high, 0, 1e-13 * strain); + } + EXPECT_NEAR(l.Value().high.strain_energy + r.Value().high.strain_energy, + 12 * strain * strain, 1e-13 * strain * strain); +} + +TEST(Hex8HighLow, NodeZeroSignedDirectionMatchesForceAndEnergyDerivatives) { + using namespace fesa::math_internal; + auto brick = UnitBrick(); + fesa::Vector high(24), low(24); + for (std::size_t a = 0; a < 8; ++a) { + high[3 * a] = 1; + low[3 * a] = (a == 1 || a == 2 || a == 5 || a == 6) + ? std::ldexp(1.0, -60) : 0; + } + const auto baseline = brick.EvaluateNonlinearHighLow(high, low); + ASSERT_TRUE(baseline.HasValue()); + const auto& base = baseline.Value(); + for (int exponent : {-65, -70, -75}) { + const double step = std::ldexp(1.0, exponent); + auto plus_high = high, minus_high = high; + auto plus_low = low, minus_low = low; + const auto plus = Add({high[0], low[0]}, HighLow{step}); + const auto minus = Subtract({high[0], low[0]}, HighLow{step}); + plus_high[0] = plus.high; plus_low[0] = plus.low; + minus_high[0] = minus.high; minus_low[0] = minus.low; + const auto positive = brick.EvaluateNonlinearHighLow(plus_high, plus_low); + const auto negative = brick.EvaluateNonlinearHighLow(minus_high, minus_low); + ASSERT_TRUE(positive.HasValue()); + ASSERT_TRUE(negative.HasValue()); + for (std::size_t i = 0; i < 24; ++i) { + const auto derivative = Divide(Subtract( + {positive.Value().high.internal_force[i], positive.Value().internal_force_low[i]}, + {negative.Value().high.internal_force[i], negative.Value().internal_force_low[i]}), + HighLow{2 * step}); + const auto error = Subtract(derivative, + {base.high.tangent(i, 0), base.tangent_low(i, 0)}); + EXPECT_LE(std::abs(error.high), 1e-25); + } + const auto gradient = Divide(Subtract( + {positive.Value().high.strain_energy, positive.Value().strain_energy_low}, + {negative.Value().high.strain_energy, negative.Value().strain_energy_low}), + HighLow{2 * step}); + const auto energy_error = Subtract(gradient, + {base.high.internal_force[0], base.internal_force_low[0]}); + // Centered energy differences retain the StVK cubic term: O(step^2). + // For this unit E=12 brick its directional coefficient is bounded by 32. + EXPECT_LE(std::abs(energy_error.high), 32 * step * step + 1e-44); + } +} + +TEST(Hex8HighLow, GeometryCertificateIncludesLowOnlyCurrentCellWidth) { + fesa::Hex8Element::Coordinates coordinates{{{0, 0, 0}, {1, 0, 0}, + {1, 1, 0}, {0, 1, 0}, + {0, 0, 1}, {1, 0, 1}, + {1, 1, 1}, {0, 1, 1}}}; + fesa::Vector high(24), positive_low(24), negative_low(24); + for (std::size_t a = 0; a < 8; ++a) { + high[3 * a] = -coordinates[a][0]; + positive_low[3 * a] = std::ldexp(coordinates[a][0], -60); + negative_low[3 * a] = -positive_low[3 * a]; + } + using fesa::hex8_internal::Certificate; + EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &positive_low), + Certificate::kPositive); + EXPECT_EQ(fesa::hex8_internal::Certify(coordinates, &high, &negative_low), + Certificate::kNonpositive); + EXPECT_NE(fesa::hex8_internal::Certify(coordinates, &high), + Certificate::kPositive); +} + +TEST(Hex8HighLow, SignedZeroMaterialCoefficientHasIdenticalResponse) { + auto positive = UnitBrick(0, 0.0), negative = UnitBrick(0, -0.0); + fesa::Vector high(24), low(24); + high[3] = 0.01; + for (bool nonlinear : {false, true}) { + const auto p = nonlinear ? positive.EvaluateNonlinearHighLow(high, low) + : positive.EvaluateLinearHighLow(high, low); + const auto n = nonlinear ? negative.EvaluateNonlinearHighLow(high, low) + : negative.EvaluateLinearHighLow(high, low); + ASSERT_TRUE(p.HasValue()); + ASSERT_TRUE(n.HasValue()); + EXPECT_EQ(p.Value().high.strain_energy, n.Value().high.strain_energy); + EXPECT_EQ(p.Value().strain_energy_low, n.Value().strain_energy_low); + for (std::size_t i = 0; i < 24; ++i) { + EXPECT_EQ(p.Value().high.internal_force[i], n.Value().high.internal_force[i]); + EXPECT_EQ(p.Value().internal_force_low[i], n.Value().internal_force_low[i]); + for (std::size_t j = 0; j < 24; ++j) { + EXPECT_EQ(p.Value().high.tangent(i,j), n.Value().high.tangent(i,j)); + EXPECT_EQ(p.Value().tangent_low(i,j), n.Value().tangent_low(i,j)); + } + } + } +} + +TEST(Hex8HighLow, FiniteRotationRetainsTinyShearForceCovariance) { + using namespace fesa::math_internal; + auto brick = UnitBrick(); + fesa::Hex8Element::Coordinates x{{{0, 0, 0}, {1, 0, 0}, + {1, 1, 0}, {0, 1, 0}, + {0, 0, 1}, {1, 0, 1}, + {1, 1, 1}, {0, 1, 1}}}; + fesa::Vector zero(24), rotation(24), shear(24), rotated_shear(24); + const double gamma = std::ldexp(1.0, -60); + for (std::size_t a = 0; a < 8; ++a) { + // Q is an exact quarter turn; Q*(gamma*y,0,0)=(0,gamma*y,0). + rotation[3*a] = HighLow{-x[a][1] - x[a][0]}.high; + rotation[3*a+1] = x[a][0] - x[a][1]; + shear[3*a] = gamma*x[a][1]; + rotated_shear[3*a+1] = gamma*x[a][1]; + } + fesa::Vector turned_high(24), turned_low(24); + for (std::size_t j = 0; j < 24; ++j) { + const auto full = Add(HighLow{rotation[j]}, HighLow{rotated_shear[j]}); + turned_high[j] = full.high; turned_low[j] = full.low; + } + const auto origin = brick.EvaluateNonlinearHighLow(zero, zero); + const auto deformed = brick.EvaluateNonlinearHighLow(shear, zero); + const auto turned = brick.EvaluateNonlinearHighLow(rotation, zero); + const auto turned_deformed = brick.EvaluateNonlinearHighLow(turned_high, turned_low); + ASSERT_TRUE(origin.HasValue()); ASSERT_TRUE(deformed.HasValue()); + ASSERT_TRUE(turned.HasValue()); ASSERT_TRUE(turned_deformed.HasValue()); + // The fixed F64 reference gradients leave coefficient-level rotation error. + // Subtract the same rigid baseline; do not project that force to zero. + for (std::size_t a = 0; a < 8; ++a) { + std::array baseline{}, rotated{}; + for (std::size_t i = 0; i < 3; ++i) { + const auto j = 3*a+i; + baseline[i] = Subtract( + {deformed.Value().high.internal_force[j], deformed.Value().internal_force_low[j]}, + {origin.Value().high.internal_force[j], origin.Value().internal_force_low[j]}); + rotated[i] = Subtract( + {turned_deformed.Value().high.internal_force[j], turned_deformed.Value().internal_force_low[j]}, + {turned.Value().high.internal_force[j], turned.Value().internal_force_low[j]}); + } + EXPECT_LE(std::abs(Add(rotated[0], baseline[1]).high), 1e-13*gamma); + EXPECT_LE(std::abs(Subtract(rotated[1], baseline[0]).high), 1e-13*gamma); + EXPECT_LE(std::abs(Subtract(rotated[2], baseline[2]).high), 1e-13*gamma); + } +} +} // namespace diff --git a/tests/unit/math/high_low_arithmetic_test.cpp b/tests/unit/math/high_low_arithmetic_test.cpp new file mode 100644 index 0000000..092410d --- /dev/null +++ b/tests/unit/math/high_low_arithmetic_test.cpp @@ -0,0 +1,128 @@ +#include "../../../src/fesa/math/high_low_arithmetic.h" + +#include + +#include +#include + +namespace { +using namespace fesa::math_internal; + +TEST(HighLowArithmetic, RetainsSubUlpAdditionsAndSignedCancellation) { + ASSERT_TRUE(ArithmeticEnvironmentSupported()); + EXPECT_FALSE(IsNormalized(HighLow{-0.0, -0.0})); + EXPECT_TRUE(IsNormalized(HighLow{-0.0})); + EXPECT_FALSE(std::signbit(HighLow{-0.0}.high)); + EXPECT_FALSE(std::signbit(Negate(HighLow{}).high)); + EXPECT_FALSE(std::signbit(Negate(HighLow{}).low)); + const double tiny = std::ldexp(1.0, -60); + const auto positive = Add(HighLow{1}, HighLow{tiny}); + EXPECT_EQ(positive.high, 1); + EXPECT_EQ(positive.low, tiny); + EXPECT_EQ(Rounded(positive), 1); + const auto difference = Subtract(positive, HighLow{1}); + EXPECT_EQ(difference.high, tiny); + EXPECT_EQ(difference.low, 0); + EXPECT_EQ(Compare(positive, HighLow{1}), 1); + EXPECT_EQ(Compare(HighLow{1, -tiny}, HighLow{1}), -1); + EXPECT_EQ(Compare(Add(positive, Negate(positive)), HighLow{}), 0); + EXPECT_EQ(Abs(HighLow{-1, -tiny}).low, tiny); +} + +TEST(HighLowArithmetic, AllProductTermsAndDivisionResidualReachState) { + const double tiny = std::ldexp(1.0, -60); + const auto product = Multiply(HighLow{1, tiny}, HighLow{1, -tiny}); + EXPECT_EQ(product.high, 1); + EXPECT_EQ(product.low, -std::ldexp(1.0, -120)); + const auto quotient = Divide(HighLow{1}, HighLow{3}); + ASSERT_TRUE(IsNormalized(quotient)); + const auto defect = Subtract(Multiply(quotient, HighLow{3}), HighLow{1}); + EXPECT_LE(std::abs(defect.high), std::ldexp(1.0, -104)); + const auto cancellation = Subtract(HighLow{1, tiny}, HighLow{1, tiny / 2}); + EXPECT_EQ(cancellation.high, tiny / 2); + EXPECT_EQ(cancellation.low, 0); +} + +TEST(HighLowArithmetic, NormalizesTiesAndCrossesPowerOfTwoBoundary) { + const double half_ulp = std::ldexp(1.0, -53); + EXPECT_TRUE(IsNormalized(HighLow{1, half_ulp})); + EXPECT_FALSE(IsNormalized(HighLow{std::nextafter(1.0, 2.0), -half_ulp})); + auto result = HighLow{1}; + for (int i = 0; i < 4; ++i) result = Add(result, HighLow{half_ulp / 2}); + EXPECT_EQ(result.high, std::nextafter(1.0, 2.0)); + EXPECT_EQ(result.low, 0); + const auto minimum = HighLow{std::numeric_limits::denorm_min()}; + EXPECT_EQ(Add(minimum, minimum).high, 2 * minimum.high); + EXPECT_EQ(Multiply(minimum, HighLow{1}).high, minimum.high); +} + +TEST(HighLowArithmetic, RejectsUnsupportedDyadicTailAndNonfiniteResults) { + const auto minimum = HighLow{std::numeric_limits::denorm_min()}; + EXPECT_FALSE(IsFinite(Multiply(minimum, HighLow{0.5}))); + EXPECT_FALSE(IsFinite(Multiply(HighLow{std::ldexp(1.0, -1022)}, + HighLow{std::nextafter(0.5, 1.0)}))); + EXPECT_FALSE(IsFinite(Divide(minimum, HighLow{2}))); + EXPECT_FALSE(IsFinite(Divide(HighLow{1}, HighLow{}))); + EXPECT_FALSE(IsFinite(Add(HighLow{std::numeric_limits::max()}, + HighLow{std::numeric_limits::max()}))); + EXPECT_EQ(Compare(InvalidHighLow(), HighLow{}), 2); + EXPECT_FALSE(IsNormalized(HighLow{1, 1})); + EXPECT_TRUE(std::isnan(Rounded(HighLow{1, 1}))); +} + +TEST(HighLowArithmetic, ExactDyadicOracleBoundsDiscardedTailAndRefinement) { + // Expected pairs were obtained independently using exact rational arithmetic. + // Add and multiply deliberately discard a term at 2^-108; neither is exact. + for (int scale : {-100, 0, 100}) { + const auto x = HighLow{std::ldexp(1.0, scale), std::ldexp(1.0, scale-54)}; + const auto y = HighLow{std::ldexp(1.0, scale-53), std::ldexp(1.0, scale-108)}; + const auto sum = Add(x, y); + EXPECT_EQ(sum.high, std::ldexp(0x1.0000000000001p+0, scale)); + EXPECT_EQ(sum.low, -std::ldexp(1.0, scale-54)); + const auto square = Multiply(x, x); + EXPECT_EQ(square.high, std::ldexp(1.0, 2*scale)); + EXPECT_EQ(square.low, std::ldexp(1.0, 2*scale-53)); + } + const auto third = Divide(HighLow{1}, HighLow{3}); + EXPECT_EQ(third.high, 0x1.5555555555555p-2); + EXPECT_EQ(third.low, 0x1.5555555555555p-56); + const auto third_residual = Subtract(Multiply(HighLow{3}, third), HighLow{1}); + EXPECT_EQ(third_residual.high, -std::ldexp(1.0, -108)); + EXPECT_EQ(third_residual.low, 0); + const auto refined = Divide(HighLow{1, std::ldexp(1.0,-60)}, + HighLow{1.25, -std::ldexp(1.0,-60)}); + EXPECT_EQ(refined.high, 0x1.999999999999ap-1); + EXPECT_EQ(refined.low, -0x1.8e147ae147ae1p-55); + const auto boundary = TwoProduct(0x1.0000000000001p+0, + 0x1.0000000000001p-970); + EXPECT_EQ(boundary.high, 0x1.0000000000002p-970); + EXPECT_EQ(boundary.low, std::numeric_limits::denorm_min()); + EXPECT_FALSE(IsFinite(TwoProduct(std::nextafter(1.0, 2.0), + std::nextafter(std::numeric_limits::min(), 1.0)))); +} + +TEST(HighLowArithmetic, RejectsChangedRoundingAndFlushToZeroEnvironment) { + struct Restore { + int rounding = std::fegetround(); +#if defined(_M_X64) || defined(__SSE2__) + unsigned csr = _mm_getcsr(); +#endif + ~Restore() { + std::fesetround(rounding); +#if defined(_M_X64) || defined(__SSE2__) + _mm_setcsr(csr); +#endif + } + } restore; + ASSERT_TRUE(ArithmeticEnvironmentSupported()); + ASSERT_EQ(std::fesetround(FE_DOWNWARD), 0); + EXPECT_FALSE(ArithmeticEnvironmentSupported()); + ASSERT_EQ(std::fesetround(FE_TONEAREST), 0); +#if defined(_M_X64) || defined(__SSE2__) + _mm_setcsr(restore.csr | 0x8000U); + EXPECT_FALSE(ArithmeticEnvironmentSupported()); + _mm_setcsr(restore.csr | 0x0040U); + EXPECT_FALSE(ArithmeticEnvironmentSupported()); +#endif +} +} // namespace