KOKO\Mimi
3684 lines
177 KiB
C++
3684 lines
177 KiB
C++
#include <gtest/gtest.h>
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cfenv>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <cstdint>
|
|
#include <filesystem>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include <initializer_list>
|
|
#include <iostream>
|
|
#include <limits>
|
|
#include <memory>
|
|
#include <optional>
|
|
#include <sstream>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#include "analysis/hht_dynamic_internal.h"
|
|
#include "analysis/hht_reduced_internal.h"
|
|
#include "analysis/hht_reduced_increment.h"
|
|
#include "analysis/hht_reduced_procedure.h"
|
|
#include "fesa/analysis/analysis_model.h"
|
|
#include "fesa/analysis/newton_iteration_log.h"
|
|
#include "fesa/assembly/parallel_for.h"
|
|
#include "fesa/elements/finite_rotation.h"
|
|
#include "elements/mitc4_drilling_internal.h"
|
|
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
|
|
#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
|
|
|
|
namespace {
|
|
|
|
// User-approved physical/numerical comparison tolerance for the local experiment.
|
|
constexpr double kNumericalTolerance = 1.0e-3;
|
|
|
|
using Frame = fesa::DynamicAcceptedFrameState;
|
|
using Assembly = fesa::hht_internal::DynamicAssembly;
|
|
using L3 = std::array<long double, 3>;
|
|
constexpr double kEndTime = 0.25;
|
|
|
|
std::string Message(const fesa::Status& status) {
|
|
return status.Diagnostics().empty()
|
|
? "Unspecified solver failure"
|
|
: status.Diagnostics().front().code + ": " +
|
|
status.Diagnostics().front().message;
|
|
}
|
|
|
|
void Require(const fesa::Status& status) {
|
|
if (!status.IsOk()) throw std::runtime_error{Message(status)};
|
|
}
|
|
|
|
template <typename T>
|
|
T Require(fesa::Result<T> result) {
|
|
if (!result.HasValue()) throw std::runtime_error{Message(result.GetStatus())};
|
|
return std::move(result.Value());
|
|
}
|
|
|
|
L3 Add(const L3& a, const L3& b) {
|
|
return {a[0] + b[0], a[1] + b[1], a[2] + b[2]};
|
|
}
|
|
L3 Scale(const L3& a, const long double factor) {
|
|
return {factor * a[0], factor * a[1], factor * a[2]};
|
|
}
|
|
L3 Cross(const L3& a, const L3& b) {
|
|
return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
|
|
a[0] * b[1] - a[1] * b[0]};
|
|
}
|
|
long double Dot(const L3& a, const L3& b) {
|
|
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
|
|
}
|
|
long double Norm(const L3& a) {
|
|
return std::hypot(a[0], std::hypot(a[1], a[2]));
|
|
}
|
|
L3 Triple(const fesa::Vector& vector, const std::size_t offset) {
|
|
return {vector[offset], vector[offset + 1U], vector[offset + 2U]};
|
|
}
|
|
L3 Coordinates(const fesa::Node& node) {
|
|
return {node.coordinates[0], node.coordinates[1], node.coordinates[2]};
|
|
}
|
|
void Store(const L3& value, fesa::Vector* vector, const std::size_t offset) {
|
|
for (std::size_t c = 0; c < 3; ++c)
|
|
(*vector)[offset + c] = static_cast<double>(value[c]);
|
|
}
|
|
|
|
// Independent subset-block exponential, with no Rodrigues/Jl scalar formulas.
|
|
// The final block gives the mixed derivative in the supplied chart directions.
|
|
L3 ExponentialDerivative(const L3& phi, const L3& initial,
|
|
const std::vector<L3>& directions = {}) {
|
|
std::array<L3, 16> term{};
|
|
std::array<L3, 16> sum{};
|
|
term[0] = sum[0] = initial;
|
|
const std::size_t blocks = std::size_t{1} << directions.size();
|
|
for (std::size_t order = 1U; order <= 64U; ++order) {
|
|
std::array<L3, 16> next{};
|
|
for (std::size_t subset = 0; subset < blocks; ++subset) {
|
|
next[subset] = Cross(phi, term[subset]);
|
|
for (std::size_t direction = 0; direction < directions.size();
|
|
++direction) {
|
|
const std::size_t bit = std::size_t{1} << direction;
|
|
if ((subset & bit) != 0U)
|
|
next[subset] = Add(next[subset],
|
|
Cross(directions[direction], term[subset ^ bit]));
|
|
}
|
|
next[subset] =
|
|
Scale(next[subset], 1.0L / static_cast<long double>(order));
|
|
sum[subset] = Add(sum[subset], next[subset]);
|
|
}
|
|
term = next;
|
|
}
|
|
return sum[blocks - 1U];
|
|
}
|
|
|
|
struct Fixture {
|
|
std::unique_ptr<fesa::Domain> domain;
|
|
fesa::DofManager dofs;
|
|
std::vector<fesa::hht_internal::RuntimeShell> shells;
|
|
double length;
|
|
};
|
|
|
|
// The finite rigid-orbit and free-spin C08 fixtures prescribe the historical
|
|
// objective potential and its hidden velocity equation explicitly.
|
|
void UseLegacyDrilling(const Fixture& fixture) {
|
|
for (const auto& shell : fixture.shells)
|
|
*shell.internal_kernel = fesa::mitc4_drilling_internal::TestAccess::WithMode(
|
|
*shell.internal_kernel, fesa::Mitc4DrillingMode::kLegacyCorotatedTwist);
|
|
}
|
|
|
|
Fixture MakeFixture(const std::size_t divisions,
|
|
const bool lifecycle_translation = false) {
|
|
fesa::ModelDefinition definition{};
|
|
definition.source_path = "hht-shell-physics.inp";
|
|
definition.source_content_identity = "internal-C08-fixture";
|
|
for (std::size_t y = 0; y <= divisions; ++y) {
|
|
for (std::size_t x = 0; x <= divisions; ++x) {
|
|
const auto index = definition.nodes.size();
|
|
const auto label = static_cast<std::int64_t>(index + 1U);
|
|
definition.nodes.push_back({{"Plate", label, std::to_string(label)},
|
|
{-1.0 + 2.0 * static_cast<double>(x) /
|
|
static_cast<double>(divisions),
|
|
-1.0 + 2.0 * static_cast<double>(y) /
|
|
static_cast<double>(divisions),
|
|
0.0},
|
|
{}});
|
|
definition.shell_node_initial_frames.push_back(
|
|
{static_cast<fesa::EntityIndex>(index),
|
|
{0.0, 0.0, 1.0},
|
|
{1.0, 0.0, 0.0},
|
|
{0.0, 1.0, 0.0}});
|
|
}
|
|
}
|
|
definition.materials = {{"Elastic", 1000.0, 0.3, {}}};
|
|
definition.densities = {{0U, 1.0, {}}};
|
|
definition.shell_sections = {{"Centered", 0.2, 0U, {}}};
|
|
for (std::size_t y = 0; y < divisions; ++y) {
|
|
for (std::size_t x = 0; x < divisions; ++x) {
|
|
const auto first =
|
|
static_cast<fesa::EntityIndex>(y * (divisions + 1U) + x);
|
|
const auto stride = static_cast<fesa::EntityIndex>(divisions + 1U);
|
|
const auto label =
|
|
static_cast<std::int64_t>(definition.shell_elements.size() + 1U);
|
|
definition.shell_elements.push_back(
|
|
{{"Plate", label, std::to_string(label)},
|
|
fesa::ShellSourceElementType::kS4,
|
|
{first, first + 1U, first + stride + 1U, first + stride},
|
|
0U,
|
|
0U,
|
|
{}});
|
|
}
|
|
}
|
|
fesa::DynamicStepDefinition step{};
|
|
step.source_step_name_raw = "C08";
|
|
step.source_step_name_canonical = "C08";
|
|
step.initial_increment = 1.0 / 64.0;
|
|
step.time_period = kEndTime;
|
|
step.increment_count = 16U;
|
|
if (lifecycle_translation) {
|
|
std::vector<fesa::EntityIndex> nodes;
|
|
for (std::size_t node = 0U; node < definition.nodes.size(); ++node)
|
|
nodes.push_back(static_cast<fesa::EntityIndex>(node));
|
|
definition.node_sets.push_back({"All", std::nullopt, nodes, {}});
|
|
step.loads = {{"All", 1, 0.2, {}}};
|
|
step.increment_count = 2U;
|
|
step.time_period = 2.0 * step.initial_increment;
|
|
}
|
|
definition.dynamic_step = step;
|
|
auto domain = std::make_unique<fesa::Domain>(
|
|
Require(fesa::Domain::Create(std::move(definition))));
|
|
auto model = Require(fesa::AnalysisModel::Create(*domain));
|
|
auto dofs = Require(fesa::DofManager::Create(model));
|
|
auto shells = Require(fesa::hht_internal::BuildRuntimeShells(*domain, dofs));
|
|
const double length = Require(fesa::hht_internal::ReferenceLength(*domain));
|
|
return {std::move(domain), std::move(dofs), std::move(shells), length};
|
|
}
|
|
|
|
Frame EmptyFrame(const Fixture& fixture) {
|
|
Frame frame{};
|
|
frame.is_initial = true;
|
|
for (auto* vector :
|
|
{&frame.generalized_coordinates, &frame.velocity, &frame.acceleration,
|
|
&frame.algorithmic_acceleration, &frame.external_force}) {
|
|
*vector = fesa::Vector(fixture.dofs.FullDofCount());
|
|
}
|
|
return frame;
|
|
}
|
|
|
|
Assembly Evaluate(const Fixture& fixture, Frame frame, const double h,
|
|
const fesa::HhtParameters& parameters, const bool third) {
|
|
const fesa::SerialParallelFor parallel;
|
|
fesa::Vector load_scale(frame.external_force.Size());
|
|
for (std::size_t i = 0; i < load_scale.Size(); ++i)
|
|
load_scale[i] = std::abs(frame.external_force[i]);
|
|
return Require(fesa::hht_internal::EvaluateAssembly(
|
|
*fixture.domain, fixture.dofs, fixture.shells, std::move(frame),
|
|
load_scale, h, parameters, parallel, third));
|
|
}
|
|
|
|
struct VolumeOracle {
|
|
fesa::Vector force;
|
|
long double kinetic = 0.0L;
|
|
long double power = 0.0L;
|
|
long double power_scale = 0.0L;
|
|
L3 momentum{};
|
|
L3 angular_momentum{};
|
|
L3 momentum_rate{};
|
|
L3 angular_momentum_rate{};
|
|
long double momentum_scale = 0.0L;
|
|
long double angular_momentum_scale = 0.0L;
|
|
long double momentum_rate_scale = 0.0L;
|
|
long double angular_momentum_rate_scale = 0.0L;
|
|
};
|
|
|
|
// Independent nodal 5^3 rectangle quadrature and exponential kinematics.
|
|
// No production mass, inertia force, momentum, or director derivative is
|
|
// reused.
|
|
VolumeOracle Integrate(const Fixture& fixture, const Frame& frame) {
|
|
constexpr std::array<long double, 5> points{
|
|
-0.906179845938663992797626878299L, -0.538469310105683091036314420700L,
|
|
0.0L, 0.538469310105683091036314420700L,
|
|
0.906179845938663992797626878299L};
|
|
constexpr std::array<long double, 5> weights{
|
|
0.236926885056189087514264040720L, 0.478628670499366468041291514836L,
|
|
0.568888888888888888888888888889L, 0.478628670499366468041291514836L,
|
|
0.236926885056189087514264040720L};
|
|
constexpr std::array<std::array<int, 2>, 4> signs{
|
|
{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}}};
|
|
VolumeOracle oracle{fesa::Vector(fixture.dofs.FullDofCount())};
|
|
std::vector<long double> force(fixture.dofs.FullDofCount(), 0.0L);
|
|
for (std::size_t element = 0;
|
|
element < fixture.domain->ShellElements().Size(); ++element) {
|
|
const auto& shell = fixture.domain->ShellElements()[element];
|
|
std::array<L3, 4> x{};
|
|
std::array<L3, 4> v{};
|
|
std::array<L3, 4> a{};
|
|
std::array<L3, 4> director{};
|
|
std::array<L3, 4> director_rate{};
|
|
std::array<L3, 4> director_acceleration{};
|
|
std::array<std::array<L3, 3>, 4> derivative{};
|
|
for (std::size_t node = 0; node < 4; ++node) {
|
|
const auto global = shell.node_indices[node];
|
|
const auto offset = 6U * global;
|
|
const L3 phi = Triple(frame.generalized_coordinates, offset + 3U);
|
|
const L3 phidot = Triple(frame.velocity, offset + 3U);
|
|
const L3 phiddot = Triple(frame.acceleration, offset + 3U);
|
|
x[node] = Add(Coordinates(fixture.domain->Nodes()[global]),
|
|
Triple(frame.generalized_coordinates, offset));
|
|
v[node] = Triple(frame.velocity, offset);
|
|
a[node] = Triple(frame.acceleration, offset);
|
|
director[node] = ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L});
|
|
director_rate[node] =
|
|
ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {phidot});
|
|
director_acceleration[node] =
|
|
Add(ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {phidot, phidot}),
|
|
ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {phiddot}));
|
|
for (std::size_t c = 0; c < 3; ++c) {
|
|
L3 axis{};
|
|
axis[c] = 1.0L;
|
|
derivative[node][c] =
|
|
ExponentialDerivative(phi, {0.0L, 0.0L, 1.0L}, {axis});
|
|
}
|
|
}
|
|
const L3 edge_x =
|
|
Add(Coordinates(fixture.domain->Nodes()[shell.node_indices[1]]),
|
|
Scale(Coordinates(fixture.domain->Nodes()[shell.node_indices[0]]),
|
|
-1.0L));
|
|
const L3 edge_y =
|
|
Add(Coordinates(fixture.domain->Nodes()[shell.node_indices[3]]),
|
|
Scale(Coordinates(fixture.domain->Nodes()[shell.node_indices[0]]),
|
|
-1.0L));
|
|
const long double jacobian = Norm(Cross(edge_x, edge_y)) * 0.025L;
|
|
for (std::size_t i = 0; i < 5; ++i) {
|
|
for (std::size_t j = 0; j < 5; ++j) {
|
|
std::array<long double, 4> shape{};
|
|
for (std::size_t node = 0; node < 4; ++node)
|
|
shape[node] = 0.25L * (1.0L + signs[node][0] * points[i]) *
|
|
(1.0L + signs[node][1] * points[j]);
|
|
for (std::size_t k = 0; k < 5; ++k) {
|
|
const long double z = 0.1L * points[k];
|
|
const long double weight =
|
|
jacobian * weights[i] * weights[j] * weights[k];
|
|
for (std::size_t node = 0; node < 4; ++node) {
|
|
const long double nodal_weight = weight * shape[node];
|
|
const L3 position = Add(x[node], Scale(director[node], z));
|
|
const L3 velocity = Add(v[node], Scale(director_rate[node], z));
|
|
const L3 acceleration =
|
|
Add(a[node], Scale(director_acceleration[node], z));
|
|
oracle.kinetic += 0.5L * nodal_weight * Dot(velocity, velocity);
|
|
oracle.power += nodal_weight * Dot(velocity, acceleration);
|
|
oracle.power_scale +=
|
|
nodal_weight * Norm(velocity) * Norm(acceleration);
|
|
oracle.momentum =
|
|
Add(oracle.momentum, Scale(velocity, nodal_weight));
|
|
oracle.angular_momentum =
|
|
Add(oracle.angular_momentum,
|
|
Scale(Cross(position, velocity), nodal_weight));
|
|
oracle.momentum_rate =
|
|
Add(oracle.momentum_rate, Scale(acceleration, nodal_weight));
|
|
oracle.angular_momentum_rate =
|
|
Add(oracle.angular_momentum_rate,
|
|
Scale(Cross(position, acceleration), nodal_weight));
|
|
oracle.momentum_scale += nodal_weight * Norm(velocity);
|
|
oracle.angular_momentum_scale +=
|
|
nodal_weight * Norm(Cross(position, velocity));
|
|
oracle.momentum_rate_scale += nodal_weight * Norm(acceleration);
|
|
oracle.angular_momentum_rate_scale +=
|
|
nodal_weight * Norm(Cross(position, acceleration));
|
|
const auto offset = 6U * shell.node_indices[node];
|
|
for (std::size_t c = 0; c < 3; ++c) {
|
|
force[offset + c] += nodal_weight * acceleration[c];
|
|
force[offset + 3U + c] +=
|
|
nodal_weight * z * Dot(derivative[node][c], acceleration);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
for (std::size_t i = 0; i < force.size(); ++i)
|
|
oracle.force[i] = static_cast<double>(force[i]);
|
|
return oracle;
|
|
}
|
|
|
|
std::array<L3, 5> OrbitChart(const double time) {
|
|
const long double t = time;
|
|
return {{{0.4L * std::sin(t), 0.6L * std::sin(0.7L * t),
|
|
0.2L * std::sin(1.3L * t)},
|
|
{0.4L * std::cos(t), 0.42L * std::cos(0.7L * t),
|
|
0.26L * std::cos(1.3L * t)},
|
|
{-0.4L * std::sin(t), -0.294L * std::sin(0.7L * t),
|
|
-0.338L * std::sin(1.3L * t)},
|
|
{-0.4L * std::cos(t), -0.2058L * std::cos(0.7L * t),
|
|
-0.4394L * std::cos(1.3L * t)},
|
|
{0.4L * std::sin(t), 0.14406L * std::sin(0.7L * t),
|
|
0.57122L * std::sin(1.3L * t)}}};
|
|
}
|
|
|
|
Frame RigidOrbit(const Fixture& fixture, const double time) {
|
|
Frame frame = EmptyFrame(fixture);
|
|
frame.time = time;
|
|
frame.is_initial = time == 0.0;
|
|
const auto chart = OrbitChart(time);
|
|
const auto& phi = chart[0];
|
|
const auto& velocity = chart[1];
|
|
const auto& acceleration = chart[2];
|
|
for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node) {
|
|
const auto offset = 6U * node;
|
|
const L3 initial = Coordinates(fixture.domain->Nodes()[node]);
|
|
Store(Add(ExponentialDerivative(phi, initial), Scale(initial, -1.0L)),
|
|
&frame.generalized_coordinates, offset);
|
|
Store(phi, &frame.generalized_coordinates, offset + 3U);
|
|
Store(ExponentialDerivative(phi, initial, {velocity}), &frame.velocity,
|
|
offset);
|
|
Store(velocity, &frame.velocity, offset + 3U);
|
|
Store(Add(ExponentialDerivative(phi, initial, {velocity, velocity}),
|
|
ExponentialDerivative(phi, initial, {acceleration})),
|
|
&frame.acceleration, offset);
|
|
Store(acceleration, &frame.acceleration, offset + 3U);
|
|
}
|
|
frame.algorithmic_acceleration = frame.acceleration;
|
|
return frame;
|
|
}
|
|
|
|
std::array<L3, 3> OrbitPhysicalRates(const std::array<L3, 5>& chart,
|
|
const L3& initial) {
|
|
const auto& q = chart[0];
|
|
const auto& v = chart[1];
|
|
const auto& a = chart[2];
|
|
const auto& j = chart[3];
|
|
const auto& s = chart[4];
|
|
const auto d2 = Add(ExponentialDerivative(q, initial, {v, v}),
|
|
ExponentialDerivative(q, initial, {a}));
|
|
const auto d3 =
|
|
Add(Add(ExponentialDerivative(q, initial, {v, v, v}),
|
|
Scale(ExponentialDerivative(q, initial, {v, a}), 3.0L)),
|
|
ExponentialDerivative(q, initial, {j}));
|
|
const auto d4 =
|
|
Add(Add(ExponentialDerivative(q, initial, {v, v, v, v}),
|
|
Scale(ExponentialDerivative(q, initial, {v, v, a}), 6.0L)),
|
|
Add(Scale(ExponentialDerivative(q, initial, {a, a}), 3.0L),
|
|
Add(Scale(ExponentialDerivative(q, initial, {v, j}), 4.0L),
|
|
ExponentialDerivative(q, initial, {s}))));
|
|
return {d2, d3, d4};
|
|
}
|
|
|
|
// For the single centered rectangle, the nodal kinetic measure gives
|
|
// m_i=0.2 and j_i=m_i*t^2/12. Differentiate the nodal force moments.
|
|
std::array<fesa::Vector, 2> RigidForceRates(const Fixture& fixture,
|
|
const double time) {
|
|
std::array<fesa::Vector, 2> rates{fesa::Vector(24U), fesa::Vector(24U)};
|
|
const auto chart = OrbitChart(time);
|
|
const auto director = OrbitPhysicalRates(chart, {0.0L, 0.0L, 1.0L});
|
|
for (std::size_t node = 0; node < 4U; ++node) {
|
|
const auto x =
|
|
OrbitPhysicalRates(chart, Coordinates(fixture.domain->Nodes()[node]));
|
|
Store(Scale(x[1], 0.2L), &rates[0], 6U * node);
|
|
Store(Scale(x[2], 0.2L), &rates[1], 6U * node);
|
|
for (std::size_t c = 0; c < 3U; ++c) {
|
|
L3 axis{};
|
|
axis[c] = 1.0L;
|
|
const auto g =
|
|
ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L}, {axis});
|
|
const auto dg =
|
|
ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L}, {axis, chart[1]});
|
|
const auto ddg = Add(ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L},
|
|
{axis, chart[1], chart[1]}),
|
|
ExponentialDerivative(chart[0], {0.0L, 0.0L, 1.0L},
|
|
{axis, chart[2]}));
|
|
rates[0][6U * node + 3U + c] = static_cast<double>(
|
|
(0.008L / 12.0L) * (Dot(dg, director[0]) + Dot(g, director[1])));
|
|
rates[1][6U * node + 3U + c] =
|
|
static_cast<double>((0.008L / 12.0L) * (Dot(ddg, director[0]) +
|
|
2.0L * Dot(dg, director[1]) +
|
|
Dot(g, director[2])));
|
|
}
|
|
}
|
|
return rates;
|
|
}
|
|
|
|
double SparseNorm(const fesa::SparseMatrix& matrix) {
|
|
double norm = 0.0;
|
|
for (const double value : matrix.Values()) norm = std::hypot(norm, value);
|
|
return norm;
|
|
}
|
|
|
|
// The manufactured F(t) is prescribed in generalized coordinates. Its null
|
|
// work vanishes on the exact orbit, but not at a Newton trial. Include every
|
|
// explicit-time term of psi=n(q)^T(f_int(q)-F(t)); the public dead-load helper
|
|
// deliberately does not implement this verification-only load convention.
|
|
void ManufacturedHidden(const Assembly& assembly,
|
|
const std::array<fesa::Vector, 2>& force_rates,
|
|
Frame* frame) {
|
|
const auto nodes = assembly.rotations.size();
|
|
for (auto* vector :
|
|
{&frame->hidden_constraint, &frame->hidden_constraint_scale,
|
|
&frame->hidden_velocity_defect, &frame->hidden_velocity_defect_scale,
|
|
&frame->hidden_acceleration_defect,
|
|
&frame->hidden_acceleration_defect_scale})
|
|
*vector = fesa::Vector(nodes);
|
|
for (std::size_t node = 0; node < nodes; ++node) {
|
|
const auto& rotation = assembly.rotations[node];
|
|
double psi = 0.0;
|
|
double psi_scale = 0.0;
|
|
double first = 0.0;
|
|
double first_scale = 0.0;
|
|
double second = assembly.hidden_second[node];
|
|
double second_scale = assembly.hidden_second_scale[node];
|
|
auto add = [](const double term, double* value, double* scale) {
|
|
*value += term;
|
|
*scale += std::abs(term);
|
|
};
|
|
for (std::size_t c = 0; c < 3U; ++c) {
|
|
const auto row = 6U * node + 3U + c;
|
|
const double n = rotation.drill_null[c];
|
|
const double g = frame->internal_force[row] - frame->external_force[row];
|
|
psi += n * g;
|
|
psi_scale += std::abs(n) * frame->residual_scale[row];
|
|
for (std::size_t d = 0; d < 3U; ++d) {
|
|
const auto column = 6U * node + 3U + d;
|
|
const double dn = rotation.null_gradient(c, d);
|
|
add(dn * g * frame->velocity[column], &first, &first_scale);
|
|
add(dn * g * frame->acceleration[column], &second, &second_scale);
|
|
add(-2.0 * dn * frame->velocity[column] * force_rates[0][row], &second,
|
|
&second_scale);
|
|
}
|
|
const auto& stiffness = assembly.internal_tangent;
|
|
for (std::size_t index = stiffness.RowOffsets()[row];
|
|
index < stiffness.RowOffsets()[row + 1U]; ++index) {
|
|
const auto column = stiffness.ColumnIndices()[index];
|
|
const double gradient = n * stiffness.Values()[index];
|
|
add(gradient * frame->velocity[column], &first, &first_scale);
|
|
add(gradient * frame->acceleration[column], &second, &second_scale);
|
|
}
|
|
add(-n * force_rates[0][row], &first, &first_scale);
|
|
add(-n * force_rates[1][row], &second, &second_scale);
|
|
}
|
|
frame->hidden_constraint[node] = psi;
|
|
frame->hidden_constraint_scale[node] = psi_scale;
|
|
frame->hidden_velocity_defect[node] = first;
|
|
frame->hidden_velocity_defect_scale[node] = first_scale;
|
|
frame->hidden_acceleration_defect[node] = second;
|
|
frame->hidden_acceleration_defect_scale[node] = second_scale;
|
|
}
|
|
}
|
|
|
|
void CompleteManufacturedConstraints(
|
|
const Fixture& fixture, const Assembly& assembly,
|
|
const std::vector<std::size_t>& initial_null_nodes,
|
|
fesa::MklPardisoSolver* positive, fesa::MklPardisoTangentSolver* tangent,
|
|
Frame* frame) {
|
|
// This read-only probe exercises current positive mass, null inventory,
|
|
// psi and psi_q regularity. psi_q includes n'*(f_int-F), including moments.
|
|
// Its autonomous hidden-rate outputs are discarded. Explicit time terms
|
|
// and initial rate consistency are evaluated separately below.
|
|
auto probe = assembly.frame;
|
|
probe.is_initial = false;
|
|
Require(fesa::hht_internal::CompleteConstraintEvidence(
|
|
*fixture.domain, fixture.dofs, assembly, fixture.length,
|
|
initial_null_nodes, *positive, *tangent, &probe));
|
|
ManufacturedHidden(assembly, RigidForceRates(fixture, frame->time), frame);
|
|
if (frame->is_initial) {
|
|
for (std::size_t node = 0; node < assembly.rotations.size(); ++node) {
|
|
for (const auto quantity :
|
|
{std::make_pair(&frame->hidden_velocity_defect,
|
|
&frame->hidden_velocity_defect_scale),
|
|
std::make_pair(&frame->hidden_acceleration_defect,
|
|
&frame->hidden_acceleration_defect_scale)}) {
|
|
const double value = std::abs((*quantity.first)[node]);
|
|
const double scale = (*quantity.second)[node];
|
|
if (scale == 0.0)
|
|
EXPECT_DOUBLE_EQ(value, 0.0);
|
|
else
|
|
EXPECT_LE(value / scale, kNumericalTolerance);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(HhtShellPhysics, NodalRigidForceRatesMatchIndependentTimeDifferences) {
|
|
const auto fixture = MakeFixture(1U);
|
|
constexpr double time = 0.37;
|
|
const auto rates = RigidForceRates(fixture, time);
|
|
for (const double delta : {0.02, 0.01, 0.005}) {
|
|
SCOPED_TRACE(delta);
|
|
std::array<fesa::Vector, 5> forces{fesa::Vector(24U), fesa::Vector(24U),
|
|
fesa::Vector(24U), fesa::Vector(24U),
|
|
fesa::Vector(24U)};
|
|
for (std::size_t sample = 0; sample < forces.size(); ++sample) {
|
|
const double sample_time =
|
|
time + (static_cast<double>(sample) - 2.0) * delta;
|
|
forces[sample] =
|
|
Integrate(fixture, RigidOrbit(fixture, sample_time)).force;
|
|
}
|
|
for (std::size_t derivative = 0; derivative < 2U; ++derivative) {
|
|
for (const bool rotary : {false, true}) {
|
|
SCOPED_TRACE(derivative);
|
|
SCOPED_TRACE(rotary);
|
|
double error = 0.0;
|
|
double scale = 0.0;
|
|
for (std::size_t row = 0; row < 24U; ++row) {
|
|
if ((row % 6U >= 3U) != rotary) continue;
|
|
const double difference =
|
|
derivative == 0U ? (forces[0][row] - 8.0 * forces[1][row] +
|
|
8.0 * forces[3][row] - forces[4][row]) /
|
|
(12.0 * delta)
|
|
: (-forces[0][row] + 16.0 * forces[1][row] -
|
|
30.0 * forces[2][row] +
|
|
16.0 * forces[3][row] - forces[4][row]) /
|
|
(12.0 * delta * delta);
|
|
error = std::hypot(error, rates[derivative][row] - difference);
|
|
scale = std::hypot(scale, difference);
|
|
}
|
|
ASSERT_GT(scale, 0.0);
|
|
EXPECT_LE(error / scale, kNumericalTolerance);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(HhtShellPhysics, ExternalDrillingHessianBalancesGradientAcceleration) {
|
|
const auto fixture = MakeFixture(1U);
|
|
const fesa::Vector3 phi{0.0, 0.0, 0.0};
|
|
const fesa::Vector3 velocity{1.0, 0.0, 1.0};
|
|
const fesa::Vector3 acceleration{0.0, -1.0 / 3.0, 0.0};
|
|
const fesa::Vector3 force{1.0, 0.0, 0.0};
|
|
const auto external = fixture.shells[0U].internal->EvaluateExternalDrillingConstraint(
|
|
0U, phi, velocity, force, {}, {});
|
|
ASSERT_TRUE(external.HasValue());
|
|
double gradient_acceleration = 0.0;
|
|
for (std::size_t component = 0U; component < 3U; ++component)
|
|
gradient_acceleration +=
|
|
external.Value().gradient[component] * acceleration[component];
|
|
EXPECT_NEAR(external.Value().directional_second, -1.0 / 6.0, kNumericalTolerance);
|
|
EXPECT_NEAR(gradient_acceleration, 1.0 / 6.0, kNumericalTolerance);
|
|
EXPECT_NEAR(external.Value().directional_second + gradient_acceleration,
|
|
0.0, kNumericalTolerance);
|
|
constexpr double step = 1.0e-4;
|
|
const auto plus = fixture.shells[0U].internal->EvaluateExternalDrillingConstraint(
|
|
0U, {step, 0.0, step}, {}, force, {}, {});
|
|
const auto middle = fixture.shells[0U].internal->EvaluateExternalDrillingConstraint(
|
|
0U, phi, {}, force, {}, {});
|
|
const auto minus = fixture.shells[0U].internal->EvaluateExternalDrillingConstraint(
|
|
0U, {-step, 0.0, -step}, {}, force, {}, {});
|
|
ASSERT_TRUE(plus.HasValue());
|
|
ASSERT_TRUE(middle.HasValue());
|
|
ASSERT_TRUE(minus.HasValue());
|
|
const double finite_difference =
|
|
(plus.Value().value - 2.0 * middle.Value().value + minus.Value().value) /
|
|
(step * step);
|
|
EXPECT_NEAR(finite_difference, external.Value().directional_second, kNumericalTolerance);
|
|
}
|
|
|
|
TEST(HhtShellPhysics, ManufacturedHiddenCountsExternalHessianOnce) {
|
|
const auto fixture = MakeFixture(1U);
|
|
auto frame = EmptyFrame(fixture);
|
|
frame.velocity[3U] = frame.velocity[5U] = 1.0;
|
|
frame.acceleration[4U] = -1.0 / 3.0;
|
|
frame.external_force[3U] = 1.0;
|
|
const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
|
|
const auto assembly = Evaluate(fixture, frame, 1.0 / 64.0, parameters, true);
|
|
auto balance = assembly.frame;
|
|
std::array<fesa::Vector, 2U> zero_rates{
|
|
fesa::Vector(fixture.dofs.FullDofCount()),
|
|
fesa::Vector(fixture.dofs.FullDofCount())};
|
|
ManufacturedHidden(assembly, zero_rates, &balance);
|
|
constexpr std::size_t node = 0U;
|
|
double expected = assembly.hidden_second[node];
|
|
const auto& rotation = assembly.rotations[node];
|
|
for (std::size_t component = 0U; component < 3U; ++component) {
|
|
const auto row = 6U * node + 3U + component;
|
|
const double force =
|
|
assembly.frame.internal_force[row] - assembly.frame.external_force[row];
|
|
for (std::size_t direction = 0U; direction < 3U; ++direction)
|
|
expected += rotation.null_gradient(component, direction) * force *
|
|
assembly.frame.acceleration[6U * node + 3U + direction];
|
|
for (std::size_t index = assembly.internal_tangent.RowOffsets()[row];
|
|
index < assembly.internal_tangent.RowOffsets()[row + 1U]; ++index)
|
|
expected += rotation.drill_null[component] *
|
|
assembly.internal_tangent.Values()[index] *
|
|
assembly.frame.acceleration[
|
|
assembly.internal_tangent.ColumnIndices()[index]];
|
|
}
|
|
const double scale = assembly.hidden_second_scale[node] + std::abs(expected);
|
|
ASSERT_GT(scale, 0.0);
|
|
EXPECT_LE(std::abs(balance.hidden_acceleration_defect[node] - expected) /
|
|
scale,
|
|
kNumericalTolerance);
|
|
}
|
|
|
|
TEST(HhtShellPhysics,
|
|
C08aRigidOrbitMatchesIndependentVolumeForcePowerAndMomentum) {
|
|
const auto fixture = MakeFixture(1U);
|
|
UseLegacyDrilling(fixture);
|
|
const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
|
|
const auto initial =
|
|
Evaluate(fixture, EmptyFrame(fixture), 1.0 / 16.0, parameters, false);
|
|
const auto scaled_stiffness = Require(fesa::hht_internal::ScaledFree(
|
|
initial.internal_tangent, fixture.dofs, fixture.length));
|
|
double maximum_rotation = 0.0;
|
|
for (std::size_t sample = 0; sample <= 16; ++sample) {
|
|
const double time = static_cast<double>(sample) / 16.0;
|
|
SCOPED_TRACE(time);
|
|
auto exact = RigidOrbit(fixture, time);
|
|
const auto oracle = Integrate(fixture, exact);
|
|
exact.external_force = oracle.force;
|
|
const auto assembly =
|
|
Evaluate(fixture, exact, 1.0 / 16.0, parameters, true);
|
|
auto difference = assembly.frame.inertia_force;
|
|
difference.Axpy(-1.0, oracle.force);
|
|
EXPECT_LE(difference.Norm() / oracle.force.Norm(), kNumericalTolerance);
|
|
EXPECT_LE(std::abs(assembly.frame.kinetic_energy -
|
|
static_cast<double>(oracle.kinetic)) /
|
|
static_cast<double>(oracle.kinetic),
|
|
kNumericalTolerance);
|
|
EXPECT_LE(std::abs(oracle.force.Dot(exact.velocity) -
|
|
static_cast<double>(oracle.power)) /
|
|
static_cast<double>(oracle.power_scale),
|
|
kNumericalTolerance);
|
|
EXPECT_LE(assembly.frame.physical_strain_energy / 800.0, kNumericalTolerance);
|
|
EXPECT_LE(assembly.frame.drilling_energy / 800.0, kNumericalTolerance);
|
|
fesa::Vector scaled_q = exact.generalized_coordinates;
|
|
fesa::Vector scaled_force = assembly.frame.internal_force;
|
|
for (std::size_t i = 0; i < scaled_q.Size(); ++i) {
|
|
if (i % 6U < 3U) {
|
|
scaled_q[i] /= fixture.length;
|
|
scaled_force[i] *= fixture.length;
|
|
}
|
|
}
|
|
const double force_scale = SparseNorm(scaled_stiffness) * scaled_q.Norm();
|
|
if (force_scale == 0.0)
|
|
EXPECT_DOUBLE_EQ(scaled_force.Norm(), 0.0);
|
|
else
|
|
EXPECT_LE(scaled_force.Norm() / force_scale, kNumericalTolerance);
|
|
double momentum_error = 0.0;
|
|
double angular_error = 0.0;
|
|
for (std::size_t c = 0; c < 3; ++c) {
|
|
momentum_error = std::hypot(momentum_error,
|
|
assembly.frame.linear_momentum[c] -
|
|
static_cast<double>(oracle.momentum[c]));
|
|
angular_error = std::hypot(
|
|
angular_error, assembly.frame.angular_momentum[c] -
|
|
static_cast<double>(oracle.angular_momentum[c]));
|
|
}
|
|
EXPECT_LE(momentum_error / static_cast<double>(oracle.momentum_scale),
|
|
kNumericalTolerance);
|
|
EXPECT_LE(
|
|
angular_error / static_cast<double>(oracle.angular_momentum_scale),
|
|
kNumericalTolerance);
|
|
auto balance = assembly.frame;
|
|
ManufacturedHidden(assembly, RigidForceRates(fixture, time), &balance);
|
|
for (std::size_t node = 0; node < 4U; ++node) {
|
|
for (const auto quantity :
|
|
{std::make_pair(&balance.hidden_constraint,
|
|
&balance.hidden_constraint_scale),
|
|
std::make_pair(&balance.hidden_velocity_defect,
|
|
&balance.hidden_velocity_defect_scale),
|
|
std::make_pair(&balance.hidden_acceleration_defect,
|
|
&balance.hidden_acceleration_defect_scale)}) {
|
|
const double value = std::abs((*quantity.first)[node]);
|
|
const double scale = (*quantity.second)[node];
|
|
if (scale == 0.0)
|
|
EXPECT_DOUBLE_EQ(value, 0.0);
|
|
else
|
|
EXPECT_LE(value / scale, kNumericalTolerance);
|
|
}
|
|
}
|
|
EXPECT_TRUE(fesa::hht_internal::CompleteEquilibrium(
|
|
*fixture.domain, fixture.dofs, assembly, &balance)
|
|
.IsOk());
|
|
const double mass_norm = SparseNorm(assembly.mass);
|
|
for (std::size_t node = 0; node < 4U; ++node) {
|
|
fesa::Vector null(24U);
|
|
for (std::size_t c = 0; c < 3U; ++c)
|
|
null[6U * node + 3U + c] = assembly.rotations[node].drill_null[c];
|
|
EXPECT_LE(assembly.mass.Multiply(null).Norm() / mass_norm, kNumericalTolerance);
|
|
}
|
|
maximum_rotation = std::max(
|
|
maximum_rotation,
|
|
static_cast<double>(Norm(Triple(exact.generalized_coordinates, 3U))));
|
|
}
|
|
std::cout << "C08a pointwise independent-volume max_rotation="
|
|
<< std::setprecision(17) << maximum_rotation << '\n';
|
|
EXPECT_GT(maximum_rotation, 0.5);
|
|
}
|
|
|
|
// At q0=0 the normalized null is e_phi3 and f_int=0, hence psi_q is the
|
|
// corresponding stiffness row. Determine all null velocities, rather than
|
|
// assigning the unobservable director-axis chart rate arbitrarily.
|
|
void InitializeNullVelocity(const Fixture& fixture, const Assembly& initial,
|
|
Frame* frame,
|
|
fesa::MklPardisoTangentSolver* solver) {
|
|
const auto nodes = fixture.domain->Nodes().size();
|
|
std::vector<fesa::CooContribution> entries;
|
|
fesa::Vector rhs(nodes);
|
|
for (std::size_t node = 0; node < nodes; ++node) {
|
|
const auto row = 6U * node + 5U;
|
|
for (auto offset = initial.internal_tangent.RowOffsets()[row];
|
|
offset < initial.internal_tangent.RowOffsets()[row + 1U]; ++offset) {
|
|
const auto column = initial.internal_tangent.ColumnIndices()[offset];
|
|
const double value = initial.internal_tangent.Values()[offset];
|
|
if (column % 6U == 5U)
|
|
entries.push_back({node, column / 6U, value, 0U, entries.size()});
|
|
else
|
|
rhs[node] -= value * frame->velocity[column];
|
|
}
|
|
}
|
|
const auto regularity =
|
|
Require(fesa::hht_internal::MakeSparse(nodes, nodes, std::move(entries)));
|
|
fesa::Vector null_velocity(nodes);
|
|
Require(solver->FactorizeFresh(regularity));
|
|
Require(solver->SolveCorrection(rhs, null_velocity));
|
|
Require(fesa::hht_internal::CheckLinearSolve(regularity, rhs, null_velocity));
|
|
for (std::size_t node = 0; node < nodes; ++node)
|
|
frame->velocity[6U * node + 5U] = null_velocity[node];
|
|
}
|
|
|
|
struct GridEvidence {
|
|
double h;
|
|
double end_time;
|
|
std::vector<Frame> history;
|
|
std::array<double, 6> solve_noise{};
|
|
std::array<double, 6> arithmetic_noise{};
|
|
std::array<double, 6> selector_history_arithmetic_noise{};
|
|
std::array<double, 2> full_filter_arithmetic_noise{};
|
|
double energy_error = 0.0;
|
|
double momentum_error = 0.0;
|
|
double angular_momentum_error = 0.0;
|
|
double maximum_rotation = 0.0;
|
|
double algorithmic_acceleration_gap = 0.0;
|
|
double hidden_velocity = 0.0;
|
|
double hidden_acceleration = 0.0;
|
|
struct RecoveredAcceleration {
|
|
fesa::Vector value;
|
|
double pre_hidden2 = 0.0;
|
|
double pre_hidden2_relative = 0.0;
|
|
double post_hidden2 = 0.0;
|
|
double post_hidden2_relative = 0.0;
|
|
double mass_delta_a_relative = 0.0;
|
|
};
|
|
// These are diagnostic-only: RunPlate never substitutes them into Frame.
|
|
std::vector<RecoveredAcceleration> recovered_acceleration;
|
|
};
|
|
|
|
// Fixed-correction-path rounding radii for the carried-A recurrence. They
|
|
// bound arithmetic in the stored operation order only; Newton sensitivity is
|
|
// covered separately by the last-solve defect below.
|
|
struct ArithmeticRadii {
|
|
fesa::Vector coordinates;
|
|
fesa::Vector velocity;
|
|
fesa::Vector acceleration;
|
|
fesa::Vector algorithmic_acceleration;
|
|
};
|
|
|
|
double Up(const double value) {
|
|
if (!std::isfinite(value) || !(value >= 0.0))
|
|
return std::numeric_limits<double>::infinity();
|
|
return std::nextafter(value, std::numeric_limits<double>::infinity());
|
|
}
|
|
|
|
double UpSum(const double left, const double right) {
|
|
return Up(left + right);
|
|
}
|
|
|
|
double UpProduct(const double left, const double right) {
|
|
// Every caller is an upper nonnegative radius/product operation. Preserve
|
|
// exact zero support under an infinite inherited radius, but turn NaN or
|
|
// nonfinite nonzero products into infinity so max-based diagnostics cannot
|
|
// silently accept an incomplete certificate.
|
|
if (std::isnan(left) || std::isnan(right) || !(left >= 0.0) || !(right >= 0.0))
|
|
return std::numeric_limits<double>::infinity();
|
|
if (left == 0.0 || right == 0.0) return 0.0;
|
|
if (!std::isfinite(left) || !std::isfinite(right))
|
|
return std::numeric_limits<double>::infinity();
|
|
const double product = left * right;
|
|
return std::isfinite(product) ? Up(product)
|
|
: std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
double RoundoffRadius(const double value) {
|
|
const double u = std::ldexp(1.0, -53);
|
|
return UpSum(UpProduct(Up(u / (1.0 - u)), std::abs(value)),
|
|
std::numeric_limits<double>::denorm_min());
|
|
}
|
|
|
|
double AddRadius(const double result, const double left, const double right) {
|
|
return UpSum(UpSum(left, right), RoundoffRadius(result));
|
|
}
|
|
|
|
double ProductRadius(const double left_value, const double left_radius,
|
|
const double right_value, const double right_radius) {
|
|
const double result = left_value * right_value;
|
|
return UpSum(UpSum(UpProduct(std::abs(right_value), left_radius),
|
|
UpProduct(std::abs(left_value), right_radius)),
|
|
UpSum(UpProduct(left_radius, right_radius), RoundoffRadius(result)));
|
|
}
|
|
|
|
// Replay the actual fixed-branch Neumaier operations of rate reconstruction.
|
|
struct RateSumRadius {
|
|
double sum = 0.0, correction = 0.0;
|
|
double sum_radius = 0.0, correction_radius = 0.0;
|
|
void Add(const double term, const double term_radius) {
|
|
if (!std::isfinite(term) || !std::isfinite(term_radius) ||
|
|
!(term_radius >= 0.0))
|
|
throw std::runtime_error{"C08 arithmetic radius input is incomplete."};
|
|
const double next = sum + term;
|
|
const double next_radius = AddRadius(next, sum_radius, term_radius);
|
|
const bool sum_first = std::abs(sum) >= std::abs(term);
|
|
const double first = (sum_first ? sum : term) - next;
|
|
const double first_radius = AddRadius(first,
|
|
sum_first ? sum_radius : term_radius, next_radius);
|
|
const double error = first + (sum_first ? term : sum);
|
|
const double error_radius = AddRadius(error, first_radius,
|
|
sum_first ? term_radius : sum_radius);
|
|
correction += error;
|
|
correction_radius = AddRadius(correction, correction_radius, error_radius);
|
|
sum = next;
|
|
sum_radius = next_radius;
|
|
}
|
|
double Value() const { return sum + correction; }
|
|
double Radius() const {
|
|
return AddRadius(Value(), sum_radius, correction_radius);
|
|
}
|
|
};
|
|
|
|
double QuotientRadius(const double numerator_value, const double numerator_radius,
|
|
const double denominator_value,
|
|
const double denominator_radius) {
|
|
const double lower = std::nextafter(
|
|
std::abs(denominator_value) - denominator_radius,
|
|
-std::numeric_limits<double>::infinity());
|
|
if (!(lower > 0.0))
|
|
throw std::runtime_error{"C08 carried-Ar arithmetic denominator radius lost sign."};
|
|
const double result = numerator_value / denominator_value;
|
|
const double input = UpSum(numerator_radius,
|
|
UpProduct(std::abs(result), denominator_radius));
|
|
return UpSum(Up(input / lower), UpProduct(RoundoffRadius(result),
|
|
UpSum(1.0, Up(denominator_radius / lower))));
|
|
}
|
|
|
|
namespace reduced = fesa::hht_reduced_internal;
|
|
|
|
fesa::Vector ProjectBasisRadii(
|
|
const std::vector<fesa::hht_reduced_internal::FixedBasisColumn>& columns,
|
|
const fesa::Vector& full, const fesa::Vector& full_radii,
|
|
fesa::Vector* replay_center = nullptr) {
|
|
if (replay_center != nullptr) *replay_center = fesa::Vector(columns.size());
|
|
fesa::Vector result(columns.size());
|
|
for (std::size_t column = 0U; column < columns.size(); ++column) {
|
|
const auto& basis = columns[column];
|
|
if (basis.support.size() != basis.support_radii.size())
|
|
throw std::runtime_error{"C08 basis projection radius inventory is invalid."};
|
|
RateSumRadius sum;
|
|
for (std::size_t i = 0U; i < basis.support.size(); ++i) {
|
|
const auto& entry = basis.support[i];
|
|
if (entry.dof >= full.Size() || entry.dof >= full_radii.Size() ||
|
|
!std::isfinite(basis.support_radii[i]))
|
|
throw std::runtime_error{"C08 basis projection radius support is invalid."};
|
|
const double term = entry.value * full[entry.dof];
|
|
const double term_radius = ProductRadius(entry.value, basis.support_radii[i],
|
|
full[entry.dof], full_radii[entry.dof]);
|
|
sum.Add(term, term_radius);
|
|
}
|
|
if (replay_center != nullptr) (*replay_center)[column] = sum.Value();
|
|
result[column] = sum.Radius();
|
|
}
|
|
return result;
|
|
}
|
|
|
|
void AccumulateLiftBasisRadii(
|
|
const std::vector<fesa::hht_reduced_internal::FixedBasisColumn>& columns,
|
|
const fesa::Vector& coordinates, const fesa::Vector& coordinate_radii,
|
|
std::vector<RateSumRadius>* sums) {
|
|
if (coordinates.Size() != columns.size() || coordinate_radii.Size() != columns.size())
|
|
throw std::runtime_error{"C08 basis lift radius inventory is invalid."};
|
|
for (std::size_t column = 0U; column < columns.size(); ++column) {
|
|
const auto& basis = columns[column];
|
|
if (basis.support.size() != basis.support_radii.size())
|
|
throw std::runtime_error{"C08 basis lift radius support inventory is invalid."};
|
|
for (std::size_t i = 0U; i < basis.support.size(); ++i) {
|
|
const auto& entry = basis.support[i];
|
|
if (entry.dof >= sums->size() || !std::isfinite(basis.support_radii[i]))
|
|
throw std::runtime_error{"C08 basis lift radius support is invalid."};
|
|
const double term = entry.value * coordinates[column];
|
|
const double term_radius = ProductRadius(entry.value, basis.support_radii[i],
|
|
coordinates[column], coordinate_radii[column]);
|
|
(*sums)[entry.dof].Add(term, term_radius);
|
|
}
|
|
}
|
|
}
|
|
|
|
struct BasisLiftReplay {
|
|
fesa::Vector center;
|
|
fesa::Vector radius;
|
|
};
|
|
|
|
BasisLiftReplay ReplayLiftBasisRadii(
|
|
const std::vector<fesa::hht_reduced_internal::FixedBasisColumn>& columns,
|
|
const fesa::Vector& coordinates, const fesa::Vector& coordinate_radii,
|
|
const std::size_t full_count) {
|
|
std::vector<RateSumRadius> sums(full_count);
|
|
AccumulateLiftBasisRadii(columns, coordinates, coordinate_radii, &sums);
|
|
BasisLiftReplay result{fesa::Vector(full_count), fesa::Vector(full_count)};
|
|
for (std::size_t full = 0U; full < full_count; ++full) {
|
|
result.center[full] = sums[full].Value();
|
|
result.radius[full] = sums[full].Radius();
|
|
}
|
|
return result;
|
|
}
|
|
|
|
BasisLiftReplay ReplayLiftBasisRadii(
|
|
const std::vector<fesa::hht_reduced_internal::FixedBasisColumn>& p_columns,
|
|
const fesa::Vector& p, const fesa::Vector& p_radii,
|
|
const std::vector<fesa::hht_reduced_internal::FixedBasisColumn>& z_columns,
|
|
const fesa::Vector& z, const fesa::Vector& z_radii,
|
|
const std::size_t full_count) {
|
|
std::vector<RateSumRadius> sums(full_count);
|
|
// Match reduced::Lift: every p support precedes every z support.
|
|
AccumulateLiftBasisRadii(p_columns, p, p_radii, &sums);
|
|
AccumulateLiftBasisRadii(z_columns, z, z_radii, &sums);
|
|
BasisLiftReplay result{fesa::Vector(full_count), fesa::Vector(full_count)};
|
|
for (std::size_t full = 0U; full < full_count; ++full) {
|
|
result.center[full] = sums[full].Value();
|
|
result.radius[full] = sums[full].Radius();
|
|
}
|
|
return result;
|
|
}
|
|
|
|
fesa::Vector LiftBasisRadii(
|
|
const std::vector<fesa::hht_reduced_internal::FixedBasisColumn>& columns,
|
|
const fesa::Vector& coordinates, const fesa::Vector& coordinate_radii,
|
|
const std::size_t full_count) {
|
|
return ReplayLiftBasisRadii(columns, coordinates, coordinate_radii, full_count).radius;
|
|
}
|
|
|
|
TEST(HhtShellPhysics, BasisCoefficientRadiusSurvivesSparseProjectionAndLift) {
|
|
// A zero stored coefficient still has construction uncertainty. Dropping
|
|
// its support radius would incorrectly erase the free angular full-row
|
|
// contribution when a nonzero private plane coordinate is lifted.
|
|
const fesa::hht_reduced_internal::FixedBasisColumn plane{
|
|
0U, 1U, true,
|
|
{{3U, 0.0}, {4U, 0.6}, {5U, 0.8}},
|
|
{0x1p-52, 0x1p-53, 0x1p-53}};
|
|
const std::vector<fesa::hht_reduced_internal::FixedBasisColumn> columns{plane};
|
|
fesa::Vector coordinate(1U);
|
|
coordinate[0U] = 4.0;
|
|
const fesa::Vector zero_coordinate_radius(1U);
|
|
const auto lifted = LiftBasisRadii(columns, coordinate,
|
|
zero_coordinate_radius, 6U);
|
|
EXPECT_GT(lifted[3U], 0.0);
|
|
EXPECT_GE(lifted[4U], 4.0 * 0x1p-53);
|
|
EXPECT_GE(lifted[5U], 4.0 * 0x1p-53);
|
|
|
|
fesa::Vector full(6U);
|
|
full[4U] = 2.0;
|
|
const fesa::Vector zero_full_radius(6U);
|
|
const auto projected = ProjectBasisRadii(columns, full, zero_full_radius);
|
|
EXPECT_GT(projected[0U], 2.0 * 0x1p-53);
|
|
}
|
|
|
|
TEST(HhtShellPhysics, BasisRadiusReplayMatchesCompensatedProjectionAndCombinedLift) {
|
|
const double large = std::ldexp(1.0, 53);
|
|
fesa::Vector full(3U);
|
|
full[0U] = large;
|
|
full[1U] = 1.0;
|
|
full[2U] = -large;
|
|
const fesa::Vector zero_full_radius(3U);
|
|
const reduced::FixedBasisColumn projection{
|
|
0U, 0U, false, {{0U, 1.0}, {1U, 1.0}, {2U, 1.0}}, {0.0, 0.0, 0.0}};
|
|
fesa::Vector projection_center(0U);
|
|
const auto projection_radius = ProjectBasisRadii({projection}, full,
|
|
zero_full_radius,
|
|
&projection_center);
|
|
reduced::FixedChart projection_chart{};
|
|
projection_chart.full_dof_count = 3U;
|
|
projection_chart.p_dofs = {0U};
|
|
projection_chart.p_columns = {projection};
|
|
const auto projected = Require(reduced::ProjectPhysical(projection_chart, full));
|
|
EXPECT_DOUBLE_EQ(projected[0U], 1.0);
|
|
EXPECT_DOUBLE_EQ(projection_center[0U], projected[0U]);
|
|
EXPECT_GT(projection_radius[0U], 0.0);
|
|
|
|
reduced::FixedChart lift_chart{};
|
|
lift_chart.full_dof_count = 1U;
|
|
lift_chart.p_dofs = {0U, 1U};
|
|
lift_chart.z_dofs = {2U};
|
|
lift_chart.p_columns = {
|
|
{0U, 0U, false, {{0U, 1.0}}, {0.0}},
|
|
{0U, 0U, false, {{0U, 1.0}}, {0.0}}};
|
|
lift_chart.z_columns = {{0U, 1U, false, {{0U, 1.0}}, {0.0}}};
|
|
fesa::Vector p(2U);
|
|
p[0U] = large;
|
|
p[1U] = 1.0;
|
|
const fesa::Vector z{1U, -large};
|
|
const auto replay = ReplayLiftBasisRadii(lift_chart.p_columns, p, fesa::Vector(2U),
|
|
lift_chart.z_columns, z, fesa::Vector(1U), 1U);
|
|
const auto lifted = Require(reduced::Lift(lift_chart, p, z));
|
|
EXPECT_DOUBLE_EQ(lifted[0U], 1.0);
|
|
EXPECT_DOUBLE_EQ(replay.center[0U], lifted[0U]);
|
|
double plain = 0.0;
|
|
for (const double value : {large, 1.0, -large}) plain += value;
|
|
EXPECT_DOUBLE_EQ(plain, 0.0);
|
|
const double target = 1.25;
|
|
const double gap = target - replay.center[0U];
|
|
const double target_radius = UpSum(replay.radius[0U],
|
|
UpSum(std::abs(gap), RoundoffRadius(gap)));
|
|
EXPECT_GE(target_radius, 0.25);
|
|
|
|
fesa::Vector small_p(2U);
|
|
small_p[0U] = 1.0;
|
|
const auto small_replay = ReplayLiftBasisRadii(
|
|
lift_chart.p_columns, small_p, fesa::Vector(2U), lift_chart.z_columns,
|
|
fesa::Vector(1U), fesa::Vector(1U), 1U);
|
|
const auto small_lifted = Require(reduced::Lift(
|
|
lift_chart, small_p, fesa::Vector(1U)));
|
|
EXPECT_DOUBLE_EQ(small_replay.center[0U], small_lifted[0U]);
|
|
EXPECT_DOUBLE_EQ(small_replay.center[0U], 1.0);
|
|
ASSERT_LT(small_replay.radius[0U], 0.25);
|
|
const double small_gap = 1.25 - small_replay.center[0U];
|
|
const double small_target_radius = UpSum(
|
|
small_replay.radius[0U],
|
|
UpSum(std::abs(small_gap), RoundoffRadius(small_gap)));
|
|
EXPECT_GE(small_target_radius, 0.25);
|
|
}
|
|
|
|
TEST(HhtShellPhysics, BasisLiftedAlgorithmicCorrectionUsesHhtFilterFactor) {
|
|
reduced::FixedChart chart{};
|
|
chart.full_dof_count = 2U;
|
|
chart.p_dofs = {0U};
|
|
chart.p_columns = {{0U, 1U, true, {{0U, 0.6}, {1U, 0.8}}, {0.0, 0.0}}};
|
|
fesa::Vector delta_a(2U);
|
|
delta_a[0U] = 1.2;
|
|
delta_a[1U] = 1.6;
|
|
const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
|
|
const double delta_ar = (1.0 + parameters.alpha) * 2.0;
|
|
fesa::Vector delta_A = delta_a;
|
|
delta_A.Scale(1.0 + parameters.alpha);
|
|
const auto projected = Require(reduced::ProjectPhysical(chart, delta_A));
|
|
EXPECT_NEAR(projected[0U], delta_ar, 1.0e-15);
|
|
const auto omitted = Require(reduced::ProjectPhysical(chart, delta_a));
|
|
const double omission_scale = std::abs(omitted[0U]) + std::abs(delta_ar);
|
|
ASSERT_GT(omission_scale, 0.0);
|
|
EXPECT_GT(std::abs(omitted[0U] - delta_ar) / omission_scale, 1.0e-12);
|
|
}
|
|
|
|
TEST(HhtShellPhysics, CarriedArRadiiIncludeAbsorbedUpdatesAndInputHistory) {
|
|
ASSERT_EQ(std::fegetround(), FE_TONEAREST);
|
|
volatile double subnormal = std::numeric_limits<double>::denorm_min();
|
|
volatile double doubled = subnormal + subnormal;
|
|
ASSERT_DOUBLE_EQ(doubled, 2.0 * std::numeric_limits<double>::denorm_min());
|
|
const double tiny = std::ldexp(1.0, -54);
|
|
const double rounded = 1.0 + tiny;
|
|
ASSERT_DOUBLE_EQ(rounded, 1.0);
|
|
const double first = AddRadius(rounded, 0.0, 0.0);
|
|
EXPECT_GE(first, tiny);
|
|
const double second = AddRadius(rounded + tiny, first, 0.0);
|
|
EXPECT_GT(second, first);
|
|
EXPECT_GE(second, 2.0 * tiny);
|
|
EXPECT_GE(ProductRadius(2.0, tiny, 4.0, tiny), 6.0 * tiny);
|
|
EXPECT_GE(QuotientRadius(1.0, tiny, 2.0, 0.0), 0.5 * tiny);
|
|
EXPECT_THROW(QuotientRadius(1.0, 0.0, 1.0, 1.0), std::runtime_error);
|
|
EXPECT_TRUE(std::isinf(AddRadius(std::numeric_limits<double>::max(),
|
|
std::numeric_limits<double>::max(), 0.0)));
|
|
}
|
|
|
|
TEST(HhtShellPhysics, PrimaryCoordinateRetainsCorrectionAbsorbedByCarriedAr) {
|
|
// The unchanged LargeDIRECT fixture has a0=-800000, h=.25 and this pp.
|
|
const double pp = -11218.75;
|
|
const double factor = 0.275625 * 0.25 * 0.25;
|
|
const double ar = -pp / factor;
|
|
const double initial_p = pp + factor * ar;
|
|
const double delta_p = 1.0e-12;
|
|
const double next_ar = ar + delta_p / factor;
|
|
EXPECT_DOUBLE_EQ(pp + factor * next_ar, initial_p);
|
|
const double next_p = initial_p + delta_p;
|
|
EXPECT_GT(next_p, initial_p);
|
|
EXPECT_DOUBLE_EQ(next_p - initial_p, delta_p);
|
|
const double scale = std::abs(pp) + std::abs(factor * next_ar);
|
|
EXPECT_LE(std::abs(next_p - (pp + factor * next_ar)) / scale, 1.0e-12);
|
|
const double inherited_radius = 2.0e-10;
|
|
const double first_radius = AddRadius(next_p, inherited_radius, 0.0);
|
|
EXPECT_GT(first_radius, inherited_radius);
|
|
EXPECT_GT(AddRadius(next_p + delta_p, first_radius, 0.0), first_radius);
|
|
}
|
|
|
|
// Diagnostic error-free transform of one finite binary64 addition: exact
|
|
// input sum minus the stored sum. It does not modify the integrator state.
|
|
double AdditionLowPart(const double left, const double right, const double sum) {
|
|
const double right_virtual = sum - left;
|
|
return (left - (sum - right_virtual)) + (right - right_virtual);
|
|
}
|
|
|
|
// TwoSum preserves the exact high+delta sum as a pair. Only adding its low
|
|
// part to the previous low part rounds the represented pair total. The owned
|
|
// history radius enters once; the unused final low part encloses high output.
|
|
void AdvanceCoordinate(const double delta, double* high, double* low,
|
|
double* pair_radius) {
|
|
const double sum = *high + delta;
|
|
const double error = AdditionLowPart(*high, delta, sum);
|
|
const double low_sum = *low + error;
|
|
*high = sum + low_sum;
|
|
*low = AdditionLowPart(sum, low_sum, *high);
|
|
*pair_radius = UpSum(*pair_radius, RoundoffRadius(low_sum));
|
|
}
|
|
|
|
TEST(HhtShellPhysics, CompensatedCoordinateRetainsAbsorbedNewtonCorrections) {
|
|
const double delta = std::ldexp(1.0, -54);
|
|
double high = 1.0, low = 0.0, radius = 2.0e-10;
|
|
const double inherited_radius = radius;
|
|
AdvanceCoordinate(delta, &high, &low, &radius);
|
|
EXPECT_DOUBLE_EQ(high, 1.0);
|
|
EXPECT_DOUBLE_EQ(low, delta);
|
|
for (unsigned i = 1U; i < 4U; ++i)
|
|
AdvanceCoordinate(delta, &high, &low, &radius);
|
|
EXPECT_DOUBLE_EQ(high, 1.0 + std::ldexp(1.0, -52));
|
|
EXPECT_DOUBLE_EQ(low, 0.0);
|
|
EXPECT_GT(radius, inherited_radius);
|
|
EXPECT_GE(UpSum(radius, std::abs(low)), inherited_radius);
|
|
high = 1.0;
|
|
const double first = 1.0 + delta;
|
|
low = AdditionLowPart(1.0, delta, first) + AdditionLowPart(first, delta, first + delta);
|
|
radius = 0.0;
|
|
AdvanceCoordinate(delta, &high, &low, &radius);
|
|
EXPECT_DOUBLE_EQ(high, 1.0 + std::ldexp(1.0, -52));
|
|
EXPECT_DOUBLE_EQ(low, -delta);
|
|
high = std::numeric_limits<double>::denorm_min();
|
|
low = 0.0;
|
|
radius = 0.0;
|
|
AdvanceCoordinate(high, &high, &low, &radius);
|
|
EXPECT_DOUBLE_EQ(high, 2.0 * std::numeric_limits<double>::denorm_min());
|
|
EXPECT_DOUBLE_EQ(low, 0.0);
|
|
high = std::numeric_limits<double>::max();
|
|
AdvanceCoordinate(high, &high, &low, &radius);
|
|
EXPECT_FALSE(std::isfinite(high));
|
|
EXPECT_FALSE(std::isfinite(radius));
|
|
}
|
|
|
|
double HhtCoordinateFactorRadius(const fesa::HhtParameters& p, const double h) {
|
|
const double first = p.beta * h;
|
|
return ProductRadius(first, RoundoffRadius(first), h, 0.0);
|
|
}
|
|
|
|
double RadiusNorm(const fesa::Vector& radii, const std::size_t offset) {
|
|
if (offset + 3U > radii.Size())
|
|
throw std::runtime_error{"C08 radius channel is outside the full vector."};
|
|
return std::hypot(radii[offset], radii[offset + 1U], radii[offset + 2U]);
|
|
}
|
|
|
|
struct GraphBasisProducts {
|
|
fesa::Matrix graph;
|
|
fesa::Matrix physical;
|
|
fesa::Matrix graph_radius;
|
|
fesa::Matrix physical_radius;
|
|
};
|
|
|
|
struct GraphInverseCertificate {
|
|
fesa::Matrix inverse_action;
|
|
fesa::Matrix inverse_bound;
|
|
fesa::Matrix residual_operator_bound;
|
|
bool complete = false;
|
|
};
|
|
|
|
GraphBasisProducts ContractGraphBasis(
|
|
const reduced::ConstraintOperators& operators,
|
|
const reduced::FixedChart& chart) {
|
|
const auto contract = [&](const std::vector<reduced::FixedBasisColumn>& columns) {
|
|
fesa::Matrix value(operators.rows.size(), columns.size());
|
|
fesa::Matrix radius(operators.rows.size(), columns.size());
|
|
for (std::size_t row = 0U; row < operators.rows.size(); ++row) {
|
|
for (std::size_t column = 0U; column < columns.size(); ++column) {
|
|
const auto& basis = columns[column];
|
|
if (basis.support.size() != basis.support_radii.size())
|
|
throw std::runtime_error{"C08 graph basis radius inventory is invalid."};
|
|
RateSumRadius sum;
|
|
for (const auto& gradient : operators.rows[row].gradient) {
|
|
for (std::size_t support = 0U; support < basis.support.size(); ++support) {
|
|
const auto& entry = basis.support[support];
|
|
if (gradient.dof != entry.dof) continue;
|
|
sum.Add(gradient.value * entry.value,
|
|
ProductRadius(gradient.value, 0.0, entry.value,
|
|
basis.support_radii[support]));
|
|
}
|
|
}
|
|
value(row, column) = sum.Value();
|
|
radius(row, column) = sum.Radius();
|
|
}
|
|
}
|
|
return std::pair{std::move(value), std::move(radius)};
|
|
};
|
|
auto [graph, graph_radius] = contract(chart.z_columns);
|
|
auto [physical, physical_radius] = contract(chart.p_columns);
|
|
return {std::move(graph), std::move(physical), std::move(graph_radius),
|
|
std::move(physical_radius)};
|
|
}
|
|
|
|
fesa::SparseMatrix SparseFromDense(const fesa::Matrix& matrix) {
|
|
std::vector<fesa::CooContribution> entries;
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row)
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column)
|
|
if (matrix(row, column) != 0.0)
|
|
entries.push_back({row, column, matrix(row, column), 0U, entries.size()});
|
|
return Require(fesa::hht_internal::MakeSparse(matrix.Rows(), matrix.Columns(),
|
|
std::move(entries)));
|
|
}
|
|
|
|
double UpperRowSum(const fesa::Matrix& matrix, const std::size_t row) {
|
|
double sum = 0.0;
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column)
|
|
sum = UpSum(sum, std::abs(matrix(row, column)));
|
|
return sum;
|
|
}
|
|
|
|
bool FiniteMatrix(const fesa::Matrix& matrix) {
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row)
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column)
|
|
if (!std::isfinite(matrix(row, column))) return false;
|
|
return true;
|
|
}
|
|
|
|
bool FiniteNonnegativeMatrix(const fesa::Matrix& matrix) {
|
|
for (std::size_t row = 0U; row < matrix.Rows(); ++row)
|
|
for (std::size_t column = 0U; column < matrix.Columns(); ++column)
|
|
if (!std::isfinite(matrix(row, column)) || !(matrix(row, column) >= 0.0))
|
|
return false;
|
|
return true;
|
|
}
|
|
|
|
bool FiniteVector(const fesa::Vector& vector) {
|
|
for (std::size_t i = 0U; i < vector.Size(); ++i)
|
|
if (!std::isfinite(vector[i])) return false;
|
|
return true;
|
|
}
|
|
|
|
bool FiniteNonnegativeVector(const fesa::Vector& vector) {
|
|
for (std::size_t i = 0U; i < vector.Size(); ++i)
|
|
if (!std::isfinite(vector[i]) || !(vector[i] >= 0.0)) return false;
|
|
return true;
|
|
}
|
|
|
|
bool HasInfiniteRadius(const fesa::Vector& vector) {
|
|
bool infinite = false;
|
|
for (std::size_t i = 0U; i < vector.Size(); ++i) {
|
|
if (std::isnan(vector[i]) || !(vector[i] >= 0.0))
|
|
throw std::runtime_error{"C08 graph arithmetic radius is invalid."};
|
|
infinite = infinite || std::isinf(vector[i]);
|
|
}
|
|
return infinite;
|
|
}
|
|
|
|
void FillInfinite(fesa::Vector* vector) {
|
|
for (std::size_t i = 0U; i < vector->Size(); ++i)
|
|
(*vector)[i] = std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
GraphInverseCertificate CertifyGraphInverse(const GraphBasisProducts& products,
|
|
fesa::TangentSolver* solver) {
|
|
const auto size = products.graph.Rows();
|
|
if (products.graph.Columns() != size || products.graph_radius.Rows() != size ||
|
|
products.graph_radius.Columns() != size)
|
|
throw std::runtime_error{"C08 graph certificate matrix inventory is invalid."};
|
|
if (!FiniteMatrix(products.graph) ||
|
|
!FiniteNonnegativeMatrix(products.graph_radius))
|
|
return {fesa::Matrix(size, size), fesa::Matrix(size, size),
|
|
fesa::Matrix(size, size), false};
|
|
if (size == 0U)
|
|
return {fesa::Matrix(0U, 0U), fesa::Matrix(0U, 0U), fesa::Matrix(0U, 0U), true};
|
|
const auto sparse = SparseFromDense(products.graph);
|
|
Require(solver->FactorizeFresh(sparse));
|
|
fesa::Matrix inverse(size, size);
|
|
for (std::size_t column = 0U; column < size; ++column) {
|
|
fesa::Vector unit(size), action(size);
|
|
unit[column] = 1.0;
|
|
Require(solver->SolveCorrection(unit, action));
|
|
Require(fesa::hht_internal::CheckLinearSolve(sparse, unit, action));
|
|
if (!FiniteVector(action))
|
|
return {fesa::Matrix(size, size), fesa::Matrix(size, size),
|
|
fesa::Matrix(size, size), false};
|
|
for (std::size_t row = 0U; row < size; ++row) inverse(row, column) = action[row];
|
|
}
|
|
fesa::Matrix hbar(size, size);
|
|
for (std::size_t row = 0U; row < size; ++row) {
|
|
for (std::size_t column = 0U; column < size; ++column) {
|
|
RateSumRadius sum;
|
|
for (std::size_t middle = 0U; middle < size; ++middle)
|
|
sum.Add(products.graph(row, middle) * inverse(middle, column),
|
|
ProductRadius(products.graph(row, middle), 0.0,
|
|
inverse(middle, column), 0.0));
|
|
const double defect = (row == column ? 1.0 : 0.0) - sum.Value();
|
|
hbar(row, column) = UpSum(std::abs(defect),
|
|
AddRadius(defect, 0.0, sum.Radius()));
|
|
}
|
|
}
|
|
double theta = 0.0;
|
|
for (std::size_t row = 0U; row < size; ++row)
|
|
theta = std::max(theta, UpperRowSum(hbar, row));
|
|
const double lower = std::nextafter(1.0 - theta,
|
|
-std::numeric_limits<double>::infinity());
|
|
if (!(lower > 0.0) || !std::isfinite(lower))
|
|
return {std::move(inverse), fesa::Matrix(size, size), std::move(hbar), false};
|
|
fesa::Matrix bound(size, size);
|
|
for (std::size_t row = 0U; row < size; ++row) {
|
|
double tail = 0.0;
|
|
for (std::size_t middle = 0U; middle < size; ++middle)
|
|
tail = UpSum(tail, UpProduct(std::abs(inverse(row, middle)),
|
|
UpperRowSum(hbar, middle)));
|
|
tail = Up(tail / lower);
|
|
for (std::size_t column = 0U; column < size; ++column)
|
|
bound(row, column) = UpSum(std::abs(inverse(row, column)), tail);
|
|
}
|
|
return {std::move(inverse), std::move(bound), std::move(hbar), true};
|
|
}
|
|
|
|
void AddPhysicalContractionRadius(const GraphBasisProducts& products,
|
|
const fesa::Vector& coordinates,
|
|
const fesa::Vector& coordinate_radii,
|
|
fesa::Vector* rhs_radius) {
|
|
if (coordinates.Size() != products.physical.Columns() ||
|
|
coordinate_radii.Size() != coordinates.Size() ||
|
|
rhs_radius == nullptr || rhs_radius->Size() != products.physical.Rows())
|
|
throw std::runtime_error{"C08 physical graph-radius inventory is invalid."};
|
|
for (std::size_t row = 0U; row < products.physical.Rows(); ++row)
|
|
for (std::size_t column = 0U; column < products.physical.Columns(); ++column)
|
|
(*rhs_radius)[row] = UpSum((*rhs_radius)[row], ProductRadius(
|
|
products.physical(row, column), products.physical_radius(row, column),
|
|
coordinates[column], coordinate_radii[column]));
|
|
}
|
|
|
|
fesa::Vector BoundGraphSolution(const GraphBasisProducts& products,
|
|
const GraphInverseCertificate& inverse,
|
|
const fesa::Vector& rhs,
|
|
const fesa::Vector& rhs_radius,
|
|
const fesa::Vector& actual_center) {
|
|
const auto size = products.graph.Rows();
|
|
if (!inverse.complete || rhs.Size() != size || rhs_radius.Size() != size ||
|
|
actual_center.Size() != size)
|
|
throw std::runtime_error{"C08 graph certificate is incomplete."};
|
|
if (!FiniteMatrix(products.graph) ||
|
|
!FiniteNonnegativeMatrix(products.graph_radius) ||
|
|
!FiniteMatrix(inverse.inverse_action) ||
|
|
!FiniteNonnegativeMatrix(inverse.inverse_bound) ||
|
|
!FiniteNonnegativeMatrix(inverse.residual_operator_bound) ||
|
|
!FiniteVector(rhs) || !FiniteNonnegativeVector(rhs_radius) ||
|
|
!FiniteVector(actual_center))
|
|
throw std::runtime_error{"C08 graph certificate has nonfinite arithmetic input."};
|
|
// `actual_center` is the accepted graph/rate solve authority. Do not
|
|
// replace it with an independently rounded solve: this certificate bounds
|
|
// the actual mapped endpoint, while its N2-L evidence remains separate.
|
|
fesa::Vector rbar(size), source(size), y(size), row_a(size);
|
|
for (std::size_t row = 0U; row < size; ++row) {
|
|
RateSumRadius sum;
|
|
for (std::size_t column = 0U; column < size; ++column)
|
|
sum.Add(products.graph(row, column) * actual_center[column],
|
|
ProductRadius(products.graph(row, column), 0.0, actual_center[column], 0.0));
|
|
const double defect = rhs[row] - sum.Value();
|
|
rbar[row] = UpSum(std::abs(defect), AddRadius(defect, 0.0, sum.Radius()));
|
|
source[row] = UpSum(rbar[row], rhs_radius[row]);
|
|
for (std::size_t column = 0U; column < size; ++column)
|
|
source[row] = UpSum(source[row],
|
|
UpProduct(products.graph_radius(row, column), std::abs(actual_center[column])));
|
|
}
|
|
for (std::size_t row = 0U; row < size; ++row)
|
|
for (std::size_t column = 0U; column < size; ++column)
|
|
y[row] = UpSum(y[row], UpProduct(inverse.inverse_bound(row, column),
|
|
source[column]));
|
|
double theta = 0.0;
|
|
for (std::size_t row = 0U; row < size; ++row) {
|
|
for (std::size_t middle = 0U; middle < size; ++middle)
|
|
for (std::size_t column = 0U; column < size; ++column)
|
|
row_a[row] = UpSum(row_a[row], UpProduct(inverse.inverse_bound(row, middle),
|
|
products.graph_radius(middle, column)));
|
|
theta = std::max(theta, row_a[row]);
|
|
}
|
|
const double lower = std::nextafter(1.0 - theta,
|
|
-std::numeric_limits<double>::infinity());
|
|
if (!(lower > 0.0) || !std::isfinite(lower))
|
|
throw std::runtime_error{"C08 graph uncertainty feedback is incomplete."};
|
|
double norm_y = 0.0;
|
|
for (std::size_t row = 0U; row < size; ++row) norm_y = std::max(norm_y, y[row]);
|
|
fesa::Vector result(size);
|
|
for (std::size_t row = 0U; row < size; ++row)
|
|
result[row] = UpSum(y[row], UpProduct(row_a[row], Up(norm_y / lower)));
|
|
return result;
|
|
}
|
|
|
|
// N7 arithmetic INDICATOR with binary64 endpoint scalar/derivative operators
|
|
// frozen. This is first-order implicit sensitivity, not an enclosure of all
|
|
// FEM/AD operator evaluation, nonlinear branch, or equilibrium history error.
|
|
// CA selector history is retained separately; no radius is reset or clamped.
|
|
ArithmeticRadii GraphArithmeticIndicator(
|
|
const Frame& frame, const Frame& old,
|
|
const reduced::FixedChart& chart,
|
|
const reduced::ConstraintOperators& operators,
|
|
const ArithmeticRadii& selector_radii,
|
|
const ArithmeticRadii& old_radii, const fesa::HhtParameters& p,
|
|
fesa::MklPardisoTangentSolver* solver) {
|
|
auto result = selector_radii;
|
|
const auto nz = chart.z_columns.size();
|
|
if (nz == 0U) return result;
|
|
const bool inherited_incomplete =
|
|
HasInfiniteRadius(selector_radii.coordinates) ||
|
|
HasInfiniteRadius(selector_radii.velocity) ||
|
|
HasInfiniteRadius(selector_radii.acceleration) ||
|
|
HasInfiniteRadius(selector_radii.algorithmic_acceleration) ||
|
|
HasInfiniteRadius(old_radii.coordinates) ||
|
|
HasInfiniteRadius(old_radii.velocity) ||
|
|
HasInfiniteRadius(old_radii.acceleration) ||
|
|
HasInfiniteRadius(old_radii.algorithmic_acceleration);
|
|
if (inherited_incomplete) {
|
|
// +Inf is an inherited arithmetic bound, not a valid zero uncertainty.
|
|
// Retain it as INCONCLUSIVE without attempting max-based graph bounds.
|
|
FillInfinite(&result.coordinates);
|
|
FillInfinite(&result.velocity);
|
|
FillInfinite(&result.acceleration);
|
|
FillInfinite(&result.algorithmic_acceleration);
|
|
return result;
|
|
}
|
|
const auto products = ContractGraphBasis(operators, chart);
|
|
const auto graph = products.graph;
|
|
const auto assembled_graph = reduced::GraphJacobian(operators, chart);
|
|
for (std::size_t row = 0U; row < graph.Rows(); ++row)
|
|
for (std::size_t column = 0U; column < graph.Columns(); ++column)
|
|
if (graph(row, column) != assembled_graph(row, column))
|
|
throw std::runtime_error{"C08 graph CA replay differs from factored G."};
|
|
const auto inverse = CertifyGraphInverse(products, solver);
|
|
if (!inverse.complete)
|
|
throw std::runtime_error{"C08 graph inverse-action certificate is incomplete."};
|
|
const fesa::Vector zero_z(chart.z_columns.size());
|
|
const auto lift_p = [&](const fesa::Vector& coordinates) {
|
|
return Require(reduced::Lift(chart, coordinates, zero_z));
|
|
};
|
|
// Stored C/H/time entries are frozen inputs. The C*Lp/C*Lz replay below
|
|
// propagates immutable basis BL and source-order contraction rounding.
|
|
const auto p_q_coordinates = Require(reduced::ProjectPhysical(
|
|
chart, frame.generalized_coordinates));
|
|
const auto p_v_coordinates = Require(reduced::ProjectPhysical(chart, frame.velocity));
|
|
const auto p_a_coordinates = Require(reduced::ProjectPhysical(chart, frame.acceleration));
|
|
const auto p_q = lift_p(p_q_coordinates);
|
|
const auto p_v = lift_p(p_v_coordinates);
|
|
const auto p_a = lift_p(p_a_coordinates);
|
|
const auto z_q = Require(reduced::ProjectAlgebraic(chart, frame.generalized_coordinates));
|
|
const auto z_v = Require(reduced::ProjectAlgebraic(chart, frame.velocity));
|
|
const auto z_a = Require(reduced::ProjectAlgebraic(chart, frame.acceleration));
|
|
const auto p_coordinate_radii = ProjectBasisRadii(chart.p_columns,
|
|
frame.generalized_coordinates, selector_radii.coordinates);
|
|
const auto p_velocity_radii = ProjectBasisRadii(chart.p_columns,
|
|
frame.velocity, selector_radii.velocity);
|
|
const auto p_acceleration_radii = ProjectBasisRadii(chart.p_columns,
|
|
frame.acceleration, selector_radii.acceleration);
|
|
const auto p_q_full_radius = LiftBasisRadii(chart.p_columns,
|
|
p_q_coordinates,
|
|
p_coordinate_radii, chart.full_dof_count);
|
|
const auto p_v_full_radius = LiftBasisRadii(chart.p_columns,
|
|
p_v_coordinates,
|
|
p_velocity_radii, chart.full_dof_count);
|
|
const auto p_a_full_radius = LiftBasisRadii(chart.p_columns,
|
|
p_a_coordinates,
|
|
p_acceleration_radii, chart.full_dof_count);
|
|
for (std::size_t full = 0U; full < chart.full_dof_count; ++full) {
|
|
result.coordinates[full] = AddRadius(p_q[full], selector_radii.coordinates[full],
|
|
p_q_full_radius[full]);
|
|
result.velocity[full] = AddRadius(p_v[full], selector_radii.velocity[full],
|
|
p_v_full_radius[full]);
|
|
result.acceleration[full] = AddRadius(p_a[full], selector_radii.acceleration[full],
|
|
p_a_full_radius[full]);
|
|
}
|
|
const auto combine_lift_radius = [&](
|
|
const fesa::Vector& target, const fesa::Vector& inherited,
|
|
const fesa::Vector& p_coordinates, const fesa::Vector& p_radii,
|
|
const fesa::Vector& z_coordinates, const fesa::Vector& z_radii) {
|
|
const auto replay = ReplayLiftBasisRadii(chart.p_columns, p_coordinates, p_radii,
|
|
chart.z_columns, z_coordinates, z_radii,
|
|
chart.full_dof_count);
|
|
const auto actual = Require(reduced::Lift(chart, p_coordinates, z_coordinates));
|
|
fesa::Vector result(chart.full_dof_count);
|
|
for (std::size_t full = 0U; full < chart.full_dof_count; ++full) {
|
|
if (actual[full] != replay.center[full])
|
|
throw std::runtime_error{"C08 lift radius replay center differs from production Lift."};
|
|
const double gap = target[full] - replay.center[full];
|
|
const double gap_radius = UpSum(std::abs(gap), RoundoffRadius(gap));
|
|
result[full] = AddRadius(target[full], inherited[full],
|
|
UpSum(replay.radius[full], gap_radius));
|
|
}
|
|
return result;
|
|
};
|
|
const auto rate_defect = [&](const std::size_t row, const bool acceleration) {
|
|
const auto& op = operators.rows[row];
|
|
RateSumRadius rhs_sum;
|
|
if (acceleration) {
|
|
RateSumRadius hessian;
|
|
for (const auto& entry : op.hessian) {
|
|
const double first = entry.value * frame.velocity[entry.row];
|
|
hessian.Add(first * frame.velocity[entry.column], ProductRadius(first,
|
|
ProductRadius(entry.value, 0.0, frame.velocity[entry.row], 0.0),
|
|
frame.velocity[entry.column], 0.0));
|
|
}
|
|
RateSumRadius time_velocity;
|
|
for (const auto& entry : op.time_gradient)
|
|
time_velocity.Add(entry.value * frame.velocity[entry.dof],
|
|
ProductRadius(entry.value, 0.0, frame.velocity[entry.dof], 0.0));
|
|
rhs_sum.Add(-hessian.Value(), hessian.Radius());
|
|
rhs_sum.Add(-2.0 * time_velocity.Value(), ProductRadius(
|
|
2.0, 0.0, time_velocity.Value(), time_velocity.Radius()));
|
|
rhs_sum.Add(-op.time_second, 0.0);
|
|
} else {
|
|
rhs_sum.Add(-op.time_first, 0.0);
|
|
}
|
|
const auto& p_channel = acceleration ? p_a : p_v;
|
|
for (const auto& entry : op.gradient)
|
|
rhs_sum.Add(-entry.value * p_channel[entry.dof],
|
|
ProductRadius(entry.value, 0.0, p_channel[entry.dof], 0.0));
|
|
const double rhs = rhs_sum.Value();
|
|
const double radius = rhs_sum.Radius();
|
|
const auto& z_channel = acceleration ? z_a : z_v;
|
|
double dot = 0.0;
|
|
double dot_radius = 0.0;
|
|
for (std::size_t column = 0U; column < nz; ++column) {
|
|
if (graph(row, column) == 0.0) continue;
|
|
const double term = graph(row, column) * z_channel[column];
|
|
dot += term;
|
|
dot_radius = AddRadius(dot, dot_radius, ProductRadius(
|
|
graph(row, column), 0.0, z_channel[column], 0.0));
|
|
}
|
|
const double residual = dot - rhs;
|
|
return UpSum(std::abs(residual), AddRadius(residual, dot_radius, radius));
|
|
};
|
|
fesa::Vector position_rhs(nz), velocity_rhs(nz), acceleration_rhs(nz);
|
|
for (std::size_t row = 0U; row < nz; ++row) {
|
|
position_rhs[row] = UpSum(std::abs(operators.rows[row].value),
|
|
RoundoffRadius(operators.rows[row].value));
|
|
for (const auto& entry : operators.rows[row].gradient)
|
|
position_rhs[row] = UpSum(position_rhs[row],
|
|
UpProduct(std::abs(entry.value), result.coordinates[entry.dof]));
|
|
}
|
|
AddPhysicalContractionRadius(products, Require(reduced::ProjectPhysical(
|
|
chart, frame.generalized_coordinates)), p_coordinate_radii, &position_rhs);
|
|
const auto bz = BoundGraphSolution(products, inverse, graph.Multiply(z_q),
|
|
position_rhs, z_q);
|
|
result.coordinates = combine_lift_radius(frame.generalized_coordinates,
|
|
selector_radii.coordinates, p_q_coordinates, p_coordinate_radii, z_q, bz);
|
|
for (std::size_t row = 0U; row < nz; ++row) {
|
|
const auto& op = operators.rows[row];
|
|
double first = op.time_first;
|
|
double first_radius = 0.0;
|
|
for (const auto& entry : op.gradient) {
|
|
first += entry.value * frame.velocity[entry.dof];
|
|
first_radius = AddRadius(first, first_radius,
|
|
ProductRadius(entry.value, 0.0, frame.velocity[entry.dof], 0.0));
|
|
velocity_rhs[row] = UpSum(velocity_rhs[row],
|
|
UpProduct(std::abs(entry.value), result.velocity[entry.dof]));
|
|
}
|
|
for (const auto& entry : op.hessian)
|
|
velocity_rhs[row] = UpSum(velocity_rhs[row], UpProduct(
|
|
UpProduct(std::abs(entry.value), std::abs(frame.velocity[entry.column])),
|
|
result.coordinates[entry.row]));
|
|
for (const auto& entry : op.time_gradient)
|
|
velocity_rhs[row] = UpSum(velocity_rhs[row],
|
|
UpProduct(std::abs(entry.value), result.coordinates[entry.dof]));
|
|
velocity_rhs[row] = UpSum(velocity_rhs[row],
|
|
UpSum(rate_defect(row, false), UpSum(std::abs(first), first_radius)));
|
|
}
|
|
AddPhysicalContractionRadius(products, Require(reduced::ProjectPhysical(
|
|
chart, frame.velocity)), p_velocity_radii, &velocity_rhs);
|
|
const auto bvz = BoundGraphSolution(products, inverse, graph.Multiply(z_v),
|
|
velocity_rhs, z_v);
|
|
result.velocity = combine_lift_radius(frame.velocity, selector_radii.velocity,
|
|
p_v_coordinates, p_velocity_radii, z_v, bvz);
|
|
for (std::size_t row = 0U; row < nz; ++row) {
|
|
const auto& op = operators.rows[row];
|
|
double second = 0.0;
|
|
double second_radius = 0.0;
|
|
for (const auto& entry : op.hessian) {
|
|
const double first = entry.value * frame.velocity[entry.row];
|
|
second += first * frame.velocity[entry.column];
|
|
second_radius = AddRadius(second, second_radius, ProductRadius(first,
|
|
ProductRadius(entry.value, 0.0, frame.velocity[entry.row], 0.0),
|
|
frame.velocity[entry.column], 0.0));
|
|
acceleration_rhs[row] = UpSum(acceleration_rhs[row], UpProduct(
|
|
UpProduct(std::abs(entry.value), std::abs(frame.acceleration[entry.column])),
|
|
result.coordinates[entry.row]));
|
|
acceleration_rhs[row] = UpSum(acceleration_rhs[row], UpProduct(2.0, UpProduct(
|
|
UpProduct(std::abs(entry.value), std::abs(frame.velocity[entry.column])),
|
|
result.velocity[entry.row])));
|
|
}
|
|
for (const auto& entry : op.third_directional_gradient)
|
|
acceleration_rhs[row] = UpSum(acceleration_rhs[row],
|
|
UpProduct(std::abs(entry.value), result.coordinates[entry.dof]));
|
|
for (const auto& entry : op.time_gradient) {
|
|
const double first = 2.0 * entry.value;
|
|
second += first * frame.velocity[entry.dof];
|
|
second_radius = AddRadius(second, second_radius, ProductRadius(first,
|
|
ProductRadius(2.0, 0.0, entry.value, 0.0), frame.velocity[entry.dof], 0.0));
|
|
acceleration_rhs[row] = UpSum(acceleration_rhs[row], UpProduct(2.0,
|
|
UpProduct(std::abs(entry.value), result.velocity[entry.dof])));
|
|
}
|
|
second += op.time_second;
|
|
second_radius = AddRadius(second, second_radius, 0.0);
|
|
for (const auto& entry : op.gradient) {
|
|
second += entry.value * frame.acceleration[entry.dof];
|
|
second_radius = AddRadius(second, second_radius,
|
|
ProductRadius(entry.value, 0.0, frame.acceleration[entry.dof], 0.0));
|
|
acceleration_rhs[row] = UpSum(acceleration_rhs[row],
|
|
UpProduct(std::abs(entry.value), result.acceleration[entry.dof]));
|
|
}
|
|
acceleration_rhs[row] = UpSum(acceleration_rhs[row],
|
|
UpSum(rate_defect(row, true), UpSum(std::abs(second), second_radius)));
|
|
}
|
|
AddPhysicalContractionRadius(products, Require(reduced::ProjectPhysical(
|
|
chart, frame.acceleration)), p_acceleration_radii, &acceleration_rhs);
|
|
const auto baz = BoundGraphSolution(products, inverse, graph.Multiply(z_a),
|
|
acceleration_rhs, z_a);
|
|
result.acceleration = combine_lift_radius(frame.acceleration,
|
|
selector_radii.acceleration, p_a_coordinates, p_acceleration_radii, z_a, baz);
|
|
for (std::size_t i = 0U; i < chart.full_dof_count; ++i) {
|
|
result.algorithmic_acceleration[i] = AddRadius(
|
|
frame.algorithmic_acceleration[i], ProductRadius(1.0 + p.alpha,
|
|
RoundoffRadius(1.0 + p.alpha), frame.acceleration[i], result.acceleration[i]),
|
|
ProductRadius(p.alpha, 0.0, old.acceleration[i], old_radii.acceleration[i]));
|
|
}
|
|
return result;
|
|
}
|
|
|
|
TEST(HhtShellPhysics, FrozenGraphIndicatorIncludesCoordinateAndRateCurvature) {
|
|
// Independent polynomial graph z=p^2: vz=2*p*vp,
|
|
// az=2*vp^2+2*p*ap. Its first-order radius coefficients are exact.
|
|
reduced::FixedChart chart{};
|
|
chart.full_dof_count = 2U;
|
|
chart.p_dofs = {0U};
|
|
chart.z_dofs = {1U};
|
|
chart.p_columns = {{0U, 0U, false, {{0U, 1.0}}, {0.0}}};
|
|
chart.z_columns = {{0U, 1U, false, {{1U, 1.0}}, {0.0}}};
|
|
reduced::ConstraintOperators operators{};
|
|
operators.rows.resize(1U);
|
|
operators.rows[0].gradient = {{0U, -2.0}, {1U, 1.0}};
|
|
operators.rows[0].hessian = {{0U, 0U, -2.0}};
|
|
Frame frame{};
|
|
frame.generalized_coordinates = fesa::Vector(2U);
|
|
frame.velocity = fesa::Vector(2U);
|
|
frame.acceleration = fesa::Vector(2U);
|
|
frame.algorithmic_acceleration = fesa::Vector(2U);
|
|
frame.generalized_coordinates[0U] = frame.generalized_coordinates[1U] = 1.0;
|
|
frame.velocity[0U] = 1.0;
|
|
frame.velocity[1U] = 2.0;
|
|
frame.acceleration[0U] = 2.0;
|
|
frame.acceleration[1U] = 6.0;
|
|
Frame old = frame;
|
|
old.acceleration = fesa::Vector(2U);
|
|
const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
|
|
for (std::size_t i = 0U; i < 2U; ++i)
|
|
frame.algorithmic_acceleration[i] = (1.0 + parameters.alpha) * frame.acceleration[i];
|
|
ArithmeticRadii selector{fesa::Vector(2U), fesa::Vector(2U),
|
|
fesa::Vector(2U), fesa::Vector(2U)};
|
|
const ArithmeticRadii zero = selector;
|
|
selector.coordinates[0U] = 1.0e-6;
|
|
selector.velocity[0U] = 1.0e-7;
|
|
selector.acceleration[0U] = 1.0e-8;
|
|
fesa::MklPardisoTangentSolver solver;
|
|
const auto indicator = GraphArithmeticIndicator(frame, old, chart, operators,
|
|
selector, zero, parameters, &solver);
|
|
EXPECT_GE(indicator.coordinates[1U], 2.0e-6);
|
|
EXPECT_GE(indicator.velocity[1U], 2.0e-6 + 2.0e-7);
|
|
EXPECT_GE(indicator.acceleration[1U], 4.0e-6 + 4.0e-7 + 2.0e-8);
|
|
EXPECT_TRUE(std::isfinite(indicator.acceleration[1U]));
|
|
EXPECT_GE(indicator.algorithmic_acceleration[1U],
|
|
(1.0 + parameters.alpha) * indicator.acceleration[1U]);
|
|
selector.coordinates[0U] = std::numeric_limits<double>::infinity();
|
|
const auto unresolved = GraphArithmeticIndicator(frame, old, chart, operators,
|
|
selector, zero, parameters, &solver);
|
|
EXPECT_TRUE(std::isinf(unresolved.coordinates[1U]));
|
|
EXPECT_TRUE(std::isinf(unresolved.velocity[1U]));
|
|
EXPECT_TRUE(std::isinf(unresolved.acceleration[1U]));
|
|
}
|
|
|
|
TEST(HhtShellPhysics, GraphBasisCertificateBoundsMatrixAndRhsPerturbations) {
|
|
const auto solve = [](const fesa::Matrix& matrix, const fesa::Vector& rhs) {
|
|
const auto sparse = SparseFromDense(matrix);
|
|
fesa::MklPardisoTangentSolver solver;
|
|
Require(solver.FactorizeFresh(sparse));
|
|
fesa::Vector value(rhs.Size());
|
|
Require(solver.SolveCorrection(rhs, value));
|
|
Require(fesa::hht_internal::CheckLinearSolve(sparse, rhs, value));
|
|
return value;
|
|
};
|
|
{
|
|
GraphBasisProducts products{fesa::Matrix(1U, 1U), fesa::Matrix(1U, 0U),
|
|
fesa::Matrix(1U, 1U), fesa::Matrix(1U, 0U)};
|
|
products.graph(0U, 0U) = 2.0;
|
|
products.graph_radius(0U, 0U) = 0.1;
|
|
fesa::MklPardisoTangentSolver solver;
|
|
const auto inverse = CertifyGraphInverse(products, &solver);
|
|
ASSERT_TRUE(inverse.complete);
|
|
const fesa::Vector rhs{1U, 1.0};
|
|
const fesa::Vector rhs_radius{1U, 0.01};
|
|
const auto center = solve(products.graph, rhs);
|
|
const auto bound = BoundGraphSolution(products, inverse, rhs, rhs_radius, center);
|
|
auto perturbed = products.graph;
|
|
perturbed(0U, 0U) = 1.95;
|
|
const auto actual = solve(perturbed, fesa::Vector{1U, 1.009});
|
|
EXPECT_LE(std::abs(actual[0U] - center[0U]), bound[0U]);
|
|
EXPECT_GT(bound[0U], 0.0);
|
|
|
|
// The actual mapped graph center is the certificate authority. Replacing
|
|
// it with a newly solved center would lose this distance when BG=Bb=0.
|
|
products.graph_radius(0U, 0U) = 0.0;
|
|
const fesa::Vector shifted_center{1U, 0.5 + std::ldexp(1.0, -30)};
|
|
const auto centered = BoundGraphSolution(products, inverse, rhs,
|
|
fesa::Vector(1U), shifted_center);
|
|
EXPECT_GE(centered[0U], std::ldexp(1.0, -30));
|
|
|
|
// A zero stored contribution preserves exact zero under an inherited
|
|
// +Inf radius, whereas NaN is invalid. The certificate must not let an
|
|
// invalid graph radius reach max-based feedback as a silent zero.
|
|
EXPECT_EQ(UpProduct(0.0, std::numeric_limits<double>::infinity()), 0.0);
|
|
EXPECT_TRUE(std::isinf(UpProduct(0.0, std::numeric_limits<double>::quiet_NaN())));
|
|
products.graph_radius(0U, 0U) = std::numeric_limits<double>::infinity();
|
|
EXPECT_FALSE(CertifyGraphInverse(products, &solver).complete);
|
|
EXPECT_THROW(BoundGraphSolution(products, inverse, rhs,
|
|
fesa::Vector(1U), center), std::runtime_error);
|
|
}
|
|
{
|
|
GraphBasisProducts products{fesa::Matrix(4U, 4U), fesa::Matrix(4U, 0U),
|
|
fesa::Matrix(4U, 4U), fesa::Matrix(4U, 0U)};
|
|
for (std::size_t i = 0U; i < 4U; ++i) {
|
|
products.graph(i, i) = 1.0;
|
|
products.graph_radius(i, i) = 0.95;
|
|
}
|
|
fesa::MklPardisoTangentSolver solver;
|
|
const auto inverse = CertifyGraphInverse(products, &solver);
|
|
ASSERT_TRUE(inverse.complete);
|
|
const fesa::Vector rhs(4U, 1.0);
|
|
const auto center = solve(products.graph, rhs);
|
|
const auto bound = BoundGraphSolution(products, inverse, rhs, fesa::Vector(4U), center);
|
|
auto perturbed = products.graph;
|
|
for (std::size_t i = 0U; i < 4U; ++i) perturbed(i, i) = 0.1;
|
|
const auto actual = solve(perturbed, rhs);
|
|
for (std::size_t i = 0U; i < 4U; ++i) {
|
|
EXPECT_LE(std::abs(actual[i] - center[i]), bound[i]);
|
|
// Without the A feedback term this would be only y=.95 and fail for 9.
|
|
EXPECT_GT(bound[i], 9.0);
|
|
}
|
|
products.graph_radius(0U, 0U) = 1.0;
|
|
const auto invalid_inverse = CertifyGraphInverse(products, &solver);
|
|
ASSERT_TRUE(invalid_inverse.complete);
|
|
EXPECT_THROW(BoundGraphSolution(products, invalid_inverse, rhs,
|
|
fesa::Vector(4U), center), std::runtime_error);
|
|
}
|
|
}
|
|
|
|
struct LastSolve {
|
|
reduced::AugmentedReducedSystem full_system;
|
|
reduced::AugmentedSolveSelection selected_system;
|
|
fesa::Vector selected_correction;
|
|
fesa::Vector full_correction;
|
|
};
|
|
|
|
// Read-only N7 probes retain the actual last Newton matrix, including all
|
|
// auxiliary rows. The endpoint RHS includes graph and both hidden defects.
|
|
void RecordNoise(const Fixture& fixture, const Assembly& accepted,
|
|
const reduced::FixedChart& chart,
|
|
const reduced::ConstraintOperators& operators,
|
|
const LastSolve& last, const fesa::HhtParameters& p,
|
|
const ArithmeticRadii& arithmetic_radii,
|
|
fesa::MklPardisoTangentSolver* solver,
|
|
GridEvidence* evidence) {
|
|
const auto solve_defect = [&](const fesa::SparseMatrix& tangent,
|
|
const fesa::Vector& correction,
|
|
const fesa::Vector& rhs) {
|
|
auto defect = tangent.Multiply(correction);
|
|
defect.Axpy(-1.0, rhs);
|
|
defect.Scale(-1.0);
|
|
fesa::Vector solution(defect.Size());
|
|
Require(solver->FactorizeFresh(tangent));
|
|
Require(solver->SolveCorrection(defect, solution));
|
|
Require(fesa::hht_internal::CheckLinearSolve(tangent, defect, solution));
|
|
return solution;
|
|
};
|
|
const auto selected_defect = solve_defect(
|
|
last.selected_system.tangent, last.selected_correction,
|
|
last.selected_system.rhs);
|
|
const auto selected_full = Require(reduced::EmbedAugmentedSolution(
|
|
last.selected_system, selected_defect, last.full_system.rhs.Size()));
|
|
fesa::Vector model_full(last.full_system.rhs.Size());
|
|
if (last.selected_system.angular_only)
|
|
model_full = solve_defect(last.full_system.tangent, last.full_correction,
|
|
last.full_system.rhs);
|
|
const double cv = p.gamma / (p.beta * evidence->h);
|
|
const double ca = 1.0 / ((1.0 + p.alpha) * p.beta * evidence->h * evidence->h);
|
|
auto endpoint = Require(reduced::AssembleAugmentedReducedSystem(
|
|
accepted, chart, operators, cv, ca, fixture.length, evidence->h));
|
|
const auto hessian = reduced::ContractConstraintHessian(
|
|
operators, accepted.frame.velocity, accepted.frame.velocity);
|
|
for (std::size_t row = 0U; row < chart.z_columns.size(); ++row) {
|
|
double first = operators.rows[row].time_first;
|
|
double second = hessian[row] + operators.rows[row].time_second;
|
|
for (const auto& entry : operators.rows[row].gradient) {
|
|
first += entry.value * accepted.frame.velocity[entry.dof];
|
|
second += entry.value * accepted.frame.acceleration[entry.dof];
|
|
}
|
|
for (const auto& entry : operators.rows[row].time_gradient)
|
|
second += 2.0 * entry.value * accepted.frame.velocity[entry.dof];
|
|
endpoint.rhs[endpoint.physical_count + row] =
|
|
-operators.rows[row].value * endpoint.row_scale[endpoint.physical_count + row];
|
|
endpoint.rhs[endpoint.physical_count + endpoint.algebraic_count + row] =
|
|
-first * endpoint.row_scale[
|
|
endpoint.physical_count + endpoint.algebraic_count + row];
|
|
endpoint.rhs[endpoint.physical_count + 2U * endpoint.algebraic_count + row] =
|
|
-second * endpoint.row_scale[
|
|
endpoint.physical_count + 2U * endpoint.algebraic_count + row];
|
|
}
|
|
fesa::Vector nonlinear_defect(endpoint.rhs.Size());
|
|
// The endpoint nonlinear probe uses its freshly assembled actual augmented
|
|
// matrix; the linear probe above belongs to the frozen last Newton matrix.
|
|
Require(solver->FactorizeFresh(endpoint.tangent));
|
|
Require(solver->SolveCorrection(endpoint.rhs, nonlinear_defect));
|
|
Require(fesa::hht_internal::CheckLinearSolve(endpoint.tangent, endpoint.rhs,
|
|
nonlinear_defect));
|
|
const auto last_full = Require(reduced::MapAugmentedCorrection(
|
|
chart, last.full_system, last.full_correction, cv, ca));
|
|
const auto selected_linear_full = Require(reduced::MapAugmentedCorrection(
|
|
chart, last.full_system, selected_full, cv, ca));
|
|
std::optional<reduced::AugmentedCorrection> model_linear_full;
|
|
if (last.selected_system.angular_only)
|
|
model_linear_full = Require(reduced::MapAugmentedCorrection(
|
|
chart, last.full_system, model_full, cv, ca));
|
|
const auto nonlinear_full = Require(reduced::MapAugmentedCorrection(
|
|
chart, endpoint, nonlinear_defect, cv, ca));
|
|
const std::array<const fesa::Vector*, 3> radii{
|
|
&arithmetic_radii.coordinates, &arithmetic_radii.velocity,
|
|
&arithmetic_radii.acceleration};
|
|
std::array<fesa::Vector, 3> selector_radii{
|
|
LiftBasisRadii(chart.p_columns,
|
|
Require(reduced::ProjectPhysical(chart, accepted.frame.generalized_coordinates)),
|
|
ProjectBasisRadii(chart.p_columns, accepted.frame.generalized_coordinates,
|
|
*radii[0U]), chart.full_dof_count),
|
|
LiftBasisRadii(chart.p_columns,
|
|
Require(reduced::ProjectPhysical(chart, accepted.frame.velocity)),
|
|
ProjectBasisRadii(chart.p_columns, accepted.frame.velocity,
|
|
*radii[1U]), chart.full_dof_count),
|
|
LiftBasisRadii(chart.p_columns,
|
|
Require(reduced::ProjectPhysical(chart, accepted.frame.acceleration)),
|
|
ProjectBasisRadii(chart.p_columns, accepted.frame.acceleration,
|
|
*radii[2U]), chart.full_dof_count)};
|
|
const std::array<const reduced::AugmentedCorrection*, 4> corrections{
|
|
&last_full, &selected_linear_full,
|
|
model_linear_full.has_value() ? &*model_linear_full : nullptr,
|
|
&nonlinear_full};
|
|
for (std::size_t node = 0U; node < fixture.domain->Nodes().size(); ++node)
|
|
for (std::size_t channel = 0U; channel < 2U; ++channel) {
|
|
const auto offset = 6U * node + 3U * channel;
|
|
const double scale = channel == 0U ? fixture.length : 1.0;
|
|
evidence->full_filter_arithmetic_noise[channel] = std::max(
|
|
evidence->full_filter_arithmetic_noise[channel], Up(UpProduct(
|
|
evidence->end_time * evidence->end_time,
|
|
RadiusNorm(arithmetic_radii.algorithmic_acceleration, offset)) / scale));
|
|
for (std::size_t rate = 0U; rate < 3U; ++rate) {
|
|
const double time_scale = rate == 0U ? 1.0 : rate == 1U
|
|
? evidence->end_time : evidence->end_time * evidence->end_time;
|
|
const auto channel_values = [&](const reduced::AugmentedCorrection& value) {
|
|
return rate == 0U ? &value.coordinate
|
|
: rate == 1U ? &value.velocity : &value.acceleration;
|
|
};
|
|
double eta = 0.0;
|
|
for (const auto* correction : corrections)
|
|
if (correction != nullptr)
|
|
eta += static_cast<double>(Norm(Triple(*channel_values(*correction),
|
|
offset)));
|
|
evidence->solve_noise[2U * rate + channel] = std::max(
|
|
evidence->solve_noise[2U * rate + channel], time_scale * eta / scale);
|
|
evidence->selector_history_arithmetic_noise[2U * rate + channel] = std::max(
|
|
evidence->selector_history_arithmetic_noise[2U * rate + channel],
|
|
Up(UpProduct(time_scale, RadiusNorm(selector_radii[rate], offset)) / scale));
|
|
evidence->arithmetic_noise[2U * rate + channel] = std::max(
|
|
evidence->arithmetic_noise[2U * rate + channel], Up(UpProduct(time_scale,
|
|
RadiusNorm(*radii[rate], offset)) / scale));
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(HhtShellPhysics, RestrictedLastSolvePreservesOmittedTranslationDefect) {
|
|
// p0 is translational and p1 is rotational. The restricted solve exactly
|
|
// satisfies its angular row, while the omitted U row remains coupled and
|
|
// must therefore contribute the separate full-model N7 solve.
|
|
reduced::FixedChart chart{};
|
|
chart.full_dof_count = 6U;
|
|
chart.p_dofs = {0U, 3U};
|
|
chart.p_columns = {{0U, 0U, false, {{0U, 1.0}}, {0.0}},
|
|
{0U, 1U, false, {{3U, 1.0}}, {0.0}}};
|
|
const auto tangent = Require(fesa::hht_internal::MakeSparse(2U, 2U,
|
|
{{0U, 0U, 1.0, 0U, 0U}, {0U, 1U, 2.0, 0U, 1U},
|
|
{1U, 0U, 3.0, 0U, 2U}, {1U, 1U, 1.0, 0U, 3U}}));
|
|
reduced::AugmentedReducedSystem full{tangent, fesa::Vector(2U), 2U, 0U,
|
|
fesa::Vector(2U, 1.0),
|
|
fesa::Vector(2U, 1.0)};
|
|
full.rhs[0U] = 0.0;
|
|
full.rhs[1U] = 1.0;
|
|
ASSERT_TRUE(reduced::ShouldRefineAngular(chart, {0.0, 1.0},
|
|
{0.0, 0.0, 0.0, 0.0}, true));
|
|
const auto selected = Require(reduced::SelectAugmentedSolve(full, chart, true));
|
|
ASSERT_TRUE(selected.angular_only);
|
|
ASSERT_EQ(selected.full_indices, std::vector<std::size_t>({1U}));
|
|
fesa::MklPardisoTangentSolver solver;
|
|
fesa::Vector angular(selected.rhs.Size());
|
|
Require(solver.FactorizeFresh(selected.tangent));
|
|
Require(solver.SolveCorrection(selected.rhs, angular));
|
|
Require(fesa::hht_internal::CheckLinearSolve(selected.tangent, selected.rhs,
|
|
angular));
|
|
const auto embedded = Require(reduced::EmbedAugmentedSolution(
|
|
selected, angular, full.rhs.Size()));
|
|
EXPECT_EQ(embedded[0U], 0.0);
|
|
EXPECT_EQ(embedded[1U], 1.0);
|
|
const auto selected_residual = selected.tangent.Multiply(angular);
|
|
EXPECT_EQ(selected_residual[0U] - selected.rhs[0U], 0.0);
|
|
auto full_defect = full.rhs;
|
|
full_defect.Axpy(-1.0, full.tangent.Multiply(embedded));
|
|
EXPECT_EQ(full_defect[0U], -2.0);
|
|
EXPECT_EQ(full_defect[1U], 0.0);
|
|
fesa::Vector full_model(full_defect.Size());
|
|
Require(solver.FactorizeFresh(full.tangent));
|
|
Require(solver.SolveCorrection(full_defect, full_model));
|
|
Require(fesa::hht_internal::CheckLinearSolve(full.tangent, full_defect,
|
|
full_model));
|
|
EXPECT_GT(std::abs(full_model[0U]), 0.0);
|
|
EXPECT_GT(std::abs(full_model[1U]), 0.0);
|
|
}
|
|
|
|
TEST(HhtShellPhysics, RestrictedRecordNoiseIncludesFullModelOmittedUDefect) {
|
|
// Keep the actual shell accepted frame/chart/hidden endpoint rows. The
|
|
// frozen last operator is a minimal full-size coupled matrix: its J rows
|
|
// are exactly solved while one omitted translational U row retains a
|
|
// nonzero full-model defect.
|
|
const auto fixture = MakeFixture(1U);
|
|
constexpr double h = 1.0 / 64.0;
|
|
const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
|
|
const auto accepted = Evaluate(fixture, EmptyFrame(fixture), h, parameters, true);
|
|
const auto chart = Require(reduced::BuildFixedChart(fixture.dofs,
|
|
accepted.rotations));
|
|
const fesa::Vector zero(fixture.dofs.FullDofCount());
|
|
const auto operators = Require(reduced::EvaluateConstraintOperators(
|
|
*fixture.domain, fixture.dofs, fixture.shells, accepted, chart, zero, zero));
|
|
const double cv = parameters.gamma / (parameters.beta * h);
|
|
const double ca = 1.0 / ((1.0 + parameters.alpha) * parameters.beta * h * h);
|
|
const std::size_t total = chart.p_columns.size() + 3U * chart.z_columns.size();
|
|
std::size_t u = chart.p_columns.size();
|
|
std::size_t angular = chart.p_columns.size();
|
|
for (std::size_t i = 0U; i < chart.p_columns.size(); ++i) {
|
|
if (reduced::PhysicalCoordinateChannel(chart, i) == 0U) u = i;
|
|
else angular = i;
|
|
}
|
|
ASSERT_LT(u, chart.p_columns.size());
|
|
ASSERT_LT(angular, chart.p_columns.size());
|
|
std::vector<fesa::CooContribution> entries;
|
|
entries.reserve(total + 2U);
|
|
for (std::size_t i = 0U; i < total; ++i)
|
|
entries.push_back({i, i, 1.0, 0U, entries.size()});
|
|
entries.push_back({u, angular, 2.0, 0U, entries.size()});
|
|
entries.push_back({angular, u, 3.0, 0U, entries.size()});
|
|
const auto synthetic_tangent = Require(fesa::hht_internal::MakeSparse(
|
|
total, total, std::move(entries)));
|
|
reduced::AugmentedReducedSystem last_system{
|
|
synthetic_tangent, fesa::Vector(total), chart.p_columns.size(),
|
|
chart.z_columns.size(), fesa::Vector(total, 1.0), fesa::Vector(total, 1.0)};
|
|
fesa::Vector full_correction(total);
|
|
full_correction[angular] = 1.0;
|
|
last_system.rhs = last_system.tangent.Multiply(full_correction);
|
|
const auto baseline_selected = Require(reduced::SelectAugmentedSolve(
|
|
last_system, chart, true));
|
|
fesa::Vector selected_correction(baseline_selected.rhs.Size());
|
|
for (std::size_t i = 0U; i < selected_correction.Size(); ++i)
|
|
selected_correction[i] = full_correction[baseline_selected.full_indices[i]];
|
|
const LastSolve baseline{last_system, baseline_selected, selected_correction,
|
|
full_correction};
|
|
|
|
auto restricted_system = last_system;
|
|
restricted_system.rhs[u] = 0.0;
|
|
const auto restricted_selected = Require(reduced::SelectAugmentedSolve(
|
|
restricted_system, chart, true));
|
|
ASSERT_EQ(restricted_selected.full_indices, baseline_selected.full_indices);
|
|
const LastSolve restricted{restricted_system, restricted_selected,
|
|
selected_correction, full_correction};
|
|
auto full_defect = restricted_system.rhs;
|
|
full_defect.Axpy(-1.0, restricted_system.tangent.Multiply(full_correction));
|
|
ASSERT_NE(full_defect[u], 0.0);
|
|
fesa::MklPardisoTangentSolver solver;
|
|
fesa::Vector full_model(full_defect.Size());
|
|
Require(solver.FactorizeFresh(restricted_system.tangent));
|
|
Require(solver.SolveCorrection(full_defect, full_model));
|
|
Require(fesa::hht_internal::CheckLinearSolve(restricted_system.tangent,
|
|
full_defect, full_model));
|
|
const auto mapped_model = Require(reduced::MapAugmentedCorrection(
|
|
chart, restricted_system, full_model, cv, ca));
|
|
EXPECT_GT(mapped_model.coordinate.Norm(), 0.0);
|
|
|
|
const ArithmeticRadii radii{zero, zero, zero, zero};
|
|
GridEvidence baseline_evidence{h, h};
|
|
GridEvidence restricted_evidence{h, h};
|
|
RecordNoise(fixture, accepted, chart, operators, baseline, parameters,
|
|
radii, &solver, &baseline_evidence);
|
|
RecordNoise(fixture, accepted, chart, operators, restricted, parameters,
|
|
radii, &solver, &restricted_evidence);
|
|
double baseline_noise = 0.0;
|
|
double restricted_noise = 0.0;
|
|
for (std::size_t i = 0U; i < baseline_evidence.solve_noise.size(); ++i) {
|
|
baseline_noise += baseline_evidence.solve_noise[i];
|
|
restricted_noise += restricted_evidence.solve_noise[i];
|
|
}
|
|
EXPECT_GT(restricted_noise, baseline_noise);
|
|
}
|
|
|
|
GridEvidence::RecoveredAcceleration RecoverAcceleration(
|
|
const Fixture& fixture, const Assembly& assembly,
|
|
fesa::MklPardisoTangentSolver* solver, const bool manufactured) {
|
|
const auto nodes = assembly.rotations.size();
|
|
GridEvidence::RecoveredAcceleration result{assembly.frame.acceleration};
|
|
std::vector<double> weights(nodes);
|
|
std::vector<fesa::CooContribution> entries;
|
|
fesa::Vector rhs(nodes);
|
|
for (std::size_t row = 0; row < nodes; ++row) {
|
|
const auto& rotation = assembly.rotations[row];
|
|
fesa::Vector null(assembly.frame.acceleration.Size());
|
|
for (std::size_t c = 0; c < 3U; ++c)
|
|
null[6U * row + 3U + c] = rotation.drill_null[c];
|
|
const auto metric_null = assembly.drilling_metric.Multiply(null);
|
|
const double stiffness = null.Dot(metric_null);
|
|
if (!std::isfinite(stiffness) || !(stiffness > 0.0))
|
|
throw std::runtime_error{"C08 recovery requires positive null metric."};
|
|
weights[row] = std::sqrt(stiffness);
|
|
}
|
|
for (std::size_t row = 0; row < nodes; ++row) {
|
|
const auto& rotation = assembly.rotations[row];
|
|
const double raw = assembly.frame.hidden_acceleration_defect[row];
|
|
const double raw_scale =
|
|
assembly.frame.hidden_acceleration_defect_scale[row];
|
|
result.pre_hidden2 = std::max(result.pre_hidden2, std::abs(raw));
|
|
result.pre_hidden2_relative =
|
|
std::max(result.pre_hidden2_relative,
|
|
raw_scale == 0.0 ? std::abs(raw) : std::abs(raw) / raw_scale);
|
|
rhs[row] = -raw / weights[row];
|
|
for (std::size_t c = 0; c < 3U; ++c) {
|
|
const auto full = 6U * row + 3U + c;
|
|
const double force = assembly.frame.internal_force[full] -
|
|
assembly.frame.external_force[full];
|
|
for (std::size_t d = 0; d < 3U; ++d)
|
|
entries.push_back({row, row,
|
|
rotation.null_gradient(c, d) * force *
|
|
rotation.drill_null[d] / weights[row] /
|
|
weights[row],
|
|
0U, entries.size()});
|
|
for (std::size_t index = assembly.internal_tangent.RowOffsets()[full];
|
|
index < assembly.internal_tangent.RowOffsets()[full + 1U]; ++index) {
|
|
const auto column = assembly.internal_tangent.ColumnIndices()[index];
|
|
if (column % 6U < 3U) continue;
|
|
const auto column_node = column / 6U;
|
|
entries.push_back(
|
|
{row, column_node,
|
|
rotation.drill_null[c] *
|
|
assembly.internal_tangent.Values()[index] *
|
|
assembly.rotations[column_node].drill_null[column % 6U - 3U] /
|
|
weights[row] / weights[column_node],
|
|
0U, entries.size()});
|
|
}
|
|
}
|
|
}
|
|
const auto regularity =
|
|
Require(fesa::hht_internal::MakeSparse(nodes, nodes, std::move(entries)));
|
|
fesa::Vector scaled_eta(nodes);
|
|
Require(solver->FactorizeFresh(regularity));
|
|
Require(solver->SolveCorrection(rhs, scaled_eta));
|
|
Require(fesa::hht_internal::CheckLinearSolve(regularity, rhs, scaled_eta));
|
|
fesa::Vector delta(result.value.Size());
|
|
for (std::size_t node = 0; node < nodes; ++node)
|
|
for (std::size_t c = 0; c < 3U; ++c)
|
|
delta[6U * node + 3U + c] = assembly.rotations[node].drill_null[c] *
|
|
scaled_eta[node] / weights[node];
|
|
result.value.Axpy(1.0, delta);
|
|
const auto mass_delta = assembly.mass.Multiply(delta);
|
|
const double mass_scale = SparseNorm(assembly.mass) * delta.Norm();
|
|
result.mass_delta_a_relative =
|
|
mass_scale == 0.0 ? mass_delta.Norm() : mass_delta.Norm() / mass_scale;
|
|
if (manufactured) {
|
|
auto post = assembly.frame;
|
|
post.acceleration = result.value;
|
|
ManufacturedHidden(assembly, RigidForceRates(fixture, post.time), &post);
|
|
for (std::size_t row = 0; row < nodes; ++row) {
|
|
const double value = std::abs(post.hidden_acceleration_defect[row]);
|
|
const double scale = post.hidden_acceleration_defect_scale[row];
|
|
result.post_hidden2 = std::max(result.post_hidden2, value);
|
|
result.post_hidden2_relative = std::max(
|
|
result.post_hidden2_relative, scale == 0.0 ? value : value / scale);
|
|
}
|
|
return result;
|
|
}
|
|
for (std::size_t row = 0; row < nodes; ++row) {
|
|
fesa::Vector gradient(result.value.Size());
|
|
const auto& rotation = assembly.rotations[row];
|
|
for (std::size_t c = 0; c < 3U; ++c) {
|
|
const auto full = 6U * row + 3U + c;
|
|
const double force = assembly.frame.internal_force[full] -
|
|
assembly.frame.external_force[full];
|
|
for (std::size_t d = 0; d < 3U; ++d) {
|
|
const auto column = 6U * row + 3U + d;
|
|
gradient[column] += rotation.null_gradient(c, d) * force;
|
|
}
|
|
for (std::size_t index = assembly.internal_tangent.RowOffsets()[full];
|
|
index < assembly.internal_tangent.RowOffsets()[full + 1U]; ++index) {
|
|
const auto column = assembly.internal_tangent.ColumnIndices()[index];
|
|
gradient[column] +=
|
|
rotation.drill_null[c] * assembly.internal_tangent.Values()[index];
|
|
}
|
|
}
|
|
const double post =
|
|
assembly.hidden_second[row] + gradient.Dot(result.value);
|
|
double post_scale = assembly.hidden_second_scale[row];
|
|
for (std::size_t column = 0; column < gradient.Size(); ++column)
|
|
post_scale += std::abs(gradient[column] * result.value[column]);
|
|
result.post_hidden2 = std::max(result.post_hidden2, std::abs(post));
|
|
result.post_hidden2_relative = std::max(
|
|
result.post_hidden2_relative,
|
|
post_scale == 0.0 ? std::abs(post) : std::abs(post) / post_scale);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
void CompleteOriginalConstraints(const Fixture& fixture, const Assembly& assembly,
|
|
const std::vector<std::size_t>& null_nodes,
|
|
fesa::MklPardisoSolver& positive,
|
|
fesa::MklPardisoTangentSolver& tangent,
|
|
Frame* frame, const double h, const char* stage,
|
|
const reduced::EmbeddingRoot* root = nullptr) {
|
|
const auto status = fesa::hht_internal::CompleteConstraintEvidence(
|
|
*fixture.domain, fixture.dofs, assembly, fixture.length, null_nodes,
|
|
positive, tangent, frame);
|
|
if (status.IsOk()) return;
|
|
std::ostringstream context;
|
|
context << std::setprecision(17) << "C08 h=" << h
|
|
<< ", increment=" << frame->increment_index << ", time=" << frame->time
|
|
<< ", stage=" << stage << ": " << Message(status);
|
|
if (root != nullptr)
|
|
context << "; graph_iterations=" << root->graph_iteration_count
|
|
<< ", graph_residual=" << root->graph_residual_infinity
|
|
<< ", graph_relative=" << root->graph_residual_relative
|
|
<< ", graph_correction=" << root->graph_correction_infinity;
|
|
std::cout << std::setprecision(17) << "C08_HIDDEN_Q=";
|
|
for (std::size_t i = 0U; i < frame->generalized_coordinates.Size(); ++i)
|
|
std::cout << frame->generalized_coordinates[i] << ',';
|
|
std::cout << "\n";
|
|
const auto chart = reduced::BuildFixedChart(fixture.dofs, assembly.rotations);
|
|
if (chart.HasValue()) {
|
|
const fesa::Vector zero(fixture.dofs.FullDofCount());
|
|
const auto operators = reduced::EvaluateConstraintOperators(
|
|
*fixture.domain, fixture.dofs, fixture.shells, assembly, chart.Value(), zero, zero);
|
|
if (operators.HasValue()) {
|
|
for (std::size_t row = 0U; row < chart.Value().null_nodes.size(); ++row) {
|
|
const auto node = chart.Value().null_nodes[row];
|
|
double projected = 0.0, contraction_scale = 0.0, residual_scale = 0.0;
|
|
for (std::size_t c = 0U; c < 3U; ++c) {
|
|
const auto full = 6U * node + 3U + c;
|
|
const double n = assembly.rotations[node].drill_null[c];
|
|
const double term = n * (frame->internal_force[full] - frame->external_force[full]);
|
|
projected += term;
|
|
contraction_scale += std::abs(term);
|
|
residual_scale += std::abs(n) * frame->residual_scale[full];
|
|
std::cout << std::setprecision(17) << "C08_HIDDEN_COMPONENT h=" << h
|
|
<< " increment=" << frame->increment_index << " node=" << node
|
|
<< " full=" << full << " n=" << n
|
|
<< " internal=" << frame->internal_force[full]
|
|
<< " external=" << frame->external_force[full]
|
|
<< " inertia=" << frame->inertia_force[full]
|
|
<< " q=" << frame->generalized_coordinates[full] << "\n";
|
|
}
|
|
std::cout << std::setprecision(17) << "C08_HIDDEN_FAILURE h=" << h
|
|
<< " increment=" << frame->increment_index << " time=" << frame->time
|
|
<< " stage=" << stage << " node=" << node
|
|
<< " raw_psi=" << operators.Value().rows[row].value
|
|
<< " raw_scale=" << operators.Value().rows[row].value_scale
|
|
<< " projected_psi=" << projected << " contraction_scale=" << contraction_scale
|
|
<< " residual_scale=" << residual_scale
|
|
<< " difference=" << projected - operators.Value().rows[row].value << "\n";
|
|
}
|
|
}
|
|
}
|
|
throw std::runtime_error{context.str()};
|
|
}
|
|
|
|
GridEvidence RunPlate(const Fixture& fixture, const double h,
|
|
const bool manufactured) {
|
|
GridEvidence evidence{h, manufactured ? 1.0 : kEndTime};
|
|
const auto p = Require(fesa::MakeHhtParameters(-0.05));
|
|
fesa::MklPardisoSolver positive;
|
|
fesa::MklPardisoTangentSolver tangent;
|
|
Frame frame = manufactured ? RigidOrbit(fixture, 0.0) : EmptyFrame(fixture);
|
|
const L3 translation{0.125L, -0.25L, 0.0625L};
|
|
const L3 omega{0.8L, 0.6L, 0.4L};
|
|
if (manufactured) {
|
|
frame.external_force = Integrate(fixture, frame).force;
|
|
} else {
|
|
for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node) {
|
|
const auto x = Coordinates(fixture.domain->Nodes()[node]);
|
|
L3 velocity = Add(translation, Cross(omega, x));
|
|
velocity[2] += 0.02L * x[0] * x[1];
|
|
Store(velocity, &frame.velocity, 6U * node);
|
|
Store({0.8L, 0.6L, 0.0L}, &frame.velocity, 6U * node + 3U);
|
|
}
|
|
}
|
|
auto assembly = Evaluate(fixture, frame, h, p, true);
|
|
if (!manufactured) {
|
|
InitializeNullVelocity(fixture, assembly, &frame, &tangent);
|
|
for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node)
|
|
EXPECT_NEAR(frame.velocity[6U * node + 5U], 0.4, kNumericalTolerance);
|
|
assembly = Evaluate(fixture, frame, h, p, true);
|
|
}
|
|
const auto null_nodes =
|
|
fesa::hht_internal::NullNodes(fixture.dofs, assembly.rotations);
|
|
if (!manufactured)
|
|
frame.acceleration = Require(fesa::hht_internal::InitializeAcceleration(
|
|
*fixture.domain, fixture.dofs, assembly, fixture.length, positive,
|
|
tangent));
|
|
frame.algorithmic_acceleration = frame.acceleration;
|
|
frame.time_integration_formulation =
|
|
fesa::DynamicTimeIntegrationFormulation::kReducedCoordinateHht;
|
|
assembly = Evaluate(fixture, frame, h, p, true);
|
|
frame = assembly.frame;
|
|
// The fresh physical assembly certifies the initialized acceleration before
|
|
// it becomes frame-zero history; the initializer's private solve is separate.
|
|
for (const double residual : frame.free_residual_metrics) {
|
|
EXPECT_TRUE(std::isfinite(residual));
|
|
EXPECT_LE(residual, kNumericalTolerance);
|
|
}
|
|
if (manufactured)
|
|
CompleteManufacturedConstraints(fixture, assembly, null_nodes, &positive,
|
|
&tangent, &frame);
|
|
else
|
|
CompleteOriginalConstraints(fixture, assembly, null_nodes, positive, tangent,
|
|
&frame, h, "initial");
|
|
Require(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs,
|
|
assembly, &frame));
|
|
const auto initial_oracle = Integrate(fixture, frame);
|
|
const double energy0 = frame.kinetic_energy + frame.physical_strain_energy +
|
|
frame.drilling_energy;
|
|
EXPECT_GT(energy0, 0.0);
|
|
EXPECT_LE(std::abs(frame.kinetic_energy -
|
|
static_cast<double>(initial_oracle.kinetic)) /
|
|
energy0,
|
|
kNumericalTolerance);
|
|
const auto initial_momentum = frame.linear_momentum;
|
|
const auto initial_angular_momentum = frame.angular_momentum;
|
|
evidence.history.push_back(frame);
|
|
evidence.recovered_acceleration.push_back({frame.acceleration});
|
|
ArithmeticRadii arithmetic_radii{fesa::Vector(frame.generalized_coordinates.Size()),
|
|
fesa::Vector(frame.generalized_coordinates.Size()),
|
|
fesa::Vector(frame.generalized_coordinates.Size()),
|
|
fesa::Vector(frame.generalized_coordinates.Size())};
|
|
const auto chart = Require(reduced::BuildFixedChart(fixture.dofs, assembly.rotations));
|
|
fesa::Vector private_p = Require(reduced::ProjectPhysical(
|
|
chart, frame.generalized_coordinates));
|
|
fesa::Vector private_v = Require(reduced::ProjectPhysical(chart, frame.velocity));
|
|
fesa::Vector private_a = Require(reduced::ProjectPhysical(chart, frame.acceleration));
|
|
fesa::Vector private_ar = private_a;
|
|
fesa::Vector private_p_radius = ProjectBasisRadii(
|
|
chart.p_columns, frame.generalized_coordinates, arithmetic_radii.coordinates);
|
|
fesa::Vector private_v_radius = ProjectBasisRadii(
|
|
chart.p_columns, frame.velocity, arithmetic_radii.velocity);
|
|
fesa::Vector private_a_radius = ProjectBasisRadii(
|
|
chart.p_columns, frame.acceleration, arithmetic_radii.acceleration);
|
|
fesa::Vector private_ar_radius = private_a_radius;
|
|
fesa::Vector seed = Require(reduced::ProjectAlgebraic(
|
|
chart, frame.generalized_coordinates));
|
|
const fesa::Vector zero(fixture.dofs.FullDofCount());
|
|
const auto force_rates = [&](const double time) {
|
|
return manufactured ? RigidForceRates(fixture, time)
|
|
: std::array<fesa::Vector, 2>{zero, zero};
|
|
};
|
|
const auto position_operators = [&](const Frame& endpoint,
|
|
const fesa::Vector& q) {
|
|
auto position = endpoint;
|
|
position.generalized_coordinates = q;
|
|
position.velocity = zero;
|
|
position.acceleration = zero;
|
|
position.algorithmic_acceleration = zero;
|
|
const auto position_assembly = Evaluate(fixture, position, h, p, false);
|
|
const auto rates = force_rates(endpoint.time);
|
|
return reduced::EvaluateConstraintOperators(*fixture.domain, fixture.dofs,
|
|
fixture.shells, position_assembly, chart, rates[0], rates[1]);
|
|
};
|
|
const auto root0 = Require(reduced::SolveEmbeddingRoot(chart, private_p, seed,
|
|
[&](const fesa::Vector& q) { return position_operators(frame, q); }, tangent));
|
|
const auto rates0 = Require(reduced::ReconstructPhysicalRates(
|
|
chart, root0, private_v, private_a, tangent));
|
|
for (std::size_t i = 0U; i < zero.Size(); ++i) {
|
|
EXPECT_DOUBLE_EQ(root0.q[i], frame.generalized_coordinates[i]);
|
|
EXPECT_DOUBLE_EQ(rates0.velocity[i], frame.velocity[i]);
|
|
}
|
|
// Preserve the original frame zero. Certify independently reconstructed a0
|
|
// by the original physical equation and N3/N4 gates, without a component
|
|
// absolute floor or silently projecting the stored initial state.
|
|
auto reconstructed_zero = frame;
|
|
reconstructed_zero.generalized_coordinates = root0.q;
|
|
reconstructed_zero.velocity = rates0.velocity;
|
|
reconstructed_zero.acceleration = rates0.acceleration;
|
|
reconstructed_zero.algorithmic_acceleration = rates0.acceleration;
|
|
const auto reconstructed_assembly = Evaluate(fixture, reconstructed_zero, h, p, true);
|
|
reconstructed_zero = reconstructed_assembly.frame;
|
|
if (manufactured)
|
|
CompleteManufacturedConstraints(fixture, reconstructed_assembly, null_nodes,
|
|
&positive, &tangent, &reconstructed_zero);
|
|
else
|
|
CompleteOriginalConstraints(fixture, reconstructed_assembly, null_nodes,
|
|
positive, tangent, &reconstructed_zero, h,
|
|
"reconstructed-initial");
|
|
Require(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs,
|
|
reconstructed_assembly, &reconstructed_zero));
|
|
for (const double residual : reconstructed_zero.free_residual_metrics)
|
|
EXPECT_LE(residual, kNumericalTolerance);
|
|
Require(reduced::CheckGraphRegularity(assembly, chart, root0.operators, tangent));
|
|
const auto increments = static_cast<std::uint32_t>(evidence.end_time / h);
|
|
for (std::uint32_t increment = 1U; increment <= increments; ++increment) {
|
|
const Frame old = frame;
|
|
// FB-08 is private-coordinate history. Basis labels must never be
|
|
// reinterpreted as full indices or be fed through a full-vector predictor.
|
|
fesa::HhtPredictor predictor{fesa::Vector(private_p.Size()),
|
|
fesa::Vector(private_p.Size())};
|
|
fesa::Vector predictor_coordinate_radius(private_p.Size());
|
|
fesa::Vector predictor_velocity_radius(private_p.Size());
|
|
for (std::size_t i = 0U; i < private_p.Size(); ++i) {
|
|
const double velocity_term = h * private_v[i];
|
|
const double h2 = h * h;
|
|
const double acceleration_term = h2 * (0.5 - p.beta) * private_ar[i];
|
|
predictor.coordinates[i] = private_p[i] + velocity_term + acceleration_term;
|
|
predictor.velocity[i] = private_v[i] + h * (1.0 - p.gamma) * private_ar[i];
|
|
const double velocity_radius = ProductRadius(h, 0.0,
|
|
private_v[i], private_v_radius[i]);
|
|
const double h2_radius = ProductRadius(h, 0.0, h, 0.0);
|
|
const double acceleration_radius = ProductRadius(
|
|
h2 * (0.5 - p.beta),
|
|
ProductRadius(h2, h2_radius, 0.5 - p.beta,
|
|
RoundoffRadius(0.5 - p.beta)),
|
|
private_ar[i], private_ar_radius[i]);
|
|
predictor_coordinate_radius[i] = AddRadius(predictor.coordinates[i],
|
|
AddRadius(private_p[i] + velocity_term, private_p_radius[i], velocity_radius),
|
|
acceleration_radius);
|
|
const double rate_factor = h * (1.0 - p.gamma);
|
|
predictor_velocity_radius[i] = AddRadius(predictor.velocity[i],
|
|
private_v_radius[i], ProductRadius(rate_factor,
|
|
ProductRadius(h, 0.0, 1.0 - p.gamma,
|
|
RoundoffRadius(1.0 - p.gamma)),
|
|
private_ar[i], private_ar_radius[i]));
|
|
}
|
|
frame.is_initial = false;
|
|
frame.frame_index = increment;
|
|
frame.increment_index = increment;
|
|
frame.time = static_cast<double>(increment) * h;
|
|
frame.time_integration_formulation =
|
|
fesa::DynamicTimeIntegrationFormulation::kReducedCoordinateHht;
|
|
if (manufactured)
|
|
frame.external_force = Integrate(fixture, RigidOrbit(fixture, frame.time)).force;
|
|
bool correction_pass = true;
|
|
// These are the last *applied* translation-channel corrections. They
|
|
// deliberately remain outside the all-channel acceptance flag: the
|
|
// approved restricted solve is eligible only after this channel itself
|
|
// has passed its original four mapped-update checks.
|
|
std::array<double, 4U> translation_corrections{};
|
|
bool translation_private_algorithmic_pass = true;
|
|
const auto private_algorithmic_equal = [](const double actual,
|
|
const double reduced_value) {
|
|
const double defect = std::abs(actual - reduced_value);
|
|
const double scale = std::abs(reduced_value);
|
|
return std::isfinite(defect) && std::isfinite(scale) &&
|
|
(scale == 0.0 ? defect == 0.0 : defect / scale <= 1.0e-12);
|
|
};
|
|
bool converged = false;
|
|
std::optional<LastSolve> last;
|
|
fesa::Vector ar(chart.p_columns.size());
|
|
fesa::Vector ar_radius(chart.p_columns.size());
|
|
fesa::Vector candidate_p(chart.p_columns.size());
|
|
fesa::Vector candidate_p_radius(chart.p_columns.size());
|
|
fesa::Vector candidate_p_low(chart.p_columns.size());
|
|
fesa::Vector candidate_p_pair_radius(chart.p_columns.size());
|
|
fesa::Vector chosen_start_p(chart.p_columns.size());
|
|
fesa::Vector last_p_before_update(chart.p_columns.size());
|
|
fesa::Vector last_requested_delta(chart.p_columns.size());
|
|
fesa::Vector addition_low_sum(chart.p_columns.size());
|
|
fesa::Vector requested_sum(chart.p_columns.size());
|
|
fesa::Vector requested_compensation(chart.p_columns.size());
|
|
fesa::Vector requested_absolute_sum(chart.p_columns.size());
|
|
bool copy_translation_rates = false;
|
|
fesa::Vector retained_velocity(chart.p_columns.size());
|
|
fesa::Vector retained_acceleration(chart.p_columns.size());
|
|
fesa::Vector retained_velocity_radius(chart.p_columns.size());
|
|
fesa::Vector retained_acceleration_radius(chart.p_columns.size());
|
|
const auto record_requested_increment = [&](const std::size_t i, const double delta) {
|
|
const double previous = requested_sum[i];
|
|
const double next = previous + delta;
|
|
requested_compensation[i] += std::abs(previous) >= std::abs(delta)
|
|
? (previous - next) + delta : (delta - next) + previous;
|
|
requested_sum[i] = next;
|
|
if (delta != 0.0)
|
|
requested_absolute_sum[i] = UpSum(requested_absolute_sum[i], std::abs(delta));
|
|
};
|
|
for (std::size_t i = 0U; i < candidate_p.Size(); ++i) {
|
|
candidate_p[i] = predictor.coordinates[i] + p.beta * h * h * ar[i];
|
|
chosen_start_p[i] = candidate_p[i];
|
|
const double velocity_term = h * private_v[i];
|
|
const double acceleration_term = h * h * (0.5 - p.beta) * private_ar[i];
|
|
const double first_sum = private_p[i] + velocity_term;
|
|
addition_low_sum[i] = AdditionLowPart(private_p[i], velocity_term, first_sum) +
|
|
AdditionLowPart(first_sum, acceleration_term, predictor.coordinates[i]);
|
|
candidate_p_low[i] = addition_low_sum[i];
|
|
record_requested_increment(i, h * private_v[i]);
|
|
record_requested_increment(i, h * h * (0.5 - p.beta) * private_ar[i]);
|
|
candidate_p_radius[i] = AddRadius(candidate_p[i],
|
|
predictor_coordinate_radius[i], ProductRadius(p.beta * h * h,
|
|
HhtCoordinateFactorRadius(p, h), ar[i], ar_radius[i]));
|
|
candidate_p_pair_radius[i] = UpSum(candidate_p_radius[i],
|
|
std::abs(candidate_p_low[i]));
|
|
}
|
|
ArithmeticRadii trial_radii{zero, zero, zero, zero};
|
|
struct Trial {
|
|
Assembly assembly;
|
|
reduced::EmbeddingRoot root;
|
|
fesa::Vector vp;
|
|
fesa::Vector ap;
|
|
fesa::Vector ar;
|
|
};
|
|
const auto evaluate_trial = [&](const fesa::Vector& trial_ar,
|
|
const fesa::Vector& physical_p) {
|
|
fesa::Vector vp(chart.p_columns.size());
|
|
fesa::Vector ap(chart.p_columns.size());
|
|
for (std::size_t i = 0U; i < physical_p.Size(); ++i) {
|
|
if (copy_translation_rates && reduced::PhysicalCoordinateChannel(chart, i) == 0U) {
|
|
vp[i] = retained_velocity[i];
|
|
ap[i] = retained_acceleration[i];
|
|
} else {
|
|
vp[i] = predictor.velocity[i] + p.gamma * h * trial_ar[i];
|
|
ap[i] = (trial_ar[i] + p.alpha * private_a[i]) / (1.0 + p.alpha);
|
|
}
|
|
}
|
|
auto root = Require(reduced::SolveEmbeddingRoot(chart, physical_p, seed,
|
|
[&](const fesa::Vector& q) { return position_operators(frame, q); }, tangent));
|
|
const auto rates = Require(reduced::ReconstructPhysicalRates(
|
|
chart, root, vp, ap, tangent));
|
|
auto candidate = frame;
|
|
candidate.generalized_coordinates = root.q;
|
|
candidate.velocity = rates.velocity;
|
|
candidate.acceleration = rates.acceleration;
|
|
for (std::size_t i = 0U; i < zero.Size(); ++i)
|
|
candidate.algorithmic_acceleration[i] =
|
|
(1.0 + p.alpha) * candidate.acceleration[i] - p.alpha * old.acceleration[i];
|
|
return Trial{Evaluate(fixture, candidate, h, p, true), std::move(root),
|
|
std::move(vp), std::move(ap), trial_ar};
|
|
};
|
|
auto trial = evaluate_trial(ar, candidate_p);
|
|
bool same_start = true;
|
|
fesa::Vector old_ar(ar.Size());
|
|
for (std::size_t i = 0U; i < ar.Size(); ++i) {
|
|
same_start = same_start && private_p[i] == predictor.coordinates[i];
|
|
old_ar[i] = (private_p[i] - predictor.coordinates[i]) / (p.beta * h * h);
|
|
}
|
|
if (!same_start) {
|
|
// Both candidates are lifted at the new endpoint time; cache the winner.
|
|
const auto old_trial = evaluate_trial(old_ar, private_p);
|
|
const auto residual_norm = [&](const Assembly& value) {
|
|
double norm = 0.0;
|
|
for (const auto full : fixture.dofs.FreeDofs())
|
|
norm = std::max(norm, std::abs(value.frame.residual[full]) *
|
|
(full % 6U < 3U ? fixture.length : 1.0));
|
|
if (!std::isfinite(norm)) throw std::runtime_error{"C08 nonfinite start residual."};
|
|
return norm;
|
|
};
|
|
if (residual_norm(old_trial.assembly) < residual_norm(trial.assembly)) {
|
|
trial = old_trial;
|
|
ar = old_ar;
|
|
candidate_p = private_p;
|
|
chosen_start_p = private_p;
|
|
candidate_p_low = fesa::Vector(ar.Size());
|
|
addition_low_sum = fesa::Vector(ar.Size());
|
|
candidate_p_radius = private_p_radius;
|
|
for (std::size_t i = 0U; i < ar.Size(); ++i)
|
|
candidate_p_pair_radius[i] = candidate_p_radius[i];
|
|
// The selected old-p object has exact zero initial displacement.
|
|
requested_sum = fesa::Vector(ar.Size());
|
|
requested_compensation = fesa::Vector(ar.Size());
|
|
requested_absolute_sum = fesa::Vector(ar.Size());
|
|
for (std::size_t i = 0U; i < ar.Size(); ++i) {
|
|
const double numerator = private_p[i] - predictor.coordinates[i];
|
|
ar_radius[i] = QuotientRadius(numerator,
|
|
AddRadius(numerator, private_p_radius[i],
|
|
predictor_coordinate_radius[i]),
|
|
p.beta * h * h, HhtCoordinateFactorRadius(p, h));
|
|
}
|
|
}
|
|
}
|
|
for (std::uint32_t iteration = 0U; iteration <= 100U; ++iteration) {
|
|
if (iteration > 0U) trial = evaluate_trial(ar, candidate_p);
|
|
assembly = trial.assembly;
|
|
frame = assembly.frame;
|
|
const auto time_rates = force_rates(frame.time);
|
|
const auto operators = Require(reduced::EvaluateConstraintOperators(
|
|
*fixture.domain, fixture.dofs, fixture.shells, assembly, chart,
|
|
time_rates[0], time_rates[1]));
|
|
// The trial radii begin with the immutable sparse p basis. The graph
|
|
// indicator appends z/vz/az uncertainty after its frozen endpoint solve.
|
|
fesa::Vector vp_radius(chart.p_columns.size());
|
|
fesa::Vector ap_radius(chart.p_columns.size());
|
|
for (std::size_t i = 0U; i < ar.Size(); ++i) {
|
|
if (copy_translation_rates &&
|
|
reduced::PhysicalCoordinateChannel(chart, i) == 0U) {
|
|
vp_radius[i] = retained_velocity_radius[i];
|
|
ap_radius[i] = retained_acceleration_radius[i];
|
|
continue;
|
|
}
|
|
vp_radius[i] = AddRadius(trial.vp[i], predictor_velocity_radius[i],
|
|
ProductRadius(p.gamma * h, ProductRadius(p.gamma, 0.0, h, 0.0),
|
|
ar[i], ar_radius[i]));
|
|
const double numerator = ar[i] + p.alpha * private_a[i];
|
|
ap_radius[i] = QuotientRadius(numerator,
|
|
AddRadius(numerator, ar_radius[i], ProductRadius(p.alpha, 0.0,
|
|
private_a[i], private_a_radius[i])),
|
|
1.0 + p.alpha, RoundoffRadius(1.0 + p.alpha));
|
|
}
|
|
trial_radii = {LiftBasisRadii(chart.p_columns, candidate_p,
|
|
candidate_p_radius, chart.full_dof_count),
|
|
LiftBasisRadii(chart.p_columns, trial.vp, vp_radius,
|
|
chart.full_dof_count),
|
|
LiftBasisRadii(chart.p_columns, trial.ap, ap_radius,
|
|
chart.full_dof_count),
|
|
LiftBasisRadii(chart.p_columns, ar, ar_radius,
|
|
chart.full_dof_count)};
|
|
if (correction_pass && frame.free_residual_metrics[0] <= 1.0e-10 &&
|
|
frame.free_residual_metrics[1] <= 1.0e-10) {
|
|
if (manufactured)
|
|
CompleteManufacturedConstraints(fixture, assembly, null_nodes,
|
|
&positive, &tangent, &frame);
|
|
else
|
|
CompleteOriginalConstraints(fixture, assembly, null_nodes, positive,
|
|
tangent, &frame, h, "accepted", &trial.root);
|
|
Require(fesa::hht_internal::CompleteEquilibrium(*fixture.domain,
|
|
fixture.dofs, assembly, &frame));
|
|
Require(reduced::CheckGraphRegularity(assembly, chart, operators, tangent));
|
|
if (!last) throw std::runtime_error{"C08 nontrivial step lacks N7 last solve."};
|
|
const auto full_indicator = GraphArithmeticIndicator(frame, old, chart,
|
|
operators, trial_radii, arithmetic_radii, p, &tangent);
|
|
RecordNoise(fixture, assembly, chart, operators, *last, p, full_indicator,
|
|
&tangent, &evidence);
|
|
auto recovery_assembly = assembly;
|
|
recovery_assembly.frame = frame;
|
|
evidence.recovered_acceleration.push_back(RecoverAcceleration(
|
|
fixture, recovery_assembly, &tangent, manufactured));
|
|
frame.iteration_count = iteration;
|
|
const auto update_equal = [](const double actual,
|
|
const std::initializer_list<double> terms) {
|
|
double value = 0.0;
|
|
double scale = 0.0;
|
|
for (const double term : terms) {
|
|
value += term;
|
|
scale += std::abs(term);
|
|
}
|
|
const double defect = std::abs(actual - value);
|
|
return std::isfinite(defect) && std::isfinite(scale) &&
|
|
(scale == 0.0 ? defect == 0.0 : defect / scale <= 1.0e-12);
|
|
};
|
|
struct UpdateTrace {
|
|
double ratio = 0.0;
|
|
double actual = 0.0;
|
|
double expected = 0.0;
|
|
double scale = 0.0;
|
|
std::size_t full = 0U;
|
|
double signed_sum = 0.0;
|
|
double compensation = 0.0;
|
|
};
|
|
std::array<UpdateTrace, 3> update_trace{};
|
|
const auto trace_update = [&](const std::size_t check, const std::size_t full,
|
|
const double actual,
|
|
const std::initializer_list<double> terms) {
|
|
double expected = 0.0;
|
|
double scale = 0.0;
|
|
for (const double term : terms) {
|
|
expected += term;
|
|
scale += std::abs(term);
|
|
}
|
|
const double defect = std::abs(actual - expected);
|
|
const double ratio = scale == 0.0
|
|
? (defect == 0.0 ? 0.0 : std::numeric_limits<double>::infinity())
|
|
: defect / scale;
|
|
if (ratio >= update_trace[check].ratio)
|
|
update_trace[check] = {ratio, actual, expected, scale, full};
|
|
return update_equal(actual, terms);
|
|
};
|
|
for (std::size_t i = 0U; i < ar.Size(); ++i) {
|
|
const auto full = chart.p_dofs[i];
|
|
SCOPED_TRACE("increment=" + std::to_string(increment) + ", full=" + std::to_string(full));
|
|
EXPECT_TRUE(update_equal(trial.root.p[i],
|
|
{predictor.coordinates[i], p.beta * h * h * ar[i]}));
|
|
EXPECT_TRUE(update_equal(trial.vp[i],
|
|
{predictor.velocity[i], p.gamma * h * ar[i]}));
|
|
trace_update(0U, full, trial.root.p[i] - private_p[i],
|
|
{h * private_v[i], h * h * (0.5 - p.beta) * private_ar[i],
|
|
p.beta * h * h * ar[i]});
|
|
const double velocity_term = h * private_v[i];
|
|
const double old_rate_term = h * h * (0.5 - p.beta) * private_ar[i];
|
|
const double new_rate_term = p.beta * h * h * ar[i];
|
|
const auto whole_pair = reduced::CheckPairIncrement(trial.root.p[i],
|
|
candidate_p_low[i], private_p[i],
|
|
{velocity_term, old_rate_term, new_rate_term},
|
|
std::abs(velocity_term) + std::abs(old_rate_term) + std::abs(new_rate_term));
|
|
EXPECT_TRUE(whole_pair.storage_pass);
|
|
EXPECT_TRUE(whole_pair.arithmetic_pass) << whole_pair.ratio;
|
|
const double stored_increment = requested_sum[i] + requested_compensation[i];
|
|
const double actual_increment = trial.root.p[i] - private_p[i];
|
|
const double trace_defect = std::abs(actual_increment - stored_increment);
|
|
const double trace_scale = requested_absolute_sum[i];
|
|
const auto requested_pair = reduced::CheckPairIncrement(trial.root.p[i],
|
|
candidate_p_low[i], private_p[i],
|
|
{requested_sum[i], requested_compensation[i]}, trace_scale);
|
|
const double trace_ratio = trace_scale == 0.0
|
|
? (trace_defect == 0.0 ? 0.0 : std::numeric_limits<double>::infinity())
|
|
: trace_defect / trace_scale;
|
|
if (!(std::isfinite(trace_ratio) && trace_ratio <= 1.0e-12))
|
|
std::cout << std::setprecision(17) << "C08_POSITION_REPRESENTATION_DIAGNOSTIC h=" << h
|
|
<< " increment=" << increment << " full=" << full
|
|
<< " old_p=" << private_p[i] << " start_p=" << chosen_start_p[i]
|
|
<< " p=" << trial.root.p[i] << " actual=" << actual_increment
|
|
<< " expected=" << stored_increment << " signed_sum=" << requested_sum[i]
|
|
<< " compensation=" << requested_compensation[i] << " scale=" << trace_scale
|
|
<< " ratio=" << trace_ratio << " last_before=" << last_p_before_update[i]
|
|
<< " last_delta=" << last_requested_delta[i]
|
|
<< " addition_low_sum=" << addition_low_sum[i]
|
|
<< " candidate_low=" << candidate_p_low[i]
|
|
<< " pair_ratio=" << requested_pair.ratio
|
|
<< " pair_storage=" << requested_pair.storage_pass
|
|
<< " subtraction_low=" << AdditionLowPart(trial.root.p[i], -private_p[i], actual_increment)
|
|
<< "\n";
|
|
EXPECT_TRUE(requested_pair.storage_pass);
|
|
EXPECT_TRUE(requested_pair.arithmetic_pass) << requested_pair.ratio;
|
|
if (trace_ratio >= update_trace[1U].ratio)
|
|
update_trace[1U] = {trace_ratio, actual_increment, stored_increment, trace_scale, full,
|
|
requested_sum[i], requested_compensation[i]};
|
|
// This subtraction is diagnostic only; N3 requires the absolute
|
|
// velocity update and the resolved position increment.
|
|
trace_update(2U, full, trial.vp[i] - private_v[i],
|
|
{h * (1.0 - p.gamma) * private_ar[i], p.gamma * h * ar[i]});
|
|
EXPECT_TRUE(update_equal((1.0 + p.alpha) * trial.ap[i],
|
|
{ar[i], p.alpha * private_a[i]}));
|
|
}
|
|
if (increment <= 5U) {
|
|
for (std::size_t check = 0U; check < update_trace.size(); ++check) {
|
|
const auto& value = update_trace[check];
|
|
std::cout << std::setprecision(17) << "C08_RESOLVED_TRACE h=" << h
|
|
<< " increment=" << increment << " check=" << check
|
|
<< " full=" << value.full << " ratio=" << value.ratio
|
|
<< " actual=" << value.actual << " expected=" << value.expected
|
|
<< " scale=" << value.scale << " signed_sum=" << value.signed_sum
|
|
<< " compensation=" << value.compensation << "\n";
|
|
}
|
|
}
|
|
for (std::size_t i = 0U; i < zero.Size(); ++i)
|
|
EXPECT_DOUBLE_EQ(frame.algorithmic_acceleration[i],
|
|
(1.0 + p.alpha) * frame.acceleration[i] - p.alpha * old.acceleration[i]);
|
|
private_p = trial.root.p;
|
|
private_v = trial.vp;
|
|
private_a = trial.ap;
|
|
private_ar = ar;
|
|
private_p_radius = candidate_p_radius;
|
|
private_v_radius = vp_radius;
|
|
private_a_radius = ap_radius;
|
|
private_ar_radius = ar_radius;
|
|
seed = trial.root.z;
|
|
arithmetic_radii = full_indicator;
|
|
// Private carried-Ar remains a sparse p-coordinate history. Its full
|
|
// diagnostic lift is already represented by the frozen graph result;
|
|
// do not overwrite a plane support with a canonical label index.
|
|
converged = true;
|
|
break;
|
|
}
|
|
if (iteration == 100U) break;
|
|
const double cv = p.gamma / (p.beta * h);
|
|
const double ca = 1.0 / ((1.0 + p.alpha) * p.beta * h * h);
|
|
auto system = Require(reduced::AssembleAugmentedReducedSystem(
|
|
assembly, chart, operators, cv, ca, fixture.length, h));
|
|
const bool angular_only = reduced::ShouldRefineAngular(
|
|
chart, frame.free_residual_metrics, translation_corrections,
|
|
translation_private_algorithmic_pass);
|
|
auto selected = Require(reduced::SelectAugmentedSolve(system, chart,
|
|
angular_only));
|
|
fesa::Vector selected_solution(selected.rhs.Size());
|
|
Require(tangent.FactorizeFresh(selected.tangent));
|
|
Require(tangent.SolveCorrection(selected.rhs, selected_solution));
|
|
Require(fesa::hht_internal::CheckLinearSolve(selected.tangent,
|
|
selected.rhs,
|
|
selected_solution));
|
|
auto solution = Require(reduced::EmbedAugmentedSolution(
|
|
selected, selected_solution, system.rhs.Size()));
|
|
const auto correction = Require(reduced::MapAugmentedCorrection(
|
|
chart, system, solution, cv, ca));
|
|
correction_pass = true;
|
|
std::array<double, 4U> next_translation_corrections{};
|
|
for (const auto full : fixture.dofs.FreeDofs()) {
|
|
const double scale = full % 6U < 3U ? fixture.length : 1.0;
|
|
const std::array<double, 4U> values{
|
|
correction.coordinate[full], h * correction.velocity[full],
|
|
h * h * correction.acceleration[full],
|
|
h * h * (1.0 + p.alpha) * correction.acceleration[full]};
|
|
for (std::size_t rate = 0U; rate < values.size(); ++rate) {
|
|
const double ratio = std::abs(values[rate]) / scale;
|
|
correction_pass = correction_pass && ratio <= 1.0e-10;
|
|
if (full % 6U < 3U) {
|
|
if (!std::isfinite(ratio))
|
|
next_translation_corrections[rate] =
|
|
std::numeric_limits<double>::infinity();
|
|
else
|
|
next_translation_corrections[rate] = std::max(
|
|
next_translation_corrections[rate], ratio);
|
|
}
|
|
}
|
|
}
|
|
for (std::size_t i = 0U; i < ar.Size(); ++i) {
|
|
const double delta = correction.p[i] / (p.beta * h * h);
|
|
const double scale = reduced::PhysicalCoordinateChannel(chart, i) == 0U ? fixture.length : 1.0;
|
|
const double algorithmic_ratio = std::abs(h * h * delta) / scale;
|
|
correction_pass = correction_pass && algorithmic_ratio <= 1.0e-10;
|
|
// Preserve the physical Newton increment even when pp + beta*h*h*Ar
|
|
// cancels large terms. Ar retains its own equivalent HHT recurrence.
|
|
if (selected.angular_only && reduced::PhysicalCoordinateChannel(chart, i) == 0U)
|
|
continue;
|
|
record_requested_increment(i, correction.p[i]);
|
|
last_p_before_update[i] = candidate_p[i];
|
|
last_requested_delta[i] = correction.p[i];
|
|
const double raw_sum = candidate_p[i] + correction.p[i];
|
|
addition_low_sum[i] += AdditionLowPart(candidate_p[i], correction.p[i], raw_sum);
|
|
AdvanceCoordinate(correction.p[i], &candidate_p[i], &candidate_p_low[i],
|
|
&candidate_p_pair_radius[i]);
|
|
candidate_p_radius[i] = UpSum(candidate_p_pair_radius[i],
|
|
std::abs(candidate_p_low[i]));
|
|
ar[i] += delta;
|
|
ar_radius[i] = AddRadius(ar[i], ar_radius[i], QuotientRadius(
|
|
correction.p[i], 0.0, p.beta * h * h, HhtCoordinateFactorRadius(p, h)));
|
|
}
|
|
translation_corrections = next_translation_corrections;
|
|
// Retain the original private k*delta_ap=delta_Ar identity. The
|
|
// sparse Lp projection below is an additional basis-lift evidence, not
|
|
// a replacement for this recurrence check.
|
|
const double private_acceleration_rate =
|
|
1.0 / ((1.0 + p.alpha) * p.beta * h * h);
|
|
fesa::Vector delta_algorithmic_acceleration = correction.acceleration;
|
|
delta_algorithmic_acceleration.Scale(1.0 + p.alpha);
|
|
const auto projected_delta_algorithmic = Require(reduced::ProjectPhysical(
|
|
chart, delta_algorithmic_acceleration));
|
|
for (std::size_t i = 0U; i < ar.Size(); ++i) {
|
|
const double delta_ar = correction.p[i] / (p.beta * h * h);
|
|
const double private_delta_a = private_acceleration_rate * correction.p[i];
|
|
EXPECT_TRUE(private_algorithmic_equal((1.0 + p.alpha) * private_delta_a,
|
|
delta_ar));
|
|
double lifted = 0.0;
|
|
double lifted_scale = std::abs(delta_ar);
|
|
for (const auto& entry : chart.p_columns[i].support) {
|
|
const double term =
|
|
entry.value * delta_algorithmic_acceleration[entry.dof];
|
|
lifted += term;
|
|
lifted_scale += std::abs(term);
|
|
}
|
|
const double lifted_defect = std::abs(lifted - delta_ar);
|
|
const bool lifted_pass = std::isfinite(lifted_defect) &&
|
|
std::isfinite(lifted_scale) &&
|
|
(lifted_scale == 0.0 ? lifted_defect == 0.0
|
|
: lifted_defect / lifted_scale <= 1.0e-12);
|
|
EXPECT_TRUE(lifted_pass);
|
|
// This computed projection is diagnostic only: it includes rounded
|
|
// plane/algebraic leakage and must not replace private FB-08 history.
|
|
EXPECT_TRUE(std::isfinite(projected_delta_algorithmic[i]));
|
|
if (reduced::PhysicalCoordinateChannel(chart, i) != 0U) continue;
|
|
const double ratio = std::abs(h * h * delta_ar) / fixture.length;
|
|
translation_private_algorithmic_pass =
|
|
translation_private_algorithmic_pass && std::isfinite(ratio) &&
|
|
ratio <= 1.0e-10 && private_algorithmic_equal(
|
|
(1.0 + p.alpha) * private_delta_a, delta_ar);
|
|
}
|
|
copy_translation_rates = selected.angular_only;
|
|
if (copy_translation_rates) {
|
|
retained_velocity = trial.vp;
|
|
retained_acceleration = trial.ap;
|
|
retained_velocity_radius = vp_radius;
|
|
retained_acceleration_radius = ap_radius;
|
|
}
|
|
seed = trial.root.z;
|
|
for (std::size_t i = 0U; i < seed.Size(); ++i)
|
|
seed[i] += correction.z[i];
|
|
last.emplace(LastSolve{std::move(system), std::move(selected),
|
|
std::move(selected_solution), std::move(solution)});
|
|
}
|
|
if (!converged)
|
|
throw std::runtime_error{"C08 reduced Newton exhausted 100 corrections at h=" +
|
|
std::to_string(h) + ", increment=" + std::to_string(increment)};
|
|
const auto physical_oracle = Integrate(fixture, frame);
|
|
EXPECT_LE(std::abs(frame.kinetic_energy -
|
|
static_cast<double>(physical_oracle.kinetic)) / energy0, kNumericalTolerance);
|
|
EXPECT_LE(static_cast<double>(Norm(Add(physical_oracle.momentum,
|
|
Scale(L3{frame.linear_momentum[0], frame.linear_momentum[1],
|
|
frame.linear_momentum[2]}, -1.0L)))) /
|
|
static_cast<double>(initial_oracle.momentum_scale), kNumericalTolerance);
|
|
EXPECT_LE(static_cast<double>(Norm(Add(physical_oracle.angular_momentum,
|
|
Scale(L3{frame.angular_momentum[0], frame.angular_momentum[1],
|
|
frame.angular_momentum[2]}, -1.0L)))) /
|
|
static_cast<double>(initial_oracle.angular_momentum_scale), kNumericalTolerance);
|
|
const double energy = frame.kinetic_energy + frame.physical_strain_energy +
|
|
frame.drilling_energy;
|
|
if (!manufactured)
|
|
evidence.energy_error =
|
|
std::max(evidence.energy_error, std::abs(energy - energy0) / energy0);
|
|
double dp = 0.0;
|
|
double dl = 0.0;
|
|
for (std::size_t c = 0; c < 3; ++c) {
|
|
dp = std::hypot(dp, frame.linear_momentum[c] - initial_momentum[c]);
|
|
dl = std::hypot(dl,
|
|
frame.angular_momentum[c] - initial_angular_momentum[c]);
|
|
}
|
|
if (!manufactured) {
|
|
evidence.momentum_error =
|
|
std::max(evidence.momentum_error,
|
|
dp / static_cast<double>(initial_oracle.momentum_scale));
|
|
evidence.angular_momentum_error = std::max(
|
|
evidence.angular_momentum_error,
|
|
dl / static_cast<double>(initial_oracle.angular_momentum_scale));
|
|
}
|
|
for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node)
|
|
evidence.maximum_rotation = std::max(
|
|
evidence.maximum_rotation,
|
|
static_cast<double>(
|
|
Norm(Triple(frame.generalized_coordinates, 6U * node + 3U))));
|
|
for (std::size_t node = 0; node < fixture.domain->Nodes().size(); ++node)
|
|
for (std::size_t channel = 0; channel < 2U; ++channel) {
|
|
const auto offset = 6U * node + 3U * channel;
|
|
const double scale = channel == 0U ? fixture.length : 1.0;
|
|
evidence.algorithmic_acceleration_gap = std::max(
|
|
evidence.algorithmic_acceleration_gap,
|
|
evidence.end_time * evidence.end_time / scale *
|
|
static_cast<double>(Norm(
|
|
Add(Triple(frame.algorithmic_acceleration, offset),
|
|
Scale(Triple(frame.acceleration, offset), -1.0L)))));
|
|
}
|
|
for (std::size_t i = 0; i < frame.hidden_constraint.Size(); ++i) {
|
|
if (frame.hidden_constraint_scale[i] == 0.0)
|
|
EXPECT_DOUBLE_EQ(frame.hidden_constraint[i], 0.0);
|
|
else
|
|
EXPECT_LE(std::abs(frame.hidden_constraint[i]) /
|
|
frame.hidden_constraint_scale[i],
|
|
kNumericalTolerance);
|
|
const double vscale = frame.hidden_velocity_defect_scale[i];
|
|
const double ascale = frame.hidden_acceleration_defect_scale[i];
|
|
if (vscale > 0.0)
|
|
evidence.hidden_velocity =
|
|
std::max(evidence.hidden_velocity,
|
|
std::abs(frame.hidden_velocity_defect[i]) / vscale);
|
|
else
|
|
EXPECT_DOUBLE_EQ(frame.hidden_velocity_defect[i], 0.0);
|
|
if (ascale > 0.0)
|
|
evidence.hidden_acceleration =
|
|
std::max(evidence.hidden_acceleration,
|
|
std::abs(frame.hidden_acceleration_defect[i]) / ascale);
|
|
else
|
|
EXPECT_DOUBLE_EQ(frame.hidden_acceleration_defect[i], 0.0);
|
|
}
|
|
evidence.history.push_back(frame);
|
|
}
|
|
std::cout << std::setprecision(17) << (manufactured ? "C08a" : "C08b")
|
|
<< " h=" << h;
|
|
if (!manufactured)
|
|
std::cout << " energy=" << evidence.energy_error
|
|
<< " P=" << evidence.momentum_error
|
|
<< " L=" << evidence.angular_momentum_error;
|
|
std::cout << " phi_max=" << evidence.maximum_rotation
|
|
<< " a_A_gap=" << evidence.algorithmic_acceleration_gap
|
|
<< " hidden_v=" << evidence.hidden_velocity
|
|
<< " hidden_a=" << evidence.hidden_acceleration << " solve_noise=";
|
|
for (const double value : evidence.solve_noise) std::cout << value << ',';
|
|
std::cout << " arithmetic_scope=frozen-operator-first-order-indicator round_noise=";
|
|
for (const double value : evidence.arithmetic_noise)
|
|
std::cout << value << ',';
|
|
std::cout << " selector_history_radii=";
|
|
for (const double value : evidence.selector_history_arithmetic_noise)
|
|
std::cout << value << ',';
|
|
std::cout << " public_A_indicator=";
|
|
for (const double value : evidence.full_filter_arithmetic_noise)
|
|
std::cout << value << ',';
|
|
std::cout << '\n';
|
|
return evidence;
|
|
}
|
|
|
|
std::array<double, 6> Differences(const GridEvidence& coarse,
|
|
const GridEvidence& fine,
|
|
const double length) {
|
|
std::array<double, 6> error{};
|
|
for (std::size_t frame = 0; frame < coarse.history.size(); ++frame) {
|
|
const auto& a = coarse.history[frame];
|
|
const auto& b = fine.history[static_cast<std::size_t>(
|
|
std::llround(coarse.h / fine.h)) * frame];
|
|
const std::array<const fesa::Vector*, 3> first{
|
|
&a.generalized_coordinates, &a.velocity, &a.acceleration};
|
|
const std::array<const fesa::Vector*, 3> second{
|
|
&b.generalized_coordinates, &b.velocity, &b.acceleration};
|
|
for (std::size_t node = 0; node < a.nodal_internal_order.size(); ++node) {
|
|
for (std::size_t channel = 0; channel < 2; ++channel) {
|
|
const auto offset = 6U * node + 3U * channel;
|
|
for (std::size_t rate = 0; rate < 3; ++rate) {
|
|
const double normalization =
|
|
(rate == 0 ? 1.0
|
|
: rate == 1 ? coarse.end_time
|
|
: coarse.end_time * coarse.end_time) /
|
|
(channel == 0U ? length : 1.0);
|
|
const double difference = static_cast<double>(
|
|
Norm(Add(Triple(*first[rate], offset),
|
|
Scale(Triple(*second[rate], offset), -1.0L))));
|
|
error[2U * rate + channel] =
|
|
std::max(error[2U * rate + channel], normalization * difference);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return error;
|
|
}
|
|
|
|
std::array<double, 6> FrameError(const Frame& actual, const Frame& exact,
|
|
const double length, const double end_time) {
|
|
std::array<double, 6> error{};
|
|
const std::array<const fesa::Vector*, 3> first{
|
|
&actual.generalized_coordinates, &actual.velocity, &actual.acceleration};
|
|
const std::array<const fesa::Vector*, 3> second{
|
|
&exact.generalized_coordinates, &exact.velocity, &exact.acceleration};
|
|
for (std::size_t node = 0; node < actual.generalized_coordinates.Size() / 6U;
|
|
++node) {
|
|
for (std::size_t channel = 0; channel < 2U; ++channel) {
|
|
const auto offset = 6U * node + 3U * channel;
|
|
for (std::size_t rate = 0; rate < 3U; ++rate) {
|
|
const double factor = (rate == 0U ? 1.0
|
|
: rate == 1U ? end_time
|
|
: end_time * end_time) /
|
|
(channel == 0U ? length : 1.0);
|
|
error[2U * rate + channel] = std::max(
|
|
error[2U * rate + channel],
|
|
factor * static_cast<double>(Norm(
|
|
Add(Triple(*first[rate], offset),
|
|
Scale(Triple(*second[rate], offset), -1.0L)))));
|
|
}
|
|
}
|
|
}
|
|
return error;
|
|
}
|
|
|
|
std::array<double, 6> RecoveredFrameError(
|
|
const Frame& actual, const fesa::Vector& recovered_acceleration,
|
|
const Frame& exact, const double length, const double end_time) {
|
|
auto recovered = actual;
|
|
recovered.acceleration = recovered_acceleration;
|
|
return FrameError(recovered, exact, length, end_time);
|
|
}
|
|
|
|
std::array<double, 6> RecoveredDifferences(const GridEvidence& coarse,
|
|
const GridEvidence& fine,
|
|
const double length) {
|
|
std::array<double, 6> error{};
|
|
for (std::size_t frame = 0; frame < coarse.history.size(); ++frame) {
|
|
const auto& a = coarse.history[frame];
|
|
const auto& b = fine.history[static_cast<std::size_t>(
|
|
std::llround(coarse.h / fine.h)) * frame];
|
|
const std::array<const fesa::Vector*, 3> first{
|
|
&a.generalized_coordinates, &a.velocity,
|
|
&coarse.recovered_acceleration[frame].value};
|
|
const std::array<const fesa::Vector*, 3> second{
|
|
&b.generalized_coordinates, &b.velocity,
|
|
&fine.recovered_acceleration[static_cast<std::size_t>(
|
|
std::llround(coarse.h / fine.h)) * frame].value};
|
|
for (std::size_t node = 0; node < a.nodal_internal_order.size(); ++node) {
|
|
for (std::size_t channel = 0; channel < 2; ++channel) {
|
|
const auto offset = 6U * node + 3U * channel;
|
|
const double scale = channel == 0U ? length : 1.0;
|
|
for (std::size_t rate = 0; rate < 3; ++rate) {
|
|
const double normalization =
|
|
(rate == 0U ? 1.0
|
|
: rate == 1U ? coarse.end_time
|
|
: coarse.end_time * coarse.end_time) /
|
|
scale;
|
|
error[2U * rate + channel] = std::max(
|
|
error[2U * rate + channel],
|
|
normalization *
|
|
static_cast<double>(
|
|
Norm(Add(Triple(*first[rate], offset),
|
|
Scale(Triple(*second[rate], offset), -1.0L)))));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return error;
|
|
}
|
|
|
|
std::array<double, 2> AccelerationSplit(const Frame& basis,
|
|
const fesa::Vector& first,
|
|
const fesa::Vector& second,
|
|
const double end_time) {
|
|
const auto dofs = basis.generalized_coordinates.Size();
|
|
if (dofs % 6U != 0U || basis.velocity.Size() != dofs ||
|
|
basis.acceleration.Size() != dofs || first.Size() != dofs ||
|
|
second.Size() != dofs) {
|
|
throw std::runtime_error{
|
|
"C08 acceleration split requires matching nodal vectors."};
|
|
}
|
|
std::array<double, 2> split{};
|
|
for (std::size_t node = 0; node < dofs / 6U; ++node) {
|
|
const auto offset = 6U * node + 3U;
|
|
auto vector3 = [](const L3& value) {
|
|
return fesa::Vector3{static_cast<double>(value[0]),
|
|
static_cast<double>(value[1]),
|
|
static_cast<double>(value[2])};
|
|
};
|
|
const auto rotation = Require(fesa::EvaluateDynamicRotation(
|
|
vector3(Triple(basis.generalized_coordinates, offset)), {0.0, 0.0, 1.0},
|
|
vector3(Triple(basis.velocity, offset)),
|
|
vector3(Triple(basis.acceleration, offset))));
|
|
const L3 null{rotation.drill_null[0], rotation.drill_null[1],
|
|
rotation.drill_null[2]};
|
|
const L3 difference =
|
|
Add(Triple(first, offset), Scale(Triple(second, offset), -1.0L));
|
|
split[0] = std::max(
|
|
split[0], end_time * end_time *
|
|
static_cast<double>(std::abs(Dot(null, difference))));
|
|
split[1] = std::max(split[1],
|
|
end_time * end_time *
|
|
static_cast<double>(Norm(Cross(null, difference))));
|
|
}
|
|
return split;
|
|
}
|
|
|
|
TEST(HhtShellPhysics, ReducedManufacturedControllerKeepsExactInitialOrbit) {
|
|
const auto fixture = MakeFixture(1U);
|
|
UseLegacyDrilling(fixture);
|
|
const auto grid = RunPlate(fixture, 1.0 / 64.0, true);
|
|
ASSERT_EQ(grid.history.size(), 65U);
|
|
const auto exact = RigidOrbit(fixture, 0.0);
|
|
for (std::size_t i = 0U; i < exact.generalized_coordinates.Size(); ++i) {
|
|
EXPECT_DOUBLE_EQ(grid.history.front().generalized_coordinates[i], exact.generalized_coordinates[i]);
|
|
EXPECT_DOUBLE_EQ(grid.history.front().velocity[i], exact.velocity[i]);
|
|
EXPECT_DOUBLE_EQ(grid.history.front().acceleration[i], exact.acceleration[i]);
|
|
EXPECT_DOUBLE_EQ(grid.history.front().algorithmic_acceleration[i], exact.acceleration[i]);
|
|
}
|
|
for (std::size_t index = 0U; index < grid.history.size(); ++index) {
|
|
const auto& frame = grid.history[index];
|
|
EXPECT_DOUBLE_EQ(frame.time, static_cast<double>(index) / 64.0);
|
|
EXPECT_EQ(frame.time_integration_formulation,
|
|
fesa::DynamicTimeIntegrationFormulation::kReducedCoordinateHht);
|
|
EXPECT_EQ(frame.generalized_coordinates.Size(), fixture.dofs.FullDofCount());
|
|
EXPECT_EQ(frame.velocity.Size(), fixture.dofs.FullDofCount());
|
|
EXPECT_EQ(frame.acceleration.Size(), fixture.dofs.FullDofCount());
|
|
for (const auto* channel : {&frame.generalized_coordinates, &frame.velocity,
|
|
&frame.acceleration, &frame.algorithmic_acceleration})
|
|
for (std::size_t i = 0U; i < channel->Size(); ++i)
|
|
EXPECT_TRUE(std::isfinite((*channel)[i]));
|
|
for (const double residual : frame.free_residual_metrics)
|
|
EXPECT_LE(residual, kNumericalTolerance);
|
|
for (const auto* defect : {&frame.hidden_constraint,
|
|
&frame.hidden_velocity_defect,
|
|
&frame.hidden_acceleration_defect})
|
|
for (std::size_t i = 0U; i < defect->Size(); ++i)
|
|
EXPECT_TRUE(std::isfinite((*defect)[i]));
|
|
}
|
|
EXPECT_EQ(grid.history.back().time_integration_formulation,
|
|
fesa::DynamicTimeIntegrationFormulation::kReducedCoordinateHht);
|
|
EXPECT_DOUBLE_EQ(grid.history.back().time, 1.0);
|
|
EXPECT_GT(grid.maximum_rotation, 0.5);
|
|
}
|
|
|
|
TEST(HhtShellPhysics, C08aRigidManufacturedOrbitReducedCoordinateTimeConvergence) {
|
|
const auto fixture = MakeFixture(1U);
|
|
UseLegacyDrilling(fixture);
|
|
std::vector<GridEvidence> grids;
|
|
std::vector<std::array<double, 6>> errors;
|
|
for (const unsigned denominator : {16U, 32U, 64U, 128U, 256U, 512U, 1024U}) {
|
|
SCOPED_TRACE(denominator);
|
|
grids.push_back(
|
|
RunPlate(fixture, 1.0 / static_cast<double>(denominator), true));
|
|
const auto& grid = grids.back();
|
|
EXPECT_GT(grid.maximum_rotation, 0.5);
|
|
std::array<double, 6> maximum{};
|
|
std::array<double, 6> first{};
|
|
std::array<double, 6> late{};
|
|
std::array<double, 6> final{};
|
|
std::array<double, 6> maximum_time{};
|
|
std::array<double, 6> recovered_maximum{};
|
|
std::array<double, 6> recovered_maximum_time{};
|
|
std::array<double, 2> recovered_acceleration_split{};
|
|
double recovery_pre_hidden2 = 0.0;
|
|
double recovery_post_hidden2 = 0.0;
|
|
double recovery_mass_delta_a = 0.0;
|
|
double angular_acceleration_null_error = 0.0;
|
|
double angular_acceleration_physical_error = 0.0;
|
|
for (std::size_t index = 0; index < grid.history.size(); ++index) {
|
|
const auto& frame = grid.history[index];
|
|
ASSERT_EQ(frame.hidden_velocity_defect.Size(), fixture.domain->Nodes().size());
|
|
ASSERT_EQ(frame.hidden_velocity_defect_scale.Size(), frame.hidden_velocity_defect.Size());
|
|
ASSERT_EQ(frame.hidden_acceleration_defect.Size(), frame.hidden_velocity_defect.Size());
|
|
ASSERT_EQ(frame.hidden_acceleration_defect_scale.Size(), frame.hidden_velocity_defect.Size());
|
|
for (std::size_t node = 0U; node < frame.hidden_velocity_defect.Size(); ++node) {
|
|
for (const auto values : {
|
|
std::make_pair(frame.hidden_velocity_defect[node], frame.hidden_velocity_defect_scale[node]),
|
|
std::make_pair(frame.hidden_acceleration_defect[node], frame.hidden_acceleration_defect_scale[node])}) {
|
|
ASSERT_TRUE(std::isfinite(values.first));
|
|
ASSERT_TRUE(std::isfinite(values.second));
|
|
ASSERT_GE(values.second, 0.0);
|
|
if (values.second == 0.0)
|
|
EXPECT_DOUBLE_EQ(values.first, 0.0);
|
|
else
|
|
EXPECT_LE(std::abs(values.first) / values.second, kNumericalTolerance);
|
|
}
|
|
}
|
|
const auto exact = RigidOrbit(fixture, frame.time);
|
|
const auto error =
|
|
FrameError(frame, exact, fixture.length, grid.end_time);
|
|
const auto recovered_error =
|
|
RecoveredFrameError(frame, grid.recovered_acceleration[index].value,
|
|
exact, fixture.length, grid.end_time);
|
|
const auto split =
|
|
AccelerationSplit(exact, grid.recovered_acceleration[index].value,
|
|
exact.acceleration, grid.end_time);
|
|
for (std::size_t channel = 0; channel < 6U; ++channel) {
|
|
if (recovered_error[channel] > recovered_maximum[channel])
|
|
recovered_maximum_time[channel] = frame.time;
|
|
recovered_maximum[channel] =
|
|
std::max(recovered_maximum[channel], recovered_error[channel]);
|
|
}
|
|
for (std::size_t part = 0; part < 2U; ++part)
|
|
recovered_acceleration_split[part] =
|
|
std::max(recovered_acceleration_split[part], split[part]);
|
|
recovery_pre_hidden2 =
|
|
std::max(recovery_pre_hidden2,
|
|
grid.recovered_acceleration[index].pre_hidden2_relative);
|
|
recovery_post_hidden2 =
|
|
std::max(recovery_post_hidden2,
|
|
grid.recovered_acceleration[index].post_hidden2_relative);
|
|
recovery_mass_delta_a =
|
|
std::max(recovery_mass_delta_a,
|
|
grid.recovered_acceleration[index].mass_delta_a_relative);
|
|
if (error[5] > maximum[5]) {
|
|
double largest = 0.0;
|
|
for (std::size_t node = 0; node < fixture.domain->Nodes().size();
|
|
++node) {
|
|
const auto offset = 6U * node + 3U;
|
|
const auto difference =
|
|
Add(Triple(frame.acceleration, offset),
|
|
Scale(Triple(exact.acceleration, offset), -1.0L));
|
|
const double magnitude = static_cast<double>(Norm(difference));
|
|
if (magnitude <= largest) continue;
|
|
largest = magnitude;
|
|
const auto phi = Triple(frame.generalized_coordinates, offset);
|
|
const auto velocity = Triple(frame.velocity, offset);
|
|
const auto acceleration = Triple(frame.acceleration, offset);
|
|
auto vector3 = [](const L3& value) {
|
|
return fesa::Vector3{static_cast<double>(value[0]),
|
|
static_cast<double>(value[1]),
|
|
static_cast<double>(value[2])};
|
|
};
|
|
const auto rotation = Require(fesa::EvaluateDynamicRotation(
|
|
vector3(phi), {0.0, 0.0, 1.0}, vector3(velocity),
|
|
vector3(acceleration)));
|
|
const L3 null{rotation.drill_null[0], rotation.drill_null[1],
|
|
rotation.drill_null[2]};
|
|
angular_acceleration_null_error =
|
|
static_cast<double>(std::abs(Dot(null, difference)));
|
|
angular_acceleration_physical_error =
|
|
static_cast<double>(Norm(Cross(null, difference)));
|
|
}
|
|
}
|
|
for (std::size_t channel = 0; channel < 6U; ++channel) {
|
|
if (error[channel] > maximum[channel])
|
|
maximum_time[channel] = frame.time;
|
|
maximum[channel] = std::max(maximum[channel], error[channel]);
|
|
if (index == 1U) first[channel] = error[channel];
|
|
if (frame.time >= grid.end_time / 4.0)
|
|
late[channel] = std::max(late[channel], error[channel]);
|
|
}
|
|
final = error;
|
|
}
|
|
errors.push_back(maximum);
|
|
std::cout << "C08a temporal_error h=" << grid.h;
|
|
for (const auto& row :
|
|
{std::make_pair(" first=", first), std::make_pair(" late=", late),
|
|
std::make_pair(" maximum=", maximum),
|
|
std::make_pair(" final=", final),
|
|
std::make_pair(" maximum_time=", maximum_time)}) {
|
|
std::cout << row.first;
|
|
for (const double value : row.second) std::cout << value << ',';
|
|
}
|
|
std::cout << " max_angular_a_null=" << angular_acceleration_null_error
|
|
<< " max_angular_a_physical="
|
|
<< angular_acceleration_physical_error << '\n';
|
|
std::cout << "C08a recovered_diagnostic h=" << grid.h << " max_time_error=";
|
|
for (const double value : recovered_maximum) std::cout << value << ',';
|
|
std::cout << " max_time=";
|
|
for (const double value : recovered_maximum_time) std::cout << value << ',';
|
|
std::cout << " angular_a_null=" << recovered_acceleration_split[0]
|
|
<< " angular_a_physical=" << recovered_acceleration_split[1]
|
|
<< " pre_hidden2=" << recovery_pre_hidden2
|
|
<< " post_hidden2=" << recovery_post_hidden2
|
|
<< " Mdeltaa=" << recovery_mass_delta_a << '\n';
|
|
}
|
|
for (std::size_t channel = 0; channel < 6U; ++channel) {
|
|
SCOPED_TRACE(channel);
|
|
EXPECT_LE(errors.back()[channel], kNumericalTolerance);
|
|
for (std::size_t grid = grids.size() - 3U; grid + 1U < grids.size(); ++grid) {
|
|
if (errors[grid][channel] == 0.0 && errors[grid + 1U][channel] == 0.0)
|
|
continue;
|
|
ASSERT_GT(errors[grid][channel], 0.0);
|
|
ASSERT_GT(errors[grid + 1U][channel], 0.0);
|
|
EXPECT_GE(std::log2(errors[grid][channel] / errors[grid + 1U][channel]),
|
|
1.8);
|
|
}
|
|
for (std::size_t grid = grids.size() - 3U; grid < grids.size(); ++grid)
|
|
EXPECT_GT(errors[grid][channel],
|
|
100.0 * std::max(grids[grid].solve_noise[channel],
|
|
grids[grid].arithmetic_noise[channel]))
|
|
<< "N7 inconclusive: temporal error does not resolve "
|
|
"solve/arithmetic noise.";
|
|
}
|
|
}
|
|
|
|
TEST(HhtShellPhysics, C08bInitialPhysicalSpectrumProbe) {
|
|
const auto fixture = MakeFixture(2U);
|
|
const auto parameters = Require(fesa::MakeHhtParameters(-0.05));
|
|
const auto assembly =
|
|
Evaluate(fixture, EmptyFrame(fixture), 1.0 / 64.0, parameters, false);
|
|
const auto print_matrix = [](const char* prefix,
|
|
const fesa::SparseMatrix& matrix) {
|
|
ASSERT_EQ(matrix.Rows(), 54U);
|
|
ASSERT_EQ(matrix.Columns(), 54U);
|
|
ASSERT_EQ(matrix.RowOffsets().size(), 55U);
|
|
for (std::size_t row = 0; row < matrix.Rows(); ++row) {
|
|
std::array<double, 54U> values{};
|
|
for (std::size_t index = matrix.RowOffsets()[row];
|
|
index < matrix.RowOffsets()[row + 1U]; ++index) {
|
|
const auto column = matrix.ColumnIndices()[index];
|
|
ASSERT_LT(column, values.size());
|
|
values[column] = matrix.Values()[index];
|
|
}
|
|
std::cout << prefix << " row=" << row
|
|
<< " values=" << std::setprecision(17);
|
|
for (const double value : values) std::cout << value << ',';
|
|
std::cout << '\n';
|
|
}
|
|
};
|
|
// Diagnostic matrix dump for offline static condensation/eigensolution only.
|
|
print_matrix("C08b_initial_M", assembly.mass);
|
|
print_matrix("C08b_initial_K", assembly.internal_tangent);
|
|
}
|
|
|
|
TEST(HhtShellPhysics,
|
|
PredictorStartCompletesPreviouslyFailingLumpedFreePlateGrid) {
|
|
const auto fixture = MakeFixture(2U);
|
|
UseLegacyDrilling(fixture);
|
|
const auto evidence = RunPlate(fixture, 1.0 / 256.0, false);
|
|
ASSERT_EQ(evidence.history.size(), 65U);
|
|
EXPECT_DOUBLE_EQ(evidence.history.back().time, kEndTime);
|
|
EXPECT_GE(evidence.maximum_rotation, 0.2);
|
|
EXPECT_LT(evidence.maximum_rotation, std::acos(-1.0));
|
|
}
|
|
|
|
TEST(HhtShellPhysics,
|
|
C08bFreeFiniteRotationPlateSelfConvergenceAndPhysicalInvariants) {
|
|
const auto fixture = MakeFixture(2U);
|
|
UseLegacyDrilling(fixture);
|
|
std::vector<GridEvidence> grids;
|
|
for (const unsigned denominator :
|
|
{64U, 128U, 256U, 512U, 1024U, 32768U, 65536U, 131072U}) {
|
|
SCOPED_TRACE(denominator);
|
|
grids.push_back(
|
|
RunPlate(fixture, 1.0 / static_cast<double>(denominator), false));
|
|
EXPECT_GE(grids.back().maximum_rotation, 0.2);
|
|
}
|
|
std::vector<std::array<double, 6>> differences;
|
|
for (std::size_t pair = 0; pair + 1U < grids.size(); ++pair) {
|
|
differences.push_back(
|
|
Differences(grids[pair], grids[pair + 1U], fixture.length));
|
|
std::cout << "C08b self_difference h=" << grids[pair].h << ' ';
|
|
for (const double value : differences.back()) std::cout << value << ',';
|
|
std::cout << '\n';
|
|
const auto recovered =
|
|
RecoveredDifferences(grids[pair], grids[pair + 1U], fixture.length);
|
|
std::array<double, 2> acceleration_split{};
|
|
double recovery_pre_hidden2 = 0.0;
|
|
double recovery_post_hidden2 = 0.0;
|
|
double recovery_mass_delta_a = 0.0;
|
|
for (std::size_t frame = 0; frame < grids[pair].history.size(); ++frame) {
|
|
const auto fine_frame = static_cast<std::size_t>(
|
|
std::llround(grids[pair].h / grids[pair + 1U].h)) * frame;
|
|
const auto split = AccelerationSplit(
|
|
grids[pair + 1U].history[fine_frame],
|
|
grids[pair].recovered_acceleration[frame].value,
|
|
grids[pair + 1U].recovered_acceleration[fine_frame].value,
|
|
grids[pair].end_time);
|
|
for (std::size_t part = 0; part < 2U; ++part)
|
|
acceleration_split[part] =
|
|
std::max(acceleration_split[part], split[part]);
|
|
}
|
|
for (const auto* grid : {&grids[pair], &grids[pair + 1U]}) {
|
|
for (const auto& recovery : grid->recovered_acceleration) {
|
|
recovery_pre_hidden2 =
|
|
std::max(recovery_pre_hidden2, recovery.pre_hidden2_relative);
|
|
recovery_post_hidden2 =
|
|
std::max(recovery_post_hidden2, recovery.post_hidden2_relative);
|
|
recovery_mass_delta_a =
|
|
std::max(recovery_mass_delta_a, recovery.mass_delta_a_relative);
|
|
}
|
|
}
|
|
std::cout << "C08b recovered_diagnostic h=" << grids[pair].h
|
|
<< " common_time_difference=";
|
|
for (const double value : recovered) std::cout << value << ',';
|
|
std::cout << " angular_a_null=" << acceleration_split[0]
|
|
<< " angular_a_physical=" << acceleration_split[1]
|
|
<< " pre_hidden2=" << recovery_pre_hidden2
|
|
<< " post_hidden2=" << recovery_post_hidden2
|
|
<< " Mdeltaa=" << recovery_mass_delta_a << '\n';
|
|
}
|
|
for (std::size_t channel = 0; channel < 6; ++channel) {
|
|
SCOPED_TRACE(channel);
|
|
EXPECT_LE(differences.back()[channel], kNumericalTolerance);
|
|
for (std::size_t pair = differences.size() - 3U;
|
|
pair + 1U < differences.size(); ++pair) {
|
|
const double first = differences[pair][channel];
|
|
const double second = differences[pair + 1U][channel];
|
|
if (first == 0.0 && second == 0.0) continue;
|
|
ASSERT_GT(first, 0.0);
|
|
ASSERT_GT(second, 0.0);
|
|
// For the preserved 1024 -> 32768 gap use the unequal-grid
|
|
// Richardson ratio; the final 2:1 triple reduces to log2(d0/d1).
|
|
const double coarse_ratio = grids[pair].h / grids[pair + 1U].h;
|
|
const double fine_ratio = grids[pair + 1U].h / grids[pair + 2U].h;
|
|
const double threshold = (std::pow(coarse_ratio, 1.8) - 1.0) /
|
|
(1.0 - std::pow(fine_ratio, -1.8));
|
|
EXPECT_GE(first / second, threshold);
|
|
}
|
|
for (std::size_t pair = 1U; pair < differences.size(); ++pair) {
|
|
const double noise = std::max(grids[pair].solve_noise[channel],
|
|
grids[pair].arithmetic_noise[channel]) +
|
|
std::max(grids[pair + 1U].solve_noise[channel],
|
|
grids[pair + 1U].arithmetic_noise[channel]);
|
|
EXPECT_GT(differences[pair][channel], 100.0 * noise)
|
|
<< "N7 inconclusive: temporal difference does not resolve "
|
|
"solve/arithmetic noise.";
|
|
}
|
|
}
|
|
for (const auto member :
|
|
{&GridEvidence::energy_error, &GridEvidence::angular_momentum_error}) {
|
|
const auto final = grids.size() - 1U;
|
|
EXPECT_LT(grids[final - 1U].*member, grids[final - 2U].*member);
|
|
EXPECT_LT(grids[final].*member, grids[final - 1U].*member);
|
|
EXPECT_LE(grids[final].*member, 0.6 * (grids[final - 1U].*member));
|
|
EXPECT_LE(grids[final].*member, 1.0e-3);
|
|
}
|
|
// In this centered, all-free fixture the exact HHT branch conserves total
|
|
// linear momentum at every h. Endpoint equilibrium and CheckUpdates above
|
|
// exercise its balance/update identities; temporal strict decrease is not a
|
|
// defined P accuracy measure when the exact error is zero.
|
|
EXPECT_LE(grids.back().momentum_error, 1.0e-3);
|
|
}
|
|
|
|
TEST(HhtReducedIncrement, FrozenNearestCoordinatesPassPairArithmeticWithoutChangingHigh) {
|
|
const auto first = reduced::CheckPairIncrement(
|
|
.071707935756794142, -2.6003877354211524e-18, .071709031890658609,
|
|
{-1.0961338644697638e-6, -2.0834558726391122e-22},
|
|
1.0961338644697645e-6);
|
|
EXPECT_TRUE(first.storage_pass);
|
|
EXPECT_TRUE(first.arithmetic_pass);
|
|
EXPECT_EQ(first.defect, 0.0);
|
|
const auto worst = reduced::CheckPairIncrement(
|
|
.07170874594859544, -3.911862569787138e-18, .07170906968934286,
|
|
{-3.23740747422167e-7, -5.266112996573883e-23},
|
|
3.237407474221672e-7);
|
|
EXPECT_TRUE(worst.storage_pass);
|
|
EXPECT_TRUE(worst.arithmetic_pass);
|
|
EXPECT_EQ(worst.defect, 0.0);
|
|
// The independently logged requested trace is not reconstructed from high.
|
|
const auto missing_low = reduced::CheckPairIncrement(
|
|
.07170874594859544, 0.0, .07170906968934286,
|
|
{-3.23740747422167e-7, -5.266112996573883e-23},
|
|
3.237407474221672e-7);
|
|
EXPECT_FALSE(missing_low.arithmetic_pass);
|
|
}
|
|
|
|
TEST(HhtReducedIncrement, StorageRequiresNearestEvenIncludingBothMidpointDirections) {
|
|
const double half_upper_spacing = std::ldexp(1.0, -53);
|
|
const double odd = std::nextafter(1.0, 2.0);
|
|
EXPECT_TRUE(reduced::CheckPairIncrement(1.0, half_upper_spacing, 1.0,
|
|
{half_upper_spacing}, half_upper_spacing).storage_pass);
|
|
EXPECT_TRUE(reduced::CheckPairIncrement(1.0, -std::ldexp(1.0, -54), 1.0,
|
|
{-std::ldexp(1.0, -54)}, std::ldexp(1.0, -54)).storage_pass);
|
|
EXPECT_TRUE(reduced::CheckPairIncrement(-1.0, -half_upper_spacing, -1.0,
|
|
{-half_upper_spacing}, half_upper_spacing).storage_pass);
|
|
const auto odd_tie = reduced::CheckPairIncrement(odd, -half_upper_spacing,
|
|
1.0, {half_upper_spacing}, half_upper_spacing);
|
|
EXPECT_TRUE(odd_tie.arithmetic_pass);
|
|
EXPECT_FALSE(odd_tie.storage_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(-odd, half_upper_spacing, -1.0,
|
|
{-half_upper_spacing}, half_upper_spacing).storage_pass);
|
|
const auto misplaced_high = reduced::CheckPairIncrement(1.0,
|
|
std::ldexp(1.0, -52), 1.0, {std::ldexp(1.0, -52)}, std::ldexp(1.0, -52));
|
|
EXPECT_TRUE(misplaced_high.arithmetic_pass);
|
|
EXPECT_FALSE(misplaced_high.storage_pass);
|
|
}
|
|
|
|
TEST(HhtReducedIncrement, WholeUpdateAndIndependentTraceRejectMissingCorrection) {
|
|
const double correction = std::ldexp(1.0, -54);
|
|
const auto valid_whole = reduced::CheckPairIncrement(1.0, 2.0 * correction,
|
|
1.0, {correction, 0.0, correction}, 2.0 * correction);
|
|
const auto valid_trace = reduced::CheckPairIncrement(1.0, 2.0 * correction,
|
|
1.0, {2.0 * correction, 0.0}, 2.0 * correction);
|
|
EXPECT_TRUE(valid_whole.arithmetic_pass && valid_whole.storage_pass);
|
|
EXPECT_TRUE(valid_trace.arithmetic_pass && valid_trace.storage_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(1.0, correction, 1.0,
|
|
{correction, 0.0, correction}, 2.0 * correction).arithmetic_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(1.0, correction, 1.0,
|
|
{2.0 * correction, 0.0}, 2.0 * correction).arithmetic_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(1.0, -2.0 * correction, 1.0,
|
|
{2.0 * correction, 0.0}, 2.0 * correction).arithmetic_pass);
|
|
}
|
|
|
|
TEST(HhtReducedIncrement, ExactThresholdRejectsTailHiddenByRoundedRatio) {
|
|
const double threshold = 1.0e-12;
|
|
const double tail = std::ldexp(1.0, -100);
|
|
const auto equality = reduced::CheckPairIncrement(0.0, 0.0, 0.0,
|
|
{threshold}, 1.0);
|
|
EXPECT_TRUE(equality.arithmetic_pass);
|
|
const auto above = reduced::CheckPairIncrement(0.0, 0.0, 0.0,
|
|
{threshold, tail}, 1.0);
|
|
EXPECT_EQ(above.ratio, threshold);
|
|
EXPECT_FALSE(above.arithmetic_pass);
|
|
const auto above_decimal = reduced::CheckPairIncrement(0.0, 0.0, 0.0,
|
|
{threshold, std::ldexp(1.0, -94)}, 1.0);
|
|
EXPECT_EQ(above_decimal.ratio, threshold);
|
|
EXPECT_FALSE(above_decimal.arithmetic_pass);
|
|
const auto below = reduced::CheckPairIncrement(0.0, 0.0, 0.0,
|
|
{threshold, -tail}, 1.0);
|
|
EXPECT_TRUE(below.arithmetic_pass);
|
|
const double largest = std::numeric_limits<double>::max();
|
|
const auto maximum = reduced::CheckPairIncrement(0.0, 0.0, 0.0,
|
|
{1.0}, largest);
|
|
EXPECT_TRUE(maximum.arithmetic_pass);
|
|
const double tiny = std::numeric_limits<double>::denorm_min();
|
|
const auto unsupported = reduced::CheckPairIncrement(0.0, 0.0, 0.0,
|
|
{tiny}, std::ldexp(1.0, -1000));
|
|
EXPECT_FALSE(unsupported.arithmetic_pass);
|
|
}
|
|
|
|
TEST(HhtReducedIncrement, ZeroScaleAndSubnormalArithmeticRemainExact) {
|
|
const double tiny = std::numeric_limits<double>::denorm_min();
|
|
const auto zero = reduced::CheckPairIncrement(1.0, 0.0, 1.0, {0.0}, 0.0);
|
|
EXPECT_TRUE(zero.arithmetic_pass && zero.storage_pass);
|
|
const auto fixed_zero = reduced::CheckPairIncrement(0.0, 0.0, 0.0, {0.0}, 0.0);
|
|
EXPECT_TRUE(fixed_zero.arithmetic_pass && fixed_zero.storage_pass);
|
|
const auto tiny_valid = reduced::CheckPairIncrement(tiny, 0.0, 0.0, {tiny}, tiny);
|
|
EXPECT_TRUE(tiny_valid.arithmetic_pass && tiny_valid.storage_pass);
|
|
EXPECT_EQ(tiny_valid.defect, 0.0);
|
|
const auto tiny_nonzero = reduced::CheckPairIncrement(1.0, tiny, 1.0, {0.0}, 0.0);
|
|
EXPECT_TRUE(tiny_nonzero.storage_pass);
|
|
EXPECT_FALSE(tiny_nonzero.arithmetic_pass);
|
|
EXPECT_EQ(tiny_nonzero.defect, tiny);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(0.0, tiny, 0.0, {tiny}, tiny).storage_pass);
|
|
}
|
|
|
|
TEST(HhtReducedIncrement, NonfiniteAndOverflowFailBeforeAcceptance) {
|
|
const double nan = std::numeric_limits<double>::quiet_NaN();
|
|
const double infinity = std::numeric_limits<double>::infinity();
|
|
for (const auto invalid : {nan, infinity}) {
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(invalid, 0.0, 0.0, {0.0}, 0.0).arithmetic_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(0.0, invalid, 0.0, {0.0}, 0.0).storage_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(0.0, 0.0, invalid, {0.0}, 0.0).arithmetic_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(0.0, 0.0, 0.0, {invalid}, 1.0).arithmetic_pass);
|
|
EXPECT_FALSE(reduced::CheckPairIncrement(0.0, 0.0, 0.0, {0.0}, invalid).arithmetic_pass);
|
|
}
|
|
const double largest = std::numeric_limits<double>::max();
|
|
const auto overflow = reduced::CheckPairIncrement(largest, largest, 0.0,
|
|
{largest}, largest);
|
|
EXPECT_FALSE(overflow.arithmetic_pass);
|
|
EXPECT_FALSE(overflow.storage_pass);
|
|
#if defined(_M_X64) || defined(__SSE2__)
|
|
const auto saved_mxcsr = _mm_getcsr();
|
|
for (const auto mask : {0x8000U, 0x0040U}) {
|
|
_mm_setcsr(saved_mxcsr | mask);
|
|
const auto unsupported = reduced::CheckPairIncrement(1.0,
|
|
std::numeric_limits<double>::denorm_min(), 1.0, {0.0}, 0.0);
|
|
_mm_setcsr(saved_mxcsr);
|
|
EXPECT_FALSE(unsupported.arithmetic_pass);
|
|
EXPECT_FALSE(unsupported.storage_pass);
|
|
}
|
|
#endif
|
|
}
|
|
|
|
TEST(HhtReducedLifecycle, AllowsEvaluationAfterCorrection100ButNeverCorrection101) {
|
|
using fesa::hht_reduced_internal::CanApplyOuterCorrection;
|
|
using fesa::hht_reduced_internal::kMaximumOuterCorrections;
|
|
EXPECT_TRUE(CanApplyOuterCorrection(0U));
|
|
EXPECT_TRUE(CanApplyOuterCorrection(99U));
|
|
EXPECT_FALSE(CanApplyOuterCorrection(100U));
|
|
EXPECT_FALSE(CanApplyOuterCorrection(101U));
|
|
// The actual executor evaluates acceptance before this correction guard.
|
|
// This verifies its budget policy, not a physical 100-correction trajectory.
|
|
for (const bool accept_at_100 : {false, true}) {
|
|
std::uint32_t corrections = 0U;
|
|
std::uint32_t last_evaluation = 0U;
|
|
bool accepted = false;
|
|
for (std::uint32_t evaluation = 0U;
|
|
evaluation <= kMaximumOuterCorrections; ++evaluation) {
|
|
last_evaluation = evaluation;
|
|
if (accept_at_100 && evaluation == 100U) {
|
|
accepted = true;
|
|
break;
|
|
}
|
|
if (!CanApplyOuterCorrection(evaluation)) break;
|
|
++corrections;
|
|
}
|
|
EXPECT_EQ(corrections, 100U);
|
|
EXPECT_EQ(last_evaluation, 100U);
|
|
EXPECT_EQ(accepted, accept_at_100);
|
|
}
|
|
}
|
|
|
|
std::string ReadLifecycleLog(const std::filesystem::path& path) {
|
|
std::ifstream stream(path);
|
|
std::ostringstream text;
|
|
text << stream.rdbuf();
|
|
return text.str();
|
|
}
|
|
|
|
class FailAfterAcceptedIncrement final : public fesa::TangentSolver {
|
|
public:
|
|
explicit FailAfterAcceptedIncrement(std::filesystem::path path)
|
|
: path_(std::move(path)) {}
|
|
fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
|
|
std::istringstream lines(ReadLifecycleLog(path_));
|
|
std::string line;
|
|
while (std::getline(lines, line)) {
|
|
if (line.find("\tincrement=1\t") != std::string::npos &&
|
|
line.find("\toutcome=accepted") != std::string::npos) {
|
|
++failures;
|
|
return fesa::Status::Failure(fesa::FailureCategory::kSolver,
|
|
{{fesa::Severity::kError, "c10-injected-backend", {}, "DYNAMIC", {},
|
|
"Injected backend failure after accepted increment one."}});
|
|
}
|
|
}
|
|
return solver_.FactorizeFresh(matrix);
|
|
}
|
|
fesa::Status SolveCorrection(const fesa::Vector& rhs,
|
|
fesa::Vector& correction) const override {
|
|
return solver_.SolveCorrection(rhs, correction);
|
|
}
|
|
std::size_t failures = 0U;
|
|
private:
|
|
std::filesystem::path path_;
|
|
fesa::MklPardisoTangentSolver solver_;
|
|
};
|
|
|
|
TEST(HhtReducedLifecycle, FailedExecutionDoesNotPublishCandidateOrContaminateRetry) {
|
|
auto fixture = MakeFixture(1U, true);
|
|
const fesa::SerialParallelFor parallel;
|
|
const auto run = [&](fesa::TangentSolver& tangent,
|
|
const std::filesystem::path& path) {
|
|
fesa::MklPardisoSolver positive;
|
|
fesa::NewtonIterationLogWriter log;
|
|
Require(log.Open(path));
|
|
auto result = fesa::ExecuteReducedHht(*fixture.domain, log, parallel,
|
|
positive, tangent);
|
|
Require(log.Close());
|
|
return result;
|
|
};
|
|
const auto directory = std::filesystem::temp_directory_path();
|
|
const auto baseline_path = directory / "fesa-c10-baseline.log";
|
|
const auto failure_path = directory / "fesa-c10-failure.log";
|
|
const auto retry_path = directory / "fesa-c10-retry.log";
|
|
fesa::MklPardisoTangentSolver baseline_solver;
|
|
auto baseline = run(baseline_solver, baseline_path);
|
|
ASSERT_TRUE(baseline.HasValue()) << Message(baseline.GetStatus());
|
|
ASSERT_EQ(baseline.Value().state.AcceptedFrames().size(), 3U);
|
|
EXPECT_GT(baseline.Value().state.AcceptedFrames()[1U].generalized_coordinates[0U],
|
|
0.0);
|
|
FailAfterAcceptedIncrement failing(failure_path);
|
|
const auto failed = run(failing, failure_path);
|
|
ASSERT_FALSE(failed.HasValue());
|
|
EXPECT_EQ(failing.failures, 1U);
|
|
EXPECT_EQ(failed.GetStatus().Diagnostics().front().code, "c10-injected-backend");
|
|
const auto partial_log = ReadLifecycleLog(failure_path);
|
|
EXPECT_NE(partial_log.find("\tincrement=2\t"), std::string::npos);
|
|
EXPECT_NE(partial_log.find("\toutcome=failed"), std::string::npos);
|
|
std::istringstream lines(partial_log);
|
|
std::string line;
|
|
while (std::getline(lines, line)) {
|
|
if (line.find("\toutcome=initial") != std::string::npos ||
|
|
line.find("\toutcome=accepted") != std::string::npos) {
|
|
EXPECT_NE(ReadLifecycleLog(baseline_path).find(line + '\n'), std::string::npos);
|
|
EXPECT_EQ(line.find("\tincrement=2\t"), std::string::npos);
|
|
}
|
|
}
|
|
fesa::MklPardisoTangentSolver retry_solver;
|
|
auto retry = run(retry_solver, retry_path);
|
|
ASSERT_TRUE(retry.HasValue()) << Message(retry.GetStatus());
|
|
const auto& expected = baseline.Value().state.AcceptedFrames();
|
|
const auto& actual = retry.Value().state.AcceptedFrames();
|
|
ASSERT_EQ(actual.size(), expected.size());
|
|
for (std::size_t frame = 0U; frame < actual.size(); ++frame) {
|
|
EXPECT_EQ(actual[frame].frame_index, expected[frame].frame_index);
|
|
EXPECT_EQ(actual[frame].iteration_count, expected[frame].iteration_count);
|
|
EXPECT_EQ(actual[frame].time, expected[frame].time);
|
|
for (const auto member : {&Frame::generalized_coordinates, &Frame::velocity,
|
|
&Frame::acceleration, &Frame::algorithmic_acceleration,
|
|
&Frame::internal_force, &Frame::inertia_force, &Frame::residual,
|
|
&Frame::reaction, &Frame::hidden_constraint,
|
|
&Frame::hidden_velocity_defect, &Frame::hidden_acceleration_defect}) {
|
|
const auto& left = actual[frame].*member;
|
|
const auto& right = expected[frame].*member;
|
|
ASSERT_EQ(left.Size(), right.Size());
|
|
for (std::size_t i = 0U; i < left.Size(); ++i)
|
|
EXPECT_EQ(left[i], right[i]);
|
|
}
|
|
}
|
|
EXPECT_EQ(ReadLifecycleLog(retry_path), ReadLifecycleLog(baseline_path));
|
|
for (const auto& path : {baseline_path, failure_path, retry_path})
|
|
std::filesystem::remove(path);
|
|
// Result<DynamicExecutionResult> contains no value after failure. Private
|
|
// histories are local to that execution; they are not directly observable.
|
|
}
|
|
|
|
} // namespace
|