From d7f6d0ce746018c30dbf0226863a34a1c844d823 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 18 Oct 2026 23:42:54 +0900 Subject: [PATCH] chore: preserve approved MITC4 tolerance experiment baseline --- docs/README.md | 1 + docs/abaqus-tolerance/README.md | 192 ++++++++++++++++++ docs/abaqus-tolerance/application-plan.md | 139 +++++++++++++ .../local-kernel-experiment.md | 132 ++++++++++++ docs/hht-implicit-dynamics/README.md | 8 + .../nonlinear-static-newton-raphson/README.md | 6 + docs/tracking/findings.md | 15 ++ docs/tracking/status.md | 7 + src/fesa/elements/nonlinear_mitc4_inertia.cpp | 15 +- src/fesa/elements/nonlinear_mitc4_shell.cpp | 5 +- .../unit/analysis/hht_shell_physics_test.cpp | 61 +++--- .../elements/nonlinear_mitc4_inertia_test.cpp | 110 +++++----- .../nonlinear_mitc4_shell_tangent_test.cpp | 49 ++--- 13 files changed, 630 insertions(+), 110 deletions(-) create mode 100644 docs/abaqus-tolerance/README.md create mode 100644 docs/abaqus-tolerance/application-plan.md create mode 100644 docs/abaqus-tolerance/local-kernel-experiment.md diff --git a/docs/README.md b/docs/README.md index 6813fe4..ef925cc 100644 --- a/docs/README.md +++ b/docs/README.md @@ -38,4 +38,5 @@ Formulation의 수식·가정은 원문으로 보존한다. 문서 안의 ready/ - [Hex8 구조해석](hex8-structural-analysis/README.md) - [Shell 열전달](shell-transient-heat-transfer/README.md) - [계산 성능](calculation-performance/README.md) +- [Abaqus 수렴 허용치 조사와 FESA 적용 계획](abaqus-tolerance/README.md) - [C++ 모듈화·Refactor 병합 현황과 이관 근거](cpp-object-oriented-modular-refactoring/README.md) diff --git a/docs/abaqus-tolerance/README.md b/docs/abaqus-tolerance/README.md new file mode 100644 index 0000000..5779df8 --- /dev/null +++ b/docs/abaqus-tolerance/README.md @@ -0,0 +1,192 @@ +# Abaqus 수렴 허용치 조사와 FESA 대응 + +조사일: 2026-10-18. FESA 조사 기준: `43867d0`. +상태: 문헌·코드 조사이며, 허용치 변경이나 성능 실험 결과가 아니다. + +## 조사 결론 + +Abaqus/Standard의 일반 비선형 평형 기준은 잔차 비율 **0.005**, 해 보정 비율 +**0.01**이다. 병진과 회전을 각각 검사한다. FESA의 현재 MITC4 HHT는 힘·모멘트 +잔차와 q/v/a/A 보정량에 **1e-10**을 사용한다. 완화 실험을 검토할 근거는 있으나, +양쪽의 정규화와 검사 대상이 달라 숫자의 비로 정확도나 실행시간을 예측할 수 없다. +일반 기준값의 근거는 [Abaqus 공식 해석 지침의 제어값 설명][S1]이다. + +비선형 반복 종료, 선형방정식 풀이 정확도, 시간적분 정확도, 결과 비교 오차는 서로 +다른 기준이다. `1e-10`을 전역 검색해 `0.005`로 치환하는 접근은 부적절하다. +이 문서에서의 적용 권고는 **실험 제안**이며 현재 승인된 FESA 판정을 대체하지 않는다. + +## 출처와 버전의 한계 + +- 아래 Abaqus 수치는 Dassault Systèmes가 작성한 **Abaqus 2025 공식 문서**의 공개 + 호스팅본에서 확인했다. 각 링크는 사용한 문서의 버전을 포함한다. +- 저장소의 [동해석 기준 입력](../../reference/shellNL-dyn-meshfine/shellNL-dyn-meshfine.inp) + 3행은 `Abaqus/CAE Learning Edition 2024` 생성 표기다. 이는 실제 solver 실행 버전이나 + 적용된 모든 제어값을 증명하지 않는다. +- 2024판의 대응 문서 주소는 조사 환경에서 HTTP 403으로 열리지 않았다. 2025값을 + 2024에서 확인한 값으로 보고하지 않는다. +- 해당 기준 폴더에는 INP와 변위 CSV만 있고 `.msg`/`.dat`는 없다. 공식 지침에 따르면 + 실제 적용 제어값은 이 두 로그에서 확인할 수 있다.[S1] +- 따라서 현재 확정 가능한 것은 **2025 문서 기본값과 FESA 코드의 차이**다. + 기준 해석과 완전히 같은 설정이라는 주장은 실제 실행 로그 또는 해당 버전 원문을 + 확보한 뒤에만 가능하다. 이번 조사에서는 Abaqus를 실행하지 않았다. + +## 1. 비선형 field 방정식의 기본값 + +다음은 `*CONTROLS, PARAMETERS=FIELD` 기본값이다. α는 field 종류다. +q̃0/q̃u 외에는 무차원 비율이며 절대 힘·길이로 해석하지 않는다.[S2] + +| 기호 | 기본값 | 역할 | +| --- | ---: | --- | +| Rn | 5e-3 | 최대 scaled 잔차 / 평균 flux | +| Cn | 1e-2 | 최대 Newton 보정 / 해당 증분의 최대 변화 | +| q̃0 | 첫 step 1e-2; 이후 이전 step 평균 | 평균 flux 초기값, 해당 flux 단위 | +| q̃u | 미지정 | 사용자가 지정하는 고정 평균 flux, 해당 flux 단위 | +| Rp | 2e-2 | 비이차 수렴에서 사용하는 대체 잔차 비율 | +| ε | 1e-5 | zero-flux 분류 기준 | +| Cε | 1e-3 | zero-flux 시 보정 비율 | +| Rl | 1e-8 | 선형 증분 판정용 잔차 비율 | +| Cf | 1 | 미소 field의 물리 차원 변환 비율 | +| εl | 1e-5 | 평균 계산의 비활성 flux 구분 | +| εd | 1e-8 | 특성길이 대비 미소 병진 증분 구분 | + +첫 데이터 줄 순서는 `Rn, Cn, q̃0, q̃u, Rp, ε, Cε, Rl`, +둘째 줄은 `Cf, εl, εd`다.[S2] + +2025 문서의 `CONVERGENCE CHECK` 기본 모드는 `STRICT`다. `MODERATE`, `RELAXED`, +`CONCEPT DESIGN`은 분모나 잔차/보정의 결합이 달라지는 별도 선택이며, 이번 기본값 +비교에서 이 모드들을 섞지 않는다.[S2] + +## 2. 같은 숫자가 같은 판정이 아닌 이유 + +Abaqus의 보통 잔차 검사는 `max_i |r_i/(f_i q̃)| <= Rn`, 보정 검사는 +`cmax <= Cn Δumax`다. 여기서 q̃는 요소가 절점에 주는 flux와 외부 flux의 크기로 +계산한 공간 평균을 다시 증분 이력에서 평균한 값이다. 단순한 최대 외력이나 +조립 후 잔차 평균이 아니다. 종속 절점이 없는 경우 fi=1이다.[S3] + +기본값 재현에는 다음 분기가 함께 필요하다.[S3] + +| 상황 | 의미 | +| --- | --- | +| 보통 비선형 | 각 활성 field에서 잔차와 실제 보정 또는 추정 다음 보정 검사 | +| zero flux | 공간 평균이 εq̃ 이하일 때 별도 잔차/보정 분기 | +| 선형 증분 | 모든 field의 잔차가 Rl 이하면 보정 검사 없이 수렴 가능 | +| 비이차 수렴 | IP 이후, 해당 조건에서 Rp와 보정 검사를 함께 적용; quasi-Newton 제외 | +| 비활성 영역 | 작은 flux를 평균에서 제외하는 규칙으로 평균 희석 방지 | +| 미소 field/증분 | 힘↔모멘트 특성길이 변환과 미소 병진 증분 예외 | + +이 분기들을 생략하면 “Abaqus 기본 수치를 사용한 부분 정책”으로 표현해야 한다. +특히 zero-flux는 **외력이 정확히 0인 경우**와 같지 않다. 평균 이력은 실패한 trial이나 +중간 Newton 반복마다 누적해서는 안 된다. FESA 도입 시 동일 trial의 평균 후보를 계산하고, +accepted increment에서만 이력을 반영하는 설계가 필요하다. + +## 3. MITC4 동해석에 관계되는 시간적분 기준 + +Abaqus의 half-increment residual은 증분 중간의 평형 불균형으로 시간간격의 정확도를 +판단하는 양이다. 증분 끝의 Newton 잔차 허용치가 아니다. 공식 `*DYNAMIC` 문서의 구분은 +다음과 같다.[S4] + +| 항목 | 문서 기본값 또는 지침 | FESA와의 관계 | +| --- | --- | --- | +| HALFINC SCALE FACTOR | transient fidelity: 비접촉 1000, 접촉 10000 | 시간 평균 힘·모멘트에 곱하는 비율 | +| HAFTOL | 사용자가 지정하는 힘 차원의 값 | DIRECT와 함께 지정 불가 | +| HAFTOL 선택 지침 | 저감쇠 탄성: 보통 정확도 대표 힘의 1–10배, 높은 정확도 0.1–1배 | 기본값이 아닌 선택 지침 | +| HAFTOL 선택 지침 | 소성 등 소산이 큰 경우: 보통 10–100배, 높은 정확도 1–10배 | 현재 FESA 탄성 shell의 직접 대응 아님 | +| NOHAF | half-increment residual 계산 생략 | DIRECT에서도 계산은 기본 수행되며 생략 가능 | +| DIRECT | 고정 시간간격 | HAFTOL·HALFINC SCALE FACTOR와 동시 지정 불가 | + +Moderate dissipation과 quasi-static의 half-increment 검사는 기본적으로 억제된다. +1000/10000은 매우 큰 Newton 허용치가 아니라 별도 시간 정확도 지표의 배율이다.[S5] + +저장소 기준 입력 1475–1477행은 `*Dynamic,direct`와 `0.1,10.,`다. +따라서 100개 고정 증분이며, 자동 시간간격 조절 배율을 바꾸어 이 모델의 증분 수를 +줄이는 계획은 현재 입력 계약과 다르다. `INC=1000`은 실제 증분 수가 아닌 상한이다. + +Abaqus에는 수렴 실패 상태도 받아들이는 `DIRECT=NO STOP`이 있지만, 고정 간격을 유지하는 +`DIRECT`와 다르다. FESA의 실패 candidate 보존 계약에 반하므로 제안하지 않는다.[S6] + +## 4. 반복·시간간격 제어값 + +아래는 허용오차 자체보다 재시도·반복 횟수에 관계되는 보조 제어값이다. +자동 시간간격과 고정 시간간격에서 적용 방식이 다르다. 이를 FESA의 100회 상한과 +그대로 교환하면 같은 해석 절차가 되지 않는다.[S1] + +| 제어값 | 기본값 | 의미 | +| --- | ---: | --- | +| I0 / IR | 4 / 8 | 발산 검사 / 수렴률 검사 시작 | +| IP | 9 | 대체 잔차 기준을 검토하는 반복 수 | +| IC | 16 | 연속 평형 반복 제어 상한; 세부 조건은 원문 참조 | +| IA | 5 | 증분 시도 횟수 | +| Df / DC | 0.25 / 0.5 | 발산 / 느린 수렴 시 자동 cutback | +| DB / DD | 0.75 / 1.5 | 다음 증분 축소 / 빠른 수렴 시 증가 | +| 동해석 DM | 1.25 | 자동 증분 증가의 최대 배율.[S2] | + +이번 FESA 실험 제안에서는 먼저 기존 고정 시간간격·반복 상한을 유지하여, +종료 판정 변경의 효과를 분리하는 편이 타당하다. 이는 조사에 따른 권고다. + +## 5. 다른 종류의 tolerance + +| 종류 | 확인값과 범위 | FESA 적용 판단 | +| --- | --- | --- | +| 반복 선형방정식 solver | 2025 기본 상대잔차 1e-6, 최대 500회.[S7] | PARDISO direct solve의 backward error와 교환 불가 | +| 선형 잔차의 정의 | 반복 solver는 L2 잔차 `||Ku-f||/||f||`; Newton 종료와 독립.[S8] | norm·분모가 다르므로 숫자만 비교 불가 | +| 구속 적합성 | hybrid 체적/축/전단 1e-5; distributing coupling 변위/회전 1e-5.[S2] | FESA hidden mass-null constraint와 다른 방정식 | +| 접촉 | 적합성 5e-3, 저압 soft-contact 0.1, CONVERT SDI=YES force error 1.[S2] | 현재 지원 subset 밖; MITC4 shear/drill 기준으로 사용 금지 | +| 열전달 | field 기준 외에 DELTMX는 사용자가 정하는 최대 온도 변화.[S9] | 현재 선형 열전달 잔차 인증의 완화 근거가 아님 | +| 고유치 | Lanczos/AMS/Subspace별 절차·설정 구분.[S10] | FEAST의 잔차·직교성 1e-10 대응값으로 임의 채택 금지 | +| 외부 결과 비교 | Abaqus Newton 기준과 다른 목적 | 기존 FESA U/RF/S 등의 비교 허용치는 유지 | + +조사 범위는 FESA 구조해석, 특히 MITC4 implicit dynamics에 관련되는 기준이다. +접촉·hybrid는 경계 설명용이며, Abaqus 전체 제품의 재료별 적분·최적화·손상·확산· +Explicit 안정시간 기준을 모두 열거한 문서는 아니다. 공개 문헌에서 확인하지 않은 +Abaqus 내부 직접분해·pivot·rank 판정 숫자는 추측하지 않았다. + +## 6. 현재 FESA 기준과 소유 위치 + +| 대상 | 현재 판정 | 근거/실험 시 취급 | +| --- | --- | --- | +| MITC4 비선형 정적 | 잔차·보정 상대 1e-2; 외력 또는 증분 scale=0에서 typed absolute 1e-3 | [공개 기준](../../include/fesa/analysis/nonlinear_static_control.h), [실제 분기](../../src/fesa/analysis/nonlinear_static_control.cpp); 동해석과 별도 | +| Shell HHT 외부 Newton 잔차 | 병진/회전 별 `max|r|/max(scale) <=1e-10` | [조립과 scale](../../src/fesa/analysis/hht_dynamic_assembly.cpp), [ResidualRecord](../../src/fesa/analysis/hht_reduced_procedure.cpp) | +| Shell HHT 외부 Newton 보정 | 병진은 L, 회전은 1로 나눈 δq·hδv·h²δa·h²δA 모두 1e-10 | 같은 파일 `CorrectionRecord`; Abaqus의 증분 대비 보정과 다름 | +| Shell HHT 내부 위치 root | hidden residual 1e-10, 회전 보정 1e-10 등; 최대 100 correction | [SolveEmbeddingRoot](../../src/fesa/analysis/hht_reduced_operators.cpp); 별도 구속 정밀도 | +| 회전만 추가 보정하는 solve 경로 | 병진 통과/회전 미통과와 병진 보정량의 1e-10에 의존 | 같은 파일 `SolveRotationOnlyCorrection`; 종료 정책과 결합됨 | +| 초기 frame0 일관성 | 초기·재구성 full 잔차 1e-10 | [procedure](../../src/fesa/analysis/hht_reduced_procedure.cpp); 양의 증분 종료와 분리 | +| HHT 선형 solve | scaled backward error 1e-12 | [CheckLinearSolve](../../src/fesa/analysis/hht_dynamic_assembly.cpp); 유지 대상 | +| HHT 전역 힘·모멘트 평형 | 각각 normalized 1e-8 | 같은 파일 `CompleteEquilibrium`; 외부 Newton 완화 후에도 남는 필수 acceptance | +| 이력 update·basis·pair 산술 | same-unit update 1e-12 및 exact-zero/정확 반올림 검사 | [procedure](../../src/fesa/analysis/hht_reduced_procedure.cpp), [pair 검사](../../src/fesa/analysis/hht_reduced_increment.h); 유지 대상 | +| hidden 위치/rate·rank·null | 별도 1e-10/1e-12 등의 물리·대수 검사 | [algebra](../../src/fesa/analysis/hht_dynamic_algebra.cpp), [기능 요약](../hht-implicit-dynamics/README.md); 일괄 완화 제외 | +| 선형 정적·고유치·열 | 각 기능별 독립 기준 | 이번 MITC4 동해석 실험의 자동 변경 대상 아님 | + +HHT의 scale은 외력과 요소별 physical/drilling force, convective inertia 및 질량×가속도 +항들의 **상쇄 전 크기**를 누적한다. Abaqus의 q̃와 동일하지 않다. 현재 보정 기준도 +`δq/Δq`가 아니라 특성길이·시간간격에 기반하므로 `1e-10 → 0.01`은 단순 완화 이상의 +의미 변경이다. 정적은 이미 상대 1e-2이므로 Abaqus의 Rn=5e-3 채택이 숫자상 더 엄격할 +수도 있다. “FESA 전체가 Abaqus보다 엄격하다”로 일반화할 수 없다. + +## 7. 성능 가설과 적용 전 확인 + +확인된 사실은 엄격한 HHT 종료 기준과, 외부 Newton 외에도 비용이 큰 내부 root· +요소 도함수·조립 경로가 존재한다는 점이다. 허용치가 실제 장시간 실행의 주원인인지는 +현재 측정하지 않았다. [기존 성능 기록](../calculation-performance/README.md)의 +6372.325621초는 과거 source의 단일 baseline으로 최신 revision의 측정값이 아니다. + +다음 실험에서 외부/내부 반복 수, 요소 평가 수, factor/solve 수, 단계별 시간과 +종료를 지연시킨 채널을 함께 측정해야 한다. 내부 구속 solve가 대부분의 시간을 쓰면 +외부 Newton 기준만 바꿔도 효과가 작을 수 있다. 기준을 넓혀 성공한 것과 해의 정확도가 +유지된 것을 구분하고, U 비교와 별도로 V/physical a·구속·시간수렴을 확인해야 한다. + +적용 단계와 승인 경계는 [적용 계획](application-plan.md)에 정리한다. + +## 출처 + +아래는 모두 Abaqus 2025 공식 문서의 공개 호스팅본이다. 실제 2024 실행값 확인은 남아 있다. + +[S1]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEANLRefMap/simaanl-c-convergecontrol.htm "Commonly Used Control Parameters" +[S2]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEKEYRefMap/simakey-r-controls.htm "*CONTROLS" +[S3]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEANLRefMap/simaanl-c-convergcriteria.htm "Convergence Criteria for Nonlinear Problems" +[S4]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEKEYRefMap/simakey-r-dynamic.htm "*DYNAMIC" +[S5]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAECAERefMap/simacae-t-simconfiguredynamicimplicit.htm "Configuring a dynamic, implicit procedure" +[S6]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEANLRefMap/simaanl-c-dynamic.htm "Implicit Dynamic Analysis Using Direct Integration" +[S7]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEKEYRefMap/simakey-r-solvercontrols.htm "*SOLVER CONTROLS" +[S8]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEANLRefMap/simaanl-c-itrsolveroverview.htm "Iterative Linear Equation Solver" +[S9]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEKEYRefMap/simakey-r-heattransfer.htm "*HEAT TRANSFER" +[S10]: https://docs.software.vt.edu/abaqusv2025/English/SIMACAEKEYRefMap/simakey-r-frequency.htm "*FREQUENCY" diff --git a/docs/abaqus-tolerance/application-plan.md b/docs/abaqus-tolerance/application-plan.md new file mode 100644 index 0000000..3898b46 --- /dev/null +++ b/docs/abaqus-tolerance/application-plan.md @@ -0,0 +1,139 @@ +# MITC4 HHT에 Abaqus 수렴 규칙을 적용하는 실행 계획 + +작성일: 2026-10-18. 조사 기준 source: 43867d0. +상태: 실제 구현·측정·검증을 위한 개정안. 이번 ready에서는 계획 문서만 수정한다. +수정안 승인 후 go는 **실제 코드 수정·빌드·FESA 실행·성능 측정·정확도 검증**을 수행한다. +이 계획의 Abaqus 외부 Newton 규칙 구현·측정은 아직 실행하지 않았다. + +2026-10-18 후속 사용자 지시로 별도 [요소 local 근사 실험](local-kernel-experiment.md)을 +먼저 적용했다. 아래 요소 수식·질량 보존 조건의 명시적 예외는 polar 1회, 회전/Jacobian +3차 급수, 질량 검사 1e-3이다. 후속 지시로 관련 관성·접선·HHT 물리 테스트 오차도 +1e-3으로 완화했지만 기존 시간수렴 차수·실제 해석 controller는 유지했다. +이는 Abaqus 기본값이 아니며 정적 비선형 MITC4에도 +polar 변경이 공유된다. 외부 수렴 규칙의 효과를 분리하려면 strict/candidate 양쪽에 같은 +local kernel을 사용하고 원래 43867d0 기준과의 차이도 별도로 기록해야 한다. + +## 목표와 확정한 방식 + +MITC4 implicit dynamics의 외부 Newton 종료에 **Abaqus/Standard 2025의 공개 STRICT field +수렴 규칙**을 적용한다. 숫자만 바꾸거나 기존 strict 종료를 계속 사용하는 방식으로 축소하지 않는다. +근거는 [조사 보고서](README.md)다. 2024 생성 INP의 실제 solver 설정 미확인은 별도 한계로 +기록하며, 그 로그가 없다는 이유로 2025 기준 후보 구현을 보류하지 않는다. + +공개 규칙을 FESA의 기존 물리량에 대응하는 실험이다. 잔차는 기존 full free endpoint +E=I+fint-F, field는 병진 u와 total rotation-vector φ, 회전 flux는 generalized moment다. +요소의 total internal+physical inertia flux와 외부 source contribution을 평균에 사용하고 +상쇄 전 개별 항 합인 기존 residual_scale과 혼동하지 않는다. 미소 field 변환의 특성길이는 +기존 ReferenceLength의 초기 shell 최대 edge 길이에 대응한다. +Abaqus의 비공개 요소·관성·회전 처리와 내부 특성길이까지 같다고 보장하지 않는다. + +대상은 MITC4 HHT 양의 증분이다. 정적·Hex8·열·고유치의 tolerance, 요소 수식·질량·HHT 적분식, +시간간격 dt=0.1/T=10/100증분 및 기존 최대100 correction은 유지한다. + +## 변경하는 판정과 유지하는 검사 + +| 변경하는 부분 | 유지하는 부분 | +| --- | --- | +| 외부 Newton 잔차·보정의 종료 판정 | frame0 초기 일관성 | +| 힘·모멘트 평균과 accepted 증분 이력 | 내부 graph/hidden constraint·rank/null | +| 일반·zero·linear·추정 보정·비이차 수렴 분기 | 선형 backward error1e-12·update1e-12·pair/rounding | +| q/φ의 증분 대비 보정으로 종료 | PrivateCorrectionEqual·physical rate 복원 | +| 기존 q/v/a/A1e-10 크기는 strict 비교 진단 | 전역 힘·모멘트 평형 각각1e-8·finite·영 구속 | +| 의미가 달라지는 후보 로그 hht-log-v2 | CLI·지원 INP·HDF5 schema2·최종 결과 보호 | + +Rn=.005, Cn=.01, Rp=.02, ε=1e-5, Cε=.001, Rl=1e-8, +Cf=1, εl=1e-5, εd=1e-8, IP=9와 첫 step q̃0=.01을 사용한다. +힘/모멘트별 일반 규칙과 예외를 함께 구현한다. 보정의 분모는 predictor 차이가 아니라 +직전 accepted endpoint부터 현재 trial까지의 증분이다. + +평균 후보는 현재 trial에서 다시 계산하고, 모든 acceptance가 끝난 증분에서 한 번 commit한다. +외부 field 수렴 후에도 유지 검사가 실패하면 정상 결과를 공개하지 않는다. 자동 완화, +숨은 strict 재시도·polishing, 기존 HDF5 재사용으로 성공을 만들지 않는다. +회전 전용 correction의 기존 strict 선택 자격은 유지하고 후보 field pass와 혼용하지 않는다. + +## 실제 실행 순서와 완료 증거 + +| 순서 | 실제 작업과 주요 대상 | 완료 증거 | +| --- | --- | --- | +| T1 | 공개 규칙·FESA 대응·분기 시점 고정. formulation-convergence.md 작성 | 수식/분기/단위/차이와 독립 수치 검토 | +| T2 | strict source·Release 실행물·DLL 보존, 기준 실행과 개발 계측 | hash·native exit·새 결과·U비교·elapsed·비용 기록 | +| T3 | analysis-private field 수렴 계산과 평균 이력 구현 | 독립 oracle의 경계·예외·rollback 단위 검사 | +| T4 | HHT 실제 종료·accepted 이력·로그 연결 | 실제 조기 종료·후속 검사·최종 파일 보호·v2 reader 검사 | +| T5 | 실제 후보 경로의 물리·시간 정확도 검사 | q/v/physical a·hidden·회전·C08a/C08b·기존 N7 판정 | +| T6 | 후보 Release 공식 측정·개발 계측·정식 HHT 비교 | 기준과 같은 조건의 elapsed/반복/비용·새 HDF5·비교 결과 | +| T7 | 전체 직접 MSVC Debug·CMake/CTest·Python·필수 reference | 전체 discovery/실행 identity·exit·raw/XML·비교 증거 | +| T8 | 결과와 한계·채택 판단 문서 갱신 | results.md와 관련 README/contracts/status/findings의 근거 일치 | + +순서는 T1→T2→T3→T4→T5→T6→T7→T8이다. 코드 의미 변경 전에 기준을 보존하고, +측정 중 다른 build/test를 실행하지 않도록 직렬 수행한다. 새 source/test는 직접 MSVC와 +CMake의 명시적 목록에 함께 등록한다. 같은 source·build tree를 여러 담당자가 동시에 수정하지 않는다. + +## 코드와 검사 범위 + +핵심 변경 위치는 src/fesa/analysis의 hht_dynamic_assembly.cpp, +hht_reduced_procedure.cpp와 관련 private 선언, 필요 시 hht_reduced_operators.cpp다. +새 field 계산은 hht_field_convergence.h/.cpp 등의 최소 private 단위로 둔다. +로그 변경은 공개 newton_iteration_log.h와 대응 구현·검사·contracts에 함께 반영한다. +범용 runtime tolerance 설정이나 새 CLI option, *CONTROLS parser는 추가하지 않는다. + +단위 검사는 상쇄되는 요소 flux, 비활성 영역, 비균일 절점 연결, 첫 평균, +분모0/NaN/Inf, 경계의 아래/같음/위, zero-flux·선형·추정·비이차 분기, 이력 commit/rollback을 포함한다. +통합 검사는 판정량 출력만 달라진 구현을 걸러낼 수 있도록 실제 후보 종료를 확인한다. +CorrectionRecord에서 strict 크기 gate를 분리하면서 PrivateCorrectionEqual까지 없애지 않도록 검사한다. + +hht_reduced_shell_test.cpp, hht_dynamic_procedure_test.cpp, +hht_time_integration_test.cpp, hht_temporal_verification_test.cpp, +hht_shell_physics_test.cpp, hht_energy_verification_test.cpp와 state/writer 검사를 수행한다. +기존 strict 상수를 사용하는 algebra fixture와 실제 후보 production 경로의 검증을 구분한다. + +## 성능 측정 방법 + +공식 elapsed는 계측 없는 baseline/candidate Release 실행을 같은 host·최적화/산술 옵션· +runtime·thread·전원 조건에서 각각1회 측정한다. Process 시작 직전부터 정상 종료·최종 HDF5 +공개까지를 포함하며 build·입력 복사·비교는 제외한다. 동일한 개발 probe를 적용한 별도 실행으로 +outer/inner correction, 요소 평가, factor/solve/RHS 수와 중복 없는 단계 시간을 수집한다. +Probe를 추가하더라도 공식 실행에서는 꺼져야 하며 수학 계산·평가 횟수를 바꾸지 않는다. + +결과는 elapsed, Δ초, 변화율, baseline/candidate speedup과 정확도를 함께 제시한다. +과거6372초·계측 CPU sample 합·여러 실행 중 최솟값을 공식 baseline으로 사용하지 않는다. +기존100초 목표는 참고이며 이번 실험이 그 달성을 약속하는 것은 아니다. + +정식 reference/shellNL-dyn-meshfine의 hash 확인 새 작업 복사본에서 실행한다. +전체729-node/101-frame schema·identity를 먼저 확인한 뒤 기존9-node U blocking/UR warning-only +비교를 수행한다. Near-zero1%·일반 행5%·RMS1%·zero-scale exact-zero를 유지한다. +후보/기준 전체 q/v/a/energy 차이는 별도 진단이며, CSV에 없는 V/a/energy의 외부 검증을 대신하지 않는다. + +## 정확도·전체 회귀·실패의 처리 + +q/v/physical a의 마지막 두 refinement 감소·차수≥1.8, finest normalized error≤1e-4, +기존 N7 해상도 판정을 유지한다. C08b 미판정이나 일반 비선형 HHT의 energy 의미를 +결과에 맞춰 다시 정의하지 않는다. 기존 물리 기준을 통과하지 못하면 수치 부적합으로 보고한다. + +관련 검사 후 전체 MSVC x64 Debug build와 unit/integration/reference 세 executable, +보조 CMake/CTest 전체 및 별도 Python suite를 실제 실행한다. 환경·명령은 +[solver-core](../solver-core/README.md#빌드와-직접-검사)를 따른다. 짧고 고유한 TMP/TEMP, +raw/XML과 discovery를 보존하고 상속 filter/sharding을 제거한다. + +정식 fine·추가 fine·pressure 등 자료가 없으면 미완료다. _02의 전용 비교 계약 미구현도 구분한다. +reference/temp로 대체하거나 fixture 기대값을 완화·삭제하지 않는다. Skip·중단·발견0건은 PASS가 아니다. +기준에서도 재현되는 실패와 후보가 만든 실패를 구분한다. 범위 안 결함은 수정·재검증하고, +유지하기로 한 물리·acceptance 변경이 필요하면 근거와 영향을 보고한다. + +실패 때문에 완결 HDF5가 없으면 성공 실행시간 비교나 reference PASS는 성립하지 않는다. +입증된 같은 수치 실패의 긴 실행을 새 근거 없이 반복하지 않는다. 가능한 독립 검사와 +결과/한계 기록은 수행하되 실행하지 못한 단계를 전체 완료로 표시하지 않는다. + +## 산출물·판정과 실행 권한 + +원시 실행물·DLL·stdout/stderr·XML·HDF5·비교 JSON·환경/입력 hash·elapsed는 +.build/tolerance-experiment//의 기준/후보별 위치에 보존한다. +추적 요약은 docs/abaqus-tolerance/results.md에 작성하고 HHT README·관련 contracts· +tracking/status와 findings 및 필요한 결정 기록을 실제 결과에 맞춘다. + +결론은 **개선 확인 / 속도 미개선 / 수치 부적합 / 검증 미완료**로 구분한다. +개선 확인에는 필수 수치 검증과 동일 조건 후보 elapsed<기준 elapsed가 모두 필요하다. +단일 표본의 개선을 통계적 보장으로 확대하지 않는다. 코드 구현 완료와 전체 검증 완료를 분리한다. + +승인 후 로컬 코드·검사·빌드·FESA 실행은 이 계획의 범위이며 단계별 승인을 다시 요구하지 않는다. +외부 Abaqus·라이선스 비용·원격 전송·release·최종 병합은 제외한다. +실험 revision의 후보 적용과 주 개발 branch의 정식 기본값 채택을 구분하며 최종 병합은 사용자가 결정한다. diff --git a/docs/abaqus-tolerance/local-kernel-experiment.md b/docs/abaqus-tolerance/local-kernel-experiment.md new file mode 100644 index 0000000..3878226 --- /dev/null +++ b/docs/abaqus-tolerance/local-kernel-experiment.md @@ -0,0 +1,132 @@ +# MITC4 local 근사 실험 + +2026-10-18 사용자 명시 지시에 따라 적용한 실험이다. Abaqus의 공식 tolerance나 +외부 Newton 수렴 규칙을 구현한 변경이 아니다. 기준 revision은 `43867d0`이며 +기존 문서와 reference 원본은 보존했다. + +후속 사용자 지시로 테스트 오차 허용치도 `1e-3`으로 변경했다. 아래 1차 결과는 그 전의 +엄격한 테스트 기준으로 실행한 이력이며, 최신 결과는 마지막 절을 따른다. + +## 적용 범위 + +| 위치 | 변경 | +| --- | --- | +| nonlinear_mitc4_shell.cpp / PolarRotation | 고정 12회 → 고정 1회 | +| nonlinear_mitc4_inertia.cpp / ExpSkew, LeftJacobian | order 1..32 → 1..3; 항등항은 별도 유지 | +| 같은 파일 / CheckPhysicalMass | 대칭성·analytic null action·full Gram·Cholesky 복원 네 비교 1e-12 → 1e-3 | + +회전 급수는 `E3=I+P+P²/2+P³/6`, Jacobian은 +`J3=I+P/2+P²/6+P³/24`다. 자동미분은 이 다항식의 D1/D2/D3을 계산한다. +질량의 양의 Cholesky pivot, reciprocal condition `2^-32`, 초기 director 단위길이 +검사는 그대로다. 외부 Newton·inner root·선형 solve·update·전역 평형과 시간수렴, +외부 reference 비교의 기존 허용치는 변경하지 않았다. Hex8도 변경하지 않았다. +공유 polar 변경은 비선형 정적 MITC4와 HHT에 모두 영향을 준다. + +## 수치적 의미와 독립 검토 + +`F=RU`일 때 한 번의 polar 반복 결과는 `R(U+U^-1)/2`이므로 일반 변형에서 +직교 회전행렬이 아니다. Drilling potential과 그 미분값이 바뀐다. 순수 강체회전은 +별개이며 실제 관련 강체운동 검사는 통과했다. + +Cubic exponential은 정확한 SO(3) 회전과 director/Jacobian 미분·null 항등식을 +일반적으로 만족하지 않는다. D3 오차도 작은 회전에서 나타날 수 있다. 변경된 식의 +AD/finite-difference 접선 일치가 원래 유한회전 물리의 정확도를 보장하지 않는다. +질량 검사 완화는 검사값을 통과시키는 범위를 10^9배 넓히며 행렬 오차를 교정하지 않는다. +독립 수치 검토에서도 요청과 production diff가 일치함을 확인했지만 정확도 통과로 +판정하지 않았다. + +## 1차 실행 증거 — 테스트 허용치 변경 전 + +실행 로그·XML·discovery는 ignored `.build/local-numerics-check/`에 있다. +MSVC x64 Debug 전체 build와 CMake/Ninja Debug 전체 build는 각각 exit 0이다. +동일 source의 unit 850·integration 72·reference 43, 총 965개가 두 경로에서 등록됐다. +CTest 표시 이름의 parameter pretty-print를 제외하고 실제 GoogleTest filter identity를 +대조해 동일함을 확인했다. 등록과 실제 실행은 다르다. + +| 실행 | 결과 | 로그 | +| --- | --- | --- | +| 변경 전 새 cubic 검사 | 1 FAIL: 기존 32차와 구분됨 | baseline-requested-test.log/.xml | +| 변경 전 관련 기존 검사 | 45 PASS, 0 FAIL | baseline-elements.log/.xml | +| 변경 후 같은 범위와 새 검사 | 46개 중 41 PASS, 5 FAIL | modified-elements.log/.xml | +| CTest 전체 선택·첫 실패 시 정지 | 43개 실행: 42 PASS, 1 FAIL; 나머지 922개 미실행, exit 8 | ctest.log/.xml | + +실패한 기존 검사와 관측 예는 다음과 같다. 기존 assertion은 완화하거나 삭제하지 않았다. + +| 검사 | 관측 | +| --- | --- | +| NonlinearMitc4Inertia.MatchesIndependentWarpedNodalMomentsEnergyAndPower | 질량 독립 oracle 상대오차 3.653585e-7, 기준 1e-12 | +| DynamicRotation.InverseIdentityAndCurrentAnalyticNullAtFiniteRotation | null 변환의 불필요 성분 9.165932e-4, 기준 1e-12 | +| DynamicRotation.ThirdDerivativesMatchAllPermutationsAndInteriorSweep | D3 oracle 오차 예 1.250000e-6, 기준 1e-12 | +| DynamicRotation.MatchesEightyDigitIndependentBlockExponentialAnchors | D3 oracle 오차 예 1.250000e-6, 기준 1e-12 | +| NonlinearMitc4Residual.ExactDrillingConstraintHasIndependentForceGradientAndHessian | constraint 차이 예 6.755203e-4, 기준 1e-11 | + +새 cubic 동작 검사, 기존 접선 유한차분 검사와 강체운동 검사는 통과했다. +이 통과로 위 실패를 상쇄하지 않는다. + +HHT 추가 선택 실행의 raw 기록은 `modified-hht.log`다. C08a 독립 volume-force/ +power/momentum oracle 검사가 실패했으며 힘 상대오차 예는 9.105479e-7(기준 1e-12)다. +`HhtDynamicProcedure.LargeDirectStepStartsFromAcceptedGeometry`는 outer correction +100회 후 `HHT exhausted 100 corrections`로 실패했다. 이 추가 HHT 선택의 변경 전 +동일 실행은 하지 않았으므로 새 회귀로 단정하는 비교 증거와 구분한다. + +이 HHT 선택은 7개 중 6개 완료(4 PASS·2 FAIL) 후 물리 실패가 확인된 상태에서 +마지막 `UniformTranslationHasNoShellStrainAndIsDeterministicAcrossPartitions` 실행을 +agent가 중단했다(native exit -1). 마지막 1개는 미판정이다. 완결 XML은 없으며 raw 로그와 +`hht-interruption.txt`를 남겼다. CTest도 전체 통과가 아니며 첫 실패 C08a에서 중단했다. +전체 Native unit/integration/reference 실행, 정식 외부 reference 비교는 수행하지 않았다. +물리 정확도 실패를 남긴 채 장시간 검증·공식 성능 측정으로 채택을 주장하지 않는다. +선택 범위·exit·source/executable hash 요약은 `verification-summary.json`에 있다. + +Python suite는 현재 Python 환경의 `No module named pytest`로 실행되지 않았다. +자동 설치하지 않았다. 전체 검사와 정식 reference 비교·Release 성능·C08a/C08b +시간수렴 통과는 확보하지 않았다. 검사 실행시간 차이는 조건을 통제한 성능 측정이 아니므로 +속도 향상률로 보고하지 않는다. + +## 1차 판정과 다음 접근 + +요청한 수치 변경은 구현됐지만 기존 정확도 검증은 실패했다. 승인된 실험값과 기존 +검증 기준을 그대로 남긴다. 정확도 손실을 수용한 solver 채택이나 성능 목표 달성을 +뜻하지 않는다. 후속 원인 분리가 필요하면 polar, cubic, mass gate를 하나씩 바꾸는 +동일 조건 비교로 각각의 영향을 측정한다. 실제 성능은 동일 local kernel을 사용한 +strict/Abaqus 외부 정책 비교와 원래 kernel 대비 비교를 분리해 측정해야 한다. + +## 후속 변경 — 테스트 오차 허용치 1e-3 + +사용자가 처음 제안한 `1e-4`를 최종 `1e-3`으로 대체했다. 다음 세 테스트 파일에 +`kNumericalTolerance=1e-3`을 두고 기존 숫자 기준이 이보다 작은 직접 오차 비교 100개 +(관성 50·접선 21·HHT 29)와 요소 유한차분 sweep의 오차 판정을 변경했다. + +- tests/unit/elements/nonlinear_mitc4_inertia_test.cpp +- tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp +- tests/unit/analysis/hht_shell_physics_test.cpp + +오차 계산식과 oracle·fixture·유한차분 간격은 그대로다. 절대오차/정규화오차 구분도 +그대로이므로 모든 `1e-3`을 상대 0.1%로 해석할 수 없다. 일부 작은 회전질량의 절대 +기대값 약 `6.67e-4`보다 허용치가 커서 해당 개별 비교는 회전질량 0도 허용한다. +정확한 equality, nonfinite/실패 검사, 누락항 negative control, 순수 산술 bookkeeping +두 검사, 시간수렴의 감소/차수 조건, 테스트 내부 해석 controller와 production 종료 +조건은 유지했다. HHT 시간수렴의 finest 오차 assertion만 `1e-4`에서 `1e-3`으로 바뀌었다. +외부 Abaqus 결과 비교나 다른 요소의 테스트는 바꾸지 않았다. + +최신 로그는 ignored `.build/test-tolerance-1e-3/`에 별도 보관했다. 두 Debug 전체 target +build는 exit 0이다. Native 및 CTest에서 같은 50개 identity를 선택해 모두 실행했으며 +각각 **46 PASS·4 FAIL·0 SKIP**(exit 1/8)이다. 전체 965개 중 나머지 915개는 이번에 +실행하지 않았다. 관련 요소 46개에 HHT 네 물리 검사를 추가한 선택이다. + +| 남은 실패 | 새 기준에서도 실패한 이유/관측 예 | +| --- | --- | +| DynamicRotation.ThirdDerivativesMatchAllPermutationsAndInteriorSweep | 독립 D3 oracle 절대오차 0.2276159 > 0.001 | +| DynamicRotation.MatchesEightyDigitIndependentBlockExponentialAnchors | 독립 D3 oracle 절대오차 0.2276159 > 0.001 | +| NonlinearMitc4Residual.ExactDrillingConstraintHasIndependentForceGradientAndHessian | gradient 절대오차 0.00329135 > 0.001 | +| HhtShellPhysics.C08aRigidOrbitMatchesIndependentVolumeForcePowerAndMomentum | 정규화 hidden defect 예 0.00123371 > 0.001 및 실제 CompleteEquilibrium 실패 | + +기존 실패 중 질량/운동량 독립 oracle와 inverse/null 검사는 완화된 기준에서 통과했다. +이는 계산값의 정확도가 개선됐다는 뜻이 아니다. D3·drilling·HHT 실패는 남겼으며 +추가 완화나 production 변경으로 숨기지 않았다. 이전 큰 시간간격 수렴 실패와 중단된 +HHT 실행은 이번에 재실행하지 않았고 이전 미판정 상태를 유지한다. + +`native-selected.log/.xml`, `ctest-selected.log/.xml`, `changed-comparisons.json`과 +`verification-summary.json`에 실행·변경 근거가 있다. 최초 CTest wrapper 호출의 shell +인용 오류(exit 255)는 `ctest-invocation-error.log`로 보존하고 CTest 직접 호출로 위 +50개를 재실행했다. 독립 검토에서도 변경 범위와 보존 조건이 일치함을 확인했다. +Production 수식·종료 조건은 직전 실험본 그대로이며 이번에는 테스트와 문서만 수정했다. diff --git a/docs/hht-implicit-dynamics/README.md b/docs/hht-implicit-dynamics/README.md index 30bb178..66200b9 100644 --- a/docs/hht-implicit-dynamics/README.md +++ b/docs/hht-implicit-dynamics/README.md @@ -2,6 +2,14 @@ ## 범위와 현재 상태 +2026-10-18 후속 사용자 지시로 polar 반복 12→1회, 관성 회전/Jacobian 급수 32→3차, +질량 대칭·null/Gram·분해 복원 검사 1e-12→1e-3의 local 근사 실험을 적용했다. +최초 엄격한 기준의 관련 요소 46개 중 5개가 실패했다. 후속 사용자 지시로 관성·접선· +HHT 물리 테스트의 오차 기준을 1e-3으로 완화했으며 최신 선택 50개는 두 빌드 경로에서 +46 PASS·4 FAIL이다. 실제 해석 루틴의 수렴 조건은 이번 테스트 수정으로 바뀌지 않았다. +아래 정식화의 유한회전 항등식과 과거 검증을 이 근사본의 통과로 승계하지 않는다. +변경 범위와 실제 검증 한계는 [실험 기록](../abaqus-tolerance/local-kernel-experiment.md)을 따른다. + MITC4 shell의 유한회전 기하비선형 동해석이다. S4/S4R은 FESA-MITC4로 해석하며 외부 solver 내부 알고리즘과의 동등성을 뜻하지 않는다. 현재 채택 경로는 **nodal-lumped-kinetic 질량, reduced-coordinate HHT, 초기 null의 고정 직교 보완기저, diff --git a/docs/nonlinear-static-newton-raphson/README.md b/docs/nonlinear-static-newton-raphson/README.md index 771156a..e80d2ae 100644 --- a/docs/nonlinear-static-newton-raphson/README.md +++ b/docs/nonlinear-static-newton-raphson/README.md @@ -2,6 +2,12 @@ ## 목적과 경계 +2026-10-18 후속 사용자 지시의 [local 근사 실험](../abaqus-tolerance/local-kernel-experiment.md)에서 +공유 nonlinear MITC4의 polar 반복을 12→1회로 변경했다. 후속 사용자 지시로 관련 접선 +테스트 오차를 1e-3으로 완화한 뒤에도 drilling constraint의 gradient 항등식 검사가 +실패했으므로 과거 정적 검증을 변경본에 승계하지 않는다. +정적 Newton tolerance 자체는 변경하지 않았다. + 등방 선형탄성 MITC4의 기하 비선형 정적 응답을 full Newton으로 계산한다. 재료 비선형·contact·arc length·line search·자동 cutback은 이 기능에 포함하지 않는다. 고정 global dead nodal force/moment와 factor-one prescribed target을 공통 load factor로 증가시킨다. diff --git a/docs/tracking/findings.md b/docs/tracking/findings.md index dcf7c06..4f942c4 100644 --- a/docs/tracking/findings.md +++ b/docs/tracking/findings.md @@ -1,5 +1,20 @@ # 해결되지 않은 문제 +## MITC4 local 근사 실험의 정확도 실패 + +- 조건·관측: 2026-10-18 사용자 지시로 polar 1회·관성 회전/Jacobian 3차·질량 검사 + 1e-3을 적용했다. 변경 전 관련 기존 45개 PASS, 변경 후 46개 중 41 PASS·5 FAIL이다. + 추가 HHT 선택에서도 C08a 독립 물리 oracle와 큰 시간간격 100회 보정 종료가 실패했다. + 후속 사용자 지시로 관련 테스트 오차 기준도 1e-3으로 변경했으며 최신 동일 선택 + 50개는 Native/CTest 각각46 PASS·4 FAIL이다. D3 오차 약0.228, drilling gradient + 약0.00329 및 C08a hidden defect/실제 평형 검사 실패가 남는다. +- 영향: 회전 미분·analytic null·질량/에너지·drilling 항등식의 기존 정확도 검증을 만족하지 + 못한다. 구현 완료를 solver correctness·시간수렴·성능 향상으로 보고할 수 없다. +- 지금 해결하지 않는 이유: 요청한 근사값의 수치 영향이며 승인된 테스트 허용치1e-3보다 + 더 완화하거나 요청을 자동 취소해 실패를 숨기지 않는다. 외부 Newton 기준 변경만으로 해당 항등식은 복원되지 않는다. +- 다음 접근: 원인 분리 실험과 정확도 손실의 별도 판정이 필요하다. 변경 범위·실제 로그· + 미검증 항목은 [local 실험 기록](../abaqus-tolerance/local-kernel-experiment.md)을 따른다. + ## 비정적 검증 공백 - 조건·관측: 2026-10-18 통합에서 장시간 HHT 검증이 사용자 요청으로 중단되고 정적·공통 diff --git a/docs/tracking/status.md b/docs/tracking/status.md index c8e4245..a425397 100644 --- a/docs/tracking/status.md +++ b/docs/tracking/status.md @@ -1,5 +1,12 @@ # 구현과 검증 상태 +2026-10-18 후속 [MITC4 local 근사 실험](../abaqus-tolerance/local-kernel-experiment.md)은 +polar 1회·관성 회전/Jacobian 3차·질량 검사 1e-3을 적용했다. 두 Debug build는 성공했으나 +최초 관련 요소 46개 중 5 FAIL이었다. 후속 사용자 지시로 관련 테스트의 오차 기준을 +1e-3으로 변경했고, 최신 Native/CTest 동일 선택 50개는 각각46 PASS·4 FAIL이다. +전체 965개 중 나머지 915개는 이번에 미실행이며 실제 해석 수렴 조건은 변경하지 않았다. +아래 과거 검증을 이 변경본에 승계하지 않는다. Abaqus 외부 Newton 정책은 미구현이다. + ## 기준 시점 2026-10-18 저장소 HEAD `945d899` 기준이다. 아래 통과 수치는 각 검증 당시 기록의 범위이며 diff --git a/src/fesa/elements/nonlinear_mitc4_inertia.cpp b/src/fesa/elements/nonlinear_mitc4_inertia.cpp index b0716d4..9a806b2 100644 --- a/src/fesa/elements/nonlinear_mitc4_inertia.cpp +++ b/src/fesa/elements/nonlinear_mitc4_inertia.cpp @@ -229,7 +229,7 @@ JetMatrix LeftJacobian(const JetVector& phi) { JetMatrix result = IdentityMatrix(); JetMatrix term = IdentityMatrix(); double factorial = 1.0; - for (std::size_t order = 1U; order <= 32U; ++order) { + for (std::size_t order = 1U; order <= 3U; ++order) { term = Multiply(term, skew); factorial *= static_cast(order + 1U); for (std::size_t row = 0U; row < 3U; ++row) { @@ -242,14 +242,13 @@ JetMatrix LeftJacobian(const JetVector& phi) { return result; } -/// @brief Uses a converged exponential series so D1/D2/D3 share one analytic -/// map. +/// @brief Uses a cubic exponential polynomial with matching D1/D2/D3. JetMatrix ExpSkew(const JetVector& phi) { const JetMatrix skew = Skew(phi); JetMatrix result = IdentityMatrix(); JetMatrix term = IdentityMatrix(); double factorial = 1.0; - for (std::size_t order = 1U; order <= 32U; ++order) { + for (std::size_t order = 1U; order <= 3U; ++order) { term = Multiply(term, skew); factorial *= static_cast(order); for (std::size_t row = 0U; row < 3U; ++row) { @@ -378,9 +377,9 @@ Status CheckPhysicalMass(const Matrix& mass, const double full_gram_error = std::hypot( gram_error, std::hypot(null_gram_error, std::sqrt(2.0) * orthogonality_error)); - if (!(symmetry_norm / mass_norm <= 1.0e-12) || - !(null_action_norm / mass_norm / 2.0 <= 1.0e-12) || - !(full_gram_error <= 1.0e-12)) { + if (!(symmetry_norm / mass_norm <= 1.0e-3) || + !(null_action_norm / mass_norm / 2.0 <= 1.0e-3) || + !(full_gram_error <= 1.0e-3)) { return fail("invalid-nonlinear-inertia-null", "Scaled mass symmetry and analytic null identities failed."); } @@ -408,7 +407,7 @@ Status CheckPhysicalMass(const Matrix& mass, block_frobenius = std::hypot(block_frobenius, block(row, column)); } } - if (!(reconstruction_error / block_frobenius <= 1.0e-12)) { + if (!(reconstruction_error / block_frobenius <= 1.0e-3)) { return fail("invalid-nonlinear-inertia-positive-mass", "Physical mass Cholesky reconstruction failed."); } diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp index d058f18..d6711cb 100644 --- a/src/fesa/elements/nonlinear_mitc4_shell.cpp +++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp @@ -1155,11 +1155,10 @@ Jet PhysicalEnergy(const std::array& coordinates, return energy; } -/// @brief Computes a fixed-iteration analytic polar factor for drill -/// corotation. +/// @brief Applies one polar iteration for drill corotation. JetMatrix3 PolarRotation(const JetMatrix3& deformation) { JetMatrix3 result = deformation; - for (std::size_t iteration = 0U; iteration < 12U; ++iteration) { + for (std::size_t iteration = 0U; iteration < 1U; ++iteration) { const JetMatrix3 coefficients = Transpose(result); Matrix primal{3, 3}; for (std::size_t row = 0; row < 3; ++row) diff --git a/tests/unit/analysis/hht_shell_physics_test.cpp b/tests/unit/analysis/hht_shell_physics_test.cpp index 13b49c9..abac0f0 100644 --- a/tests/unit/analysis/hht_shell_physics_test.cpp +++ b/tests/unit/analysis/hht_shell_physics_test.cpp @@ -33,6 +33,9 @@ namespace { +// User-approved physical/numerical comparison tolerance for the local experiment. +constexpr double kNumericalTolerance = 1.0e-3; + using Frame = fesa::DynamicAcceptedFrameState; using Assembly = fesa::hht_internal::DynamicAssembly; using L3 = std::array; @@ -520,7 +523,7 @@ void CompleteManufacturedConstraints( if (scale == 0.0) EXPECT_DOUBLE_EQ(value, 0.0); else - EXPECT_LE(value / scale, 1.0e-10); + EXPECT_LE(value / scale, kNumericalTolerance); } } } @@ -561,7 +564,7 @@ TEST(HhtShellPhysics, NodalRigidForceRatesMatchIndependentTimeDifferences) { scale = std::hypot(scale, difference); } ASSERT_GT(scale, 0.0); - EXPECT_LE(error / scale, 1.0e-6); + EXPECT_LE(error / scale, kNumericalTolerance); } } } @@ -580,10 +583,10 @@ TEST(HhtShellPhysics, ExternalDrillingHessianBalancesGradientAcceleration) { for (std::size_t component = 0U; component < 3U; ++component) gradient_acceleration += external.Value().gradient[component] * acceleration[component]; - EXPECT_NEAR(external.Value().directional_second, -1.0 / 6.0, 1.0e-14); - EXPECT_NEAR(gradient_acceleration, 1.0 / 6.0, 1.0e-14); + EXPECT_NEAR(external.Value().directional_second, -1.0 / 6.0, kNumericalTolerance); + EXPECT_NEAR(gradient_acceleration, 1.0 / 6.0, kNumericalTolerance); EXPECT_NEAR(external.Value().directional_second + gradient_acceleration, - 0.0, 1.0e-14); + 0.0, kNumericalTolerance); constexpr double step = 1.0e-4; const auto plus = fixture.shells[0U].internal->EvaluateExternalDrillingConstraint( 0U, {step, 0.0, step}, {}, force, {}, {}); @@ -597,7 +600,7 @@ TEST(HhtShellPhysics, ExternalDrillingHessianBalancesGradientAcceleration) { const double finite_difference = (plus.Value().value - 2.0 * middle.Value().value + minus.Value().value) / (step * step); - EXPECT_NEAR(finite_difference, external.Value().directional_second, 1.0e-8); + EXPECT_NEAR(finite_difference, external.Value().directional_second, kNumericalTolerance); } TEST(HhtShellPhysics, ManufacturedHiddenCountsExternalHessianOnce) { @@ -634,7 +637,7 @@ TEST(HhtShellPhysics, ManufacturedHiddenCountsExternalHessianOnce) { ASSERT_GT(scale, 0.0); EXPECT_LE(std::abs(balance.hidden_acceleration_defect[node] - expected) / scale, - 1.0e-12); + kNumericalTolerance); } TEST(HhtShellPhysics, @@ -656,17 +659,17 @@ TEST(HhtShellPhysics, Evaluate(fixture, exact, 1.0 / 16.0, parameters, true); auto difference = assembly.frame.inertia_force; difference.Axpy(-1.0, oracle.force); - EXPECT_LE(difference.Norm() / oracle.force.Norm(), 1.0e-12); + EXPECT_LE(difference.Norm() / oracle.force.Norm(), kNumericalTolerance); EXPECT_LE(std::abs(assembly.frame.kinetic_energy - static_cast(oracle.kinetic)) / static_cast(oracle.kinetic), - 1.0e-12); + kNumericalTolerance); EXPECT_LE(std::abs(oracle.force.Dot(exact.velocity) - static_cast(oracle.power)) / static_cast(oracle.power_scale), - 1.0e-10); - EXPECT_LE(assembly.frame.physical_strain_energy / 800.0, 1.0e-10); - EXPECT_LE(assembly.frame.drilling_energy / 800.0, 1.0e-10); + kNumericalTolerance); + EXPECT_LE(assembly.frame.physical_strain_energy / 800.0, kNumericalTolerance); + EXPECT_LE(assembly.frame.drilling_energy / 800.0, kNumericalTolerance); fesa::Vector scaled_q = exact.generalized_coordinates; fesa::Vector scaled_force = assembly.frame.internal_force; for (std::size_t i = 0; i < scaled_q.Size(); ++i) { @@ -679,7 +682,7 @@ TEST(HhtShellPhysics, if (force_scale == 0.0) EXPECT_DOUBLE_EQ(scaled_force.Norm(), 0.0); else - EXPECT_LE(scaled_force.Norm() / force_scale, 1.0e-10); + EXPECT_LE(scaled_force.Norm() / force_scale, kNumericalTolerance); double momentum_error = 0.0; double angular_error = 0.0; for (std::size_t c = 0; c < 3; ++c) { @@ -691,10 +694,10 @@ TEST(HhtShellPhysics, static_cast(oracle.angular_momentum[c])); } EXPECT_LE(momentum_error / static_cast(oracle.momentum_scale), - 1.0e-12); + kNumericalTolerance); EXPECT_LE( angular_error / static_cast(oracle.angular_momentum_scale), - 1.0e-12); + kNumericalTolerance); auto balance = assembly.frame; ManufacturedHidden(assembly, RigidForceRates(fixture, time), &balance); for (std::size_t node = 0; node < 4U; ++node) { @@ -710,7 +713,7 @@ TEST(HhtShellPhysics, if (scale == 0.0) EXPECT_DOUBLE_EQ(value, 0.0); else - EXPECT_LE(value / scale, 1.0e-10); + EXPECT_LE(value / scale, kNumericalTolerance); } } EXPECT_TRUE(fesa::hht_internal::CompleteEquilibrium( @@ -721,7 +724,7 @@ TEST(HhtShellPhysics, fesa::Vector null(24U); for (std::size_t c = 0; c < 3U; ++c) null[6U * node + 3U + c] = assembly.rotations[node].drill_null[c]; - EXPECT_LE(assembly.mass.Multiply(null).Norm() / mass_norm, 1.0e-12); + EXPECT_LE(assembly.mass.Multiply(null).Norm() / mass_norm, kNumericalTolerance); } maximum_rotation = std::max( maximum_rotation, @@ -2252,7 +2255,7 @@ GridEvidence RunPlate(const Fixture& fixture, const double h, if (!manufactured) { InitializeNullVelocity(fixture, assembly, &frame, &tangent); for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node) - EXPECT_NEAR(frame.velocity[6U * node + 5U], 0.4, 1.0e-12); + EXPECT_NEAR(frame.velocity[6U * node + 5U], 0.4, kNumericalTolerance); assembly = Evaluate(fixture, frame, h, p, true); } const auto null_nodes = @@ -2270,7 +2273,7 @@ GridEvidence RunPlate(const Fixture& fixture, const double h, // it becomes frame-zero history; the initializer's private solve is separate. for (const double residual : frame.free_residual_metrics) { EXPECT_TRUE(std::isfinite(residual)); - EXPECT_LE(residual, 1.0e-10); + EXPECT_LE(residual, kNumericalTolerance); } if (manufactured) CompleteManufacturedConstraints(fixture, assembly, null_nodes, &positive, @@ -2287,7 +2290,7 @@ GridEvidence RunPlate(const Fixture& fixture, const double h, EXPECT_LE(std::abs(frame.kinetic_energy - static_cast(initial_oracle.kinetic)) / energy0, - 1.0e-12); + kNumericalTolerance); const auto initial_momentum = frame.linear_momentum; const auto initial_angular_momentum = frame.angular_momentum; evidence.history.push_back(frame); @@ -2356,7 +2359,7 @@ GridEvidence RunPlate(const Fixture& fixture, const double h, Require(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs, reconstructed_assembly, &reconstructed_zero)); for (const double residual : reconstructed_zero.free_residual_metrics) - EXPECT_LE(residual, 1.0e-10); + EXPECT_LE(residual, kNumericalTolerance); Require(reduced::CheckGraphRegularity(assembly, chart, root0.operators, tangent)); const auto increments = static_cast(evidence.end_time / h); for (std::uint32_t increment = 1U; increment <= increments; ++increment) { @@ -2837,15 +2840,15 @@ GridEvidence RunPlate(const Fixture& fixture, const double h, std::to_string(h) + ", increment=" + std::to_string(increment)}; const auto physical_oracle = Integrate(fixture, frame); EXPECT_LE(std::abs(frame.kinetic_energy - - static_cast(physical_oracle.kinetic)) / energy0, 1.0e-10); + static_cast(physical_oracle.kinetic)) / energy0, kNumericalTolerance); EXPECT_LE(static_cast(Norm(Add(physical_oracle.momentum, Scale(L3{frame.linear_momentum[0], frame.linear_momentum[1], frame.linear_momentum[2]}, -1.0L)))) / - static_cast(initial_oracle.momentum_scale), 1.0e-10); + static_cast(initial_oracle.momentum_scale), kNumericalTolerance); EXPECT_LE(static_cast(Norm(Add(physical_oracle.angular_momentum, Scale(L3{frame.angular_momentum[0], frame.angular_momentum[1], frame.angular_momentum[2]}, -1.0L)))) / - static_cast(initial_oracle.angular_momentum_scale), 1.0e-10); + static_cast(initial_oracle.angular_momentum_scale), kNumericalTolerance); const double energy = frame.kinetic_energy + frame.physical_strain_energy + frame.drilling_energy; if (!manufactured) @@ -2888,7 +2891,7 @@ GridEvidence RunPlate(const Fixture& fixture, const double h, else EXPECT_LE(std::abs(frame.hidden_constraint[i]) / frame.hidden_constraint_scale[i], - 1.0e-10); + kNumericalTolerance); const double vscale = frame.hidden_velocity_defect_scale[i]; const double ascale = frame.hidden_acceleration_defect_scale[i]; if (vscale > 0.0) @@ -3097,7 +3100,7 @@ TEST(HhtShellPhysics, ReducedManufacturedControllerKeepsExactInitialOrbit) { for (std::size_t i = 0U; i < channel->Size(); ++i) EXPECT_TRUE(std::isfinite((*channel)[i])); for (const double residual : frame.free_residual_metrics) - EXPECT_LE(residual, 1.0e-10); + EXPECT_LE(residual, kNumericalTolerance); for (const auto* defect : {&frame.hidden_constraint, &frame.hidden_velocity_defect, &frame.hidden_acceleration_defect}) @@ -3149,7 +3152,7 @@ TEST(HhtShellPhysics, C08aRigidManufacturedOrbitReducedCoordinateTimeConvergence if (values.second == 0.0) EXPECT_DOUBLE_EQ(values.first, 0.0); else - EXPECT_LE(std::abs(values.first) / values.second, 1.0e-10); + EXPECT_LE(std::abs(values.first) / values.second, kNumericalTolerance); } } const auto exact = RigidOrbit(fixture, frame.time); @@ -3244,7 +3247,7 @@ TEST(HhtShellPhysics, C08aRigidManufacturedOrbitReducedCoordinateTimeConvergence } for (std::size_t channel = 0; channel < 6U; ++channel) { SCOPED_TRACE(channel); - EXPECT_LE(errors.back()[channel], 1.0e-4); + EXPECT_LE(errors.back()[channel], kNumericalTolerance); for (std::size_t grid = grids.size() - 3U; grid + 1U < grids.size(); ++grid) { if (errors[grid][channel] == 0.0 && errors[grid + 1U][channel] == 0.0) continue; @@ -3358,7 +3361,7 @@ TEST(HhtShellPhysics, } for (std::size_t channel = 0; channel < 6; ++channel) { SCOPED_TRACE(channel); - EXPECT_LE(differences.back()[channel], 1.0e-4); + EXPECT_LE(differences.back()[channel], kNumericalTolerance); for (std::size_t pair = differences.size() - 3U; pair + 1U < differences.size(); ++pair) { const double first = differences[pair][channel]; diff --git a/tests/unit/elements/nonlinear_mitc4_inertia_test.cpp b/tests/unit/elements/nonlinear_mitc4_inertia_test.cpp index 79f5761..5002648 100644 --- a/tests/unit/elements/nonlinear_mitc4_inertia_test.cpp +++ b/tests/unit/elements/nonlinear_mitc4_inertia_test.cpp @@ -18,6 +18,9 @@ namespace { +// User-approved physical/numerical comparison tolerance for the local experiment. +constexpr double kNumericalTolerance = 1.0e-3; + fesa::Node Node(const std::int64_t label, const std::array& coordinates) { return {{"Shell-Instance", label, std::to_string(label)}, @@ -88,17 +91,17 @@ TEST(NonlinearMitc4Inertia, IntegratesTheRectangleNodalLumpedMassAndNull) { for (std::size_t row = 0U; row < 4U; ++row) { for (std::size_t column = 0U; column < 4U; ++column) { EXPECT_NEAR(mass(6U * row, 6U * column), row == column ? 0.2 : 0.0, - 1.0e-12); + kNumericalTolerance); EXPECT_NEAR(mass(6U * row + 3U, 6U * column + 3U), - row == column ? 0.2 * 0.2 * 0.2 / 12.0 : 0.0, 1.0e-12); + row == column ? 0.2 * 0.2 * 0.2 / 12.0 : 0.0, kNumericalTolerance); } - EXPECT_NEAR(mass(6U * row + 5U, 6U * row + 5U), 0.0, 1.0e-12); + EXPECT_NEAR(mass(6U * row + 5U, 6U * row + 5U), 0.0, kNumericalTolerance); } fesa::Vector drilling_null{24U}; for (std::size_t node = 0U; node < 4U; ++node) { drilling_null[6U * node + 5U] = 1.0; } - EXPECT_LE(mass.Multiply(drilling_null).Norm(), 1.0e-12); + EXPECT_LE(mass.Multiply(drilling_null).Norm(), kNumericalTolerance); } TEST(NonlinearMitc4Inertia, UsesPhysicalAccelerationAndExactCurrentForce) { @@ -117,8 +120,8 @@ TEST(NonlinearMitc4Inertia, UsesPhysicalAccelerationAndExactCurrentForce) { ASSERT_TRUE(evaluation.HasValue()); const fesa::Vector expected = evaluation.Value().mass.Multiply(acceleration); for (std::size_t index = 0U; index < expected.Size(); ++index) { - EXPECT_NEAR(evaluation.Value().force[index], expected[index], 1.0e-11); - EXPECT_NEAR(evaluation.Value().convective_force[index], 0.0, 1.0e-11); + EXPECT_NEAR(evaluation.Value().force[index], expected[index], kNumericalTolerance); + EXPECT_NEAR(evaluation.Value().convective_force[index], 0.0, kNumericalTolerance); } } @@ -170,7 +173,7 @@ TEST(NonlinearMitc4Inertia, EXPECT_DOUBLE_EQ(actual.Value().convective_force[full], 0.0); } ASSERT_GT(scale, 0.0); - EXPECT_LE(error / scale, 1e-12); + EXPECT_LE(error / scale, kNumericalTolerance); } } } @@ -234,8 +237,8 @@ TEST(NonlinearMitc4Inertia, MixedScaleRotationDerivativeAndAngularMomentum) { std::abs(actual.Value().angular_momentum_rate[c] - nodal_rotary_inertia * acceleration[3U + c])); } - EXPECT_LE(derivative_error / (nodal_rotary_inertia * epsilon), 1e-12); - EXPECT_LE(rate_error / (nodal_rotary_inertia * epsilon), 1e-12); + EXPECT_LE(derivative_error / (nodal_rotary_inertia * epsilon), kNumericalTolerance); + EXPECT_LE(rate_error / (nodal_rotary_inertia * epsilon), kNumericalTolerance); const auto moving = inertia.Value().Evaluate( {}, acceleration, fesa::Vector{24U}); ASSERT_TRUE(moving.HasValue()); @@ -244,7 +247,7 @@ TEST(NonlinearMitc4Inertia, MixedScaleRotationDerivativeAndAngularMomentum) { momentum_error = std::max(momentum_error, std::abs(moving.Value().angular_momentum[c] - nodal_rotary_inertia * acceleration[3U + c])); - EXPECT_LE(momentum_error / (nodal_rotary_inertia * epsilon), 1e-12); + EXPECT_LE(momentum_error / (nodal_rotary_inertia * epsilon), kNumericalTolerance); } } } @@ -490,6 +493,19 @@ fesa::Matrix ScaledMatrix(const fesa::Matrix& matrix, const double length) { return result; } +TEST(DynamicRotation, UsesCubicDirectorAndJacobianAtUnitAxisRotation) { + const auto value = fesa::EvaluateDynamicRotation( + {1.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, {0.0, 1.0, 0.0}, {}); + ASSERT_TRUE(value.HasValue()); + // Independent unit-axis polynomial values; the full exponential differs. + EXPECT_NEAR(value.Value().director.X(), 0.0, kNumericalTolerance); + EXPECT_NEAR(value.Value().director.Y(), -5.0 / 6.0, kNumericalTolerance); + EXPECT_NEAR(value.Value().director.Z(), 0.5, kNumericalTolerance); + EXPECT_NEAR(value.Value().angular_velocity.X(), 0.0, kNumericalTolerance); + EXPECT_NEAR(value.Value().angular_velocity.Y(), 5.0 / 6.0, kNumericalTolerance); + EXPECT_NEAR(value.Value().angular_velocity.Z(), 11.0 / 24.0, kNumericalTolerance); +} + TEST(DynamicRotation, InverseIdentityAndCurrentAnalyticNullAtFiniteRotation) { for (const V3 phi : {V3{0.5, 0.0, 0.0}, V3{0.2, -0.3, 0.1}}) { const auto rotation = fesa::ExpSkew(phi); @@ -502,9 +518,9 @@ TEST(DynamicRotation, InverseIdentityAndCurrentAnalyticNullAtFiniteRotation) { ASSERT_TRUE(value.HasValue()); const V3 mapped = Multiply(jacobian, value.Value().drill_null); EXPECT_GT(mapped[column], 0.0); - EXPECT_NEAR(value.Value().drill_null.Norm(), 1.0, 1e-12); + EXPECT_NEAR(value.Value().drill_null.Norm(), 1.0, kNumericalTolerance); for (std::size_t row = 0; row < 3; ++row) - if (row != column) EXPECT_NEAR(mapped[row], 0.0, 1e-12); + if (row != column) EXPECT_NEAR(mapped[row], 0.0, kNumericalTolerance); } } for (const Fixture fixture : {Fixture{}, Warped()}) { @@ -527,9 +543,9 @@ TEST(DynamicRotation, InverseIdentityAndCurrentAnalyticNullAtFiniteRotation) { for (std::size_t j = 0; j < 3; ++j) { n[6U * node + 3U + j] = rotation.Value().drill_null[j]; } - EXPECT_LE(mass.Multiply(n).Norm() / MatrixNorm(mass), 1.0e-12); + EXPECT_LE(mass.Multiply(n).Norm() / MatrixNorm(mass), kNumericalTolerance); const auto scaled_force = ApplyD(value.Value().force, fixture.Length()); - EXPECT_LE(std::abs(n.Dot(scaled_force)) / scaled_force.Norm(), 1.0e-10); + EXPECT_LE(std::abs(n.Dot(scaled_force)) / scaled_force.Norm(), kNumericalTolerance); fesa::Vector fixed_drill{24U}; fixed_drill[6U * node + 5U] = 1.0; EXPECT_GT(mass.Multiply(fixed_drill).Norm() / MatrixNorm(mass), 1.0e-8); @@ -630,7 +646,7 @@ TEST(NonlinearMitc4Inertia, ? analytical.Norm() + fd.Norm() : MatrixNorm(scaled); const double metric = scale > 0.0 ? error / scale : error; - const double tolerance = analytical.Norm() > 0.0 ? 1.0e-6 : 1.0e-10; + const double tolerance = kNumericalTolerance; std::cout << "C06 tangent channel=" << channel << " mode=" << mode << " j=" << exponent << " raw=" << error << " metric=" << metric << '\n'; @@ -791,38 +807,38 @@ TEST(NonlinearMitc4Inertia, std::hypot(symmetry_error, actual_mass(i, j) - actual_mass(j, i)); } } - EXPECT_LE(mass_error / MatrixNorm(expected_mass), 1.0e-12); - EXPECT_LE(symmetry_error / MatrixNorm(expected_mass), 1.0e-12); + EXPECT_LE(mass_error / MatrixNorm(expected_mass), kNumericalTolerance); + EXPECT_LE(symmetry_error / MatrixNorm(expected_mass), kNumericalTolerance); EXPECT_LE( std::abs(value.Value().kinetic_energy - oracle.energy) / oracle.energy, - 1.0e-12); + kNumericalTolerance); EXPECT_LE(std::abs(0.5 * state.v.Dot(value.Value().mass.Multiply(state.v)) - oracle.energy) / oracle.energy, - 1.0e-12); + kNumericalTolerance); EXPECT_LE(std::abs(state.v.Dot(value.Value().force) - oracle.power) / std::abs(oracle.power), - 1.0e-10); + kNumericalTolerance); EXPECT_LE(Difference(ApplyD(value.Value().force, fixture.Length()), ApplyD(oracle.force, fixture.Length()), 1.0) .Norm() / ApplyD(oracle.force, fixture.Length()).Norm(), - 1.0e-12); + kNumericalTolerance); EXPECT_LE((value.Value().momentum - oracle.momentum).Norm() / oracle.momentum.Norm(), - 1.0e-12); + kNumericalTolerance); EXPECT_LE( (value.Value().angular_momentum - oracle.angular_momentum).Norm() / oracle.angular_momentum.Norm(), - 1.0e-12); + kNumericalTolerance); EXPECT_LE((value.Value().momentum_rate - oracle.momentum_rate).Norm() / oracle.momentum_rate.Norm(), - 1.0e-12); + kNumericalTolerance); EXPECT_LE( (value.Value().angular_momentum_rate - oracle.angular_momentum_rate) .Norm() / oracle.angular_momentum_rate.Norm(), - 1.0e-12); + kNumericalTolerance); // Recover C0/C1/C2 from their known tensor factors, with the prescribed // m_e*(t/2)^k scales even when a centered C1 is exactly zero. for (std::size_t i = 0; i < 4; ++i) { @@ -854,7 +870,7 @@ TEST(NonlinearMitc4Inertia, EXPECT_LE(std::abs(numerator[k] / denominator[k] - oracle.moments[k][i][j]) / scale, - 1.0e-12); + kNumericalTolerance); scale *= fixture.thickness / 2.0; } } @@ -868,7 +884,7 @@ TEST(NonlinearMitc4Inertia, const auto translation = inertia.Value().Evaluate({}, uniform.v, uniform.a); ASSERT_TRUE(translation.HasValue()); EXPECT_NEAR(translation.Value().kinetic_energy / oracle.total_mass, 0.25, - 1.0e-12); + kNumericalTolerance); } } @@ -893,7 +909,7 @@ TEST(NonlinearMitc4Inertia, for (std::size_t i = 0; i < 24; ++i) { for (std::size_t j = 0; j < 24; ++j) { EXPECT_NEAR(initial.Value().mass(i, j), - eigen.Value().lumped_global24(i, j), 1.0e-12); + eigen.Value().lumped_global24(i, j), kNumericalTolerance); if (i / 6U != j / 6U) EXPECT_EQ(initial.Value().mass(i, j), 0.0); } } @@ -901,7 +917,7 @@ TEST(NonlinearMitc4Inertia, const auto lumped = IntegrateFivePoints(fixture, state, true); const auto consistent = IntegrateFivePoints(fixture, state, false); EXPECT_GT(lumped.energy - consistent.energy, 0.0); - EXPECT_NEAR(lumped.total_mass, consistent.total_mass, 1.0e-12); + EXPECT_NEAR(lumped.total_mass, consistent.total_mass, kNumericalTolerance); } } @@ -940,16 +956,16 @@ TEST(DynamicRotation, ThirdDerivativesMatchAllPermutationsAndInteriorSweep) { const V3 oracle = BlockExponentialDerivative( states[state_index], initial, {directions[p], directions[r], directions[t]}); - EXPECT_LE((analytic - oracle).Norm(), 1.0e-12); + EXPECT_LE((analytic - oracle).Norm(), kNumericalTolerance); EXPECT_LE((analytic - Third(value.Value(), directions[t], directions[p], directions[r])) .Norm(), - 1.0e-12); + kNumericalTolerance); const double error = (analytic - fd).Norm(); const double scale = analytic.Norm() > 0.0 ? analytic.Norm() + fd.Norm() : 1.0; const double metric = error / scale; - const double tolerance = analytic.Norm() > 0.0 ? 1.0e-6 : 1.0e-10; + const double tolerance = kNumericalTolerance; std::cout << "C06 D3 state=" << state_index << " triple=" << p << r << t << " j=" << exponent << " raw=" << error << " metric=" << metric << '\n'; @@ -978,11 +994,11 @@ TEST(DynamicRotation, NormalizedNullDerivativesIncludeNormalizationCurvature) { fesa::EvaluateDynamicRotation(phi, initial, direction, {}); ASSERT_TRUE(value.HasValue()); const V3 first = Multiply(value.Value().null_gradient, direction); - EXPECT_NEAR(value.Value().drill_null.Norm(), 1.0, 1.0e-12); - EXPECT_NEAR(value.Value().drill_null.Dot(first), 0.0, 1.0e-12); + EXPECT_NEAR(value.Value().drill_null.Norm(), 1.0, kNumericalTolerance); + EXPECT_NEAR(value.Value().drill_null.Dot(first), 0.0, kNumericalTolerance); EXPECT_NEAR( value.Value().drill_null.Dot(value.Value().null_directional_second), - -first.Dot(first), 1.0e-12); + -first.Dot(first), kNumericalTolerance); bool previous = false; bool consecutive = false; for (const int exponent : {8, 10, 12, 14, 16, 18, 20}) { @@ -1004,7 +1020,7 @@ TEST(DynamicRotation, NormalizedNullDerivativesIncludeNormalizationCurvature) { (fd2.Norm() + value.Value().null_directional_second.Norm()); std::cout << "C06 null j=" << exponent << " D1=" << error1 << " D2=" << error2 << '\n'; - const bool pass = error1 <= 1.0e-6 && error2 <= 1.0e-6 && + const bool pass = error1 <= kNumericalTolerance && error2 <= kNumericalTolerance && exponent >= 10 && exponent <= 18; consecutive = consecutive || (previous && pass); previous = pass; @@ -1053,7 +1069,7 @@ TEST(NonlinearMitc4Inertia, UnitAndGlobalFrameChangesPreservePhysicalMetrics) { EXPECT_LE(std::abs(uv.Value().kinetic_energy / energy_scale - value.Value().kinetic_energy) / value.Value().kinetic_energy, - 1.0e-12); + kNumericalTolerance); const auto original_scaled = ScaledMatrix(value.Value().mass, fixture.Length()); const auto converted_scaled = ScaledMatrix(uv.Value().mass, units.Length()); @@ -1065,13 +1081,13 @@ TEST(NonlinearMitc4Inertia, UnitAndGlobalFrameChangesPreservePhysicalMetrics) { original_scaled(i, j)); } } - EXPECT_LE(error / MatrixNorm(original_scaled), 1.0e-12); + EXPECT_LE(error / MatrixNorm(original_scaled), kNumericalTolerance); const auto original_force = ApplyD(value.Value().force, fixture.Length()); auto converted_force = ApplyD(uv.Value().force, units.Length()); converted_force.Scale(1.0 / energy_scale); EXPECT_LE(Difference(original_force, converted_force, 1.0).Norm() / original_force.Norm(), - 1.0e-10); + kNumericalTolerance); const auto qx = fesa::ExpSkew(V3{0.4, 0, 0}); const auto qz = fesa::ExpSkew(V3{0, 0, 0.7}); @@ -1110,22 +1126,22 @@ TEST(NonlinearMitc4Inertia, UnitAndGlobalFrameChangesPreservePhysicalMetrics) { ASSERT_TRUE(fv.HasValue()); EXPECT_LE(std::abs(fv.Value().kinetic_energy - value.Value().kinetic_energy) / value.Value().kinetic_energy, - 1.0e-12); + kNumericalTolerance); EXPECT_LE((fv.Value().momentum - rotate(value.Value().momentum)).Norm() / value.Value().momentum.Norm(), - 1.0e-12); + kNumericalTolerance); EXPECT_LE( (fv.Value().angular_momentum - rotate(value.Value().angular_momentum)) .Norm() / value.Value().angular_momentum.Norm(), - 1.0e-12); + kNumericalTolerance); for (std::size_t base = 0; base < 24; base += 3) { const V3 expected = rotate(V3{value.Value().force[base], value.Value().force[base + 1U], value.Value().force[base + 2U]}); const V3 actual{fv.Value().force[base], fv.Value().force[base + 1U], fv.Value().force[base + 2U]}; - EXPECT_LE((actual - expected).Norm() / expected.Norm(), 1.0e-10); + EXPECT_LE((actual - expected).Norm() / expected.Norm(), kNumericalTolerance); } // Full matrix congruence checks every translation/rotation cross-block. double covariance_error = 0.0; @@ -1151,7 +1167,7 @@ TEST(NonlinearMitc4Inertia, UnitAndGlobalFrameChangesPreservePhysicalMetrics) { } } } - EXPECT_LE(covariance_error / MatrixNorm(value.Value().mass), 1.0e-12); + EXPECT_LE(covariance_error / MatrixNorm(value.Value().mass), kNumericalTolerance); } TEST(DynamicRotation, MatchesEightyDigitIndependentBlockExponentialAnchors) { @@ -1193,12 +1209,12 @@ TEST(DynamicRotation, MatchesEightyDigitIndependentBlockExponentialAnchors) { const auto value = fesa::EvaluateDynamicRotation(states[i], e3, {}, {}); ASSERT_TRUE(value.HasValue()); EXPECT_LE((Third(value.Value(), e1, e1, e1) - expected[i][0]).Norm(), - 1.0e-12); + kNumericalTolerance); EXPECT_LE((Third(value.Value(), e1, e2, e3) - expected[i][1]).Norm(), - 1.0e-12); + kNumericalTolerance); EXPECT_LE( (Third(value.Value(), mixed, mixed, mixed) - expected[i][2]).Norm(), - 1.0e-12); + kNumericalTolerance); } } diff --git a/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp index 331abee..5256a2b 100644 --- a/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp +++ b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp @@ -16,6 +16,9 @@ namespace { +// User-approved physical/numerical comparison tolerance for the local experiment. +constexpr double kNumericalTolerance = 1.0e-3; + using Coordinates = std::array; fesa::Node Node(const std::int64_t label, const Coordinates coordinates) { @@ -188,10 +191,10 @@ fesa::Matrix WorkConjugateTangent(const fesa::Matrix& tangent, void ExpectRelative(const double actual, const double expected) { if (expected == 0.0) { - EXPECT_LE(std::abs(actual), 1.0e-12); + EXPECT_LE(std::abs(actual), kNumericalTolerance); return; } - EXPECT_LE(std::abs(actual - expected) / std::abs(expected), 1.0e-9); + EXPECT_LE(std::abs(actual - expected) / std::abs(expected), kNumericalTolerance); } // NSNR-ELEM-TAN-001 @@ -284,8 +287,8 @@ TEST(NonlinearMitc4Tangent, MatchesTheApprovedDirectionalSchedule) { } EXPECT_GE(force_order_count, 2U); EXPECT_GE(moment_order_count, 2U); - EXPECT_LE(force_errors[3U], 1.0e-6); - EXPECT_LE(moment_errors[3U], 1.0e-6); + EXPECT_LE(force_errors[3U], kNumericalTolerance); + EXPECT_LE(moment_errors[3U], kNumericalTolerance); } // NSNR-ELEM-TAN-004 and NSNR-ELEM-TAN-005 @@ -304,7 +307,7 @@ TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) { WorkConjugate(equilibrium.Value().internal_tangent.Multiply( ToVector(rigid_translation)), length); - EXPECT_LE(null_action.Norm(), 1.0e-10); + EXPECT_LE(null_action.Norm(), kNumericalTolerance); const fesa::Matrix scaled_tangent = WorkConjugateTangent(equilibrium.Value().internal_tangent, length); @@ -316,7 +319,7 @@ TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) { } ASSERT_GT(FrobeniusNorm(scaled_tangent), 0.0); EXPECT_LE(FrobeniusNorm(antisymmetric) / FrobeniusNorm(scaled_tangent), - 1.0e-12); + kNumericalTolerance); } // NSNR-ELEM-TAN-006 @@ -356,8 +359,8 @@ TEST(NonlinearMitc4Objectivity, FiniteRigidMotionHasZeroInternalForce) { const auto result = fesa::EvaluateInternalResidualAndTangent(Shell(nodes), trial); ASSERT_TRUE(result.HasValue()); - EXPECT_LE(result.Value().internal_residual.Norm(), 1.0e-10); - EXPECT_LE(result.Value().physical_strain_energy, 1.0e-12); + EXPECT_LE(result.Value().internal_residual.Norm(), kNumericalTolerance); + EXPECT_LE(result.Value().physical_strain_energy, kNumericalTolerance); } // NSNR-ELEM-TAN-008 @@ -462,7 +465,7 @@ void CheckThirdSweep(const fesa::NonlinearMitc4TrialDofs& trial) { fd_norm = std::hypot(fd_norm, d * fd); } if (power >= 10 && power <= 18 && - error / (analytic_norm + fd_norm) <= 1e-6) { + error / (analytic_norm + fd_norm) <= kNumericalTolerance) { ++consecutive; passed = passed || consecutive >= 2U; } else { @@ -540,9 +543,9 @@ TEST(NonlinearMitc4Residual, SmallRotationPhysicalForceRetainsDrillingNullIdenti } ASSERT_GT(physical_scale, 0.0); ASSERT_GT(total_scale, 0.0); - EXPECT_LE(std::abs(physical) / physical_scale, 1.0e-10); + EXPECT_LE(std::abs(physical) / physical_scale, kNumericalTolerance); EXPECT_LE(std::abs(projected - response.Value().drilling_constraint[node]) / total_scale, - 1.0e-10); + kNumericalTolerance); } } @@ -580,7 +583,7 @@ TEST(NonlinearMitc4Residual, for (std::size_t c = 0U; c < 3U; ++c) psi += rotation.Value().drill_null[c] * exact.Value().internal_force[6U * node + 3U + c]; - EXPECT_NEAR(exact.Value().drilling_constraint[node], psi, 1e-11); + EXPECT_NEAR(exact.Value().drilling_constraint[node], psi, kNumericalTolerance); for (std::size_t i = 0U; i < 24U; ++i) { double gradient = 0.; for (std::size_t c = 0U; c < 3U; ++c) { @@ -591,7 +594,7 @@ TEST(NonlinearMitc4Residual, exact.Value().internal_force[6U * node + 3U + c]; } EXPECT_NEAR(exact.Value().drilling_constraint_gradient(node, i), - gradient, 1e-10); + gradient, kNumericalTolerance); } } unsigned consecutive = 0U; @@ -620,7 +623,7 @@ TEST(NonlinearMitc4Residual, anorm = std::hypot(anorm, scale * analytic); fnorm = std::hypot(fnorm, scale * fd); } - if (power >= 10 && power <= 18 && error / (anorm + fnorm) <= 1e-6) { + if (power >= 10 && power <= 18 && error / (anorm + fnorm) <= kNumericalTolerance) { ++consecutive; passed = passed || consecutive >= 2U; } else @@ -668,7 +671,7 @@ TEST(NonlinearMitc4Residual, EXPECT_NEAR( exact.Value().drilling_constraint_third_directional_gradient[node] [row], - finite_difference, 2.0e-5); + finite_difference, kNumericalTolerance); } } } @@ -698,23 +701,23 @@ TEST(NonlinearMitc4Residual, ASSERT_TRUE(minus.HasValue()); EXPECT_NEAR(exact.Value().directional_first, (plus.Value().value - minus.Value().value) / (2.0 * h), - 2.0e-8); + kNumericalTolerance); EXPECT_NEAR(exact.Value().directional_second, (plus.Value().value - 2.0 * exact.Value().value + minus.Value().value) / (h * h), - 3.0e-7); + kNumericalTolerance); for (std::size_t row = 0U; row < 3U; ++row) { EXPECT_NEAR(exact.Value().gradient_directional_first[row], (plus.Value().gradient[row] - minus.Value().gradient[row]) / (2.0 * h), - 3.0e-8); + kNumericalTolerance); EXPECT_NEAR( exact.Value().gradient_directional_second[row], (plus.Value().gradient[row] - 2.0 * exact.Value().gradient[row] + minus.Value().gradient[row]) / (h * h), - 4.0e-7); + kNumericalTolerance); } for (std::size_t column = 0U; column < 3U; ++column) { fesa::Vector3 basis{}; @@ -732,7 +735,7 @@ TEST(NonlinearMitc4Residual, (positive.Value().gradient[row] - negative.Value().gradient[row]) / (2.0 * h), - 2.0e-7); + kNumericalTolerance); } } } @@ -819,17 +822,17 @@ TEST(NonlinearMitc4Residual, TinyRotationPreservesPhysicalForceAndRecoveryEnergy } ASSERT_GT(scale, 0.0); ASSERT_GT(quadratic_energy, 0.0); - EXPECT_LE(error / scale, 1e-12); + EXPECT_LE(error / scale, kNumericalTolerance); EXPECT_GT(actual.Value().physical_strain_energy, 0.0); EXPECT_LE(std::abs(actual.Value().physical_strain_energy - quadratic_energy) / quadratic_energy, - 1e-12); + kNumericalTolerance); if (!warped) { // Uniform director tilt has .5 * (5/6) * G * t * area * theta^2. constexpr double shear_energy = 1.6e-33; EXPECT_LE(std::abs(actual.Value().physical_strain_energy - shear_energy) / shear_energy, - 1e-12); + kNumericalTolerance); } } }