From adb809bc01c862f252a96121e2ee39059257ced6 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Mon, 19 Oct 2026 02:09:22 +0900 Subject: [PATCH] =?UTF-8?q?=EC=97=84=EA=B2=A9=20=EC=88=98=EB=A0=B4=20?= =?UTF-8?q?=EA=B8=B0=EC=A4=80=20=EC=8B=A4=ED=96=89=20=EB=B3=B4=EC=A1=B4?= =?UTF-8?q?=EA=B3=BC=20=ED=95=9C=EC=A0=95=20HHT=20=EA=B0=9C=EB=B0=9C=20?= =?UTF-8?q?=EA=B3=84=EC=B8=A1=20=EA=B8=B0=EB=A1=9D?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 1 + docs/abaqus-tolerance/results.md | 158 ++++++++++++++++++ msvc/Fesa.Common.props | 4 + msvc/fesa_solver.vcxproj | 1 + msvc/fesa_unit_tests.vcxproj | 1 + src/fesa/CMakeLists.txt | 4 + src/fesa/analysis/hht_development_probe.h | 141 ++++++++++++++++ src/fesa/analysis/hht_dynamic_assembly.cpp | 5 + src/fesa/analysis/hht_reduced_operators.cpp | 5 + src/fesa/analysis/hht_reduced_procedure.cpp | 12 ++ tests/CMakeLists.txt | 4 + .../analysis/hht_development_probe_test.cpp | 114 +++++++++++++ 12 files changed, 450 insertions(+) create mode 100644 docs/abaqus-tolerance/results.md create mode 100644 src/fesa/analysis/hht_development_probe.h create mode 100644 tests/unit/analysis/hht_development_probe_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 5cd78db..1dbc94b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,6 +3,7 @@ project(FESA VERSION 0.1.0 LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) +option(FESA_HHT_DEVELOPMENT_PROBE "Build private HHT development counters" OFF) set(CMAKE_MSVC_RUNTIME_LIBRARY "MultiThreaded$<$:Debug>DLL") # Pin the original backend ABI/link/thread choices rather than package defaults. set(MKL_INTERFACE ilp64) diff --git a/docs/abaqus-tolerance/results.md b/docs/abaqus-tolerance/results.md new file mode 100644 index 0000000..bfe029a --- /dev/null +++ b/docs/abaqus-tolerance/results.md @@ -0,0 +1,158 @@ +# 외부 수렴 정책 실험 결과 — T2 기준 보존 + +기록일: 2026-10-19. **정상 strict 기준은 확보했으나 전체 실험·물리 검증은 미완료**다. +후보 정책은 아직 구현·측정하지 않았다. 가속률이나 채택 가능성을 주장하지 않는다. + +## 기준과 원시 증거 + +측정 source는 `78cb33920742a91b4801788fbcc953d3306b0884`다. +`d7f6d0c`와 제품 source·두 빌드 설정의 diff는 exit 0이다. 원래 조사 기준 +`43867d0`와 달리 사용자 승인 local kernel인 polar 1회, ExpSkew/Jacobian 3차, +질량 gate와 관련 테스트 오차 `1e-3`을 사용했다. 외부 Newton은 잔차와 +q/v/a/A 보정 `1e-10`, 최대 100회인 기존 strict 경로다. + +원시 자료는 **부모 작업공간**의 다음 ignored 위치에 보존했다. T2 worktree의 +`.build`를 정리해도 아래 자료는 별도로 남는다. + +```text +.build/tolerance-experiment/go-20261018/baseline/ + source/strict-source-78cb339.zip, source-hashes.json, revision.txt + runtime/ # 공식 EXE·DLL·PDB + runtime/tests/ # 같은 Release 검사 EXE·DLL·PDB + build/ # 실제 명령·전체 로그·각 exit·선택 CMake compile 명령 + run/ # 새 INP·CSV 복사본·로그·새 HDF5 + execution.json, runtime-hashes.json, test-runtime-hashes.json + input-hashes-before.json, input-hashes-after.json, environment-before.json + reference-comparison.json, official-log-summary.json + development-one-increment/ # 최초 ON 진단, 보완 전 + development-off-control/ # 보존된 공식 OFF EXE의 1증분 대조 + development-final/ # 보완 후 최종 ON 진단·카운터 관계 검증 +``` + +| 보존 대상 | SHA256 | +| --- | --- | +| strict source ZIP | `0B655E0257720AD9CEF31274931578FE17C4456B9988C158A61DF88D93792BED` | +| 공식 fesa.exe | `3CDB5D9D35F6AD95781044BD311ECB35C9FB121780EF81FC535434BF3515172C` | +| 정식 INP 및 공식 작업 복사 | `B8922243CA8DDFE20B1DD39F3123207F789A1A75F6B10DBB6B0BC0B5A40DFFF5` | +| 정식 변위 CSV 및 복사 | `422ECF17609E14FD119D60A3D5AF7A0D1BC4A95E0F556608A7D791E2718C56C5` | +| 새 공식 HDF5 | `55EF1FFAAE061311F3C1AE705824658B668A818D5F1F7D27BB538CA66EE4792C` | + +원본 reference는 실행·수정하지 않았고 전후 hash가 같다. 기존 결과를 재사용하지 않았다. + +## 계측 없는 공식 실행 + +`MSBuild FESA.sln /m:2 /p:Configuration=Release /p:Platform=x64`는 native exit 0이다. +VS 18.11/MSVC 14.51 x64, `/O2 /Ob2 /fp:precise /MD`, staged MKL 2026.1, +TBB 2023.1, HDF5 2.1.1, Intel runtime 2026.1을 사용했다. +기존 Hex8 검사 source의 C4334 경고 1개가 있었고 오류는 없었다. + +호스트는 KOKO, Ryzen 7 5800H(8코어/16논리), 전원은 고성능이다. +관련 thread 환경변수와 PATH를 `environment-before.json`에 기록했고 변경하지 않았다. +실제 worker 수를 고정·계측한 결과는 아니다. 후보 공식 측정도 같은 환경을 사용해야 한다. +측정 중 다른 build/test/solver를 실행하지 않았다. 워밍업 없이 정확히 1회 실행했다. + +| 항목 | 실제 결과 | +| --- | --- | +| 명령 | 보존된 `runtime/fesa.exe run/shellNL-dyn-meshfine.inp` | +| 입력 | dt=.1, T=10, 100증분, 729 nodes, 676 elements | +| 시작/종료 UTC | 2026-10-18 15:01:20.6895159 / 16:47:11.4199462 | +| native / wrapper exit | 0 / 0 | +| 공식 elapsed | **6350.7895521초** | +| 프로세스 CPU | 18021.3125초; elapsed와 다른 값 | +| 새 HDF5 | 126254676 bytes, schema2, 729 nodes / 101 frames | +| peak working set | 종료 후 API가 null을 반환하여 미확보 | +| 정식 비교 | comparator exit 0, 사전검사 통과, U blocking 실패 0, UR 경고 200 | +| U 최대 family relative RMS | 0.002282424502275574 | + +정식 비교는 기존 9-node U blocking/UR warning 정책 그대로 수행했다. +전체 441774 변위 scalar의 identity/finite 확인 후 202 families를 평가했고 +near-zero 1%, 일반 상대 5%, RMS 1%, zero-scale exact-zero를 유지했다. +Python의 로컬 uv cache를 child `PYTHONPATH`로 읽어 numpy 2.4.6/h5py 3.16.0을 사용했다. +설치·다운로드·캐시 쓰기를 하지 않았다. 실제 argv와 exit는 원시 로그에 있다. + +공식 hht-log-v1에서 양의 accepted 100건, 평가 204건이며 accepted iteration 합은 +**outer applied correction 204회**다. 96증분은 2회, 4증분은 3회다. +잔차가 통과했지만 보정이 실패한 평가는 27건이다. 보정 채널별 실패 평가 건수는 +u의 q/v/a/A가 102/103/104/104, 회전은 100/100/102/102다. +이는 중복을 허용한 로그 채널 진단이며 배타적 실행시간이나 후보의 절약 가능 회수가 아니다. + +## 별도 1증분 개발 계측 + +기존 4개 물리·미분 검사 실패가 남은 상태에서 동일 100증분을 장시간 반복하지 않았다. +승인된 별도 개발 입력 복사에서 T만 `.1`로 줄였고 dt=.1과 local kernel·strict 종료는 유지했다. +이 입력 SHA256은 `79D1C09FC5AFB792312F1908E13D71ABC228D98D10B438B91412A9BFA4C95A50`다. + +계측은 기본 컴파일 OFF다. Native는 `/p:FesaHhtDevelopmentProbe=true`, CMake는 +`-DFESA_HHT_DEVELOPMENT_PROBE=ON`으로 켜고, ON 실행에서만 +`FESA_HHT_PROBE_OUTPUT=<새 sidecar 경로>`를 설정한다. 정상 OFF 객체에는 해당 환경 키가 없다. +추가 요소 평가·수치 보정·판정 완화 없이 기존 호출 앞의 계수와 serial scope의 정수 ns 타이머만 쓴다. +단일 procedure 개발 실행용이며 동시 다중 procedure 계측을 지원한다는 주장은 없다. +병렬 요소 callback은 atomic 계수만 쓰고 타이머/단계 배열은 serial thread만 수정한다. + +최종 개발 실행 native exit 0, 개발 elapsed 77.5268139초다. 이는 공식 시간이 아니다. +최초 ON, OFF 대조, 보완 후 최종 ON은 모두 별도 파일로 보존했다. + +| 최종 sidecar 계수 | 값과 범위 | +| --- | --- | +| outer_applied_correction / inner_applied_correction | 3 / 1; 실제 좌표 갱신만 포함 | +| assembly_calls | 14 | +| assembly_internal_evaluation | 9464 = 14×676 | +| assembly_internal_failure | 3; 대체 시작 trial 평가의 실패이며 최종 해석 실패가 아님 | +| assembly_inertia_evaluation | 9461 = 9464−3; internal 실패 뒤 생략된 관성 호출 제외 | +| constraint_internal_evaluation | 6760; cache miss의 실제 추가 EvaluateTrial만 포함 | +| logical_factor_call / logical_solve_call / logical_rhs | 1274 / 7516 / 7516 | + +Solver 수는 HHT 경계 forwarding wrapper의 논리 호출이다. 실패 시도도 세며 vendor의 +phase11/22/33 전체 호출수로 해석하지 않는다. inner correction도 stagnation 때문에 +계산만 하고 적용하지 않은 solve는 제외하지만 logical solve에는 포함한다. +계수 관계와 공식 첫 증분의 outer 3회 일치는 `counter-check.json`에서 assert했다. + +| 중복 없는 최종 단계 시간 | 초 | +| --- | ---: | +| procedure_other | 30.6919325 | +| inner_graph | 26.6134860 | +| assembly | 13.8891082 | +| outer_linear | 0.1445448 | +| physical_rates | 0.2650750 | +| factor | 5.0312963 | +| solve | 0.5715257 | +| 합계 = procedure inclusive | 77.2069685 | + +각 `_exclusive_ns`는 자식 시간을 뺀 값이고 합계가 `procedure_inclusive_ns`와 정확히 같음을 +assert했다. ExecuteReducedHht 경계이므로 CLI 입력·최종 HDF5 finalize·sidecar I/O를 제외한다. +요소 kernel 시간은 병렬 wall-time을 포함한 assembly 범위이며 각 요소 CPU 시간의 합이 아니다. +전체 100증분의 inner/요소/factor/solve 비용과 후반 채널 분포는 **계측하지 않았다**. +첫 증분 계수·시간을 100배하거나 후보 전체 비용으로 외삽하지 않는다. + +## 정확 byte parity와 검증 한계 + +공식 frame0/1과 개발 결과의 40개 numeric dataset을 identity/shape/type/finite 확인 후 비교했다. +ON 최초·ON 최종·보존된 공식 OFF EXE의 짧은 입력 대조 모두 22개 dataset이 byte 일치하고 +frame1의 18개 dataset은 불일치했다. **정확 byte parity는 미통과**다. +최종 ON의 q/v/physical a 최대 절대차는 각각 6.776263578034403e-21, +2.168404344971009e-19, 6.938893903907228e-18이고 energy는 8.881784197001252e-16이다. +거의 영인 값의 scalar 상대오차가 클 수 있으므로 최대 절대차만으로 통과 판정하지 않는다. + +동일 보존 OFF EXE의 대조에서도 q/v/a 최대 절대차는 1.3552527156068805e-20, +3.552713678800501e-15, 5.684341886080802e-14였다. T 단축·반복 실행·기본 thread 조건· +컴파일/계측 영향을 완전히 분리하지 못했으며 원인을 하나로 단정하지 않는다. +새 parity tolerance나 floor를 도입하지 않았다. dataset별 차이 수·최대 절대/상대차·첫 위치와 +원본/실제 bytes는 각 `parity-details.json`에 보존했다. NaN을 동일 처리하지 않고 finite를 요구했다. +독립 source 검토에서 산술·평가·진단 순서 변경이 없다고 확인한 사실은 실행 parity 통과와 별개다. + +계측 검사는 no-op에서 2FAIL/native1을 확인한 뒤 카운트·중첩 배타 시간·solver forwarding/ +실패 보존·실제 inner 1회/0회 갱신의 4개 검사를 native0으로 통과했다. +실제 inner 계수 검사는 CMake OFF operators 객체를 연결하면 1FAIL/native1, ON이면 4PASS/native0다. +Standalone 준비 중 link/include/MKL runtime 오류와 짧은 clock tick assertion 실패도 별도 보존했다. +짧은 tick은 0ns를 허용하고 scope 진입 수와 시간 분할의 exact equality를 검사한다. +Standalone에 필요한 기존 설치의 sequential DLL은 테스트 전용이며 공식 runtime에 추가하지 않았다. + +Native ON CLI 최종 빌드는 exit0(경고0), CMake Release 기본 OFF의 변경 production 3개와 +새 검사 객체 compile은 exit0, ON 새 검사 객체 compile도 exit0이다. source 등록을 양쪽에 맞췄다. +OFF에서는 새 검사 identity가 없고 ON에서는 4개다. 전체 CMake link/CTest나 전체 Debug 검사는 +이번 T2에서 실행하지 않았고 후속 전체 회귀를 대신하지 않는다. +기존 local-kernel의 4FAIL은 [기록](local-kernel-experiment.md)을 유지하며 재실행/해결한 것으로 보고하지 않는다. +U 통과는 V/a/energy의 외부 검증이나 기존 물리 실패의 해결을 뜻하지 않는다. + +T2 판정은 **기준 보존 완료·개발 진단 한정 확보·전체 검증 미완료**다. +후보 측정, 전체 100증분 개발 비용, 정확 parity 원인 분리와 필수 전체 물리/회귀 검증은 남았다. diff --git a/msvc/Fesa.Common.props b/msvc/Fesa.Common.props index d8e4e27..34bf8c1 100644 --- a/msvc/Fesa.Common.props +++ b/msvc/Fesa.Common.props @@ -3,6 +3,7 @@ + false $([System.IO.Path]::GetFullPath('$(MSBuildThisFileDirectory)..\')) $(FesaRoot).build\native\$(Configuration)\ $(FesaRoot.Replace('\','/').TrimEnd('/')) @@ -20,6 +21,9 @@ $(FesaExternalOptions) /external:I"$(FESA_MKL_ROOT)\include" + + FESA_HHT_DEVELOPMENT_PROBE=1;%(PreprocessorDefinitions) + stdcpp17 false diff --git a/msvc/fesa_solver.vcxproj b/msvc/fesa_solver.vcxproj index 897f7f5..3ab073a 100644 --- a/msvc/fesa_solver.vcxproj +++ b/msvc/fesa_solver.vcxproj @@ -14,6 +14,7 @@ + diff --git a/msvc/fesa_unit_tests.vcxproj b/msvc/fesa_unit_tests.vcxproj index d0826de..80e95ee 100644 --- a/msvc/fesa_unit_tests.vcxproj +++ b/msvc/fesa_unit_tests.vcxproj @@ -20,6 +20,7 @@ + diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index 475748b..804c94a 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -7,6 +7,7 @@ add_library( analysis/eigen_reduction.cpp analysis/hht_time_integration.cpp analysis/hht_dynamic_algebra.cpp + analysis/hht_development_probe.h analysis/hht_dynamic_assembly.cpp analysis/hht_reduced_operators.cpp analysis/hht_reduced_procedure.cpp @@ -126,6 +127,9 @@ target_compile_options( PRIVATE $<$:/W4 /WX /EHsc> ) +if(FESA_HHT_DEVELOPMENT_PROBE) + target_compile_definitions(fesa_solver PRIVATE FESA_HHT_DEVELOPMENT_PROBE=1) +endif() # Keep the CPU/OS dispatch on the baseline ISA; only its private kernel uses AVX2. if(CMAKE_CXX_COMPILER_ID STREQUAL "MSVC" AND diff --git a/src/fesa/analysis/hht_development_probe.h b/src/fesa/analysis/hht_development_probe.h new file mode 100644 index 0000000..803203c --- /dev/null +++ b/src/fesa/analysis/hht_development_probe.h @@ -0,0 +1,141 @@ +#ifndef FESA_ANALYSIS_HHT_DEVELOPMENT_PROBE_H_ +#define FESA_ANALYSIS_HHT_DEVELOPMENT_PROBE_H_ + +// This header is compiled out of normal builds, including environment/file I/O. +#if defined(FESA_HHT_DEVELOPMENT_PROBE) +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/solvers/linear/linear_solver.h" +#include "fesa/solvers/linear/tangent_solver.h" + +namespace fesa::hht_development_probe { +enum Stage { kProcedure, kInner, kAssembly, kOuterSolve, kRates, + kFactor, kSolve, kStageCount }; +enum Counter { kOuterCorrection, kInnerCorrection, kInternalEvaluation, + kInertiaEvaluation, kFactorCall, kSolveCall, kRhs, + kConstraintInternalEvaluation, kInternalFailure, kCounterCount }; +struct Statistics { + std::array, kCounterCount> counts{}; + // Timers/stage attribution are serial procedure-thread data. Only the three + // element counters below and internal failures may be touched by workers. + std::array inclusive_ns{}; + std::array exclusive_ns{}; + std::array stage_calls{}; + std::array, kStageCount> by_stage{}; +}; +inline Statistics* active = nullptr; +class Scope { + public: + explicit Scope(Stage stage) : stage_(stage), parent_(current_), + start_(std::chrono::steady_clock::now()) { + current_ = this; + if (active != nullptr) ++active->stage_calls[stage_]; + } + ~Scope() { + const auto elapsed = std::chrono::duration_cast( + std::chrono::steady_clock::now() - start_).count(); + if (active != nullptr) { + active->inclusive_ns[stage_] += elapsed; + active->exclusive_ns[stage_] += elapsed - children_; + } + if (parent_ != nullptr) parent_->children_ += elapsed; + current_ = parent_; + } + static Stage Current() { return current_ == nullptr ? kProcedure : current_->stage_; } + private: + Stage stage_; + Scope* parent_; + std::chrono::steady_clock::time_point start_; + std::int64_t children_ = 0; + inline static thread_local Scope* current_ = nullptr; +}; +inline void Count(Counter counter, std::uint64_t count = 1U) { + if (active == nullptr) return; + active->counts[counter].fetch_add(count, std::memory_order_relaxed); + // Element callbacks are parallel; their totals use only the atomics above. + if (counter != kInternalEvaluation && counter != kInertiaEvaluation && + counter != kConstraintInternalEvaluation && counter != kInternalFailure) + active->by_stage[Scope::Current()][counter] += count; +} +class Session { + public: + Session() : previous_(active) { active = &statistics_; } + ~Session() { + active = previous_; +#if defined(_WIN32) + char* environment = nullptr; + std::size_t size = 0U; + if (_dupenv_s(&environment, &size, "FESA_HHT_PROBE_OUTPUT") != 0) return; + const std::string path = environment == nullptr ? "" : environment; + std::free(environment); +#else + const char* environment = std::getenv("FESA_HHT_PROBE_OUTPUT"); + const std::string path = environment == nullptr ? "" : environment; +#endif + if (path.empty()) return; + std::ofstream output(path, std::ios::out | std::ios::binary); + constexpr const char* counters[] = {"outer_applied_correction", "inner_applied_correction", + "assembly_internal_evaluation", "assembly_inertia_evaluation", + "logical_factor_call", "logical_solve_call", "logical_rhs", + "constraint_internal_evaluation", "assembly_internal_failure"}; + constexpr const char* stages[] = {"procedure_other", "inner_graph", "assembly", + "outer_linear", "physical_rates", "factor", "solve"}; + for (int c = 0; c < kCounterCount; ++c) + output << counters[c] << '=' << statistics_.counts[c].load() << '\n'; + output << "procedure_inclusive_ns=" << statistics_.inclusive_ns[kProcedure] << '\n'; + for (int s = 0; s < kStageCount; ++s) { + output << stages[s] << "_calls=" << statistics_.stage_calls[s] << '\n'; + output << stages[s] << "_exclusive_ns=" << statistics_.exclusive_ns[s] << '\n'; + for (int c = 0; c < kCounterCount; ++c) + if (statistics_.by_stage[s][c] != 0U) + output << stages[s] << '_' << counters[c] << '=' + << statistics_.by_stage[s][c] << '\n'; + } + } + private: + Statistics statistics_; + Statistics* previous_; +}; +class PositiveSolver final : public LinearSolver { + public: + explicit PositiveSolver(LinearSolver& solver) : solver_(solver) {} + Status Factorize(const SparseMatrix& matrix) override { + Count(kFactorCall); Scope scope{kFactor}; return solver_.Factorize(matrix); + } + Status Solve(const Vector& rhs, Vector& solution) const override { + Count(kSolveCall); Count(kRhs); Scope scope{kSolve}; + return solver_.Solve(rhs, solution); + } + private: + LinearSolver& solver_; +}; +class CorrectionSolver final : public TangentSolver { + public: + explicit CorrectionSolver(TangentSolver& solver) : solver_(solver) {} + Status FactorizeFresh(const SparseMatrix& matrix) override { + Count(kFactorCall); Scope scope{kFactor}; return solver_.FactorizeFresh(matrix); + } + Status SolveCorrection(const Vector& rhs, Vector& solution) const override { + Count(kSolveCall); Count(kRhs); Scope scope{kSolve}; + return solver_.SolveCorrection(rhs, solution); + } + private: + TangentSolver& solver_; +}; +} // namespace fesa::hht_development_probe +#define FESA_HHT_PROBE_COUNT(counter) \ + ::fesa::hht_development_probe::Count(::fesa::hht_development_probe::counter) +#define FESA_HHT_PROBE_SCOPE(stage) \ + ::fesa::hht_development_probe::Scope probe_scope{ \ + ::fesa::hht_development_probe::stage} +#else +#define FESA_HHT_PROBE_COUNT(counter) ((void)0) +#define FESA_HHT_PROBE_SCOPE(stage) ((void)0) +#endif +#endif // FESA_ANALYSIS_HHT_DEVELOPMENT_PROBE_H_ diff --git a/src/fesa/analysis/hht_dynamic_assembly.cpp b/src/fesa/analysis/hht_dynamic_assembly.cpp index 17012ab..bec9a3f 100644 --- a/src/fesa/analysis/hht_dynamic_assembly.cpp +++ b/src/fesa/analysis/hht_dynamic_assembly.cpp @@ -12,6 +12,7 @@ #endif #include "analysis/hht_dynamic_internal.h" +#include "analysis/hht_development_probe.h" #include "math/dense_blas_internal.h" #include "fesa/assembly/parallel_for.h" #include "fesa/elements/finite_rotation.h" @@ -382,6 +383,7 @@ Result EvaluateAssembly( const Vector& load_scale, const double time_increment, const HhtParameters& parameters, const ParallelFor& parallel_for, const bool directional_third) { + FESA_HHT_PROBE_SCOPE(kAssembly); const std::size_t count = dofs.FullDofCount(); for (const auto* values : std::array{ &frame.generalized_coordinates, &frame.velocity, &frame.acceleration, @@ -446,9 +448,11 @@ Result EvaluateAssembly( } const auto before = retain_constraints ? ReadTrialEvaluationContext() : TrialEvaluationContext{}; + FESA_HHT_PROBE_COUNT(kInternalEvaluation); auto f = shells[e].internal->EvaluateTrial( trial, directional_third ? &v : nullptr); if (!f.HasValue()) { + FESA_HHT_PROBE_COUNT(kInternalFailure); failures[e] = f.GetStatus(); return; } @@ -456,6 +460,7 @@ Result EvaluateAssembly( const auto after = ReadTrialEvaluationContext(); if (ReusableTrialContext(before, after)) evaluation_contexts[e] = after; } + FESA_HHT_PROBE_COUNT(kInertiaEvaluation); auto m = shells[e].inertia->Evaluate(trial, v, a); if (!m.HasValue()) { failures[e] = m.GetStatus(); diff --git a/src/fesa/analysis/hht_reduced_operators.cpp b/src/fesa/analysis/hht_reduced_operators.cpp index 4f0b94b..be4781f 100644 --- a/src/fesa/analysis/hht_reduced_operators.cpp +++ b/src/fesa/analysis/hht_reduced_operators.cpp @@ -1,4 +1,5 @@ #include "analysis/hht_reduced_internal.h" +#include "analysis/hht_development_probe.h" #include #include @@ -513,6 +514,7 @@ Result EvaluateConstraintOperators( shell_index++, shells.size(), q, v); std::optional fresh; if (saved == nullptr) { + FESA_HHT_PROBE_COUNT(kConstraintInternalEvaluation); auto internal = shell.internal->EvaluateTrial(Trial(q), &v); if (!internal.HasValue()) @@ -713,6 +715,7 @@ Result SolveEmbeddingRoot(const FixedChart& chart, const Vector& p, const Vector& seed, const ConstraintEvaluator& evaluate, TangentSolver& solver) { + FESA_HHT_PROBE_SCOPE(kInner); if (!evaluate || p.Size() != chart.p_dofs.size() || seed.Size() != chart.z_dofs.size()) return Result::Failure( @@ -830,6 +833,7 @@ Result SolveEmbeddingRoot(const FixedChart& chart, prior_q = q.Value(); has_prior_z = true; z = std::move(next_z); + FESA_HHT_PROBE_COUNT(kInnerCorrection); } return Result::Failure(Failure("Unreachable graph loop.")); } @@ -876,6 +880,7 @@ Result EvaluateEmbeddingFirstDerivatives( Result ReconstructPhysicalRates( const FixedChart& chart, const EmbeddingRoot& root, const Vector& vp, const Vector& ap, TangentSolver& solver) { + FESA_HHT_PROBE_SCOPE(kRates); if (!ValidBasis(chart) || root.q.Size() != chart.full_dof_count || root.operators.rows.size() != chart.z_dofs.size() || vp.Size() != chart.p_dofs.size() || ap.Size() != chart.p_dofs.size() || diff --git a/src/fesa/analysis/hht_reduced_procedure.cpp b/src/fesa/analysis/hht_reduced_procedure.cpp index efcd372..408b8e4 100644 --- a/src/fesa/analysis/hht_reduced_procedure.cpp +++ b/src/fesa/analysis/hht_reduced_procedure.cpp @@ -1,4 +1,5 @@ #include "analysis/hht_reduced_procedure.h" +#include "analysis/hht_development_probe.h" #include "analysis/hht_reduced_integrator.h" #include "equilibrium_algorithm.h" @@ -396,8 +397,16 @@ Result ExecuteReducedHht( TangentSolver& tangent_solver) { const ReducedCoordinateShellHhtIntegrator integrator; NewtonStepAlgorithm algorithm; +#if defined(FESA_HHT_DEVELOPMENT_PROBE) + hht_development_probe::Session probe_session; + hht_development_probe::PositiveSolver measured_positive{positive_solver}; + hht_development_probe::CorrectionSolver measured_tangent{tangent_solver}; + return ExecuteReducedHht(domain, log, parallel_for, measured_positive, + measured_tangent, integrator, algorithm); +#else return ExecuteReducedHht(domain, log, parallel_for, positive_solver, tangent_solver, integrator, algorithm); +#endif } Result ExecuteReducedHht( @@ -405,6 +414,7 @@ Result ExecuteReducedHht( const ParallelFor& parallel_for, LinearSolver& positive_solver, TangentSolver& tangent_solver, const ShellHhtIntegrator& integrator, EquilibriumAlgorithm& algorithm) { + FESA_HHT_PROBE_SCOPE(kProcedure); if (domain.Dynamic() == nullptr) return Result::Failure(Status::Failure( FailureCategory::kModel, {{Severity::kError, "invalid-dynamic-step", {}, @@ -1016,6 +1026,7 @@ Result ExecuteReducedHht( }; Vector selected_solution{selected.Value().rhs.Size()}; if (selected_solution.Size() > 0U) { + FESA_HHT_PROBE_SCOPE(kOuterSolve); status = tangent_solver.FactorizeFresh(selected.Value().tangent); if (!status.IsOk()) return linear_failure(status); @@ -1091,6 +1102,7 @@ Result ExecuteReducedHht( seed = trial.Value().root.z; for (std::size_t i = 0U; i < seed.Size(); ++i) seed[i] += correction.Value().z[i]; + FESA_HHT_PROBE_COUNT(kOuterCorrection); return StepResult::Success(std::nullopt); }, kMaximumOuterCorrections); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5845708..1e77eac 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -14,6 +14,7 @@ add_executable( unit/analysis/hht_energy_verification_test.cpp unit/analysis/hht_temporal_verification_test.cpp unit/analysis/hht_dynamic_algebra_test.cpp + unit/analysis/hht_development_probe_test.cpp unit/analysis/hht_reduced_shell_test.cpp unit/analysis/hht_dynamic_procedure_test.cpp unit/analysis/nonlinear_analysis_state_test.cpp @@ -99,6 +100,9 @@ target_include_directories( PRIVATE "${PROJECT_SOURCE_DIR}/src/fesa" ) +if(FESA_HHT_DEVELOPMENT_PROBE) + target_compile_definitions(fesa_unit_tests PRIVATE FESA_HHT_DEVELOPMENT_PROBE=1) +endif() target_link_libraries( fesa_unit_tests diff --git a/tests/unit/analysis/hht_development_probe_test.cpp b/tests/unit/analysis/hht_development_probe_test.cpp new file mode 100644 index 0000000..3fa5cea --- /dev/null +++ b/tests/unit/analysis/hht_development_probe_test.cpp @@ -0,0 +1,114 @@ +#include "analysis/hht_development_probe.h" + +#include + +#if defined(FESA_HHT_DEVELOPMENT_PROBE) +#include "fesa/math/sparse_matrix.h" +#include "fesa/math/vector.h" +#include "analysis/hht_reduced_internal.h" +namespace fesa::hht_development_probe { +class IdentityTangentSolver final : public TangentSolver { + public: + Status FactorizeFresh(const SparseMatrix&) override { return Status::Ok(); } + Status SolveCorrection(const Vector& rhs, Vector& result) const override { + result = rhs; return Status::Ok(); + } +}; +TEST(HhtDevelopmentProbe, CountsActualInnerUpdatesForLinearGraphAndZeroRoot) { + hht_reduced_internal::FixedChart chart{}; + chart.full_dof_count = 1U; + chart.z_dofs = {0U}; + chart.z_columns = {{0U, 1U, false, {{0U, 1.0}}, {0.0}}}; + IdentityTangentSolver solver; + const auto evaluate = [](const Vector& q) { + hht_reduced_internal::ConstraintOperators operators; + operators.rows.resize(1U); + operators.rows[0U].value = q[0U]; + operators.rows[0U].gradient = {{0U, 1.0}}; + return Result::Success( + std::move(operators)); + }; + Statistics statistics; + active = &statistics; + const auto first = hht_reduced_internal::SolveEmbeddingRoot( + chart, Vector{0U}, Vector{1U, 1.0}, evaluate, solver); + active = nullptr; + ASSERT_TRUE(first.HasValue()); + EXPECT_EQ(statistics.counts[kInnerCorrection].load(), 1U); + EXPECT_DOUBLE_EQ(first.Value().z[0U], 0.0); + Statistics zero_statistics; + active = &zero_statistics; + const auto zero = hht_reduced_internal::SolveEmbeddingRoot( + chart, Vector{0U}, Vector{1U, 0.0}, evaluate, solver); + active = nullptr; + ASSERT_TRUE(zero.HasValue()); + EXPECT_EQ(zero_statistics.counts[kInnerCorrection].load(), 0U); +} +class FailingSolver final : public LinearSolver, public TangentSolver { + public: + Status Factorize(const SparseMatrix&) override { + ++factor_calls; return Status::Failure(FailureCategory::kSolver, {}); + } + Status FactorizeFresh(const SparseMatrix& matrix) override { return Factorize(matrix); } + Status Solve(const Vector&, Vector&) const override { + ++solve_calls; return Status::Failure(FailureCategory::kSolver, {}); + } + Status SolveCorrection(const Vector& rhs, Vector& result) const override { + return Solve(rhs, result); + } + int factor_calls = 0; + mutable int solve_calls = 0; +}; +TEST(HhtDevelopmentProbe, ForwardsEveryLogicalCallAndPreservesFailure) { + Statistics statistics; + active = &statistics; + FailingSolver backend; + PositiveSolver positive{backend}; + CorrectionSolver tangent{backend}; + auto matrix = SparseMatrix::FromCoo(0U, 0U, {}, SparsePattern{{0U}, {}}); + ASSERT_TRUE(matrix.HasValue()); + const Vector rhs{0U}; + Vector result{0U}; + EXPECT_FALSE(positive.Factorize(matrix.Value()).IsOk()); + EXPECT_FALSE(positive.Solve(rhs, result).IsOk()); + EXPECT_FALSE(tangent.FactorizeFresh(matrix.Value()).IsOk()); + EXPECT_FALSE(tangent.SolveCorrection(rhs, result).IsOk()); + active = nullptr; + EXPECT_EQ(backend.factor_calls, 2); + EXPECT_EQ(backend.solve_calls, 2); + EXPECT_EQ(statistics.counts[kFactorCall].load(), 2U); + EXPECT_EQ(statistics.counts[kSolveCall].load(), 2U); + EXPECT_EQ(statistics.counts[kRhs].load(), 2U); +} +TEST(HhtDevelopmentProbe, CountsOnlyWhileActiveAndPreservesMultiplicity) { + Statistics statistics; + active = &statistics; + Count(kOuterCorrection); + Count(kInternalEvaluation, 7U); + active = nullptr; + Count(kOuterCorrection); + EXPECT_EQ(statistics.counts[kOuterCorrection].load(), 1U); + EXPECT_EQ(statistics.counts[kInternalEvaluation].load(), 7U); + EXPECT_EQ(statistics.counts[kInertiaEvaluation].load(), 0U); +} + +TEST(HhtDevelopmentProbe, NestedTimesHaveNoDoubleCounting) { + Statistics statistics; + active = &statistics; + { + Scope parent{kProcedure}; + { Scope child{kAssembly}; } + } + active = nullptr; + EXPECT_EQ(statistics.stage_calls[kProcedure], 1U); + EXPECT_EQ(statistics.stage_calls[kAssembly], 1U); + EXPECT_GE(statistics.inclusive_ns[kProcedure], 0); + EXPECT_GE(statistics.inclusive_ns[kAssembly], 0); + EXPECT_EQ(statistics.exclusive_ns[kAssembly], + statistics.inclusive_ns[kAssembly]); + EXPECT_EQ(statistics.exclusive_ns[kProcedure] + + statistics.exclusive_ns[kAssembly], + statistics.inclusive_ns[kProcedure]); +} +} // namespace fesa::hht_development_probe +#endif