diff --git a/docs/hex8-structural-analysis/PRD.md b/docs/hex8-structural-analysis/PRD.md index f2b2df9..74de1a4 100644 --- a/docs/hex8-structural-analysis/PRD.md +++ b/docs/hex8-structural-analysis/PRD.md @@ -434,7 +434,7 @@ HHT 식·기준 consistent M·입력·고정 schedule·100회·Si/exact-zero· - [x] Solid Domain/parser/3N DOF/ordered assembly/amplitude 검증 - [x] 선형·비선형 정적 accepted state/반력/잔차 검증 - [x] Generalized consistent-mass frequency·최저 N·original-pencil 검증 -- [ ] 선형·비선형 HHT u/v/physical a·시간수렴·수지 검증 +- [x] 선형·비선형 HHT u/v/physical a·시간수렴·수지 내부 검증 - [ ] 공식 HDF5/CLI/self-check/실패 final 보존 검증 - [ ] 정적2·동적2 원본 필수 family 비교 PASS - [ ] 사용자 원본 최저10 modal Value 제공 및 필수 비교 PASS @@ -446,9 +446,14 @@ Local kernel과 입력·3N mapping은 독립 검토 및 연결 후 관련180개 native CTest exit0이며 독립 수치 검토를 통과했다. 두 정적 원본과442절점 고유치 원본은 직접 API로 실행했지만 원본 결과값 비교를 수행한 것은 아니다. 연결 source `28be39f`의 전체 MSVC Debug build는 exit0이며 이 build 자체를 전체 CTest 통과로 해석하지 않는다. -HHT 내부 사례는 진행 중이다. 별도 구현 `f661864`의 관련76개 검사와 독립 수치 검토는 -통과했지만, 두 동적 원본은 세 번째 증분100회 제한에서 엄격 잔차 기준을 만족하지 못했다. -최신 잔차 비율은 선형 `1.77739e-8`, 비선형 `8.73572e-9`이며 기준은 `1e-10`이다. -동적 검증은 미완료이며 원인과 추가 검증은 [lessons](lessons.md)에 기록한다. +과거 단일 Float64 구현 `f661864`는 관련76개 검사와 독립 수치 검토를 통과했지만, +두 동적 원본의 세 번째 증분에서 strict gate를 만족하지 못했다. 이 실패는 +[lessons](lessons.md)에 보존하며 현재 방식의 결과와 구분한다. +승인된 전체 high+low 구현 `c49cae4`는 소규모32/32 및 별도 성공한 LIN/TL·SDOF/MDOF +시간수렴 사례와 독립 numerical/physics 검토를 확보했다. Root 연결 후 관련49/49와 +추가 직렬/TBB low 일치1/1 검사도 native build/CTest0/0이다. 두 원본의 초기 상태와 +증분1–3은 unchanged schedule로 관측했고, 원래 모든 free 성분의 전체 pair 물리 gate를 +독립 재계산했다. 이 관측은 네 번째 증분 전에 종료했으며 전체101 endpoint 완료가 아니다. +상세 범위와 한계는 [수치 검토](numerical-review.md), [물리 검토](physics-review.md)에 있다. 공식 결과·외부 비교·최신 전체 회귀도 남아 있다. 원본 frequency 결과 부재는 구현 진행을 막지 않지만 기능 검증 완료를 막는다. diff --git a/docs/hex8-structural-analysis/lessons.md b/docs/hex8-structural-analysis/lessons.md index 10b402e..5518a9d 100644 --- a/docs/hex8-structural-analysis/lessons.md +++ b/docs/hex8-structural-analysis/lessons.md @@ -259,3 +259,16 @@ trial에서 계산한다. 독립 reviewer가 이 캐시와 native44/44의 원시 독립 numerical reviewer는 실제 shared arithmetic·kernel·geometry의 한정된 범위를 PASS로 판정했다. 재현 로그는 `.harness/hex8-evidence/precision/`에 보존한다. 이 결과는 동적 전체 상태·시간수렴·원본 endpoint·공식 출력의 완료를 뜻하지 않는다. 현재 checklist는 미완료다. + +## 전체 상태의 제한적 원본 관측과 연결 + +전체 high+low T6 `c49cae4`는 원래 시간 간격·INC를 바꾸지 않고 두 원본의 초기 상태와 +증분1–3을 accept했다. Test-only typed stop으로 다음 증분 전에 종료해 prefix 상태를 +관측했으며 native0·2/2를 전체 원본 완료로 보고하지 않는다. 메인이 전체 pair·full M/Mff, +모든1275 free r/Si·zero branch·HHT update·weighted old-r·실제 solve 보정을 exact Fraction으로 +검산했다. 물리 비율 최대값은 선형3.12680e-20, 비선형1.3043731e-16이다. + +Root 연결의 related49/49와 두 brick의 LIN/TL 전체 low 직렬/TBB 일치1/1은 native0이다. +최초 broad filter가 긴 전체 원본 smoke를 함께 선택한 실행은 의도적으로 중단했고 +CTest exit−1로 보존했다. 완료하지 않은 실행을 수렴 실패로 해석하거나 현재 PASS에 합산하지 +않는다. 공식 출력·101 endpoint·원본 비교·전체 회귀는 후속 검증이다. diff --git a/docs/hex8-structural-analysis/numerical-review.md b/docs/hex8-structural-analysis/numerical-review.md index e874602..9f8a452 100644 --- a/docs/hex8-structural-analysis/numerical-review.md +++ b/docs/hex8-structural-analysis/numerical-review.md @@ -20,6 +20,37 @@ T6 consumer·backend·worker 경계와 rollback은 해당 구현이 고정된 동일 graph·선형 force/energy/recovery·비선형 trial Hessian·기존 F64 Create 지원을 유지하며 native44/44 원시 로그를 확인했다. 이 캐시 변경도 bounded PASS다. +## 2026-10-10 — 전체 동적 상태와 원본 세 증분의 제한적 검토 + +동일 독립 reviewer가 고정된 T6 `c49cae4`의 predictor/rates, 기준 M의 전체 pair 작용, +ordered force/tangent 누적, full residual·scale, 실제 old residual을 포함한 weighted 식, +fresh solve와 실제 상태 변경의 네 보정 채널, backend/worker 이후 실행 전제 검사, +accepted low의 복사·검증·rollback을 검토했다. 원시73개 inventory와 owned/shared source, +binary·patch의 hash를 확인했고 해당 consumer 범위 bounded PASS, actionable finding0이다. +Current32/32 검사와 과거 temporal4/5 실행의 성공 사례를 구분한다. + +이후 변경하지 않은 두 원본의 초기 상태와 증분1–3을 관측했다. 원래 h=0.1,T=10,INC=100, +1326 full/1275 free/51 constrained를 유지하고, 네 번째 증분 local 계산 전 test-only 예외로 +종료했다. 관측 native0·2/2는 이 종료 조건의 성공이며 전체 해석 완료가 아니다. +메인은51개 raw artifact의 hash/크기를 확인한 뒤 두 binary64 성분을 exact dyadic Fraction으로 +읽어 모든 상태의 normalization, 실제 Mff의 full M restriction, 원래 모든 free 성분의 +물리 r/Si와 exact-zero, HHT 세 update 식과 weighted old-r 항등식을 독립 재계산했다. + +Exact CSR M 작용으로 재계산한 물리 비율의 최댓값은 선형3.12680e-20, +비선형1.3043731e-16으로 승인1e-10 이내다. Update·weighted identity는1e-12 이내이며 +초기1275개 free row는 모두 exact0이다. 실제 PARDISO 보정의 선형6회·비선형9회 기록으로 +Newton 시작값인 이전 accepted u부터 최종 u를 재구성하고, fresh/observed 네 보정 채널도 +각각 확인했다. 비선형 최대 보정 비율1.01046e-13은 승인1e-10 이내다. + +Storage graph 대조의 별도 절대 오차 budget과 물리 gate를 혼동하지 않는다. Exact M 작용과 +저장된 pair Ma의 상대 차이는 최대2.88e-32였다. Strict physical gate에 floor를 추가하거나 +high만으로 full 상태를 대신하지 않았다. 재현 자료는 로컬 +`.harness/hex8-evidence/precision/{prefix3-frozen,prefix-independent.json,check_prefix.py}`다. +Root 연결 후 관련49/49 검사의 native build/CTest는0/0이다. 이어 두 brick의 LIN/TL 직렬·TBB +전체 low 필드 일치와 실제 nonzero low를 확인한 focused1/1 검사도0/0이다. +최초 넓은 filter가 전체 원본 smoke를 선택해 중단한 실행은 CTest exit−1로 보존하며 +수렴 실패 또는 통과로 해석하지 않는다. 전체101 endpoint·공식 HDF5·원본 비교는 아직 없다. + ## HEX8-NR-001 — 영 회전 경계와 필수 원본 충돌 검토 당시 PRD008/I/O는 BC/CLOAD DOF4–6을 값0이어도 거부했다. 그러나 모든 필수 INP에 diff --git a/docs/hex8-structural-analysis/physics-review.md b/docs/hex8-structural-analysis/physics-review.md index 38af2e1..1bb2181 100644 --- a/docs/hex8-structural-analysis/physics-review.md +++ b/docs/hex8-structural-analysis/physics-review.md @@ -35,7 +35,30 @@ temporal 실행은4/5 PASS·CTest exit8이며 실패는 이전 principal trace a owned/shared source·원본 두 INP·binary·patch hash를 대조했다. Reviewer는 관련22개 raw와 7개 owned source·patch를 별도로 대조했다. -전체 high+low endpoint 배열의 독립 gate 재계산과 원본 prefix/101 endpoint·필수 reference -비교는 아직 수행하지 않았다. Rounded observable 불변량 검사와 전체 pair의 production +첫32개 검토 이후 원본 prefix의 전체 pair를 아래와 같이 별도 재계산했다. 전체101 endpoint· +필수 reference 비교는 아직 없다. Rounded observable 불변량 검사와 전체 pair의 production gate 통과를 구분한다. Alpha−0.05의 장시간 에너지, warped dynamic mesh, near-incompressibility 및 큰 strain StVK 적용은 이번 evidence로 보증하지 않는다. + +## 원본 초기 상태와 세 증분의 독립 물리 계산 + +변경하지 않은 LIN/TL 두 원본의 초기 상태와 증분1–3에 대해 reviewer가 모든 high+low를 +exact Fraction으로 읽었다. 기준 CSR M·reference X·실제 full force를 사용해 T/W/balance, +p/pdot/L/Ldot, reference/current/free-residual moments와 global force/moment 수지를 +총240개 대조했다. 실제 pair Add의 discarded-tail bound를 전파한 별도 절대 graph budget +내였으며 재생 graph의 high/low도240/240 일치했다. 구속 u/v/a/A와 support work는 exact0다. +Source442절점·300연결·51구속·하중·원래 h/T/INC와 test-only 관측 종료 조건도 대조했다. +선형 U와 독립0.5uᵀf의 최대 절대차는1.30e-29였다. + +Captured F64 M의 ordered coefficient에는7752개의 비대칭 entry가 있고 +최대|Mij−Mji|=2.84217e-14다. 이에 따른 exact Σv×Mv의 작은 비영 값은512eps same-unit +coefficient-roundoff budget 이내이며 이 비대칭으로 정확히 설명된다. Pair 산술 오류와 +계수 오차를 구분하고 M을 대칭화하거나 모멘트를0으로 보정하지 않는다. Ldot=x×Ma는 +같은 endpoint의 계약 observable로 확인했다. LIN의 current internal moment는 O(u²)로 +비영이며 TL의 reference moment도 별도0 기대량이 아니다. + +Alpha−0.05의 frame별 energy balance는 관측량이지 exact conservation gate가 아니다. +Prefix에서 나타난 balance와 진동에 따른 trapezoid work 감소를 그대로 보존한다. +이 범위에서 물리 blocker는 없다. 재현 자료는 로컬 +`.harness/hex8-evidence/precision/physics-prefix-review/`와 `prefix3-frozen/`다. +TL의 local polynomial U 검산은 별도 frozen reference cache에 대해 추가 확인한다. diff --git a/docs/tracking/findings.md b/docs/tracking/findings.md index 206cfc7..dbc5472 100644 --- a/docs/tracking/findings.md +++ b/docs/tracking/findings.md @@ -3,7 +3,7 @@ ## Hex8 동해석의 엄격 잔차와 내부 상태 정밀도 - 조건·현상: Hex8 동적 구현 `f661864`에서 `reference/hexa/hexa_dynlinear.inp`와 - `hexa_dynnonlinear.inp`를 실행하면 세 번째 증분100회 제한에서 실패한다. 최신 원본 + `hexa_dynnonlinear.inp`를 실행하면 세 번째 증분100회 제한에서 실패했다. 이 과거 원본 physical residual 비율은 `1.77739e-8`·`8.73572e-9`, 기준은 `1e-10`이다. - 영향: 원본101개 endpoint의 필수 검증과 동해석 완료 조건을 충족하지 못했다. 현재 전체 요구와 checklist는 [기능 PRD](../hex8-structural-analysis/PRD.md)에 있다. @@ -12,9 +12,10 @@ Float64 해의 부재를 증명하지 않으며, 수렴 기준·원본을 변경할 근거로 쓰지 않는다. - 후속 방향: [실패 근거](../hex8-structural-analysis/lessons.md)를 보존하고, 전체 물리 상태의 high+low 계산과 Float64 출력의 반올림 의미를 정하는 - [승인된 변경안](../hex8-structural-analysis/precision-proposal.md)을 구현한다. - 사용자가 2026-10-09 선택 A로 승인했다. 산술·요소의 첫 native 검사는 통과했으나 - 전체 동적 state·시간수렴·원본 endpoint 검증은 남아 있다. Acceptance 기준은 유지한다. + [승인된 변경안](../hex8-structural-analysis/precision-proposal.md)을 구현했다. + 사용자가 2026-10-09 선택 A로 승인했고 소규모 상태·시간수렴·수지와 두 원본의 세 증분 + 관측을 검증했다. 독립 full-pair 재계산에서 원래 모든 free의 strict gate를 통과했다. + 전체101 endpoint·공식 출력·필수 원본 비교는 남아 있으며 acceptance 기준은 유지한다. ## HHT shell의 시간수렴 diff --git a/include/fesa/analysis/solid_dynamic_analysis_state.h b/include/fesa/analysis/solid_dynamic_analysis_state.h new file mode 100644 index 0000000..d6cdfb9 --- /dev/null +++ b/include/fesa/analysis/solid_dynamic_analysis_state.h @@ -0,0 +1,203 @@ +#ifndef FESA_ANALYSIS_SOLID_DYNAMIC_ANALYSIS_STATE_H_ +#define FESA_ANALYSIS_SOLID_DYNAMIC_ANALYSIS_STATE_H_ + +#include +#include +#include +#include +#include +#include + +#include "fesa/fem/dof_manager.h" +#include "fesa/results/solid_result_records.h" + +namespace fesa { + +/// @brief Dimensionally explicit translation-only Newton evidence. +struct SolidDynamicIterationEvidence { + std::uint32_t iteration_index{0}; + double residual_numerator{0}; // Force infinity norm. + double residual_scale{0}; // Uncancelled force infinity scale. + double residual_metric{0}; // Maximum componentwise force ratio. + std::size_t worst_full_dof{0}; + std::array worst_force_terms{}; // r,Ma,fint,F,uncancelled scale. + double weighted_residual_numerator{0}; // Original HHT force infinity norm. + double recombined_residual_numerator{ + 0}; // Actual solved weighted force norm. + double old_residual_numerator{0}; // Previous physical free force norm. + double weighted_equivalence_metric{0}; // Uncancelled force identity ratio. + std::array + worst_weighted_terms{}; // Original,recombined,old r,scale. + bool zero_scale_exact{true}; + std::array correction_metrics{}; // du/L,h*dv/L,h²*da/L,h²*dA/L. + double update_metric{0}; // Maximum uncancelled update equality ratio. + bool residual_pass{false}; + bool correction_pass{false}; + bool update_pass{false}; +}; + +/// @brief Failure-only shadow dot products; these never change force or +/// acceptance. +struct SolidDynamicLinearDotEvidence { + EntityIndex element_index{0}; + std::size_t local_row{0}; + double absolute_product_sum{0}; + double kernel_force{0}; + double compensated_force{ + 0}; // Shadow using rounded relative differences only. +}; + +/// @brief Failure-only coordinate arithmetic evidence, with length units. +struct SolidDynamicCoordinateEvidence { + EntityIndex element_index{0}; + std::size_t local_row{0}; + double displacement{0}; + double previous_displacement{0}; + double fresh_correction{0}; + double reference_displacement{0}; + double previous_reference_displacement{0}; + double relative_high{0}; + double previous_relative_high{0}; +}; + +/// @brief Failure-only arithmetic shadow; never substitutes physical forces. +struct SolidDynamicForceSumEvidence { + double assembled_internal{0}; + double compensated_internal{0}; + double absolute_local_internal_sum{0}; + double assembled_inertia{0}; + double compensated_inertia{0}; + double absolute_inertia_product_sum{0}; +}; + +/// @brief Internal exact-candidate recurrence evidence, separate from gates. +struct SolidDynamicNewtonRepeatEvidence { + std::size_t history_entry_count{0}; + std::uint32_t matched_entry_iteration{0}; + std::size_t changed_translations{0}; + std::size_t changed_holding_mask{0}; + std::size_t changed_correction_metrics{0}; + bool repeat{false}; +}; + +/// @brief Keeps failed attempt evidence outside the accepted physical history. +struct SolidDynamicIterationRecord { + std::uint32_t increment_index{0}; + double time{0}; + std::string stage; + std::string outcome; + std::string failure_code; + SolidDynamicIterationEvidence evidence; + std::vector failure_linear_dot; + std::vector failure_coordinates; + std::optional failure_force_sum; + std::optional newton_repeat; +}; + +/// @brief Retains normalized physical low parts; never serialized as extra +/// result fields. Corresponding endpoint fields contain the nearest Float64 +/// high parts. Both parts are copied and validated together at commit. +struct SolidDynamicFrameLowParts { + Vector displacement{0}; + Vector velocity{0}; + Vector acceleration{0}; + Vector algorithmic_acceleration{0}; + Vector inertia_force{0}; + Vector internal_force{0}; + Vector external_force{0}; + Vector residual{0}; + Vector reaction{0}; + Vector residual_scale{0}; + double strain_energy{0}; + double kinetic_energy{0}; + double external_work{0}; + double energy_balance_error{0}; + std::array linear_momentum{}; + std::array angular_momentum{}; + std::array linear_momentum_rate{}; + std::array angular_momentum_rate{}; + std::array global_equilibrium{}; + std::array linearization_current_internal_moment{}; + std::array reference_internal_moment{}; + std::array current_free_residual_moment{}; +}; + +/// @brief Owns one complete Solid endpoint; physical a differs from auxiliary +/// A. +struct SolidDynamicAcceptedFrameState { + std::uint32_t frame_index{0}; + std::uint32_t increment_index{0}; + std::uint32_t iteration_count{0}; + double time{0}; + bool is_initial{false}; + std::vector nodal_internal_order; + std::vector nodal_source_ids; + Vector displacement{0}; + Vector velocity{0}; + Vector acceleration{0}; + Vector algorithmic_acceleration{0}; + Vector inertia_force{0}; + Vector internal_force{0}; + Vector external_force{0}; + Vector residual{0}; + Vector reaction{0}; + Vector residual_scale{0}; + std::vector solid_results; + double strain_energy{0}; + double kinetic_energy{0}; + double external_work{0}; + double energy_balance_error{0}; + std::array linear_momentum{}; + std::array angular_momentum{}; + std::array linear_momentum_rate{}; + std::array angular_momentum_rate{}; + std::array global_equilibrium{}; + // Internal evidence, not additional official output fields. + std::array linearization_current_internal_moment{}; + std::array reference_internal_moment{}; + std::array current_free_residual_moment{}; + std::vector iterations; + SolidDynamicFrameLowParts low; +}; + +/// @brief Owns only validated accepted Solid history; rejects partial +/// candidates. +class SolidDynamicAnalysisState { + public: + /// @brief Creates an empty state for exactly one zero-BC fixed dynamic + /// schedule. + static Result Create(const Domain& domain, + const DofManager& dofs); + /// @brief Validates shape, identity, schedule, finite values and Newton + /// evidence. + /// @note A failed candidate never changes accepted frames or current state. + Status CommitAcceptedFrame(const SolidDynamicAcceptedFrameState& frame); + /// @brief Returns an independent accepted endpoint copy for trial + /// preparation. + Result BeginTrialFromCommitted() const; + /// @brief Returns accepted endpoints in source increment order. + const std::vector& AcceptedFrames() + const noexcept; + /// @brief Returns the last endpoint; throws when no frame is accepted. + const SolidDynamicAcceptedFrameState& CommittedState() const; + /// @brief Reports whether schedule and inventory have been initialized. + bool IsInitialized() const noexcept; + /// @brief Appends attempt evidence without committing physical state. + void RecordIteration(SolidDynamicIterationRecord record); + /// @brief Returns attempt evidence, including a terminal failure if present. + const std::vector& IterationRecords() + const noexcept; + + private: + std::size_t full_dof_count_{0}; + std::uint32_t increment_count_{0}; + double time_increment_{0}; + std::vector constrained_dofs_; + std::vector source_ids_; + std::vector point_inventory_; + std::vector accepted_frames_; + std::vector iteration_records_; + bool initialized_{false}; +}; +} // namespace fesa +#endif // FESA_ANALYSIS_SOLID_DYNAMIC_ANALYSIS_STATE_H_ diff --git a/include/fesa/analysis/solid_hht_dynamic_procedure.h b/include/fesa/analysis/solid_hht_dynamic_procedure.h new file mode 100644 index 0000000..b4a5e56 --- /dev/null +++ b/include/fesa/analysis/solid_hht_dynamic_procedure.h @@ -0,0 +1,44 @@ +#ifndef FESA_ANALYSIS_SOLID_HHT_DYNAMIC_PROCEDURE_H_ +#define FESA_ANALYSIS_SOLID_HHT_DYNAMIC_PROCEDURE_H_ + +#include "fesa/analysis/solid_dynamic_analysis_state.h" + +namespace fesa { +class LinearSolver; +class ParallelFor; +class TangentSolver; + +/// @brief Owns the validated 3N dynamic result independently of serialization. +struct SolidDynamicExecutionResult { + DofManager dofs; + SolidDynamicAnalysisState state; +}; + +/// @brief Integrates linear or TL Solid mechanics with constant consistent +/// mass. +class SolidHhtDynamicProcedureExecutor { + public: + /// @brief Retains caller-owned adapters which must outlive execution. + SolidHhtDynamicProcedureExecutor(const ParallelFor& parallel_for, + LinearSolver& mass_solver, + TangentSolver& tangent_solver); + /// @brief Returns complete accepted history without opening output files. + Result Execute(const Domain& domain) const; + /// @brief Executes into an empty initialized state, preserving accepted + /// history on every failure. Domain must match the state inventory and + /// schedule. + Status ExecuteInto(const Domain& domain, + SolidDynamicAnalysisState& state) const; + /// @brief Internal verification only; alpha0 verifies conservative dynamics. + /// @note This does not extend the approved INP parameter subset. + Status ExecuteIntoForVerification(const Domain& domain, + SolidDynamicAnalysisState& state, + double alpha) const; + + private: + const ParallelFor& parallel_for_; + LinearSolver& mass_solver_; + TangentSolver& tangent_solver_; +}; +} // namespace fesa +#endif // FESA_ANALYSIS_SOLID_HHT_DYNAMIC_PROCEDURE_H_ diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index bd85b6e..58d1627 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -17,7 +17,9 @@ add_library( analysis/nonlinear_static_control.cpp analysis/nonlinear_static_procedure.cpp analysis/prepared_analysis.cpp + analysis/solid_dynamic_analysis_state.cpp analysis/solid_frequency_procedure.cpp + analysis/solid_hht_dynamic_procedure.cpp analysis/solid_static_procedure.cpp analysis/solid_static_state.cpp analysis/thermal_analysis_state.cpp diff --git a/src/fesa/analysis/solid_dynamic_analysis_state.cpp b/src/fesa/analysis/solid_dynamic_analysis_state.cpp new file mode 100644 index 0000000..d5e2cab --- /dev/null +++ b/src/fesa/analysis/solid_dynamic_analysis_state.cpp @@ -0,0 +1,280 @@ +#include "fesa/analysis/solid_dynamic_analysis_state.h" + +#include +#include +#include +#include + +#include "fesa/elements/element_factory.h" +#include "fesa/elements/hex8_element.h" +#include "math/high_low_arithmetic.h" + +namespace fesa { +namespace { +/// @brief Rejects a candidate without changing observable history. +Status Failure(const std::string& message) { + return Status::Failure(FailureCategory::kSolver, + {{Severity::kError, + "solid-dynamic-invalid-state", + {}, + "DYNAMIC", + {}, + message}}); +} +/// @brief Checks a physical record without replacing nonfinite quantities. +template +bool Finite(const std::array& values) { + return std::all_of(values.begin(), values.end(), + [](double x) { return std::isfinite(x); }); +} +bool ValidPair(double high, double low) { + const math_internal::HighLow value(high, low); + return math_internal::IsFinite(value) && math_internal::IsNormalized(value); +} +bool ValidPair(const Vector& high, const Vector& low, std::size_t size) { + if (high.Size() != size || low.Size() != size) return false; + for (std::size_t i = 0; i < size; ++i) + if (!ValidPair(high[i], low[i])) return false; + return true; +} +template +bool ValidPair(const std::array& high, + const std::array& low) { + for (std::size_t i = 0; i < N; ++i) + if (!ValidPair(high[i], low[i])) return false; + return true; +} +/// @brief Preserves canonical and raw source identity independently of indices. +bool SameId(const SourceEntityId& a, const SourceEntityId& b) { + return a.instance_name == b.instance_name && + a.source_label == b.source_label && + a.source_label_text == b.source_label_text; +} +} // namespace + +Result SolidDynamicAnalysisState::Create( + const Domain& domain, const DofManager& dofs) { + const auto* step = domain.Dynamic(); + if (!domain.IsSolid() || step == nullptr || + dofs.Profile() != DofProfile::kSolidStructural || + dofs.FullDofCount() != 3 * domain.Nodes().size() || + !dofs.ValidateInvariants().IsOk() || step->increment_count == 0 || + !std::isfinite(step->initial_increment) || step->initial_increment <= 0 || + !std::isfinite(step->time_period) || + static_cast(step->increment_count) * step->initial_increment != + step->time_period || + step->increment_count > step->maximum_increments) { + return Result::Failure( + Failure("Solid HHT requires a valid 3N exact fixed schedule.")); + } + for (std::size_t i = 0; i < dofs.PrescribedValues().Size(); ++i) + if (dofs.PrescribedValues()[i] != 0) + return Result::Failure( + Failure("Dynamic prescribed translations must be exactly zero.")); + SolidDynamicAnalysisState state; + state.full_dof_count_ = dofs.FullDofCount(); + state.increment_count_ = step->increment_count; + state.time_increment_ = step->initial_increment; + state.constrained_dofs_ = dofs.ConstrainedDofs(); + for (const auto& node : domain.Nodes()) + state.source_ids_.push_back(node.source_id); + ElementFactory factory; + for (std::size_t e = 0; e < domain.Elements().Size(); ++e) { + auto runtime = factory.Create(domain.Elements()[e], domain); + if (!runtime.HasValue()) + return Result::Failure(runtime.GetStatus()); + const auto* solid = dynamic_cast(runtime.Value().get()); + if (solid == nullptr) + return Result::Failure(Failure( + "Dynamic Solid inventory contains a different element family.")); + auto zero = solid->EvaluateLinear(Vector(24)); + if (!zero.HasValue()) + return Result::Failure(zero.GetStatus()); + state.point_inventory_.insert(state.point_inventory_.end(), + zero.Value().points.begin(), + zero.Value().points.end()); + } + state.initialized_ = true; + return Result::Success(std::move(state)); +} + +Status SolidDynamicAnalysisState::CommitAcceptedFrame( + const SolidDynamicAcceptedFrameState& f) { + if (!math_internal::ArithmeticEnvironmentSupported()) + return Failure( + "Paired state commit requires nearest rounding and gradual " + "underflow."); + const std::size_t index = accepted_frames_.size(); + if (!initialized_ || index > increment_count_ || f.frame_index != index || + f.increment_index != index || f.is_initial != (index == 0) || + f.time != static_cast(index) * time_increment_ || + f.iteration_count > 100 || (index == 0 && f.iteration_count != 0) || + f.nodal_internal_order.size() != source_ids_.size() || + f.nodal_source_ids.size() != source_ids_.size() || + f.solid_results.size() != point_inventory_.size()) { + return Failure( + "Incomplete or out-of-order Solid endpoint inventory/schedule."); + } + for (std::size_t n = 0; n < source_ids_.size(); ++n) + if (f.nodal_internal_order[n] != n || + !SameId(f.nodal_source_ids[n], source_ids_[n])) + return Failure("Solid nodal source identity differs from its inventory."); + const std::array high{ + &f.displacement, &f.velocity, + &f.acceleration, &f.algorithmic_acceleration, + &f.inertia_force, &f.internal_force, + &f.external_force, &f.residual, + &f.reaction, &f.residual_scale}; + const std::array low{ + &f.low.displacement, &f.low.velocity, + &f.low.acceleration, &f.low.algorithmic_acceleration, + &f.low.inertia_force, &f.low.internal_force, + &f.low.external_force, &f.low.residual, + &f.low.reaction, &f.low.residual_scale}; + for (std::size_t i = 0; i < high.size(); ++i) + if (!ValidPair(*high[i], *low[i], full_dof_count_)) + return Failure( + "Solid endpoint pairs are incomplete, unnormalized or " + "nonfinite."); + using namespace math_internal; + for (std::size_t i = 0; i < full_dof_count_; ++i) { + const bool constrained = std::binary_search(constrained_dofs_.begin(), + constrained_dofs_.end(), i); + const HighLow inertia(f.inertia_force[i], f.low.inertia_force[i]); + const HighLow internal(f.internal_force[i], f.low.internal_force[i]); + const HighLow external(f.external_force[i], f.low.external_force[i]); + const HighLow residual(f.residual[i], f.low.residual[i]); + const HighLow scale(f.residual_scale[i], f.low.residual_scale[i]); + const HighLow reaction(f.reaction[i], f.low.reaction[i]); + const auto expected_residual = Subtract(Add(inertia, internal), external); + const auto expected_scale = + Add(Add(Abs(inertia), Abs(internal)), Abs(external)); + if (!IsFinite(expected_residual) || !IsFinite(expected_scale) || + Compare(residual, expected_residual) != 0 || + Compare(scale, expected_scale) != 0 || + Compare(reaction, constrained ? residual : HighLow(0)) != 0) + return Failure( + "Physical reaction/residual/uncancelled force scale disagree."); + if (constrained && (f.displacement[i] != 0 || f.low.displacement[i] != 0 || + f.velocity[i] != 0 || f.low.velocity[i] != 0 || + f.acceleration[i] != 0 || f.low.acceleration[i] != 0 || + f.algorithmic_acceleration[i] != 0 || + f.low.algorithmic_acceleration[i] != 0)) + return Failure( + "Prescribed dynamic translations and rates must be exact zero."); + if (!constrained) { + if (Compare(scale, HighLow(0)) < 0) + return Failure("Negative physical force scale."); + if (Compare(scale, HighLow(0)) == 0) { + if (Compare(residual, HighLow(0)) != 0) + return Failure("Physical zero-scale residual must be exactly zero."); + } else { + const auto ratio = Divide(Abs(residual), scale); + if (!IsFinite(ratio) || !IsNormalized(ratio) || + Compare(ratio, HighLow(1e-10)) > 0) + return Failure("Physical free endpoint equilibrium is incomplete."); + } + } + if (index == 0 && + (f.displacement[i] != 0 || f.velocity[i] != 0 || + f.low.displacement[i] != 0 || f.low.velocity[i] != 0 || + Compare(HighLow(f.acceleration[i], f.low.acceleration[i]), + HighLow(f.algorithmic_acceleration[i], + f.low.algorithmic_acceleration[i])) != 0)) + return Failure( + "Initial Solid displacement/velocity must be zero and A=a."); + } + for (std::size_t i = 0; i < point_inventory_.size(); ++i) { + const auto& row = f.solid_results[i]; + const auto& expected = point_inventory_[i]; + if (row.element != expected.element || + row.integration_point != expected.integration_point || + row.natural_coordinates != expected.natural_coordinates || + row.reference_coordinates != expected.reference_coordinates || + !Finite(row.cauchy_stress) || !Finite(row.strain)) + return Failure( + "Solid point recovery has invalid identity or finite values."); + } + if (!ValidPair(f.strain_energy, f.low.strain_energy) || + !ValidPair(f.kinetic_energy, f.low.kinetic_energy) || + !ValidPair(f.external_work, f.low.external_work) || + !ValidPair(f.energy_balance_error, f.low.energy_balance_error)) + return Failure("Nonfinite or unnormalized Solid energy/work pair."); + if ((index == 0 && + (f.external_work != 0 || f.low.external_work != 0 || + f.energy_balance_error != 0 || f.low.energy_balance_error != 0)) || + Compare(HighLow(f.kinetic_energy, f.low.kinetic_energy), HighLow(0)) < 0) + return Failure( + "Initial work/error or kinetic energy violates its definition."); + if (!ValidPair(f.linear_momentum, f.low.linear_momentum) || + !ValidPair(f.angular_momentum, f.low.angular_momentum) || + !ValidPair(f.linear_momentum_rate, f.low.linear_momentum_rate) || + !ValidPair(f.angular_momentum_rate, f.low.angular_momentum_rate) || + !ValidPair(f.linearization_current_internal_moment, + f.low.linearization_current_internal_moment) || + !ValidPair(f.reference_internal_moment, + f.low.reference_internal_moment) || + !ValidPair(f.current_free_residual_moment, + f.low.current_free_residual_moment)) + return Failure("Nonfinite or unnormalized physical momentum/moment pair."); + if (!ValidPair(f.global_equilibrium, f.low.global_equilibrium) || + f.iterations.empty()) + return Failure("Missing finite equilibrium or iteration evidence."); + for (std::size_t i = 0; i < f.iterations.size(); ++i) { + const auto& e = f.iterations[i]; + if (e.iteration_index != i || !Finite(e.correction_metrics) || + !Finite(e.worst_force_terms) || !Finite(e.worst_weighted_terms)) + return Failure("Incomplete or nonfinite Newton correction evidence."); + for (double x : {e.residual_numerator, e.residual_scale, e.residual_metric, + e.weighted_residual_numerator, e.update_metric, + e.old_residual_numerator, e.recombined_residual_numerator, + e.weighted_equivalence_metric}) + if (!std::isfinite(x) || x < 0) + return Failure("Invalid raw Newton residual/update evidence."); + for (double x : e.correction_metrics) + if (x < 0) return Failure("Invalid correction scale."); + } + const auto& last = f.iterations.back(); + if (f.iterations.size() != static_cast(f.iteration_count) + 1 || + !last.residual_pass || !last.correction_pass || !last.update_pass || + !last.zero_scale_exact || last.residual_metric > 1e-10 || + last.update_metric > 1e-12 || last.weighted_equivalence_metric > 1e-12 || + std::any_of(last.correction_metrics.begin(), + last.correction_metrics.end(), + [](double x) { return x > 1e-10; })) + return Failure( + "Candidate does not carry complete accepted convergence evidence."); + accepted_frames_.push_back(f); + return Status::Ok(); +} + +Result +SolidDynamicAnalysisState::BeginTrialFromCommitted() const { + if (accepted_frames_.empty()) + return Result::Failure( + Failure("No accepted Solid endpoint exists.")); + return Result::Success( + accepted_frames_.back()); +} +const std::vector& +SolidDynamicAnalysisState::AcceptedFrames() const noexcept { + return accepted_frames_; +} +const SolidDynamicAcceptedFrameState& +SolidDynamicAnalysisState::CommittedState() const { + if (accepted_frames_.empty()) + throw std::logic_error("No accepted Solid endpoint."); + return accepted_frames_.back(); +} +bool SolidDynamicAnalysisState::IsInitialized() const noexcept { + return initialized_; +} +void SolidDynamicAnalysisState::RecordIteration( + SolidDynamicIterationRecord record) { + iteration_records_.push_back(std::move(record)); +} +const std::vector& +SolidDynamicAnalysisState::IterationRecords() const noexcept { + return iteration_records_; +} +} // namespace fesa diff --git a/src/fesa/analysis/solid_hht_dynamic_procedure.cpp b/src/fesa/analysis/solid_hht_dynamic_procedure.cpp new file mode 100644 index 0000000..db5b9dc --- /dev/null +++ b/src/fesa/analysis/solid_hht_dynamic_procedure.cpp @@ -0,0 +1,870 @@ +#include "fesa/analysis/solid_hht_dynamic_procedure.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/hht_time_integration.h" +#include "fesa/assembly/load_assembler.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/constraints/essential_constraint_policy.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/hex8_element.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/solvers/linear/linear_solver.h" +#include "fesa/solvers/linear/tangent_solver.h" +#include "solid_hht_force_identity.h" +#include "solid_hht_high_low.h" + +namespace fesa { +namespace { +using math_internal::Abs; +using math_internal::Add; +using math_internal::Compare; +using math_internal::Divide; +using math_internal::HighLow; +using math_internal::IsFinite; +using math_internal::IsNormalized; +using math_internal::Multiply; +using math_internal::Negate; +using math_internal::Rounded; +using math_internal::Subtract; + +HighLow Pair(const Vector& high, const Vector& low, std::size_t i) { + return {high[i], low[i]}; +} +void Store(HighLow value, double* high, double* low) { + *high = value.high; + *low = value.low; +} +void Store(HighLow value, Vector* high, Vector* low, std::size_t i) { + Store(value, &(*high)[i], &(*low)[i]); +} +bool Valid(HighLow value) { return IsFinite(value) && IsNormalized(value); } +HighLow Maximum(HighLow left, HighLow right) { + if (!Valid(left) || !Valid(right)) return math_internal::InvalidHighLow(); + return Compare(left, right) >= 0 ? left : right; +} +void InitializeLowVectors(SolidDynamicAcceptedFrameState* f, std::size_t size) { + for (auto* v : {&f->low.displacement, &f->low.velocity, &f->low.acceleration, + &f->low.algorithmic_acceleration, &f->low.inertia_force, + &f->low.internal_force, &f->low.external_force, + &f->low.residual, &f->low.reaction, &f->low.residual_scale}) + *v = Vector(size); +} +struct PairedTangent { + std::vector values; +}; +/// @brief Reports source-context procedure failure without committing a trial. +Status Failure(const Domain& domain, const char* code, + const std::string& message) { + return Status::Failure( + FailureCategory::kSolver, + {{Severity::kError, code, + domain.Dynamic() == nullptr ? SourceLocation{} + : domain.Dynamic()->dynamic_location, + "DYNAMIC", + domain.Dynamic() == nullptr ? "" + : domain.Dynamic()->source_step_name_raw, + message}}); +} +Status CheckArithmeticEnvironment(const Domain& domain) { + return math_internal::ArithmeticEnvironmentSupported() + ? Status::Ok() + : Failure(domain, "solid-dynamic-arithmetic-environment", + "Paired HHT requires nearest rounding and gradual " + "underflow, including after backend callbacks."); +} +/// @brief Owns immutable kernels and the stable full-space translation scatter. +struct RuntimeSolid { + std::unique_ptr owner; + const Hex8Element* solid; + std::vector scatter; +}; +/// @brief Creates checked kernels in declaration order, with no mutable element +/// state. +Result> BuildSolids(const Domain& domain, + const DofManager& dofs) { + std::vector solids; + ElementFactory factory; + for (std::size_t e = 0; e < domain.Elements().Size(); ++e) { + auto owner = factory.Create(domain.Elements()[e], domain); + if (!owner.HasValue()) + return Result>::Failure(owner.GetStatus()); + const auto* solid = dynamic_cast(owner.Value().get()); + if (solid == nullptr) + return Result>::Failure( + Failure(domain, "solid-dynamic-element-family", + "HHT requires exclusively Hex8 kernels.")); + auto scatter = dofs.ElementScatter(solid->DofLayout()); + if (!scatter.HasValue()) + return Result>::Failure(scatter.GetStatus()); + solids.push_back( + {std::move(owner.Value()), solid, std::move(scatter.Value())}); + } + return Result>::Success(std::move(solids)); +} +/// @brief Preserves all 576 structural entries in stable element/local order. +void AppendMatrix(const Matrix& local, const RuntimeSolid& solid, std::size_t e, + std::vector* coo) { + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) + coo->push_back( + {solid.scatter[i], solid.scatter[j], local(i, j), e, 24 * i + j}); +} +/// @brief Assembles constant reference consistent mass once without lumping. +Result AssembleMass(const std::vector& solids, + const DofManager& dofs) { + std::vector coo; + coo.reserve(solids.size() * 576); + for (std::size_t e = 0; e < solids.size(); ++e) { + auto mass = solids[e].solid->ComputeMass(); + if (!mass.HasValue()) + return Result::Failure(mass.GetStatus()); + AppendMatrix(mass.Value().values, solids[e], e, &coo); + } + return SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(), + std::move(coo), dofs.GetSparsePattern()); +} +/// @brief Evaluates independent local trials, then deterministically reduces +/// recovery. +Result EvaluateSolids(const Domain& domain, + const std::vector& solids, + const DofManager& dofs, + const ParallelFor& parallel, + SolidDynamicAcceptedFrameState* f) { + const auto environment = CheckArithmeticEnvironment(domain); + if (!environment.IsOk()) return Result::Failure(environment); + if (!detail::ValidSolidHhtPair(f->displacement, f->low.displacement, + dofs.FullDofCount())) + return Result::Failure( + Failure(domain, "solid-dynamic-invalid-paired-state", + "Local evaluation requires complete normalized pairs.")); + std::vector>> responses( + solids.size()); + parallel.Execute(solids.size(), [&](std::size_t e) { + Vector local(24), low(24); + for (std::size_t i = 0; i < 24; ++i) { + local[i] = f->displacement[solids[e].scatter[i]]; + low[i] = f->low.displacement[solids[e].scatter[i]]; + } + responses[e] = domain.Dynamic()->nlgeom + ? solids[e].solid->EvaluateNonlinearHighLow(local, low) + : solids[e].solid->EvaluateLinearHighLow(local, low); + }); + const auto reduction_environment = CheckArithmeticEnvironment(domain); + if (!reduction_environment.IsOk()) + return Result::Failure(reduction_environment); + const auto& pattern = dofs.GetSparsePattern(); + PairedTangent tangent; + tangent.values.resize(pattern.column_indices.size()); + f->internal_force = Vector(dofs.FullDofCount()); + f->low.internal_force = Vector(dofs.FullDofCount()); + f->strain_energy = 0; + f->low.strain_energy = 0; + f->solid_results.clear(); + for (std::size_t e = 0; e < solids.size(); ++e) { + if (!responses[e]) + return Result::Failure( + Failure(domain, "solid-dynamic-local-execution", + "Local worker did not evaluate its index.")); + if (!responses[e]->HasValue()) + return Result::Failure(responses[e]->GetStatus()); + const auto& r = responses[e]->Value(); + for (std::size_t i = 0; i < 24; ++i) { + const auto row = solids[e].scatter[i]; + const auto force = + Add(Pair(f->internal_force, f->low.internal_force, row), + Pair(r.high.internal_force, r.internal_force_low, i)); + if (!Valid(force)) + return Result::Failure( + Failure(domain, "solid-dynamic-nonfinite-assembly", + "Paired internal force reduction failed.")); + Store(force, &f->internal_force, &f->low.internal_force, row); + for (std::size_t j = 0; j < 24; ++j) { + const auto column = solids[e].scatter[j]; + const auto begin = + pattern.column_indices.begin() + pattern.row_offsets[row]; + const auto end = + pattern.column_indices.begin() + pattern.row_offsets[row + 1]; + const auto found = std::lower_bound(begin, end, column); + if (found == end || *found != column) + return Result::Failure( + Failure(domain, "solid-dynamic-sparse-pattern", + "Local tangent is outside the fixed sparse pattern.")); + const auto k = + static_cast(found - pattern.column_indices.begin()); + tangent.values[k] = + Add(tangent.values[k], + HighLow(r.high.tangent(i, j), r.tangent_low(i, j))); + if (!Valid(tangent.values[k])) + return Result::Failure( + Failure(domain, "solid-dynamic-nonfinite-assembly", + "Paired tangent reduction failed.")); + } + } + const auto energy = Add(HighLow(f->strain_energy, f->low.strain_energy), + HighLow(r.high.strain_energy, r.strain_energy_low)); + if (!Valid(energy)) + return Result::Failure( + Failure(domain, "solid-dynamic-nonfinite-assembly", + "Paired strain energy reduction failed.")); + Store(energy, &f->strain_energy, &f->low.strain_energy); + f->solid_results.insert(f->solid_results.end(), r.high.points.begin(), + r.high.points.end()); + } + return Result::Success(std::move(tangent)); +} +/// @brief Returns the finite positive reference nodal bounding-box diagonal. +Result ReferenceLength(const Domain& domain) { + auto lo = domain.Nodes().front().coordinates; + auto hi = lo; + for (const auto& node : domain.Nodes()) + for (std::size_t c = 0; c < 3; ++c) { + lo[c] = (std::min)(lo[c], node.coordinates[c]); + hi[c] = (std::max)(hi[c], node.coordinates[c]); + } + const double length = + std::hypot(std::hypot(hi[0] - lo[0], hi[1] - lo[1]), hi[2] - lo[2]); + if (!std::isfinite(length) || length <= 0) + return Result::Failure( + Failure(domain, "solid-dynamic-reference-length", + "Reference length is not finite positive.")); + return Result::Success(length); +} +/// @brief Forms M/(beta*h²)+(1+alpha)*Kt in the unchanged full sparse pattern. +Result EffectiveTangent(const SparseMatrix& m, + const PairedTangent& k, + const DofManager& dofs, double h, + const HhtParameters& p) { + std::vector coo; + coo.reserve(m.Values().size()); + const auto denominator = + Multiply(Multiply(HighLow(p.beta), HighLow(h)), HighLow(h)); + const auto coefficient = Add(HighLow(1), HighLow(p.alpha)); + for (std::size_t i = 0; i < m.Rows(); ++i) + for (std::size_t j = m.RowOffsets()[i]; j < m.RowOffsets()[i + 1]; ++j) { + const auto value = Add(Divide(HighLow(m.Values()[j]), denominator), + Multiply(coefficient, k.values[j])); + if (!Valid(value)) + return Result::Failure( + detail::SolidHhtPairedKinematicsFailure()); + coo.push_back({i, m.ColumnIndices()[j], Rounded(value), 0, j}); + } + return SparseMatrix::FromCoo(m.Rows(), m.Columns(), std::move(coo), + dofs.GetSparsePattern()); +} + +using Paired3 = std::array; +void CrossAdd(const Paired3& x, const Paired3& force, Paired3* sum) { + for (std::size_t c = 0; c < 3; ++c) { + const auto j = (c + 1) % 3; + const auto k = (c + 2) % 3; + (*sum)[c] = Add((*sum)[c], Subtract(Multiply(x[j], force[k]), + Multiply(x[k], force[j]))); + } +} +bool Store3(const Paired3& values, std::array* high, + std::array* low) { + for (std::size_t c = 0; c < 3; ++c) { + if (!Valid(values[c])) return false; + Store(values[c], &(*high)[c], &(*low)[c]); + } + return true; +} + +/// @brief Recovers full physical endpoint quantities before any publication or +/// commit. The constant Float64 mass coefficients are lifted exactly once. +Status RecoverEndpoint(const Domain& domain, const DofManager& dofs, + const SparseMatrix& mass, + const SolidDynamicAcceptedFrameState* old, + HighLow initial_energy, + SolidDynamicAcceptedFrameState* f, + SolidDynamicIterationEvidence* e) { + auto inertia = detail::MultiplySolidHhtMassHighLow(mass, f->acceleration, + f->low.acceleration); + auto momentum = + detail::MultiplySolidHhtMassHighLow(mass, f->velocity, f->low.velocity); + if (!inertia.HasValue()) return inertia.GetStatus(); + if (!momentum.HasValue()) return momentum.GetStatus(); + f->inertia_force = std::move(inertia.Value().high); + f->low.inertia_force = std::move(inertia.Value().low); + const auto size = dofs.FullDofCount(); + for (auto* v : {&f->residual, &f->low.residual, &f->reaction, + &f->low.reaction, &f->residual_scale, &f->low.residual_scale}) + *v = Vector(size); + HighLow numerator, scale_maximum, metric; + for (std::size_t i = 0; i < size; ++i) { + const auto ma = Pair(f->inertia_force, f->low.inertia_force, i); + const auto internal = Pair(f->internal_force, f->low.internal_force, i); + const auto external = Pair(f->external_force, f->low.external_force, i); + const auto residual = Subtract(Add(ma, internal), external); + const auto scale = Add(Add(Abs(ma), Abs(internal)), Abs(external)); + if (!Valid(residual) || !Valid(scale)) + return Failure(domain, "solid-dynamic-nonfinite-residual", + "Physical paired force arithmetic failed."); + Store(residual, &f->residual, &f->low.residual, i); + Store(scale, &f->residual_scale, &f->low.residual_scale, i); + } + for (auto i : dofs.ConstrainedDofs()) { + f->reaction[i] = f->residual[i]; + f->low.reaction[i] = f->low.residual[i]; + } + for (auto i : dofs.FreeDofs()) { + const auto residual = Pair(f->residual, f->low.residual, i); + const auto scale = Pair(f->residual_scale, f->low.residual_scale, i); + numerator = Maximum(numerator, Abs(residual)); + scale_maximum = Maximum(scale_maximum, scale); + if (Compare(scale, HighLow(0)) == 0) { + e->zero_scale_exact = + e->zero_scale_exact && Compare(residual, HighLow(0)) == 0; + } else { + const auto ratio = Divide(Abs(residual), scale); + if (!Valid(ratio)) + return Failure(domain, "solid-dynamic-nonfinite-residual", + "Physical paired force ratio is unsupported."); + if (Compare(ratio, metric) > 0) { + metric = ratio; + e->worst_full_dof = i; + e->worst_force_terms = { + Rounded(residual), + Rounded(Pair(f->inertia_force, f->low.inertia_force, i)), + Rounded(Pair(f->internal_force, f->low.internal_force, i)), + Rounded(Pair(f->external_force, f->low.external_force, i)), + Rounded(scale)}; + } + } + } + if (!Valid(numerator) || !Valid(scale_maximum) || !Valid(metric)) + return Failure(domain, "solid-dynamic-nonfinite-residual", + "Physical paired evidence is nonfinite."); + e->residual_numerator = Rounded(numerator); + e->residual_scale = Rounded(scale_maximum); + e->residual_metric = Rounded(metric); + e->residual_pass = + e->zero_scale_exact && Compare(metric, HighLow(1e-10)) <= 0; + + HighLow kinetic; + HighLow work = old == nullptr + ? HighLow(0) + : HighLow(old->external_work, old->low.external_work); + for (std::size_t i = 0; i < size; ++i) { + kinetic = + Add(kinetic, + Multiply(Pair(f->velocity, f->low.velocity, i), + Pair(momentum.Value().high, momentum.Value().low, i))); + if (old != nullptr) { + const auto average = + Multiply(HighLow(.5), + Add(Pair(old->external_force, old->low.external_force, i), + Pair(f->external_force, f->low.external_force, i))); + const auto delta = + Subtract(Pair(f->displacement, f->low.displacement, i), + Pair(old->displacement, old->low.displacement, i)); + work = Add(work, Multiply(average, delta)); + } + } + kinetic = Multiply(HighLow(.5), kinetic); + const auto error = Subtract( + Subtract(Add(kinetic, HighLow(f->strain_energy, f->low.strain_energy)), + initial_energy), + work); + if (!Valid(kinetic) || !Valid(work) || !Valid(error)) + return Failure(domain, "solid-dynamic-nonfinite-energy", + "Paired energy or endpoint work is unsupported."); + Store(kinetic, &f->kinetic_energy, &f->low.kinetic_energy); + Store(work, &f->external_work, &f->low.external_work); + Store(error, &f->energy_balance_error, &f->low.energy_balance_error); + + Paired3 linear{}, angular{}, linear_rate{}, angular_rate{}, applied_moment{}, + linear_equilibrium{}, current_internal_moment{}, reference_moment{}, + free_residual_moment{}; + for (std::size_t n = 0; n < domain.Nodes().size(); ++n) { + Paired3 x{}, reference{}, mv{}, ma{}, applied{}, internal{}, + free_residual{}, u{}; + for (std::size_t c = 0; c < 3; ++c) { + const auto i = 3 * n + c; + reference[c] = HighLow(domain.Nodes()[n].coordinates[c]); + u[c] = Pair(f->displacement, f->low.displacement, i); + x[c] = Add(reference[c], u[c]); + mv[c] = Pair(momentum.Value().high, momentum.Value().low, i); + ma[c] = Pair(f->inertia_force, f->low.inertia_force, i); + internal[c] = Pair(f->internal_force, f->low.internal_force, i); + applied[c] = Add(Pair(f->external_force, f->low.external_force, i), + Pair(f->reaction, f->low.reaction, i)); + free_residual[c] = Subtract(Pair(f->residual, f->low.residual, i), + Pair(f->reaction, f->low.reaction, i)); + linear[c] = Add(linear[c], mv[c]); + linear_rate[c] = Add(linear_rate[c], ma[c]); + linear_equilibrium[c] = + Add(linear_equilibrium[c], Subtract(applied[c], ma[c])); + } + CrossAdd(x, mv, &angular); + CrossAdd(x, ma, &angular_rate); + CrossAdd(x, applied, &applied_moment); + CrossAdd(u, internal, ¤t_internal_moment); + CrossAdd(reference, internal, &reference_moment); + CrossAdd(x, free_residual, &free_residual_moment); + } + if (!Store3(linear, &f->linear_momentum, &f->low.linear_momentum) || + !Store3(angular, &f->angular_momentum, &f->low.angular_momentum) || + !Store3(linear_rate, &f->linear_momentum_rate, + &f->low.linear_momentum_rate) || + !Store3(angular_rate, &f->angular_momentum_rate, + &f->low.angular_momentum_rate) || + !Store3(current_internal_moment, + &f->linearization_current_internal_moment, + &f->low.linearization_current_internal_moment) || + !Store3(reference_moment, &f->reference_internal_moment, + &f->low.reference_internal_moment) || + !Store3(free_residual_moment, &f->current_free_residual_moment, + &f->low.current_free_residual_moment)) + return Failure( + domain, "solid-dynamic-nonfinite-momentum", + "Paired momentum or current/reference moment is unsupported."); + for (std::size_t c = 0; c < 3; ++c) { + const auto moment = Subtract(applied_moment[c], angular_rate[c]); + if (!Valid(linear_equilibrium[c]) || !Valid(moment)) + return Failure(domain, "solid-dynamic-nonfinite-equilibrium", + "Paired global equilibrium is unsupported."); + Store(linear_equilibrium[c], &f->global_equilibrium[c], + &f->low.global_equilibrium[c]); + Store(moment, &f->global_equilibrium[3 + c], + &f->low.global_equilibrium[3 + c]); + } + return Status::Ok(); +} + +/// @brief Uses full paired same-unit uncancelled formula terms. +HighLow UpdateMetric(HighLow left, std::initializer_list terms) { + HighLow sum; + auto scale = Abs(left); + for (const auto x : terms) { + sum = Add(sum, x); + scale = Add(scale, Abs(x)); + } + const auto defect = Abs(Subtract(left, sum)); + if (!Valid(defect) || !Valid(scale)) return math_internal::InvalidHighLow(); + if (Compare(scale, HighLow(0)) == 0) + return Compare(defect, HighLow(0)) == 0 ? HighLow(0) + : math_internal::InvalidHighLow(); + return Divide(defect, scale); +} +HighLow CheckUpdates(const SolidDynamicAcceptedFrameState& old, + const SolidDynamicAcceptedFrameState& f, + const HhtParameters& p, double h) { + HighLow metric; + const auto step = HighLow(h); + const auto hh = Multiply(step, step); + for (std::size_t i = 0; i < f.displacement.Size(); ++i) { + const auto A = + Pair(f.algorithmic_acceleration, f.low.algorithmic_acceleration, i); + const auto oldA = + Pair(old.algorithmic_acceleration, old.low.algorithmic_acceleration, i); + const auto oldv = Pair(old.velocity, old.low.velocity, i); + metric = Maximum( + metric, + UpdateMetric( + Pair(f.displacement, f.low.displacement, i), + {Pair(old.displacement, old.low.displacement, i), + Multiply(step, oldv), + Multiply(Multiply(hh, Subtract(HighLow(.5), HighLow(p.beta))), + oldA), + Multiply(Multiply(hh, HighLow(p.beta)), A)})); + metric = Maximum( + metric, + UpdateMetric( + Pair(f.velocity, f.low.velocity, i), + {oldv, + Multiply(Multiply(step, Subtract(HighLow(1), HighLow(p.gamma))), + oldA), + Multiply(Multiply(step, HighLow(p.gamma)), A)})); + metric = Maximum( + metric, + UpdateMetric(A, {Multiply(Add(HighLow(1), HighLow(p.alpha)), + Pair(f.acceleration, f.low.acceleration, i)), + Negate(Multiply(HighLow(p.alpha), + Pair(old.acceleration, + old.low.acceleration, i)))})); + } + return metric; +} +std::array FreshCorrectionMetrics(const Vector& du, double length, + const HhtParameters& p) { + HighLow q; + for (std::size_t i = 0; i < du.Size(); ++i) + q = Maximum(q, Divide(Abs(HighLow(du[i])), HighLow(length))); + return { + q, Divide(Multiply(HighLow(p.gamma), q), HighLow(p.beta)), + Divide(q, Multiply(Add(HighLow(1), HighLow(p.alpha)), HighLow(p.beta))), + Divide(q, HighLow(p.beta))}; +} +std::array ObservedCorrectionMetrics( + const SolidDynamicAcceptedFrameState& before, + const SolidDynamicAcceptedFrameState& after, double h, double length) { + std::array metric{}; + const auto hh = Multiply(HighLow(h), HighLow(h)); + for (std::size_t i = 0; i < after.displacement.Size(); ++i) { + metric[0] = Maximum( + metric[0], + Divide(Abs(Subtract( + Pair(after.displacement, after.low.displacement, i), + Pair(before.displacement, before.low.displacement, i))), + HighLow(length))); + metric[1] = Maximum( + metric[1], + Divide(Multiply(HighLow(h), + Abs(Subtract( + Pair(after.velocity, after.low.velocity, i), + Pair(before.velocity, before.low.velocity, i)))), + HighLow(length))); + metric[2] = Maximum( + metric[2], + Divide(Multiply(hh, Abs(Subtract(Pair(after.acceleration, + after.low.acceleration, i), + Pair(before.acceleration, + before.low.acceleration, i)))), + HighLow(length))); + metric[3] = Maximum( + metric[3], + Divide(Multiply(hh, Abs(Subtract( + Pair(after.algorithmic_acceleration, + after.low.algorithmic_acceleration, i), + Pair(before.algorithmic_acceleration, + before.low.algorithmic_acceleration, i)))), + HighLow(length))); + } + return metric; +} +bool CorrectionPass(const std::array& metrics) { + return std::all_of(metrics.begin(), metrics.end(), [](HighLow x) { + return Valid(x) && Compare(x, HighLow(0)) >= 0 && + Compare(x, HighLow(1e-10)) <= 0; + }); +} +Status SetRates(const SolidDynamicAcceptedFrameState& old, + const detail::SolidHhtPairedPredictor& predictor, double h, + const HhtParameters& p, + SolidDynamicAcceptedFrameState* candidate) { + auto rates = detail::EvaluateSolidHhtHighLowRates( + candidate->displacement, candidate->low.displacement, predictor, + old.acceleration, old.low.acceleration, h, p); + if (!rates.HasValue()) return rates.GetStatus(); + candidate->velocity = std::move(rates.Value().velocity); + candidate->low.velocity = std::move(rates.Value().velocity_low); + candidate->acceleration = std::move(rates.Value().acceleration); + candidate->low.acceleration = std::move(rates.Value().acceleration_low); + candidate->algorithmic_acceleration = + std::move(rates.Value().algorithmic_acceleration); + candidate->low.algorithmic_acceleration = + std::move(rates.Value().algorithmic_acceleration_low); + return Status::Ok(); +} +} // namespace + +SolidHhtDynamicProcedureExecutor::SolidHhtDynamicProcedureExecutor( + const ParallelFor& parallel_for, LinearSolver& mass_solver, + TangentSolver& tangent_solver) + : parallel_for_(parallel_for), + mass_solver_(mass_solver), + tangent_solver_(tangent_solver) {} + +Result SolidHhtDynamicProcedureExecutor::Execute( + const Domain& domain) const { + auto model = AnalysisModel::Create(domain); + if (!model.HasValue()) + return Result::Failure(model.GetStatus()); + auto dofs = DofManager::Create(model.Value()); + if (!dofs.HasValue()) + return Result::Failure(dofs.GetStatus()); + auto state = SolidDynamicAnalysisState::Create(domain, dofs.Value()); + if (!state.HasValue()) + return Result::Failure(state.GetStatus()); + auto status = ExecuteInto(domain, state.Value()); + if (!status.IsOk()) + return Result::Failure(status); + return Result::Success( + {std::move(dofs.Value()), std::move(state.Value())}); +} +Status SolidHhtDynamicProcedureExecutor::ExecuteInto( + const Domain& domain, SolidDynamicAnalysisState& state) const { + return ExecuteIntoForVerification(domain, state, -0.05); +} +Status SolidHhtDynamicProcedureExecutor::ExecuteIntoForVerification( + const Domain& domain, SolidDynamicAnalysisState& state, + double alpha) const { + SolidDynamicIterationRecord record; + auto fail = [&](const Status& failure) { + record.stage = "failed"; + record.outcome = "failed"; + record.failure_code = failure.Diagnostics().empty() + ? "solid-dynamic-failure" + : failure.Diagnostics().front().code; + state.RecordIteration(record); + return failure; + }; + if (!state.IsInitialized() || !state.AcceptedFrames().empty() || + !domain.IsSolid() || domain.Dynamic() == nullptr) + return fail(Failure( + domain, "solid-dynamic-invalid-execution", + "Execution requires an initialized empty Solid dynamic state.")); + auto environment = CheckArithmeticEnvironment(domain); + if (!environment.IsOk()) return fail(environment); + auto model = AnalysisModel::Create(domain); + if (!model.HasValue()) return fail(model.GetStatus()); + auto dofs_result = DofManager::Create(model.Value()); + if (!dofs_result.HasValue()) return fail(dofs_result.GetStatus()); + const auto& dofs = dofs_result.Value(); + auto solids = BuildSolids(domain, dofs); + if (!solids.HasValue()) return fail(solids.GetStatus()); + auto mass = AssembleMass(solids.Value(), dofs); + if (!mass.HasValue()) return fail(mass.GetStatus()); + EssentialConstraintPolicy partition; + auto mass_partition = partition.Partition(mass.Value(), dofs); + if (!mass_partition.HasValue()) return fail(mass_partition.GetStatus()); + if (dofs.FreeDofCount() != 0) { + const auto status = mass_solver_.Factorize(mass_partition.Value().k_ff); + if (!status.IsOk()) return fail(status); + environment = CheckArithmeticEnvironment(domain); + if (!environment.IsOk()) return fail(environment); + } + auto parameters = MakeHhtParameters(alpha); + if (!parameters.HasValue()) return fail(parameters.GetStatus()); + const auto p = parameters.Value(); + auto length = ReferenceLength(domain); + if (!length.HasValue()) return fail(length.GetStatus()); + const double h = domain.Dynamic()->initial_increment; + SolidDynamicAcceptedFrameState initial; + initial.is_initial = true; + for (std::size_t n = 0; n < domain.Nodes().size(); ++n) { + initial.nodal_internal_order.push_back(static_cast(n)); + initial.nodal_source_ids.push_back(domain.Nodes()[n].source_id); + } + for (auto* v : {&initial.displacement, &initial.velocity, + &initial.acceleration, &initial.algorithmic_acceleration}) + *v = Vector(dofs.FullDofCount()); + InitializeLowVectors(&initial, dofs.FullDofCount()); + auto force = LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs, 0); + if (!force.HasValue()) return fail(force.GetStatus()); + initial.external_force = std::move(force.Value()); + auto tangent = + EvaluateSolids(domain, solids.Value(), dofs, parallel_for_, &initial); + if (!tangent.HasValue()) return fail(tangent.GetStatus()); + SolidDynamicIterationEvidence evidence; + bool initial_converged = false; + for (std::uint32_t refinement = 0; refinement <= 100; ++refinement) { + initial.algorithmic_acceleration = initial.acceleration; + initial.low.algorithmic_acceleration = initial.low.acceleration; + evidence = {}; + evidence.correction_pass = evidence.update_pass = true; + auto status = RecoverEndpoint(domain, dofs, mass.Value(), nullptr, + HighLow(0), &initial, &evidence); + if (!status.IsOk()) return fail(status); + if (evidence.residual_pass) { + initial_converged = true; + break; + } + record.stage = "initial-mass"; + record.outcome = "residual-refinement"; + record.evidence = evidence; + record.evidence.iteration_index = refinement; + state.RecordIteration(record); + if (refinement == 100) break; + Vector rhs(dofs.FreeDofCount()), a(dofs.FreeDofCount()); + for (std::size_t j = 0; j < dofs.FreeDofCount(); ++j) { + const auto i = dofs.FreeDofs()[j]; + rhs[j] = Rounded(Negate(Pair(initial.residual, initial.low.residual, i))); + } + status = mass_solver_.Solve(rhs, a); + if (!status.IsOk()) return fail(status); + environment = CheckArithmeticEnvironment(domain); + if (!environment.IsOk()) return fail(environment); + if (a.Size() != dofs.FreeDofCount()) + return fail( + Failure(domain, "solid-dynamic-invalid-backend-result", + "Initial mass solve returned incompatible acceleration.")); + for (std::size_t j = 0; j < a.Size(); ++j) { + if (!std::isfinite(a[j])) + return fail( + Failure(domain, "solid-dynamic-invalid-backend-result", + "Initial mass solve returned nonfinite acceleration.")); + const auto i = dofs.FreeDofs()[j]; + const auto value = + Add(Pair(initial.acceleration, initial.low.acceleration, i), + HighLow(a[j] == 0 ? 0.0 : a[j])); + if (!Valid(value)) return fail(detail::SolidHhtPairedKinematicsFailure()); + Store(value, &initial.acceleration, &initial.low.acceleration, i); + } + } + if (!initial_converged) + return fail(Failure(domain, "solid-dynamic-initial-mass-nonconvergence", + "Physical initial equilibrium did not converge in100 " + "mass refinements.")); + initial.iterations.push_back(evidence); + auto status = state.CommitAcceptedFrame(initial); + if (!status.IsOk()) { + record.evidence = evidence; + return fail(status); + } + record.stage = "accepted"; + record.outcome = "initial"; + record.evidence = evidence; + state.RecordIteration(record); + const auto initial_energy = + Add(HighLow(initial.strain_energy, initial.low.strain_energy), + HighLow(initial.kinetic_energy, initial.low.kinetic_energy)); + for (std::uint32_t increment = 1; + increment <= domain.Dynamic()->increment_count; ++increment) { + const auto& old = state.CommittedState(); + record = {}; + record.increment_index = increment; + record.time = static_cast(increment) * h; + auto predictor = detail::PredictSolidHhtHighLow(old, h, p); + if (!predictor.HasValue()) return fail(predictor.GetStatus()); + auto candidate = old; + candidate.frame_index = candidate.increment_index = increment; + candidate.iteration_count = 0; + candidate.time = record.time; + candidate.is_initial = false; + candidate.iterations.clear(); + // Begin from certified accepted geometry, not an unchecked predictor. + force = LoadAssembler::AssembleFullNodalLoad(model.Value(), dofs, + candidate.time); + if (!force.HasValue()) return fail(force.GetStatus()); + candidate.external_force = std::move(force.Value()); + candidate.low.external_force = Vector(dofs.FullDofCount()); + std::array correction{}, fresh_correction{}; + bool accepted = false; + for (std::uint32_t iteration = 0; iteration <= 100; ++iteration) { + record.stage = "evaluation"; + record.outcome = "paired-candidate"; + record.evidence = {}; + record.evidence.iteration_index = iteration; + status = SetRates(old, predictor.Value(), h, p, &candidate); + if (!status.IsOk()) return fail(status); + tangent = EvaluateSolids(domain, solids.Value(), dofs, parallel_for_, + &candidate); + if (!tangent.HasValue()) return fail(tangent.GetStatus()); + status = RecoverEndpoint(domain, dofs, mass.Value(), &old, initial_energy, + &candidate, &record.evidence); + if (!status.IsOk()) return fail(status); + auto MA = detail::MultiplySolidHhtMassHighLow( + mass.Value(), candidate.algorithmic_acceleration, + candidate.low.algorithmic_acceleration); + if (!MA.HasValue()) return fail(MA.GetStatus()); + Vector rhs(dofs.FreeDofCount()); + HighLow identity_metric, original_norm, solved_norm, old_norm; + for (std::size_t j = 0; j < dofs.FreeDofCount(); ++j) { + const auto i = dofs.FreeDofs()[j]; + detail::SolidHhtPairedForceIdentity identity; + if (!detail::EvaluateSolidHhtHighLowForceIdentity( + mass.Value(), i, Pair(MA.Value().high, MA.Value().low, i), + candidate, old, p.alpha, &identity)) + return fail( + Failure(domain, "solid-dynamic-nonfinite-residual", + "Paired weighted force arithmetic is unsupported.")); + rhs[j] = Rounded(Negate(identity.solved)); + identity_metric = Maximum(identity_metric, identity.metric); + original_norm = Maximum(original_norm, Abs(identity.original)); + solved_norm = Maximum(solved_norm, Abs(identity.solved)); + old_norm = + Maximum(old_norm, Abs(Pair(old.residual, old.low.residual, i))); + if (i == record.evidence.worst_full_dof) + record.evidence.worst_weighted_terms = { + Rounded(identity.original), Rounded(identity.solved), + Rounded(Pair(old.residual, old.low.residual, i)), + Rounded(identity.scale)}; + } + if (!Valid(identity_metric) || !Valid(original_norm) || + !Valid(solved_norm) || !Valid(old_norm) || + Compare(identity_metric, HighLow(1e-12)) > 0) + return fail(Failure(domain, + "solid-dynamic-weighted-algebra-inconsistency", + "Paired original and solved forces disagree beyond " + "their uncancelled scale.")); + record.evidence.weighted_equivalence_metric = Rounded(identity_metric); + record.evidence.weighted_residual_numerator = Rounded(original_norm); + record.evidence.recombined_residual_numerator = Rounded(solved_norm); + record.evidence.old_residual_numerator = Rounded(old_norm); + for (std::size_t c = 0; c < 4; ++c) { + if (!Valid(correction[c]) || !Valid(fresh_correction[c])) + return fail(detail::SolidHhtPairedKinematicsFailure()); + record.evidence.correction_metrics[c] = Rounded(correction[c]); + } + record.evidence.correction_pass = + CorrectionPass(correction) && CorrectionPass(fresh_correction); + const auto update = CheckUpdates(old, candidate, p, h); + if (!Valid(update)) + return fail(detail::SolidHhtPairedKinematicsFailure()); + record.evidence.update_metric = Rounded(update); + record.evidence.update_pass = Compare(update, HighLow(1e-12)) <= 0; + candidate.iterations.push_back(record.evidence); + candidate.iteration_count = iteration; + if (record.evidence.residual_pass && record.evidence.correction_pass && + record.evidence.update_pass) { + status = state.CommitAcceptedFrame(candidate); + if (!status.IsOk()) return fail(status); + record.stage = "accepted"; + record.outcome = "increment"; + state.RecordIteration(record); + accepted = true; + break; + } + state.RecordIteration(record); + if (iteration == 100) + return fail(Failure( + domain, "solid-dynamic-nonconvergence", + "Paired HHT endpoint did not converge in100 Newton solves.")); + auto effective = + EffectiveTangent(mass.Value(), tangent.Value(), dofs, h, p); + if (!effective.HasValue()) return fail(effective.GetStatus()); + auto reduced = partition.Partition(effective.Value(), dofs); + if (!reduced.HasValue()) return fail(reduced.GetStatus()); + Vector du(dofs.FreeDofCount()); + if (dofs.FreeDofCount() != 0) { + status = tangent_solver_.FactorizeFresh(reduced.Value().k_ff); + if (!status.IsOk()) return fail(status); + environment = CheckArithmeticEnvironment(domain); + if (!environment.IsOk()) return fail(environment); + status = tangent_solver_.SolveCorrection(rhs, du); + if (!status.IsOk()) return fail(status); + environment = CheckArithmeticEnvironment(domain); + if (!environment.IsOk()) return fail(environment); + } + if (du.Size() != dofs.FreeDofCount()) + return fail(Failure(domain, "solid-dynamic-invalid-backend-result", + "Tangent solve returned incompatible correction.")); + for (std::size_t j = 0; j < du.Size(); ++j) + if (!std::isfinite(du[j])) + return fail(Failure(domain, "solid-dynamic-invalid-backend-result", + "Tangent solve returned nonfinite correction.")); + fresh_correction = FreshCorrectionMetrics(du, length.Value(), p); + const auto before = candidate; + for (std::size_t j = 0; j < du.Size(); ++j) { + const auto i = dofs.FreeDofs()[j]; + const auto value = + Add(Pair(candidate.displacement, candidate.low.displacement, i), + HighLow(du[j] == 0 ? 0.0 : du[j])); + if (!Valid(value)) + return fail(detail::SolidHhtPairedKinematicsFailure()); + Store(value, &candidate.displacement, &candidate.low.displacement, i); + } + status = SetRates(old, predictor.Value(), h, p, &candidate); + if (!status.IsOk()) return fail(status); + correction = + ObservedCorrectionMetrics(before, candidate, h, length.Value()); + } + if (!accepted) + return fail(Failure(domain, "solid-dynamic-nonconvergence", + "No complete paired endpoint was accepted.")); + } + return Status::Ok(); +} +} // namespace fesa diff --git a/src/fesa/analysis/solid_hht_force_identity.h b/src/fesa/analysis/solid_hht_force_identity.h new file mode 100644 index 0000000..c091fe1 --- /dev/null +++ b/src/fesa/analysis/solid_hht_force_identity.h @@ -0,0 +1,182 @@ +#ifndef FESA_ANALYSIS_SOLID_HHT_FORCE_IDENTITY_H_ +#define FESA_ANALYSIS_SOLID_HHT_FORCE_IDENTITY_H_ + +#include +#include +#include + +#include "fesa/analysis/solid_dynamic_analysis_state.h" +#include "fesa/math/sparse_matrix.h" +#include "solid_hht_high_low.h" + +namespace fesa::detail { +struct SolidHhtForceIdentity { + double original{0}; + double solved{0}; + double scale{0}; + double defect{0}; + double metric{0}; +}; + +struct SolidHhtPairedForceIdentity { + math_internal::HighLow original; + math_internal::HighLow solved; + math_internal::HighLow scale; + math_internal::HighLow defect; + math_internal::HighLow metric; +}; + +/// @brief Keeps both physical parts and the actual old residual in the approved +/// weighted identity. Its scale includes each uncancelled mass product. +inline bool EvaluateSolidHhtHighLowForceIdentity( + const SparseMatrix& mass, std::size_t row, + math_internal::HighLow algorithmic_inertia, + const SolidDynamicAcceptedFrameState& next, + const SolidDynamicAcceptedFrameState& old, double alpha, + SolidHhtPairedForceIdentity* result) { + using namespace math_internal; + if (!ArithmeticEnvironmentSupported() || row >= mass.Rows() || + !IsFinite(algorithmic_inertia) || !IsNormalized(algorithmic_inertia) || + !IsFinite(HighLow(alpha))) + return false; + for (const auto* frame : {&next, &old}) + if (frame->acceleration.Size() != mass.Columns() || + frame->low.acceleration.Size() != mass.Columns() || + frame->algorithmic_acceleration.Size() != mass.Columns() || + frame->low.algorithmic_acceleration.Size() != mass.Columns() || + frame->internal_force.Size() != mass.Rows() || + frame->low.internal_force.Size() != mass.Rows() || + frame->external_force.Size() != mass.Rows() || + frame->low.external_force.Size() != mass.Rows() || + frame->residual.Size() != mass.Rows() || + frame->low.residual.Size() != mass.Rows()) + return false; + const HighLow a(alpha); + const auto c = Add(HighLow(1), a); + const HighLow f(next.internal_force[row], next.low.internal_force[row]); + const HighLow F(next.external_force[row], next.low.external_force[row]); + const HighLow old_f(old.internal_force[row], old.low.internal_force[row]); + const HighLow old_F(old.external_force[row], old.low.external_force[row]); + const HighLow r(next.residual[row], next.low.residual[row]); + const HighLow old_r(old.residual[row], old.low.residual[row]); + result->original = + Subtract(Add(algorithmic_inertia, Multiply(c, Subtract(f, F))), + Multiply(a, Subtract(old_f, old_F))); + result->solved = Subtract(Multiply(c, r), Multiply(a, old_r)); + HighLow mass_scale; + HighLow relation_action; + for (std::size_t k = mass.RowOffsets()[row]; k < mass.RowOffsets()[row + 1]; + ++k) { + const auto j = mass.ColumnIndices()[k]; + const HighLow coefficient(mass.Values()[k]); + const HighLow A(next.algorithmic_acceleration[j], + next.low.algorithmic_acceleration[j]); + const HighLow physical_a(next.acceleration[j], next.low.acceleration[j]); + const HighLow old_a(old.acceleration[j], old.low.acceleration[j]); + const auto relation = + Add(Subtract(A, Multiply(c, physical_a)), Multiply(a, old_a)); + mass_scale = + Add(mass_scale, + Add(Add(Abs(Multiply(coefficient, A)), + Multiply(c, Abs(Multiply(coefficient, physical_a)))), + Multiply(Abs(a), Abs(Multiply(coefficient, old_a))))); + relation_action = + Add(relation_action, Abs(Multiply(coefficient, relation))); + if (!IsFinite(mass_scale) || !IsFinite(relation_action)) return false; + } + result->scale = Add( + Add(Add(mass_scale, relation_action), Multiply(c, Add(Abs(f), Abs(F)))), + Add(Add(Multiply(Abs(a), Add(Abs(old_f), Abs(old_F))), + Multiply(c, Abs(r))), + Abs(Multiply(a, old_r)))); + result->defect = Abs(Subtract(result->original, result->solved)); + if (!IsFinite(result->original) || !IsFinite(result->solved) || + !IsFinite(result->scale) || !IsFinite(result->defect)) + return false; + if (Compare(result->scale, HighLow(0)) == 0) { + if (Compare(result->defect, HighLow(0)) != 0) return false; + result->metric = HighLow(0); + } else { + result->metric = Divide(result->defect, result->scale); + } + return IsFinite(result->metric) && IsNormalized(result->metric); +} + +/// @brief Evaluates one original/equivalent HHT row, retaining old equilibrium +/// error. Checks each scalar before it can enter an evidence maximum. +inline bool EvaluateSolidHhtForceIdentity( + const SparseMatrix& mass, std::size_t row, double algorithmic_inertia, + const SolidDynamicAcceptedFrameState& next, + const SolidDynamicAcceptedFrameState& old, double alpha, + SolidHhtForceIdentity* result) { + for (double x : + {alpha, algorithmic_inertia, next.internal_force[row], + next.external_force[row], old.internal_force[row], + old.external_force[row], next.residual[row], old.residual[row]}) + if (!std::isfinite(x)) return false; + result->original = + algorithmic_inertia + + (1 + alpha) * (next.internal_force[row] - next.external_force[row]) - + alpha * (old.internal_force[row] - old.external_force[row]); + result->solved = (1 + alpha) * next.residual[row] - alpha * old.residual[row]; + double mass_action_scale = 0; + double relation_defect_action = 0; + for (std::size_t entry = mass.RowOffsets()[row]; + entry < mass.RowOffsets()[row + 1]; ++entry) { + const auto column = mass.ColumnIndices()[entry]; + const double coefficient = mass.Values()[entry]; + const double a = next.acceleration[column]; + const double A = next.algorithmic_acceleration[column]; + const double a_old = old.acceleration[column]; + const double relation_defect = A - (1 + alpha) * a + alpha * a_old; + const double action_scale = std::abs(coefficient * A) + + (1 + alpha) * std::abs(coefficient * a) + + std::abs(alpha) * std::abs(coefficient * a_old); + const double defect_action = std::abs(coefficient * relation_defect); + for (double x : {a, A, a_old, relation_defect, action_scale, defect_action}) + if (!std::isfinite(x)) return false; + mass_action_scale += action_scale; + relation_defect_action += defect_action; + } + result->scale = mass_action_scale + relation_defect_action + + (1 + alpha) * (std::abs(next.internal_force[row]) + + std::abs(next.external_force[row])) + + std::abs(alpha) * (std::abs(old.internal_force[row]) + + std::abs(old.external_force[row])) + + (1 + alpha) * std::abs(next.residual[row]) + + std::abs(alpha * old.residual[row]); + result->defect = std::abs(result->original - result->solved); + for (double x : {result->original, result->solved, mass_action_scale, + relation_defect_action, result->scale, result->defect}) + if (!std::isfinite(x)) return false; + result->metric = + result->scale == 0 + ? (result->defect == 0 ? 0 : std::numeric_limits::infinity()) + : result->defect / result->scale; + return std::isfinite(result->metric); +} +/// @brief Finite checked failure-shadow accumulation, separate from physics. +inline bool AccumulateSolidHhtForceShadow(double term, double product_error, + double* sum, double* correction, + double* absolute_sum) { + if (!std::isfinite(term) || !std::isfinite(product_error) || + !std::isfinite(*sum) || !std::isfinite(*correction) || + !std::isfinite(*absolute_sum)) + return false; + const double next = *sum + term; + const double absolute_next = *absolute_sum + std::abs(term); + const double error_next = *correction + product_error; + const double lost = std::abs(*sum) >= std::abs(term) ? (*sum - next) + term + : (term - next) + *sum; + const double correction_next = error_next + lost; + if (!std::isfinite(next) || !std::isfinite(absolute_next) || + !std::isfinite(error_next) || !std::isfinite(lost) || + !std::isfinite(correction_next)) + return false; + *sum = next; + *correction = correction_next; + *absolute_sum = absolute_next; + return true; +} +} // namespace fesa::detail +#endif // FESA_ANALYSIS_SOLID_HHT_FORCE_IDENTITY_H_ diff --git a/src/fesa/analysis/solid_hht_high_low.h b/src/fesa/analysis/solid_hht_high_low.h new file mode 100644 index 0000000..45a0794 --- /dev/null +++ b/src/fesa/analysis/solid_hht_high_low.h @@ -0,0 +1,190 @@ +#ifndef FESA_ANALYSIS_SOLID_HHT_HIGH_LOW_H_ +#define FESA_ANALYSIS_SOLID_HHT_HIGH_LOW_H_ + +#include "fesa/analysis/hht_time_integration.h" +#include "fesa/analysis/solid_dynamic_analysis_state.h" +#include "fesa/math/sparse_matrix.h" +#include "math/high_low_arithmetic.h" + +namespace fesa::detail { +/// @brief Keeps both physical parts of the fixed-increment predictor. +struct SolidHhtPairedPredictor { + Vector displacement{0}; + Vector displacement_low{0}; + Vector velocity{0}; + Vector velocity_low{0}; +}; + +/// @brief Owns physical and auxiliary rates without collapsing either part. +struct SolidHhtPairedRates { + Vector velocity{0}; + Vector velocity_low{0}; + Vector acceleration{0}; + Vector acceleration_low{0}; + Vector algorithmic_acceleration{0}; + Vector algorithmic_acceleration_low{0}; +}; + +struct SolidHhtPairedVector { + Vector high{0}; + Vector low{0}; +}; + +inline Status SolidHhtPairedKinematicsFailure() { + return Status::Failure(FailureCategory::kSolver, + {{Severity::kError, + "solid-dynamic-invalid-paired-kinematics", + {}, + "DYNAMIC", + {}, + "Incomplete, unnormalized or nonfinite paired HHT " + "kinematics."}}); +} + +inline bool ValidSolidHhtPair(const Vector& high, const Vector& low, + std::size_t size) { + if (high.Size() != size || low.Size() != size) return false; + for (std::size_t i = 0; i < size; ++i) { + const math_internal::HighLow value(high[i], low[i]); + if (!math_internal::IsFinite(value) || !math_internal::IsNormalized(value)) + return false; + } + return true; +} + +inline bool ValidSolidHhtParameters(double h, const HhtParameters& p) { + if (!math_internal::ArithmeticEnvironmentSupported()) return false; + const auto expected = MakeHhtParameters(p.alpha); + if (!expected.HasValue() || p.beta != expected.Value().beta || + p.gamma != expected.Value().gamma) + return false; + const math_internal::HighLow step(h); + const auto square = math_internal::Multiply(step, step); + return math_internal::IsFinite(step) && + math_internal::Compare(step, math_internal::HighLow(0)) > 0 && + math_internal::IsFinite(square) && + math_internal::Compare(square, math_internal::HighLow(0)) > 0; +} + +/// @brief Uses the authoritative constant Float64 mass coefficients without +/// collapsing either physical vector part or the ordered row accumulation. +inline Result MultiplySolidHhtMassHighLow( + const SparseMatrix& mass, const Vector& high, const Vector& low) { + using namespace math_internal; + if (!ArithmeticEnvironmentSupported() || mass.Rows() != mass.Columns() || + !ValidSolidHhtPair(high, low, mass.Columns())) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + SolidHhtPairedVector result{Vector(mass.Rows()), Vector(mass.Rows())}; + for (std::size_t row = 0; row < mass.Rows(); ++row) { + HighLow sum; + for (std::size_t k = mass.RowOffsets()[row]; k < mass.RowOffsets()[row + 1]; + ++k) { + const auto col = mass.ColumnIndices()[k]; + sum = Add(sum, Multiply(HighLow(mass.Values()[k]), + HighLow(high[col], low[col]))); + if (!IsFinite(sum) || !IsNormalized(sum)) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + } + result.high[row] = sum.high; + result.low[row] = sum.low; + } + return Result::Success(std::move(result)); +} + +/// @brief Applies the unchanged Newmark predictor to the full physical state. +inline Result PredictSolidHhtHighLow( + const SolidDynamicAcceptedFrameState& old, double h, + const HhtParameters& p) { + using namespace math_internal; + const auto size = old.displacement.Size(); + if (!ValidSolidHhtPair(old.displacement, old.low.displacement, size) || + !ValidSolidHhtPair(old.velocity, old.low.velocity, size) || + !ValidSolidHhtPair(old.algorithmic_acceleration, + old.low.algorithmic_acceleration, size) || + !ValidSolidHhtParameters(h, p)) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + SolidHhtPairedPredictor result{Vector(size), Vector(size), Vector(size), + Vector(size)}; + const auto hh = Multiply(HighLow(h), HighLow(h)); + const auto u_coefficient = + Multiply(hh, Subtract(HighLow(.5), HighLow(p.beta))); + const auto v_coefficient = + Multiply(HighLow(h), Subtract(HighLow(1), HighLow(p.gamma))); + if (!IsFinite(u_coefficient) || !IsFinite(v_coefficient)) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + for (std::size_t i = 0; i < size; ++i) { + const HighLow u(old.displacement[i], old.low.displacement[i]); + const HighLow v(old.velocity[i], old.low.velocity[i]); + const HighLow A(old.algorithmic_acceleration[i], + old.low.algorithmic_acceleration[i]); + const auto up = + Add(Add(u, Multiply(HighLow(h), v)), Multiply(u_coefficient, A)); + const auto vp = Add(v, Multiply(v_coefficient, A)); + if (!IsFinite(up) || !IsNormalized(up) || !IsFinite(vp) || + !IsNormalized(vp)) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + result.displacement[i] = up.high; + result.displacement_low[i] = up.low; + result.velocity[i] = vp.high; + result.velocity_low[i] = vp.low; + } + return Result::Success(std::move(result)); +} + +/// @brief Recovers A, v and physical a from the paired displacement difference. +inline Result EvaluateSolidHhtHighLowRates( + const Vector& high, const Vector& low, + const SolidHhtPairedPredictor& predictor, const Vector& old_acceleration, + const Vector& old_acceleration_low, double h, const HhtParameters& p) { + using namespace math_internal; + const auto size = high.Size(); + if (!ValidSolidHhtPair(high, low, size) || + !ValidSolidHhtPair(predictor.displacement, predictor.displacement_low, + size) || + !ValidSolidHhtPair(predictor.velocity, predictor.velocity_low, size) || + !ValidSolidHhtPair(old_acceleration, old_acceleration_low, size) || + !ValidSolidHhtParameters(h, p)) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + SolidHhtPairedRates result{Vector(size), Vector(size), Vector(size), + Vector(size), Vector(size), Vector(size)}; + const auto denominator = + Multiply(Multiply(HighLow(p.beta), HighLow(h)), HighLow(h)); + const auto velocity_coefficient = Multiply(HighLow(p.gamma), HighLow(h)); + const auto physical_coefficient = Add(HighLow(1), HighLow(p.alpha)); + if (!IsFinite(denominator) || Compare(denominator, HighLow(0)) <= 0 || + !IsFinite(velocity_coefficient) || !IsFinite(physical_coefficient)) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + for (std::size_t i = 0; i < size; ++i) { + const auto A = Divide(Subtract(HighLow(high[i], low[i]), + HighLow(predictor.displacement[i], + predictor.displacement_low[i])), + denominator); + const auto v = + Add(HighLow(predictor.velocity[i], predictor.velocity_low[i]), + Multiply(velocity_coefficient, A)); + const auto a = Divide( + Add(A, Multiply(HighLow(p.alpha), + HighLow(old_acceleration[i], old_acceleration_low[i]))), + physical_coefficient); + if (!IsFinite(A) || !IsNormalized(A) || !IsFinite(v) || !IsNormalized(v) || + !IsFinite(a) || !IsNormalized(a)) + return Result::Failure( + SolidHhtPairedKinematicsFailure()); + result.algorithmic_acceleration[i] = A.high; + result.algorithmic_acceleration_low[i] = A.low; + result.velocity[i] = v.high; + result.velocity_low[i] = v.low; + result.acceleration[i] = a.high; + result.acceleration_low[i] = a.low; + } + return Result::Success(std::move(result)); +} +} // namespace fesa::detail +#endif // FESA_ANALYSIS_SOLID_HHT_HIGH_LOW_H_ diff --git a/src/fesa/analysis/solid_hht_representability.h b/src/fesa/analysis/solid_hht_representability.h new file mode 100644 index 0000000..207db94 --- /dev/null +++ b/src/fesa/analysis/solid_hht_representability.h @@ -0,0 +1,217 @@ +#ifndef FESA_ANALYSIS_SOLID_HHT_REPRESENTABILITY_H_ +#define FESA_ANALYSIS_SOLID_HHT_REPRESENTABILITY_H_ +#include +#include +#include +#include +#include + +#include "fesa/analysis/solid_dynamic_analysis_state.h" +#include "fesa/math/sparse_matrix.h" +namespace fesa::detail { +struct SolidHhtRepresentableNeighbor { + bool found{false}; + std::size_t full_dof{0}; + double displacement{0}; + double delta{0}; + double predicted_metric{0}; +}; +/// @brief Predicts a strictly improving one-ULP neighbor using separate +/// physical inertia/internal-force actions and every original free score. +/// This selects a trial only; callers must recompute actual endpoint gates. +inline bool SelectSolidHhtRepresentableNeighbor( + const SparseMatrix& mass, const SparseMatrix& tangent, + const std::vector& free_dofs, + const SolidDynamicAcceptedFrameState& frame, double h, double alpha, + double beta, SolidHhtRepresentableNeighbor* selected, + std::vector* ranked = nullptr) { + *selected = {}; + if (ranked != nullptr) ranked->clear(); + const auto size = mass.Rows(); + if (tangent.Rows() != size || mass.Columns() != size || + tangent.Columns() != size || mass.RowOffsets() != tangent.RowOffsets() || + mass.ColumnIndices() != tangent.ColumnIndices() || + frame.displacement.Size() != size || frame.inertia_force.Size() != size || + frame.internal_force.Size() != size || + frame.external_force.Size() != size || frame.residual.Size() != size || + frame.residual_scale.Size() != size) + return false; + const double inertia_factor = 1 / ((1 + alpha) * beta * h * h); + if (!std::isfinite(inertia_factor) || inertia_factor <= 0) return false; + const auto score = [](double ma, double f, double F, double* metric) { + const double scale = std::abs(ma) + std::abs(f) + std::abs(F); + const double r = ma + f - F; + if (!std::isfinite(scale) || !std::isfinite(r)) return false; + *metric = scale == 0 + ? (r == 0 ? 0 : std::numeric_limits::infinity()) + : std::abs(r) / scale; + return std::isfinite(*metric); + }; + std::vector is_free(size, false); + std::vector scores(size, 0); + double best = 0; + for (auto i : free_dofs) { + if (i >= size || !std::isfinite(frame.displacement[i]) || + !score(frame.inertia_force[i], frame.internal_force[i], + frame.external_force[i], &scores[i])) + return false; + is_free[i] = true; + best = (std::max)(best, scores[i]); + } + struct Action { + std::size_t row; + double inertia; + double internal; + }; + std::vector> columns(size); + for (auto i : free_dofs) + for (std::size_t k = mass.RowOffsets()[i]; k < mass.RowOffsets()[i + 1]; + ++k) { + const auto j = mass.ColumnIndices()[k]; + if (!is_free[j]) continue; + const double inertia = mass.Values()[k] * inertia_factor; + const double internal = tangent.Values()[k]; + if (!std::isfinite(inertia) || !std::isfinite(internal)) return false; + columns[j].push_back({i, inertia, internal}); + } + const double original_metric = best; + auto predicted = scores; + for (auto j : free_dofs) + for (const double direction : {-std::numeric_limits::infinity(), + std::numeric_limits::infinity()}) { + const double proposed = std::nextafter(frame.displacement[j], direction); + const double delta = proposed - frame.displacement[j]; + if (!std::isfinite(proposed) || !std::isfinite(delta) || delta == 0) + continue; + predicted = scores; + for (const auto& action : columns[j]) { + const auto i = action.row; + const double ma = frame.inertia_force[i] + action.inertia * delta; + const double f = frame.internal_force[i] + action.internal * delta; + if (!score(ma, f, frame.external_force[i], &predicted[i])) return false; + } + double metric = 0; + for (auto i : free_dofs) metric = (std::max)(metric, predicted[i]); + if (ranked != nullptr && metric < original_metric) + ranked->push_back({true, j, proposed, delta, metric}); + // Stable full-DOF/direction ordering resolves ties without a floor. + if (metric < best) { + best = metric; + *selected = {true, j, proposed, delta, metric}; + } + } + if (ranked != nullptr) + std::stable_sort(ranked->begin(), ranked->end(), + [](const auto& a, const auto& b) { + return a.predicted_metric < b.predicted_metric; + }); + return true; +} +/// @brief Keeps complete eligible failed Newton entries for a fixed increment. +class SolidHhtNewtonHistory { + public: + bool Capture(const Vector& u, const std::vector& holding_mask, + const std::array& incoming_metrics, + std::uint32_t iteration, bool physical_pass) { + if (!Valid(u, holding_mask, incoming_metrics)) return false; + if (!physical_pass && + std::all_of(incoming_metrics.begin(), incoming_metrics.end(), + [](double metric) { return metric <= 1e-10; })) + entries_.push_back({u, holding_mask, incoming_metrics, iteration}); + return true; + } + bool Match(const Vector& proposed_u, + const std::vector& next_holding_mask, + const std::array& fresh_metrics, bool physical_pass, + SolidDynamicNewtonRepeatEvidence* evidence) const { + *evidence = {}; + evidence->history_entry_count = entries_.size(); + if (!Valid(proposed_u, next_holding_mask, fresh_metrics)) return false; + // A physically converged prospective endpoint must reach the normal gate. + if (physical_pass || + !std::all_of(fresh_metrics.begin(), fresh_metrics.end(), + [](double metric) { return metric <= 1e-10; })) + return true; + for (const auto& entry : entries_) { + evidence->changed_translations = 0; + evidence->changed_holding_mask = 0; + evidence->changed_correction_metrics = 0; + for (std::size_t i = 0; i < proposed_u.Size(); ++i) + evidence->changed_translations += proposed_u[i] != entry.u[i]; + for (std::size_t i = 0; i < next_holding_mask.size(); ++i) + evidence->changed_holding_mask += + next_holding_mask[i] != entry.holding_mask[i]; + for (std::size_t i = 0; i < 4; ++i) + evidence->changed_correction_metrics += + fresh_metrics[i] != entry.incoming_metrics[i]; + // Numeric correction magnitudes remain evidence. Their finite/PASS + // controls, rather than their exact magnitudes, affect the transition. + if (evidence->changed_translations == 0 && + evidence->changed_holding_mask == 0) { + evidence->matched_entry_iteration = entry.iteration; + evidence->repeat = true; + break; + } + } + return true; + } + std::size_t Size() const { return entries_.size(); } + + private: + struct Entry { + Vector u; + std::vector holding_mask; + std::array incoming_metrics; + std::uint32_t iteration; + }; + bool Valid(const Vector& u, const std::vector& mask, + const std::array& metrics) const { + if (mask.size() > u.Size() || + (!entries_.empty() && + (u.Size() != entries_.front().u.Size() || + mask.size() != entries_.front().holding_mask.size()))) + return false; + for (std::size_t i = 0; i < u.Size(); ++i) + if (!std::isfinite(u[i])) return false; + for (double metric : metrics) + if (!std::isfinite(metric) || metric < 0) return false; + return true; + } + std::vector entries_; +}; +/// @brief Checks the fresh complete Newton trial before representable fallback. +/// Historical scores are deliberately absent: a lower half of a two-cycle +/// still needs actual proposed-step evaluation. False keeps common Newton +/// apply. +template +Result TrySolidHhtActualStep( + bool residual_pass, const std::array& fresh_corrections, + Newton newton, Neighbor neighbor) { + if (residual_pass || !std::all_of(fresh_corrections.begin(), + fresh_corrections.end(), [](double metric) { + return std::isfinite(metric) && + metric <= 1e-10; + })) + return Result::Success(false); + auto keep_newton = newton(); + if (!keep_newton.HasValue()) + return Result::Failure(keep_newton.GetStatus()); + if (keep_newton.Value()) return Result::Success(false); + return neighbor(); +} +/// @brief Tries ranked predictions until actual endpoint evaluation improves. +/// Evaluator failure is terminal; a rejected prediction does not discard later +/// candidates, and success means only that a local trial improved. +template +Result TrySolidHhtRankedNeighbors( + const std::vector& ranked, + Evaluate evaluate) { + for (const auto& neighbor : ranked) { + auto actual = evaluate(neighbor); + if (!actual.HasValue()) return Result::Failure(actual.GetStatus()); + if (actual.Value()) return Result::Success(true); + } + return Result::Success(false); +} +} // namespace fesa::detail +#endif // FESA_ANALYSIS_SOLID_HHT_REPRESENTABILITY_H_ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3da1acd..f1a270e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -19,6 +19,8 @@ add_executable( unit/analysis/nonlinear_static_control_test.cpp unit/analysis/newton_iteration_log_test.cpp unit/analysis/solid_frequency_procedure_test.cpp + unit/analysis/solid_hht_dynamic_procedure_test.cpp + unit/analysis/solid_hht_temporal_verification_test.cpp unit/analysis/solid_static_procedure_test.cpp unit/analysis/solid_static_state_test.cpp unit/analysis/thermal_time_grid_test.cpp diff --git a/tests/unit/analysis/solid_hht_dynamic_procedure_test.cpp b/tests/unit/analysis/solid_hht_dynamic_procedure_test.cpp new file mode 100644 index 0000000..78c8fd4 --- /dev/null +++ b/tests/unit/analysis/solid_hht_dynamic_procedure_test.cpp @@ -0,0 +1,1512 @@ +#include "fesa/analysis/solid_hht_dynamic_procedure.h" + +#include + +#include +#include +#include +#include +#include + +#include "analysis/solid_hht_force_identity.h" +#include "analysis/solid_hht_high_low.h" +#include "analysis/solid_hht_representability.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/hex8_element.h" +#include "fesa/io/abaqus/domain_mapper.h" +#include "fesa/io/abaqus/input_reader.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/solvers/linear/mkl_pardiso_solver.h" +#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" +#include "solid_hht_test_fixture.h" + +namespace { +fesa::Result RunDynamics( + fesa::ModelDefinition d) { + auto domain = fesa::Domain::Create(std::move(d)); + if (!domain.HasValue()) + return fesa::Result::Failure( + domain.GetStatus()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + return fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .Execute(domain.Value()); +} +fesa::Status Injected() { + return fesa::Status::Failure(fesa::FailureCategory::kSolver, + {{fesa::Severity::kError, + "injected-backend", + {}, + "DYNAMIC", + {}, + "Injected backend failure."}}); +} +class ControlledTangent : public fesa::TangentSolver { + public: + enum Mode { kFail, kNan, kZero, kInvert }; + explicit ControlledTangent(Mode mode, unsigned successful_solves = 0) + : mode_(mode), successful_solves_(successful_solves) {} + fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override { + size_ = matrix.Rows(); + return real_.FactorizeFresh(matrix); + } + fesa::Status SolveCorrection(const fesa::Vector& rhs, + fesa::Vector& du) const override { + if (calls_++ < successful_solves_) return real_.SolveCorrection(rhs, du); + if (mode_ == kFail) return Injected(); + du = fesa::Vector(size_); + if (mode_ == kNan) du[0] = std::numeric_limits::quiet_NaN(); + if (mode_ == kInvert) du[0] = 2.0; + return fesa::Status::Ok(); + } + + private: + Mode mode_; + unsigned successful_solves_; + mutable unsigned calls_{0}; + std::size_t size_{0}; + fesa::MklPardisoTangentSolver real_; +}; +class FailedMass : public fesa::LinearSolver { + public: + explicit FailedMass(bool fail_factorization) + : factorization_(fail_factorization) {} + fesa::Status Factorize(const fesa::SparseMatrix&) override { + return factorization_ ? Injected() : fesa::Status::Ok(); + } + fesa::Status Solve(const fesa::Vector&, fesa::Vector&) const override { + return Injected(); + } + + private: + bool factorization_; +}; +fesa::ModelDefinition OneFree(bool nonlinear = false, double h = .01, + unsigned count = 20) { + auto d = solid_hht_test::Brick(h, count, nonlinear); + d.dynamic_step->boundaries.push_back({"N0", 2, 3, 0, {}}); + for (unsigned node = 1; node < 8; ++node) + d.dynamic_step->boundaries.push_back( + {"N" + std::to_string(node), 1, 3, 0, {}}); + d.dynamic_step->loads = {{"N0", 1, .01, {}}}; + return d; +} +} // namespace + +namespace { +// Independent uniform-acceleration brick initial state: rho*V/8=1/4 at each +// node. High a=1 and low a=2^-60 therefore give Ma=(1/4,2^-62). +fesa::SolidDynamicAcceptedFrameState PairedUniformInitial( + const fesa::Domain& domain) { + fesa::SolidDynamicAcceptedFrameState frame; + frame.is_initial = true; + for (std::size_t n = 0; n < domain.Nodes().size(); ++n) { + frame.nodal_internal_order.push_back(static_cast(n)); + frame.nodal_source_ids.push_back(domain.Nodes()[n].source_id); + } + for (auto* part : + {&frame.displacement, &frame.velocity, + &frame.acceleration, &frame.algorithmic_acceleration, + &frame.inertia_force, &frame.internal_force, + &frame.external_force, &frame.residual, + &frame.reaction, &frame.residual_scale, + &frame.low.displacement, &frame.low.velocity, + &frame.low.acceleration, &frame.low.algorithmic_acceleration, + &frame.low.inertia_force, &frame.low.internal_force, + &frame.low.external_force, &frame.low.residual, + &frame.low.reaction, &frame.low.residual_scale}) + *part = fesa::Vector(24); + for (std::size_t n = 0; n < 8; ++n) { + const auto i = 3 * n; + frame.acceleration[i] = frame.algorithmic_acceleration[i] = 1; + frame.low.acceleration[i] = frame.low.algorithmic_acceleration[i] = + std::ldexp(1.0, -60); + frame.inertia_force[i] = frame.external_force[i] = .25; + frame.low.inertia_force[i] = frame.low.external_force[i] = + std::ldexp(1.0, -62); + frame.residual_scale[i] = .5; + frame.low.residual_scale[i] = std::ldexp(1.0, -61); + } + frame.linear_momentum_rate[0] = 2; + frame.low.linear_momentum_rate[0] = std::ldexp(1.0, -59); + frame.angular_momentum_rate = {0, 1, -1}; + frame.low.angular_momentum_rate = {0, std::ldexp(1.0, -60), + -std::ldexp(1.0, -60)}; + fesa::ElementFactory factory; + auto runtime = factory.Create(domain.Elements()[0], domain); + const auto* solid = + dynamic_cast(runtime.Value().get()); + const auto zero = solid->EvaluateLinear(fesa::Vector(24)); + frame.solid_results.assign(zero.Value().points.begin(), + zero.Value().points.end()); + fesa::SolidDynamicIterationEvidence evidence; + evidence.residual_pass = evidence.correction_pass = evidence.update_pass = + true; + frame.iterations.push_back(evidence); + return frame; +} +} // namespace + +TEST(SolidHhtPairedState, InitialLowPartsSurviveTrialCopyAndRejectedCommit) { + auto domain = fesa::Domain::Create(solid_hht_test::Brick(.01, 1)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + auto initial = PairedUniformInitial(domain.Value()); + ASSERT_TRUE(state.Value().CommitAcceptedFrame(initial).IsOk()); + auto trial = state.Value().BeginTrialFromCommitted(); + ASSERT_TRUE(trial.HasValue()); + EXPECT_EQ(trial.Value().low.acceleration[0], std::ldexp(1.0, -60)); + EXPECT_EQ(trial.Value().low.residual_scale[0], std::ldexp(1.0, -61)); + EXPECT_EQ(trial.Value().low.linear_momentum_rate[0], std::ldexp(1.0, -59)); + EXPECT_EQ(trial.Value().low.angular_momentum_rate[1], std::ldexp(1.0, -60)); + trial.Value().low.acceleration[0] = 0; + EXPECT_EQ(state.Value().CommittedState().low.acceleration[0], + std::ldexp(1.0, -60)); + auto bad = initial; + bad.frame_index = bad.increment_index = 1; + bad.time = .01; + bad.is_initial = false; + bad.low.internal_force[0] = std::numeric_limits::quiet_NaN(); + EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk()); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 1U); + EXPECT_EQ(state.Value().CommittedState().low.acceleration[0], + std::ldexp(1.0, -60)); + EXPECT_EQ(state.Value().CommittedState().low.inertia_force[0], + std::ldexp(1.0, -62)); + EXPECT_EQ(state.Value().CommittedState().low.angular_momentum_rate[2], + -std::ldexp(1.0, -60)); +} + +TEST(SolidHhtPairedState, InvalidLowShapeFiniteNormalizationAndIdentitiesFail) { + auto domain = fesa::Domain::Create(solid_hht_test::Brick(.01, 1)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + const auto initial = PairedUniformInitial(domain.Value()); + for (unsigned mutation = 0; mutation < 8; ++mutation) { + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + auto bad = initial; + switch (mutation) { + case 0: + bad.low.velocity = fesa::Vector(23); + break; + case 1: + bad.low.acceleration[0] = std::numeric_limits::infinity(); + break; + case 2: + bad.low.acceleration[0] = std::ldexp(1.0, -51); + break; + case 3: + bad.low.algorithmic_acceleration[0] = 0; + break; + case 4: + bad.low.residual[0] = std::ldexp(1.0, -70); + break; + case 5: + bad.low.residual_scale[0] = 0; + break; + case 6: + bad.low.reaction[0] = std::ldexp(1.0, -70); + break; + case 7: + bad.low.angular_momentum[0] = std::numeric_limits::quiet_NaN(); + break; + } + EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk()) << mutation; + EXPECT_TRUE(state.Value().AcceptedFrames().empty()) << mutation; + EXPECT_FALSE(state.Value().BeginTrialFromCommitted().HasValue()) + << mutation; + } +} + +TEST(SolidHhtPairedKinematics, SubUlpPhysicalDisplacementRecoversNonzeroRates) { + const double epsilon = std::ldexp(1.0, -60); + auto p = fesa::MakeHhtParameters(0); + ASSERT_TRUE(p.HasValue()); + fesa::SolidDynamicAcceptedFrameState old; + old.displacement = fesa::Vector(1, 1); + old.low.displacement = fesa::Vector(1, epsilon); + old.velocity = old.low.velocity = fesa::Vector(1); + old.acceleration = old.low.acceleration = fesa::Vector(1); + old.algorithmic_acceleration = old.low.algorithmic_acceleration = + fesa::Vector(1); + auto predictor = fesa::detail::PredictSolidHhtHighLow(old, 1, p.Value()); + ASSERT_TRUE(predictor.HasValue()); + EXPECT_EQ(predictor.Value().displacement[0], 1); + EXPECT_EQ(predictor.Value().displacement_low[0], epsilon); + auto rates = fesa::detail::EvaluateSolidHhtHighLowRates( + fesa::Vector(1, 1), fesa::Vector(1, 2 * epsilon), predictor.Value(), + old.acceleration, old.low.acceleration, 1, p.Value()); + ASSERT_TRUE(rates.HasValue()); + // beta=1/4 and gamma=1/2 give exact dyadic physical oracles. + EXPECT_EQ(rates.Value().algorithmic_acceleration[0], 4 * epsilon); + EXPECT_EQ(rates.Value().acceleration[0], 4 * epsilon); + EXPECT_EQ(rates.Value().velocity[0], 2 * epsilon); + EXPECT_EQ(rates.Value().algorithmic_acceleration_low[0], 0); + EXPECT_EQ(rates.Value().acceleration_low[0], 0); + EXPECT_EQ(rates.Value().velocity_low[0], 0); + auto rounded_predictor = + fesa::PredictHht(old.displacement, old.velocity, + old.algorithmic_acceleration, 1, p.Value()); + auto rounded_rates = + fesa::EvaluateHhtRates(fesa::Vector(1, 1), rounded_predictor.Value(), + old.acceleration, 1, p.Value()); + ASSERT_TRUE(rounded_rates.HasValue()); + EXPECT_EQ(rounded_rates.Value().acceleration[0], 0); +} + +TEST(SolidHhtPairedKinematics, PhysicalAndAuxiliaryRatesKeepOldPhysicalTerm) { + const double epsilon = std::ldexp(1.0, -60); + auto p = fesa::MakeHhtParameters(-.25); + ASSERT_TRUE(p.HasValue()); + // h=8, beta=25/64 and gamma=3/4: denominator25, gamma*h=6. + fesa::detail::SolidHhtPairedPredictor predictor{ + fesa::Vector(1, 1), fesa::Vector(1), fesa::Vector(1), fesa::Vector(1)}; + auto rates = fesa::detail::EvaluateSolidHhtHighLowRates( + fesa::Vector(1, 1), fesa::Vector(1, 25 * epsilon), predictor, + fesa::Vector(1, 4 * epsilon), fesa::Vector(1), 8, p.Value()); + ASSERT_TRUE(rates.HasValue()); + EXPECT_EQ(rates.Value().algorithmic_acceleration[0], epsilon); + EXPECT_EQ(rates.Value().acceleration[0], 0); + EXPECT_EQ(rates.Value().velocity[0], 6 * epsilon); +} + +TEST(SolidHhtPairedKinematics, InvalidLowInputsFailWithoutChangingPredictor) { + auto p = fesa::MakeHhtParameters(0); + ASSERT_TRUE(p.HasValue()); + fesa::detail::SolidHhtPairedPredictor predictor{ + fesa::Vector(1, 1), fesa::Vector(1), fesa::Vector(1), fesa::Vector(1)}; + for (const double low : + {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity(), std::ldexp(1.0, -51)}) { + EXPECT_FALSE(fesa::detail::EvaluateSolidHhtHighLowRates( + fesa::Vector(1, 1), fesa::Vector(1, low), predictor, + fesa::Vector(1), fesa::Vector(1), 1, p.Value()) + .HasValue()); + EXPECT_EQ(predictor.displacement[0], 1); + EXPECT_EQ(predictor.displacement_low[0], 0); + } + EXPECT_FALSE(fesa::detail::EvaluateSolidHhtHighLowRates( + fesa::Vector(1, 1), fesa::Vector(0), predictor, + fesa::Vector(1), fesa::Vector(1), 1, p.Value()) + .HasValue()); + auto invalid_parameters = p.Value(); + invalid_parameters.beta = .5; + EXPECT_FALSE(fesa::detail::EvaluateSolidHhtHighLowRates( + fesa::Vector(1, 1), fesa::Vector(1), predictor, + fesa::Vector(1), fesa::Vector(1), 1, invalid_parameters) + .HasValue()); +} + +TEST(SolidHhtPairedKinematics, ConsistentMassKeepsSubUlpActionAndCancellation) { + fesa::SparsePattern pattern; + pattern.row_offsets = {0, 2, 4}; + pattern.column_indices = {0, 1, 0, 1}; + auto mass = fesa::SparseMatrix::FromCoo( + 2, 2, + {{0, 0, 2, 0, 0}, {0, 1, 1, 0, 1}, {1, 0, 1, 0, 2}, {1, 1, 2, 0, 3}}, + pattern); + ASSERT_TRUE(mass.HasValue()); + const double epsilon = std::ldexp(1.0, -60); + fesa::Vector high(2), low(2); + high[0] = 1; + high[1] = -2; + low[0] = epsilon; + low[1] = epsilon; + auto result = + fesa::detail::MultiplySolidHhtMassHighLow(mass.Value(), high, low); + ASSERT_TRUE(result.HasValue()); + EXPECT_EQ(result.Value().high[0], 3 * epsilon); + EXPECT_EQ(result.Value().low[0], 0); + EXPECT_EQ(result.Value().high[1], -3); + EXPECT_EQ(result.Value().low[1], 3 * epsilon); + EXPECT_EQ(mass.Value().Multiply(high)[0], 0); +} + +TEST(SolidHhtPairedKinematics, WeightedIdentityPreservesNonzeroOldResidualLow) { + fesa::SparsePattern pattern; + pattern.row_offsets = {0, 1}; + pattern.column_indices = {0}; + auto mass = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 1, 0, 0}}, pattern); + ASSERT_TRUE(mass.HasValue()); + const double epsilon = std::ldexp(1.0, -60); + fesa::SolidDynamicAcceptedFrameState old, next; + for (auto* frame : {&old, &next}) + for (auto* v : + {&frame->acceleration, &frame->low.acceleration, + &frame->algorithmic_acceleration, + &frame->low.algorithmic_acceleration, &frame->internal_force, + &frame->low.internal_force, &frame->external_force, + &frame->low.external_force, &frame->residual, &frame->low.residual}) + *v = fesa::Vector(1); + old.acceleration[0] = 4 * epsilon; + old.internal_force[0] = 3; + old.external_force[0] = 5; + old.residual[0] = -2; + old.low.residual[0] = 4 * epsilon; + next.algorithmic_acceleration[0] = epsilon; + next.internal_force[0] = 7; + next.external_force[0] = 2; + next.residual[0] = 5; + fesa::detail::SolidHhtPairedForceIdentity result; + ASSERT_TRUE(fesa::detail::EvaluateSolidHhtHighLowForceIdentity( + mass.Value(), 0, fesa::math_internal::HighLow(epsilon), next, old, -.25, + &result)); + EXPECT_EQ(result.original.high, 3.25); + EXPECT_EQ(result.original.low, epsilon); + EXPECT_EQ(result.solved.high, 3.25); + EXPECT_EQ(result.solved.low, epsilon); + EXPECT_EQ(result.defect.high, 0); + EXPECT_EQ(result.metric.high, 0); + // Dropping only the old residual low changes the approved solved force. + old.low.residual[0] = 0; + ASSERT_TRUE(fesa::detail::EvaluateSolidHhtHighLowForceIdentity( + mass.Value(), 0, fesa::math_internal::HighLow(epsilon), next, old, -.25, + &result)); + EXPECT_EQ(result.solved.low, 0); + EXPECT_EQ(result.defect.high, epsilon); +} + +namespace { +class RestoreArithmeticEnvironment { + public: + RestoreArithmeticEnvironment() : rounding_(std::fegetround()) { +#if defined(_M_X64) || defined(__SSE2__) + mxcsr_ = _mm_getcsr(); +#endif + } + ~RestoreArithmeticEnvironment() { + std::fesetround(rounding_); +#if defined(_M_X64) || defined(__SSE2__) + _mm_setcsr(mxcsr_); +#endif + } + + private: + int rounding_; +#if defined(_M_X64) || defined(__SSE2__) + unsigned mxcsr_; +#endif +}; +void MutateArithmeticEnvironment(bool flush) { +#if defined(_M_X64) || defined(__SSE2__) + if (flush) { + _mm_setcsr(_mm_getcsr() | 0x8040U); + return; + } +#else + (void)flush; +#endif + std::fesetround(FE_UPWARD); +} +class EnvironmentMass final : public fesa::LinearSolver { + public: + explicit EnvironmentMass(bool flush) : flush_(flush) {} + fesa::Status Factorize(const fesa::SparseMatrix& matrix) override { + return real_.Factorize(matrix); + } + fesa::Status Solve(const fesa::Vector& rhs, + fesa::Vector& result) const override { + const auto status = real_.Solve(rhs, result); + if (status.IsOk()) MutateArithmeticEnvironment(flush_); + return status; + } + + private: + bool flush_; + fesa::MklPardisoSolver real_; +}; +class EnvironmentTangent final : public fesa::TangentSolver { + public: + explicit EnvironmentTangent(bool flush) : flush_(flush) {} + fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override { + return real_.FactorizeFresh(matrix); + } + fesa::Status SolveCorrection(const fesa::Vector& rhs, + fesa::Vector& result) const override { + const auto status = real_.SolveCorrection(rhs, result); + if (status.IsOk() && calls_++ >= 2) MutateArithmeticEnvironment(flush_); + return status; + } + + private: + bool flush_; + mutable unsigned calls_{0}; + fesa::MklPardisoTangentSolver real_; +}; +} // namespace + +TEST(SolidHhtPairedState, BackendEnvironmentMutationFailsAndPreservesPairs) { + const auto baseline = RunDynamics(OneFree(false, .01, 3)); + ASSERT_TRUE(baseline.HasValue()); + const auto& expected = baseline.Value().state.AcceptedFrames()[1]; + for (bool flush : {false, true}) { + for (bool mutate_mass : {false, true}) { + SCOPED_TRACE(flush); + SCOPED_TRACE(mutate_mass); + auto domain = fesa::Domain::Create(OneFree(false, .01, 3)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver real_mass; + fesa::MklPardisoTangentSolver real_tangent; + EnvironmentMass changed_mass(flush); + EnvironmentTangent changed_tangent(flush); + auto status = fesa::Status::Ok(); + { + RestoreArithmeticEnvironment restore; + status = mutate_mass ? fesa::SolidHhtDynamicProcedureExecutor( + parallel, changed_mass, real_tangent) + .ExecuteInto(domain.Value(), state.Value()) + : fesa::SolidHhtDynamicProcedureExecutor( + parallel, real_mass, changed_tangent) + .ExecuteInto(domain.Value(), state.Value()); + EXPECT_FALSE(fesa::math_internal::ArithmeticEnvironmentSupported()); + } + ASSERT_FALSE(status.IsOk()); + ASSERT_FALSE(status.Diagnostics().empty()); + EXPECT_EQ(status.Diagnostics().front().code, + "solid-dynamic-arithmetic-environment"); + if (mutate_mass) { + EXPECT_TRUE(state.Value().AcceptedFrames().empty()); + } else { + ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U); + const auto& actual = state.Value().CommittedState(); + EXPECT_EQ(actual.displacement[0], expected.displacement[0]); + EXPECT_EQ(actual.low.displacement[0], expected.low.displacement[0]); + EXPECT_EQ(actual.acceleration[0], expected.acceleration[0]); + EXPECT_EQ(actual.low.acceleration[0], expected.low.acceleration[0]); + EXPECT_EQ(actual.external_work, expected.external_work); + EXPECT_EQ(actual.low.external_work, expected.low.external_work); + } + } + } +} + +TEST(SolidHhtDynamics, AllConstrainedRetainsPhysicalReactionAndEightPoints) { + auto d = solid_hht_test::Brick(); + d.dynamic_step->boundaries = {{"ALL", 1, 3, 0, {}}}; + d.dynamic_step->loads = {{"N0", 2, 7, {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto result = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus()); + ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 21U); + for (const auto& f : result.Value().state.AcceptedFrames()) { + EXPECT_EQ(f.displacement.Size(), 24U); + EXPECT_EQ(f.solid_results.size(), 8U); + EXPECT_DOUBLE_EQ(f.reaction[1], -7); + EXPECT_DOUBLE_EQ(f.acceleration[1], 0); + EXPECT_DOUBLE_EQ(f.external_work, 0); + EXPECT_EQ(f.iteration_count, 0U); + } +} + +TEST(SolidHhtDynamics, + ActualEndpointsUseDifferentCurvesKnotsInteriorsAndHolds) { + auto d = solid_hht_test::Brick(.125, 8); + d.dynamic_step->boundaries = {{"ALL", 1, 3, 0, {}}}; + d.amplitudes = {{"A", "A", {}, {{0, 0}, {.375, 1}, {.75, 2}}}, + {"B", "B", {}, {{0, 3}, {.25, -1}, {.5, 1}}}}; + d.dynamic_step->loads = {{"N0", 1, 2, {}, 0, "A", {}}, + {"N0", 1, 3, {}, 1, "B", {}}, + {"N1", 3, 7, {}}}; + auto result = RunDynamics(std::move(d)); + ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus()); + for (const auto& f : result.Value().state.AcceptedFrames()) { + const double a = f.time <= .375 ? f.time / .375 + : f.time <= .75 ? 1 + (f.time - .375) / .375 + : 2; + const double b = f.time <= .25 ? 3 - 16 * f.time + : f.time <= .5 ? -1 + 8 * (f.time - .25) + : 1; + EXPECT_NEAR(f.external_force[0], 2 * a + 3 * b, 1e-14); + EXPECT_DOUBLE_EQ(f.reaction[0], -f.external_force[0]); + EXPECT_DOUBLE_EQ(f.external_force[5], 7); + EXPECT_DOUBLE_EQ(f.time, static_cast(f.increment_index) * .125); + } +} + +TEST(SolidHhtDynamics, RampStartsAtZeroAndAuxiliaryAccelerationIsDistinct) { + auto d = OneFree(); + d.amplitudes = {{"A", "A", {}, {{0, 0}, {.1, 1}, {.2, 1}}}}; + d.dynamic_step->loads[0].amplitude_index = 0; + d.dynamic_step->loads[0].amplitude_name_raw = "A"; + auto result = RunDynamics(std::move(d)); + ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus()); + const auto& frames = result.Value().state.AcceptedFrames(); + EXPECT_DOUBLE_EQ(frames.front().acceleration[0], 0); + EXPECT_GT(std::abs(frames[1].acceleration[0] - + frames[1].algorithmic_acceleration[0]), + 1e-5); + for (std::size_t n = 1; n < frames.size(); ++n) { + const auto& f = frames[n]; + const auto& old = frames[n - 1]; + EXPECT_NEAR(f.algorithmic_acceleration[0], + .95 * f.acceleration[0] + .05 * old.acceleration[0], 1e-14); + EXPECT_LT(f.iterations.back().update_metric, 1e-12); + EXPECT_LE(f.iterations.back().residual_metric, 1e-10); + for (std::size_t i = 1; i < 24; ++i) { + EXPECT_DOUBLE_EQ(f.displacement[i], 0); + EXPECT_DOUBLE_EQ(f.velocity[i], 0); + EXPECT_DOUBLE_EQ(f.acceleration[i], 0); + EXPECT_DOUBLE_EQ(f.algorithmic_acceleration[i], 0); + } + } +} + +TEST(SolidHhtDynamics, + OneFreeTranslationHasAnalyticConsistentInertialSupportReactions) { + auto result = RunDynamics(OneFree()); + ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus()); + const auto d = solid_hht_test::Brick(); + constexpr double lambda = 4.8; + constexpr double mu = 4.8; + const double roundoff = 512 * std::numeric_limits::epsilon(); + const auto& frames = result.Value().state.AcceptedFrames(); + for (std::size_t n = 0; n < frames.size(); ++n) { + const auto& f = frames[n]; + for (unsigned b = 0; b < 8; ++b) { + std::array product{}, sign{}; + double scalar_mass = 2.0 / 216; + for (unsigned c = 0; c < 3; ++c) { + const bool same = d.nodes[b].coordinates[c] == 0; + product[c] = same ? 1.0 / 3 : 1.0 / 6; + sign[c] = same ? -1 : 1; + scalar_mass *= same ? 2 : 1; + } + for (unsigned c = 0; c < 3; ++c) { + double stiffness; + double stiffness_term_scale; + if (c == 0) { + stiffness = (lambda + 2 * mu) * (-sign[0]) * product[1] * product[2] + + mu * (-sign[1]) * product[0] * product[2] + + mu * (-sign[2]) * product[0] * product[1]; + stiffness_term_scale = (lambda + 2 * mu) * product[1] * product[2] + + mu * product[0] * product[2] + + mu * product[0] * product[1]; + } else { + const unsigned other = c == 1 ? 2 : 1; + stiffness = -.25 * (lambda * sign[c] + mu * sign[0]) * product[other]; + stiffness_term_scale = .25 * (lambda + mu) * product[other]; + } + const double inertia = c == 0 ? scalar_mass * f.acceleration[0] : 0; + const double internal = stiffness * f.displacement[0]; + const double applied = b == 0 && c == 0 ? .01 : 0; + const double reaction = + b == 0 && c == 0 ? 0 : inertia + internal - applied; + const double scale = + std::abs(inertia) + + stiffness_term_scale * std::abs(f.displacement[0]) + + std::abs(applied); + EXPECT_NEAR(f.reaction[3 * b + c], reaction, roundoff * scale); + } + } + if (n != 0) { + const auto& old = frames[n - 1]; + const double mass = 2.0 / 27; + const double weighted = + mass * f.algorithmic_acceleration[0] + + .95 * (f.internal_force[0] - f.external_force[0]) + + .05 * (old.internal_force[0] - old.external_force[0]); + const double endpoint_identity = + .95 * f.residual[0] + .05 * old.residual[0]; + const double scale = std::abs(mass * f.algorithmic_acceleration[0]) + + .95 * (std::abs(f.internal_force[0]) + + std::abs(f.external_force[0])) + + .05 * (std::abs(old.internal_force[0]) + + std::abs(old.external_force[0])); + EXPECT_NEAR(weighted, endpoint_identity, roundoff * scale); + EXPECT_NEAR(f.iterations.back().weighted_residual_numerator, + std::abs(weighted), roundoff * scale); + } + } +} + +TEST(SolidHhtDynamics, + BackendNanNonconvergenceAndInversionPreserveAcceptedInitial) { + for (const auto mode : + {ControlledTangent::kFail, ControlledTangent::kNan, + ControlledTangent::kZero, ControlledTangent::kInvert}) { + SCOPED_TRACE(mode); + auto d = OneFree(true); + // Moving node0 past the opposite x corner creates a genuinely folded cell. + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + ASSERT_TRUE(state.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + ControlledTangent tangent(mode); + auto status = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + EXPECT_FALSE(status.IsOk()); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 1U); + EXPECT_EQ(state.Value().CommittedState().increment_index, 0U); + EXPECT_DOUBLE_EQ(state.Value().CommittedState().displacement[0], 0); + EXPECT_NEAR(state.Value().CommittedState().acceleration[0], + .01 / (2.0 / 27.0), 1e-14); + ASSERT_FALSE(state.Value().IterationRecords().empty()); + EXPECT_EQ(state.Value().IterationRecords().back().stage, "failed"); + EXPECT_FALSE(state.Value().IterationRecords().back().failure_code.empty()); + if (mode == ControlledTangent::kZero) { + const auto& failed = state.Value().IterationRecords().back(); + EXPECT_EQ(failed.evidence.iteration_index, 100U); + EXPECT_EQ(failed.failure_code, "solid-dynamic-nonconvergence"); + // Paired physical evidence replaces the retired Float64 shadow sum. + EXPECT_FALSE(failed.failure_force_sum.has_value()); + const double expected_inertia = + .01 * (-(.5 - .275625) / .275625 - .05) / .95; + EXPECT_NEAR(failed.evidence.worst_force_terms[1], expected_inertia, + 32 * std::numeric_limits::epsilon() * + std::abs(expected_inertia)); + EXPECT_DOUBLE_EQ(failed.evidence.worst_force_terms[2], 0); + EXPECT_DOUBLE_EQ(failed.evidence.worst_force_terms[3], .01); + EXPECT_FALSE(failed.evidence.residual_pass); + EXPECT_TRUE(failed.evidence.update_pass); + EXPECT_TRUE(failed.evidence.correction_pass); + for (double value : failed.evidence.worst_force_terms) + EXPECT_TRUE(std::isfinite(value)); + for (std::size_t i = 0; + i < state.Value().CommittedState().displacement.Size(); ++i) + EXPECT_DOUBLE_EQ(state.Value().CommittedState().low.displacement[i], 0); + } + } +} + +TEST(SolidHhtDynamics, + LateBackendFailurePreservesActualCommittedMotionAndHistory) { + auto domain = fesa::Domain::Create(OneFree()); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + ASSERT_TRUE(state.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + ControlledTangent tangent(ControlledTangent::kFail, 2); + auto status = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + EXPECT_FALSE(status.IsOk()); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U); + EXPECT_EQ(state.Value().CommittedState().increment_index, 1U); + EXPECT_GT(state.Value().CommittedState().displacement[0], 0); + const auto committed = state.Value().CommittedState(); + auto bad = committed; + bad.frame_index = 2; + bad.increment_index = 2; + bad.time = .02; + bad.solid_results[0].cauchy_stress[0] = + std::numeric_limits::infinity(); + EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk()); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U); + EXPECT_DOUBLE_EQ(state.Value().CommittedState().displacement[0], + committed.displacement[0]); + EXPECT_DOUBLE_EQ(state.Value().CommittedState().strain_energy, + committed.strain_energy); + EXPECT_EQ(state.Value().CommittedState().solid_results[0].cauchy_stress, + committed.solid_results[0].cauchy_stress); +} + +TEST(SolidHhtDynamics, + NonfiniteWeightedEvidenceCannotContaminateAcceptedHistory) { + auto domain = fesa::Domain::Create(OneFree()); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + ControlledTangent tangent(ControlledTangent::kFail, 2); + EXPECT_FALSE(fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()) + .IsOk()); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U); + const auto committed = state.Value().CommittedState(); + for (double invalid : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity()}) { + for (unsigned field = 0; field < 8; ++field) { + auto bad = committed; + bad.frame_index = bad.increment_index = 2; + bad.time = .02; + auto& evidence = bad.iterations.back(); + if (field == 0) evidence.weighted_residual_numerator = invalid; + if (field == 1) evidence.recombined_residual_numerator = invalid; + if (field == 2) evidence.old_residual_numerator = invalid; + if (field == 3) evidence.weighted_equivalence_metric = invalid; + if (field >= 4) evidence.worst_weighted_terms[field - 4] = invalid; + EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk()); + EXPECT_EQ(state.Value().AcceptedFrames().size(), 2U); + EXPECT_EQ(state.Value().CommittedState().solid_results[0].cauchy_stress, + committed.solid_results[0].cauchy_stress); + EXPECT_DOUBLE_EQ(state.Value().CommittedState().displacement[0], + committed.displacement[0]); + } + } +} + +TEST(SolidHhtDynamics, + AlphaZeroLinearUnforcedTailConservesEnergyAndTrapezoidWork) { + auto d = OneFree(false, .01, 100); + d.amplitudes = {{"A", "A", {}, {{0, 0}, {.1, 1}, {.2, 0}, {1, 0}}}}; + d.dynamic_step->loads[0].amplitude_index = 0; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto status = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteIntoForVerification(domain.Value(), state.Value(), 0); + ASSERT_TRUE(status.IsOk()) << solid_hht_test::Message(status); + const auto& frames = state.Value().AcceptedFrames(); + const double tail = frames[20].kinetic_energy + frames[20].strain_energy; + ASSERT_GT(tail, 1e-7); + double work = 0; + for (std::size_t n = 1; n < frames.size(); ++n) { + work += .5 * + (frames[n - 1].external_force[0] + frames[n].external_force[0]) * + (frames[n].displacement[0] - frames[n - 1].displacement[0]); + EXPECT_NEAR(frames[n].external_work, work, 1e-16); + EXPECT_NEAR(frames[n].energy_balance_error, 0, 2e-12 * tail); + if (n >= 20) + EXPECT_NEAR(frames[n].kinetic_energy + frames[n].strain_energy, tail, + 2e-12 * tail); + } +} + +TEST(SolidHhtDynamics, InitialMassBackendFailuresPreserveEmptyHistory) { + for (bool factorization : {false, true}) { + auto domain = fesa::Domain::Create(OneFree()); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + ASSERT_TRUE(state.HasValue()); + fesa::SerialParallelFor parallel; + FailedMass mass(factorization); + fesa::MklPardisoTangentSolver tangent; + const auto status = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + EXPECT_FALSE(status.IsOk()); + EXPECT_TRUE(state.Value().AcceptedFrames().empty()); + EXPECT_FALSE(state.Value().BeginTrialFromCommitted().HasValue()); + const auto& records = state.Value().IterationRecords(); + ASSERT_EQ(records.size(), factorization ? 1U : 2U); + if (!factorization) { + EXPECT_EQ(records.front().stage, "initial-mass"); + EXPECT_EQ(records.front().outcome, "residual-refinement"); + EXPECT_EQ(records.front().evidence.iteration_index, 0U); + EXPECT_FALSE(records.front().evidence.residual_pass); + } + EXPECT_EQ(records.back().stage, "failed"); + EXPECT_EQ(records.back().failure_code, "injected-backend"); + } +} + +TEST(SolidHhtDynamics, AllFreeNoForceIsExactlyAtRestForLinearAndTl) { + for (bool nonlinear : {false, true}) { + auto result = RunDynamics(solid_hht_test::Brick(.01, 3, nonlinear)); + ASSERT_TRUE(result.HasValue()) + << solid_hht_test::Message(result.GetStatus()); + for (const auto& frame : result.Value().state.AcceptedFrames()) { + EXPECT_EQ(frame.iteration_count, 0U); + EXPECT_DOUBLE_EQ(frame.kinetic_energy, 0); + EXPECT_DOUBLE_EQ(frame.strain_energy, 0); + for (std::size_t i = 0; i < frame.displacement.Size(); ++i) { + EXPECT_DOUBLE_EQ(frame.displacement[i], 0); + EXPECT_DOUBLE_EQ(frame.velocity[i], 0); + EXPECT_DOUBLE_EQ(frame.acceleration[i], 0); + EXPECT_DOUBLE_EQ(frame.algorithmic_acceleration[i], 0); + EXPECT_DOUBLE_EQ(frame.residual[i], 0); + EXPECT_DOUBLE_EQ(frame.reaction[i], 0); + } + } + } +} + +void CheckOriginalSourceDeck(const std::string& name) { + auto root = std::filesystem::path{__FILE__}.parent_path(); + for (unsigned i = 0; i < 3; ++i) root = root.parent_path(); + { + auto syntax = + fesa::AbaqusInputReader{}.Read(root / "reference/hexa" / name); + ASSERT_TRUE(syntax.HasValue()); + auto domain = fesa::AbaqusDomainMapper{}.Map(syntax.Value()); + ASSERT_TRUE(domain.HasValue()) + << solid_hht_test::Message(domain.GetStatus()); + fesa::TbbParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + ASSERT_TRUE(state.HasValue()); + auto status = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + if (!status.IsOk()) { + const auto& r = state.Value().IterationRecords().back(); + std::cout << name << " accepted=" << state.Value().AcceptedFrames().size() + << " increment=" << r.increment_index + << " iteration=" << r.evidence.iteration_index + << " metric=" << r.evidence.residual_metric + << " update=" << r.evidence.update_metric + << " worst=" << r.evidence.worst_full_dof << '\n'; + std::cout.precision(17); + std::cout << "Worst terms r,Ma,fint,F,S: "; + for (const double term : r.evidence.worst_force_terms) + std::cout << term << ' '; + std::cout << '\n'; + std::cout << "Original,recombined,old residual,equivalence scale: "; + for (double term : r.evidence.worst_weighted_terms) + std::cout << term << ' '; + std::cout << " equivalence_metric=" + << r.evidence.weighted_equivalence_metric << '\n'; + for (const auto& coordinate : r.failure_coordinates) { + std::cout + << "Coordinate element=" << coordinate.element_index + << " row=" << coordinate.local_row + << " u=" << coordinate.displacement + << " previous_u=" << coordinate.previous_displacement + << " fresh_du=" << coordinate.fresh_correction + << " node0_u=" << coordinate.reference_displacement + << " previous_node0_u=" + << coordinate.previous_reference_displacement + << " relative_high=" << coordinate.relative_high + << " previous_relative_high=" << coordinate.previous_relative_high + << " u_changed=" + << (coordinate.displacement != coordinate.previous_displacement) + << " relative_changed=" + << (coordinate.relative_high != coordinate.previous_relative_high) + << '\n'; + } + if (r.failure_force_sum) { + const auto& sums = *r.failure_force_sum; + std::cout << "Failure force sums " + "internal,compensated,abs_local,inertia,compensated,abs_" + "products: " + << sums.assembled_internal << ' ' << sums.compensated_internal + << ' ' << sums.absolute_local_internal_sum << ' ' + << sums.assembled_inertia << ' ' << sums.compensated_inertia + << ' ' << sums.absolute_inertia_product_sum << '\n'; + } + double shadow_delta = 0, local_scale = 0; + for (const auto& dot : r.failure_linear_dot) { + std::cout << "Local dot element=" << dot.element_index + << " row=" << dot.local_row + << " sum_abs_products=" << dot.absolute_product_sum + << " kernel_force=" << dot.kernel_force + << " rounded_relative_compensated_force=" + << dot.compensated_force + << " delta=" << dot.compensated_force - dot.kernel_force + << '\n'; + shadow_delta += dot.compensated_force - dot.kernel_force; + local_scale += dot.absolute_product_sum; + } + std::cout << "Assembled rounded-relative shadow delta=" << shadow_delta + << " local_product_scale=" << local_scale << '\n'; + for (const auto& rec : state.Value().IterationRecords()) { + if (rec.newton_repeat && + (rec.newton_repeat->repeat || rec.evidence.iteration_index < 6 || + rec.evidence.iteration_index > 94)) { + const auto& probe = *rec.newton_repeat; + std::cout << "Source exact repeat " << rec.increment_index << ' ' + << rec.evidence.iteration_index + << " history_entries=" << probe.history_entry_count + << " matched_entry=" << probe.matched_entry_iteration + << " u_changed=" << probe.changed_translations + << " holding_mask_changed=" << probe.changed_holding_mask + << " correction_changed=" + << probe.changed_correction_metrics << '\n'; + } + if (rec.evidence.iteration_index < 6 || + rec.evidence.iteration_index > 94) + std::cout << "Source trace " << rec.increment_index << ' ' + << rec.evidence.iteration_index << ' ' << rec.outcome + << " ratio=" << rec.evidence.residual_metric + << " correction=" << rec.evidence.correction_metrics[0] + << '\n'; + } + } + ASSERT_TRUE(status.IsOk()) + << name << ' ' << solid_hht_test::Message(status); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 101U); + EXPECT_EQ(dofs.Value().FullDofCount(), 1326U); + const auto& final = state.Value().CommittedState(); + EXPECT_DOUBLE_EQ(final.time, 10); + EXPECT_EQ(final.solid_results.size(), 2400U); + std::cout << name + << " accepted=101 final_u_norm=" << final.displacement.Norm() + << " final_force_metric=" + << final.iterations.back().residual_metric + << " final_update_metric=" + << final.iterations.back().update_metric << '\n'; + } +} + +TEST(SolidHhtSourceSmoke, ReadOnlyOriginalLinearAcceptsAll101Endpoints) { + CheckOriginalSourceDeck("hexa_dynlinear.inp"); +} +TEST(SolidHhtSourceSmoke, ReadOnlyOriginalNonlinearAcceptsAll101Endpoints) { + CheckOriginalSourceDeck("hexa_dynnonlinear.inp"); +} + +TEST(SolidHhtDynamics, + WeightedIdentityRetainsNonzeroOldResidualAndEffectiveTangent) { + const fesa::SparsePattern pattern{{0, 2, 4}, {0, 1, 0, 1}}; + auto mass = fesa::SparseMatrix::FromCoo( + 2, 2, + {{0, 0, 2, 0, 0}, {0, 1, .3, 0, 1}, {1, 0, .3, 0, 2}, {1, 1, 1.5, 0, 3}}, + pattern); + ASSERT_TRUE(mass.HasValue()); + constexpr double alpha = -.05, beta = .275625, h = .2; + fesa::SolidDynamicAcceptedFrameState old; + old.acceleration = fesa::Vector(2); + old.acceleration[0] = .4; + old.acceleration[1] = -.7; + old.internal_force = fesa::Vector(2); + old.internal_force[0] = 1; + old.internal_force[1] = -2; + old.external_force = fesa::Vector(2); + old.external_force[0] = .1; + old.external_force[1] = .2; + old.residual = mass.Value().Multiply(old.acceleration); + for (unsigned i = 0; i < 2; ++i) + old.residual[i] += old.internal_force[i] - old.external_force[i]; + ASSERT_NE(old.residual[0], 0); + ASSERT_NE(old.residual[1], 0); + const std::array u{.12, -.03}, predictor{.09, -.01}, + direction{.7, -.4}; + for (double cubic : {0.0, 11.0}) { + const auto evaluate = [&](double epsilon) { + fesa::SolidDynamicAcceptedFrameState next; + next.acceleration = fesa::Vector(2); + next.algorithmic_acceleration = fesa::Vector(2); + next.internal_force = fesa::Vector(2); + next.external_force = fesa::Vector(2); + next.external_force[0] = .3; + next.external_force[1] = -.8; + const double x = u[0] + epsilon * direction[0]; + const double y = u[1] + epsilon * direction[1]; + for (unsigned i = 0; i < 2; ++i) { + const double q = i == 0 ? x : y; + next.algorithmic_acceleration[i] = (q - predictor[i]) / (beta * h * h); + next.acceleration[i] = + (next.algorithmic_acceleration[i] + alpha * old.acceleration[i]) / + (1 + alpha); + } + next.internal_force[0] = 7 * x + 2 * y + cubic * x * x * x; + next.internal_force[1] = 2 * x + 5 * y + cubic * y * y * y; + next.residual = mass.Value().Multiply(next.acceleration); + const auto MA = mass.Value().Multiply(next.algorithmic_acceleration); + std::array identities{}; + for (unsigned i = 0; i < 2; ++i) { + next.residual[i] += next.internal_force[i] - next.external_force[i]; + EXPECT_TRUE(fesa::detail::EvaluateSolidHhtForceIdentity( + mass.Value(), i, MA[i], next, old, alpha, &identities[i])); + EXPECT_LE(identities[i].metric, 1e-12); + EXPECT_NEAR( + identities[i].original, identities[i].solved, + 512 * std::numeric_limits::epsilon() * identities[i].scale); + EXPECT_GT( + std::abs(identities[i].solved - (1 + alpha) * next.residual[i]), + .01); + } + return identities; + }; + constexpr double epsilon = 1e-6; + const auto left = evaluate(-epsilon), right = evaluate(epsilon); + const std::array tangent_direction{ + (2 * direction[0] + .3 * direction[1]) / (beta * h * h) + + (1 + alpha) * ((7 + 3 * cubic * u[0] * u[0]) * direction[0] + + 2 * direction[1]), + (.3 * direction[0] + 1.5 * direction[1]) / (beta * h * h) + + (1 + alpha) * (2 * direction[0] + + (5 + 3 * cubic * u[1] * u[1]) * direction[1])}; + for (unsigned i = 0; i < 2; ++i) + EXPECT_NEAR((right[i].solved - left[i].solved) / (2 * epsilon), + tangent_direction[i], 1e-9 * std::abs(tangent_direction[i])); + } +} + +TEST(SolidHhtDynamics, + WeightedIdentityRejectsNonfiniteAndOverflowBeforeNormReduction) { + const fesa::SparsePattern pattern{{0, 1}, {0}}; + auto mass = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 2, 0, 0}}, pattern); + ASSERT_TRUE(mass.HasValue()); + fesa::SolidDynamicAcceptedFrameState old, next; + for (auto* frame : {&old, &next}) + for (auto* v : + {&frame->acceleration, &frame->algorithmic_acceleration, + &frame->internal_force, &frame->external_force, &frame->residual}) + *v = fesa::Vector(1); + fesa::detail::SolidHhtForceIdentity identity; + EXPECT_TRUE(fesa::detail::EvaluateSolidHhtForceIdentity( + mass.Value(), 0, 0, next, old, -.05, &identity)); + EXPECT_EQ(identity.scale, 0); + EXPECT_EQ(identity.metric, 0); + for (double invalid : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity()}) { + EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity( + mass.Value(), 0, invalid, next, old, -.05, &identity)); + for (unsigned field = 0; field < 8; ++field) { + auto bad = next; + fesa::Vector* vectors[]{&bad.acceleration, &bad.algorithmic_acceleration, + &bad.internal_force, &bad.external_force, + &bad.residual}; + if (field < 5) (*vectors[field])[0] = invalid; + auto previous = old; + if (field == 5) previous.acceleration[0] = invalid; + if (field == 6) previous.residual[0] = invalid; + if (field == 7) previous.external_force[0] = invalid; + EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity( + mass.Value(), 0, 0, bad, previous, -.05, &identity)); + } + } + const double huge = std::numeric_limits::max(); + next.internal_force[0] = next.external_force[0] = huge; + EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity( + mass.Value(), 0, 0, next, old, -.05, &identity)); + next.internal_force[0] = next.external_force[0] = 0; + next.algorithmic_acceleration[0] = huge; + EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity( + mass.Value(), 0, 0, next, old, -.05, &identity)); +} + +TEST(SolidHhtDynamics, + FullFreeNewtonRecoversAllComponentsAfterCoupledTrialPerturbation) { + class CoupledKick final : public fesa::TangentSolver { + public: + fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override { + rows_ = matrix.Rows(); + full_size_all_calls = full_size_all_calls && rows_ == 24; + if (kicked) solved_after_kick = true; + return real_.FactorizeFresh(matrix); + } + fesa::Status SolveCorrection(const fesa::Vector& rhs, + fesa::Vector& du) const override { + auto status = real_.SolveCorrection(rhs, du); + if (status.IsOk() && !kicked) { + // Deliberately perturb one transverse correction; the next full + // physical evaluation must recover every coupled original component. + du[1] += 1e-6; + kicked = true; + } + return status; + } + mutable bool kicked{false}; + bool full_size_all_calls{true}; + bool solved_after_kick{false}; + + private: + std::size_t rows_{0}; + fesa::MklPardisoTangentSolver real_; + } tangent; + auto d = solid_hht_test::Brick(.01, 1); + d.dynamic_step->loads = {{"ALL", 1, .25, {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + auto status = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + EXPECT_TRUE(tangent.kicked); + EXPECT_TRUE(tangent.full_size_all_calls); + EXPECT_TRUE(tangent.solved_after_kick); + ASSERT_TRUE(status.IsOk()) << solid_hht_test::Message(status); + const double roundoff = 512 * std::numeric_limits::epsilon(); + for (const auto& frame : state.Value().AcceptedFrames()) { + EXPECT_LE(frame.iterations.back().residual_metric, 1e-10); + EXPECT_LE(frame.iterations.back().update_metric, 1e-12); + for (unsigned n = 0; n < 8; ++n) + EXPECT_NEAR(frame.displacement[3 * n], .5 * frame.time * frame.time, + roundoff * .5 * frame.time * frame.time); + } +} + +TEST(SolidHhtFloat64History, + RepresentableNeighborUsesHeldCouplingAndOriginalFullForceMetric) { + const fesa::SparsePattern pattern{{0, 2, 4}, {0, 1, 0, 1}}; + auto mass = fesa::SparseMatrix::FromCoo(2, 2, + {{0, 0, 1000, 0, 0}, + {0, 1, 0, 0, 1}, + {1, 0, 0, 0, 2}, + {1, 1, 1000, 0, 3}}, + pattern); + auto tangent = fesa::SparseMatrix::FromCoo(2, 2, + {{0, 0, 1e6, 0, 0}, + {0, 1, 1e6, 0, 1}, + {1, 0, 1e6, 0, 2}, + {1, 1, 2e6, 0, 3}}, + pattern); + ASSERT_TRUE(mass.HasValue()); + ASSERT_TRUE(tangent.HasValue()); + fesa::SolidDynamicAcceptedFrameState frame; + for (auto* v : + {&frame.displacement, &frame.inertia_force, &frame.internal_force, + &frame.external_force, &frame.residual, &frame.residual_scale}) + *v = fesa::Vector(2); + frame.displacement[0] = 1.5; + frame.displacement[1] = .5; + frame.inertia_force[0] = -.5; + frame.internal_force[0] = .500000000102; + frame.inertia_force[1] = -50; + frame.internal_force[1] = 50; + for (unsigned i = 0; i < 2; ++i) { + frame.residual[i] = frame.inertia_force[i] + frame.internal_force[i]; + frame.residual_scale[i] = + std::abs(frame.inertia_force[i]) + std::abs(frame.internal_force[i]); + } + const double initial = std::abs(frame.residual[0]) / frame.residual_scale[0]; + ASSERT_GT(initial, 1e-10); + ASSERT_EQ(frame.residual[1], 0); // Equilibrated coupled coordinate is held. + fesa::detail::SolidHhtRepresentableNeighbor selected; + ASSERT_TRUE(fesa::detail::SelectSolidHhtRepresentableNeighbor( + mass.Value(), tangent.Value(), {0, 1}, frame, 1, -.05, .275625, + &selected)); + ASSERT_TRUE(selected.found); + EXPECT_EQ(selected.full_dof, 1U); + EXPECT_EQ(selected.displacement, + std::nextafter(.5, -std::numeric_limits::infinity())); + EXPECT_EQ(selected.delta, selected.displacement - frame.displacement[1]); + EXPECT_LT(selected.predicted_metric, initial); + EXPECT_LE(selected.predicted_metric, 1e-10); + // Independent separate inertia/internal-force update checks the original + // physical component score, including the previously equilibrated row. + double exact_metric = 0; + for (unsigned i = 0; i < 2; ++i) { + const double ma = frame.inertia_force[i] + + (i == 1 ? 1000 : 0) * selected.delta / (.95 * .275625); + const double f = + frame.internal_force[i] + (i == 0 ? 1e6 : 2e6) * selected.delta; + exact_metric = (std::max)(exact_metric, + std::abs(ma + f) / (std::abs(ma) + std::abs(f))); + } + EXPECT_EQ(selected.predicted_metric, exact_metric); + auto alternative_tangent = fesa::SparseMatrix::FromCoo(2, 2, + {{0, 0, 4e5, 0, 0}, + {0, 1, 1e6, 0, 1}, + {1, 0, 1e6, 0, 2}, + {1, 1, 2e6, 0, 3}}, + pattern); + ASSERT_TRUE(alternative_tangent.HasValue()); + std::vector ranked; + ASSERT_TRUE(fesa::detail::SelectSolidHhtRepresentableNeighbor( + mass.Value(), alternative_tangent.Value(), {0, 1}, frame, 1, -.05, + .275625, &selected, &ranked)); + ASSERT_EQ(ranked.size(), 2U); + EXPECT_EQ(ranked[0].full_dof, 0U); + EXPECT_EQ(ranked[1].full_dof, 1U); + EXPECT_LT(ranked[0].predicted_metric, ranked[1].predicted_metric); + EXPECT_LT(ranked[1].predicted_metric, initial); + // Independent actual-score fixture models a rounded force evaluation whose + // best linear prediction worsens the physical ratio, while another improves. + // Neither prediction is accepted as an actual force result. + const std::array actual_scores{initial * 1.1, initial * .8}; + std::vector tried; + auto accepted = + fesa::detail::TrySolidHhtRankedNeighbors(ranked, [&](const auto& trial) { + tried.push_back(trial.full_dof); + return fesa::Result::Success(actual_scores[trial.full_dof] < + initial); + }); + ASSERT_TRUE(accepted.HasValue()); + EXPECT_TRUE(accepted.Value()); + EXPECT_EQ(tried, (std::vector{0, 1})); + tried.clear(); + auto invalid = + fesa::detail::TrySolidHhtRankedNeighbors(ranked, [&](const auto& trial) { + tried.push_back(trial.full_dof); + return fesa::Result::Failure( + fesa::Status::Failure(fesa::FailureCategory::kSolver, {})); + }); + EXPECT_FALSE(invalid.HasValue()); + EXPECT_EQ(tried, (std::vector{0})); +} + +TEST(SolidHhtDynamics, + RepresentableNeighborKeepsExactZeroAndRejectsInvalidArithmetic) { + const fesa::SparsePattern pattern{{0, 1}, {0}}; + auto mass = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 1, 0, 0}}, pattern); + auto tangent = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 2, 0, 0}}, pattern); + ASSERT_TRUE(mass.HasValue()); + ASSERT_TRUE(tangent.HasValue()); + fesa::SolidDynamicAcceptedFrameState frame; + for (auto* v : + {&frame.displacement, &frame.inertia_force, &frame.internal_force, + &frame.external_force, &frame.residual, &frame.residual_scale}) + *v = fesa::Vector(1); + fesa::detail::SolidHhtRepresentableNeighbor selected; + EXPECT_TRUE(fesa::detail::SelectSolidHhtRepresentableNeighbor( + mass.Value(), tangent.Value(), {0}, frame, .1, -.05, .275625, &selected)); + EXPECT_FALSE(selected.found); + for (double invalid : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity()}) { + auto bad = frame; + bad.internal_force[0] = invalid; + EXPECT_FALSE(fesa::detail::SelectSolidHhtRepresentableNeighbor( + mass.Value(), tangent.Value(), {0}, bad, .1, -.05, .275625, &selected)); + } + EXPECT_FALSE(fesa::detail::SelectSolidHhtRepresentableNeighbor( + mass.Value(), tangent.Value(), {0}, frame, 0, -.05, .275625, &selected)); +} + +TEST(SolidHhtFloat64History, + ActualNewtonStallAtLowHalfOfTwoCycleTriesFallback) { + // The current physical force ratio improved from the previous iterate, but + // the next complete Newton trial returns to the larger ratio. Historical + // progress must not suppress actual evaluation or representable fallback. + constexpr double previous_metric = 1.6e-10; + constexpr double current_metric = 1.4e-10; + constexpr double actual_neighbor_metric = 1.2e-10; + ASSERT_LT(current_metric, previous_metric); + fesa::Vector current(2), previous(2); + current[0] = .5; + previous[0] = 1; + const std::vector mask(2, false); + const std::array metrics{1e-20, 2e-20, 3e-20, 4e-20}; + fesa::detail::SolidHhtNewtonHistory history; + ASSERT_TRUE(history.Capture(previous, mask, metrics, 0, false)); + ASSERT_TRUE(history.Capture(current, mask, metrics, 1, false)); + std::vector evaluated; + auto selected = fesa::detail::TrySolidHhtActualStep( + false, {1e-20, 2e-20, 3e-20, 4e-20}, + [&] { + evaluated.push_back("complete-newton"); + fesa::SolidDynamicNewtonRepeatEvidence repeat; + EXPECT_TRUE(history.Match(previous, mask, metrics, false, &repeat)); + return fesa::Result::Success(!repeat.repeat); + }, + [&] { + evaluated.push_back("actual-neighbor"); + return fesa::Result::Success(actual_neighbor_metric < + current_metric); + }); + ASSERT_TRUE(selected.HasValue()); + EXPECT_TRUE(selected.Value()); + EXPECT_EQ(evaluated, + (std::vector{"complete-newton", "actual-neighbor"})); + evaluated.clear(); + auto productive = fesa::detail::TrySolidHhtActualStep( + false, {1e-20, 2e-20, 3e-20, 4e-20}, + [&] { + evaluated.push_back("complete-newton"); + return fesa::Result::Success(true); + }, + [&] { + evaluated.push_back("actual-neighbor"); + return fesa::Result::Success(true); + }); + ASSERT_TRUE(productive.HasValue()); + EXPECT_FALSE(productive.Value()); // Keep complete Newton common application. + EXPECT_EQ(evaluated, (std::vector{"complete-newton"})); +} + +TEST(SolidHhtFloat64History, + ExactNewtonLoopHistoryDetectsThreeStateRecurrence) { + // Three distinct entry states form a closed deterministic Newton cycle. + // The prospective state is neither the current nor its immediate predecessor. + fesa::detail::SolidHhtNewtonHistory history; + fesa::Vector first(3), second(3), third(3); + first[0] = 1; + second[0] = 2; + third[0] = 3; + const std::vector first_mask{true, false}, second_mask{false, true}; + const std::array first_metrics{1e-20, 2e-20, 3e-20, 4e-20}; + const std::array second_metrics{2e-20, 3e-20, 4e-20, 5e-20}; + ASSERT_TRUE(history.Capture(first, first_mask, first_metrics, 4, false)); + ASSERT_TRUE(history.Capture(second, second_mask, second_metrics, 5, false)); + ASSERT_TRUE(history.Capture(third, first_mask, second_metrics, 6, false)); + EXPECT_EQ(history.Size(), 3U); + fesa::SolidDynamicNewtonRepeatEvidence e; + ASSERT_TRUE(history.Match(first, first_mask, first_metrics, false, &e)); + EXPECT_TRUE(e.repeat); + EXPECT_EQ(e.matched_entry_iteration, 4U); + EXPECT_EQ(e.history_entry_count, 3U); + EXPECT_EQ(e.changed_translations, 0U); + EXPECT_EQ(e.changed_holding_mask, 0U); + EXPECT_EQ(e.changed_correction_metrics, 0U); + ASSERT_TRUE(history.Match(first, second_mask, first_metrics, false, &e)); + EXPECT_FALSE(e.repeat); + auto different_metrics = first_metrics; + different_metrics[3] = std::nextafter(first_metrics[3], 1.0); + ASSERT_TRUE(history.Match(first, first_mask, different_metrics, false, &e)); + EXPECT_TRUE( + e.repeat); // Both finite/PASS controls produce the same transition. + EXPECT_EQ(e.changed_correction_metrics, 1U); + auto failing_metrics = first_metrics; + failing_metrics[2] = 1.1e-10; + ASSERT_TRUE(history.Match(first, first_mask, failing_metrics, false, &e)); + EXPECT_FALSE(e.repeat); + first[2] = std::numeric_limits::denorm_min(); + ASSERT_TRUE(history.Match(first, first_mask, first_metrics, false, &e)); + EXPECT_FALSE(e.repeat); + first[2] = 0; + ASSERT_TRUE(history.Match(first, first_mask, first_metrics, true, &e)); + EXPECT_FALSE(e.repeat); // Actual physical equilibrium keeps ordinary Newton. + EXPECT_FALSE( + history.Match(fesa::Vector(2), first_mask, first_metrics, false, &e)); + EXPECT_FALSE(history.Match(first, {false}, first_metrics, false, &e)); + for (double bad : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity()}) { + first[2] = bad; + EXPECT_FALSE(history.Match(first, first_mask, first_metrics, false, &e)); + EXPECT_FALSE(history.Capture(first, first_mask, first_metrics, 7, false)); + first[2] = 0; + auto bad_metrics = first_metrics; + bad_metrics[0] = bad; + EXPECT_FALSE(history.Match(first, first_mask, bad_metrics, false, &e)); + EXPECT_FALSE(history.Capture(first, first_mask, bad_metrics, 7, false)); + } +} + +TEST(SolidHhtFloat64History, + EquivalentFailedNewtonEntrySurvivesRepresentableMove) { + fesa::detail::SolidHhtNewtonHistory history; + fesa::Vector failed(3); + failed[0] = .5; + const std::vector mask{false, true}; + const std::array correction{1e-20, 2e-20, 3e-20, 4e-20}; + ASSERT_TRUE(history.Capture(failed, mask, correction, 23, true)); + EXPECT_EQ(history.Size(), 0U); + auto large = correction; + large[0] = 1.1e-10; + ASSERT_TRUE(history.Capture(failed, mask, large, 23, false)); + EXPECT_EQ(history.Size(), 0U); + ASSERT_TRUE(history.Capture(failed, mask, correction, 23, false)); + fesa::Vector actual = failed; + std::vector applied; + auto first = fesa::detail::TrySolidHhtActualStep( + false, correction, + [&] { + fesa::SolidDynamicNewtonRepeatEvidence e; + EXPECT_TRUE(history.Match(failed, mask, correction, false, &e)); + return fesa::Result::Success(!e.repeat); + }, + [&] { + actual[0] = std::nextafter(actual[0], 0.0); + applied.push_back("accepted-ulp"); + return fesa::Result::Success(true); + }); + ASSERT_TRUE(first.HasValue()); + ASSERT_TRUE(first.Value()); + EXPECT_NE(actual[0], failed[0]); + auto new_metrics = correction; + new_metrics[0] = std::nextafter(new_metrics[0], 1.0); + // The ordinary Newton proposal returns to the exact previously failed state. + auto second = fesa::detail::TrySolidHhtActualStep( + false, new_metrics, + [&] { + applied.push_back("returning-newton-evaluated"); + fesa::SolidDynamicNewtonRepeatEvidence e; + EXPECT_TRUE(history.Match(failed, mask, new_metrics, false, &e)); + EXPECT_TRUE(e.repeat); + EXPECT_EQ(e.matched_entry_iteration, 23U); + return fesa::Result::Success(!e.repeat); + }, + [&] { + actual[0] = std::nextafter(actual[0], 0.0); + applied.push_back("another-accepted-ulp"); + return fesa::Result::Success(true); + }); + ASSERT_TRUE(second.HasValue()); + EXPECT_TRUE(second.Value()); + EXPECT_EQ(applied, (std::vector{"accepted-ulp", + "returning-newton-evaluated", + "another-accepted-ulp"})); + auto fresh = failed; + fresh[2] = std::numeric_limits::denorm_min(); + fesa::SolidDynamicNewtonRepeatEvidence e; + ASSERT_TRUE(history.Match(fresh, mask, new_metrics, false, &e)); + EXPECT_FALSE(e.repeat); + ASSERT_TRUE(history.Match(failed, {true, false}, new_metrics, false, &e)); + EXPECT_FALSE(e.repeat); + ASSERT_TRUE(history.Match(failed, mask, new_metrics, true, &e)); + EXPECT_FALSE(e.repeat); +} + +TEST(SolidHhtDynamics, UnsupportedPairedProductsPreserveEmptyHistory) { + // Preserve the original extreme finite model. Its coefficient products + // lie outside the approved error-free paired multiplication range. + const double q = std::ldexp(1.0, -960); + const double delta = std::ldexp(1.0, -1010); + const double force1 = 4 * q, force2 = 8 * delta; + auto d = OneFree(false, 1, 2); + d.materials = {{"M", 1e-300, 0, {}}}; + d.densities = {{0, 27, {}}}; + d.amplitudes = {{"A", "A", {}, {{0, 0}, {1, force1}, {2, force2}}}}; + d.dynamic_step->loads = {{"N0", 1, 1, {}, 0, "A", {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + ASSERT_TRUE(state.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + const auto status = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteIntoForVerification(domain.Value(), state.Value(), 0); + ASSERT_FALSE(status.IsOk()); + ASSERT_FALSE(status.Diagnostics().empty()); + EXPECT_NE(status.Diagnostics().front().message.find("supported"), + std::string::npos); + EXPECT_TRUE(state.Value().AcceptedFrames().empty()); + EXPECT_FALSE(state.Value().BeginTrialFromCommitted().HasValue()); +} + +TEST(SolidHhtDynamics, FailureShadowRecoversCancellationAndRejectsOverflow) { + double sum = 0, correction = 0, absolute_sum = 0; + for (double term : {1e16, 1.0, -1e16}) + ASSERT_TRUE(fesa::detail::AccumulateSolidHhtForceShadow( + term, 0, &sum, &correction, &absolute_sum)); + EXPECT_EQ(sum, 0); + EXPECT_EQ(sum + correction, 1); // Independent exactly representable sum. + sum = 0; + correction = 0; + absolute_sum = 0; + ASSERT_TRUE(fesa::detail::AccumulateSolidHhtForceShadow( + 1e308, 0, &sum, &correction, &absolute_sum)); + EXPECT_FALSE(fesa::detail::AccumulateSolidHhtForceShadow( + -1e308, 0, &sum, &correction, &absolute_sum)); + EXPECT_EQ(sum, 1e308); // No invalid partial shadow is returned or committed. + EXPECT_EQ(absolute_sum, 1e308); + EXPECT_FALSE(fesa::detail::AccumulateSolidHhtForceShadow( + 0, std::numeric_limits::quiet_NaN(), &sum, &correction, + &absolute_sum)); +} diff --git a/tests/unit/analysis/solid_hht_temporal_verification_test.cpp b/tests/unit/analysis/solid_hht_temporal_verification_test.cpp new file mode 100644 index 0000000..bc9f5ad --- /dev/null +++ b/tests/unit/analysis/solid_hht_temporal_verification_test.cpp @@ -0,0 +1,653 @@ +#include + +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/analysis/solid_hht_dynamic_procedure.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/elements/hex8_element.h" +#include "fesa/solvers/linear/mkl_pardiso_solver.h" +#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" +#include "solid_hht_test_fixture.h" + +namespace { +using M3 = std::array, 3>; +struct AffineMotion { + M3 b{}, velocity{}, acceleration{}, piola{}; +}; +AffineMotion Motion(double time, bool nonlinear) { + const double q = 1 - std::cos(2 * time); + const double dq = 2 * std::sin(2 * time); + const double ddq = 4 * std::cos(2 * time); + AffineMotion m; + constexpr double lambda = 4.8; + constexpr double mu = 4.8; + if (!nonlinear) { + const M3 d{{{.015, .008, -.004}, {.006, -.01, .003}, {.005, .004, .012}}}; + double trace = 0; + for (unsigned i = 0; i < 3; ++i) trace += q * d[i][i]; + for (unsigned i = 0; i < 3; ++i) + for (unsigned j = 0; j < 3; ++j) { + m.b[i][j] = q * d[i][j]; + m.velocity[i][j] = dq * d[i][j]; + m.acceleration[i][j] = ddq * d[i][j]; + m.piola[i][j] = + (i == j ? lambda * trace : 0) + mu * q * (d[i][j] + d[j][i]); + } + } else { + const double s = .004 * q, ds = .004 * dq, dds = .004 * ddq; + const double theta = .5 * q, dt = .5 * dq, ddt = .5 * ddq; + const double c = std::cos(theta), sn = std::sin(theta); + const M3 r{{{c, -sn, 0}, {sn, c, 0}, {0, 0, 1}}}; + const M3 jr{{{-sn, -c, 0}, {c, -sn, 0}, {0, 0, 0}}}; + const M3 jjr{{{-c, sn, 0}, {-sn, -c, 0}, {0, 0, 0}}}; + const double eg = .5 * ((1 + s) * (1 + s) - 1); + const double stress = (3 * lambda + 2 * mu) * eg; + for (unsigned i = 0; i < 3; ++i) + for (unsigned j = 0; j < 3; ++j) { + m.b[i][j] = (1 + s) * r[i][j] - (i == j ? 1 : 0); + m.velocity[i][j] = ds * r[i][j] + (1 + s) * dt * jr[i][j]; + m.acceleration[i][j] = dds * r[i][j] + 2 * ds * dt * jr[i][j] + + (1 + s) * (ddt * jr[i][j] + dt * dt * jjr[i][j]); + m.piola[i][j] = (1 + s) * r[i][j] * stress; + } + } + return m; +} +std::array Apply(const M3& matrix, const std::array& x) { + std::array result{}; + for (unsigned i = 0; i < 3; ++i) + for (unsigned j = 0; j < 3; ++j) result[i] += matrix[i][j] * x[j]; + return result; +} +fesa::ModelDefinition TwoBricks(double h, unsigned count, bool nonlinear) { + auto d = solid_hht_test::Brick(h, count, nonlinear); + const std::array, 4> added{ + {{2, 0, 0}, {2, 1, 0}, {2, 0, 1}, {2, 1, 1}}}; + for (unsigned j = 0; j < 4; ++j) { + unsigned n = 8 + j; + d.nodes.push_back({{"I", n + 1, std::to_string(n + 1)}, added[j], {}}); + d.node_sets.push_back({"N" + std::to_string(n), {}, {n}, {}}); + } + d.solid_elements.emplace_back( + fesa::SourceEntityId{"I", 9, "9"}, + std::array{1, 8, 9, 2, 5, 10, 11, 6}, 0, 0, + fesa::SourceLocation{}); + return d; +} +// Closed polynomial brick integrals: integral Na*X and integral grad Na. +// These do not call the production stiffness, mass, quadrature or response. +std::vector ExactForce(const fesa::ModelDefinition& d, double t, + bool nonlinear) { + const auto m = Motion(t, nonlinear); + std::vector f(3 * d.nodes.size()); + for (unsigned e = 0; e < 2; ++e) { + const auto& nodes = d.solid_elements[e].node_indices; + for (unsigned a = 0; a < 8; ++a) { + const auto& x = d.nodes[nodes[a]].coordinates; + std::array integral_x{}, gradient{}; + for (unsigned c = 0; c < 3; ++c) { + const double origin = c == 0 ? static_cast(e) : 0; + const double local = x[c] - origin; + integral_x[c] = 2 * (origin / 8 + (1 + local) / 24); + gradient[c] = (2 * local - 1) / 4; + } + const auto inertia = Apply(m.acceleration, integral_x); + const auto internal = Apply(m.piola, gradient); + for (unsigned c = 0; c < 3; ++c) + f[3 * nodes[a] + c] += inertia[c] + internal[c]; + } + } + return f; +} +std::vector ExactMassProduct(const fesa::Domain& d, + const fesa::Vector& v) { + std::vector result(v.Size()); + for (unsigned e = 0; e < 2; ++e) { + const auto& nodes = d.SolidElements()[e].node_indices; + for (unsigned a = 0; a < 8; ++a) + for (unsigned b = 0; b < 8; ++b) { + double coefficient = 2.0 / 216; + for (unsigned c = 0; c < 3; ++c) + coefficient *= d.Nodes()[nodes[a]].coordinates[c] == + d.Nodes()[nodes[b]].coordinates[c] + ? 2 + : 1; + for (unsigned c = 0; c < 3; ++c) + result[3 * nodes[a] + c] += coefficient * v[3 * nodes[b] + c]; + } + } + return result; +} +void AddCross(const std::array& x, const std::array& f, + std::array* result) { + (*result)[0] += x[1] * f[2] - x[2] * f[1]; + (*result)[1] += x[2] * f[0] - x[0] * f[2]; + (*result)[2] += x[0] * f[1] - x[1] * f[0]; +} +std::array ManufacturedRun(double h, unsigned count, + bool nonlinear, + bool compare_parallel = false) { + auto d = TwoBricks(h, count, nonlinear); + for (unsigned i = 0; i < 3 * d.nodes.size(); ++i) { + fesa::TabularAmplitudeDefinition amplitude; + amplitude.name_raw = "A" + std::to_string(i); + amplitude.name_canonical = amplitude.name_raw; + for (unsigned n = 0; n <= count; ++n) { + const double t = static_cast(n) * h; + amplitude.points.push_back({t, ExactForce(d, t, nonlinear)[i]}); + } + d.amplitudes.push_back(amplitude); + d.dynamic_step->loads.push_back({"N" + std::to_string(i / 3), + static_cast(i % 3 + 1), + 1, + {}, + i, + amplitude.name_raw, + {}}); + } + auto domain = fesa::Domain::Create(std::move(d)); + EXPECT_TRUE(domain.HasValue()); + if (!domain.HasValue()) return {1, 1, 1}; + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + auto status = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + EXPECT_TRUE(status.IsOk()) << solid_hht_test::Message(status); + if (!status.IsOk()) { + const auto& e = state.Value().IterationRecords().back().evidence; + std::cout << "Failed force metric " << e.residual_metric << " update " + << e.update_metric << '\n'; + return {1, 1, 1}; + } + if (count == 30 || compare_parallel) { + // Both element workers have distinct index-owned buffers; compare every + // accepted field, not merely the final displacement norm. + auto parallel_state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + fesa::TbbParallelFor tbb; + fesa::MklPardisoSolver parallel_mass; + fesa::MklPardisoTangentSolver parallel_tangent; + auto parallel_status = + fesa::SolidHhtDynamicProcedureExecutor(tbb, parallel_mass, + parallel_tangent) + .ExecuteInto(domain.Value(), parallel_state.Value()); + EXPECT_TRUE(parallel_status.IsOk()) + << solid_hht_test::Message(parallel_status); + if (parallel_status.IsOk()) { + const auto& serial_frames = state.Value().AcceptedFrames(); + const auto& parallel_frames = parallel_state.Value().AcceptedFrames(); + EXPECT_EQ(serial_frames.size(), parallel_frames.size()); + bool has_nonzero_low = false; + for (std::size_t n = 0; n < serial_frames.size(); ++n) { + const auto& sf = serial_frames[n]; + const auto& pf = parallel_frames[n]; + EXPECT_DOUBLE_EQ(sf.kinetic_energy, pf.kinetic_energy); + EXPECT_DOUBLE_EQ(sf.strain_energy, pf.strain_energy); + EXPECT_DOUBLE_EQ(sf.external_work, pf.external_work); + EXPECT_EQ(sf.angular_momentum, pf.angular_momentum); + EXPECT_EQ(sf.low.strain_energy, pf.low.strain_energy); + EXPECT_EQ(sf.low.kinetic_energy, pf.low.kinetic_energy); + EXPECT_EQ(sf.low.external_work, pf.low.external_work); + EXPECT_EQ(sf.low.energy_balance_error, pf.low.energy_balance_error); + EXPECT_EQ(sf.low.linear_momentum, pf.low.linear_momentum); + EXPECT_EQ(sf.low.angular_momentum, pf.low.angular_momentum); + EXPECT_EQ(sf.low.linear_momentum_rate, pf.low.linear_momentum_rate); + EXPECT_EQ(sf.low.angular_momentum_rate, pf.low.angular_momentum_rate); + EXPECT_EQ(sf.low.global_equilibrium, pf.low.global_equilibrium); + EXPECT_EQ(sf.low.linearization_current_internal_moment, + pf.low.linearization_current_internal_moment); + EXPECT_EQ(sf.low.reference_internal_moment, + pf.low.reference_internal_moment); + EXPECT_EQ(sf.low.current_free_residual_moment, + pf.low.current_free_residual_moment); + for (auto member : {&fesa::SolidDynamicFrameLowParts::displacement, + &fesa::SolidDynamicFrameLowParts::velocity, + &fesa::SolidDynamicFrameLowParts::acceleration, + &fesa::SolidDynamicFrameLowParts::algorithmic_acceleration, + &fesa::SolidDynamicFrameLowParts::inertia_force, + &fesa::SolidDynamicFrameLowParts::internal_force, + &fesa::SolidDynamicFrameLowParts::external_force, + &fesa::SolidDynamicFrameLowParts::residual, + &fesa::SolidDynamicFrameLowParts::reaction, + &fesa::SolidDynamicFrameLowParts::residual_scale}) { + const auto& sl = sf.low.*member; + const auto& pl = pf.low.*member; + EXPECT_EQ(sl.Size(), pl.Size()); + if (sl.Size() != pl.Size()) continue; + for (std::size_t i = 0; i < sl.Size(); ++i) { + EXPECT_EQ(sl[i], pl[i]); + has_nonzero_low = has_nonzero_low || sl[i] != 0; + } + } + for (std::size_t i = 0; i < sf.displacement.Size(); ++i) { + EXPECT_DOUBLE_EQ(sf.displacement[i], pf.displacement[i]); + EXPECT_DOUBLE_EQ(sf.velocity[i], pf.velocity[i]); + EXPECT_DOUBLE_EQ(sf.acceleration[i], pf.acceleration[i]); + EXPECT_DOUBLE_EQ(sf.algorithmic_acceleration[i], + pf.algorithmic_acceleration[i]); + EXPECT_DOUBLE_EQ(sf.residual[i], pf.residual[i]); + } + for (std::size_t point = 0; point < sf.solid_results.size(); ++point) { + EXPECT_EQ(sf.solid_results[point].cauchy_stress, + pf.solid_results[point].cauchy_stress); + EXPECT_EQ(sf.solid_results[point].strain, + pf.solid_results[point].strain); + } + } + EXPECT_TRUE(has_nonzero_low); + } + } + std::array errors{}, scales{}; + const auto exact = Motion(static_cast(count) * h, nonlinear); + const auto& final = state.Value().CommittedState(); + for (unsigned node = 0; node < domain.Value().Nodes().size(); ++node) { + const auto u = Apply(exact.b, domain.Value().Nodes()[node].coordinates); + const auto v = + Apply(exact.velocity, domain.Value().Nodes()[node].coordinates); + const auto a = + Apply(exact.acceleration, domain.Value().Nodes()[node].coordinates); + for (unsigned c = 0; c < 3; ++c) { + const std::array truth{u[c], v[c], a[c]}; + const std::array actual{final.displacement[3 * node + c], + final.velocity[3 * node + c], + final.acceleration[3 * node + c]}; + for (unsigned j = 0; j < 3; ++j) { + errors[j] += (actual[j] - truth[j]) * (actual[j] - truth[j]); + scales[j] += truth[j] * truth[j]; + } + } + } + for (unsigned j = 0; j < 3; ++j) errors[j] = std::sqrt(errors[j] / scales[j]); + for (const auto& f : state.Value().AcceptedFrames()) { + const auto mv = ExactMassProduct(domain.Value(), f.velocity); + const auto ma = ExactMassProduct(domain.Value(), f.acceleration); + std::array p{}, angular{}, pdot{}, ldot{}, vv{}; + double kinetic = 0, energy_scale = 0, momentum_scale = 0; + double rate_scale = 0, moment_scale = 0, angular_scale = 0, vv_scale = 0; + double free_residual_sum = 0; + for (unsigned n = 0; n < domain.Value().Nodes().size(); ++n) { + auto x = domain.Value().Nodes()[n].coordinates; + std::array nmv{}, nma{}, v{}; + for (unsigned c = 0; c < 3; ++c) { + x[c] += f.displacement[3 * n + c]; + v[c] = f.velocity[3 * n + c]; + nmv[c] = mv[3 * n + c]; + nma[c] = ma[3 * n + c]; + p[c] += nmv[c]; + pdot[c] += nma[c]; + kinetic += .5 * v[c] * nmv[c]; + energy_scale += .5 * std::abs(v[c] * nmv[c]); + momentum_scale += std::abs(nmv[c]); + rate_scale += std::abs(nma[c]) + std::abs(f.internal_force[3 * n + c]) + + std::abs(f.external_force[3 * n + c]) + + std::abs(f.reaction[3 * n + c]); + free_residual_sum += + std::abs(f.residual[3 * n + c] - f.reaction[3 * n + c]); + } + double arm = 0, velocity_scale = 0, node_mv = 0, node_rate = 0; + for (unsigned c = 0; c < 3; ++c) { + arm += + std::abs(x[c]) + std::abs(domain.Value().Nodes()[n].coordinates[c]); + velocity_scale += std::abs(v[c]); + node_mv += std::abs(nmv[c]); + node_rate += std::abs(nma[c]) + std::abs(f.internal_force[3 * n + c]) + + std::abs(f.external_force[3 * n + c]) + + std::abs(f.reaction[3 * n + c]); + } + angular_scale += arm * node_mv; + moment_scale += arm * node_rate; + vv_scale += velocity_scale * node_mv; + AddCross(x, nmv, &angular); + AddCross(x, nma, &ldot); + AddCross(v, nmv, &vv); + } + const double roundoff = 512 * std::numeric_limits::epsilon(); + EXPECT_NEAR(f.kinetic_energy, kinetic, roundoff * energy_scale); + for (unsigned c = 0; c < 3; ++c) { + EXPECT_NEAR(f.linear_momentum[c], p[c], roundoff * momentum_scale); + EXPECT_NEAR(f.angular_momentum[c], angular[c], roundoff * angular_scale); + EXPECT_NEAR(f.linear_momentum_rate[c], pdot[c], roundoff * rate_scale); + EXPECT_NEAR(f.angular_momentum_rate[c], ldot[c], roundoff * moment_scale); + EXPECT_NEAR(vv[c], 0, roundoff * vv_scale); + EXPECT_NEAR(f.global_equilibrium[c], 0, + free_residual_sum + roundoff * rate_scale); + const double identity = f.reference_internal_moment[c] + + f.linearization_current_internal_moment[c] - + f.current_free_residual_moment[c]; + EXPECT_NEAR(f.global_equilibrium[3 + c], identity, + roundoff * moment_scale); + if (nonlinear) + EXPECT_NEAR(f.reference_internal_moment[c] + + f.linearization_current_internal_moment[c], + 0, roundoff * moment_scale); + else + EXPECT_NEAR(f.reference_internal_moment[c], 0, roundoff * moment_scale); + } + } + return errors; +} +} // namespace + +TEST(SolidHhtTemporal, PairedSerialAndTbbPreserveEveryLowComponent) { + for (bool nonlinear : {false, true}) { + SCOPED_TRACE(nonlinear); + ManufacturedRun(.02, 2, nonlinear, true); + } +} + +TEST(SolidHhtTemporal, FreeConstantForceHasAnalyticalInitialAcceleration) { + auto d = solid_hht_test::Brick(); + d.dynamic_step->boundaries = {{"ALL", 2, 3, 0, {}}}; + d.dynamic_step->loads = {{"ALL", 1, 0.25, {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto result = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus()); + const auto& initial = result.Value().state.AcceptedFrames().front(); + // Brick mass row sum rho V/8 = 1/4; rigid translation has zero stiffness. + for (unsigned node = 0; node < 8; ++node) { + EXPECT_NEAR(initial.acceleration[3 * node], 1, 1e-13); + EXPECT_DOUBLE_EQ(initial.displacement[3 * node], 0); + } +} + +TEST(SolidHhtTemporal, + AllFreeUniformForceRetainsExactRigidTranslationBehavior) { + auto d = solid_hht_test::Brick(); + d.dynamic_step->loads = {{"ALL", 1, .25, {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto model = fesa::AnalysisModel::Create(domain.Value()); + auto dofs = fesa::DofManager::Create(model.Value()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + auto status = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + if (!status.IsOk()) { + const auto& e = state.Value().IterationRecords().back().evidence; + std::cout << "Rigid-force failure worst DOF " << e.worst_full_dof + << " force ratio " << e.residual_metric + << " terms r,Ma,f,F,scale "; + for (auto term : e.worst_force_terms) std::cout << term << ' '; + std::cout << '\n'; + for (const auto& rec : state.Value().IterationRecords()) { + if (rec.evidence.iteration_index < 10) { + std::cout << "Trace " << rec.increment_index << ' ' + << rec.evidence.iteration_index << ' ' << rec.outcome + << " ratio=" << rec.evidence.residual_metric + << " worst=" << rec.evidence.worst_full_dof + << " correction=" << rec.evidence.correction_metrics[0] + << '\n'; + } + } + } + ASSERT_TRUE(status.IsOk()) << solid_hht_test::Message(status); + EXPECT_TRUE(std::any_of( + state.Value().IterationRecords().begin(), + state.Value().IterationRecords().end(), + [](const auto& record) { return record.outcome == "paired-candidate"; })); + const double roundoff = 512 * std::numeric_limits::epsilon(); + for (const auto& f : state.Value().AcceptedFrames()) + for (unsigned n = 0; n < 8; ++n) { + EXPECT_NEAR(f.displacement[3 * n], .5 * f.time * f.time, + roundoff * .5 * f.time * f.time); + EXPECT_NEAR(f.velocity[3 * n], f.time, roundoff * f.time); + EXPECT_NEAR(f.acceleration[3 * n], 1, roundoff); + for (unsigned c = 1; c < 3; ++c) { + EXPECT_NEAR(f.displacement[3 * n + c], 0, + roundoff * .5 * f.time * f.time); + EXPECT_NEAR(f.velocity[3 * n + c], 0, roundoff * f.time); + EXPECT_NEAR(f.acceleration[3 * n + c], 0, roundoff); + } + } +} + +TEST(SolidHhtTemporal, + PairedNewtonKeepsAllOriginalComponentsAcrossAxesAndScales) { + const double roundoff = 512 * std::numeric_limits::epsilon(); + for (unsigned axis = 0; axis < 3; ++axis) + for (double acceleration : {1e-4, 1.0, 1000.0}) { + SCOPED_TRACE(axis); + SCOPED_TRACE(acceleration); + auto d = solid_hht_test::Brick(.01, 8); + d.dynamic_step->loads = { + {"ALL", static_cast(axis + 1), .25 * acceleration, {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + auto state = + fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value()); + ASSERT_TRUE(state.HasValue()); + const auto status = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .ExecuteInto(domain.Value(), state.Value()); + if (!status.IsOk()) { + const auto& records = state.Value().IterationRecords(); + for (const auto& record : records) { + if (record.increment_index != records.back().increment_index || + (record.evidence.iteration_index > 8 && + record.evidence.iteration_index < 95)) + continue; + std::cout << "Rigid paired Newton trace " + << record.evidence.iteration_index << ' ' << record.outcome + << " metric=" << record.evidence.residual_metric + << " worst=" << record.evidence.worst_full_dof << " terms="; + for (double term : record.evidence.worst_force_terms) + std::cout << term << ' '; + std::cout << " corrections="; + for (double metric : record.evidence.correction_metrics) + std::cout << metric << ' '; + std::cout << '\n'; + } + } + ASSERT_TRUE(status.IsOk()) << solid_hht_test::Message(status); + for (const auto& frame : state.Value().AcceptedFrames()) { + EXPECT_LE(frame.iterations.back().residual_metric, 1e-10); + EXPECT_LE(frame.iterations.back().update_metric, 1e-12); + for (unsigned node = 0; node < 8; ++node) + for (unsigned c = 0; c < 3; ++c) { + const double a = c == axis ? acceleration : 0; + EXPECT_NEAR(frame.displacement[3 * node + c], + .5 * a * frame.time * frame.time, + roundoff * .5 * acceleration * frame.time * frame.time); + EXPECT_NEAR(frame.velocity[3 * node + c], a * frame.time, + roundoff * acceleration * frame.time); + EXPECT_NEAR(frame.acceleration[3 * node + c], a, + roundoff * acceleration); + } + } + } +} + +TEST(SolidHhtTemporal, + ActualTwoBrickLinearAndRotatingSmallStrainTlAreSecondOrder) { + for (bool nonlinear : {false, true}) { + SCOPED_TRACE(nonlinear); + const auto coarse = ManufacturedRun(.04, 30, nonlinear); + const auto medium = ManufacturedRun(.02, 60, nonlinear); + const auto fine = ManufacturedRun(.01, 120, nonlinear); + for (unsigned j = 0; j < 3; ++j) { + const double p1 = std::log(coarse[j] / medium[j]) / std::log(2.0); + const double p2 = std::log(medium[j] / fine[j]) / std::log(2.0); + std::cout << "Hex8 " << (nonlinear ? "TL" : "linear") << " u/v/a channel " + << j << " errors " << coarse[j] << ' ' << medium[j] << ' ' + << fine[j] << " orders " << p1 << ' ' << p2 << '\n'; + EXPECT_GT(fine[j], + 1e-8); // Meaningful truncation interval, away from noise. + EXPECT_LT(coarse[j], .1); + EXPECT_GE(p1, 1.8); + EXPECT_GE(p2, 1.8); + } + } +} + +TEST(SolidHhtTemporal, IndependentThreeDofBrickOscillatorHasSecondOrderUVA) { + // Node0's three translations are free. Exact polynomial integrals yield + // M=(rho/27)I, Kii=(lambda+4mu)/9, Kij=(lambda+mu)/12. + std::array, 3> errors{}; + for (unsigned level = 0; level < 3; ++level) { + const double h = .02 / static_cast(1U << level); + const unsigned count = 32U << level; + auto d = solid_hht_test::Brick(h, count); + for (unsigned n = 1; n < 8; ++n) + d.dynamic_step->boundaries.push_back( + {"N" + std::to_string(n), 1, 3, 0, {}}); + d.dynamic_step->loads = { + {"N0", 1, .01, {}}, {"N0", 2, .02, {}}, {"N0", 3, .03, {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto result = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) + << solid_hht_test::Message(result.GetStatus()); + const auto& f = result.Value().state.CommittedState(); + const double m = 2.0 / 27; + const double diagonal = (4.8 + 4 * 4.8) / 9; + const double off = (4.8 + 4.8) / 12; + const std::array k{diagonal + 2 * off, diagonal - off}; + for (unsigned c = 0; c < 3; ++c) { + std::array truth{}; + const std::array force{.02, .01 * (c + 1) - .02}; + for (unsigned mode = 0; mode < 2; ++mode) { + const double w = std::sqrt(k[mode] / m); + truth[0] += force[mode] / k[mode] * (1 - std::cos(w * f.time)); + truth[1] += force[mode] / (m * w) * std::sin(w * f.time); + truth[2] += force[mode] / m * std::cos(w * f.time); + } + const std::array actual{f.displacement[c], f.velocity[c], + f.acceleration[c]}; + for (unsigned j = 0; j < 3; ++j) + errors[level][j] += std::pow(actual[j] - truth[j], 2); + } + } + for (unsigned j = 0; j < 3; ++j) { + const double p1 = + .5 * std::log(errors[0][j] / errors[1][j]) / std::log(2.0); + const double p2 = + .5 * std::log(errors[1][j] / errors[2][j]) / std::log(2.0); + std::cout << "Independent MDOF channel " << j << " orders " << p1 << ' ' + << p2 << '\n'; + EXPECT_GE(p1, 1.8); + EXPECT_GE(p2, 1.8); + EXPECT_GT(errors[2][j], 1e-15); + } +} + +TEST(SolidHhtTemporal, IndependentScalarBrickOscillatorHasSecondOrderUVA) { + std::array, 3> errors{}; + const double m = 2.0 / 27; + const double k = (4.8 + 4 * 4.8) / 9; + const double w = std::sqrt(k / m); + for (unsigned level = 0; level < 3; ++level) { + const double h = .02 / static_cast(1U << level); + auto d = solid_hht_test::Brick(h, 32U << level); + d.dynamic_step->boundaries.push_back({"N0", 2, 3, 0, {}}); + for (unsigned n = 1; n < 8; ++n) + d.dynamic_step->boundaries.push_back( + {"N" + std::to_string(n), 1, 3, 0, {}}); + d.dynamic_step->loads = {{"N0", 1, .01, {}}}; + auto domain = fesa::Domain::Create(std::move(d)); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + auto result = + fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent) + .Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) + << solid_hht_test::Message(result.GetStatus()); + const auto& f = result.Value().state.CommittedState(); + const std::array truth{.01 / k * (1 - std::cos(w * f.time)), + .01 / (m * w) * std::sin(w * f.time), + .01 / m * std::cos(w * f.time)}; + const std::array actual{f.displacement[0], f.velocity[0], + f.acceleration[0]}; + for (unsigned j = 0; j < 3; ++j) + errors[level][j] = std::abs(actual[j] - truth[j]); + EXPECT_NEAR( + f.kinetic_energy, .5 * m * f.velocity[0] * f.velocity[0], + 512 * std::numeric_limits::epsilon() * f.kinetic_energy); + EXPECT_NEAR(f.strain_energy, .5 * k * f.displacement[0] * f.displacement[0], + 512 * std::numeric_limits::epsilon() * f.strain_energy); + } + for (unsigned j = 0; j < 3; ++j) { + const double p1 = std::log(errors[0][j] / errors[1][j]) / std::log(2.0); + const double p2 = std::log(errors[1][j] / errors[2][j]) / std::log(2.0); + std::cout << "Independent SDOF channel " << j << " orders " << p1 << ' ' + << p2 << '\n'; + EXPECT_GE(p1, 1.8); + EXPECT_GE(p2, 1.8); + EXPECT_GT(errors[2][j], 1e-8); + } +} + +TEST(SolidHhtTemporal, + BendingMeshRefinementQuantifiesFullIntegrationShearLocking) { + constexpr double thickness = .125; + constexpr double curvature = .01; + constexpr double lambda = 4.8; + constexpr double mu = 4.8; + const double continuum = .5 * (lambda + 2 * mu) * curvature * curvature * + thickness * thickness * thickness / 12; + double previous_error = 0; + for (unsigned nx : {1U, 2U, 4U, 8U}) { + double numerical = 0; + for (unsigned e = 0; e < nx; ++e) { + auto d = solid_hht_test::Brick(); + fesa::Hex8Element::Coordinates x{}; + fesa::ElementDofLayout layout{ + {"I", e + 1, std::to_string(e + 1)}, + {0, 1, 2, 3, 4, 5, 6, 7}, + {fesa::DofComponent::kUx, fesa::DofComponent::kUy, + fesa::DofComponent::kUz}}; + fesa::Vector u(24); + for (unsigned a = 0; a < 8; ++a) { + x[a] = {static_cast(e) / nx + d.nodes[a].coordinates[0] / nx, + thickness * (d.nodes[a].coordinates[1] - .5), + d.nodes[a].coordinates[2]}; + u[3 * a] = -curvature * x[a][0] * x[a][1]; + u[3 * a + 1] = .5 * curvature * x[a][0] * x[a][0]; + } + auto element = fesa::Hex8Element::Create(e, layout, x, d.materials[0], 2); + ASSERT_TRUE(element.HasValue()); + auto response = element.Value().EvaluateLinear(u); + ASSERT_TRUE(response.HasValue()); + numerical += response.Value().strain_energy; + } + const double error = numerical / continuum - 1; + const double exact_parasitic = + mu / (lambda + 2 * mu) / std::pow(nx * thickness, 2); + EXPECT_NEAR(error, exact_parasitic, 2e-11); + if (previous_error != 0) EXPECT_NEAR(previous_error / error, 4, 1e-10); + std::cout << "Bending nx=" << nx + << " energy/continuum=" << numerical / continuum << '\n'; + previous_error = error; + } +} diff --git a/tests/unit/analysis/solid_hht_test_fixture.h b/tests/unit/analysis/solid_hht_test_fixture.h new file mode 100644 index 0000000..d7b13cf --- /dev/null +++ b/tests/unit/analysis/solid_hht_test_fixture.h @@ -0,0 +1,52 @@ +#ifndef FESA_TESTS_SOLID_HHT_TEST_FIXTURE_H_ +#define FESA_TESTS_SOLID_HHT_TEST_FIXTURE_H_ + +#include +#include + +#include "fesa/elements/hex8_definition.h" +#include "fesa/model/domain.h" +#include "fesa/properties/solid_section.h" + +namespace solid_hht_test { +inline fesa::ModelDefinition Brick(double h = 0.01, unsigned count = 20, + bool nonlinear = false) { + fesa::ModelDefinition d; + d.source_path = "solid-hht-test.inp"; + const std::array, 8> x{{{0, 0, 0}, + {1, 0, 0}, + {1, 1, 0}, + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, + {1, 1, 1}, + {0, 1, 1}}}; + for (unsigned i = 0; i < 8; ++i) { + d.nodes.push_back({{"I", i + 1, std::to_string(i + 1)}, x[i], {}}); + d.node_sets.push_back({"N" + std::to_string(i), {}, {i}, {}}); + } + d.node_sets.push_back({"ALL", {}, {0, 1, 2, 3, 4, 5, 6, 7}, {}}); + d.materials = {{"M", 12.0, 0.25, {}}}; + d.densities = {{0, 2.0, {}}}; + d.solid_sections.emplace_back(fesa::SourceEntityId{"I", 1, "1"}, "E", 0, + fesa::SourceLocation{}); + d.solid_elements.emplace_back( + fesa::SourceEntityId{"I", 8, "8"}, + std::array{0, 1, 2, 3, 4, 5, 6, 7}, 0, 0, + fesa::SourceLocation{}); + fesa::DynamicStepDefinition step; + step.source_step_name_raw = "Dynamic"; + step.source_step_name_canonical = "DYNAMIC"; + step.initial_increment = h; + step.time_period = static_cast(count) * h; + step.increment_count = count; + step.maximum_increments = count; + step.nlgeom = nonlinear; + d.dynamic_step = step; + return d; +} +inline std::string Message(const fesa::Status& s) { + return s.Diagnostics().empty() ? "" : s.Diagnostics()[0].message; +} +} // namespace solid_hht_test +#endif