From dd54b0ae1e68a313ac46933b886c865c85ec4432 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Mon, 19 Oct 2026 02:49:53 +0900 Subject: [PATCH] Connect HHT field convergence to production acceptance and v2 logs --- docs/contracts.md | 24 ++ include/fesa/analysis/newton_iteration_log.h | 22 +- src/fesa/analysis/hht_dynamic_assembly.cpp | 36 +- src/fesa/analysis/hht_dynamic_internal.h | 12 +- src/fesa/analysis/hht_reduced_procedure.cpp | 110 +++++- src/fesa/analysis/newton_iteration_log.cpp | 90 +++++ .../analysis/hht_dynamic_algebra_test.cpp | 35 ++ .../analysis/hht_dynamic_procedure_test.cpp | 359 ++++++++++++++++++ .../analysis/newton_iteration_log_test.cpp | 146 +++++++ 9 files changed, 811 insertions(+), 23 deletions(-) diff --git a/docs/contracts.md b/docs/contracts.md index b8b2a0c..868b542 100644 --- a/docs/contracts.md +++ b/docs/contracts.md @@ -105,3 +105,27 @@ Procedure는 family별 completed outcome을 반환하며 로그 close와 최종 Hex8/HHT Python 비교기는 입력·reference·HDF5를 읽고 지정한 JSON report만 쓴다. Solver를 기동하거나 원본을 수정하지 않는다. 입력 부재·hash 충돌·schema/identity/finite 불일치는 실패이며 합성 fixture 통과는 실제 모델 성공을 뜻하지 않는다. + +## Shell HHT iteration 로그 독자 + +Frame0과 직접 v1 writer는 기존 `hht-log-v1` 의미를 유지한다. 양의 reduced-coordinate +MITC4 HHT 증분은 `hht-log-v2`, `policy_id=abaqus2025-field-on-fesa-v1`을 기록하므로 +한 파일 안의 각 record에서 schema를 읽는다. v1 숫자를 v2 field 값으로 재해석하지 않는다. +V2의 고정 키 순서·분기·단위·평균·threshold는 +[field 수렴 정식화](abaqus-tolerance/formulation-convergence.md#로그-소비-계약)를 따른다. +`strict_` 진단은 기존 힘/모멘트 residual과 무차원 q/v/a/A correction이며 종료 기준이 아니다. + +UTF-8 LF 한 줄과 tab 구분 key=value를 사용한다. 첫 `=`만 구분자이며 값의 backslash, +tab, LF, CR은 각각 `\\`, `\t`, `\n`, `\r`로 escape한다. 미평가 numeric/bool은 +literal `not-evaluated`, 평가 bool은 `0`/`1`이다. 독자는 중복·필수키 누락·알 수 없는 +branch/policy/schema·NaN/Inf를 오류로 처리한다. 숫자는 locale-independent binary64 +round-trip 표현이며 step은 raw source 이름이다. 실패 전 평가된 부분값만 보존한다. +`field_pass=1`은 후속 acceptance 통과를 뜻하지 않는다. `accepted=1`도 최종 HDF5 성공을 +뜻하지 않으며 CLI exit와 완결 결과 identity를 함께 확인한다. 실패 증분은 accepted +frame/private/flux history에 반영되지 않는다. + +현재 acceptance transaction은 기존 accepted frame 전체를 가진 state를 지역 staged copy로 +검증·commit한 뒤 로그 write/flush 성공 시 move로 공개한다. private/flux history도 같은 +성공 경계에서만 교체한다. 이 최소 구현은 증분마다 누적 frame 수에 비례하는 copy 시간과 +추가 메모리를 사용하므로 전체 실행의 복사 비용은 frame 수에 대해 이차로 증가할 수 있다. +후보의 실제 elapsed 측정에는 이 비용을 포함하며 출력 실패 안전성을 위해 이를 생략하지 않는다. diff --git a/include/fesa/analysis/newton_iteration_log.h b/include/fesa/analysis/newton_iteration_log.h index c59ab27..6ce4a7b 100644 --- a/include/fesa/analysis/newton_iteration_log.h +++ b/include/fesa/analysis/newton_iteration_log.h @@ -5,6 +5,7 @@ #include #include #include +#include #include #include "fesa/analysis/nonlinear_static_control.h" @@ -52,6 +53,17 @@ struct NewtonIterationLogRecord { bool combined_passed; }; +/// @brief Vendor-neutral field evidence in full global u/total-phi coordinates. +/// @note Missing quantities serialize as literal not-evaluated. Force uses +/// force/length units; generalized moment uses moment/rad units. +struct DynamicFieldLogEvidence { + std::optional residual_numerator, residual_scale, spatial_flux; + std::optional time_average_flux, correction_numerator, increment_norm; + std::optional estimated_correction, residual_threshold, correction_threshold; + std::string branch = "not-evaluated"; + std::optional residual_pass, correction_pass, field_pass; +}; + /// @brief Carries two-channel HHT endpoint convergence evidence. /// @note Array entries are translation and rotation channels, respectively. /// Corrections are dimensionless N2 metrics; rate corrections already include @@ -71,6 +83,14 @@ struct DynamicNewtonIterationLogRecord { std::array algorithmic_correction{}; bool residual_pass = false; bool correction_pass = false; + // Absent policy preserves the historical v1 direct-writer/frame-zero contract. + std::string policy_id; + std::string step_name; + std::array fields; + std::optional field_pass; + bool accepted = false; + bool strict_residual_evaluated = false; + bool strict_correction_evaluated = false; }; /// @brief Writes deterministic, line-oriented Newton iteration diagnostics. @@ -90,7 +110,7 @@ class NewtonIterationLogWriter { /// failure with code `log-output-failure`. Status Write(const NewtonIterationLogRecord& record); - /// @brief Writes one complete HHT record using schema `hht-log-v1`. + /// @brief Writes HHT v1 evidence, or v2 when the field policy is identified. /// @return Success after the complete line is written; invalid numeric /// evidence or stream failure permanently fails this writer instance. Status Write(const DynamicNewtonIterationLogRecord& record); diff --git a/src/fesa/analysis/hht_dynamic_assembly.cpp b/src/fesa/analysis/hht_dynamic_assembly.cpp index bec9a3f..f4a5e1a 100644 --- a/src/fesa/analysis/hht_dynamic_assembly.cpp +++ b/src/fesa/analysis/hht_dynamic_assembly.cpp @@ -191,8 +191,34 @@ Result ReferenceLength(const Domain& domain) { return Result::Success(length); } +Result> MeasureEndpointFields( + const DofManager& dofs, const DynamicAcceptedFrameState& trial, + const DynamicAcceptedFrameState& accepted, const std::optional& previous) { + using Fields = std::array; + const auto count = dofs.FullDofCount(); + if (trial.generalized_coordinates.Size() != count || trial.residual.Size() != count || + accepted.generalized_coordinates.Size() != count || (previous && previous->Size() != count)) + return Result::Failure(Failure("Invalid full endpoint field inventory.")); + Fields fields; + for (auto& field : fields) if (previous) field.correction = 0.0; + for (const auto full : dofs.FreeDofs()) { + auto& field = fields[full % 6U < 3U ? 0U : 1U]; + const double residual = std::abs(trial.residual[full]); + const double increment = std::abs(trial.generalized_coordinates[full] - accepted.generalized_coordinates[full]); + const double correction = previous ? std::abs(trial.generalized_coordinates[full] - (*previous)[full]) : 0.0; + if (!std::isfinite(residual) || !std::isfinite(increment) || !std::isfinite(correction)) + return Result::Failure(Failure("Nonfinite full endpoint field arithmetic.")); + field.active = true; + field.residual = (std::max)(field.residual, residual); + field.increment = (std::max)(field.increment, increment); + if (previous) field.correction = (std::max)(*field.correction, correction); + } + return Result::Success(std::move(fields)); +} + Result> AssembleLoads(const Domain& domain, - const DofManager& dofs) { + const DofManager& dofs, + std::array, 2>* field_flux) { const auto index = SourceTargetIndex::FromDomain(domain); const SourceTargetResolver resolver(index); const LoadContext context{domain, dofs, resolver}; @@ -203,6 +229,9 @@ Result> AssembleLoads(const Domain& domain, if (!contributions.HasValue()) return Result>::Failure(contributions.GetStatus()); for (const auto& entry : contributions.Value()) { + if (field_flux != nullptr) + (*field_flux)[entry.full_dof_index % 6U < 3U ? 0U : 1U].push_back( + {entry.value, 0.0}); result[0U][entry.full_dof_index] += entry.value; result[1U][entry.full_dof_index] += std::abs(entry.value); if (!std::isfinite(result[0U][entry.full_dof_index]) || @@ -489,6 +518,7 @@ Result EvaluateAssembly( Vector third(count); Vector third_scale(count); std::array, 6> entries; + std::array, 2> field_flux; for (auto& values : entries) values.reserve(shells.size() * 24U * 24U); double momentum_scale = 0.0; double angular_scale = 0.0; @@ -546,6 +576,8 @@ Result EvaluateAssembly( const Vector mass_action = EntryProducts(mass_values, acceleration_values); for (std::size_t row = 0U; row < 24U; ++row) { const auto global = shells[e].scatter[row]; + field_flux[row % 6U < 3U ? 0U : 1U].push_back( + {f.internal_force[row], m.force[row]}); frame.internal_force[global] += f.internal_force[row]; frame.inertia_force[global] += m.force[row]; third[global] += f.force_directional_second[row]; @@ -708,7 +740,7 @@ Result EvaluateAssembly( std::move(third), std::move(third_scale), std::move(rotations), momentum_scale, angular_scale, std::move(hidden_second), std::move(hidden_second_scale), std::move(hidden_gradient), - std::move(trial_constraints)}); + std::move(trial_constraints), std::move(field_flux)}); } Status CompleteEquilibrium(const Domain& domain, const DofManager& dofs, diff --git a/src/fesa/analysis/hht_dynamic_internal.h b/src/fesa/analysis/hht_dynamic_internal.h index 2911f41..e18dd62 100644 --- a/src/fesa/analysis/hht_dynamic_internal.h +++ b/src/fesa/analysis/hht_dynamic_internal.h @@ -16,6 +16,7 @@ #include "fesa/math/sparse_matrix.h" #include "elements/nonlinear_mitc4_configuration_internal.h" #include "math/dense_blas_internal.h" +#include "analysis/hht_field_convergence.h" namespace fesa { class LinearSolver; @@ -95,6 +96,7 @@ struct DynamicAssembly { // Canonical scalar drilling-constraint derivative, indexed by global node. std::optional hidden_gradient; std::optional trial_constraints; + std::array, 2> field_flux; }; /// @brief Builds checked shell evaluators in stable source order. @@ -104,9 +106,17 @@ Result> BuildRuntimeShells(const Domain& domain, /// @brief Returns the maximum initial active-element edge length. Result ReferenceLength(const Domain& domain); +/// @brief Measures free full endpoint residuals and actual u/total-phi changes. +/// @note Previous is the selected trial in this increment, not a tangent lift. +Result> MeasureEndpointFields( + const DofManager& dofs, const DynamicAcceptedFrameState& trial, + const DynamicAcceptedFrameState& accepted, + const std::optional& previous); + /// @brief Assembles dead loads and their uncancelled absolute row scales. Result> AssembleLoads(const Domain& domain, - const DofManager& dofs); + const DofManager& dofs, + std::array, 2>* field_flux = nullptr); /// @brief Evaluates current kinetic and internal quantities without commit. Result EvaluateAssembly( diff --git a/src/fesa/analysis/hht_reduced_procedure.cpp b/src/fesa/analysis/hht_reduced_procedure.cpp index 408b8e4..b109d1c 100644 --- a/src/fesa/analysis/hht_reduced_procedure.cpp +++ b/src/fesa/analysis/hht_reduced_procedure.cpp @@ -150,7 +150,7 @@ void ResidualRecord(const DofManager& dofs, frame.free_residual_metrics[1U] <= 1e-10; } -void CorrectionRecord(const Vector& coordinate, const Vector& velocity, +bool CorrectionRecord(const Vector& coordinate, const Vector& velocity, const Vector& acceleration, const Vector& full_algorithmic, const Vector& reduced_algorithmic, @@ -189,12 +189,14 @@ void CorrectionRecord(const Vector& coordinate, const Vector& velocity, if (reduced_algorithmic.Size() != chart.p_columns.size() || reduced_physical.Size() != chart.p_columns.size()) { record->correction_pass = false; - return; + return false; } + bool private_pass = true; for (std::size_t i = 0U; i < chart.p_dofs.size(); ++i) - record->correction_pass = record->correction_pass && - PrivateCorrectionEqual(chart, i, full_algorithmic, - reduced_algorithmic, reduced_physical, alpha); + private_pass = private_pass && PrivateCorrectionEqual( + chart, i, full_algorithmic, reduced_algorithmic, reduced_physical, alpha); + record->correction_pass = record->correction_pass && private_pass; + return private_pass; } bool ExactlyEqual(const Vector& left, const Vector& right) { @@ -433,7 +435,8 @@ Result ExecuteReducedHht( auto length_result = ReferenceLength(domain); if (!length_result.HasValue()) return Result::Failure(length_result.GetStatus()); const double length = length_result.Value(); - auto loads = AssembleLoads(domain, dofs); + std::array, 2> load_flux; + auto loads = AssembleLoads(domain, dofs, &load_flux); if (!loads.HasValue()) return Result::Failure(loads.GetStatus()); auto parameters_result = MakeHhtParameters(domain.Dynamic()->alpha); if (!parameters_result.HasValue()) return Result::Failure(parameters_result.GetStatus()); @@ -699,9 +702,11 @@ Result ExecuteReducedHht( status = log.Write(record); if (status.IsOk()) status = log.Flush(); if (!status.IsOk()) return Result::Failure(status); + std::array flux_history{}; struct ReducedStepCandidate { DynamicAcceptedFrameState frame; ReducedHistory history; + std::array flux_history; }; for (std::uint32_t increment = 1U; increment <= domain.Dynamic()->increment_count; ++increment) { @@ -711,9 +716,14 @@ Result ExecuteReducedHht( domain, failure, stage, increment, record.time); record.stage = stage; record.outcome = "failed"; - record.failure_code = contextual.Diagnostics().empty() - ? "dynamic-solver-failure" - : contextual.Diagnostics().front().code; + record.accepted = false; + const auto cause = std::find_if( + failure.Diagnostics().begin(), failure.Diagnostics().end(), + [](const Diagnostic& diagnostic) { + return diagnostic.code != "dynamic-procedure-context"; + }); + record.failure_code = cause == failure.Diagnostics().end() + ? "dynamic-solver-failure" : cause->code; const Status written = log.Write(record); return written.IsOk() ? contextual : written; }; @@ -744,6 +754,10 @@ Result ExecuteReducedHht( record.increment_index = increment; record.time = trial_frame.time; record.correction_pass = true; + record.policy_id = hht_field_convergence::kPolicyId; + record.step_name = domain.Dynamic()->source_step_name_raw; + std::optional previous_selected_q; + std::array, 2> previous_residual, older_residual; bool translation_private_algorithmic_pass = true; bool copy_translation_rates = false; Vector retained_velocity{history.p.Size()}; @@ -755,6 +769,9 @@ Result ExecuteReducedHht( }; ReducedShellHhtStep iteration_step( [&](std::uint32_t iteration) -> StepResult { + record.fields = {}; + record.field_pass.reset(); + record.strict_residual_evaluated = false; record.iteration_index = iteration; record.stage = "evaluation"; record.outcome = "candidate"; @@ -896,7 +913,53 @@ Result ExecuteReducedHht( if (!operators.HasValue()) return fail(operators.GetStatus()); ResidualRecord(dofs, candidate, &record); - if (record.residual_pass && record.correction_pass) { + record.strict_residual_evaluated = true; + auto measured_fields = MeasureEndpointFields(dofs, candidate, old, previous_selected_q); + if (!measured_fields.HasValue()) return fail(measured_fields.GetStatus()); + auto field_inputs = std::move(measured_fields.Value()); + for (std::size_t c = 0; c < 2; ++c) { + field_inputs[c].flux = trial.Value().assembly.field_flux[c]; + field_inputs[c].flux.insert(field_inputs[c].flux.end(), + load_flux[c].begin(), load_flux[c].end()); + field_inputs[c].previous_residual = previous_residual[c]; + field_inputs[c].older_residual = older_residual[c]; + } + const auto field_evaluation = hht_field_convergence::Evaluate( + field_inputs, flux_history, iteration, length); + for (std::size_t c = 0; c < 2; ++c) { + auto& e = record.fields[c]; + const auto& in = field_inputs[c]; + const auto& out = field_evaluation.fields[c]; + if (in.active) { + e.residual_numerator = in.residual; + e.increment_norm = in.increment; + e.correction_numerator = in.correction; + } + e.residual_scale = out.residual_scale; + e.spatial_flux = out.spatial_flux; + e.time_average_flux = out.time_average_flux; + e.estimated_correction = out.estimated_correction; + e.residual_threshold = out.residual_threshold; + e.correction_threshold = out.correction_threshold; + e.residual_pass = out.residual_pass; + e.correction_pass = out.correction_pass; + e.field_pass = out.field_pass; + using Branch = hht_field_convergence::Branch; + switch (out.branch) { + case Branch::kNotEvaluated: e.branch = "not-evaluated"; break; + case Branch::kLinear: e.branch = "linear"; break; + case Branch::kZeroFlux: e.branch = "zero-flux"; break; + case Branch::kNonquadratic: e.branch = "nonquadratic"; break; + case Branch::kNormal: e.branch = "normal"; break; + } + } + if (!field_evaluation.status.IsOk()) return fail(field_evaluation.status); + record.field_pass = field_evaluation.field_pass; + previous_selected_q = candidate.generalized_coordinates; + older_residual = previous_residual; + for (std::size_t c = 0; c < 2; ++c) + previous_residual[c] = field_inputs[c].residual; + if (field_evaluation.field_pass) { status = CompleteConstraintEvidence( domain, dofs, trial.Value().assembly, length, null_nodes, positive_solver, tangent_solver, &candidate); @@ -932,7 +995,7 @@ Result ExecuteReducedHht( return StepResult::Success(ReducedStepCandidate{ std::move(candidate), {trial_p, trial_v, trial_a, trial_algorithmic, - trial.Value().root.z}}); + trial.Value().root.z}, field_evaluation.pending}); } if (!CanApplyOuterCorrection(iteration)) { double full_hidden = 0.0; @@ -1059,11 +1122,15 @@ Result ExecuteReducedHht( correction.Value().p[i] / (parameters.beta * h * h); reduced_physical[i] = acceleration_rate * correction.Value().p[i]; } - CorrectionRecord(correction.Value().coordinate, + const bool private_correction_pass = CorrectionRecord(correction.Value().coordinate, correction.Value().velocity, correction.Value().acceleration, full_algorithmic, reduced_algorithmic, reduced_physical, chart.Value(), dofs, length, h, parameters.alpha, &record); + record.strict_correction_evaluated = true; + if (!private_correction_pass) + return fail(Failure(domain, "reduced-dynamic-update-inconsistency", + "The full/private correction acceleration arithmetic is inconsistent.")); translation_private_algorithmic_pass = true; for (std::size_t i = 0U; i < chart.Value().p_dofs.size(); ++i) if (PhysicalCoordinateChannel(chart.Value(), i) == 0U) @@ -1115,16 +1182,21 @@ Result ExecuteReducedHht( iteration_step.TakeCandidate()); }, [&](ReducedStepCandidate&& accepted_candidate) -> Status { - status = state.CommitAcceptedFrame(accepted_candidate.frame); - if (!status.IsOk()) - return report_failure(status, "accepted"); - history = std::move(accepted_candidate.history); + // Stage the complete accepted history. Output failure discards all + // three snapshots; publication below consists only of noexcept moves. + auto staged_state = state; + status = staged_state.CommitAcceptedFrame(accepted_candidate.frame); + if (!status.IsOk()) return report_failure(status, "accepted"); record.stage = "accepted"; record.outcome = "accepted"; + record.accepted = true; status = log.Write(record); - if (status.IsOk()) - status = log.Flush(); - return status; + if (status.IsOk()) status = log.Flush(); + if (!status.IsOk()) return status; + state = std::move(staged_state); + history = std::move(accepted_candidate.history); + flux_history = accepted_candidate.flux_history; + return Status::Ok(); }); if (!status.IsOk()) return Result::Failure(status); } diff --git a/src/fesa/analysis/newton_iteration_log.cpp b/src/fesa/analysis/newton_iteration_log.cpp index 895364e..517515b 100644 --- a/src/fesa/analysis/newton_iteration_log.cpp +++ b/src/fesa/analysis/newton_iteration_log.cpp @@ -261,10 +261,100 @@ Result SerializeRecord(const NewtonIterationLogRecord& record, return Result::Success(std::move(line)); } +// Candidate policy has separate names and units; v1 evidence is never redefined. +Result SerializeFieldRecord( + const DynamicNewtonIterationLogRecord& record, + const std::filesystem::path& path) { + if (record.policy_id != "abaqus2025-field-on-fesa-v1" || + record.time < 0.0 || record.stage.empty() || record.outcome.empty()) + return Result::Failure( + LogOutputFailure(path, "Invalid HHT field-policy log identity.")); + std::vector fields; + const auto number = [&](const std::string& key, + const std::optional& value) -> Status { + if (!value.has_value()) { + AppendField(&fields, key, "not-evaluated"); + return Status::Ok(); + } + auto formatted = FormatFiniteNumber(*value, path); + if (!formatted.HasValue()) return formatted.GetStatus(); + if (*value < 0.0) + return LogOutputFailure(path, "Negative HHT field norm or scale."); + AppendField(&fields, key, formatted.Value()); + return Status::Ok(); + }; + const auto boolean = [&](const std::string& key, const std::optional& value) { + AppendField(&fields, key, value.has_value() ? (*value ? "1" : "0") : "not-evaluated"); + }; + AppendField(&fields, "record", "dynamic-newton-iteration"); + AppendField(&fields, "schema", "hht-log-v2"); + AppendField(&fields, "policy_id", record.policy_id); + AppendField(&fields, "step", EscapeText(record.step_name)); + AppendField(&fields, "increment", std::to_string(record.increment_index)); + auto status = number("time", record.time); + if (!status.IsOk()) return Result::Failure(status); + AppendField(&fields, "iteration", std::to_string(record.iteration_index)); + AppendField(&fields, "stage", EscapeText(record.stage)); + AppendField(&fields, "outcome", EscapeText(record.outcome)); + AppendField(&fields, "failure_code", EscapeText(record.failure_code)); + boolean("field_pass", record.field_pass); + boolean("accepted", record.accepted); + for (std::size_t c = 0; c < 2; ++c) { + const auto& e = record.fields[c]; + const std::string prefix = c == 0 ? "force_" : "moment_"; + const std::array>, 9> numbers{{ + {"residual_numerator", e.residual_numerator}, {"residual_scale", e.residual_scale}, + {"spatial_flux", e.spatial_flux}, {"time_average_flux", e.time_average_flux}, + {"correction_numerator", e.correction_numerator}, {"increment_norm", e.increment_norm}, + {"estimated_correction", e.estimated_correction}, {"residual_threshold", e.residual_threshold}, + {"correction_threshold", e.correction_threshold}}}; + for (const auto& entry : numbers) { + status = number(prefix + entry.first, entry.second); + if (!status.IsOk()) return Result::Failure(status); + } + AppendField(&fields, prefix + "residual_unit", c == 0 ? "force" : "moment"); + AppendField(&fields, prefix + "correction_unit", c == 0 ? "length" : "rad"); + if (e.branch != "normal" && e.branch != "zero-flux" && e.branch != "linear" && + e.branch != "nonquadratic" && e.branch != "not-evaluated") + return Result::Failure(LogOutputFailure(path, "Invalid HHT field branch.")); + AppendField(&fields, prefix + "branch", e.branch); + boolean(prefix + "residual_pass", e.residual_pass); + boolean(prefix + "correction_pass", e.correction_pass); + boolean(prefix + "field_pass", e.field_pass); + } + for (std::size_t c = 0; c < 2; ++c) { + const std::string suffix = c == 0 ? "u" : "phi"; + const auto strict = [&](double value, bool evaluated) { + return evaluated ? std::optional{value} : std::nullopt; + }; + const std::array>, 6> numbers{{ + {"residual_" + suffix + "_numerator", strict(record.residual_numerator[c], record.strict_residual_evaluated)}, + {"residual_" + suffix + "_scale", strict(record.residual_scale[c], record.strict_residual_evaluated)}, + {"coordinate_correction_" + suffix, strict(record.coordinate_correction[c], record.strict_correction_evaluated)}, + {"velocity_correction_" + suffix, strict(record.velocity_correction[c], record.strict_correction_evaluated)}, + {"acceleration_correction_" + suffix, strict(record.acceleration_correction[c], record.strict_correction_evaluated)}, + {"algorithmic_correction_" + suffix, strict(record.algorithmic_correction[c], record.strict_correction_evaluated)}}}; + for (const auto& entry : numbers) { + status = number("strict_" + entry.first, entry.second); + if (!status.IsOk()) return Result::Failure(status); + } + } + boolean("strict_residual_pass", record.strict_residual_evaluated ? std::optional{record.residual_pass} : std::nullopt); + boolean("strict_correction_pass", record.strict_correction_evaluated ? std::optional{record.correction_pass} : std::nullopt); + std::string line; + for (const auto& field : fields) { + if (!line.empty()) line += '\t'; + line += field; + } + line += '\n'; + return Result::Success(std::move(line)); +} + /// @brief Serializes dynamic evidence without static load-factor semantics. Result SerializeRecord( const DynamicNewtonIterationLogRecord& record, const std::filesystem::path& log_path) { + if (!record.policy_id.empty()) return SerializeFieldRecord(record, log_path); const auto time = FormatFiniteNumber(record.time, log_path); if (!time.HasValue()) return Result::Failure(time.GetStatus()); if (record.time < 0.0 || record.stage.empty() || record.outcome.empty()) { diff --git a/tests/unit/analysis/hht_dynamic_algebra_test.cpp b/tests/unit/analysis/hht_dynamic_algebra_test.cpp index f77f540..a6da082 100644 --- a/tests/unit/analysis/hht_dynamic_algebra_test.cpp +++ b/tests/unit/analysis/hht_dynamic_algebra_test.cpp @@ -669,4 +669,39 @@ TEST(HhtDynamicAlgebra, RejectsAdmissibleNullRankChangeBeforeAcceptingFrame) { EXPECT_EQ(frame.hidden_constraint[0], 23.0); } +TEST(HhtEndpointFields, MeasuresActualFullChartAcrossNonlinearRootAndPartialBoundary) { + const auto fixture = MakeFixture(1U, true); + auto trial = EmptyTrial(fixture), old = EmptyTrial(fixture); + trial.residual = fesa::Vector{24U}; + fesa::Vector previous{24U}; + // Independent nonlinear graph q_phiZ=p^2: p .1 -> .2. The + // previous tangent lift is .02, while the actual endpoint change is .03. + previous[0] = .1; previous[5] = .01; + trial.generalized_coordinates[0] = .2; + trial.generalized_coordinates[5] = .04; + trial.residual[0] = 3; trial.residual[5] = -5; + // Large constrained values and spatial support moments are not free + // generalized endpoint evidence and cannot change either scalar maximum. + trial.generalized_coordinates[3] = 1000; + trial.residual[3] = 2000; + trial.reaction = fesa::Vector{24U}; trial.reaction[5] = 3000; + const auto result = fesa::hht_internal::MeasureEndpointFields(fixture.dofs, trial, old, previous); + ASSERT_TRUE(result.HasValue()); + EXPECT_TRUE(result.Value()[0].active); EXPECT_TRUE(result.Value()[1].active); + EXPECT_DOUBLE_EQ(result.Value()[0].residual, 3); + EXPECT_DOUBLE_EQ(result.Value()[1].residual, 5); + EXPECT_DOUBLE_EQ(result.Value()[0].increment, .2); + EXPECT_DOUBLE_EQ(result.Value()[1].increment, .04); + EXPECT_DOUBLE_EQ(*result.Value()[0].correction, .2 - .1); + EXPECT_DOUBLE_EQ(*result.Value()[1].correction, .04 - .01); + EXPECT_NE(*result.Value()[1].correction, 2 * .1 * (.2 - .1)); + const auto initial = fesa::hht_internal::MeasureEndpointFields(fixture.dofs, trial, old, std::nullopt); + ASSERT_TRUE(initial.HasValue()); + EXPECT_FALSE(initial.Value()[0].correction.has_value()); + EXPECT_FALSE(initial.Value()[1].correction.has_value()); + previous[0] = -std::numeric_limits::max(); + trial.generalized_coordinates[0] = std::numeric_limits::max(); + EXPECT_FALSE(fesa::hht_internal::MeasureEndpointFields(fixture.dofs, trial, old, previous).HasValue()); +} + } // namespace diff --git a/tests/unit/analysis/hht_dynamic_procedure_test.cpp b/tests/unit/analysis/hht_dynamic_procedure_test.cpp index 9e87653..1e97734 100644 --- a/tests/unit/analysis/hht_dynamic_procedure_test.cpp +++ b/tests/unit/analysis/hht_dynamic_procedure_test.cpp @@ -8,9 +8,16 @@ #include #include #include +#include +#include #include #include "fesa/analysis/newton_iteration_log.h" +#include "fesa/analysis/analysis_model.h" +#include "analysis/hht_dynamic_internal.h" +#include "analysis/hht_reduced_integrator.h" +#include "analysis/hht_reduced_procedure.h" +#include "analysis/equilibrium_algorithm.h" #include "fesa/assembly/parallel_for.h" #include "fesa/solvers/linear/mkl_pardiso_solver.h" #include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" @@ -569,3 +576,355 @@ TEST(HhtDynamicProcedure, } } // namespace + +namespace { +std::map LogFields(const std::string& line) { + std::map fields; + std::istringstream input(line); + std::string token; + while (std::getline(input, token, '\t')) { + const auto equal = token.find('='); + if (equal != std::string::npos) fields.emplace(token.substr(0, equal), token.substr(equal + 1)); + } + return fields; +} +} + +namespace { +fesa::ModelDefinition SymmetricExtension() { + auto definition = MakeDefinition(false); + definition.dynamic_step->initial_increment = 0.1; + definition.dynamic_step->time_period = 0.1; + definition.dynamic_step->increment_count = 1U; + definition.node_sets.push_back({"Left", std::nullopt, {0U, 3U}, {"dynamic-state.inp", 30U}}); + definition.node_sets.push_back({"Right", std::nullopt, {1U, 2U}, {"dynamic-state.inp", 31U}}); + definition.dynamic_step->boundaries = { + {"All", 2, 6, 0.0, {"dynamic-state.inp", 32U}}}; + definition.dynamic_step->loads = { + {"Left", 1, -1.0, {"dynamic-state.inp", 33U}}, + {"Right", 1, 1.0, {"dynamic-state.inp", 34U}}}; + return definition; +} +} + +TEST(HhtDynamicProcedure, FieldPolicyActuallyTerminatesSymmetricExtensionBeforeStrict) { + auto definition = SymmetricExtension(); + fesa::SerialParallelFor parallel; + auto run = RunRigidTranslation(definition, parallel, "field-policy-extension"); + ASSERT_TRUE(run.HasValue()) << (run.HasValue() ? "" : FailureText(run.GetStatus())); + const auto& frame = run.Value().execution.state.CommittedState(); + ASSERT_GT(frame.iteration_count, 0U); + ASSERT_LE(frame.iteration_count, 100U); + const auto& text = run.Value().log; + const auto accepted = text.rfind("record=dynamic-newton-iteration"); + ASSERT_NE(accepted, std::string::npos); + const auto line = text.substr(accepted); + EXPECT_NE(line.find("\tschema=hht-log-v2\t"), std::string::npos); + EXPECT_NE(line.find("\tpolicy_id=abaqus2025-field-on-fesa-v1\t"), std::string::npos); + EXPECT_NE(line.find("\tstep=Transient Load\t"), std::string::npos); + EXPECT_NE(line.find("\tfield_pass=1\taccepted=1\t"), std::string::npos); + EXPECT_NE(line.find("\tstrict_correction_pass=0"), std::string::npos); + for (const auto value : frame.global_equilibrium) EXPECT_LE(std::abs(value), 1e-8); + auto fields = LogFields(line.substr(0, line.find('\n'))); + double rho = 0, increment = 0, flux = 4.0; + for (std::size_t node = 0; node < 4; ++node) { + rho = (std::max)(rho, std::abs(frame.residual[6 * node])); + increment = (std::max)(increment, std::abs(frame.generalized_coordinates[6 * node])); + for (std::size_t c = 0; c < 3; ++c) + flux += std::abs(frame.internal_force[6 * node + c] + frame.inertia_force[6 * node + c]); + } + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_residual_numerator")), rho); + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_increment_norm")), increment); + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_spatial_flux")), flux / 24.0); + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_time_average_flux")), flux / 24.0); + EXPECT_LE(rho, .005 * std::stod(fields.at("force_residual_scale"))); + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_residual_threshold")), .005); + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_correction_threshold")), .01); + EXPECT_EQ(fields.at("force_residual_pass"), "1"); + EXPECT_EQ(fields.at("force_correction_pass"), "1"); + const double c = std::stod(fields.at("force_correction_numerator")); + const auto& estimated = fields.at("force_estimated_correction"); + EXPECT_TRUE(c <= .01 * increment || + (estimated != "not-evaluated" && std::stod(estimated) <= .01 * increment)); + EXPECT_EQ(fields.at("force_branch"), "normal"); + EXPECT_EQ(fields.at("moment_branch"), "not-evaluated"); + // Symmetric monotone extension has q=(-D,+D,+D,-D). This independent + // endpoint oracle includes the selected iteration-zero start, not the lift. + std::istringstream records(text); + std::string row; + double previous_D = 0; + bool first = true; + while (std::getline(records, row)) { + const auto values = LogFields(row); + if (values.at("schema") != "hht-log-v2") continue; + const double D = std::stod(values.at("force_increment_norm")); + if (first) { + EXPECT_EQ(values.at("force_correction_numerator"), "not-evaluated"); + EXPECT_EQ(values.at("force_estimated_correction"), "not-evaluated"); + EXPECT_EQ(values.at("field_pass"), "0"); + } else { + EXPECT_NEAR(std::stod(values.at("force_correction_numerator")), + std::abs(D - previous_D), 1e-16); + } + previous_D = D; + first = false; + } +} + +TEST(HhtDynamicProcedure, FieldPassDoesNotBypassPhysicalAcceptanceAndLogsActualCause) { + auto definition = MakeDefinition(false); + definition.materials[0].youngs_modulus = 1e9; + definition.dynamic_step->increment_count = 1; + definition.dynamic_step->time_period = definition.dynamic_step->initial_increment; + definition.node_sets[0].node_indices = {0U, 3U}; + definition.node_sets[1].node_indices = {2U}; + definition.dynamic_step->boundaries = {{"All", 1, 6, 0, {"dynamic-state.inp", 11U}}}; + definition.dynamic_step->loads = {{"Loaded", 3, -40000, {"dynamic-state.inp", 12U}}}; + auto domain = fesa::Domain::Create(definition); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + fesa::NewtonIterationLogWriter log; + const auto path = std::filesystem::temp_directory_path() / "fesa-field-post-gate.log"; + ASSERT_TRUE(log.Open(path).IsOk()); + auto result = fesa::HhtDynamicProcedureExecutor(parallel, mass, tangent).Execute(domain.Value(), log); + ASSERT_TRUE(log.Close().IsOk()); + std::ifstream input(path, std::ios::binary); + const std::string text{std::istreambuf_iterator{input}, {}}; + input.close(); + std::filesystem::remove(path); + ASSERT_FALSE(result.HasValue()); + const auto row = text.substr(text.rfind("record=")); + const auto fields = LogFields(row.substr(0, row.find('\n'))); + EXPECT_EQ(fields.at("field_pass"), "1"); + EXPECT_EQ(fields.at("accepted"), "0"); + EXPECT_EQ(fields.at("outcome"), "failed"); + const auto cause = std::find_if(result.GetStatus().Diagnostics().begin(), + result.GetStatus().Diagnostics().end(), [](const auto& d) { + return d.code != "dynamic-procedure-context"; + }); + ASSERT_NE(cause, result.GetStatus().Diagnostics().end()); + EXPECT_EQ(fields.at("failure_code"), cause->code) << FailureText(result.GetStatus()); + EXPECT_EQ(text.find("\taccepted=1"), std::string::npos); +} + +TEST(HhtDynamicProcedure, AcceptedFluxHistoryCountsEachIncrementExactlyOnce) { + auto definition = SymmetricExtension(); + definition.dynamic_step->time_period = .2; + definition.dynamic_step->increment_count = 2; + fesa::SerialParallelFor parallel; + auto run = RunRigidTranslation(definition, parallel, "field-history"); + ASSERT_TRUE(run.HasValue()) << (run.HasValue() ? "" : FailureText(run.GetStatus())); + const auto& frames = run.Value().execution.state.AcceptedFrames(); + ASSERT_EQ(frames.size(), 3U); + double sum = 0, maximum_sum = 0; + std::size_t n = 0; + std::istringstream input(run.Value().log); + std::string row; + while (std::getline(input, row)) { + const auto fields = LogFields(row); + if (fields.at("schema") != "hht-log-v2" || fields.at("accepted") != "1") continue; + ++n; + const auto& frame = frames[n]; + std::vector flux; + for (std::size_t node = 0; node < 4; ++node) + for (std::size_t c = 0; c < 3; ++c) + flux.push_back(std::abs(frame.internal_force[6 * node + c] + frame.inertia_force[6 * node + c])); + for (int node = 0; node < 4; ++node) { + flux.push_back(1); flux.push_back(0); flux.push_back(0); + } + const double m = *std::max_element(flux.begin(), flux.end()); + const double H = n == 1 ? 0 : maximum_sum / static_cast(n - 1); + double selected_sum = 0; + std::size_t count = 0; + for (double x : flux) + if (m < .1 * H || x >= 1e-5 * H) { selected_sum += x; ++count; } + ASSERT_GT(count, 0U); + const double spatial = selected_sum / static_cast(count); + sum += spatial; maximum_sum += m; + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_spatial_flux")), spatial); + EXPECT_DOUBLE_EQ(std::stod(fields.at("force_time_average_flux")), sum / static_cast(n)); + EXPECT_EQ(std::stoul(fields.at("iteration")), frame.iteration_count); + } + EXPECT_EQ(n, 2U); +} + +TEST(HhtDynamicProcedure, LoadFluxPreservesSourceCancellationAndExplicitZeroEntries) { + auto definition = MakeDefinition(true); + definition.dynamic_step->loads = { + {"Loaded", 1, 3, {"dynamic-state.inp", 30U}}, + {"Loaded", 1, -3, {"dynamic-state.inp", 31U}}, + {"Loaded", 2, 0, {"dynamic-state.inp", 32U}}}; + auto domain = fesa::Domain::Create(definition); + 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()); + std::array, 2> flux; + const auto loads = fesa::hht_internal::AssembleLoads(domain.Value(), dofs.Value(), &flux); + ASSERT_TRUE(loads.HasValue()); + EXPECT_DOUBLE_EQ(loads.Value()[0][0], 0); + EXPECT_DOUBLE_EQ(loads.Value()[1][0], 6); + ASSERT_EQ(flux[0].size(), 9U); + ASSERT_EQ(flux[1].size(), 9U); + const std::array expected{3, 0, 0, -3, 0, 0, 0, 0, 0}; + for (std::size_t i = 0; i < expected.size(); ++i) { + EXPECT_DOUBLE_EQ(flux[0][i].internal, expected[i]); + EXPECT_DOUBLE_EQ(flux[0][i].inertia, 0); + EXPECT_DOUBLE_EQ(flux[1][i].internal, 0); + } +} + +namespace { +// All-fixed fixture has no private coordinates. This adapter exercises the +// real production acceptance/output transaction without a test solver path. +class EmptyCoordinateIntegrator final : public fesa::hht_reduced_internal::ShellHhtIntegrator { + public: + fesa::hht_reduced_internal::ReducedShellPredictor Predict( + const fesa::hht_reduced_internal::ReducedHistory& history, double, + const fesa::HhtParameters&) const override { + EXPECT_EQ(history.p.Size(), 0U); + return {}; + } + void EvaluateRates(const fesa::hht_reduced_internal::ReducedShellRateRequest&, + fesa::Vector& v, fesa::Vector& a) const override { + EXPECT_EQ(v.Size(), 0U); EXPECT_EQ(a.Size(), 0U); + } + fesa::Result ReconstructPhysicalRates( + const fesa::hht_reduced_internal::FixedChart& chart, + const fesa::hht_reduced_internal::EmbeddingRoot& root, + const fesa::Vector& v, const fesa::Vector& a, + fesa::TangentSolver& solver) const override { + return fesa::hht_reduced_internal::ReconstructPhysicalRates(chart, root, v, a, solver); + } +}; +class CloseAtCandidate final : public fesa::EquilibriumAlgorithm { + public: + explicit CloseAtCandidate(fesa::NewtonIterationLogWriter& log) : log_(log) {} + fesa::Status SolveCurrentStep(std::uint32_t first, std::uint32_t last, + const Evaluation& evaluate) override { + EXPECT_EQ(first, 0U); EXPECT_EQ(last, 100U); + auto value = evaluate(first); + if (!value.HasValue()) return value.GetStatus(); + EXPECT_TRUE(value.Value()); + return log_.Close(); + } + private: + fesa::NewtonIterationLogWriter& log_; +}; +} + +TEST(HhtDynamicProcedure, FailedAcceptanceOutputPublishesNoCandidateAndFreshRetryIsUnchanged) { + auto domain = fesa::Domain::Create(MakeDefinition()); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + fesa::NewtonIterationLogWriter log; + const auto path = std::filesystem::temp_directory_path() / "fesa-field-output-fail.log"; + ASSERT_TRUE(log.Open(path).IsOk()); + EmptyCoordinateIntegrator integrator; + CloseAtCandidate algorithm(log); + const auto result = fesa::ExecuteReducedHht(domain.Value(), log, parallel, mass, tangent, integrator, algorithm); + ASSERT_FALSE(result.HasValue()); + EXPECT_EQ(result.GetStatus().Category(), fesa::FailureCategory::kOutput); + std::ifstream input(path, std::ios::binary); + const std::string text{std::istreambuf_iterator{input}, {}}; + input.close(); std::filesystem::remove(path); + EXPECT_EQ(text.find("accepted=1"), std::string::npos); + auto first = RunRigidTranslation(MakeDefinition(), parallel, "field-retry"); + auto second = RunRigidTranslation(MakeDefinition(), parallel, "field-retry-repeat"); + ASSERT_TRUE(first.HasValue()); ASSERT_TRUE(second.HasValue()); + ExpectSameDynamicHistory(first.Value(), second.Value()); +} + +namespace { +class FailNextSelectedTrial final : public fesa::TangentSolver { + public: + fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override { + if (outer_solved_) return fesa::Status::Failure(fesa::FailureCategory::kSolver, + {{fesa::Severity::kError, "injected-next-trial", {}, "DYNAMIC", {}, "Fail before the next selected field evaluation."}}); + large_ = matrix.Rows() > 4; + return solver_.FactorizeFresh(matrix); + } + fesa::Status SolveCorrection(const fesa::Vector& rhs, fesa::Vector& solution) const override { + auto status = solver_.SolveCorrection(rhs, solution); + if (status.IsOk() && large_) outer_solved_ = true; + return status; + } + private: + fesa::MklPardisoTangentSolver solver_; + bool large_ = false; + mutable bool outer_solved_ = false; +}; +} + +TEST(HhtDynamicProcedure, FailureBeforeNextEvaluationDoesNotReusePreviousFieldEvidence) { + auto domain = fesa::Domain::Create(MakeFiniteRotationDefinition()); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + FailNextSelectedTrial tangent; + fesa::NewtonIterationLogWriter log; + const auto path = std::filesystem::temp_directory_path() / "fesa-field-next-fail.log"; + ASSERT_TRUE(log.Open(path).IsOk()); + const auto result = fesa::HhtDynamicProcedureExecutor(parallel, mass, tangent).Execute(domain.Value(), log); + ASSERT_TRUE(log.Close().IsOk()); + std::ifstream input(path, std::ios::binary); + const std::string text{std::istreambuf_iterator{input}, {}}; + input.close(); std::filesystem::remove(path); + ASSERT_FALSE(result.HasValue()); + const auto row = text.substr(text.rfind("record=")); + const auto fields = LogFields(row.substr(0, row.find('\n'))); + EXPECT_EQ(fields.at("failure_code"), "injected-next-trial") << FailureText(result.GetStatus()); + EXPECT_EQ(fields.at("iteration"), "1"); + EXPECT_EQ(fields.at("force_residual_numerator"), "not-evaluated"); + EXPECT_EQ(fields.at("moment_residual_scale"), "not-evaluated"); + EXPECT_EQ(fields.at("field_pass"), "not-evaluated"); + EXPECT_EQ(fields.at("accepted"), "0"); +} + +namespace { +class EvaluateEndpoint100 final : public fesa::EquilibriumAlgorithm { + public: + fesa::Status SolveCurrentStep(std::uint32_t first, std::uint32_t last, + const Evaluation& evaluate) override { + EXPECT_EQ(first, 0U); EXPECT_EQ(last, 100U); + auto value = evaluate(last); + if (!value.HasValue()) return value.GetStatus(); + EXPECT_TRUE(value.Value()); + return fesa::Status::Ok(); + } +}; +} + +TEST(HhtDynamicProcedure, FieldAcceptanceEvaluatesProductionEndpoint100BeforeCorrectionBudget) { + auto domain = fesa::Domain::Create(MakeDefinition()); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor parallel; + fesa::MklPardisoSolver mass; + fesa::MklPardisoTangentSolver tangent; + fesa::NewtonIterationLogWriter log; + const auto path = std::filesystem::temp_directory_path() / "fesa-field-endpoint100.log"; + ASSERT_TRUE(log.Open(path).IsOk()); + EmptyCoordinateIntegrator integrator; + EvaluateEndpoint100 algorithm; + // Tests the real endpoint callback at the inclusive boundary. No physical + // trajectory is forced through 100 corrections; all fixed coordinates stay zero. + const auto result = fesa::ExecuteReducedHht(domain.Value(), log, parallel, mass, tangent, integrator, algorithm); + ASSERT_TRUE(log.Close().IsOk()); + std::ifstream input(path, std::ios::binary); + const std::string text{std::istreambuf_iterator{input}, {}}; + input.close(); std::filesystem::remove(path); + ASSERT_TRUE(result.HasValue()) << (result.HasValue() ? "" : FailureText(result.GetStatus())); + ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 3U); + EXPECT_EQ(result.Value().state.AcceptedFrames()[1].iteration_count, 100U); + EXPECT_EQ(result.Value().state.AcceptedFrames()[2].iteration_count, 100U); + const auto row = text.substr(text.rfind("record=")); + const auto fields = LogFields(row.substr(0, row.find('\n'))); + EXPECT_EQ(fields.at("iteration"), "100"); + EXPECT_EQ(fields.at("field_pass"), "1"); + EXPECT_EQ(fields.at("accepted"), "1"); +} diff --git a/tests/unit/analysis/newton_iteration_log_test.cpp b/tests/unit/analysis/newton_iteration_log_test.cpp index afebbd9..bd97b6c 100644 --- a/tests/unit/analysis/newton_iteration_log_test.cpp +++ b/tests/unit/analysis/newton_iteration_log_test.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -562,4 +563,149 @@ TEST(DynamicNewtonIterationLog, ExpectLogOutputFailure(unopened.Write(CompleteDynamicRecord())); } +TEST(DynamicFieldIterationLog, KeepsOrderedUnitsMissingEvidenceAndEscapedSourceName) { + TempDirectory directory{"field-v2"}; + const auto path = directory.Path() / "log.txt"; + auto record = CompleteDynamicRecord(); + record.policy_id = "abaqus2025-field-on-fesa-v1"; + record.step_name = "\xED\x95\x9C\xEA\xB8\x80=Step\\tail\tline\nend\r"; + record.stage = "evaluation"; + record.outcome = "candidate"; + record.fields[0].residual_numerator = 1.2345678901234567; + record.fields[0].increment_norm = 0; + fesa::NewtonIterationLogWriter writer; + ASSERT_TRUE(writer.Open(path).IsOk()); + ASSERT_TRUE(writer.Write(record).IsOk()); + ASSERT_TRUE(writer.Close().IsOk()); + const auto text = ReadLogBytes(path); + EXPECT_EQ(std::count(text.begin(), text.end(), '\n'), 1); + EXPECT_EQ(text.find('\r'), std::string::npos); + const auto fields = SplitTabs(text.substr(0, text.size() - 1)); + const std::vector base{"record", "schema", "policy_id", "step", "increment", "time", "iteration", "stage", "outcome", "failure_code", "field_pass", "accepted"}; + for (std::size_t i = 0; i < base.size(); ++i) + EXPECT_EQ(fields[i].substr(0, fields[i].find('=')), base[i]); + const std::vector suffixes{ + "residual_numerator", "residual_scale", "spatial_flux", "time_average_flux", + "correction_numerator", "increment_norm", "estimated_correction", + "residual_threshold", "correction_threshold", "residual_unit", "correction_unit", + "branch", "residual_pass", "correction_pass", "field_pass"}; + for (std::size_t c = 0; c < 2; ++c) + for (std::size_t i = 0; i < suffixes.size(); ++i) + EXPECT_EQ(fields[base.size() + c * suffixes.size() + i].substr( + 0, fields[base.size() + c * suffixes.size() + i].find('=')), + (c == 0 ? "force_" : "moment_") + suffixes[i]); + EXPECT_EQ(FieldValue(fields, "step"), "step=\xED\x95\x9C\xEA\xB8\x80=Step\\\\tail\\tline\\nend\\r"); + EXPECT_EQ(FieldValue(fields, "force_correction_numerator"), "force_correction_numerator=not-evaluated"); + EXPECT_EQ(FieldValue(fields, "force_increment_norm"), "force_increment_norm=0"); + EXPECT_EQ(FieldValue(fields, "strict_residual_pass"), "strict_residual_pass=not-evaluated"); + EXPECT_EQ(FieldValue(fields, "moment_correction_unit"), "moment_correction_unit=rad"); + EXPECT_DOUBLE_EQ(std::stod(ValueAfterEquals(FieldValue(fields, "force_residual_numerator"))), + *record.fields[0].residual_numerator); +} + +TEST(DynamicFieldIterationLog, InvalidFieldEvidenceFailsWithoutWritingPartialRecord) { + TempDirectory directory{"field-v2-invalid"}; + const auto path = directory.Path() / "log.txt"; + for (int kind = 0; kind < 4; ++kind) { + SCOPED_TRACE(kind); + fesa::NewtonIterationLogWriter writer; + ASSERT_TRUE(writer.Open(path).IsOk()); + auto record = CompleteDynamicRecord(); + record.policy_id = "abaqus2025-field-on-fesa-v1"; + ASSERT_TRUE(writer.Write(record).IsOk()); + ASSERT_TRUE(writer.Flush().IsOk()); + const auto prefix = ReadLogBytes(path); + if (kind == 0) record.fields[0].branch = "unknown"; + if (kind == 1) record.fields[0].estimated_correction = std::numeric_limits::infinity(); + if (kind == 2) record.fields[1].residual_scale = -1; + if (kind == 3) record.policy_id = "unknown"; + ExpectLogOutputFailure(writer.Write(record)); + ExpectLogOutputFailure(writer.Flush()); + EXPECT_EQ(ReadLogBytes(path), prefix); + } +} + +// Independent strict consumer: production has a writer, no public log reader. +bool ReadFieldV2(const std::string& bytes) { + if (bytes.empty() || bytes.back() != '\n' || bytes.find('\r') != std::string::npos || + std::count(bytes.begin(), bytes.end(), '\n') != 1) return false; + std::map values; + for (const auto& field : SplitTabs(bytes.substr(0, bytes.size() - 1))) { + const auto equal = field.find('='); + if (equal == std::string::npos || !values.emplace(field.substr(0, equal), field.substr(equal + 1)).second) + return false; + } + for (const auto* key : {"record", "schema", "policy_id", "step", "increment", "time", "iteration", "stage", "outcome", "failure_code", "field_pass", "accepted"}) + if (!values.count(key)) return false; + if (values["record"] != "dynamic-newton-iteration" || values["schema"] != "hht-log-v2" || + values["policy_id"] != "abaqus2025-field-on-fesa-v1") return false; + const auto numeric = [](const std::string& text, bool missing) { + if (missing && text == "not-evaluated") return true; + double value = 0; + const auto parsed = std::from_chars(text.data(), text.data() + text.size(), value); + return parsed.ec == std::errc{} && parsed.ptr == text.data() + text.size() && std::isfinite(value) && value >= 0; + }; + const auto boolean = [](const std::string& text, bool missing) { + return text == "0" || text == "1" || (missing && text == "not-evaluated"); + }; + for (const auto* key : {"increment", "iteration"}) { + std::uint64_t value = 0; + const auto& text = values[key]; + const auto parsed = std::from_chars(text.data(), text.data() + text.size(), value); + if (parsed.ec != std::errc{} || parsed.ptr != text.data() + text.size()) return false; + } + if (!numeric(values["time"], false) || !boolean(values["field_pass"], true) || !boolean(values["accepted"], false)) return false; + for (const auto* prefix : {"force_", "moment_"}) { + for (const auto* suffix : {"residual_numerator", "residual_scale", "spatial_flux", "time_average_flux", "correction_numerator", "increment_norm", "estimated_correction", "residual_threshold", "correction_threshold"}) { + const auto found = values.find(std::string(prefix) + suffix); + if (found == values.end() || !numeric(found->second, true)) return false; + } + for (const auto* suffix : {"residual_pass", "correction_pass", "field_pass"}) { + const auto found = values.find(std::string(prefix) + suffix); + if (found == values.end() || !boolean(found->second, true)) return false; + } + const auto branch = values.find(std::string(prefix) + "branch"); + if (branch == values.end() || (branch->second != "normal" && branch->second != "zero-flux" && + branch->second != "linear" && branch->second != "nonquadratic" && branch->second != "not-evaluated")) return false; + if (values[std::string(prefix) + "residual_unit"] != (std::string(prefix) == "force_" ? "force" : "moment") || + values[std::string(prefix) + "correction_unit"] != (std::string(prefix) == "force_" ? "length" : "rad")) return false; + } + for (const auto& entry : values) + for (std::size_t i = 0; i < entry.second.size(); ++i) + if (entry.second[i] == '\\') { + if (++i == entry.second.size() || (entry.second[i] != '\\' && entry.second[i] != 't' && entry.second[i] != 'n' && entry.second[i] != 'r')) return false; + } + return true; +} + +TEST(DynamicFieldIterationLog, IndependentReaderRejectsMalformedVersionedRecords) { + TempDirectory directory{"field-reader"}; + const auto path = directory.Path() / "log.txt"; + fesa::NewtonIterationLogWriter writer; + ASSERT_TRUE(writer.Open(path).IsOk()); + auto record = CompleteDynamicRecord(); + record.policy_id = "abaqus2025-field-on-fesa-v1"; + ASSERT_TRUE(writer.Write(record).IsOk()); + ASSERT_TRUE(writer.Close().IsOk()); + const auto valid = ReadLogBytes(path); + ASSERT_TRUE(ReadFieldV2(valid)); + const auto replaced = [&](const std::string& from, const std::string& to) { + auto text = valid; + const auto start = text.find(from); + EXPECT_NE(start, std::string::npos); + text.replace(start, from.size(), to); + return text; + }; + EXPECT_FALSE(ReadFieldV2(replaced("schema=hht-log-v2", "schema=hht-log-v1"))); + EXPECT_FALSE(ReadFieldV2(replaced("\tincrement=2", "\tincrement=-1"))); + EXPECT_FALSE(ReadFieldV2(replaced("\ttime=0.25", "\ttime=nan"))); + EXPECT_FALSE(ReadFieldV2(replaced("\ttime=0.25", "\ttime=inf"))); + EXPECT_FALSE(ReadFieldV2(replaced("force_residual_numerator=not-evaluated", "force_residual_numerator=1junk"))); + EXPECT_FALSE(ReadFieldV2(replaced("force_branch=not-evaluated", "force_branch=unknown"))); + EXPECT_FALSE(ReadFieldV2(replaced("\tfield_pass=not-evaluated", "\tfield_pass=true"))); + EXPECT_FALSE(ReadFieldV2(replaced("\tstep=", "\tstep=bad\\x"))); + EXPECT_FALSE(ReadFieldV2(replaced("\taccepted=0", "\taccepted=0\taccepted=0"))); + EXPECT_FALSE(ReadFieldV2(replaced("\tforce_increment_norm=not-evaluated", ""))); +} + } // namespace