From f0dc42a280e710e2423e7390b7a331b61f67de9c Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Tue, 22 Sep 2026 21:38:42 +0900 Subject: [PATCH] feat: implement consistent-mass HHT shell dynamics --- HANDOFF.md | 2 + docs/hht-implicit-dynamics/ADR.md | 10 + .../implementation-numerical-review.md | 49 + docs/hht-implicit-dynamics/io-contract.md | 2 +- .../numerical-controls.md | 2 +- .../physics-evaluation.md | 47 + .../fesa/analysis/dynamic_analysis_state.h | 113 ++ include/fesa/analysis/hht_dynamic_procedure.h | 29 + include/fesa/analysis/hht_time_integration.h | 50 + include/fesa/analysis/newton_iteration_log.h | 27 + include/fesa/analysis/prepared_analysis.h | 17 +- .../fesa/elements/nonlinear_mitc4_inertia.h | 106 ++ include/fesa/elements/nonlinear_mitc4_shell.h | 17 +- .../io/hdf5/dynamic_hdf5_results_writer.h | 20 + include/fesa/results/dynamic_results_writer.h | 28 + scripts/verification/compare_hht_dynamic.py | 406 +++++++ src/fesa/CMakeLists.txt | 7 + src/fesa/analysis/dynamic_analysis_state.cpp | 406 +++++++ src/fesa/analysis/hht_dynamic_algebra.cpp | 732 ++++++++++++ src/fesa/analysis/hht_dynamic_assembly.cpp | 483 ++++++++ src/fesa/analysis/hht_dynamic_internal.h | 102 ++ src/fesa/analysis/hht_dynamic_procedure.cpp | 322 +++++ src/fesa/analysis/hht_time_integration.cpp | 97 ++ src/fesa/analysis/newton_iteration_log.cpp | 83 ++ src/fesa/analysis/prepared_analysis.cpp | 38 +- src/fesa/app/fesa_application.cpp | 13 +- src/fesa/elements/nonlinear_mitc4_inertia.cpp | 871 ++++++++++++++ src/fesa/elements/nonlinear_mitc4_shell.cpp | 427 +++++-- .../io/hdf5/dynamic_hdf5_results_writer.cpp | 1017 ++++++++++++++++ tests/CMakeLists.txt | 8 + tests/test_compare_hht_dynamic.py | 302 +++++ .../analysis/dynamic_analysis_state_test.cpp | 305 +++++ .../analysis/hht_dynamic_algebra_test.cpp | 330 ++++++ .../analysis/hht_dynamic_procedure_test.cpp | 125 ++ .../analysis/hht_energy_verification_test.cpp | 223 ++++ .../hht_temporal_verification_test.cpp | 507 ++++++++ .../analysis/hht_time_integration_test.cpp | 97 ++ .../analysis/newton_iteration_log_test.cpp | 150 +++ .../elements/nonlinear_mitc4_inertia_test.cpp | 1051 +++++++++++++++++ .../nonlinear_mitc4_shell_tangent_test.cpp | 162 +++ .../hdf5/dynamic_hdf5_results_writer_test.cpp | 636 ++++++++++ 41 files changed, 9294 insertions(+), 125 deletions(-) create mode 100644 docs/hht-implicit-dynamics/implementation-numerical-review.md create mode 100644 docs/hht-implicit-dynamics/physics-evaluation.md create mode 100644 include/fesa/analysis/dynamic_analysis_state.h create mode 100644 include/fesa/analysis/hht_dynamic_procedure.h create mode 100644 include/fesa/analysis/hht_time_integration.h create mode 100644 include/fesa/elements/nonlinear_mitc4_inertia.h create mode 100644 include/fesa/io/hdf5/dynamic_hdf5_results_writer.h create mode 100644 include/fesa/results/dynamic_results_writer.h create mode 100644 scripts/verification/compare_hht_dynamic.py create mode 100644 src/fesa/analysis/dynamic_analysis_state.cpp create mode 100644 src/fesa/analysis/hht_dynamic_algebra.cpp create mode 100644 src/fesa/analysis/hht_dynamic_assembly.cpp create mode 100644 src/fesa/analysis/hht_dynamic_internal.h create mode 100644 src/fesa/analysis/hht_dynamic_procedure.cpp create mode 100644 src/fesa/analysis/hht_time_integration.cpp create mode 100644 src/fesa/elements/nonlinear_mitc4_inertia.cpp create mode 100644 src/fesa/io/hdf5/dynamic_hdf5_results_writer.cpp create mode 100644 tests/test_compare_hht_dynamic.py create mode 100644 tests/unit/analysis/dynamic_analysis_state_test.cpp create mode 100644 tests/unit/analysis/hht_dynamic_algebra_test.cpp create mode 100644 tests/unit/analysis/hht_dynamic_procedure_test.cpp create mode 100644 tests/unit/analysis/hht_energy_verification_test.cpp create mode 100644 tests/unit/analysis/hht_temporal_verification_test.cpp create mode 100644 tests/unit/analysis/hht_time_integration_test.cpp create mode 100644 tests/unit/elements/nonlinear_mitc4_inertia_test.cpp create mode 100644 tests/unit/io/hdf5/dynamic_hdf5_results_writer_test.cpp diff --git a/HANDOFF.md b/HANDOFF.md index bf8b950..027ab1b 100644 --- a/HANDOFF.md +++ b/HANDOFF.md @@ -214,3 +214,5 @@ PowerShell helper의 `-B`는 named parameter로 잘못 bind될 수 있어 위 ar - 2026-09-22 初回記録: 本文のsnapshot作成。Reference16/100受理、最新Debug focused実行中。Comparator41testsのmain再実行PASSを確認。今後も作業を継続する。 - 2026-09-22 후속 확인: 독립 physics reviewer가 `physics-evaluation.md` 초안을 작성했다. 최신 Debug focused는 68 tests 중67 PASS/1 C04a INCONCLUSIVE SKIP으로 종료했다. Main은 다음 명령에서 종료코드와 summary를 직접 확인한다. C08 및 최종 reference 검증은 계속 PENDING이다. + +- 2026-09-22 후속 작업: main이 Debug focused session81971 종료코드0을 직접 확인했다. 68 tests/67 PASS/1 INCONCLUSIVE SKIP, 총117.526s. Comparator41tests도 main 확인 완료. 구현 묶음을 로컬 커밋한 뒤 C08 파일을 등록한다. 신규 private test는 아직 이 Debug 실행에 포함되지 않았다. diff --git a/docs/hht-implicit-dynamics/ADR.md b/docs/hht-implicit-dynamics/ADR.md index 4a22476..7c56df2 100644 --- a/docs/hht-implicit-dynamics/ADR.md +++ b/docs/hht-implicit-dynamics/ADR.md @@ -78,3 +78,13 @@ Feature 범위·상태는 [PRD](PRD.md), 상세 acceptance intent는 [requiremen - 승인 근거: 사용자의 “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는 포함하지 않는다. 실패·미실행 검증을 통과로 표시하지 않는다. + + +## HHT-ADR-010 — Newton 시작 형상과 영변형 산술 + +- 상태: `accepted-implementation-detail`; reference 통과를 의미하지 않는다. +- Owner: 메인 assistant; 독립 수치 검토는 numerical-review-agent. +- 근거: 승인된 전체 구현 범위 안의 수치 시작값 선택. `.harness/build/hht-numerical/zero-trial-red.log`와 original reference 실행에서 먼저 실패를 재현했다. +- 결정: 모든 increment의 Newton 초기 좌표를 `q_n`으로 둔다. Constant-A 외삽은 지정 DIRECT h=.1에서 최초 trial geometry를 무효화했으므로 N2의 시작값을 명시적으로 갱신한다. HHT predictor/rate/접선, alpha/beta/gamma, 고정 dt, Newton100 및 모든 tolerance는 유지한다. 유효한 시작점에서 나온 correction도 chart/geometry 밖이면 즉시 실패하며 retry·clamp·cutback은 없다. +- 추가 수정: 초기 형상에서 차감한 큰 covariant metric 대신 `delta g`로 Green strain과 `F=I+delta g⊗Gdual`을 계산한다. 수학적으로 같은 식이며 영변형에서 생긴 artificial drilling force를 제거한다. 강성·정적 결과 영향은 기존 regression으로 확인한다. +- 검증: decimal 좌표 영변형 exact-zero 회귀, finite/warped scalar drilling constraint·gradient·Hessian, 작은각 D3 threshold 검사가 통과했다. 상세 실행/물리 acceptance는 최종 implementation evidence에서 구분한다. diff --git a/docs/hht-implicit-dynamics/implementation-numerical-review.md b/docs/hht-implicit-dynamics/implementation-numerical-review.md new file mode 100644 index 0000000..c9a91ec --- /dev/null +++ b/docs/hht-implicit-dynamics/implementation-numerical-review.md @@ -0,0 +1,49 @@ +# HHT 관성 구현 독립 수치 검토 + +- Owner: `numerical-review-agent`; date: `2026-09-22`. +- 상태: `source-delta-reviewed; runtime-verification-pending`. 수식/API 중간 검토이며 전체 구현 gate 통과 판정이 아니다. +- 범위: [관성 API](../../include/fesa/elements/nonlinear_mitc4_inertia.h), [관성 구현](../../src/fesa/elements/nonlinear_mitc4_inertia.cpp), [관성 tests](../../tests/unit/elements/nonlinear_mitc4_inertia_test.cpp), [비선형 내력](../../src/fesa/elements/nonlinear_mitc4_shell.cpp)의 directional third-derivative, [초기화·null algebra](../../src/fesa/analysis/hht_dynamic_algebra.cpp), [동적 조립](../../src/fesa/analysis/hht_dynamic_assembly.cpp), [HHT procedure](../../src/fesa/analysis/hht_dynamic_procedure.cpp). +- 기준: [정식화](formulation.md), [관성 상세](formulation-inertia.md), [N1–N7](numerical-controls.md), [C01–C12](verification-cases.md). Source/tests/상위 계약은 읽기 전용이다. + +## 첫 검토 findings + +아래는 작업 중 최초로 관찰한 내용이다. Owner 수정 뒤의 대조·실행 evidence는 별도로 갱신한다. + +| ID | 심각도·상태 | 관찰·수정과 필요한 검증 | +| --- | --- | --- | +| HHT-IMP-NR-001 | **P1; source corrected, execution pending** | `Inverse(JetMatrix)`의 cyclic minor에 이미 포함된 cofactor 부호를 다시 적용하던 오류가 제거됐다. 새 test가 독립 `Jl^-T`와 현재 null, finite rotation의 `M n=0`을 대조한다. 실행 evidence 확인 전이다. | +| HHT-IMP-NR-002 | **P2; corrected, focused regression passed** | `convective_force=force-M*a` cancellation을 직접 `∫rho Hᵀ B_x[v,v]dV0` 누적으로 수정했다. 같은 q/v에서 a를 `1e15`배 변경해도 c가 동일한지 검사하는 test의 통과 로그를 확인했다. | +| HHT-IMP-NR-003 | **P2; corrected, focused sweep passed** | 최초 raw mixed-unit direction/`1e-5` test를 D scaling/단위 direction, N5 `1e-6`·연속한 두 내부 sweep 점으로 수정했다. Position/velocity/acceleration partial 및 두 h의 dependent-rate effective inertia tangent, pure translation/rotation/mixed 방향과 해당 실행 통과를 확인했다. | +| HHT-IMP-NR-004 | **미완료 coverage; 실행 대기** | 새 N4 guard는 analytic null의 20차원 직교 여공간에서 scaled symmetry/null, positive Cholesky와 `rcond_inf>=2^-32`를 검사한다. 물리 DOF를 numeric rank cutoff로 삭제하지 않는다. C05/C06의 finite rotation·warped/oracle·D3·단위변환 test와 실행 evidence는 owner 완료 후 확인한다. | +| HHT-IMP-NR-005 | **P2; corrected, threshold regression passed** | 작은 각에서 기존 직접 삼각함수 D3의 cancellation 위험을 안정 series로 수정했다. `|phi|<.5`의 squared-norm 9차 coefficient recurrence와 최종 회전행렬 extras 복사를 대조했으며 작은각 경계 양쪽을 포함한 4개 C++ sweep의 통과를 확인했다. | +| HHT-IMP-NR-006 | **P2; closed, exact scale restored** | 기존 `Σe |D²f_e[v,v]|`는 element 내부 coordinate-pair 상쇄를 잃어 N3와 달랐다. 수정본은 아래 exact scalar ψ Hessian을 global `(q_j,q_k)` identity로 합친 **뒤** `Σjk |ψ_qq,jk v_j v_k|`를 누적하고 algebra가 이를 소비한다. N3/[I/O 계약](io-contract.md) 의미와 일치한다. Finite/warped ψ·gradient·Hessian 독립 test 통과도 확인했다. | +| HHT-IMP-NR-007 | **P2; source corrected, execution pending** | 최초 `BuildPhysicalSpace`에 없던 N4 `PᵀZ` 및 `[P Z]` Gram guard가 추가됐다. Nodal sparse block의 실제 free dimension과 global Frobenius 오차, normalized null Gram을 검사하며 `1e-12` 기준을 유지한다. 수정본 실행 evidence는 대기다. | +| HHT-IMP-NR-008 | **수치 절차 변경; 문서 정합성 필요** | Reference의 큰 DIRECT step에서 constant-A 초기 추정이 geometry를 벗어나, Newton 시작 좌표를 `q_n`으로 변경했다. HHT predictor/종속 rates/endpoint 식은 그대로이며 cutback/retry도 없어 수식 변경은 아니다. 다만 N2가 기존 시작값을 명시하므로 새 deterministic 시작값과 실패 동작을 N2에 명시해야 한다. 기존 문구 상태로 계약 일치를 주장하지 않는다. | + +NR-001의 독립 scalar 반례: `phi=(.5,0,0)`, `d0=eZ`, Jl의 yz block은 `[[sin(.5)/.5,−(1−cos(.5))/.5],[(1−cos(.5))/.5,sin(.5)/.5]]`다. 관찰한 inverse 식을 적용하면 unit null의 yz 성분이 `(-.6816387600,.7316888689)`이고 `||B n||=||(Jl n)×d||=.4744471057`이다. Python 표준 산술로 실제 계산했다. 이는 zero rotation에서 숨는 부호 오류이며 같은 잘못된 n을 미분한 finite-difference test만으로는 검출되지 않는다. + +## 일치한 수식과 남은 대조 + +최초 snapshot의 `H`, `B_x`, `C_x` 조립과 `I=∫rho Hᵀ(Ha+B_x[v,v])`, `DqI`, `DvI=2∫rho HᵀB_x[v,p]`, `M`의 표시된 index contraction은 정식화와 일치했다. `P/L` 및 `Pdot/Ldot`는 같은 현재 x와 물리적 acceleration을 사용하고 `Ldot=∫rho x×xddot`의 `v×v=0` 항등식과 맞는다. Referenced density/measure, 3³ mass quadrature, current geometry 검사와 fixed-global rate 의미도 대조했다. 이것은 runtime 통과 판정이 아니다. + +Internal Jet의 product 6항 및 Compose의 `f''' x_i x_v²+f''(2x_iv x_v+x_i x_vv)+f' x_ivv`를 독립 대조했다. Atan2의 bounded-ratio 미분은 quadrant 상수의 미분이 0인 점과 맞으며 기존 first/second reduction을 유지한다. 반환 `force_directional_second=D²f_int[v,v]`, `drilling_metric=Σkd gradγ⊗gradγ` 및 physical/drilling force·energy 분리는 정식화와 일치한다. 실제 procedure의 hidden-constraint 연결은 이 kernel 범위 밖이며 별도 검증이 필요하다. + +NR-005의 독립 scalar 검산은 `phi=(θ,0,0), d0=eZ`에 대해 `D³d_x[eX,eX,eZ]=cosθ/θ−2sinθ/θ²+2(1−cosθ)/θ³`와 안정 전개 `−θ/4+θ³/36−θ⁵/960+θ⁷/50400`를 Python binary64로 비교했다. `θ=1e-4`에서 각각 `−7.7441382018e-5`, `−2.4999999972e-5`, 대칭 상대오차 `0.512`이고 `θ=1e-4*(1+2^-10)`에서 `0.339`다. 이는 algebraic scalar cancellation의 재현이며 수정 전 C++ Jet 출력 자체를 실행했다고 주장하지 않는다. 관성 rotation kernel의 32항 matrix exponential series에는 같은 직접 뺄셈이 없다. `|phi|