diff --git a/docs/hht-implicit-dynamics/ADR.md b/docs/hht-implicit-dynamics/ADR.md new file mode 100644 index 0000000..4a22476 --- /dev/null +++ b/docs/hht-implicit-dynamics/ADR.md @@ -0,0 +1,80 @@ +# HHT Implicit Dynamics 결정 기록 + +기준: [프로젝트 지침](../../AGENTS.md), [프로젝트 PRD](../PRD.md), [공통 ADR](../ADR.md#adr-028). +Feature 범위·상태는 [PRD](PRD.md), 상세 acceptance intent는 [requirements](requirements.md)를 따른다. + +## HHT-ADR-001 — Requirement와 Research 작업 범위 + +- 상태: `accepted-for-current-task` +- 날짜: 2026-09-22 +- Owner: 메인 assistant +- 승인 근거: 사용자가 HHT implicit dynamics의 requirement/research gate를 요청하고 `reference/shellNL-dyn-meshfine`을 지정했다. +- 결정: `dev` commit `23b37ca644563b51472654af3b22811f077f15f6`에서 분리한 `hht-implicit-dynamics` worktree에서 요구사항·연구 문서를 작성한다. 원본 reference와 solver source는 이번 작업에서 수정하지 않는다. +- 이유·영향: 기존 논문 재현 실험과 eigen 질량 진단의 변경을 현행 dev나 새로운 동적 물리 계약으로 오인하지 않도록 기준을 구별한다. 새로운 dynamics 구현은 이 두 gate 후의 별도 작업이다. + +## HHT-ADR-002 — Reference 입력과 비교 정책 보존 + +- 상태: `accepted-existing-contract` +- Owner: 메인 assistant +- 승인 근거: 사용자 지정 reference 및 프로젝트 AGENTS/TOLERANCE. +- 결정: 원본 `NLGEOM=YES`, `DIRECT`, 시간증분 `0.1`, 종료시간 `10`을 보존하고 101 frames / 9 nodes 전체 U/UR 이력을 대상으로 한다. U blocking, UR warning-only와 공통 frame별 수치 정책을 적용한다. Eigen 첫 구현의 RMS 5% 예외는 적용하지 않는다. +- 이유·영향: 동적 이력의 일부만 선택하거나 reference를 변경하여 일치를 만드는 것을 방지한다. CSV가 제공하지 않는 속도·가속도·반력·에너지에는 독립 검증 evidence가 필요하다. 최종 family/schema/mapping은 후속 owner가 정의한다. + +## HHT-ADR-003 — 물리·입출력 결정의 경계 + +- 상태: `partially-resolved`; 질량 종류는 HHT-ADR-007로 확정, 나머지는 PRD의 현재 상태 참조. +- Owner: 메인 assistant; 근거는 Research, 정식화는 Formulation, 독립 판단은 Numerical Review. +- 결정 대상: `requirements.md`의 `HHT-OPEN-001–006`. +- 이유: 현재 dev의 small-rotation mass kernel과 nonlinear internal-force evaluator만으로 비선형 동적 관성의 일관성이 성립하지 않는다. HHT의 coefficient convention, 초기 가속도와 load onset, drilling 영질량 처리, 회전 운동학 및 공식 동적 출력 의미가 함께 연결되어야 한다. +- 영향: 이번 gate에서 연구 후보식을 implementation contract로 승인하지 않는다. 공개 문헌으로 확인 가능한 입력 기본값은 근거를 기록하여 해소하되, 질량·관성의 물리 범위 및 영속 결과 schema는 후속 검토와 필요한 사용자 결정을 거친다. 비교 결과에 맞춘 parameter/mass 조정은 결정 근거로 사용하지 않는다. + +## HHT-ADR-004 — 확인된 공개 입력 기본값 + +- 상태: `researched-input-interpretation`; FESA 정식화·전체 입력 지원 승인과 구별. +- 날짜: 2026-09-22 +- Owner: 메인 assistant; 근거 수집은 research-agent. +- 근거: [Research §3 및 S4–S7](research.md)의 공식 Abaqus 2025 공개 문서와 원본 INP 대조. +- 해석: 무접촉 기본 application은 TRANSIENT FIDELITY/HHT-TF이고 `alpha=-0.05`, `beta=0.275625`, `gamma=0.55`다. 명시 amplitude가 없는 집중하중은 STEP이며 첫 증분부터 전체 크기다. DIRECT는 이 모델에서 고정 시간간격이고, 기본 INITIAL 의미는 초기 가속도 계산이다. +- 영향: 이 공개 의미는 CSV fitting으로 정하지 않으며 기본값을 정하기 위한 반복 사용자 질문도 필요하지 않다. FESA 최초 old-force history, `t0-/t0+`와 frame 0, mass-null 초기화는 명시적으로 정식화해야 한다. 공개 의미만으로 reference를 생성한 solver의 내부 초기화·관성 알고리즘을 확인했다고 주장하지 않는다. + +## HHT-ADR-005 — Formulation 작업과 독립 검토 + +- 상태: `accepted-for-current-task` +- Owner: 메인 assistant +- 승인 근거: 사용자의 후속 formulation gate 진행 요청 및 프로젝트의 수식 변경 독립 검토 원칙. +- 결정: 기존 Feature worktree에서 시간적분·관성·초기화·잔차·접선·회복량을 정식화하고 독립 검토한다. 유한회전과 drilling 영질량의 영향을 검산하며 질량 종류 등 물리 선택은 명시적인 제안으로 제시한다. +- 영향: 기존 requirement/research 산출물은 해당 작업 시점의 기록으로 유지한다. 독립 수식 검산은 허용하지만 solver 구현·reference 수정·Abaqus 실행·영속 schema 확정은 포함하지 않는다. + +## HHT-ADR-006 — 물리적 가속도와 보조 가속도를 분리하는 HHT 후보 + +- 상태: `recommended-after-formulation-review`; 물리·구현 범위의 사용자 승인과 구별. +- Owner: 메인 assistant; 수식 owner는 formulation-agent, 독립 검토는 numerical-review-agent. +- 근거: [정식화](formulation.md), [관성 상세](formulation-inertia.md), [독립 검토](numerical-review.md), [Brüls–Arnold (2008) 원문 식 (5)–(8)](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf). +- 추천: `a`는 물리적 좌표 가속도, `A`는 Newmark에 쓰는 보조량으로 정의하고 `A_(n+1)=(1+alpha)a_(n+1)-alpha a_n`와 endpoint 동적 평형을 함께 푼다. `M(q)` 및 convective inertia의 정확한 미분을 접선에 포함한다. +- 이유: 단순 current-mass 대입은 variable-mass 제조해에서 1차로 떨어졌다. 두 시점 질량의 가중은 회전하는 질량 영공간을 소실시킬 수 있다. 추천 후보는 제한된 variable-mass 및 moving-null 사례에서 2차 수렴을 보이고 일정 질량 HHT로 환원한다. +- 영향: `A`를 physical acceleration으로 출력하지 않는다. 초기화는 전값 STEP force의 `0+` 동적 평형과 대수 일관성을 함께 만족하도록 제안한다. 원문 nonsingular-mass 수렴 결과를 실제 shell의 drilling DAE에 자동 적용하지 않으며 구현 후 별도 시간 수렴·관성·반력·에너지 검증이 필요하다. + +## HHT-ADR-007 — HHT 일관질량 선택 + +- 상태: `accepted-user-decision` +- Owner: 메인 assistant; 수학 정의는 formulation-agent. +- 승인 근거: 사용자가 질량 선택 질문에 “일관질량 (추천): shell 속도 보간의 운동에너지를 그대로 적분”을 선택했다. +- 결정: [관성 상세 §2](formulation-inertia.md)의 **consistent kinetic mass**를 HHT 기본 및 이번 구현 범위로 사용한다. Scalar-moment nodal lumped 후보는 비교 연구로 보존하며 이번 HHT 지원 요구에서 제외한다. +- 이유: Consistent는 선택한 shell 속도 보간의 운동에너지를 그대로 적분한다. Lumped는 총질량과 scalar thickness moment를 유지하지만 일반 속도장·orbital inertia의 운동에너지가 달라진다. 이전 eigen reference와의 근접성은 채택 근거로 쓰지 않는다. +- 영향: 선택한 유한회전 운동학에서 관성력·접선을 계산하며 drilling에 임의 질량을 주지 않는다. 기존 정적·eigen mass 계약과 thesis 실험 branch에는 영향을 주지 않는다. 이 결정만으로 다른 HHT 후보식·공개 parameter·영속 결과 schema까지 승인되었다고 확대하지 않는다. + +## HHT-ADR-008 — Numerical Review와 추가 수치 연구 + +- 상태: `accepted-for-current-task`; 수치 검증 설계와 실제 solver 통과를 구별. +- Owner: 메인 assistant; 독립 수치 기준·검증 설계는 numerical-review-agent, 추가 근거는 research-agent. +- 승인 근거: 사용자의 “numerical research 단계도 진행하세요” 요청과 프로젝트의 독립 수치 검토 원칙. +- 결정: [수치 기준](numerical-controls.md), [검증 사례](verification-cases.md), [추가 연구](numerical-research.md), [reference 비교 계약](reference-model.md)을 작성한다. User-selected consistent mass에 대해 Newton·초기 대수조건·질량 영공간·conditioning·관성 접선·시간수렴의 반증 가능한 기준을 정의한다. +- 이유·영향: 외부 reference 공통 tolerance를 유지하면서 내부 계산 정확도와 별도의 물리 coverage를 검증한다. 제한된 수식 probe의 실행 결과를 전체 shell·구현·reference 검증으로 승계하지 않는다. I/O field binding과 실제 구현·build/test·reference 비교는 후속 작업이다. + +## HHT-ADR-009 — I/O·구현·Physics Evaluation 수행 + +- 상태: `accepted-for-implementation`; 검증 결과의 통과와 구별. +- Owner: 메인 assistant; I/O 정의·구현·독립 physics review는 각 할당 담당자. +- 승인 근거: 사용자의 “I/O definition, implementation, physics evaluation gate 모두 진행하세요” 및 후속 “깃에 커밋 하면서 진행하세요”. +- 결정: HHT-ADR-006의 검토된 physical a/auxiliary A 분리, 0+ full STEP 초기화, HHT-ADR-007의 consistent mass를 실제 구현 대상으로 채택한다. 공개 입력은 [I/O 계약](io-contract.md)의 단일 NLGEOM shell DIRECT·기본 계수·zero 초기조건/고정 BC·constant force CLOAD subset이다. 별도 dynamic state와 dynamic-only HDF5 schema v1을 추가하고 기존 static/eigen 계약은 보존한다. +- 영향: [구현 계획](implementation-plan.md)에 따라 MSVC Debug build/CTest, 독립 수치·물리 검증과 지정 원본 reference 비교를 수행한다. 확인된 문서·검증된 구현 묶음을 현 branch에서 순차 커밋한다. 외부 solver 실행·reference 보정·push/release는 포함하지 않는다. 실패·미실행 검증을 통과로 표시하지 않는다. diff --git a/docs/hht-implicit-dynamics/FEATURED_AGENTS.md b/docs/hht-implicit-dynamics/FEATURED_AGENTS.md new file mode 100644 index 0000000..a6efd0e --- /dev/null +++ b/docs/hht-implicit-dynamics/FEATURED_AGENTS.md @@ -0,0 +1,40 @@ +# HHT Implicit Dynamics 작업 지침 + +- Feature ID: `hht-implicit-dynamics` +- Owner: 메인 assistant +- 목적: 지정 shell reference를 대상으로 HHT implicit 동적 해석의 요구사항·정식화·수치 검증 계약을 마련한다. +- 현재 작업 범위: I/O definition·implementation·physics evaluation 전체 수행 및 검증 묶음별 Git commit. +- 현재 상태와 유일한 checklist: [PRD](PRD.md#현재-checklist). + +기준: [프로젝트 지침](../../AGENTS.md), [프로젝트 PRD](../PRD.md), +[공통 ADR](../ADR.md#adr-028), [ARCHITECTURE](../ARCHITECTURE.md), +[TOLERANCE](../TOLERANCE.md), [Feature 문서 관리](../FEATURE_DOCUMENTATION.md). + +## 지역 지침 + +- 기준 코드는 `dev` commit `23b37ca644563b51472654af3b22811f077f15f6`다. `hht-implicit-dynamics` worktree에서 구현·검증한다. +- 사용자가 지정한 `reference/shellNL-dyn-meshfine`의 INP와 CSV를 그대로 조사한다. 원본 reference는 보존하고 FESA solver 구현·실행·비교를 수행한다. Reference 재생성·외부 solver 실행은 포함하지 않는다. +- 요구사항은 관찰 가능한 동작과 검증 의도를, 연구는 일차 출처의 사실·계산·추론·남은 공백을 소유한다. 연구식을 승인된 구현 정식화로 취급하지 않는다. +- 이전 논문 재현 실험의 drilling·질량·접선 변경과 dev 일관질량 eigen 진단은 이 기능의 물리 선택이나 구현 승인으로 이어지지 않는다. +- 사용자가 [HHT-ADR-007](ADR.md)에서 consistent kinetic mass를 선택했다. 이번 HHT 구현 범위는 일관질량이며 lumped는 비교 연구·기존 eigen 회귀의 경계에서만 다룬다. +- `NLGEOM=YES`, `DIRECT` 시간 의미와 전체 reference 이력을 보존한다. Abaqus 입력 subset 지원은 내부 요소·시간적분 구현 동등성 약속이 아니다. +- 공통 tolerance와 shell U blocking / UR warning-only 분류를 유지한다. Eigen 전용 RMS 5% 예외를 동적 해석에 확대하지 않는다. +- 미결정 물리 범위와 근거 부족을 명시하며 결과에 맞추어 alpha·질량·하중 이력을 변경하지 않는다. 공개 문헌으로 확정되는 입력 의미에는 불필요한 재승인 절차를 추가하지 않는다. + +## 문서 소유권 + +[구현 계획](implementation-plan.md)은 메인 소유다. [I/O 계약](io-contract.md)은 io-definition-agent, 코드·시험은 메인이 지정한 input/inertia/state-output/procedure 소유권을 따른다. 최종 물리 평가는 독립 physics-evaluation-agent가 맡는다. + +| 문서 | Owner | 책임 | +| --- | --- | --- | +| [FEATURED_AGENTS.md](FEATURED_AGENTS.md) | 메인 assistant | 지역 지침과 할당 | +| [PRD.md](PRD.md) | 메인 assistant | 요청·범위·현재 checklist·현행 코드와의 차이 | +| [ADR.md](ADR.md) | 메인 assistant | 결정·승인 근거와 미결정 경계 | +| [requirements.md](requirements.md) | requirement-agent | 원자적 요구사항·acceptance intent·reference inventory | +| [research.md](research.md), [numerical-research.md](numerical-research.md) | research-agent | HHT·공개 입력 의미·수치 근거·검증 후보·근거 공백 | +| [formulation.md](formulation.md), [formulation-inertia.md](formulation-inertia.md) | formulation-agent | HHT 시간적분·초기화·shell 관성·잔차와 접선·검증 불변식 | +| [numerical-review.md](numerical-review.md), [numerical-review-probes.md](numerical-review-probes.md) | numerical-review-agent | 정식화와 수치 계약의 독립 검토·반례·재현 가능한 수식 검산 | +| [numerical-controls.md](numerical-controls.md), [verification-cases.md](verification-cases.md) | numerical-review-agent | 내부 수치 기준과 반증 가능한 검증 사례 | +| [reference-model.md](reference-model.md) | 메인 assistant; numerical-review-agent 검토 | 원본 inventory·시간·identity·수치 비교 계약; I/O binding은 후속 | + +담당자는 할당 문서만 수정하고 peer를 재위임하지 않는다. 현재 I/O·구현·solver 물리 검증을 수행하며 빈 보고서는 생성하지 않는다. 수식 검산은 solver build/test나 지정 reference의 동적 해석 통과 evidence와 구별한다. diff --git a/docs/hht-implicit-dynamics/PRD.md b/docs/hht-implicit-dynamics/PRD.md new file mode 100644 index 0000000..8223c0c --- /dev/null +++ b/docs/hht-implicit-dynamics/PRD.md @@ -0,0 +1,94 @@ +# HHT Implicit Dynamics PRD + +- Feature ID: `hht-implicit-dynamics` +- 기준: `dev` commit `23b37ca644563b51472654af3b22811f077f15f6` +- 문서 상태: `io-contract-defined; implementation-in-progress; physics-pending`, 2026-09-22 +- 기준 문서: [프로젝트 지침](../../AGENTS.md), [프로젝트 PRD](../PRD.md), [공통 ADR](../ADR.md#adr-028). + +## 사용자 요청과 범위 + +사용자는 HHT implicit 적분법을 사용한 dynamic 해석 구현을 준비하기 위해 requirement/research, formulation, numerical research를 차례로 요청했고, 검증 모델로 `reference/shellNL-dyn-meshfine`을 지정했다. 현재 산출물은 요구사항·근거·정식화·독립 수치 검토와 검증 계약이다. 사용자의 후속 I/O definition·implementation·physics evaluation 전체 진행 요청에 따라 [I/O 계약](io-contract.md)과 [구현 계획](implementation-plan.md)을 작성하고 실제 solver 구현·검증을 진행한다. + +[시간적분 정식화](formulation.md)와 [shell 관성 상세](formulation-inertia.md)의 후속으로 [독립 수치 검토](numerical-review.md), [수치 기준](numerical-controls.md), [검증 사례](verification-cases.md), [추가 연구](numerical-research.md)를 작성한다. 사용자가 **consistent kinetic mass**를 선택했으며 이번 HHT 구현 범위는 일관질량이다. 후속 전체 gate 진행 요청에 따른 구현 대상·공개 subset·출력 schema는 [HHT-ADR-009](ADR.md)에 기록한다. + +최소 구현 목표 후보는 기존 FESA-MITC4 shell의 선형 탄성 재료·기하비선형 동적 해석이다. 단일 step의 `*Dynamic,direct`, 초기 정지 상태, 고정 외곽 경계와 집중하중에 필요한 지원을 우선 정의한다. 질량·회전 관성·HHT 생략값·공식 동적 결과 범위는 근거를 확인하여 후속 정식화와 연결한다. 임의 amplitude, 시간변화 경계, 물리적 damping, contact, 소성, explicit, 다중 step/restart 또는 다른 요소의 dynamics는 이번 요청에서 지원을 약속하지 않는다. + +## 지정 reference와 acceptance intent + +원본: [INP](../../reference/shellNL-dyn-meshfine/shellNL-dyn-meshfine.inp), [변위 CSV](<../../reference/shellNL-dyn-meshfine/shellNL-dyn-meshfine displacements.csv>). 상세 inventory와 요구사항 ID는 [requirements](requirements.md)가 소유한다. + +- 평면 shell 729 nodes / 676 S4 elements, 두께 `0.5`, `E=2.1e11`, `nu=0.3`, `rho=2700`. +- `NLGEOM=YES`, `DIRECT`, 시간증분 `0.1`, 종료시간 `10`, 중심 node `2`에 global 1–3방향 각각 `-100000` 집중하중. +- CSV는 101 frames × 9 nodes = 909 rows, U/UR 6성분을 제공한다. 초기 frame을 포함한 모든 시간 이력이 비교 대상이다. +- source instance/node와 increment/time을 결정적으로 대응시킨다. HDF5 전체 node 결과에서 승인된 reference subset으로 projection한 뒤 누락·추가·중복·nonfinite를 검사한다. 최근접 시간 선택, 시간 보간 또는 평균으로 대응 오류를 숨기지 않는다. +- [공통 TOLERANCE](../TOLERANCE.md)의 frame별 U blocking / UR warning-only family를 사용한다. [Reference 비교 계약](reference-model.md)은 101 frames × 2 quantities의 202 families와 각 27 scalars, source/time mapping을 정의한다. Eigen 전용 예외는 적용하지 않는다. 실제 HDF5 field binding은 후속 I/O 계약이다. +- V/A, 동적 반력, 에너지 reference는 없다. 독립 oscillator·시간 수렴·관성 평형·회전 및 에너지 검증이 필요하다. 원본 모델의 작은 변위만으로 유한회전 동역학을 검증했다고 주장하지 않는다. + +## 현행 dev와의 차이 + +다음은 소스 정적 조사 결과이며, 현재 dynamic 실행 결과가 아니다. + +| 경계 | 현재 코드 evidence | 후속 작업에서 정의할 내용 | +| --- | --- | --- | +| 입력·semantic model | [domain_mapping_context.cpp](../../src/fesa/io/abaqus/domain_mapping_context.cpp)의 `ParseStepBlock` 분기는 STATIC/FREQUENCY를 처리한다. `FinalizeDensities`는 FREQUENCY 외 DENSITY를 거부한다. [Domain](../../include/fesa/model/domain.h)의 step 정의도 동적 절차가 없다. | DYNAMIC/DIRECT 및 density의 승인 subset, 실제 시간·기본값·초기조건·진단 | +| 절차 선택 | [prepared_analysis](../../src/fesa/analysis/prepared_analysis.cpp)는 선형 정적·비선형 정적·고유치만 선택한다. | 정적 lifecycle에 누적하지 않는 별도 dynamic procedure와 결과 경계 | +| 상태 | [NonlinearAnalysisState](../../include/fesa/analysis/nonlinear_analysis_state.h)는 load factor, q, triad/director, 내외력과 정적 recovery를 저장한다. v/a 및 실제 시간 이력은 없다. | trial/committed q/v/a·시간·이력의 소유권과 일괄 commit/실패 보존 | +| 내력·접선 | [NonlinearMitc4Shell](../../include/fesa/elements/nonlinear_mitc4_shell.h)은 total rotation-vector trial의 내력·접선·recovery 평가를 제공한다. | 이 좌표계에 맞는 관성력·HHT residual·유효 접선과 일관된 회전 운동학 | +| 질량 | [Mitc4Mass](../../include/fesa/elements/mitc4_shell.h)는 consistent/lumped global24를 제공한다. [현행 frequency](../../src/fesa/analysis/linear_frequency_procedure.cpp)는 lumped nodal mass를 사용한다. | 동적 질량 선택·구성 의존성·회전 관성과 mass-null 자유도 처리; 정적/eigen 계약의 무단 변경 금지 | +| 선형계 backend | [비선형 PARDISO adapter](../../src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp)는 이미 비대칭 full CSR 접선을 처리한다. | 동적 관성 접선의 비대칭 성분 보존과 기존 경계 재사용 검토 | +| 에너지 | [비선형 요소](../../src/fesa/elements/nonlinear_mitc4_shell.cpp)는 physical+drilling potential의 내력·Hessian을 계산하지만 반환 에너지는 physical strain energy다. | 운동에너지·수치 drilling 에너지·외력 일을 구분하는 회복 계약 | +| 공식 결과 | 정적/eigen 전용 state/writer가 존재한다. | HDF5의 step/time/frame 및 q/v/a·반력·에너지 공개 범위, readback·comparison | + +Density 허용과 HHT update 식 추가만으로 기능이 완성되는 것은 아니다. Semantic model부터 동적 상태, 방정식, recovery, 공식 결과 및 검증까지 연결해야 한다. + +## 결정과 후속 인계 + +결정 기록은 [ADR](ADR.md), 결정 항목의 원본은 [requirements의 HHT-OPEN-001–006](requirements.md#범위-경계와-미결정-사항)이다. 질량 종류는 사용자 결정으로 확정했다. 수식과 수치 검증 계약을 연결하되 공개 parameter·결과 범위 및 I/O binding은 후속으로 확정한다. 문서 검토만으로 구현이 검증되었다고 표시하지 않는다. + +이번 [Research](research.md)에서 확인한 공개 입력 의미는 무접촉 기본 HHT-TF 계수 `alpha=-0.05`, `beta=0.275625`, `gamma=0.55`, amplitude 생략 시 STEP 하중, DIRECT 고정 증분 및 기본 초기 가속도 계산이다. 최초 HHT force history와 초기 frame의 결합, FESA 질량·유한회전 관성 및 공개 결과 범위는 여전히 후속 계약이 필요하다. 원 HHT의 표준 결합 계수 범위와 vendor가 허용하는 전체 parameter 범위를 혼동하지 않는다. + +Reference의 `max|U3|/thickness≈9.84e-5`, 기존 동일 mesh eigen reference에 대한 `dt/T1=7.4567`은 작은 응답과 성긴 시간간격을 보여준다. 이 CSV와의 일치만으로 큰 회전 관성, V/A 또는 시간 정확도를 검증할 수 없다. 단위계는 입력의 일관된 단위를 따르며 SI로 단정하지 않는다. + +## Formulation 결과 + +[정식화](formulation.md)의 추천 후보는 유한회전 kinetic energy의 정확한 관성력과 물리적 가속도 `a` / 적분용 보조 가속도 `A`의 분리다. `A_(n+1)=(1+alpha)a_(n+1)-alpha a_n`를 Newmark update에 사용하고 endpoint 물리 평형을 푼다. 일정 질량에서는 기존 HHT 식으로 환원한다. 비선형 내력은 현행 dev의 objective drilling을 포함한 에너지 미분을 보존하며 이전 논문 실험의 접선·drilling 변경을 도입하지 않는다. + +[독립 수식 검토](numerical-review.md)와 [재현 가능한 검산](numerical-review-probes.md)은 일정 질량 HHT 식에 `M(q)`만 대입하면 일반적으로 1차 오차가 남음을 확인했다. 추천 후보는 smooth variable-mass 제조해와 영질량 방향이 회전하는 제한된 모델에서 U/V/물리적 a의 2차 수렴을 보였다. 메인도 공개된 재현 코드를 실행하여 표의 결과를 확인했다. 이는 실제 MITC4 shell의 동적 해석 또는 보편적인 singular-DAE 증명의 대체물이 아니다. + +| 기존 open 항목 | 이번 정식화에서 구체화한 후보 | 남은 결정/검증 | +| --- | --- | --- | +| HHT-OPEN-001 | 음수 alpha 관례, 결합 beta/gamma, 물리적 a와 보조 A 분리 | 추천 HHT 확장 채택 및 공개 parameter subset | +| HHT-OPEN-002 | frame 0을 `0+`로 정의, 전값 force와 동적·대수 초기 평형으로 a0 계산, A0=a0; 수치 기준 N2–N4 정의 | 초기화 계약 채택과 I/O 진단 binding, 구현 후 검증 | +| HHT-OPEN-003 | **사용자 선택: consistent kinetic mass**; lumped는 비교 연구로 보존하고 이번 HHT 지원에서 제외 | 질량 종류 결정 완료; 임의 drilling inertia 없이 회전관성·영공간의 실제 구현 검증 필요 | +| HHT-OPEN-004 | 기존 total rotation-vector chart에서 exact 관성·convective 항·비대칭 접선, moving-null 조건; N4–N7 및 C05–C08 정의 | 구현 후 실제 shell의 regularity·시간 수렴·회전 검증 | +| HHT-OPEN-005 | 물리적 a/각가속도와 A, endpoint 반력, T/Uphys/Udrill/W의 수학적 의미 구별 | 공식 공개 quantity·HDF5 schema 및 mapping의 I/O 결정 | +| HHT-OPEN-006 | 고정 100 increments, full Newton, 실패 trial 폐기, 동일 state의 전체 이력 commit; N2의 100 corrections/수렴 기준 정의 | I/O 실패 lifecycle과 reference 밖 입력 범위, 구현 후 실패 보존 검증 | + +현재 후보에 대한 독립 검토와 후속 구현 acceptance를 구별한다. 실제 shell 검증의 미실행은 현재 수식에서 오류를 발견했다는 뜻이 아니며, 구현 이후 반드시 확보해야 할 evidence다. 물리 선택과 수치·입출력 계약은 구현 착수 전 확정할 항목이다. + +## Numerical Review 결과 + +[내부 수치 기준 N1–N7](numerical-controls.md)은 변위·회전의 scale를 분리하고, 관성·내력·외력의 상쇄 전 크기로 endpoint 잔차를 판정한다. Newton 잔차와 correction/rate 판정은 각각 `1e-10`, 선형 solve backward error는 `1e-12`다. 해석적으로 정의한 질량 영공간과 양의 물리 질량을 구별하며, 국소 positive block·대수 regularity에는 scale된 `rcond_inf>=2^-32`를 요구한다. 이 값은 외부 reference tolerance의 변경이 아니다. + +[검증 사례 C01–C12](verification-cases.md)는 oscillator·STEP 초기화·variable mass·moving null, 운동에너지 적분·관성 접선, 반력·에너지·실제 shell 유한회전, 실패 보존·결정론·전체 reference 이력을 포함한다. D3 회전 미분에는 독립 block-exponential oracle를 지정했다. 모든 신규 solver 사례는 **NOT RUN**이다. + +추가 full Cartesian moving-null 수식 probe는 전체 시점 최대 오차에서 q/v/물리적 a의 최종 refinement 차수 `1.99985/1.99992/1.99834`를 보였고 독립 reviewer가 재현했다. 이는 기하학적 normalization이나 각도 축약 없이 계산한 제한된 수식 evidence다. 전체 N1–N7 구현 검사와 실제 shell·reference 통과를 의미하지 않는다. 자세한 실행 범위는 [수식 검산](numerical-review-probes.md)에 기록한다. + +## 현재 checklist + +- [x] 사용자 지정 모델·출력과 현행 dev 기준 확인. +- [x] 격리 worktree 및 Feature 핵심 문서 작성. +- [x] Requirement 초안: 원자적 요구사항과 acceptance intent 작성. +- [x] Research 초안: 일차·공식 근거의 열람 범위, 입력 생략 의미, 검증 후보와 근거 공백 정리. +- [x] 문서 간 범위·비교 정책·로컬 링크 및 reference 원본 보존 정적 검토. +- [x] 사용자 질량 선택 반영: 이번 HHT 기본·구현 범위는 consistent kinetic mass. +- [x] 나머지 물리·공개 parameter·공식 결과 범위 확정: HHT-ADR-009와 I/O 계약. +- [x] Formulation 후보 작성: HHT·정확한 유한회전 관성·초기화·잔차·접선·회복. +- [x] 독립 수식 검토와 variable-mass / moving-null 제조해 검산; 메인 재현 확인. +- [x] 추가 수치 근거 조사와 full Cartesian moving-null 독립 검산·재현. +- [x] Numerical Review의 내부 기준·검증 사례·reference 수치 비교 계약 검토 완료. +- [x] I/O 계약과 실제 field/schema binding 정의·정식화 대조. +- [ ] HDF5 구현과 readback 검증. +- [ ] Solver 구현·구성된 build/test·독립 물리 및 reference 검증. + +문서 gate의 완료와 solver 기능의 완료를 구별한다. 현재 구현 단계에서는 내부 수치 검증, solver build/test 및 지정 reference 해석·비교를 실행한다. 통과 여부와 physics 평가는 실제 evidence로 갱신한다. diff --git a/docs/hht-implicit-dynamics/formulation-inertia.md b/docs/hht-implicit-dynamics/formulation-inertia.md new file mode 100644 index 0000000..669f99a --- /dev/null +++ b/docs/hht-implicit-dynamics/formulation-inertia.md @@ -0,0 +1,209 @@ +# 유한회전 관성·초기화·HHT 상세 + +[주 정식화](formulation.md)의 구현 수식 상세다. 상태와 미결정 항목은 주 문서를 따른다. 아래 전개는 명시한 FESA 운동학에서의 자체 유도이며 Abaqus 내부 절차의 기술이 아니다. + +사용자 선택에 따라 이번 HHT 기본 및 구현 질량은 consistent `T_C/M_C`다. Lumped 식은 비교용 이론·기존 eigen 회귀 참고로 보존하며 이번 HHT 구현·검증 요구를 추가하지 않는다. 다른 physics·알고리즘·공개 output의 승인은 이 선택과 별개다. + +## 1. 회전과 관성의 미분 + +`P=[phi]_×`, `theta=|phi|`에 대해 + +\[ +J_l=I+\frac{1-\cos\theta}{\theta^2}P ++\frac{\theta-\sin\theta}{\theta^3}P^2,\qquad B=Dd=-[d]_\times J_l. +\] + +0에서 scalar coefficient는 analytic series를 사용한다. `v_phi,a_phi`는 fixed-chart 시간미분이다. 방향 `p,s,r`에 대한 director 미분을 +`B p=Dd[p]`, `C[p,s]=D²d[p,s]`, `D3[p,s,r]=D³d[p,s,r]`라 정의하면 + +\[ +C[p,s]=(DJ_l[s]p)\times d+(J_lp)\times((J_ls)\times d), +\] +\[ +\begin{aligned} +D3[p,s,r]={}&(D^2J_l[s,r]p)\times d+(DJ_l[s]p)\times B r\\ +&+(DJ_l[r]p)\times B s+(J_lp)\times C[s,r]. +\end{aligned} +\] + +방향 교환에 대한 대칭성이 필요하다. 정확한 미분을 scalar coefficient의 비특이 series까지 구현하기 위한 동등한 정의는 + +\[ +D^k\exp_P[H_1,\ldots,H_k] +=\sum_{m=k}^{\infty}\frac1{m!} +\sum_{\sigma\in S_k}\ \sum_{r_0+\cdots+r_k=m-k} +P^{r_0}H_{\sigma(1)}P^{r_1}\cdots H_{\sigma(k)}P^{r_k},\quad k=1,2,3. +\] + +여기서 `H=[direction]_×`이고 오른쪽에 `d0`를 곱하면 필요한 director derivative다. 이는 수렴하는 해석적 정의이며 production에서 무한합을 그대로 계산하라는 뜻이 아니다. Rodrigues/Fréchet/정확한 미분으로 같은 도함수를 계산하고 finite difference는 독립 검증에 사용한다. + +물질점에 대해 `H=D_q x`, `B_x[p,s]=D²_q x[p,s]`, `C_x[p,s,r]=D³_q x[p,s,r]`로 쓰면 + +\[ +H p=\sum_iN_i(p_{ui}+zB_i p_{\phi i}),\quad +B_x[p,s]=\sum_iN_i z C_i[p_{\phi i},s_{\phi i}], +\] +\[ +C_x[p,s,r]=\sum_iN_i zD3_i[p_{\phi i},s_{\phi i},r_{\phi i}],\quad +\dot x=Hv,\quad \ddot x=Ha+B_x[v,v]. +\] + +Consistent 관성 covector와 partial derivative action은 다음으로 직접 조립 가능하다: + +\[ +I_C=\int\rho H^T\{Ha+B_x[v,v]\}\,dV_0, +\quad M_C=\int\rho H^TH\,dV_0, +\] +\[ +\begin{aligned} +D_qI_C[p]=\int\rho\big(&B_x[p,\cdot]^T\{Ha+B_x[v,v]\}\\ +&+H^T\{B_x[p,a]+C_x[p,v,v]\}\big)dV_0,\\ +D_vI_C[p]&=2\int\rho H^TB_x[v,p]\,dV_0,\qquad D_aI_C[p]=M_Cp. +\end{aligned} +\] + +첫 항 `B_x[p,·]^T ẍ`를 빼면 단순 current-mass tangent에 불과하다. `B_x[p,a]`, `C_x[p,v,v]`, `2B_x[v,δv]`도 생략하지 않는다. `T`의 Euler–Lagrange 전개에서 `D H`의 대칭성에 의해 `d(H^T)/dt`와 `∂q T` 항이 상쇄되어 위 식을 얻는다. 원형상 밀도·measure는 시간에 따라 갱신하지 않는다. + +비교용 lumped 식에는 `x_i(z)=X_i+u_i+z d_i`, `H_i=D_q x_i`를 사용하여 위 적분을 `Σ_i ∫rho N_i (·)_i dV_0`로 바꾼다. `N_i`는 q와 무관하므로 동일한 derivative action이 성립한다. 별도의 경험적 gyroscopic 또는 Coriolis 보정항을 추가하지 않는다. + +## 2. 채택 질량·비교식과 적분 + +주 문서의 scalar moment로 + +\[ +(M_C)_{ij}=\begin{bmatrix} +C_{ij}^{(0)}I_3&C_{ij}^{(1)}B_j\\ +C_{ij}^{(1)}B_i^T&C_{ij}^{(2)}B_i^TB_j +\end{bmatrix},\quad +(M_L)_{ii}=\begin{bmatrix}m_iI_3&s_iB_i\\s_iB_i^T&j_iB_i^TB_i\end{bmatrix}, +\quad (M_L)_{ij}=0\ (i\ne j). +\] + +`rho`의 차원은 `mass/length³`, `(C0,C1,C2)` 및 `(m,s,j)`의 차원은 `(mass,mass·length,mass·length²)`다. `M_uu/M_uphi/M_phiphi`는 각각 같은 세 차원, residual translation은 force, rotation은 force·length이다. `v_u,a_u`는 `length/time,length/time²`; `v_phi,a_phi,omega,omega_dot`는 `1/time,1/time²`. 각도 단위는 rad이며 SI를 추정하지 않는다. + +Moments는 현행 [eigen2 mass 계약 M2–M3](../mitc4-shell-feast-eigen2/formulation-mass-stabilization.md#m2-consistent-scalar-moments와-적분)의 `3×3×3` Gauss–Legendre를 기준으로 하며, 그 scalar-moment lumping은 비교식에만 사용한다. 긴 유도는 선행 [mass 계약 M3](../mitc4-shell-feast-eigen/formulation-mass.md#m3-consistent-matrix와-충분한-quadrature)에 있다. 점은 `−sqrt(3/5),0,+sqrt(3/5)`, weight는 `5/9,8/9,5/9`, 순서는 증가하는 `xi→eta→zeta`다. Geometry가 bilinear degenerated shell이고 rho/t가 일정할 때 `Ni Nj z^k J0`의 축별 차수는 최대 4여서 이 moment 적분은 정확하다. `B_i`와 그 도함수는 nodal 값이므로 현재 관성과 tangent도 같은 moment를 사용한다. + +모든 27 mass 점에서 finite positive `J0`와 기존 geometry predicate를 확인한다. Current trial도 내력 8점·tying·center·recovery와 추가 inertia 평가점에서 허용 geometry를 만족해야 한다. `abs(J)`, 점 생략, 강제 평면화는 없다. 물리 내력·물리 에너지와 section recovery는 기존 `2×2×2` 적분, drill potential은 기존 element-level 합을 보존한다. + +Lumped positivity 조건은 `m_i>0`, `j_i−s_i²/m_i>0`다: + +\[ +2T_{L,i}=m_i\left|v_{ui}+(s_i/m_i)B_iv_{\phi i}\right|^2 ++(j_i-s_i^2/m_i)|B_iv_{\phi i}|^2. +\] + +두 질량식은 총질량 `Σm_i=∫rho dV0`, uniform translation energy, scalar thickness moments를 공유한다. 일반 속도장과 rigid angular motion의 orbital inertia는 다르다. 같은 현재 nodal `w_i=v_ui+z dDot_i`에 대해 + +\[ +T_L-T_C=\frac14\int\rho\sum_{i,j}N_iN_j|w_i-w_j|^2dV_0\ge0. +\] + +따라서 consistent와 lumped의 동적 응답이나 spectrum 동일성을 가정하지 않는다. 평판의 `s_i=0`, `j_i=t²m_i/12`를 warped geometry 일반식으로 사용하지 않는다. 사용자는 보간된 shell 속도의 운동에너지를 직접 적분하는 consistent mass를 선택했으며 기존 eigen의 lumped 계약은 변경하지 않는다. + +## 3. 객관성·운동량·회복 + +공간에서 강체속도를 가하면 `v_ui=V+Omega×(Xi+ui)`, `v_phii=J_l^−1 Omega`이고 `dDot_i=Omega×d_i`다. Consistent는 `xDot=V+Omega×x`를 재현하고, lumped는 각 `x_i(z)`의 같은 강체속도를 적분한다. 따라서 각자의 운동에너지 measure로 회전 관성을 검증해야 한다. + +Consistent의 물리적 총운동량·각운동량은 + +\[ +P=\int\rho\dot x\,dV_0,\quad L_O=\int\rho(x-O)\times\dot x\,dV_0. +\] + +Lumped는 `Σ_i ∫rho N_i xDot_i dV0`, `Σ_i ∫rho N_i(x_i−O)×xDot_i dV0`다. 이 정의에서 시간미분한 `Pdot,Ldot`와 외력·support reaction의 총 force/moment를 대조한다. Kinematic 변환은 algorithmic `A`가 아닌 물리적 `a`로 평가한다. Objective `U_phys+U_drill`의 전체 내력은 강체 가상운동에 일을 하지 않아야 한다. Numerical drill의 개별 nodal moment를 physical stress/resultant로 오인하지 않는다. + +고정 proper global frame change 아래 모든 벡터·구속·하중을 함께 변환하면 scalar `T`, `U`와 virtual work가 불변이어야 한다. 기준축 일부만 회전시킨 비교는 이 검증에 해당하지 않는다. 중간 회전 branch를 통과하는 rechart는 지원하지 않는다. + +## 4. Mass-null과 초기 일관성 + +`B_i`는 rank 2이며 instantaneous drilling 방향은 + +\[ +n_i(q)=\begin{bmatrix}0\\J_l(\phi_i)^{-1}d_i\end{bmatrix},\quad +B_iJ_l^{-1}d_i=0. +\] + +Consistent와 lumped 모두 이 방향에 `M n_i=0`, `n_i^T I=0`이다. 이 null은 일반적으로 `[0;d0]` 또는 고정 global `UR3`가 아니다. 필요한 geometry/rank 조건 아래 물리적 5 DOF/node 외 추가 null이 없는지 검사한다. 구속 후에는 `M_ff`와 free admissible variation 안에서의 null basis `N(q)`를 사용하며, constrained component를 단순 삭제한 nodal null을 그대로 basis라 부르지 않는다. + +`g=f_int−F`라 두면 연속 대수 조건은 + +\[ +\psi(q,t):=N(q)^Tg_f(q,t)=0. +\] + +Physical energy는 director-parallel virtual work가 0이다. 본 force-only subset에서 `psi`는 numerical drill channel의 stationarity를 포함한다. Drilling stiffness는 algebraic coordinate를 결정할 수 있지만 mass rank를 바꾸지는 않는다. `N`의 q 미분도 포함하여 + +\[ +\psi_q v_f+\psi_t=0,\qquad +\psi_q a_f+\psi_{qq}[v_f,v_f]+2\psi_{qt}v_f+\psi_{tt}=0 +\] + +를 초기 연속 일관성으로 요구한다. 첫 시간은 `0+`이며 이후 constant load의 우측 도함수 `Fdot=Fddot=0`를 사용한다. STEP jump를 impulse로 해석하여 초기 velocity를 바꾸지 않는다. + +`q0=v0=0`, 고정 zero boundary에서는 다음 절차로 `a0`를 정할 수 있다. + +1. `psi(q0,0+)=0`와 1차 조건을 확인한다. 여기서 물리적 initial acceleration을 풀기 전부터 불일관하면 실패다. +2. `b=F_f0−f_int,f(q0)−c_f(q0,v0)`에 `N0^T b=0`이 필요하다. 양의 질량 부분공간에서 `M_ff a_P=b`, `N0^T a_P=0`인 particular solution을 구한다. +3. `a_f0=a_P+N0 eta`를 놓고 + +\[ +(\psi_qN0)\eta=-\psi_q a_P-\psi_{qq}[v0,v0]-2\psi_{qt}v0-\psi_{tt} +\] + +를 푼다. 이 subset에서 우변 뒤 세 항은 0이다. `psi_qN0`가 nonsingular여야 초기 acceleration이 유일하게 정해진다. 단순 Moore–Penrose solve 결과를 그대로 가속도라고 채택하거나 massless acceleration을 0으로 놓지 않는다. + +초기점에서는 `g0`가 병진 force를 포함하므로 `D N` 항을 형식적으로 삭제하지 않는다. 다만 zero rotation에서 실제 vanished 항은 식으로 확인할 수 있다. Null basis의 normalization/부호는 물리적 판정에 영향을 주지 않으며 source identity와 결정적 순서를 유지한다. + +Bounds/regularity는 다음을 포함한다: 양의 `rho,t`; reference/current valid geometry; `||phi||_results.h5`다. `U/UR`, 물리적 `v/a`, `omega/omega_dot`, algorithmic `A`, 반력과 에너지의 **공개 여부 및 영속 schema는 미결정**이다. Algorithmic `A`를 physical acceleration이라는 이름으로 공개해서는 안 된다. 출력 선택과 무관하게 검증에는 서로 구별되는 내부 quantity가 필요하다. + +## 7. 검증 인계와 미해결 결정 + +| 항목 | 필수 evidence / 남은 경계 | +| --- | --- | +| HHT 계수·initialization | 표준 일정 질량 HHT와 algebraic reduction; SDOF 자유진동/STEP onset 및 `a_0=F/m`; `alpha=0` 한계 | +| 관성·회전 | Consistent mass의 총질량, uniform translation, 유한회전 angular motion, `d·ḋ=0`, `d·d̈+|ḋ|²=0`, objectivity, kinetic-energy derivative·tangent directional derivative | +| mass-null | 구성 의존 null, `ψ=0`과 초기 1·2차 일관성; coupled massless DAE에서 q/v/a 시간 수렴, finite-rotation drill coupling. 이 source theorem의 직접 적용은 미확인 | +| 정확도 | Smooth `M(q)` manufactured case의 q/v/a를 h,h/2,h/4로 검산; constant-M와 variable-M를 분리. STEP의 시작 불연속과 이후 smooth 구간을 분리 | +| 에너지·반력 | 관성 구속 평형, generalized/spatial moment work, algorithmic/endpoint reaction 구별, `T/U_phys/U_drill/W`와 momentum evidence | +| 구현 lifecycle | 동일 state tangent, 실패 rollback, time/frame 순서, deterministic assembly와 반복 결과 | +| 외부 reference | 지정 원본 101 frames × 9 nodes, frame/instance/node/component의 exact mapping; U blocking / UR warning-only 및 공통 TOLERANCE. V/A/energy correctness는 해당 CSV에 coverage 없음 | + +수식 후보의 [독립 검토](numerical-review.md)는 완료되었으며 표시된 수식에서 새 오류는 발견되지 않았다. 질량 선택은 consistent로 확정되었다. 구현 전 남은 계약은 (1) HHT 확장·초기화의 FESA 채택과 나머지 physics 범위, (2) initial algebraic regularity와 trial 실패 판정의 수치 기준, (3) 공식 output 의미의 공개 범위/I/O다. 이 선택 후 구현할 수 있는 수식 후보이며, 실제 유한회전 shell의 drilling DAE·관성 tangent·q/v/a 시간 수렴은 **구현 후 기능 acceptance에 필요한 검증**이다. 부모 문헌의 nonsingular-mass 조건을 shell mass-null에 자동 적용하지 않되, 모든 shell DAE의 보편적 수렴 증명을 코드 작성의 선행 조건으로 요구하지 않는다. 현재 문서 완료나 제한된 probe를 전체 solver 승인·검증 통과로 표시하지 않는다. + +이번 owner 작업은 기존 정식화·source 대조, 일차 문헌 대조와 수식 유도 및 문서 정적 검토다. [독립 시간정확도 검산](numerical-review-probes.md)은 variable-mass SDOF와 moving-null 예제에서 추천 후보의 약 2차 q/v/physical-a 수렴을 관찰했고 naive 확대의 약 1차 오차를 구별했다. 이는 shell solver 검증이나 일반 drilling DAE 증명을 대신하지 않는다. Build/test·solver/reference 실행은 문서 작업 범위가 아니므로 미실행이다. 독립 검산과 Numerical Review 결과는 해당 owner의 보고서가 소유하며 본 문서가 스스로 승인하지 않는다. diff --git a/docs/hht-implicit-dynamics/implementation-plan.md b/docs/hht-implicit-dynamics/implementation-plan.md new file mode 100644 index 0000000..60d65e0 --- /dev/null +++ b/docs/hht-implicit-dynamics/implementation-plan.md @@ -0,0 +1,43 @@ +# HHT Implementation Plan + +**Goal:** 일관질량 HHT shell 동적 해석을 원본 입력부터 공식 HDF5·전체 reference 비교·독립 물리 평가까지 구현한다. + +**Architecture:** Domain의 별도 dynamic step과 별도 동적 state/procedure/writer를 연결한다. 기존 MITC4 내력·objective drilling은 보존하고 유한회전 consistent 관성을 정확히 미분한다. 물리적 a와 보조 A를 분리하고 candidate 검증 후에만 frame·최종 파일을 반영한다. + +**Tech Stack:** C++17, MSVC x64 Debug, CMake/Ninja/CTest, MKL/TBB/HDF5. + +**Spec:** [정식화](formulation.md), [관성](formulation-inertia.md), [내부 기준](numerical-controls.md), [검증 사례](verification-cases.md), [I/O](io-contract.md), [Reference](reference-model.md). + +## 실행 범위와 소유권 + +사용자가 I/O definition, implementation, physics evaluation을 모두 요청했으므로 기존 정식화의 추천 a/A 분리와 0+ 초기화를 구현 대상으로 채택한다. 메인은 통합·procedure·검증을 맡고 독립적인 kernel/input/writer 묶음만 구현 담당자에게 할당한다. 현재 checklist는 [PRD](PRD.md)에만 유지한다. Commit·push는 이번 요청에 포함되지 않는다. + +| 묶음 | 파일·책임 | 검증 | +| --- | --- | --- | +| 입력 | `model_types.h`, `domain.h/.cpp`, Abaqus mapping: `DynamicStepDefinition`, `StepDefinition::Dynamic/IsDynamic`, `Domain::Dynamic` | 원본 INP, STEP·DIRECT·density·zero BC·force subset, unsupported·invalid input; static/eigen 회귀 | +| 관성 | `elements/nonlinear_mitc4_inertia.h/.cpp`: `Create(nodes,frames,section,material,density)`, `Evaluate(trial,v24,a24)` → M/I/c/DqI/DvI/T/P/L/rates | C05–C06 독립 volume·미분·null·좌표/단위 oracle | +| 상태·출력 | `analysis/dynamic_analysis_state.h/.cpp`, `results/dynamic_results_writer.h`, `io/hdf5/dynamic_hdf5_results_writer.h/.cpp` | frame0 포함 commit/rollback, finite/identity/schedule 검증, schema readback·atomic failure | +| 시간적분 | `analysis/hht_time_integration.h/.cpp`와 tests | C01–C04 정확해·HHT 환원·moving-null·초기화·소산 | +| 동적 절차 | `analysis/hht_dynamic_procedure.h/.cpp` 및 필요한 internal assembly/constraint helper | 초기 positive-range/null solve, 정확한 유효 접선, N1–N7, C07–C11 shell·관성반력·에너지·실패 | +| 통합 | `prepared_analysis.h/.cpp`, application, CMake 등록 | 원본 CLI 실행·진단·최종 HDF5; 기존 procedures 보존 | +| 비교·물리 | `tests/reference/hht_dynamic_reference_*`, 보고서와 ignored 실행 evidence | C12 202 comparison families, 101 frames 전부, 독립 physics evaluation | + +## 구현·검증 순서 + +각 묶음은 관찰 가능한 실패 사례를 먼저 추가하고 구현 후 관련 CMake target/CTest filter로 확인한다. 아직 없는 API의 컴파일 실패와 수치 oracle 실패를 구별해 기록한다. 코드 소유권이 겹치는 변경은 메인이 직렬로 통합한다. Build는 메인이 단일 경로에서 실행한다. + +1. I/O 계약의 필드·진단·source/time 연결을 정의하고 기존 문서와 대조한다. +2. Input/관성 kernel/동적 state를 구현한다. 관성의 director D3와 convective 항은 production finite difference로 대체하지 않는다. +3. HHT predictor와 물리적 endpoint 평형을 연결한다. `A=(q-q_p)/(beta*h*h)`, `v=v_p+gamma*(q-q_p)/(beta*h)`, `a=(A+alpha*a_n)/(1+alpha)`이며 tangent는 `DqI+gamma/(beta*h)*DvI+M/((1+alpha)*beta*h*h)+Kint`다. +4. N2 잔차·correction/rate AND 판정, N3 초기 대수조건, N4 analytic null·positive mass·regularity 및 실패 보존을 검증한다. Backend 성공만으로 frame을 commit하지 않는다. +5. Official HDF5 writer와 deterministic comparator를 연결한다. Reference의 source increment로 대응하고 시간 보간·행 평균·값 보정은 하지 않는다. +6. 수정 영향에 맞는 검사를 수행하고 최종 `python scripts/execute.py`로 전체 MSVC Debug build/CTest를 한 번 수행한다. 필수 원본 dynamic reference는 별도 비교 evidence를 남긴다. +7. 독립 physics reviewer가 실제 결과의 평형·운동량·에너지·coverage를 평가한다. 실패는 현 owner가 수정하며 tolerance를 결과에 맞추어 변경하지 않는다. + +## Review focus + +- a/A 혼동 또는 current M만의 naive HHT 확대: C01–C04와 output readback으로 검출한다. +- 질량 영공간을 작은 질량·fixed UR3 삭제로 처리: C04–C06의 초기·moving-null 사례로 검출한다. +- Initial full STEP force와 초기 acceleration 누락: C02와 원본 frame0의 U/V/A 의미를 대조한다. +- 실패 trial·불완전 HDF5 노출: state/writer 실패 주입 및 기존 최종 파일 보존을 검사한다. +- 101-frame identity·reference subset 오류: C12의 누락·추가·중복·nonfinite·wrong-time mutation을 사용한다. diff --git a/docs/hht-implicit-dynamics/io-contract.md b/docs/hht-implicit-dynamics/io-contract.md new file mode 100644 index 0000000..d6c664b --- /dev/null +++ b/docs/hht-implicit-dynamics/io-contract.md @@ -0,0 +1,153 @@ +# HHT I/O 계약 + +- 상태: `io-contract-defined; implementation-and-readback-pending`, 2026-09-22. +- Owner: `io-definition-agent`. 메인 인계로 새 procedure의 아래 공개 subset/schema를 확정한다. 물리·수치 기준의 원본은 [formulation](formulation.md), [inertia](formulation-inertia.md), [numerical-controls](numerical-controls.md)다. +- 기준: [AGENTS](../../AGENTS.md), [프로젝트 PRD](../PRD.md), [공통 ADR](../ADR.md), [ARCHITECTURE](../ARCHITECTURE.md), [Feature PRD](PRD.md), [Feature ADR](ADR.md). 기존 문서의 `I/O pending`은 본 계약 이전 상태다. + +## 1. 입력과 호환성 + +기존 [비선형 I/O](../nonlinear-static-newton-raphson/io.md)의 lexical parsing, source target resolution, part/identity instance, node/set, S4/S4R → FESA-MITC4, centered homogeneous shell, isotropic elastic material, geometry와 chart 검증을 상속한다. 아래 절차 차이만 dynamic에 적용한다. 기존 static/frequency 입력·결과의 의미나 schema는 바꾸지 않는다. Source S4/S4R 표기는 보존하며 Abaqus 내부 요소·관성·시간 알고리즘 동등성을 주장하지 않는다. 공개 Abaqus 관례의 근거는 이미 조사한 [research](research.md)의 관측 가능한 입력 의미에 한정한다. + +| 구문 | 이번 지원 의미 / 검증 | +| --- | --- | +| `*STEP, NAME=, NLGEOM=YES, INC=` | 한 step만 허용. NAME은 기존 규칙, 생략 시 `Step-1`; NLGEOM=YES 필수. INC 생략값 100, 명시값 양의 uint32. INC는 시간 increment 상한이며 Newton 반복 상한이 아님 | +| `*DYNAMIC, DIRECT` | step의 첫 procedure. DIRECT는 값 없는 flag, 다른 parameter는 거부. ALPHA/BETA/GAMMA, APPLICATION, INITIAL, EXPLICIT, 자동 증분은 미지원 | +| data `h,T` 또는 `h,T,` | 정확히 한 row, 두 finite positive 값. h≤T. `N=round(T/h)`가 uint32 범위의 양의 정수이고 `double(N)*h == T`, N≤INC일 때만 허용. 중간 누적 덧셈·끝 증분 단축·time tolerance 없음. 그 밖의 nonempty field/row 거부 | +| `*DENSITY` 아래 `rho` 또는 `rho,` | 활성 material마다 정확히 한 finite positive density; temperature/dependency/추가 field 미지원. 질량 단위는 `force*time^2/length`, density는 `force*time^2/length^4`. Frequency의 기존 허용 유지, static에 새 허용하지 않음 | +| `*BOUNDARY` | 모델 앞 또는 procedure 뒤/하중 앞. 기존 target, DOF 1–6 범위; magnitude 생략 또는 exact 0. 전체 시간에 q/v/a/A constrained 성분 exact 0. 초기 q/v와 충돌하는 값, AMPLITUDE/OP/TYPE/변경 경계 거부 | +| partial zero BC | component selector 의미를 그대로 보존. 전 DOF clamp/free 필수 지원, 부분구속은 N4 analytic null과 free 공간 교집합 사용. Rank-changing branch/regularity 실패는 오류이며 부분구속을 full clamp로 바꾸지 않음 | +| `*CLOAD` | step 내 global dead force DOF 1–3만 허용. finite magnitude; 동일 resolved node/DOF의 여러 row는 source order로 합산. DOF 4–6은 값이 0이어도 미지원. AMPLITUDE/FOLLOWER/OP 및 모델-level 하중 미지원 | +| `*PREPRINT, ECHO=NO, MODEL=NO, HISTORY=NO, CONTACT=NO` | 정확한 flags, data 없음: `ignored-output-request` warning. 출력 제어 효과 없음 | +| `*RESTART, WRITE, FREQUENCY=0` | data 없는 이 조합만 warning no-op. Restart 생성·읽기 지원 아님 | +| `*OUTPUT, FIELD, VARIABLE=PRESELECT, FREQUENCY=1` | data 없음. HISTORY도 같은 조합 허용, FIELD/HISTORY 둘 중 하나. FREQUENCY 생략도 1 의미로 수용. 요청마다 warning no-op; FESA 필수 전 frame 출력을 바꾸지 않음 | + +`*NODE OUTPUT`, `*ELEMENT OUTPUT` 등 child selectors와 다른 OUTPUT parameter/frequency는 미지원이다. 원본 [INP](../../reference/shellNL-dyn-meshfine/shellNL-dyn-meshfine.inp)는 위 root 요청만 사용한다. 초기조건 keyword, amplitude, damping, distributed/follower load, contact, 비탄성, 다른 dynamic 요소, 다중 step/restart는 fail closed다. 기존 shell section의 source integration-point field `5`는 기존 allowlisted warning 의미를 유지하며 FESA 관성 적분 선택이 아니다. + +공개 계수는 기본 `alpha=-0.05,beta=0.275625,gamma=0.55`만 제공한다. Internal 검증용 alpha=0은 formulation의 결합계수로 실행할 수 있으나 INP 지원 확대가 아니다. 질량은 `consistent-kinetic`; 사용자 선택을 lumped fallback으로 바꾸지 않는다. 초기 frame은 `0+`, q0=v0=0, F0=F*, A0=a0이며 physical a0와 null 일관성은 N3에 따라 구한다. 하중은 이후 모든 endpoint에도 전값이다. + +## 2. Identity, 시간, state 경계 + +Canonical storage slot은 `/steps/Step-1`이다. 실제 source name을 별도 `source_step_name`에 저장한다. Reference의 실제 name은 `Step-1`이다. Instance는 기존 ASCII case-insensitive canonical identity와 raw label text를 보존하고, logical label은 positive int64다. Internal node/element ID와 equation ID를 source label 대신 사용하지 않는다. + +N increments의 완성 결과는 frame 0..N, increment 0..N을 각각 한 번 포함한다. Frame n의 `time=double(n)*h`, frame0=+0, `is_initial=(n==0)`이며 n≥1에서 시간이 strictly increasing해야 한다. 종료시각은 T와 exact 일치한다. Frame0 iteration_count=0; 이후 0..100 corrections, 초기 추정 수렴/empty free solve는 0 허용. INC에 frame0는 포함하지 않는다. Original h=.1,T=10,INC=1000은 N=100/101 frames다. + +Domain이 dynamic semantic 정의를 소유한다. AnalysisModel은 non-owning view, DofManager는 full/free/constraint mapping을 소유한다. State vector의 node-major 6성분은 `[u1,u2,u3,phi1,phi2,phi3]`이며 node inventory는 Domain 선언 순서다. Rotation은 fixed-initial global total rotation-vector chart다. Trial은 committed state/history와 alias되지 않는다. Candidate의 time, q/v/a/A, orientation, force, recovery와 evidence가 모두 검증된 뒤 한 번에 commit한다. Failure는 마지막 committed frame과 history를 보존하며 failed/trial frame을 노출하지 않는다. + +### 구현 인계용 C++ 경계 + +이는 backend-neutral 경계이며 vendor type을 포함하지 않는다. Default-constructed record는 구성 중 candidate일 뿐 유효한 완료 state가 아니다. + +```cpp +struct DynamicStepDefinition { + std::string source_step_name_raw, source_step_name_canonical; + SourceLocation step_location, dynamic_location, data_location; + double initial_increment = 0.0, time_period = 0.0; + std::uint32_t increment_count = 0, maximum_increments = 100; + double alpha = -0.05, beta = 0.275625, gamma = 0.55; + std::vector boundaries; + std::vector loads; +}; + +struct DynamicAcceptedFrameState { + std::uint32_t frame_index = 0, increment_index = 0, iteration_count = 0; + double time = 0.0; + bool is_initial = false; + std::vector nodal_internal_order; + std::vector nodal_source_ids; + Vector generalized_coordinates, velocity, acceleration; + Vector algorithmic_acceleration; // Each 6*N; physical a != auxiliary A. + std::vector, 3>> nodal_triad; + std::vector> nodal_director; + Vector angular_velocity, angular_acceleration; // Each 3*N, spatial. + Vector inertia_force, internal_force, external_force, residual, reaction; + Vector residual_scale; // Each 6*N, N2 uncancelled per-row scale. + std::vector shell_results; + double kinetic_energy = 0.0, physical_strain_energy = 0.0; + double drilling_energy = 0.0, external_work = 0.0, energy_defect = 0.0; + std::array global_equilibrium{}; + std::array linear_momentum{}, angular_momentum{}; + std::array linear_momentum_rate{}, angular_momentum_rate{}; + std::array free_residual_metrics{}; // u, phi channels. + Vector hidden_constraint, hidden_velocity_defect, hidden_acceleration_defect; + Vector hidden_constraint_scale, hidden_velocity_defect_scale; + Vector hidden_acceleration_defect_scale; // Each analytic null count nz. +}; +``` + +`ModelDefinition.dynamic_step:optional`, `StepDefinition::IsDynamic()`, `StepDefinition::Dynamic()` 및 `Domain::Dynamic()` accessor를 추가한다. Dynamic과 static/frequency 정의는 상호배타적이다. `DynamicAnalysisState::Create(const Domain&,const DofManager&):Result`는 immutable inventory/schedule을 검증한다. 기본 empty 구성이 가능하되 Create 전 commit/실행/Write를 거부한다. `BeginTrialFromCommitted():Result`는 deep copy, `CommitAcceptedFrame(const DynamicAcceptedFrameState&):Status`는 강한 실패 보존, `DiscardTrial():void`는 caller-owned trial 폐기를 위한 no-op이다. `AcceptedFrames():const vector&`, `CommittedState():const DynamicAcceptedFrameState&`를 제공한다. Frame0는 Create 후 완성하여 같은 commit 검사를 통과해야 한다. 첫 commit 전 CommittedState 읽기는 금지다. Full q6N이 authority이며 displacement/rotation state를 중복 소유하지 않는다. + +`DynamicExecutionResult { DynamicAnalysisState state; vector diagnostics; }`를 PreparedAnalysisOutcome에 추가하고 별도 dynamic executor를 dispatch한다. `DynamicResultsWriter::Write(const filesystem::path&, const Domain&, const DynamicAnalysisState&, const vector&):Status` 및 `DynamicHdf5ResultsWriter`를 추가한다. Existing static writer에 time 또는 static load_factor로 dynamic을 위장하지 않는다. + +Momentum/rates는 전체 두께 volume kinetic kinematics로 구한 physical quantities, 원점 O=(0,0,0)에 대한 global angular momentum이다. Chart generalized momentum을 spatial angular momentum으로 표시하지 않는다. `global_equilibrium`은 N6의 applied+support minus momentum-rate이며 current arms와 spatial moments를 사용한다. Hidden vectors는 N3의 raw ψ/1차/2차 defects, `_scale`는 그 항별 uncancelled norm scale다. n_z=0은 empty vector, empty free channel metric 저장값 0은 해당 constraint inventory와 함께 N/A로 해석한다. NaN sentinel을 사용하지 않는다. Hidden 1·2차 defect는 진단/시간수렴 대상이며 Newton 수렴 기준으로 임의 승격하지 않는다. + +## 3. 공식 HDF5 schema v1 + +공식 경로는 input lexical sibling `_results.h5`, CLI/lock/log/atomic file lifecycle은 기존 I/O §7.1을 상속한다. Dynamic log는 `nsnr-log-v1` 또는 target_factor를 사용하지 않는다. 기존 log stream에 dynamic diagnostic을 기록할 경우 `record=dynamic-newton-iteration`, `schema=hht-log-v1`, increment/time/iteration/stage와 N2 metric의 raw numerator/scale/branch/pass를 명시한다. 정적 전용 tolerance를 재사용하지 않는다. + +Schema identity는 `(feature_id,analysis_procedure,schema_version)`다. `/metadata` attribute: `feature_id=hht-implicit-dynamics`, `analysis_procedure=hht-implicit-dynamics`, `schema_version:uint64=1`, `mass_formulation=consistent-kinetic`, `internal_formulation=FESA-MITC4-total-lagrangian`, `initial_state=zero-q-v-step-load-0+`, `acceleration_convention=physical-endpoint-and-auxiliary-separated`, `source_step_name`, `step_name=Step-1`, `alpha/beta/gamma:float64`, `initial_increment/time_period:float64`, `increment_count/maximum_increments:uint64`. 기존 solver_version/source_input_identity/unit_system_label도 필수, 단위 label은 `consistent-input-units`이며 SI로 추정하지 않는다. + +모든 numeric data는 finite-first, float64=IEEE binary64 little-endian, integer는 fixed width little-endian, string은 UTF-8 variable length다. Mandatory datasets는 contiguous/uncompressed. 아래 N=nodes, E=elements, F=Nincrements+1이다. Paths의 ``은 leading-zero 없는 decimal ordinal이다. + +Model schema는 기존 I/O §7.3의 `/model/nodes`, `/model/elements`, `/model/shell/materials`, `/model/shell/sections`, `/model/shell/nodal_director`, `/model/shell/nodal_frame`, `/model/nodal_constraint_mask`, `/model/shell/midsurface_locations`, `/model/shell/section_positions`의 정확한 dtype/shape/identity를 재사용한다. `/model/prescribed_displacement_start/target`은 dynamic에서 쓰지 않는다. 추가 `/model/prescribed_displacement:[N,6] float64`는 전값 0, `/model/material_density:[material_count]` compound는 `internal_material_id:uint64,rho:float64`로 materials 행 순서다. 활성 material의 rho는 finite positive다. + +`/steps/Step-1/accepted_frames:[F]` compound fields는 `frame_index:uint64,increment_index:uint64,time:float64,iteration_count:uint64,is_initial:uint8`. Table의 시간과 모든 frame attrs가 일치해야 한다. Load factor는 없다. + +각 `/steps/Step-1/frames/` 아래 필수 dataset은 다음과 같다. + +| 상대 path | shape, dtype | component / 의미 / 차원 | +| --- | --- | --- | +| `nodal/displacement` | [N,6] float64 | UX,UY,UZ,URX,URY,URZ; q; length×3,1×3 (UR는 rad) | +| `nodal/velocity` | [N,6] float64 | V1,V2,V3,PHIDOT1,PHIDOT2,PHIDOT3; chart v; length/time×3,1/time×3 | +| `nodal/acceleration` | [N,6] float64 | A1,A2,A3,PHIDDOT1,PHIDDOT2,PHIDDOT3; physical chart a; length/time²×3,1/time²×3 | +| `nodal/algorithmic_acceleration` | [N,6] float64 | AUX1,AUX2,AUX3,AUX4,AUX5,AUX6; auxiliary A; a와 같은 차원 | +| `nodal/angular_velocity` | [N,3] float64 | OMEGA1,OMEGA2,OMEGA3; J_l(phi) phidot; 1/time | +| `nodal/angular_acceleration` | [N,3] float64 | OMEGADOT1,OMEGADOT2,OMEGADOT3; J_l phiddot+DJ_l[phidot]phidot; 1/time² | +| `nodal/triad`, `nodal/director` | [N,3,3], [N,3] float64 | 기존 current triad/director 의미, axis row A,B,D / global XYZ | +| `nodal/generalized_inertia_force` | [N,6] float64 | I=M(q)a+c(q,v), physical endpoint | +| `nodal/generalized_internal_force` | [N,6] float64 | physical+objective drilling internal force | +| `nodal/generalized_external_force` | [N,6] float64 | F*, frame0에도 full value | +| `nodal/generalized_residual` | [N,6] float64 | E=I+f_int-F, free와 constrained 전체 | +| `nodal/generalized_reaction` | [N,6] float64 | constrained selector를 E에 적용, free 성분 exact zero | +| `nodal/reaction` | [N,6] float64 | RF1,RF2,RF3,RM1,RM2,RM3; support-on-structure; moment=J_l^-T generalized reaction | +| `element/shell/local_frame` | [E,4,3,3] float64 | 기존 GP1..4, local axes E1,E2,E3 | +| `element/shell/generalized_strain` | [E,4,8] float64 | 기존 E11,E22,G12,K11,K22,K12,G13,G23 | +| `element/shell/section_resultant` | [E,4,8] float64 | 기존 N11,N22,N12,M11,M22,M12,Q13,Q23 | +| `element/shell/stress` | [E,4,3,3] float64 | 기존 BOTTOM,MIDDLE,TOP × S11,S22,S12 | +| `global/energy` | [5] float64 | KINETIC_ENERGY,PHYSICAL_STRAIN_ENERGY,DRILLING_ENERGY,EXTERNAL_WORK,ENERGY_DEFECT; 모두 force*length | +| `global/equilibrium` | [6] float64 | FORCE_1..3,MOMENT_1..3; N6 관성 포함 balance error | + +Generalized force 계열 component_names=`F1,F2,F3,MPHI1,MPHI2,MPHI3`, units=`force,force,force,force*length,force*length,force*length`, coordinate_system=`fixed-initial-global-chart`. Spatial U/omega/reaction은 `global-cartesian`; q/v/a/A는 translation global 및 rotation chart임을 `rotation_representation=total-rotation-vector-chart`와 `quantity_semantics`로 명시한다. Shell units, local_frame/source_element_type/midsurface/section_position 참조 attrs는 기존 I/O §7.5 그대로다. Numerical drill force를 physical shell stress/resultant에 섞지 않는다. + +모든 result dataset에 `schema_version`, `component_names`, `component_unit_dimensions`, `coordinate_system`, `location`, `step_name`, `source_step_name`, `frame_index`, `increment_index`, `time`, `is_initial` attr를 쓴다. Scalar numeric attrs는 위 dtype, flags uint8이다. Location은 `nodal`, shell 기존 location, `global`/`global-balance-error`다. 모든 nodal row는 `/model/nodes`, element row는 `/model/elements` 행을 참조한다. + +Energy의 W=F*ᵀ(qn-q0), defect=T+Uphys+Udrill-E0-W다. Initial physical energy E0는 frame0의 첫 세 성분 합이며 W0/defect0=0. A를 kinetic energy나 physical inertia에 쓰지 않는다. Algorithmic dissipation을 physical damping energy로 이름 붙이지 않는다. Raw momentum/hidden diagnostics는 위 internal frame의 필수 evidence이며 이번 persistent schema에는 별도 필수 dataset을 추가하지 않는다. HDF5에서 계산할 수 없는 상세 검증은 실행 보고서가 internal evidence를 기록한다. + +`/diagnostics`는 기존 diagnostic compound의 severity/code/source_file/source_line/keyword/entity_identity/message와 정렬을 재사용한다. Writer는 전체 frames/node/element inventory, dtype/rank/component/attrs/time/finite/constraint, complete final frame를 검사하고 write/flush/close/reopen/self-check 후 원자 replace한다. 기존 final file은 실패 시 보존한다. Live log의 이미 flush된 prefix는 rollback하지 않는다. No-op warning이 있어도 필수 dataset을 생략할 수 없다. + +## 4. 진단과 비교 binding + +기존 `Diagnostic`의 severity/code/location/keyword/entity_identity/message를 사용하며 `SortDiagnostics` 순서 유지. 위치는 가능한 source keyword/data row, runtime에서는 step+increment와 node/element source identity를 message/entity_identity에 포함한다. Parser/model/output/backend의 기존 FailureCategory/CLI exit mapping을 유지하고 새 exit code를 만들지 않는다. + +| code | 발생 조건 | +| --- | --- | +| `unsupported-dynamic-parameter` | DIRECT 외 dynamic parameter, 잘못된 flag, 계수 입력 | +| `invalid-dynamic-data` | h/T/nonfinite/arity/범위 오류 | +| `invalid-dynamic-schedule` | N*h≠T, N>INC, uint32 초과, 시간 비증가 | +| `unsupported-dynamic-boundary` | 비영/시간변화/BC option | +| `unsupported-dynamic-load` | CLOAD4–6, follower/amplitude/distributed load | +| `invalid-density-data` | 누락/중복/nonpositive/비상수 density; 기존 code 재사용 | +| `unsupported-analysis-procedure`, `unsupported-multiple-step`, `unsupported-keyword` | 잘못된 procedure/NLGEOM, 다중 step, allowlist 밖 의미 | +| `ignored-output-request` | 위 정확한 allowlisted no-op만 warning | +| `dynamic-initial-consistency`, `dynamic-mass-regularity` | N3 초기화/숨은 조건 또는 N4 rank/positive block/regularity 실패 | +| `dynamic-nonconvergence`, `dynamic-linear-solve` | N2 correction 한도 또는 solve/backward-error 실패 | +| `dynamic-invalid-state`, `dynamic-output-schema` | chart/geometry/nonfinite/inventory/update 또는 저장/재독 계약 위반 | + +기존 geometry/backend code가 더 구체적이면 보존하며 위 code로 감싸 원래 source context를 버리지 않는다. 질량 singularity 자체는 오류가 아니다. Formulation에 정의된 허용 null을 제거하거나 epsilon mass를 주지 않는다. + +[Reference 계약](reference-model.md)의 source step/instance/node/increment/time/component authority를 그대로 적용한다. Binding: `/model/nodes`의 `(instance_name,source_label)`와 frame table을 조인하고 `nodal/displacement[:,0:3]→U-U1/U2/U3`, `[:,3:6]→UR-UR1/UR2/UR3`. 먼저 전체 729-node 101-frame identity/finite 검사를 수행하고 고정 9-node subset으로 projection한다. Decimal CSV time=n/10 검증과 binary HDF5 time=double(n)*double(.1) 검증을 각각 수행한다. Nearest-time/보간/평균/row-index matching은 금지다. Wrong procedure/mass/coefficient/초기 convention, missing/extra/duplicate/nonfinite는 tolerance 전에 실패다. Frame별 U blocking/UR warning-only 202 families와 공통 [TOLERANCE](../TOLERANCE.md)는 변경하지 않는다. + +## 5. 인계 및 확인 상태 + +입력 원본의 procedure/density/BC/CLOAD/output 요청과 기존 semantic/prepared_analysis/state/writer 경계를 대조했다. 이번 변경은 이 문서뿐이며 source/test/reference artifact는 수정하지 않았다. 문서 링크와 diff를 정적 검사한다. Build/test·dynamic HDF5 readback·reference 비교는 문서 owner의 작업이 아니므로 미실행이다. + +상위 물리 선택 충돌은 없다. 구현에 아직 필요한 upstream quantity는 별도 Udrill energy, exact inertia/kinetic energy, spatial rates, physical momentum/rates, N2 uncancelled per-row scale 및 N3 hidden defects다. 기존 nonlinear recovery가 Uphys만 반환한다는 이유로 Udrill=0 또는 누락으로 바꾸면 안 된다. 메인은 이 문서를 main report에서 연결하고 공개 subset/schema 결정 근거를 Feature ADR에 기록한다. 이 계약 정의만으로 solver/physics gate가 통과한 것은 아니다. diff --git a/docs/hht-implicit-dynamics/numerical-controls.md b/docs/hht-implicit-dynamics/numerical-controls.md new file mode 100644 index 0000000..889610c --- /dev/null +++ b/docs/hht-implicit-dynamics/numerical-controls.md @@ -0,0 +1,98 @@ +# HHT 내부 수치 판정 + +[Numerical Review](numerical-review.md)의 구현용 계약이다. IEEE binary64를 기준으로 하며 [정식화](formulation.md), [관성 상세](formulation-inertia.md), [검증 사례](verification-cases.md)를 함께 적용한다. 사용자 결정 [HHT-ADR-007](ADR.md)에 따라 이번 HHT 구현은 **consistent kinetic mass**를 사용한다. 아래 lumped 대조는 기존 eigen/대안 evidence의 경계이며 새 HHT lumped 지원 요구가 아니다. 외부 reference는 [TOLERANCE](../TOLERANCE.md)만 따른다. + +## N1. Scale와 공통 규칙 + +`ell`은 활성 요소의 **초기** 최대 edge 길이이며 finite positive여야 한다. `D=diag_i(ell I3,I3)`, `x=D^-1 q`, `w=D^-1 v`, `b=D^-1 a`, `B=D^-1 A`, `r=D^T E`, `Khat=D^T K D`, `Mhat=D^T M D`를 사용한다. x는 무차원, w는 시간^-1, b/B는 시간^-2, r/Khat은 energy, Mhat은 mass·length²다. D는 분석 내 고정이다. Element 독립 비교에는 해당 element의 ell을 양쪽에 동일 적용한다. + +별도 표시 없는 vector norm은 infinity norm, matrix norm은 infinity norm이다. Frobenius/2-norm은 식에 표시한다. Euclidean/Frobenius norm에는 scaled sum-of-squares, 조립에는 stable source 순서의 reduction을 사용한다. 원 force와 moment, 원 translation과 rotation을 직접 합하지 않는다. N2의 u/phi 채널도 분리한다. + +모든 검사는 finite-first다. 항등식 ratio의 분자=분모=0은 metric 0, 분모=0·분자≠0은 실패다. 필요한 positive mass, 길이, stiffness, conditioning scale가 0이면 이 zero-pass 대상이 아니다. `max(1,physical_scale)`, absolute-force floor, epsilon mass, rank 잘라내기, raw 값 clamp를 금지한다. Nonzero 값이 필요한 변환에서 underflow-to-zero/overflow이면 진단하고 실패한다. + +## N2. Full Newton 종료와 선형 solve + +매 trial의 원래 free endpoint 식 `E=Ma+c+f_phys+f_drill-F`를 다시 계산한다. 각 free row i의 uncancelled scale는 + +`s_i = sum_element,j |(M_e)_ij a_ej| + sum_element |c_ei| + sum_element |f_phys,ei| + sum_element |f_drill,ei| + sum_load |F_load,i|`. + +`Ma+c`나 physical/drill을 먼저 상쇄한 값을 분모로 사용하지 않는다. q/v/a는 같은 trial이며 원 mass와 physical a를 쓴다. 채널 `J=u,phi`에 대해 `S_J=max_(i∈free J) s_i`, `rho_J=max|E_i|/S_J`. 반력과 constrained row는 Newton norm에서 제외한다. 외력만 0이거나 acceleration turning point여도 남은 항이 scale를 제공한다. 모든 항이 0인 branch는 E의 exact zero만 통과한다. Zero vector 근방에서 임의 분모를 만들어 통과시키지 않는다. Empty free channel은 N/A다. + +| ID | 판정 | 고정 값 | +| --- | --- | --- | +| N2-R | Force와 moment의 `rho_u`, `rho_phi` | 각각 `<=1e-10` | +| N2-Q | 마지막 적용 correction의 `||Δu||/ell`, `||Δphi||` | 각각 `<=1e-10` | +| N2-V | `h||Δv_u||/ell`, `h||Δv_phi||` | 각각 `<=1e-10` | +| N2-a | `h²||Δa_u||/ell`, `h²||Δa_phi||` | 각각 `<=1e-10` | +| N2-A | `h²||ΔA_u||/ell`, `h²||ΔA_phi||` | 각각 `<=1e-10` | +| N2-L | Scaled solve `Khat_ff Δx=-r_f`: `||Khat Δx+r||/(|| |Khat| |Δx|+|r| ||)` | `<=1e-12`; finite-first, 원 unsymmetric operator로 검사 | + +모든 적용 가능한 gate의 **AND**다. a와 A의 서로 다른 chain-rule 계수를 반영한다. h-based rate gate는 증가량의 운동학 일관성을 검사하며 physical acceleration의 시간 정확도를 인증하지 않는다. 해당 정확도는 N7의 독립 검증으로 확인한다. + +Newton 초기 추정은 constant-acceleration `q_trial=q_p+beta h² A_n`으로 둔다. 이는 Newton 시작값이며 시간 scheme 변경이 아니다. 모든 trial에서 고정 q/v/a/A 성분은 exact 0이다. Residual이 이미 통과하는 초기 추정은 correction=0으로 명시하여 동일 gate를 평가한다. 그 외에는 `Keff`를 매회 새로 조립하고 correction을 수정 없이 적용한 후 residual을 다시 검사한다. Initial guess가 chart/geometry 밖이면 실패이며 clamp·line search·cutback으로 바꾸지 않는다. + +최대 **100개의 Newton correction**을 허용한다. 100번째 correction 뒤의 candidate가 모든 gate를 만족하면 통과, 아니면 실패하며 101번째 solve는 없다. 기존 정적 100회/no-cutback 제한을 동일 실패 상한으로 채택한 것이고, 동적 문제가 100회 안에 수렴한다는 보장은 아니다. Backend success나 작은 correction 하나만으로 성공하지 않는다. Singular/invalid/nonfinite/linear backward-error failure는 즉시 실패다. + +## N3. Update·초기화·대수조건 + +HHT 기본 `alpha=-.05,beta=.275625,gamma=.55`, `A0=a0`, `A1=(1+alpha)a1-alpha*a0`. `q,v`는 A로 적분하고 endpoint I에는 a를 넣는다. 각 update 식의 결함을 **동일 차원의 우변 항 절댓값 합**으로 나누어 `<=1e-12`를 요구한다. 큰 q 두 값을 빼서 잃은 작은 increment를 숨기지 않도록 `q1-qn`과 저장된 predictor/correction의 일치도 함께 검사한다. Fixed boundary의 exact zero는 이 ratio보다 우선한다. + +Scaled analytic null basis Z는 `Mhat_ff Z=0`, `Z^TZ=I`이며 ψ=`Z^T(D^Tg)_f`다. Physical inertia의 null action도 N4로 검사한다. 초기 ψ는 `|| |Z|^T (D^T s)_f ||`로, 초기 velocity/acceleration hidden residual은 각각 `|| |ψ_x| |w|+|ψ_t| ||`, `|| |ψ_x| |b|+|ψ_xx|[|w|,|w|]+2|ψ_xt| |w|+|ψ_tt| ||`로 나누어 `<=1e-10`을 요구한다. Matrix/tensor contraction 내부 상쇄도 분모에서 제거한다. 초기 acceleration solve는 physical range equation과 null consistency equation 둘 다 통과해야 한다. 비영 ψ에 대해 초기 q/v를 몰래 보정하지 않는다. + +초기 q0=v0=0, STEP의 `0+` force 전값과 우측 Fdot=Fddot=0을 사용한다. `M a_P=F-f_int-c`, `Z^T D^-1 a_P=0` 후 `a0=a_P+D Z eta`의 η를 differentiated ψ 식으로 결정한다. Null acceleration을 0으로 고정하거나 pseudoinverse 결과만 채택하지 않는다. c0=0은 v0=0에서 유도되는 값이지 일반 초기조건 생략 규칙이 아니다. + +Accepted endpoint의 ψ도 `<=1e-10`으로 검사하되 denominator는 `|| |Z|^T (D^T s)_f ||`로 계산한다. 매 frame의 1·2차 hidden-constraint 결함은 같은 항 norm 합으로 정규화해 **raw value와 함께 진단**한다. 이산 v/a가 미분 제약을 exact zero로 만족한다는 주장은 하지 않는다. Hidden defect에 N2 Newton tolerance를 강제하지 않으며, 기능 acceptance는 N7/C04/C08의 시간수렴으로 검증한다. Nonfinite 또는 rank/regularity 실패는 진단 수준에 머무르지 않고 즉시 실패한다. + +## N4. 질량·null·대수 regularity + +Analytic expected null은 각 node의 `[0;Jl^-1 d]`와 **실제 free component 공간의 교집합**이다. Constrained 성분 중 하나라도 정확히 비영인 null vector는 그 node의 admissible null이 아니다. Tiny nonzero를 0으로 취급하지 않는다. Rank-changing partial BC 구성은 regular branch를 벗어나므로 실패한다. 처음 예상한 rank를 numerical eigenvalue cutoff로 다시 정의하지 않는다. 지정 reference의 625 free nodes는 expected positive dimension 3125/null 625이며 실측 rank가 아니다. Z의 nodal 방향은 위 analytic 식을 positive norm으로 정규화하여 연속적으로 정의한다. Trial마다 eigen/SVD basis나 최대성분 부호 규칙을 새로 선택하여 ψ 미분에 불연속을 넣지 않는다. + +| 대상 | 판정·수치 | +| --- | --- | +| `Mhat` symmetry, analytic null action | `||Mhat-Mhat^T||F/||Mhat||F <=1e-12`; `||Mhat Z||F/(||Mhat||F ||Z||F)<=1e-12`. Empty Z는 N/A | +| Basis·complement | `||Z^TZ-I||F/sqrt(nz)<=1e-12`, `||P^TZ||F<=1e-12`, `[P Z]` Gram도 `1e-12`; deterministic basis sign/source ordering | +| Positive mass | Consistent element physical20와 lumped node physical5의 projected scaled block에 Cholesky strict-positive pivot; assembled free physical mass factor도 SPD. Small positive pivot을 없애지 않음 | +| Local conditioning | 각 nonempty element/node positive block의 `rcond_inf=1/(||B||inf ||B^-1||inf) >=2^-32` (약 `2.3283064e-10`). Small-block triangular solves로 계산하고 global dense inverse는 만들지 않음 | +| Lumped scalar positivity | finite `m>0,j>0`, `j-(s/m)*s>0`; `m*j-s*s`의 overflow/underflow에 판정을 맡기지 않음. Independent centered thickness integral과 `1e-12` 대조 | +| Mass/inertia identities | Scalar Ck entry error / `[m_e(t/2)^k]`, 총질량·uniform translation energy·직접 volume kinetic energy·Cholesky reconstruction: 각각 `1e-12`. Nonlinear physical inertia null action과 kinetic power identity는 `1e-10` | + +초기 및 각 accepted state에서 `A_alg=ψ_x Z`의 null-coordinate regularity를 검사한다(A_alg는 HHT 보조가속도 A와 다른 행렬이다). Scale는 numerical drill Hessian의 단순 현재 diagonal 대신 `k_i=sum_e kd_e ||Dχ_e[D Z_i]||²>0`를 사용한다. `W=diag(sqrt(k_i))`, `B_alg=W^-1 A_alg W^-1`. Empty null system은 N/A. 모든 k_i finite positive, factorization nonsingular와 `rcond_inf(B_alg)>=2^-32`를 요구한다. Small synthetic/local blocks는 inverse-action norm을 직접 계산한다. 실제 assembled sparse matrix는 1/∞ norm condition estimator를 사용하고 추정값·알고리즘을 기록한다; 추정치는 정확한 forward-error certificate가 아니므로 같은 ψ/solve/시간수렴 검사를 대체하지 않는다. 추정 불능이나 기준 미달은 실패다. + +Newton 도중 `Keff` 자체는 대칭·SPD일 필요가 없다. Positive mass 및 regular null constraint를 `Keff`의 양정성 검사로 대체하지 않는다. 오래된 N 또는 초기 M을 현재 geometry에 고정하지 않는다. Free support가 전부 fixed인 유효 모델은 0×0 free solve를 건너뛰고 반력을 회복한다. + +## N5. Rotation·geometry·관성 tangent + +기존 [finite_rotation](../../src/fesa/elements/finite_rotation.cpp)와 [nonlinear shell](../../src/fesa/elements/nonlinear_mitc4_shell.cpp)의 predicates를 보존한다: `pi-|phi|>sqrt(DBL_EPSILON)`, proper Q/director의 `1e-12` 일관성, director opposition margin `1e-8`, positive valid Jacobian·polar/swing/twist branch. Mass의 추가 27점도 정식화의 유효성 검사를 적용한다. Small-angle series 전환은 기존 `|phi|=1e-4` 경계 양쪽의 D1/D2/D3 연속성 검증을 포함한다. 허용 geometry를 넓히거나 norm을 재정규화하는 승인이 아니다. + +`x±s p`, `||p||2=1`의 perturbation에서 **모든 상태는 원 q±s Dp부터 재평가**한다. `s=2^-j, j=8,10,12,14,16,18,20`의 고정 sweep을 사용한다. Keff 검증에는 각 q perturbation에 종속 v/a/A까지 재계산하고, partial DqI/DvI는 나머지 독립변수를 고정한다. 각 u-only, phi-only, mixed direction과 nonzero v/a를 검사한다. + +| 검사 | Metric와 한계 | +| --- | --- | +| Nonzero tangent action | `||Khat p-[r(x+sp)-r(x-sp)]/(2s)||2 / (||Khat p||2+||FD||2) <=1e-6` | +| Exact zero action | FD mismatch / `[||Khat||F ||p||2] <=1e-10`; Khat도 0이면 N1 exact-zero branch | +| D3 director | 중앙차분 `D2d(phi±s p)[r,t]`와 D3 비교, nonzero ratio `1e-6`; zero analytic direction은 `||p||||r||||t||`로 나누어 `1e-10` | +| Symmetry | M, conservative internal Hessian만 scaled Frobenius `<=1e-12`; 일반 dynamic Keff에는 symmetry gate 없음 | + +`j=10..18` 안에서 **연속한 두 sweep 점**이 통과해야 한다. 모든 raw error/step을 보고하고 가장 작은 s 하나만 골라 통과시키지 않는다. Truncation 감소와 roundoff plateau를 구별하며, analytic polynomial/zero directions에 인위적인 2차 FD 기울기를 요구하지 않는다. 전체 sweep이 branch 밖이면 fixture 실패이며 한쪽 derivative로 바꾸지 않는다. + +## N6. 불변식·평형·에너지 + +동일 물리 문제의 proper global frame change·단위변환·direct kinetic-energy oracle는 positive physical scale에 대해 비교한다. Matrix congruence와 질량 항등식 `1e-12`, nonlinear objectivity·energy/power·derivative-recovered rates `1e-10`. Zero-energy rigid mode는 `E_ref=sum_e E_e t_e area_e`로 나눈 physical/drill energy `<=1e-10`와 element force action / `[||Khat_initial||F ||x||2]`를 함께 확인한다. E_ref는 테스트의 알려진 물리 stiffness scale이며 Newton residual 분모에 넣는 absolute floor가 아니다. + +전체 force balance는 `sum(F_applied+F_reaction)-Pdot`, moment는 같은 현재 위치와 공간 moment를 쓴 `sum((x-O)×F+M)-Ldot`다. 분모는 같은 식에 들어간 절점/적분점별 force·moment norm의 **상쇄 전 합**이며 채널별 normalized error `<=1e-8`. 기존 비선형 global balance 기준과 같은 값이지만 새 inertial terms, current arms, `Jl^-T`를 포함해 새로 검증한다. Free residual N2와 함께 통과해야 한다. + +Alpha=0 선형 무감쇠 oscillator에서 `max|E_n-E0|/E0<=1e-10`; nonlinear/alpha<0에는 적용하지 않는다. 일반 nonlinear `T+U_phys+U_drill-W`는 부호·최대 절댓값·시간수렴을 보고하며 매 frame 감소·정확 보존을 요구하지 않는다. 자세한 수치 bounds는 C07/C08에 고정한다. HHT algorithmic dissipation을 physical damping으로 기록하지 않는다. + +## N7. 수렴·결정론·실행 근거 + +Temporal 측정량·grid·amplitude scale는 [검증 사례](verification-cases.md)가 고정한다. Exact solution이 있는 C01–C04/C08의 q/v/physical a는 마지막 두 refinement에서 error 감소 및 `log2(e_h/e_h2)>=1.8`, finest normalized max error `<=1e-4`를 요구한다. 이는 지배적인 2차항의 asymptotic window를 허용하는 내부 oracle 조건이고 외부 reference tolerance가 아니다. A에는 physical a의 2차 기준을 적용하지 않는다. 모든 error가 정확히 0이면 exact branch를 보고하고 관측 차수를 만들지 않는다. + +Temporal 검증에서만 다음 **read-only defect solve**를 추가한다. 마지막 Newton 선형계의 계산된 correction Δx와 잔차 `e_L=Khat_k Δx+r_k`에 대해 같은 행렬로 `Khat_k z_L=-e_L`을 푼다. 별도로 accepted endpoint에서 **새로 평가한** `Khat_* z_N=-r_*`를 푼다. z_L은 선형 solve의 결함, z_N은 비선형 endpoint의 남은 defect correction이며 서로 대체하지 않는다. 두 해의 N2-L backward check도 필요하다. 이 결과를 q/history에 적용하거나 Newton100 제한을 우회하는 polishing step으로 쓰지 않는다. + +채널/node별 `eta_x=||Δx_last||2+||z_L||2+||z_N||2`로 두고, `eta_w=gamma/(beta h)*eta_x`, `eta_b=eta_x/[(1+alpha)beta h²]`, `eta_B=eta_x/(beta h²)`로 전파한다. 각 항은 같은 scaled coordinate channel의 norm이며, D를 적용한 뒤 case별 q/v/a amplitude로 나누어 temporal error와 같은 단위·정규화로 비교한다. 최소 산술 해상도도 기록한다: `gamma32=32u/(1-32u)`, `eta_x,round=gamma32*(||x_*||2+||x_p||2)`, `eta_b,round=eta_x,round/[(1+alpha)beta h²]`; v/A도 해당 chain-rule 계수를 사용한다. 이 해상도는 position cancellation을 보여 주는 보수적 산술 지표이며 모든 element evaluation의 엄밀한 forward-error bound라는 주장은 하지 않는다. + +시간오차는 각 grid의 전체 시점에서 위 estimate/산술 지표의 최댓값보다 **100배 커야** order evidence로 사용한다. 미달이면 roundoff/iteration dominated로 **inconclusive**, 통과가 아니다. 정해진 coarse grids 또는 고정밀 독립 oracle로 다시 검토하며 실패한 판정값을 완화하지 않는다. 실제 shell self-convergence에는 C08의 별도 기준을 적용한다. 일반 production analysis/reference 비교에 이 두 추가 solve를 매 frame 수행하라는 요구는 아니다. + +동일 binary/runtime/thread/input 세 실행의 source/frame/diagnostic order와 sparse structure는 exact 일치, scaled scalar state 차이는 `1e-12` 이하를 요구한다. Failed trial은 committed q/v/a/A·history·frame count·기존 최종 HDF5를 바꾸지 않아야 한다. Reference fixed DIRECT는 100 increments이며 원 n*h schedule과 identity를 검사한다. + +수치값의 근거: `u=2^-53`, `gamma4096≈4.55e-13`에 대해 작은 dense/27점 항등식 `1e-12`, 그보다 100배 큰 assembled residual `1e-10`을 사용한다. `rcond>=2^-32`는 `gamma64/rcond≈3.06e-5`의 local worst-case 규모를 제한하며 보편 forward accuracy 보장이 아니다. FD `1e-6`/zero `1e-10`, objectivity `1e-10`, balance `1e-8`은 기존 [비선형 검토](../nonlinear-static-newton-raphson/numerical-review.md) 및 [eigen2 controls](../mitc4-shell-feast-eigen2/numerical-review-controls.md)의 해당 식을 현재 관성/좌표에 대조해 채택했다. 기존 정적 Newton 외력 상대 1%·typed absolute `1e-3`, FEAST controls, eigen 외부 RMS 예외는 채택하지 않았다. 새로운 runtime/solver 사례의 실행 상태는 아직 **NOT RUN**이며 숫자를 실행 통과처럼 보고하지 않는다. diff --git a/docs/hht-implicit-dynamics/numerical-research.md b/docs/hht-implicit-dynamics/numerical-research.md new file mode 100644 index 0000000..230db88 --- /dev/null +++ b/docs/hht-implicit-dynamics/numerical-research.md @@ -0,0 +1,152 @@ +# HHT 수치 제어·검증 설계 추가 근거 + +- Feature: `hht-implicit-dynamics` +- Owner: `research-agent` +- Date: `2026-09-22` +- Status: `Evidence and verification candidates; numerical thresholds not selected` +- 기준: [프로젝트 지침](../../AGENTS.md), [PRD](../PRD.md), [ADR](../ADR.md#adr-028), [기존 연구](research.md). + +[관성 정식화](formulation-inertia.md)와 [독립 probe](numerical-review-probes.md)를 읽고 +남은 질문만 조사했다. 기존 SDOF/moving-null 계산은 재실행하지 않았다. +**사실**은 열람한 일차 자료, **도출**은 명시한 수학적 계산, **후보**는 검증 제안이다. +수치 threshold·질량 채택·물리 범위·입출력 계약은 이 문서에서 결정하지 않는다. +현재 질량 선택은 사용자 결정으로 consistent mass가 확정되었다. Shell 속도 보간의 +운동에너지를 그대로 적분하며, lumped는 이번 HHT 구현의 지원 요구에서 제외한다. + +## 1. Physical acceleration, startup와 DAE 적용 경계 + +**사실 [N1].** Brüls–Arnold (2008) author postprint pp.3–4, Eq. (5)–(8)은 +물리적 `q_ddot`와 보조 acceleration-like variable을 구분하고 보조 초기값을 +물리적 초기 가속도와 같게 둔다. 이 문서의 `A`가 논문의 보조 `a`에 대응한다. +`alpha_m=0`, `alpha_f=-alpha`를 대입하면 현재 후보의 +`A_(n+1)=(1+alpha)a_(n+1)-alpha a_n`이다. 같은 논문은 variable M일 때 이를 +force-level 가중 평형과 동일시할 수 없다고 명시한다. 이 질문에는 기존 출처가 충분하다. + +**사실 [N1]: 정리의 가정은 정확히 구분한다.** Printed p.6, “Equivalent model”은 +`[[M,-L],[-h_qddot,I-h_y]]` block의 비특이성과 index-3 조건 +`Phi_q g*_lambda`의 비특이성을 사용한다. 여기의 `h`는 논문 Eq. (4)의 algebraic +output function이며 시간간격이 아니다. Controller/algebraic coupling을 제거한 +특수형에서는 첫 조건이 M의 비특이성으로 줄어든다. 논문 전체가 무조건 “M만 +비특이이면 충분”하다고 읽어서는 안 된다. “Initial value condition”의 초기 오차 조건과 +“Constrained case”의 정칙성 하에서 얻은 수렴 결과를 singular drilling에 자동 적용할 수 없다. + +**도출/수치 제어에의 함의.** 현재 shell의 `psi=N(q)^T g=0`는 논문의 holonomic +`Phi=0`와 동일한 모델임이 아직 증명되지 않았다. 따라서 다음을 분리해서 관찰한다. + +| 관찰 대상 | 필요한 검증 후보 | 남은 한계 | +| --- | --- | --- | +| 첫 시점 | `a0`는 우측 물리 평형·초기 대수조건으로 결정한 뒤 `A0=a0`; 첫 step과 이후 오차를 분리 | 작은 Newton residual로 불일관 초기값이 정당화되지 않음 | +| 물리 a와 A | 동일 시각 정확해에 대한 q/v/physical a 차수와 `A-a`를 별도 측정 | A를 물리 출력·각운동량 미분에 대입하면 안 됨 | +| Current null | 매 configuration에서 `MN`, `N^T I`, `psi`, 양의 부분공간 mass 및 `psi_q N` 정칙성 | 이전 시점 null, 고정 UR3, 전체 M의 역행렬로 대체 불가 | +| Hidden constraints | `psi_q v+psi_t`와 `psi_q a+psi_qq[v,v]+2psi_qt v+psi_tt`의 최대·종단 결함과 dt 수렴 | 위치 평형만 만족해도 두 결함이 0인 것은 아님 | +| 반례 coverage | 현재 움직이는 영공간 예제를 full redundant coordinates로 풀어 기존 reduced-angle probe와 대조 | 기존 probe는 위치조건을 정확히 강제했으므로 full-coordinate Newton 구현 검증이 아님 | +| 고정 간격 | 기존 후보에서 dt를 줄여 같은 물리 종단시각 비교 | adaptive dt·startup 재설정은 별도 scope이며 이번 근거로 추가하지 않음 | + +**도출.** Chart 내부의 smooth q에서 projector와 null basis는 같은 물리 영공간을 나타내야 +한다. Null vector의 부호·정규화 변경 자체와 rank 손실을 구분해야 한다. +`psi_q`에는 `DN`이 들어가므로 basis를 고정한 finite difference는 별도 식을 검증하게 된다. +투영 오차·hidden defect가 큰 결과를 발견하면 정확한 projection/condensation의 채택 여부는 +formulation owner에게 돌린다. 여기서는 추가 projection을 자동 제안·채택하지 않는다. + +## 2. Double precision, conditioning 및 미분 검증 + +**사실 [N2–N4].** IEEE double의 정상수 significand는 53 bits다. LAPACK DGECON은 +선택한 1/∞ norm의 reciprocal condition number **추정값**을 반환한다. DGESVX는 +equilibration 여부를 구별하고, factor pivot이 정확히 0인 경우와 reciprocal condition이 +machine precision보다 작은 경우를 서로 다른 INFO로 보고한다. 후자는 해·오차 추정을 +계속 계산할 수 있는 warning이다. 이 공개 library 동작은 FESA failure threshold를 +자동 승인하지 않으며 vendor routine 채택 요구도 아니다. + +**도출.** Binary64에서 1의 다음 표현수 간격 `eps=2^-52`, round-to-nearest 정상 범위의 +unit roundoff `u=2^-53`를 구분한다. 문헌의 “machine epsilon” 정의와 library 반환 의미를 +확인하지 않고 두 값을 혼용하면 안 된다. 작은 perturbation의 forward error는 +대략 `condition × backward error`로 증폭되므로 residual 하나로 해의 정확도를 보장하지 못한다. + +| 질문 | 근거에 따른 분리 및 검사 후보 | +| --- | --- | +| Physical mass-null인가? | `n=[0;J_l^-1 d]`, admissible free-space intersection 같은 분석적 운동학을 원본으로 삼는다. `MN`의 수치 결함과 양의 complement의 conditioning을 각각 검사한다. | +| 작은 양의 eigenvalue인가? | 단위·두께 때문에 작은 inertia를 rcond 하나로 drilling null로 바꾸지 않는다. 예상 추가 null/geometry 이상과 positive-space accuracy loss를 구분한다. | +| 어느 rcond인가? | Full singular M의 rcond는 positive inertia 품질을 설명하지 못한다. Positive-space mass, 초기화 algebraic block, effective Newton matrix를 각각 이름 붙이고 norm·scaling을 기록한다. | +| 단위를 바꾸면? | 병진·회전 단위가 섞인 원행렬 condition은 단위 선택에 의존한다. `q=D qhat`, `Rhat=D^T R/Escale`이면 `Khat=D^T K D/Escale`라는 virtual-work 일관 scaling을 대조 후보로 삼는다. 기준 scale의 수치 선택은 NR 소유다. | +| Library equilibration이면 충분한가? | Row/column equilibration은 solver conditioning 도구다. 물리적 nondimensionalization이나 예상 null proof와 동일하지 않다. 원·scaled 행렬 중 무엇을 추정했는지 보존한다. | + +**사실 [N5].** GSL의 central differentiation은 truncation과 roundoff를 함께 추정하여 +step을 고르며, 5-point와 3-point 차이를 오차 추정에 쓴다. 따라서 가장 작은 step을 +항상 정답으로 간주하지 않는 독립 derivative 검사 설계가 공식 수치 library의 접근과 부합한다. + +**도출/후보.** Dimensionless direction p에 대해 scaled 잔차의 중심차분과 `Khat p`를 +비교하면 일반적인 smooth·비상쇄 구간에서 오차형은 `C_t delta²+C_r u/delta`다. +최적 step의 `u^(1/3)` 크기는 이 가정에서의 점근적 안내일 뿐 공통 고정값이 아니다. +Zero derivative, 큰 residual cancellation, branch 접근에서는 계수와 error floor가 달라진다. +여러 delta를 sweep하여 감소구간·최저 plateau·roundoff 증가구간을 기록하고, 모든 sample이 +같은 허용 chart와 geometry 안에 있는지 확인한다. Unit-aware scale을 쓰고 임의 +`max(1, norm)`로 작은 문제를 덮지 않는다. + +Inertia partial derivative 검사에서는 q/v/physical a 중 하나만 바꾸고 나머지는 고정한다. +Effective tangent 검사에서는 같은 trial q에 연동되는 v, A, physical a를 모두 갱신한다. +두 검사는 서로 다른 미분이므로 섞지 않는다. 작은 dt의 acceleration 재구성에는 +위치 차분 오차가 `1/dt²`로 증폭되므로 시간 차수 측정에서 너무 미세한 dt의 plateau를 +수학적 차수 저하로 오인하지 않는다. 기준 oracle 정밀도와 nonlinear solve 오차도 분리한다. + +## 3. 3차 회전 미분과 운동량·에너지 case 후보 + +**사실 [N6].** Higham–Relton의 Theorem 3.5/Algorithm 3.6은 고차 Fréchet derivative를 +재귀 block matrix function의 우상단 block으로 구한다. Matrix exponential은 smoothness +조건을 만족한다. `P=[phi]_cross` 및 세 방향 행렬에 적용하면 24×24 exponential 하나의 +3×3 우상단 block이 3차 혼합 미분을 준다. + +**후보.** 이 방식에 `d0`를 곱한 결과를 Rodrigues/회전 Jacobian의 해석 미분과 비교하면 +동일 scalar coefficient 도함수를 재사용하는 자기검증을 피할 수 있다. “Exact derivative”는 +수학적 항등식을 뜻하며 부동소수점 exp 계산까지 오차가 없다는 뜻은 아니다. +Production 의존성 추가나 24×24 exp의 runtime 사용은 요구하지 않는다. + +**자체 도출: 영회전의 유한 다항 oracle.** `H_p=[p]_cross` 등으로 쓰면 + +\[ +D^3d(0)[p,s,r]=\frac16\sum_{\sigma\in S_3} +H_{\sigma(p)}H_{\sigma(s)}H_{\sigma(r)}d_0. +\] + +이는 exponential Taylor series의 cubic 항에서 바로 얻는다. 방향 permutation 대칭, +각 방향의 선형성, 동일 방향 3회 `H_p³ d0`를 함께 검사할 수 있다. 비영 phi의 +block-exponential oracle와 `D²d`의 방향 중심차분을 추가해 영점 특례만 맞는 오류를 검출한다. +Chart 경계 안쪽, 작은 phi, 서로 비평행한 방향을 포함하되 경계를 넘어 차분하지 않는다. + +| ID | 실행 가능한 검증 설계 후보 | 판정에 필요한 한계 | +| --- | --- | --- | +| NR-E01 | 임의 유효 shell 상태에서 spatial rigid velocity를 부여하여 kinetic energy·총 P·총 L을 독립 volume integration과 대조 | Rigid velocity는 순간 운동학 검사다. 회전하는 탄성 shell이 자유운동에서도 강체를 유지한다는 가정은 하지 않음. Consistent/lumped 각각 자신의 measure를 사용 | +| NR-E02 | `I(q,v,a)`를 q/v/a별 미분하고 전체 effective tangent도 별도 directional check | 3차 director 항과 `D(H^T)` 누락을 분리 검출 | +| NR-E03 | 작은 smooth 회전 구간의 자유 shell에 일관된 initial velocity를 부여하거나 짧은 smooth loading 후 unload, P/L/총에너지·work와 dt 수렴 관찰 | 초기조건/시간하중이 CLI 범위 밖이면 kernel/procedure verification 후보로만 유지. 승인 input subset 확대 금지 | +| NR-E04 | 같은 spatial problem에 alpha=0 및 승인된 음수 alpha를 적용하고 저주파 해상 case와 높은 modal content case를 분리 | 일반 비선형 HHT의 총에너지 매 step 감소 또는 정확한 각운동량 보존을 무조건 요구하지 않음 | +| NR-E05 | 문헌 free-flying single-layer plate를 후속 benchmark 후보로 검토 | 아래 N7은 large rotation/energy 사례의 출처다. 현재 fixed chart를 넘어서는 tumble·장시간 운동은 그대로 채택할 수 없음 | + +**사실 [N7].** Vu-Quoc–Tan (2003), §5.1.2/Figs.9–12는 free-flying single-layer plate의 +큰 전체 운동과 P/L/에너지 및 dissipative/nondissipative integrator 결과를 비교한다. +해당 요소는 rotation DOF가 없는 solid-shell이므로 그 mass·응력 모델과 여기의 director shell은 +다르다. Author-uploaded 본문의 관련 절을 확인했으나 그림의 하중 위치·수치 history를 +독립 추출하지 않았다. 따라서 현 단계에서는 case family 근거이며 완성된 reference artifact가 아니다. + +**추론.** 현재 shell 운동학에서는 `Pdot/Ldot`에 physical a를 사용하고 numerical drill +potential을 포함한 에너지·내력 일을 일관되게 대조해야 한다. N7의 보존형 integrator가 +보이는 exact conservation을 HHT의 요구조건으로 옮길 수 없다. [기존 연구 S8](research.md#sources)의 +비선형 Newmark 보존 한계 근거는 그대로 충분하다. 추가 shell 연구는 기하·하중·시간 구간이 +승인된 chart와 입력 범위 안에 있는지 결정한 뒤에만 구체화할 필요가 있다. + +## 출처·열람 한계와 인계 + +`accessed: 2026-09-22`. 문서 정적 검사만 수행하며 solver/build/test/reference 생성은 수행하지 않았다. + +| ID | Stable citation / locator | 열람 범위 | +| --- | --- | --- | +| N1 | Brüls, O.; Arnold, M. (2008), *The generalized-alpha scheme as a linear multistep integrator: Toward a general mechatronic simulator*, J. Computational and Nonlinear Dynamics 3, 041007, [DOI:10.1115/1.2960475](https://doi.org/10.1115/1.2960475), [author PDF](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf) | Printed pp.3–6 Eq.5–8/Equivalent model, Initial value condition/Constrained case. Shell singular-mass 정리 아님 | +| N2 | Oracle, *Numerical Computation Guide*, [§2.2.3 Double Format](https://docs.oracle.com/cd/E37069_01/html/E39019/z4000ac019269.html) | 53-bit significand 설명. IEEE 표준 전체 적합성 시험을 수행한 것은 아님 | +| N3 | LAPACK, [DGECON](https://netlib.org/lapack/explore-html/d4/daf/group__gecon_ga4f9b830e19e12c7f082ddb497a57af18.html), Purpose/NORM/ANORM/RCOND | Normwise reciprocal condition estimate | +| N4 | LAPACK, [DGESVX](https://netlib.org/lapack/explore-html/d5/dbe/group__gesvx_ga82173a93234afc15d70b64233b3e5bc8.html), Description steps 1–6/INFO | Equilibration, singular pivot와 working-precision warning 구별 | +| N5 | GNU Scientific Library 2.8, [Numerical Differentiation](https://www.gnu.org/software/gsl/doc/html/diff.html), `gsl_deriv_central` | Truncation/roundoff 기반 step/error estimate. FESA tolerance 원본 아님 | +| N6 | Higham, N.J.; Relton, S.D. (2014), *Higher Order Fréchet Derivatives of Matrix Functions and the Level-2 Condition Number*, SIAM J. Matrix Anal. Appl. 35, 1019–1037, [DOI:10.1137/130945259](https://doi.org/10.1137/130945259), [author PDF](https://eprints.maths.manchester.ac.uk/2066/1/paper.pdf), Theorem 3.5/Algorithm 3.6 | 2013 author version printed pp.5–6; oracle identity와 계산법 | +| N7 | Vu-Quoc, L.; Tan, X.G. (2003), *Optimal solid shells for non-linear analyses of multilayer composites. II. Dynamics*, CMAME 192, 1017–1059, [DOI:10.1016/S0045-7825(02)00336-5](https://doi.org/10.1016/S0045-7825(02)00336-5), [author-uploaded text](https://www.researchgate.net/publication/222546721_Optimal_solid_shells_for_non-linear_analyses_of_multilayer_composites_II_Dynamics), §5.1.2/Figs.9–12 | Relevant text 확인; 하중 figure·곡선 digitization 미수행, exact benchmark reproduction 미정 | + +Numerical Review에 인계할 결정은 scaling·conditioning/hidden-defect 판정·FD sweep과 +시간 수렴 구간·case 최소 조합이다. Formulation에 남는 사항은 singular-mass regularity의 +실제 shell 적용과 full-coordinate hidden constraint 검증이다. 사용자 결정에 따라 이번 +HHT 구현은 consistent mass를 사용한다. 위 lumped 비교 근거는 연구 이력이며 현재 지원 요구가 아니다. diff --git a/docs/hht-implicit-dynamics/numerical-review-probes.md b/docs/hht-implicit-dynamics/numerical-review-probes.md new file mode 100644 index 0000000..6dce1a2 --- /dev/null +++ b/docs/hht-implicit-dynamics/numerical-review-probes.md @@ -0,0 +1,162 @@ +# HHT 독립 시간정확도 검산 + +[검토 보고서](numerical-review.md)의 상세 evidence다. 2026-09-22 Python 표준 라이브러리만 사용했다. FESA solver, Abaqus 또는 reference 비교 실행이 아니다. 아래 오차는 승인 tolerance를 새로 정하는 값이 아니라 정확해와의 수렴 관찰값이다. + +## 모델과 결과 + +1. Smooth SDOF: `T=.5(1+q²)v²`, `c=qv²`, `q_exact=.4+sin(t)`, `f=(1+q_exact²)a_exact+q_exact v_exact²`, `0≤t≤1`. 초기 `q=.4,v=1,a=0`. 입력은 일관된 무차원 검산 단위다. `alpha=-.05,beta=.275625,gamma=.55`. +2. Moving-null: `x(q)=q/|q|`, `T=.5|Dx(q)v|²`, `U=.5(|q|-1)²`, `q∈R²\{0}`. `r=|q|`, `d=q/r`, `M=(I-ddᵀ)/r²`, `c=-2(q·v)Mv/r²`, `f_int=(r-1)d`. `ker M=span(q)`가 회전한다. 정확해 `q=(cos t,sin t)`, `v=(-sin t,cos t)`, `a=-q`; 초기 위치·속도·가속도 대수조건은 일관된다. Endpoint 반경 대수조건 `r=1`은 각도 parameterization으로 정확하게 강제하며 tangent 평형만 Newton으로 푼다. 이 모델은 shell/drilling 구현이나 일반 singular-DAE 정리의 대체물이 아니다. + +| SDOF 방법 | `h=1/512` U 오차 | `h=1/1024` U 오차 | 관측 차수 | +| --- | ---: | ---: | ---: | +| Naive `M1 A1+(1+alpha)g1-alpha g0=0` | 1.0272820368e-5 | 5.2200068232e-6 | 0.976709 | +| Weighted mass `((1+alpha)M1-alpha M0)A1+(1+alpha)g1-alpha g0=0` | 3.2914014847e-7 | 8.2299941973e-8 | 1.999739 | +| Endpoint physical `M1 a1+g1=0`, `A1=(1+alpha)a1-alpha a0` | 2.7561411953e-7 | 6.8902116013e-8 | 2.000030 | + +`g=c-f`. 마지막 방법의 physical acceleration 오차는 `2.3176116737e-7 → 5.7545108678e-8`로 줄지만 algorithmic `A` 오차는 `5.2612378749e-5 → 2.6344436092e-5`다. 따라서 `A`를 endpoint physical acceleration으로 보고하면 안 된다. + +Moving-null 마지막 방법은 같은 두 h에서 U 오차 `4.0606413493e-7 → 1.0152120526e-7`, V 오차 `4.1505360817e-7 → 1.0376522277e-7`, physical a 오차 `4.1511310820e-7 → 1.0262747089e-7`를 보였다. `max|q·v|`는 `8.4282e-11 → 1.0594e-11`, `max|q·a+v²|`는 `7.9368e-7 → 1.9445e-7`다. Endpoint 위치 대수조건을 지켜도 미분된 속도·가속도 조건이 매 frame 정확히 0이 되는 것은 아니다. 이 결함의 시간수렴과 초기 일관성을 별도로 검증해야 한다. + +## 독립 수학 대조 + +Smooth 물리식 `M(q)a+g(q,v,t)=0`에서 Newmark 내부 가속도의 정확해 대응은 `A=a+alpha*h*adot+O(h²)`다. Naive 식의 일차 결함은 + +`alpha*h*(M*adot+gdot)=-alpha*h*Mdot*a`. + +따라서 `Mdot*a≠0`인 일반 문제에서 상수 M 증명을 그대로 사용할 수 없다. Weighted mass는 이 항을 상쇄하지만 `alpha∈(-1,0)`이면 서로 다른 두 시점의 양의 준정부호 질량의 양의 결합이다. 그 영공간은 두 영공간의 교집합이므로, 영공간이 회전하면 rank가 늘어날 수 있다. 단순 SDOF의 2차 결과가 이를 정당화하지 않는다. + +Physical/algorithmic 가속도 분리는 [Brüls–Arnold 원문](https://orbi.uliege.be/bitstream/2268/27088/1/pa_BrulsArnold2008_TheGeneralizedAlphaSchemeAsALinearMultistepIntegrator.pdf)의 Eq. (5), (7), (8)에 `alpha_m=0,alpha_f=-alpha`를 대입한 것이다. 원문의 p.6 equivalent-model/index-3 가정은 비특이 mass 및 제약 정칙성을 사용한다. Singular drilling mass에 그 증명을 직접 적용하지 않는다. 이 문서의 제조해와 결함 계산은 독립 검산이다. + +## 재현 코드와 명령 + +Workspace root에서 다음 명령으로 아래 Python block만 실행한다. 외부 package, 파일 생성 또는 solver 실행은 없다. + +```powershell +python -c 'from pathlib import Path; s=Path("docs/hht-implicit-dynamics/numerical-review-probes.md").read_text(encoding="utf-8"); exec(s.split(chr(96)*3+"python\n",1)[1].split(chr(96)*3,1)[0])' +``` + +```python +import math + +alpha = -.05 +beta = (1-alpha)**2/4 +gamma = .5-alpha + +def exact(t): + return .4+math.sin(t), math.cos(t), -math.sin(t) + +def force(t): + q, v, a = exact(t) + return (1+q*q)*a+q*v*v + +def scalar(h, variant): + q, v, ap = exact(0) + A = ap + for n in range(round(1/h)): + t = (n+1)*h + m0 = 1+q*q + g0 = q*v*v-force(n*h) + oldap = ap + pred = q+h*v+h*h*(.5-beta)*A + vpred = v+h*(1-gamma)*A + def calc(x): + an = (x-pred)/(beta*h*h) + vn = vpred+gamma*h*an + m = 1+x*x + if variant == 'physical': + pn = (an+alpha*oldap)/(1+alpha) + r = m*pn+x*vn*vn-force(t) + dr = 2*x*pn+m/((1+alpha)*beta*h*h)+vn*vn+2*x*vn*gamma/(beta*h) + else: + me = m if variant == 'naive' else (1+alpha)*m-alpha*m0 + dm = 2*x if variant == 'naive' else (1+alpha)*2*x + r = me*an+(1+alpha)*(x*vn*vn-force(t))-alpha*g0 + dr = dm*an+me/(beta*h*h)+(1+alpha)*(vn*vn+2*x*vn*gamma/(beta*h)) + pn = an + return r, dr, vn, an, pn + x = pred+beta*h*h*A + for _ in range(20): + r, dr, *_ = calc(x) + dx = -r/dr + x += dx + if abs(dx) < 1e-14: + break + else: + raise RuntimeError('scalar Newton failed') + q = x + _, _, v, A, ap = calc(q) + qe, ve, ae = exact(1) + return abs(q-qe), abs(v-ve), abs(ap-ae), abs(A-ae) + +def moving_null(h): + q, v, ap = [1.,0.], [0.,1.], [-1.,0.] + A = ap[:] + theta = 0. + max_cv = max_ca = 0. + dot = lambda x,y: sum(x[i]*y[i] for i in range(2)) + for _ in range(round(1/h)): + pred = [q[i]+h*v[i]+h*h*(.5-beta)*A[i] for i in range(2)] + vp = [v[i]+h*(1-gamma)*A[i] for i in range(2)] + old = ap[:] + theta += h + def calc(angle): + d = [math.cos(angle),math.sin(angle)] + tangent = [-d[1],d[0]] + an = [(d[i]-pred[i])/(beta*h*h) for i in range(2)] + vn = [vp[i]+gamma*h*an[i] for i in range(2)] + pn = [(an[i]+alpha*old[i])/(1+alpha) for i in range(2)] + dv, tv = dot(d,vn), dot(tangent,vn) + residual = dot(pn,tangent)-2*dv*tv + derivative = 1/((1+alpha)*beta*h*h)-dot(pn,d)-2*(tv*tv+dv*(-dv+gamma/(beta*h))) + return residual, derivative, d, vn, pn, an + for _ in range(15): + r, dr, *_ = calc(theta) + dx = -r/dr + theta += dx + if abs(dx) < 1e-14: + break + else: + raise RuntimeError('moving-null Newton failed') + _, _, q, v, ap, A = calc(theta) + max_cv = max(max_cv,abs(dot(q,v))) + max_ca = max(max_ca,abs(dot(q,ap)+dot(v,v))) + qe = [math.cos(1),math.sin(1)] + ve = [-math.sin(1),math.cos(1)] + ae = [-x for x in qe] + errors = [math.hypot(*[x[i]-y[i] for i in range(2)]) for x,y in ((q,qe),(v,ve),(ap,ae))] + return *errors, max_cv, max_ca + +for variant in ('naive','weighted_mass','physical','moving_null'): + print(variant) + previous = None + for k in range(5,11): + h = 2**(-k) + errors = moving_null(h) if variant == 'moving_null' else scalar(h,variant) + order = math.log(previous/errors[0],2) if previous else float('nan') + print(f'h={h:.9g} errors='+','.join(f'{x:.10e}' for x in errors)+f' pU={order:.6f}') + previous = errors[0] +``` + +Newton correction cutoff `1e-14`는 이 독립 double-precision probe의 반복 종료값이다. Solver convergence 또는 reference acceptance tolerance로 채택하지 않는다. 미세 h의 acceleration 오차에는 displacement 차분의 roundoff가 포함되므로 표의 관측값을 정확한 수학 상수로 해석하지 않는다. + +## 후속 full Cartesian moving-null 검산 + +추가 Numerical 작업에서 메인이 [full_coordinate_probe.py](../../.harness/build/hht-numerical/full_coordinate_probe.py)를 작성했고 수치 검토자가 residual/Jacobian 식을 읽고 아래 명령을 재실행했다. Exit `0`, 2026-09-22. 이 script는 ignored build evidence이며 FESA source가 아니다. 위 angle-reduced probe와 달리 q의 두 Cartesian 성분을 직접 Newton으로 풀며 radial projection·normalization을 하지 않는다. 같은 `x=q/|q|` 모델의 exact 2×2 Jacobian을 사용한다. + +```powershell +python .harness/build/hht-numerical/full_coordinate_probe.py +``` + +오차는 각 h의 **모든 시점**에서 q/v/physical a의 2-norm 최대다. 앞 절의 종단 오차와 다른 측정량이다. + +| h | max q error | max v error | max physical a error | 직전 h 대비 q/v/a 차수 | +| --- | ---: | ---: | ---: | --- | +| 1/128 | 6.4948801161e-6 | 6.6395874231e-6 | 1.2690547915e-5 | 1.99932 / 1.99962 / 2.01024 | +| 1/256 | 1.6240799623e-6 | 1.6601117645e-6 | 3.1667175260e-6 | 1.99968 / 1.99981 / 2.00269 | +| 1/512 | 4.0606299832e-7 | 4.1505196280e-7 | 7.9259171465e-7 | 1.99985 / 1.99992 / 1.99834 | + +마지막 grid의 최대 `||q|-1|=2.22045e-16`, `|q·v|=8.39836e-11`, `|q·a+v²|=7.91118e-7`. 상쇄 전 `Ma`의 entry products·c·f_int를 사용한 최대 normwise residual은 `7.94824e-11`, Newton 최대 3 corrections, 마지막 correction 최대 `4.00220e-16`이었다. 이 probe의 종료는 correction `<=1e-13`이며 N2의 모든 AND gate를 구현한 solver가 아니다. 위 관찰을 N2–N7 전체 통과나 실제 shell 통과로 확대하지 않는다. + +별도 trial `q=(.85,.52)`, predictor `(.83,.50)`, velocity predictor `(-.35,.71)`, old physical a `(-.7,-.3)`, h=.1, direction `(.6,-.8)`의 exact effective Jacobian을 중앙차분했다. N5의 j=8,10,12,14,16,18,20 순서에서 상대 오차는 `[3.13656e-5,1.96035e-6,1.22522e-7,7.65779e-9,4.77929e-10,3.25134e-11,1.04250e-11]`이다. 연속한 내부 sweep 점 j=12/14 등이 `1e-6`을 만족한다. 이는 이 2D trial/direction에 대한 미분 evidence이며 shell D3 oracle를 대체하지 않는다. + +Full-coordinate 결과는 새로운 [C04b](verification-cases.md)의 moving-null 구조와 max-time 측정을 지지한다. C04b의 전체 prescribed grids, 초기화 실패 mutation, regularity·scale·반복오차 분리 등 모든 판정을 실행한 것은 아니다. 더 미세한 h에서 residual이 증가할 수 있으므로 h=1/1024까지 같은 residual 기준으로 통과했다는 주장을 하지 않는다. diff --git a/docs/hht-implicit-dynamics/numerical-review.md b/docs/hht-implicit-dynamics/numerical-review.md new file mode 100644 index 0000000..a059ed9 --- /dev/null +++ b/docs/hht-implicit-dynamics/numerical-review.md @@ -0,0 +1,64 @@ +# HHT 독립 수치 검토·검증 계약 + +- Feature: `hht-implicit-dynamics`; owner: `numerical-review-agent`; date: `2026-09-22`. +- 상태: `numerical-contract-defined; io-binding-and-implementation-verification-pending`. 기존 정식화 검토에 후속 Numerical gate의 내부 판정·검증 사례·reference mapping 검토를 추가했다. Solver 기능 통과 또는 I/O schema 승인이 아니다. +- 기준: [프로젝트 지침](../../AGENTS.md), [Feature PRD](PRD.md), [요구사항](requirements.md), [연구 근거](research.md), [공통 tolerance](../TOLERANCE.md). +- 상세: [내부 수치 기준](numerical-controls.md), [검증 사례](verification-cases.md), [시간정확도 검산](numerical-review-probes.md), [추가 일차 근거](numerical-research.md), [reference mapping](reference-model.md). Solver source/tests/reference/상위 계약은 수정하지 않았다. + +## 결론 + +[주 정식화](formulation.md)와 [관성 상세](formulation-inertia.md)를 독립 대조했다. **현재 physical/algorithmic 가속도 분리 후보의 표시된 수식에서 새 오류는 발견하지 않았다.** 사용자는 [HHT-ADR-007](ADR.md)에서 consistent kinetic mass를 선택했다. 내부 수치값과 반증 가능한 사례를 구체화했으며 Numerical 계약은 후속 I/O/구현에 인계할 수 있다. 실제 finite-rotation shell의 drilling DAE·시간수렴·reference 검증은 구현 후 필수 evidence다. 이 미실행을 수식 모순이나 코드 작성 금지로 해석하지 않는다. + +| 확정한 내부 판정 | 핵심 기준 | +| --- | --- | +| Newton·선형계 | Force/moment별 상쇄 전 기여 scale의 residual `1e-10`, correction 및 h-scaled v/a/A 각각 `1e-10`의 AND; linear backward error `1e-12`; 최대100 corrections, fixed DIRECT/no cutback | +| Mass·DAE | Analytic `Jl^-1 d`의 free-space intersection으로 expected rank; positive-space Cholesky, scaled local/대수 rcond `>=2^-32`, null/Gram/moment `1e-12`; 초기 hidden consistency `1e-10` | +| 미분·불변식 | D3/관성/전체 tangent step sweep의 nonzero error `1e-6`, exact zero action `1e-10`; nonlinear objectivity `1e-10`, force/moment balance `1e-8`; dynamic Keff 대칭 가정 없음 | +| 시간·coverage | 지정 smooth oracle의 q/v/physical a 최대오차와 마지막 두 refinement 차수 `>=1.8`; A는 별도량. Full Cartesian moving-null, coupled massless startup, finite-rotation shell, energy, rollback, 전체 reference 사례 | + +숫자의 norm·zero branch·단위·실행 상태는 [numerical-controls](numerical-controls.md)가 원본이다. 위 요약의 값만 떼어 다른 metric에 적용하지 않는다. 기존 정적 외력 상대 1%/absolute `1e-3`나 eigen 외부 RMS 예외를 동적 Newton에 승계하지 않았다. + +## 위험과 수정 반영 + +| ID | 판단 | 근거·필요 조치 | +| --- | --- | --- | +| HHT-NR-001 | Naive variable-mass HHT의 일반 2차 주장 **FAIL; 현 후보에서 배제 완료** | 여기서 `g=c+f_int-F`. `M(q1)A1+(1+alpha)g1-alpha*g0`에 `-alpha*h*Mdot*a` 결함. Smooth SDOF의 마지막 refinement 차수 `0.976709`. 관성 상세 §5.1이 이를 배제한다. | +| HHT-NR-002 | 두 시점 mass 가중을 moving-null 해결책으로 채택 **FAIL; 현 후보에서 배제 완료** | 양의 가중 합의 nullspace는 교집합이다. 회전하는 drill null이 사라지면 인공 관성이 생긴다. SDOF 2차 관찰은 shell 채택 근거가 아니다. 관성 상세 §5.3이 이를 채택하지 않는다. | +| HHT-NR-003 | Endpoint physical `a`와 auxiliary `A` 분리 **SUPPORTED CANDIDATE** | `A1=(1+alpha)a1-alpha*a0`, Newmark에는 A, endpoint balance에는 a. Smooth variable-mass와 한 moving-null 제조해에서 U/V/a의 2차 경향 확인. 일반 shell DAE 수렴 증명·구현 통과는 아니다. | +| HHT-NR-004 | `A`를 physical acceleration으로 출력 **FAIL; 현 후보에서 구분 완료** | 위 후보에서도 A의 endpoint 오차는 1차다. 주 정식화 §§4–6이 a/A와 endpoint reaction을 구분한다. 영속 결과 선택은 I/O owner의 후속 계약이다. | + +## 수식 대조 결과 + +| 검토 항목 | 결과와 evidence | +| --- | --- | +| DOF·좌표·회전 | `[u,phi]`, total fixed-global chart `|phi|