feat(thermal): execute certified backward Euler histories
This commit is contained in:
KOKO\Mimi
@@ -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
|
||||
Reference in New Issue
Block a user