feat: integrate certified transient heat procedure
This commit is contained in:
KOKO\Mimi
@@ -0,0 +1,39 @@
|
||||
#ifndef FESA_ANALYSIS_TRANSIENT_HEAT_PROCEDURE_H_
|
||||
#define FESA_ANALYSIS_TRANSIENT_HEAT_PROCEDURE_H_
|
||||
|
||||
#include <vector>
|
||||
|
||||
#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<Diagnostic> 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<ThermalExecutionResult> Execute(const Domain& domain) const;
|
||||
|
||||
private:
|
||||
const ParallelFor& parallel_for_;
|
||||
LinearSolver& linear_solver_;
|
||||
};
|
||||
|
||||
} // namespace fesa
|
||||
#endif // FESA_ANALYSIS_TRANSIENT_HEAT_PROCEDURE_H_
|
||||
@@ -0,0 +1,503 @@
|
||||
#include "thermal_transient_algebra.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
|
||||
#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<bool> 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<I> exact_a;
|
||||
std::vector<I> exact_b;
|
||||
};
|
||||
/// @brief Forms the reduced operator and raw-old-state backward Euler RHS.
|
||||
Result<Equations> Build(const ThermalSystem& s, const ThermalPartition& p,
|
||||
const Vector& old, double delta) {
|
||||
const auto n = p.free.size(), m = old.Size();
|
||||
std::vector<std::size_t> 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<CooContribution> coo;
|
||||
Vector b{n};
|
||||
std::vector<I> 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<Equations>::Failure(a.GetStatus());
|
||||
if (!Finite(b))
|
||||
return Result<Equations>::Failure(Failure("Nonfinite backward Euler RHS."));
|
||||
return Result<Equations>::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<double>(m), nd = static_cast<double>(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<I> 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<double>(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<I> 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<double>(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<std::ptrdiff_t>(s.c.RowOffsets()[col]);
|
||||
const auto end = s.c.ColumnIndices().begin() +
|
||||
static_cast<std::ptrdiff_t>(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<std::size_t>(
|
||||
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<double>(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<ThermalStepCandidate> 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<ThermalStepCandidate>::Failure(status);
|
||||
auto equations = Build(s, p, old, delta);
|
||||
if (!equations.HasValue())
|
||||
return Result<ThermalStepCandidate>::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<ThermalStepCandidate>::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<ThermalStepCandidate>::Failure(status);
|
||||
if (free.Size() != n || !Finite(free))
|
||||
return Result<ThermalStepCandidate>::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<ThermalStepCandidate>::Failure(status);
|
||||
return Result<ThermalStepCandidate>::Success(std::move(candidate));
|
||||
}
|
||||
|
||||
Status RunThermalHistory(const ThermalAssembly& assembly,
|
||||
const DofManager& dofs, ThermalAnalysisState& state,
|
||||
LinearSolver& solver, const ThermalRecovery& recovery,
|
||||
std::vector<ThermalStepEvidence>* 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<std::uint64_t>(state.AcceptedFrames().size());
|
||||
i <= count; ++i) {
|
||||
const auto frame = state.TimeGrid().Frame(static_cast<std::uint32_t>(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<std::uint32_t>(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
|
||||
@@ -0,0 +1,105 @@
|
||||
#ifndef FESA_ANALYSIS_THERMAL_TRANSIENT_ALGEBRA_H_
|
||||
#define FESA_ANALYSIS_THERMAL_TRANSIENT_ALGEBRA_H_
|
||||
|
||||
#include <functional>
|
||||
#include <optional>
|
||||
#include <vector>
|
||||
|
||||
#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<std::size_t> free;
|
||||
std::vector<std::size_t> constrained;
|
||||
Vector prescribed{0U};
|
||||
};
|
||||
|
||||
/// @brief Owns only the matrix backing the solver's current reusable
|
||||
/// factorization.
|
||||
struct ThermalFactorization {
|
||||
std::optional<SparseMatrix> matrix;
|
||||
const SparseMatrix* k = nullptr;
|
||||
const SparseMatrix* c = nullptr;
|
||||
const LinearSolver* solver = nullptr;
|
||||
std::vector<std::size_t> 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<ThermalStepCandidate> 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<Result<std::vector<ThermalElementFlux>>(
|
||||
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<ThermalStepEvidence>* evidence = nullptr);
|
||||
|
||||
} // namespace fesa::thermal_transient_internal
|
||||
#endif // FESA_ANALYSIS_THERMAL_TRANSIENT_ALGEBRA_H_
|
||||
@@ -0,0 +1,90 @@
|
||||
#ifndef FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_
|
||||
#define FESA_ANALYSIS_THERMAL_TRANSIENT_BOUNDS_H_
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
|
||||
#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<double>::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_
|
||||
@@ -0,0 +1,39 @@
|
||||
#include "fesa/analysis/transient_heat_procedure.h"
|
||||
|
||||
#include <utility>
|
||||
|
||||
#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<ThermalExecutionResult> TransientHeatProcedureExecutor::Execute(
|
||||
const Domain& domain) const {
|
||||
auto model = AnalysisModel::Create(domain);
|
||||
if (!model.HasValue())
|
||||
return Result<ThermalExecutionResult>::Failure(model.GetStatus());
|
||||
auto dofs = DofManager::CreateThermal(model.Value());
|
||||
if (!dofs.HasValue())
|
||||
return Result<ThermalExecutionResult>::Failure(dofs.GetStatus());
|
||||
auto assembly = AssembleThermal(model.Value(), dofs.Value(), parallel_for_);
|
||||
if (!assembly.HasValue())
|
||||
return Result<ThermalExecutionResult>::Failure(assembly.GetStatus());
|
||||
auto state = ThermalAnalysisState::Create(model.Value(), dofs.Value(),
|
||||
assembly.Value());
|
||||
if (!state.HasValue())
|
||||
return Result<ThermalExecutionResult>::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<ThermalExecutionResult>::Failure(completed);
|
||||
return Result<ThermalExecutionResult>::Success(
|
||||
{std::move(state.Value()), domain.Warnings()});
|
||||
}
|
||||
|
||||
} // namespace fesa
|
||||
@@ -0,0 +1,472 @@
|
||||
#include "fesa/analysis/transient_heat_procedure.h"
|
||||
|
||||
#include <gtest/gtest.h>
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <functional>
|
||||
#include <limits>
|
||||
#include <sstream>
|
||||
|
||||
#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<void(int)> before_solve;
|
||||
fesa::MklPardisoSolver solver;
|
||||
};
|
||||
std::vector<std::uint64_t> Snapshot(const fesa::ThermalAnalysisState& state) {
|
||||
std::vector<std::uint64_t> 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<double>::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<int>(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<double>::quiet_NaN() : 1e20;
|
||||
}
|
||||
std::vector<std::uint64_t> before;
|
||||
solver.before_solve = [&](int count) {
|
||||
if (count == 2) before = Snapshot(state.Value());
|
||||
};
|
||||
int recoveries = 0;
|
||||
std::vector<ti::ThermalStepEvidence> 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<std::vector<fesa::ThermalElementFlux>>::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<double>::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<ti::ThermalStepEvidence> 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<double>::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<double>::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<double>::infinity(),
|
||||
std::numeric_limits<double>::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<double>::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));
|
||||
}
|
||||
Reference in New Issue
Block a user