diff --git a/include/fesa/analysis/transient_heat_procedure.h b/include/fesa/analysis/transient_heat_procedure.h new file mode 100644 index 0000000..59f7a1d --- /dev/null +++ b/include/fesa/analysis/transient_heat_procedure.h @@ -0,0 +1,39 @@ +#ifndef FESA_ANALYSIS_TRANSIENT_HEAT_PROCEDURE_H_ +#define FESA_ANALYSIS_TRANSIENT_HEAT_PROCEDURE_H_ + +#include + +#include "fesa/analysis/thermal_analysis_state.h" +#include "fesa/core/status.h" + +namespace fesa { +class LinearSolver; +class ParallelFor; + +/// @brief Owns complete thermal history before any result output is opened. +struct ThermalExecutionResult { + ThermalAnalysisState state; + std::vector diagnostics; +}; + +/// @brief Executes fixed-increment backward Euler with checked thermal +/// recovery. +class TransientHeatProcedureExecutor final { + public: + /// @brief Retains caller-owned adapters, which must outlive Execute calls. + TransientHeatProcedureExecutor(const ParallelFor& parallel_for, + LinearSolver& linear_solver); + /// @brief Returns a complete independently owned state or a structured + /// failure. + /// @note Uses raw initial values; prescribed temperatures start at increment + /// one. No output files are opened. A failed trial is never committed to + /// history. + Result Execute(const Domain& domain) const; + + private: + const ParallelFor& parallel_for_; + LinearSolver& linear_solver_; +}; + +} // namespace fesa +#endif // FESA_ANALYSIS_TRANSIENT_HEAT_PROCEDURE_H_ diff --git a/src/fesa/analysis/thermal_transient_algebra.cpp b/src/fesa/analysis/thermal_transient_algebra.cpp new file mode 100644 index 0000000..7552e83 --- /dev/null +++ b/src/fesa/analysis/thermal_transient_algebra.cpp @@ -0,0 +1,503 @@ +#include "thermal_transient_algebra.h" + +#include +#include +#include +#include +#include + +#include "thermal_transient_bounds.h" + +namespace fesa::thermal_transient_internal { +namespace { +using bounds::Abs; +using bounds::Add; +using bounds::Div; +using bounds::Gamma; +using bounds::I; +using bounds::Magnitude; +using bounds::Mul; +using bounds::Roundoff; +using bounds::Sub; +using bounds::Tiny; +using bounds::Zero; + +/// @brief Reports certification failure without publishing a candidate. +Status Failure(const std::string& message) { + return Status::Failure(FailureCategory::kSolver, {{Severity::kError, + "thermal-increment-failed", + {}, + "HEAT TRANSFER", + "", + message}}); +} +/// @brief Checks every scalar while preserving signed temperatures and powers. +bool Finite(const Vector& vector) { + for (std::size_t i = 0; i < vector.Size(); ++i) + if (!std::isfinite(vector[i])) return false; + return true; +} +/// @brief Checks aligned full operators and their coefficient-error storage. +bool SamePattern(const SparseMatrix& a, const SparseMatrix& b) { + return a.Rows() == b.Rows() && a.Columns() == b.Columns() && + a.RowOffsets() == b.RowOffsets() && + a.ColumnIndices() == b.ColumnIndices(); +} +/// @brief Encloses a componentwise maximum for infinity norms. +I Max(I a, I b) { + return {std::max(a.lower, b.lower), std::max(a.upper, b.upper)}; +} +/// @brief Encloses the exact infinity norm of stored scalar entries. +I Norm(const Vector& v) { + I result; + for (std::size_t i = 0; i < v.Size(); ++i) + result = Max(result, Magnitude(v[i])); + return result; +} +/// @brief Bounds construction error against the independently enclosed value. +double Error(I exact, double rounded) { + return Abs(Sub(exact, I{rounded})).upper; +} +/// @brief Checks optional entrywise coefficient bounds without inventing +/// defaults. +bool Nonnegative(const SparseMatrix* error, const SparseMatrix& matrix) { + if (!error) return true; + if (!error->Validate().IsOk() || !SamePattern(*error, matrix)) return false; + for (const double x : error->Values()) + if (x < 0) return false; + return true; +} +/// @brief Validates the private generic algebra boundary before indexing. +Status Validate(const ThermalSystem& s, const ThermalPartition& p, + const Vector& old, double delta) { + const auto m = s.k.Rows(); + if (!s.k.Validate().IsOk() || !s.c.Validate().IsOk() || + !SamePattern(s.k, s.c) || s.k.Columns() != m || s.q.Size() != m || + old.Size() != m || !Finite(s.q) || !Finite(old) || + !Finite(p.prescribed) || !std::isfinite(delta) || delta <= 0 || + p.prescribed.Size() != p.constrained.size() || + p.free.size() + p.constrained.size() != m || + !Nonnegative(s.k_error, s.k) || !Nonnegative(s.c_error, s.c)) + return Failure( + "Invalid finite thermal operators, partition, or actual increment."); + if (s.q_error) { + if (s.q_error->Size() != m || !Finite(*s.q_error)) + return Failure("Invalid power error bounds."); + for (std::size_t i = 0; i < m; ++i) + if ((*s.q_error)[i] < 0) return Failure("Negative power error bound."); + } + std::vector used(m, false); + for (const auto* rows : {&p.free, &p.constrained}) { + if (!std::is_sorted(rows->begin(), rows->end())) + return Failure("Partition order is invalid."); + for (const auto row : *rows) { + if (row >= m || used[row]) + return Failure("Partition is not disjoint and exhaustive."); + used[row] = true; + } + } + return Status::Ok(); +} + +/// @brief Keeps rounded solve data and separate exact-algebra enclosures. +struct Equations { + SparseMatrix a; + Vector b; + std::vector exact_a; + std::vector exact_b; +}; +/// @brief Forms the reduced operator and raw-old-state backward Euler RHS. +Result Build(const ThermalSystem& s, const ThermalPartition& p, + const Vector& old, double delta) { + const auto n = p.free.size(), m = old.Size(); + std::vector free(m, n), constrained(m, p.constrained.size()); + for (std::size_t i = 0; i < n; ++i) free[p.free[i]] = i; + for (std::size_t i = 0; i < p.constrained.size(); ++i) + constrained[p.constrained[i]] = i; + SparsePattern pattern{{0}, {}}; + std::vector coo; + Vector b{n}; + std::vector exact_a, exact_b(n); + for (std::size_t f = 0; f < n; ++f) { + const auto row = p.free[f]; + b[f] = s.q[row]; + exact_b[f] = I{s.q[row]}; + for (auto j = s.c.RowOffsets()[row]; j < s.c.RowOffsets()[row + 1]; ++j) { + const auto col = s.c.ColumnIndices()[j]; + const double c = s.c.Values()[j], k = s.k.Values()[j]; + const I cd = Div(I{c}, I{delta}); + if (free[col] != n) { + coo.push_back({f, free[col], c / delta + k, 0, coo.size()}); + pattern.column_indices.push_back(free[col]); + exact_a.push_back(Add(cd, I{k})); + b[f] += (c / delta) * old[col]; + exact_b[f] = Add(exact_b[f], Mul(cd, I{old[col]})); + } else { + const double g = p.prescribed[constrained[col]]; + // Preserve raw old constrained capacity; never project frame zero. + b[f] += (c / delta) * (old[col] - g) - k * g; + exact_b[f] = Add(exact_b[f], + Sub(Mul(cd, Sub(I{old[col]}, I{g})), Mul(I{k}, I{g}))); + } + } + pattern.row_offsets.push_back(coo.size()); + } + auto a = SparseMatrix::FromCoo(n, n, std::move(coo), pattern); + if (!a.HasValue()) return Result::Failure(a.GetStatus()); + if (!Finite(b)) + return Result::Failure(Failure("Nonfinite backward Euler RHS.")); + return Result::Success({std::move(a.Value()), std::move(b), + std::move(exact_a), std::move(exact_b)}); +} + +/// @brief Certifies a residual using the reduced denominator and exact-zero +/// rule. +bool Certificate(double numerator, I denominator, double limit, double& beta) { + if (!std::isfinite(numerator) || !denominator.IsFinite()) return false; + if (denominator.lower <= 0) { + beta = 0; + return Zero(denominator) && numerator == 0; + } + beta = Div(I{numerator}, I{denominator.lower}).upper; + return std::isfinite(beta) && beta <= limit; +} + +/// @brief Checks reduced/full residuals, physical heat balance, and BE +/// dissipation. +/// @note Coefficient errors enter forward evidence only, never the solve +/// tolerance. +Status Check(const ThermalSystem& s, const ThermalPartition& p, + const Vector& old, double delta, const Equations& eq, + ThermalStepCandidate& candidate) { + const auto m = old.Size(), n = p.free.size(); + const double md = static_cast(m), nd = static_cast(n); + auto& ev = candidate.evidence; + const I beta_limit = + Add(Mul(I{nd}, Gamma(nd + 1)), Add(Gamma(2 * nd + 4), Gamma(4 * nd + 8))); + ev.beta_limit = beta_limit.lower; + I a_norm, t_norm, b_norm = Norm(eq.b); + Vector red{n}, red_error{n}; + std::vector construction(n); + for (std::size_t f = 0; f < n; ++f) { + I row_norm, scale, ea, identity; + double dot = 0; + for (auto j = eq.a.RowOffsets()[f]; j < eq.a.RowOffsets()[f + 1]; ++j) { + const double a = eq.a.Values()[j], + t = candidate.temperature[p.free[eq.a.ColumnIndices()[j]]]; + dot += a * t; + row_norm = Add(row_norm, Magnitude(a)); + scale = Add(scale, Mul(Magnitude(a), Magnitude(t))); + const I error{Error(eq.exact_a[j], a)}; + ea = Add(ea, error); + identity = Add(identity, Mul(Sub(eq.exact_a[j], I{a}), I{t})); + } + a_norm = Max(a_norm, row_norm); + t_norm = Max(t_norm, Magnitude(candidate.temperature[p.free[f]])); + red[f] = dot - eq.b[f]; + scale = Add(scale, Magnitude(eq.b[f])); + const double nu = + static_cast(eq.a.RowOffsets()[f + 1] - eq.a.RowOffsets()[f]); + red_error[f] = Roundoff(2 * nu + 2, scale).upper; + ev.free_residual_cert = std::max( + ev.free_residual_cert, Add(Magnitude(red[f]), I{red_error[f]}).upper); + ev.operator_construction_error = + std::max(ev.operator_construction_error, ea.upper); + const double eb = Error(eq.exact_b[f], eq.b[f]); + ev.rhs_construction_error = std::max(ev.rhs_construction_error, eb); + construction[f] = Sub(identity, Sub(eq.exact_b[f], I{eq.b[f]})); + } + const I denominator = Add(Mul(a_norm, t_norm), b_norm); + if (n && !Certificate(ev.free_residual_cert, denominator, ev.beta_limit, + ev.beta_cert)) + return Failure( + "Reduced residual is not certified on the reduced-system scale."); + + Vector change{m}, full_error{m}; + for (std::size_t i = 0; i < m; ++i) + change[i] = candidate.temperature[i] - old[i]; + if (!Finite(change)) return Failure("Nonfinite temperature change."); + I c_change, k_new, q_abs, c_temperatures, underflow_sum; + std::vector k_columns(m); + double energy = 0; + for (std::size_t i = 0; i < m; ++i) { + double cd = 0, kt = 0, ksum = 0; + I cs, ks, ts, cabs, nullscale, nullerror, coefficient; + for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) { + const auto col = s.c.ColumnIndices()[j]; + const double c = s.c.Values()[j], k = s.k.Values()[j], + t = candidate.temperature[col]; + cd += c * change[col]; + kt += k * t; + ksum += k; + cs = Add(cs, Mul(Magnitude(c), Magnitude(change[col]))); + ks = Add(ks, Mul(Magnitude(k), Magnitude(t))); + ts = Add(ts, Mul(Magnitude(c), Add(Magnitude(t), Magnitude(old[col])))); + cabs = Add(cabs, Magnitude(c)); + nullscale = Add(nullscale, Magnitude(k)); + k_columns[col] = Add(k_columns[col], I{k}); + if (s.k_error) { + nullerror = Add(nullerror, I{s.k_error->Values()[j]}); + coefficient = + Add(coefficient, Mul(I{s.k_error->Values()[j]}, Magnitude(t))); + } + if (s.c_error) + coefficient = + Add(coefficient, + Div(Mul(I{s.c_error->Values()[j]}, Magnitude(change[col])), + I{delta})); + } + energy += cd; + candidate.residual[i] = (cd / delta + kt) - s.q[i]; + const double nu = + static_cast(s.c.RowOffsets()[i + 1] - s.c.RowOffsets()[i]); + const I scale = Add(Div(cs, I{delta}), Add(ks, Magnitude(s.q[i]))); + // Propagate underflow in the capacity dot product through division, and + // underflow in the old/new subtraction through each capacity coefficient. + I under = Add(Div(Tiny(2 * nu + 2, cs), I{delta}), + Tiny(2 * nu + 4, Add(ks, Magnitude(s.q[i])))); + under = Add(under, Div(Mul(cabs, Tiny(1, ts)), I{delta})); + full_error[i] = Add(Add(Mul(Gamma(4 * nu + 6), scale), + Div(Mul(Gamma(1), ts), I{delta})), + under) + .upper; + if (s.q_error) coefficient = Add(coefficient, I{(*s.q_error)[i]}); + candidate.certified_residual[i] = + Add(Add(Magnitude(candidate.residual[i]), I{full_error[i]}), + coefficient) + .upper; + const I nullbound = Add(nullerror, Roundoff(nu, nullscale)); + if (!nullbound.IsFinite() || std::abs(ksum) > nullbound.upper) + return Failure("Conductivity does not certify its constant null mode."); + c_change = Add(c_change, cs); + k_new = Add(k_new, ks); + q_abs = Add(q_abs, Magnitude(s.q[i])); + c_temperatures = Add(c_temperatures, ts); + underflow_sum = Add(underflow_sum, under); + } + if (!Finite(candidate.residual) || !Finite(candidate.certified_residual) || + !std::isfinite(energy)) + return Failure("Nonfinite full residual or residual enclosure."); + for (std::size_t f = 0; f < n; ++f) { + const auto row = p.free[f]; + const I direct = + Add(I{candidate.residual[row]}, I{-full_error[row], full_error[row]}); + const I transferred = + Add(Add(I{red[f]}, I{-red_error[f], red_error[f]}), construction[f]); + const I intersection{std::max(direct.lower, transferred.lower), + std::min(direct.upper, transferred.upper)}; + if (!direct.IsFinite() || !transferred.IsFinite() || + !intersection.IsFinite()) + return Failure( + "Independent full residual enclosures have empty intersection."); + ev.full_free_residual_cert = + std::max(ev.full_free_residual_cert, Abs(intersection).upper); + const I bound = + Add(Abs(construction[f]), Add(I{red_error[f]}, I{full_error[row]})); + ev.construction_identity_bound = + std::max(ev.construction_identity_bound, bound.upper); + if (!bound.IsFinite() || + Abs(Sub(I{candidate.residual[row]}, I{red[f]})).lower > bound.upper) + return Failure( + "Full/reduced equation construction identity is not enclosed."); + } + if (n && !Certificate(ev.full_free_residual_cert, denominator, ev.beta_limit, + ev.full_beta_cert)) + return Failure( + "Full free residual is not certified on the reduced-system scale."); + + double qsum = 0, rsum = 0; + I rabs; + for (std::size_t i = 0; i < m; ++i) qsum += s.q[i]; + for (std::size_t i = 0; i < p.constrained.size(); ++i) { + candidate.reaction[i] = candidate.residual[p.constrained[i]]; + rsum += candidate.reaction[i]; + rabs = Add(rabs, Magnitude(candidate.reaction[i])); + } + I column_norm; + for (const auto x : k_columns) column_norm = Add(column_norm, Abs(x)); + const I se = Add(c_change, Mul(I{delta}, Add(q_abs, rabs))); + I be = Mul(Mul(Mul(I{delta}, I{nd}), beta_limit), denominator); + be = Add(be, Mul(Mul(I{delta}, column_norm), Norm(candidate.temperature))); + be = Add(be, Roundoff(4 * md + 12, se)); + be = Add(be, Mul(Mul(I{delta}, Gamma(2 * md + 6)), + Add(Div(c_change, I{delta}), Add(k_new, q_abs)))); + be = Add(be, Mul(Gamma(1), c_temperatures)); + be = Add(be, Mul(I{delta}, underflow_sum)); + ev.energy_change = energy; + ev.energy_defect = std::abs(energy - delta * (qsum + rsum)); + ev.energy_bound = be.upper; + if (!be.IsFinite() || !std::isfinite(ev.energy_defect) || + ev.energy_defect > be.upper) + return Failure( + "Physical heat balance failed independently of the free residual."); + + if (p.constrained.empty() && Zero(q_abs)) { + ev.dissipation_checked = true; + double cn = 0, co = 0, cd = 0, kn = 0; + I sc, subscale, t_one, d_square, o_square, skew_square, quadratic_underflow; + for (std::size_t i = 0; i < m; ++i) { + const double t = candidate.temperature[i]; + double row_cn = 0, row_co = 0, row_cd = 0, row_kn = 0; + I scale_cn, scale_co, scale_cd, scale_kn; + t_one = Add(t_one, Magnitude(t)); + d_square = Add(d_square, Mul(Magnitude(change[i]), Magnitude(change[i]))); + o_square = Add(o_square, Mul(Magnitude(old[i]), Magnitude(old[i]))); + for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) { + const auto col = s.c.ColumnIndices()[j]; + const double c = s.c.Values()[j], k = s.k.Values()[j]; + row_cn += c * candidate.temperature[col]; + row_co += c * old[col]; + row_cd += c * change[col]; + row_kn += k * candidate.temperature[col]; + scale_cn = Add( + scale_cn, Mul(Magnitude(c), Magnitude(candidate.temperature[col]))); + scale_co = Add(scale_co, Mul(Magnitude(c), Magnitude(old[col]))); + scale_cd = Add(scale_cd, Mul(Magnitude(c), Magnitude(change[col]))); + scale_kn = Add( + scale_kn, Mul(Magnitude(k), Magnitude(candidate.temperature[col]))); + const I nc = Mul(Mul(Magnitude(t), Magnitude(c)), + Magnitude(candidate.temperature[col])); + const I oc = + Mul(Mul(Magnitude(old[i]), Magnitude(c)), Magnitude(old[col])); + const I dc = Mul(Mul(Magnitude(change[i]), Magnitude(c)), + Magnitude(change[col])); + sc = Add(sc, Mul(I{.5}, Add(nc, Add(oc, dc)))); + sc = Add(sc, Mul(I{delta}, Mul(Mul(Magnitude(t), Magnitude(k)), + Magnitude(candidate.temperature[col])))); + subscale = + Add(subscale, + Mul(Mul(Add(Magnitude(t), Magnitude(old[i])), Magnitude(c)), + Magnitude(change[col]))); + const auto begin = s.c.ColumnIndices().begin() + + static_cast(s.c.RowOffsets()[col]); + const auto end = s.c.ColumnIndices().begin() + + static_cast(s.c.RowOffsets()[col + 1]); + const auto at = std::lower_bound(begin, end, i); + const double ct = (at != end && *at == i) + ? s.c.Values()[static_cast( + at - s.c.ColumnIndices().begin())] + : 0; + const I skew = c == ct ? I{0} : Abs(Sub(I{c}, I{ct})); + skew_square = Add(skew_square, Mul(skew, skew)); + } + cn += t * row_cn; + co += old[i] * row_co; + cd += change[i] * row_cd; + kn += t * row_kn; + const double operations = + 2 * + static_cast(s.c.RowOffsets()[i + 1] - s.c.RowOffsets()[i]); + I c_under = + Add(Mul(Magnitude(t), Tiny(operations, scale_cn)), + Add(Mul(Magnitude(old[i]), Tiny(operations, scale_co)), + Mul(Magnitude(change[i]), Tiny(operations, scale_cd)))); + c_under = Add(c_under, Tiny(6, Add(scale_cn, Add(scale_co, scale_cd)))); + const I k_under = + Add(Mul(Magnitude(t), Tiny(operations, scale_kn)), Tiny(2, scale_kn)); + quadratic_underflow = + Add(quadratic_underflow, + Add(Mul(I{.5}, c_under), Mul(I{delta}, k_under))); + } + I bc = Mul(Mul(Mul(I{delta}, t_one), beta_limit), denominator); + bc = Add(bc, Roundoff(6 * md + 20, sc)); + // Frobenius bounds the requested spectral skew norm without dense storage. + if (!Zero(skew_square)) + bc = Add(bc, Mul(I{.5}, Mul(thermal_internal::Sqrt(d_square), + Mul(thermal_internal::Sqrt(skew_square), + thermal_internal::Sqrt(o_square))))); + bc = Add(bc, Mul(Gamma(1), subscale)); + bc = Add(bc, quadratic_underflow); + ev.dissipation_defect = std::abs(.5 * ((cn - co) + cd) + delta * kn); + ev.dissipation_bound = bc.upper; + if (!bc.IsFinite() || !std::isfinite(ev.dissipation_defect) || + ev.dissipation_defect > bc.upper || !std::isfinite(cn) || + !std::isfinite(co) || Sub(I{cn}, I{co}).lower > Mul(I{2}, bc).upper) + return Failure( + "Insulated backward Euler capacity-norm dissipation failed."); + } + return Status::Ok(); +} +} // namespace + +Result SolveThermalIncrement( + const ThermalSystem& s, const ThermalPartition& p, const Vector& old, + double delta, LinearSolver& solver, ThermalFactorization& cache) { + auto status = Validate(s, p, old, delta); + if (!status.IsOk()) return Result::Failure(status); + auto equations = Build(s, p, old, delta); + if (!equations.HasValue()) + return Result::Failure(equations.GetStatus()); + const auto n = p.free.size(), m = old.Size(); + ThermalStepCandidate candidate{ + Vector{m}, Vector{m}, Vector{p.constrained.size()}, Vector{m}, {}}; + for (std::size_t i = 0; i < p.constrained.size(); ++i) + candidate.temperature[p.constrained[i]] = p.prescribed[i]; + if (n) { + if (!cache.matrix || cache.k != &s.k || cache.c != &s.c || + cache.solver != &solver || cache.free != p.free || + cache.delta != delta) { + cache = ThermalFactorization{}; + cache.matrix = equations.Value().a; + status = solver.Factorize(*cache.matrix); + if (!status.IsOk()) { + cache = ThermalFactorization{}; + return Result::Failure(status); + } + cache.k = &s.k; + cache.c = &s.c; + cache.solver = &solver; + cache.free = p.free; + cache.delta = delta; + } + Vector free{n}; + status = solver.Solve(equations.Value().b, free); + if (!status.IsOk()) return Result::Failure(status); + if (free.Size() != n || !Finite(free)) + return Result::Failure( + Failure("Solver returned invalid temperatures.")); + for (std::size_t i = 0; i < n; ++i) + candidate.temperature[p.free[i]] = free[i]; + } + status = Check(s, p, old, delta, equations.Value(), candidate); + if (!status.IsOk()) return Result::Failure(status); + return Result::Success(std::move(candidate)); +} + +Status RunThermalHistory(const ThermalAssembly& assembly, + const DofManager& dofs, ThermalAnalysisState& state, + LinearSolver& solver, const ThermalRecovery& recovery, + std::vector* evidence) { + if (!dofs.ValidateInvariants().IsOk() || + dofs.Profile() != DofProfile::kThermal || + dofs.FullDofCount() != state.CommittedState().temperature.Size() || + assembly.ElementOrder() != state.ElementOrder() || !recovery) + return Failure("History does not match the thermal operators and mapping."); + ThermalSystem system{assembly.Conductivity(), assembly.Capacity(), + assembly.HeatInput(), &assembly.ConductivityError(), + &assembly.CapacityError(), &assembly.HeatInputError()}; + ThermalPartition partition{dofs.FreeDofs(), dofs.ConstrainedDofs(), + dofs.PrescribedValues()}; + ThermalFactorization cache; + const std::uint64_t count = state.TimeGrid().Count(); + for (std::uint64_t i = + static_cast(state.AcceptedFrames().size()); + i <= count; ++i) { + const auto frame = state.TimeGrid().Frame(static_cast(i)); + if (!frame.HasValue()) return frame.GetStatus(); + auto candidate = SolveThermalIncrement(system, partition, + state.CommittedState().temperature, + frame.Value().delta, solver, cache); + if (!candidate.HasValue()) return candidate.GetStatus(); + auto flux = recovery(candidate.Value().temperature); + if (!flux.HasValue()) return flux.GetStatus(); + auto completed = state.MakeCompletedFrame( + static_cast(i), std::move(candidate.Value().temperature), + std::move(candidate.Value().reaction), std::move(flux.Value())); + if (!completed.HasValue()) return completed.GetStatus(); + const auto committed = state.CommitAcceptedFrame(completed.Value()); + if (!committed.IsOk()) return committed; + if (evidence) evidence->push_back(candidate.Value().evidence); + } + return Status::Ok(); +} +} // namespace fesa::thermal_transient_internal diff --git a/src/fesa/analysis/thermal_transient_algebra.h b/src/fesa/analysis/thermal_transient_algebra.h new file mode 100644 index 0000000..5b05ce7 --- /dev/null +++ b/src/fesa/analysis/thermal_transient_algebra.h @@ -0,0 +1,105 @@ +#ifndef FESA_ANALYSIS_THERMAL_TRANSIENT_ALGEBRA_H_ +#define FESA_ANALYSIS_THERMAL_TRANSIENT_ALGEBRA_H_ + +#include +#include +#include + +#include "fesa/analysis/thermal_analysis_state.h" +#include "fesa/solvers/linear/linear_solver.h" + +namespace fesa::thermal_transient_internal { + +/// @brief References fixed operators for one algebra/history invocation. +/// @note All referenced objects outlive calls; optional errors bound +/// coefficients, not quadrature/model error. K and C use the same full CSR +/// pattern. +struct ThermalSystem { + const SparseMatrix& k; + const SparseMatrix& c; + const Vector& q; + const SparseMatrix* k_error = nullptr; + const SparseMatrix* c_error = nullptr; + const Vector* q_error = nullptr; +}; + +/// @brief Declares an explicit generic partition, separate from public BC +/// syntax. +struct ThermalPartition { + std::vector free; + std::vector constrained; + Vector prescribed{0U}; +}; + +/// @brief Owns only the matrix backing the solver's current reusable +/// factorization. +struct ThermalFactorization { + std::optional matrix; + const SparseMatrix* k = nullptr; + const SparseMatrix* c = nullptr; + const LinearSolver* solver = nullptr; + std::vector free; + double delta = 0; +}; + +/// @brief Records internal scalar numerical evidence, never mandatory output +/// fields. +struct ThermalStepEvidence { + // Stored reduced system certificate; retain observed beta for forward + // budgets. + double beta_cert = 0; + double beta_limit = 0; + double free_residual_cert = 0; + double full_free_residual_cert = 0; + double full_beta_cert = 0; + // Stored-data arithmetic construction, separate from assembly E_C/E_K. + double operator_construction_error = 0; + double rhs_construction_error = 0; + double construction_identity_bound = 0; + double energy_change = 0; + double energy_defect = 0; + double energy_bound = 0; + bool dissipation_checked = false; + double dissipation_defect = 0; + double dissipation_bound = 0; +}; + +/// @brief Holds an unpublished finite/equilibrium-checked algebra candidate. +struct ThermalStepCandidate { + Vector temperature{0U}; + Vector residual{0U}; + Vector reaction{0U}; + // Entrywise |full r| + evaluation error + E_C|dT|/delta + E_K|T| + E_Q. + // Forward verification evidence never relaxes solve acceptance. + Vector certified_residual{0U}; + ThermalStepEvidence evidence; +}; + +/// @brief Solves/checks one actual-delta increment without modifying old state. +/// @note An empty free partition skips factorization and solve. Cache reuse is +/// restricted to identical operator addresses, partition, solver and actual +/// delta. +Result SolveThermalIncrement( + const ThermalSystem& system, const ThermalPartition& partition, + const Vector& old_temperature, double delta, LinearSolver& solver, + ThermalFactorization& factorization); + +/// @brief Provides an internal recovery seam for the actual history commit +/// path. +using ThermalRecovery = std::function>( + const Vector& temperature)>; + +/// @brief Advances an existing raw/accepted history, preserving it on every +/// failure. +/// @note Recovery runs after algebra validation and before atomic commit. +/// Optional evidence contains successful increments only; no count-sized +/// preallocation. +/// @pre State, assembly and dofs were created from the same fixed model. The +/// assembly's original Domain remains alive and unmoved throughout this call. +Status RunThermalHistory(const ThermalAssembly& assembly, + const DofManager& dofs, ThermalAnalysisState& state, + LinearSolver& solver, const ThermalRecovery& recovery, + std::vector* evidence = nullptr); + +} // namespace fesa::thermal_transient_internal +#endif // FESA_ANALYSIS_THERMAL_TRANSIENT_ALGEBRA_H_ diff --git a/src/fesa/analysis/thermal_transient_bounds.h b/src/fesa/analysis/thermal_transient_bounds.h new file mode 100644 index 0000000..0a08de0 --- /dev/null +++ b/src/fesa/analysis/thermal_transient_bounds.h @@ -0,0 +1,90 @@ +#ifndef FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_ +#define FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_ + +#include +#include +#include + +#include "elements/thermal_interval.h" + +namespace fesa::thermal_transient_internal::bounds { +using I = thermal_internal::Interval; +/// @brief Recognizes proven zero, distinct from an interval containing zero. +inline bool Zero(I x) { return x.lower == 0 && x.upper == 0; } +/// @brief Encloses addition, retaining proven exact sums and zero identities. +inline I Add(I a, I b) { + if (Zero(a)) return b; + if (Zero(b)) return a; + if (a.lower == a.upper && b.lower == b.upper) { + const double sum = a.lower + b.lower; + const double bv = sum - a.lower, av = sum - bv; + const double error = (a.lower - av) + (b.lower - bv); + // TwoSum proves exact additions (including exact cancellation) without + // creating spurious construction uncertainty at prescribed stationary DOFs. + if (std::isfinite(sum) && std::isfinite(bv) && std::isfinite(av) && + error == 0) + return I{sum}; + } + return a + b; +} +/// @brief Recognizes singleton positive powers of two for exact scaling. +inline bool PowerOfTwo(I x, int& exponent) { + if (x.lower != x.upper || x.lower <= 0 || !std::isfinite(x.lower)) + return false; + const double fraction = std::frexp(x.lower, &exponent); + --exponent; + return fraction == .5; +} +/// @brief Encloses products, including loss in subnormal power-of-two scaling. +inline I Mul(I a, I b) { + if (Zero(a) || Zero(b)) return I{0}; + if (a.lower == a.upper && a.lower < 0) { + const I value = Mul(I{-a.lower}, b); + return {-value.upper, -value.lower}; + } + if (b.lower == b.upper && b.lower < 0) { + const I value = Mul(a, I{-b.lower}); + return {-value.upper, -value.lower}; + } + int exponent = 0; + if (PowerOfTwo(a, exponent)) return thermal_internal::Scale(b, exponent); + if (PowerOfTwo(b, exponent)) return thermal_internal::Scale(a, exponent); + return a * b; +} +/// @brief Encloses division without discarding subnormal quotient uncertainty. +inline I Div(I a, I b) { + if (b.ContainsZero()) return a / b; + if (Zero(a)) return I{0}; + int exponent = 0; + if (PowerOfTwo(b, exponent)) return thermal_internal::Scale(a, -exponent); + return a / b; +} +/// @brief Encloses the absolute value of a signed interval. +inline I Abs(I a) { + return {a.ContainsZero() ? 0 : std::min(std::abs(a.lower), std::abs(a.upper)), + std::max(std::abs(a.lower), std::abs(a.upper))}; +} +/// @brief Changes interval sign exactly. +inline I Neg(I a) { return {-a.upper, -a.lower}; } +/// @brief Preserves exact cancellation through the addition certificate. +inline I Sub(I a, I b) { return Add(a, Neg(b)); } +/// @brief Lifts an exact binary64 magnitude into an interval. +inline I Magnitude(double a) { return I{std::abs(a)}; } +/// @brief Encloses the approved binary64 gamma operation-count bound. +inline I Gamma(double count) { + const I product = Mul(I{count}, I{0x1p-53}); + return Div(product, Sub(I{1}, product)); +} +/// @brief Bounds absolute underflow by one minimum subnormal per operation. +/// @note A proven zero expression performs no inexact operation. +inline I Tiny(double count, I scale) { + return Zero(scale) + ? I{0} + : Mul(I{count}, I{std::numeric_limits::denorm_min()}); +} +/// @brief Combines relative gamma error with absolute underflow error. +inline I Roundoff(double count, I scale) { + return Add(Mul(Gamma(count), scale), Tiny(count, scale)); +} +} // namespace fesa::thermal_transient_internal::bounds +#endif // FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_ diff --git a/src/fesa/analysis/transient_heat_procedure.cpp b/src/fesa/analysis/transient_heat_procedure.cpp new file mode 100644 index 0000000..b8cfdf1 --- /dev/null +++ b/src/fesa/analysis/transient_heat_procedure.cpp @@ -0,0 +1,39 @@ +#include "fesa/analysis/transient_heat_procedure.h" + +#include + +#include "thermal_transient_algebra.h" + +namespace fesa { + +TransientHeatProcedureExecutor::TransientHeatProcedureExecutor( + const ParallelFor& parallel_for, LinearSolver& linear_solver) + : parallel_for_{parallel_for}, linear_solver_{linear_solver} {} + +Result TransientHeatProcedureExecutor::Execute( + const Domain& domain) const { + auto model = AnalysisModel::Create(domain); + if (!model.HasValue()) + return Result::Failure(model.GetStatus()); + auto dofs = DofManager::CreateThermal(model.Value()); + if (!dofs.HasValue()) + return Result::Failure(dofs.GetStatus()); + auto assembly = AssembleThermal(model.Value(), dofs.Value(), parallel_for_); + if (!assembly.HasValue()) + return Result::Failure(assembly.GetStatus()); + auto state = ThermalAnalysisState::Create(model.Value(), dofs.Value(), + assembly.Value()); + if (!state.HasValue()) + return Result::Failure(state.GetStatus()); + const Status completed = thermal_transient_internal::RunThermalHistory( + assembly.Value(), dofs.Value(), state.Value(), linear_solver_, + [&](const Vector& temperature) { + return assembly.Value().RecoverFlux(temperature, parallel_for_); + }); + if (!completed.IsOk()) + return Result::Failure(completed); + return Result::Success( + {std::move(state.Value()), domain.Warnings()}); +} + +} // namespace fesa diff --git a/tests/unit/analysis/transient_heat_procedure_test.cpp b/tests/unit/analysis/transient_heat_procedure_test.cpp new file mode 100644 index 0000000..b31af42 --- /dev/null +++ b/tests/unit/analysis/transient_heat_procedure_test.cpp @@ -0,0 +1,472 @@ +#include "fesa/analysis/transient_heat_procedure.h" + +#include + +#include +#include +#include +#include +#include + +#include "analysis/thermal_transient_algebra.h" +#include "analysis/thermal_transient_bounds.h" +#include "fesa/solvers/linear/mkl_pardiso_solver.h" + +namespace { +namespace ti = fesa::thermal_transient_internal; +std::string Messages(const fesa::Status& status) { + std::ostringstream out; + for (const auto& d : status.Diagnostics()) out << d.message << "\n"; + return out.str(); +} +fesa::ModelDefinition Plate(double initial = 0, double period = .9, + double dt = .3) { + fesa::ModelDefinition m{}; + m.nodes = {{{"plate", 1, "1"}, {0, 0, 0}, {}}, + {{"plate", 2, "2"}, {1, 0, 0}, {}}, + {{"plate", 3, "3"}, {1, 1, 0}, {}}, + {{"plate", 4, "4"}, {0, 1, 0}, {}}}; + m.thermal_materials = {{{"plate", 1, "1"}, "heat", 1, 1, 1, {}}}; + m.thermal_sections = {{{"plate", 1, "1"}, "section", 2, 0, {}}}; + m.thermal_elements = {{{"plate", 1, "1"}, {{0, 1, 2, 3}}, 0, {}}}; + for (fesa::EntityIndex i = 0; i < 4; ++i) + m.thermal_initial_conditions.push_back({i, initial, {}}); + m.thermal_step = fesa::ThermalStepDefinition{"heat", dt, period, 100}; + return m; +} +fesa::SparseMatrix TwoMatrix(double a, double b) { + return fesa::SparseMatrix::FromCoo(2, 2, + {{0, 0, a, 0, 0}, + {0, 1, b, 0, 1}, + {1, 0, b, 0, 2}, + {1, 1, a, 0, 3}}, + {{0, 2, 4}, {0, 1, 0, 1}}) + .Value(); +} +class CountingSolver : public fesa::LinearSolver { + public: + fesa::Status Factorize(const fesa::SparseMatrix& a) override { + ++factorizations; + return solver.Factorize(a); + } + fesa::Status Solve(const fesa::Vector& b, fesa::Vector& x) const override { + ++solves; + if (before_solve) before_solve(solves); + if (solves == fail_on) + return fesa::Status::Failure(fesa::FailureCategory::kSolver, {}); + const auto status = solver.Solve(b, x); + if (status.IsOk() && solves == corrupt_on) x[0] = corrupt_value; + return status; + } + int factorizations = 0; + mutable int solves = 0; + int fail_on = 0, corrupt_on = 0; + double corrupt_value = 0; + std::function before_solve; + fesa::MklPardisoSolver solver; +}; +std::vector Snapshot(const fesa::ThermalAnalysisState& state) { + std::vector result; + const auto scalar = [&](double x) { + std::uint64_t bits = 0; + std::memcpy(&bits, &x, sizeof(bits)); + result.push_back(bits); + }; + for (const auto& frame : state.AcceptedFrames()) { + result.insert(result.end(), + {frame.frame_index, frame.increment_index, frame.is_initial, + frame.completed.has_value()}); + scalar(frame.time); + scalar(frame.time_increment); + for (const auto node : frame.node_order) result.push_back(node); + for (std::size_t i = 0; i < frame.temperature.Size(); ++i) + scalar(frame.temperature[i]); + if (frame.completed) { + for (const auto dof : frame.completed->reaction_dofs) { + result.push_back(dof.node); + result.push_back(dof.source_dof); + } + for (std::size_t i = 0; i < frame.completed->reaction_heat_rate.Size(); + ++i) + scalar(frame.completed->reaction_heat_rate[i]); + for (const auto e : frame.completed->element_order) result.push_back(e); + for (const auto& e : frame.completed->global_heat_flux) + for (const auto& q : e) + for (const auto x : q.Components()) scalar(x); + } + } + return result; +} +} // namespace + +TEST(TransientHeatProcedure, RawConstrainedCapacityProducesSignedTwoDofJump) { + const auto c = TwoMatrix(2, 1), k = TwoMatrix(1, -1); + const fesa::Vector q{2}; + ti::ThermalPartition p{{0}, {1}, fesa::Vector{1, 1}}; + CountingSolver solver; + ti::ThermalFactorization cache; + auto first = ti::SolveThermalIncrement({k, c, q}, p, fesa::Vector{2}, .5, + solver, cache); + ASSERT_TRUE(first.HasValue()) << Messages(first.GetStatus()); + EXPECT_DOUBLE_EQ(first.Value().temperature[0], -.2); + EXPECT_DOUBLE_EQ(first.Value().reaction[0], 4.8); + auto second = ti::SolveThermalIncrement( + {k, c, q}, p, first.Value().temperature, .5, solver, cache); + ASSERT_TRUE(second.HasValue()) << Messages(second.GetStatus()); + EXPECT_NEAR(second.Value().temperature[0], .04, + 16 * std::numeric_limits::epsilon() * .04); + EXPECT_DOUBLE_EQ(second.Value().reaction[0], 1.44); + EXPECT_EQ(solver.factorizations, 1); + EXPECT_EQ(solver.solves, 2); +} + +TEST(TransientHeatProcedure, AllConstrainedBypassesZeroSizeFactorization) { + const auto c = TwoMatrix(2, 1), k = TwoMatrix(1, -1); + const fesa::Vector q{2}; + ti::ThermalPartition p{{}, {0, 1}, fesa::Vector{2, 1}}; + CountingSolver solver; + ti::ThermalFactorization cache; + auto result = ti::SolveThermalIncrement({k, c, q}, p, fesa::Vector{2}, .5, + solver, cache); + ASSERT_TRUE(result.HasValue()) << Messages(result.GetStatus()); + EXPECT_EQ(result.Value().reaction[0], 6); + EXPECT_EQ(result.Value().reaction[1], 6); + EXPECT_EQ(solver.factorizations, 0); + EXPECT_EQ(solver.solves, 0); +} + +TEST(TransientHeatProcedure, + CancellationCannotUseFullSystemScaleToPassCertificate) { + const auto c = TwoMatrix(2, 1), k = TwoMatrix(1, -1); + const fesa::Vector q{2}; + const double large = std::ldexp(1., 56); + ti::ThermalPartition p{{0}, {1}, fesa::Vector{1, 2 * large}}; + fesa::Vector old{2}; + old[0] = 1; + old[1] = large; + CountingSolver solver; + ti::ThermalFactorization cache; + auto valid = ti::SolveThermalIncrement({k, c, q}, p, old, .5, solver, cache); + ASSERT_TRUE(valid.HasValue()) << Messages(valid.GetStatus()); + EXPECT_DOUBLE_EQ(valid.Value().temperature[0], .8); + // The stable old-minus-g RHS retains b=4. Inject the naive lost-b solution; + // it must fail, even though normalization by the huge full state could pass. + solver.corrupt_on = 2; + EXPECT_FALSE(ti::SolveThermalIncrement({k, c, q}, p, old, .5, solver, cache) + .HasValue()); +} + +TEST(TransientHeatProcedure, + ExactZeroAndUniformInsulatedStatesAndActualDeltaCache) { + struct Schedule { + double period, dt; + std::size_t steps; + int factors; + double last; + }; + const Schedule schedules[] = {{.9, .3, 3, 1, .3}, + {std::nextafter(.9, 1.), .3, 4, 2, + std::fma(-3., .3, std::nextafter(.9, 1.))}, + {1., .3, 4, 2, std::fma(-3., .3, 1.)}, + {.2, .3, 1, 1, .2}}; + for (const auto schedule : schedules) + for (const double initial : {0., 2.}) { + auto domain = + fesa::Domain::Create(Plate(initial, schedule.period, schedule.dt)); + ASSERT_TRUE(domain.HasValue()) << Messages(domain.GetStatus()); + CountingSolver solver; + fesa::SerialParallelFor parallel; + auto result = + fesa::TransientHeatProcedureExecutor{parallel, solver}.Execute( + domain.Value()); + ASSERT_TRUE(result.HasValue()) << Messages(result.GetStatus()); + const auto& state = result.Value().state; + EXPECT_TRUE(state.IsComplete()); + EXPECT_EQ(state.AcceptedFrames().size(), schedule.steps + 1); + EXPECT_EQ(state.CommittedState().time, schedule.period); + EXPECT_EQ(state.CommittedState().time_increment, schedule.last); + EXPECT_EQ(solver.factorizations, schedule.factors); + EXPECT_EQ(solver.solves, static_cast(schedule.steps)); + EXPECT_FALSE(state.AcceptedFrames()[0].completed.has_value()); + for (const auto& frame : state.AcceptedFrames()) { + for (std::size_t i = 0; i < 20; ++i) { + if (initial == 0) + EXPECT_EQ(frame.temperature[i], 0); + else + EXPECT_NEAR(frame.temperature[i], initial, + 4096 * 0x1p-53 * initial); + } + if (frame.completed) + EXPECT_EQ(frame.completed->reaction_heat_rate.Size(), 0U); + } + } +} + +TEST(TransientHeatProcedure, + MklMatchesIndependentAffineDiscreteModeAndDissipation) { + // Q1 unit-square mass/stiffness give lambda=12 for (1,-1,-1,1). + // Constant thickness is its exact P2 null mode, independently of production. + auto definition = Plate(0, .25, .25); + const double mode[] = {1, -1, -1, 1}; + for (fesa::EntityIndex i = 0; i < 4; ++i) + definition.thermal_initial_conditions[i] = {i, mode[i], {}}; + auto domain = fesa::Domain::Create(std::move(definition)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + fesa::SerialParallelFor parallel; + auto assembly = fesa::AssembleThermal(model.Value(), dofs.Value(), parallel); + ASSERT_TRUE(assembly.HasValue()); + auto state = fesa::ThermalAnalysisState::Create(model.Value(), dofs.Value(), + assembly.Value()); + ASSERT_TRUE(state.HasValue()); + CountingSolver solver; + ti::ThermalFactorization cache; + const auto& a = assembly.Value(); + auto result = ti::SolveThermalIncrement( + {a.Conductivity(), a.Capacity(), a.HeatInput(), &a.ConductivityError(), + &a.CapacityError(), &a.HeatInputError()}, + {dofs.Value().FreeDofs(), {}, fesa::Vector{0}}, + state.Value().CommittedState().temperature, .25, solver, cache); + ASSERT_TRUE(result.HasValue()) << Messages(result.GetStatus()); + const auto& candidate = result.Value(); + // lambda_min(Msurface)>=1/36 and lambda_min(Z0)>=1/30 by exact + // tensor eigenvalues/Gershgorin; h/2=1. Thus ||A^-1||inf<=sqrt(20)*1080*dt. + double residual = 0; + for (std::size_t i = 0; i < 20; ++i) + residual = std::max(residual, candidate.certified_residual[i]); + const double bound = std::sqrt(20.) * 1080 * .25 * residual; + for (std::size_t i = 0; i < 20; ++i) + EXPECT_NEAR(candidate.temperature[i], mode[i / 5] / 4, bound); + EXPECT_LE(candidate.evidence.beta_cert, candidate.evidence.beta_limit); + EXPECT_LE(candidate.evidence.energy_defect, candidate.evidence.energy_bound); + EXPECT_TRUE(candidate.evidence.dissipation_checked); + EXPECT_LE(candidate.evidence.dissipation_defect, + candidate.evidence.dissipation_bound); + auto flux = a.RecoverFlux(candidate.temperature, parallel); + ASSERT_TRUE(flux.HasValue()); + // q=(.5,0,0); each exact B row has absolute sum bounded by 16 for this h=2 + // square. + for (const auto& q : flux.Value()[0]) { + EXPECT_NEAR(q.X(), .5, 16 * bound + 4096 * 0x1p-53 * .5); + EXPECT_NEAR(q.Y(), 0, 16 * bound); + EXPECT_NEAR(q.Z(), 0, 16 * bound); + } +} + +TEST(TransientHeatProcedure, + LaterSolveRecoveryAndPayloadFailuresPreserveEveryAcceptedBit) { + for (int failure = 0; failure < 5; ++failure) { + auto definition = Plate(2, .6, .3); + definition.thermal_model_boundaries = {{0, 11, 1, {}}}; + auto domain = fesa::Domain::Create(std::move(definition)); + ASSERT_TRUE(domain.HasValue()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + fesa::SerialParallelFor parallel; + auto assembly = + fesa::AssembleThermal(model.Value(), dofs.Value(), parallel); + ASSERT_TRUE(assembly.HasValue()); + auto state = fesa::ThermalAnalysisState::Create(model.Value(), dofs.Value(), + assembly.Value()); + ASSERT_TRUE(state.HasValue()); + CountingSolver solver; + if (failure == 0) solver.fail_on = 2; + if (failure >= 3) { + solver.corrupt_on = 2; + solver.corrupt_value = + failure == 3 ? std::numeric_limits::quiet_NaN() : 1e20; + } + std::vector before; + solver.before_solve = [&](int count) { + if (count == 2) before = Snapshot(state.Value()); + }; + int recoveries = 0; + std::vector evidence; + const auto status = ti::RunThermalHistory( + assembly.Value(), dofs.Value(), state.Value(), solver, + [&](const fesa::Vector& temperature) { + ++recoveries; + if (recoveries == 2 && failure == 1) + return fesa::Result>::Failure( + fesa::Status::Failure(fesa::FailureCategory::kSolver, {})); + auto flux = assembly.Value().RecoverFlux(temperature, parallel); + if (recoveries == 2 && failure == 2) + flux.Value()[0][0] = + fesa::Vector3{std::numeric_limits::infinity(), 0, 0}; + return flux; + }, + &evidence); + EXPECT_FALSE(status.IsOk()); + ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U) << Messages(status); + EXPECT_FALSE(state.Value().IsComplete()); + EXPECT_EQ(evidence.size(), 1U); + EXPECT_EQ(Snapshot(state.Value()), before); + EXPECT_EQ(state.Value().AcceptedFrames()[0].temperature[0], 2); + EXPECT_EQ(state.Value().CommittedState().temperature[0], 1); + EXPECT_EQ(solver.factorizations, 1); + } +} + +TEST(TransientHeatProcedure, + SignedPowerConstrainedReactionAndDisconnectedComponentBalance) { + auto definition = Plate(0, .5, .25); + // A separate translated patch keeps an independent no-BC heat balance. + for (fesa::EntityIndex i = 0; i < 4; ++i) { + auto node = definition.nodes[i]; + node.source_id = {"other", i + 1, std::to_string(i + 1)}; + node.coordinates[0] += 3; + definition.nodes.push_back(node); + definition.thermal_initial_conditions.push_back({i + 4, 0, {}}); + } + definition.thermal_elements.push_back( + {{"other", 2, "2"}, {{4, 5, 6, 7}}, 0, {}}); + definition.thermal_model_boundaries = {{0, 11, 0, {}}}; + definition.thermal_step->loads = { + {0, 11, 3, {}}, {1, 15, -3, {}}, {4, 11, 5, {}}, {6, 15, -5, {}}}; + auto domain = fesa::Domain::Create(std::move(definition)); + ASSERT_TRUE(domain.HasValue()) << Messages(domain.GetStatus()); + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::CreateThermal(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + fesa::SerialParallelFor parallel; + auto assembly = fesa::AssembleThermal(model.Value(), dofs.Value(), parallel); + ASSERT_TRUE(assembly.HasValue()); + auto state = fesa::ThermalAnalysisState::Create(model.Value(), dofs.Value(), + assembly.Value()); + ASSERT_TRUE(state.HasValue()); + CountingSolver solver; + std::vector evidence; + const auto status = ti::RunThermalHistory( + assembly.Value(), dofs.Value(), state.Value(), solver, + [&](const fesa::Vector& t) { + return assembly.Value().RecoverFlux(t, parallel); + }, + &evidence); + ASSERT_TRUE(status.IsOk()) << Messages(status); + ASSERT_EQ(evidence.size(), 2U); + const double weight[] = {1. / 24, 1. / 6, 1. / 12, 1. / 6, 1. / 24}; + for (std::size_t step = 1; step < 3; ++step) { + const auto& frame = state.Value().AcceptedFrames()[step]; + const auto& old = state.Value().AcceptedFrames()[step - 1]; + double balance[2] = {0, 0}, scale[2] = {0, 0}; + for (std::size_t i = 0; i < 40; ++i) { + const double term = + weight[i % 5] * (frame.temperature[i] - old.temperature[i]); + balance[i / 20] += term; + scale[i / 20] += std::abs(term); + } + const double reaction = frame.completed->reaction_heat_rate[0]; + EXPECT_NEAR(balance[0], .25 * reaction, + evidence[step - 1].energy_bound + 4096 * 0x1p-53 * scale[0]); + EXPECT_NEAR(balance[1], 0, + evidence[step - 1].energy_bound + 4096 * 0x1p-53 * scale[1]); + EXPECT_FALSE(evidence[step - 1].dissipation_checked); + } + // Adding a constrained power changes only its full-residual reaction by -P. + auto unloaded = Plate(0, .5, .25); + unloaded.thermal_model_boundaries = {{0, 11, 0, {}}}; + unloaded.thermal_step->loads = {{1, 15, -3, {}}}; + auto other_domain = fesa::Domain::Create(std::move(unloaded)); + ASSERT_TRUE(other_domain.HasValue()); + CountingSolver other_solver; + auto other = + fesa::TransientHeatProcedureExecutor{parallel, other_solver}.Execute( + other_domain.Value()); + ASSERT_TRUE(other.HasValue()); + for (std::size_t step = 1; step < 3; ++step) { + const auto& first = state.Value().AcceptedFrames()[step]; + const auto& second = other.Value().state.AcceptedFrames()[step]; + for (std::size_t i = 0; i < 20; ++i) + EXPECT_NEAR(first.temperature[i], second.temperature[i], + 4096 * 0x1p-53 * std::abs(second.temperature[i])); + EXPECT_NEAR(first.completed->reaction_heat_rate[0] - + second.completed->reaction_heat_rate[0], + -3., 4096 * 0x1p-53 * 3); + } +} + +TEST(TransientHeatProcedure, + ConstructionEnclosuresRetainExactIdentitiesAndSubnormalLoss) { + namespace b = ti::bounds; + const double tiny = std::numeric_limits::denorm_min(); + const auto exact = [&](b::I value, double expected) { + EXPECT_EQ(value.lower, expected); + EXPECT_EQ(value.upper, expected); + }; + exact(b::Add(b::I{0x1p56}, b::I{-0x1p56}), 0); + exact(b::Add(b::I{0x1p-1022}, b::I{-std::nextafter(0x1p-1022, 0.)}), tiny); + exact(b::Div(b::I{2}, b::I{.5}), 4); + exact(b::Mul(b::I{-1}, b::I{-0x1p-100}), 0x1p-100); + exact(b::Mul(b::I{0x1p-500}, b::I{0x1p400}), 0x1p-100); + exact(b::Add(b::I{-.8}, b::I{1}), 0x1.9999999999998p-3); + // Exact rational targets lie strictly between the listed adjacent binary64 + // numbers: 1+2^-53 and 2^-1075. Neither may become a singleton rounded value. + const auto midpoint = b::Add(b::I{1}, b::I{0x1p-53}); + EXPECT_LE(midpoint.lower, 1); + EXPECT_GE(midpoint.upper, std::nextafter(1., 2.)); + for (const auto half : + {b::Mul(b::I{tiny}, b::I{.5}), b::Div(b::I{tiny}, b::I{2})}) { + EXPECT_LE(half.lower, 0); + EXPECT_GE(half.upper, tiny); + EXPECT_NE(half.lower, half.upper); + } + EXPECT_FALSE( + b::Mul(b::I{std::numeric_limits::max()}, b::I{2}).IsFinite()); + exact(b::Roundoff(80, b::I{0}), 0); + EXPECT_GT(b::Roundoff(80, b::I{tiny}).upper, 0); +} + +TEST(TransientHeatProcedure, + GenericAlgebraRejectsInvalidInputsAndUncertifiedSubnormalSolve) { + const auto c = TwoMatrix(2, 1), k = TwoMatrix(1, -1); + const fesa::Vector q{2}; + const ti::ThermalPartition p{{0}, {1}, fesa::Vector{1, 1}}; + CountingSolver solver; + ti::ThermalFactorization cache; + for (const double delta : {0., -1., std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN()}) + EXPECT_FALSE(ti::SolveThermalIncrement({k, c, q}, p, fesa::Vector{2}, delta, + solver, cache) + .HasValue()); + EXPECT_FALSE(ti::SolveThermalIncrement({k, c, q}, + {{0}, {0}, fesa::Vector{1, 1}}, + fesa::Vector{2}, .5, solver, cache) + .HasValue()); + EXPECT_FALSE(ti::SolveThermalIncrement({k, c, q}, p, fesa::Vector{1}, .5, + solver, cache) + .HasValue()); + auto invalid = TwoMatrix(2, -1); + EXPECT_FALSE(ti::SolveThermalIncrement({invalid, c, q}, p, fesa::Vector{2}, + .5, solver, cache) + .HasValue()); + ti::ThermalPartition small{ + {0}, {1}, fesa::Vector{1, std::numeric_limits::denorm_min()}}; + EXPECT_FALSE(ti::SolveThermalIncrement({k, c, q}, small, fesa::Vector{2}, .5, + solver, cache) + .HasValue()); +} + +TEST(TransientHeatProcedure, + SerialAndTbbExecutionsHaveIdenticalAcceptedHistories) { + auto definition = Plate(2, .6, .3); + definition.thermal_step->loads = {{0, 11, 3, {}}, {2, 15, -3, {}}}; + definition.thermal_model_boundaries = {{1, 15, -1, {}}}; + auto domain = fesa::Domain::Create(std::move(definition)); + ASSERT_TRUE(domain.HasValue()); + fesa::SerialParallelFor serial; + fesa::TbbParallelFor parallel; + CountingSolver a, b; + auto left = + fesa::TransientHeatProcedureExecutor{serial, a}.Execute(domain.Value()); + auto right = + fesa::TransientHeatProcedureExecutor{parallel, b}.Execute(domain.Value()); + ASSERT_TRUE(left.HasValue()) << Messages(left.GetStatus()); + ASSERT_TRUE(right.HasValue()) << Messages(right.GetStatus()); + EXPECT_EQ(Snapshot(left.Value().state), Snapshot(right.Value().state)); +}