KOKO\Mimi
1517 lines
66 KiB
C++
1517 lines
66 KiB
C++
#include "fesa/analysis/solid_hht_dynamic_procedure.h"
|
|
|
|
#include <gtest/gtest.h>
|
|
|
|
#include <cfenv>
|
|
#include <cmath>
|
|
#include <filesystem>
|
|
#include <iostream>
|
|
#include <limits>
|
|
|
|
#include "analysis/solid_hht_force_identity.h"
|
|
#include "analysis/solid_hht_high_low.h"
|
|
#include "analysis/solid_hht_representability.h"
|
|
#include "fesa/assembly/parallel_for.h"
|
|
#include "fesa/elements/element_factory.h"
|
|
#include "fesa/elements/hex8_element.h"
|
|
#include "fesa/io/abaqus/domain_mapper.h"
|
|
#include "fesa/io/abaqus/input_reader.h"
|
|
#include "fesa/math/sparse_matrix.h"
|
|
#include "fesa/solvers/linear/mkl_pardiso_solver.h"
|
|
#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h"
|
|
#include "solid_hht_test_fixture.h"
|
|
|
|
#include "fesa/analysis/analysis_model.h"
|
|
|
|
namespace {
|
|
fesa::Result<fesa::SolidDynamicExecutionResult> RunDynamics(
|
|
fesa::ModelDefinition d) {
|
|
auto domain = fesa::Domain::Create(std::move(d));
|
|
if (!domain.HasValue())
|
|
return fesa::Result<fesa::SolidDynamicExecutionResult>::Failure(
|
|
domain.GetStatus());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
fesa::MklPardisoTangentSolver tangent;
|
|
return fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.Execute(domain.Value());
|
|
}
|
|
fesa::Status Injected() {
|
|
return fesa::Status::Failure(fesa::FailureCategory::kSolver,
|
|
{{fesa::Severity::kError,
|
|
"injected-backend",
|
|
{},
|
|
"DYNAMIC",
|
|
{},
|
|
"Injected backend failure."}});
|
|
}
|
|
class ControlledTangent : public fesa::TangentSolver {
|
|
public:
|
|
enum Mode { kFail, kNan, kZero, kInvert };
|
|
explicit ControlledTangent(Mode mode, unsigned successful_solves = 0)
|
|
: mode_(mode), successful_solves_(successful_solves) {}
|
|
fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
|
|
size_ = matrix.Rows();
|
|
return real_.FactorizeFresh(matrix);
|
|
}
|
|
fesa::Status SolveCorrection(const fesa::Vector& rhs,
|
|
fesa::Vector& du) const override {
|
|
if (calls_++ < successful_solves_) return real_.SolveCorrection(rhs, du);
|
|
if (mode_ == kFail) return Injected();
|
|
du = fesa::Vector(size_);
|
|
if (mode_ == kNan) du[0] = std::numeric_limits<double>::quiet_NaN();
|
|
if (mode_ == kInvert) du[0] = 2.0;
|
|
return fesa::Status::Ok();
|
|
}
|
|
|
|
private:
|
|
Mode mode_;
|
|
unsigned successful_solves_;
|
|
mutable unsigned calls_{0};
|
|
std::size_t size_{0};
|
|
fesa::MklPardisoTangentSolver real_;
|
|
};
|
|
class FailedMass : public fesa::LinearSolver {
|
|
public:
|
|
explicit FailedMass(bool fail_factorization)
|
|
: factorization_(fail_factorization) {}
|
|
fesa::Status Factorize(const fesa::SparseMatrix&) override {
|
|
return factorization_ ? Injected() : fesa::Status::Ok();
|
|
}
|
|
fesa::Status Solve(const fesa::Vector&, fesa::Vector&) const override {
|
|
return Injected();
|
|
}
|
|
|
|
private:
|
|
bool factorization_;
|
|
};
|
|
fesa::ModelDefinition OneFree(bool nonlinear = false, double h = .01,
|
|
unsigned count = 20) {
|
|
auto d = solid_hht_test::Brick(h, count, nonlinear);
|
|
d.dynamic_step->boundaries.push_back({"N0", 2, 3, 0, {}});
|
|
for (unsigned node = 1; node < 8; ++node)
|
|
d.dynamic_step->boundaries.push_back(
|
|
{"N" + std::to_string(node), 1, 3, 0, {}});
|
|
d.dynamic_step->loads = {{"N0", 1, .01, {}}};
|
|
return d;
|
|
}
|
|
} // namespace
|
|
|
|
namespace {
|
|
// Independent uniform-acceleration brick initial state: rho*V/8=1/4 at each
|
|
// node. High a=1 and low a=2^-60 therefore give Ma=(1/4,2^-62).
|
|
fesa::SolidDynamicAcceptedFrameState PairedUniformInitial(
|
|
const fesa::Domain& domain) {
|
|
fesa::SolidDynamicAcceptedFrameState frame;
|
|
frame.is_initial = true;
|
|
for (std::size_t n = 0; n < domain.Nodes().size(); ++n) {
|
|
frame.nodal_internal_order.push_back(static_cast<fesa::EntityIndex>(n));
|
|
frame.nodal_source_ids.push_back(domain.Nodes()[n].source_id);
|
|
}
|
|
for (auto* part :
|
|
{&frame.displacement, &frame.velocity,
|
|
&frame.acceleration, &frame.algorithmic_acceleration,
|
|
&frame.inertia_force, &frame.internal_force,
|
|
&frame.external_force, &frame.residual,
|
|
&frame.reaction, &frame.residual_scale,
|
|
&frame.low.displacement, &frame.low.velocity,
|
|
&frame.low.acceleration, &frame.low.algorithmic_acceleration,
|
|
&frame.low.inertia_force, &frame.low.internal_force,
|
|
&frame.low.external_force, &frame.low.residual,
|
|
&frame.low.reaction, &frame.low.residual_scale})
|
|
*part = fesa::Vector(24);
|
|
for (std::size_t n = 0; n < 8; ++n) {
|
|
const auto i = 3 * n;
|
|
frame.acceleration[i] = frame.algorithmic_acceleration[i] = 1;
|
|
frame.low.acceleration[i] = frame.low.algorithmic_acceleration[i] =
|
|
std::ldexp(1.0, -60);
|
|
frame.inertia_force[i] = frame.external_force[i] = .25;
|
|
frame.low.inertia_force[i] = frame.low.external_force[i] =
|
|
std::ldexp(1.0, -62);
|
|
frame.residual_scale[i] = .5;
|
|
frame.low.residual_scale[i] = std::ldexp(1.0, -61);
|
|
}
|
|
frame.linear_momentum_rate[0] = 2;
|
|
frame.low.linear_momentum_rate[0] = std::ldexp(1.0, -59);
|
|
frame.angular_momentum_rate = {0, 1, -1};
|
|
frame.low.angular_momentum_rate = {0, std::ldexp(1.0, -60),
|
|
-std::ldexp(1.0, -60)};
|
|
fesa::ElementFactory factory;
|
|
auto runtime = factory.Create(domain.Elements()[0], domain);
|
|
const auto* solid =
|
|
dynamic_cast<const fesa::Hex8Element*>(runtime.Value().get());
|
|
const auto zero = solid->EvaluateLinear(fesa::Vector(24));
|
|
frame.solid_results.assign(zero.Value().points.begin(),
|
|
zero.Value().points.end());
|
|
fesa::SolidDynamicIterationEvidence evidence;
|
|
evidence.residual_pass = evidence.correction_pass = evidence.update_pass =
|
|
true;
|
|
frame.iterations.push_back(evidence);
|
|
return frame;
|
|
}
|
|
} // namespace
|
|
|
|
TEST(SolidHhtPairedState, InitialLowPartsSurviveTrialCopyAndRejectedCommit) {
|
|
auto domain = fesa::Domain::Create(solid_hht_test::Brick(.01, 1));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
auto initial = PairedUniformInitial(domain.Value());
|
|
ASSERT_TRUE(state.Value().CommitAcceptedFrame(initial).IsOk());
|
|
auto trial = state.Value().BeginTrialFromCommitted();
|
|
ASSERT_TRUE(trial.HasValue());
|
|
EXPECT_EQ(trial.Value().low.acceleration[0], std::ldexp(1.0, -60));
|
|
EXPECT_EQ(trial.Value().low.residual_scale[0], std::ldexp(1.0, -61));
|
|
EXPECT_EQ(trial.Value().low.linear_momentum_rate[0], std::ldexp(1.0, -59));
|
|
EXPECT_EQ(trial.Value().low.angular_momentum_rate[1], std::ldexp(1.0, -60));
|
|
trial.Value().low.acceleration[0] = 0;
|
|
EXPECT_EQ(state.Value().CommittedState().low.acceleration[0],
|
|
std::ldexp(1.0, -60));
|
|
auto bad = initial;
|
|
bad.frame_index = bad.increment_index = 1;
|
|
bad.time = .01;
|
|
bad.is_initial = false;
|
|
bad.low.internal_force[0] = std::numeric_limits<double>::quiet_NaN();
|
|
EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk());
|
|
ASSERT_EQ(state.Value().AcceptedFrames().size(), 1U);
|
|
EXPECT_EQ(state.Value().CommittedState().low.acceleration[0],
|
|
std::ldexp(1.0, -60));
|
|
EXPECT_EQ(state.Value().CommittedState().low.inertia_force[0],
|
|
std::ldexp(1.0, -62));
|
|
EXPECT_EQ(state.Value().CommittedState().low.angular_momentum_rate[2],
|
|
-std::ldexp(1.0, -60));
|
|
}
|
|
|
|
TEST(SolidHhtPairedState, InvalidLowShapeFiniteNormalizationAndIdentitiesFail) {
|
|
auto domain = fesa::Domain::Create(solid_hht_test::Brick(.01, 1));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
const auto initial = PairedUniformInitial(domain.Value());
|
|
for (unsigned mutation = 0; mutation < 8; ++mutation) {
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
auto bad = initial;
|
|
switch (mutation) {
|
|
case 0:
|
|
bad.low.velocity = fesa::Vector(23);
|
|
break;
|
|
case 1:
|
|
bad.low.acceleration[0] = std::numeric_limits<double>::infinity();
|
|
break;
|
|
case 2:
|
|
bad.low.acceleration[0] = std::ldexp(1.0, -51);
|
|
break;
|
|
case 3:
|
|
bad.low.algorithmic_acceleration[0] = 0;
|
|
break;
|
|
case 4:
|
|
bad.low.residual[0] = std::ldexp(1.0, -70);
|
|
break;
|
|
case 5:
|
|
bad.low.residual_scale[0] = 0;
|
|
break;
|
|
case 6:
|
|
bad.low.reaction[0] = std::ldexp(1.0, -70);
|
|
break;
|
|
case 7:
|
|
bad.low.angular_momentum[0] = std::numeric_limits<double>::quiet_NaN();
|
|
break;
|
|
}
|
|
EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk()) << mutation;
|
|
EXPECT_TRUE(state.Value().AcceptedFrames().empty()) << mutation;
|
|
EXPECT_FALSE(state.Value().BeginTrialFromCommitted().HasValue())
|
|
<< mutation;
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtPairedKinematics, SubUlpPhysicalDisplacementRecoversNonzeroRates) {
|
|
const double epsilon = std::ldexp(1.0, -60);
|
|
auto p = fesa::MakeHhtParameters(0);
|
|
ASSERT_TRUE(p.HasValue());
|
|
fesa::SolidDynamicAcceptedFrameState old;
|
|
old.displacement = fesa::Vector(1, 1);
|
|
old.low.displacement = fesa::Vector(1, epsilon);
|
|
old.velocity = old.low.velocity = fesa::Vector(1);
|
|
old.acceleration = old.low.acceleration = fesa::Vector(1);
|
|
old.algorithmic_acceleration = old.low.algorithmic_acceleration =
|
|
fesa::Vector(1);
|
|
auto predictor = fesa::detail::PredictSolidHhtHighLow(old, 1, p.Value());
|
|
ASSERT_TRUE(predictor.HasValue());
|
|
EXPECT_EQ(predictor.Value().displacement[0], 1);
|
|
EXPECT_EQ(predictor.Value().displacement_low[0], epsilon);
|
|
auto rates = fesa::detail::EvaluateSolidHhtHighLowRates(
|
|
fesa::Vector(1, 1), fesa::Vector(1, 2 * epsilon), predictor.Value(),
|
|
old.acceleration, old.low.acceleration, 1, p.Value());
|
|
ASSERT_TRUE(rates.HasValue());
|
|
// beta=1/4 and gamma=1/2 give exact dyadic physical oracles.
|
|
EXPECT_EQ(rates.Value().algorithmic_acceleration[0], 4 * epsilon);
|
|
EXPECT_EQ(rates.Value().acceleration[0], 4 * epsilon);
|
|
EXPECT_EQ(rates.Value().velocity[0], 2 * epsilon);
|
|
EXPECT_EQ(rates.Value().algorithmic_acceleration_low[0], 0);
|
|
EXPECT_EQ(rates.Value().acceleration_low[0], 0);
|
|
EXPECT_EQ(rates.Value().velocity_low[0], 0);
|
|
auto rounded_predictor =
|
|
fesa::PredictHht(old.displacement, old.velocity,
|
|
old.algorithmic_acceleration, 1, p.Value());
|
|
auto rounded_rates =
|
|
fesa::EvaluateHhtRates(fesa::Vector(1, 1), rounded_predictor.Value(),
|
|
old.acceleration, 1, p.Value());
|
|
ASSERT_TRUE(rounded_rates.HasValue());
|
|
EXPECT_EQ(rounded_rates.Value().acceleration[0], 0);
|
|
}
|
|
|
|
TEST(SolidHhtPairedKinematics, PhysicalAndAuxiliaryRatesKeepOldPhysicalTerm) {
|
|
const double epsilon = std::ldexp(1.0, -60);
|
|
auto p = fesa::MakeHhtParameters(-.25);
|
|
ASSERT_TRUE(p.HasValue());
|
|
// h=8, beta=25/64 and gamma=3/4: denominator25, gamma*h=6.
|
|
fesa::detail::SolidHhtPairedPredictor predictor{
|
|
fesa::Vector(1, 1), fesa::Vector(1), fesa::Vector(1), fesa::Vector(1)};
|
|
auto rates = fesa::detail::EvaluateSolidHhtHighLowRates(
|
|
fesa::Vector(1, 1), fesa::Vector(1, 25 * epsilon), predictor,
|
|
fesa::Vector(1, 4 * epsilon), fesa::Vector(1), 8, p.Value());
|
|
ASSERT_TRUE(rates.HasValue());
|
|
EXPECT_EQ(rates.Value().algorithmic_acceleration[0], epsilon);
|
|
EXPECT_EQ(rates.Value().acceleration[0], 0);
|
|
EXPECT_EQ(rates.Value().velocity[0], 6 * epsilon);
|
|
}
|
|
|
|
TEST(SolidHhtPairedKinematics, InvalidLowInputsFailWithoutChangingPredictor) {
|
|
auto p = fesa::MakeHhtParameters(0);
|
|
ASSERT_TRUE(p.HasValue());
|
|
fesa::detail::SolidHhtPairedPredictor predictor{
|
|
fesa::Vector(1, 1), fesa::Vector(1), fesa::Vector(1), fesa::Vector(1)};
|
|
for (const double low :
|
|
{std::numeric_limits<double>::quiet_NaN(),
|
|
std::numeric_limits<double>::infinity(), std::ldexp(1.0, -51)}) {
|
|
EXPECT_FALSE(fesa::detail::EvaluateSolidHhtHighLowRates(
|
|
fesa::Vector(1, 1), fesa::Vector(1, low), predictor,
|
|
fesa::Vector(1), fesa::Vector(1), 1, p.Value())
|
|
.HasValue());
|
|
EXPECT_EQ(predictor.displacement[0], 1);
|
|
EXPECT_EQ(predictor.displacement_low[0], 0);
|
|
}
|
|
EXPECT_FALSE(fesa::detail::EvaluateSolidHhtHighLowRates(
|
|
fesa::Vector(1, 1), fesa::Vector(0), predictor,
|
|
fesa::Vector(1), fesa::Vector(1), 1, p.Value())
|
|
.HasValue());
|
|
auto invalid_parameters = p.Value();
|
|
invalid_parameters.beta = .5;
|
|
EXPECT_FALSE(fesa::detail::EvaluateSolidHhtHighLowRates(
|
|
fesa::Vector(1, 1), fesa::Vector(1), predictor,
|
|
fesa::Vector(1), fesa::Vector(1), 1, invalid_parameters)
|
|
.HasValue());
|
|
}
|
|
|
|
TEST(SolidHhtPairedKinematics, MassMultiplyKeepsSubUlpActionAndCancellation) {
|
|
fesa::SparsePattern pattern;
|
|
pattern.row_offsets = {0, 2, 4};
|
|
pattern.column_indices = {0, 1, 0, 1};
|
|
auto mass = fesa::SparseMatrix::FromCoo(
|
|
2, 2,
|
|
{{0, 0, 2, 0, 0}, {0, 1, 1, 0, 1}, {1, 0, 1, 0, 2}, {1, 1, 2, 0, 3}},
|
|
pattern);
|
|
ASSERT_TRUE(mass.HasValue());
|
|
const double epsilon = std::ldexp(1.0, -60);
|
|
fesa::Vector high(2), low(2);
|
|
high[0] = 1;
|
|
high[1] = -2;
|
|
low[0] = epsilon;
|
|
low[1] = epsilon;
|
|
auto result =
|
|
fesa::detail::MultiplySolidHhtMassHighLow(mass.Value(), high, low);
|
|
ASSERT_TRUE(result.HasValue());
|
|
EXPECT_EQ(result.Value().high[0], 3 * epsilon);
|
|
EXPECT_EQ(result.Value().low[0], 0);
|
|
EXPECT_EQ(result.Value().high[1], -3);
|
|
EXPECT_EQ(result.Value().low[1], 3 * epsilon);
|
|
EXPECT_EQ(mass.Value().Multiply(high)[0], 0);
|
|
}
|
|
|
|
TEST(SolidHhtPairedKinematics, WeightedIdentityPreservesNonzeroOldResidualLow) {
|
|
fesa::SparsePattern pattern;
|
|
pattern.row_offsets = {0, 1};
|
|
pattern.column_indices = {0};
|
|
auto mass = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 1, 0, 0}}, pattern);
|
|
ASSERT_TRUE(mass.HasValue());
|
|
const double epsilon = std::ldexp(1.0, -60);
|
|
fesa::SolidDynamicAcceptedFrameState old, next;
|
|
for (auto* frame : {&old, &next})
|
|
for (auto* v :
|
|
{&frame->acceleration, &frame->low.acceleration,
|
|
&frame->algorithmic_acceleration,
|
|
&frame->low.algorithmic_acceleration, &frame->internal_force,
|
|
&frame->low.internal_force, &frame->external_force,
|
|
&frame->low.external_force, &frame->residual, &frame->low.residual})
|
|
*v = fesa::Vector(1);
|
|
old.acceleration[0] = 4 * epsilon;
|
|
old.internal_force[0] = 3;
|
|
old.external_force[0] = 5;
|
|
old.residual[0] = -2;
|
|
old.low.residual[0] = 4 * epsilon;
|
|
next.algorithmic_acceleration[0] = epsilon;
|
|
next.internal_force[0] = 7;
|
|
next.external_force[0] = 2;
|
|
next.residual[0] = 5;
|
|
fesa::detail::SolidHhtPairedForceIdentity result;
|
|
ASSERT_TRUE(fesa::detail::EvaluateSolidHhtHighLowForceIdentity(
|
|
mass.Value(), 0, fesa::math_internal::HighLow(epsilon), next, old, -.25,
|
|
&result));
|
|
EXPECT_EQ(result.original.high, 3.25);
|
|
EXPECT_EQ(result.original.low, epsilon);
|
|
EXPECT_EQ(result.solved.high, 3.25);
|
|
EXPECT_EQ(result.solved.low, epsilon);
|
|
EXPECT_EQ(result.defect.high, 0);
|
|
EXPECT_EQ(result.metric.high, 0);
|
|
// Dropping only the old residual low changes the approved solved force.
|
|
old.low.residual[0] = 0;
|
|
ASSERT_TRUE(fesa::detail::EvaluateSolidHhtHighLowForceIdentity(
|
|
mass.Value(), 0, fesa::math_internal::HighLow(epsilon), next, old, -.25,
|
|
&result));
|
|
EXPECT_EQ(result.solved.low, 0);
|
|
EXPECT_EQ(result.defect.high, epsilon);
|
|
}
|
|
|
|
namespace {
|
|
class RestoreArithmeticEnvironment {
|
|
public:
|
|
RestoreArithmeticEnvironment() : rounding_(std::fegetround()) {
|
|
#if defined(_M_X64) || defined(__SSE2__)
|
|
mxcsr_ = _mm_getcsr();
|
|
#endif
|
|
}
|
|
~RestoreArithmeticEnvironment() {
|
|
std::fesetround(rounding_);
|
|
#if defined(_M_X64) || defined(__SSE2__)
|
|
_mm_setcsr(mxcsr_);
|
|
#endif
|
|
}
|
|
|
|
private:
|
|
int rounding_;
|
|
#if defined(_M_X64) || defined(__SSE2__)
|
|
unsigned mxcsr_;
|
|
#endif
|
|
};
|
|
void MutateArithmeticEnvironment(bool flush) {
|
|
#if defined(_M_X64) || defined(__SSE2__)
|
|
if (flush) {
|
|
_mm_setcsr(_mm_getcsr() | 0x8040U);
|
|
return;
|
|
}
|
|
#else
|
|
(void)flush;
|
|
#endif
|
|
std::fesetround(FE_UPWARD);
|
|
}
|
|
class EnvironmentMass final : public fesa::LinearSolver {
|
|
public:
|
|
explicit EnvironmentMass(bool flush) : flush_(flush) {}
|
|
fesa::Status Factorize(const fesa::SparseMatrix& matrix) override {
|
|
return real_.Factorize(matrix);
|
|
}
|
|
fesa::Status Solve(const fesa::Vector& rhs,
|
|
fesa::Vector& result) const override {
|
|
const auto status = real_.Solve(rhs, result);
|
|
if (status.IsOk()) MutateArithmeticEnvironment(flush_);
|
|
return status;
|
|
}
|
|
|
|
private:
|
|
bool flush_;
|
|
fesa::MklPardisoSolver real_;
|
|
};
|
|
class EnvironmentTangent final : public fesa::TangentSolver {
|
|
public:
|
|
explicit EnvironmentTangent(bool flush) : flush_(flush) {}
|
|
fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
|
|
return real_.FactorizeFresh(matrix);
|
|
}
|
|
fesa::Status SolveCorrection(const fesa::Vector& rhs,
|
|
fesa::Vector& result) const override {
|
|
const auto status = real_.SolveCorrection(rhs, result);
|
|
if (status.IsOk() && calls_++ >= 2) MutateArithmeticEnvironment(flush_);
|
|
return status;
|
|
}
|
|
|
|
private:
|
|
bool flush_;
|
|
mutable unsigned calls_{0};
|
|
fesa::MklPardisoTangentSolver real_;
|
|
};
|
|
} // namespace
|
|
|
|
TEST(SolidHhtPairedState, BackendEnvironmentMutationFailsAndPreservesPairs) {
|
|
const auto baseline = RunDynamics(OneFree(false, .01, 3));
|
|
ASSERT_TRUE(baseline.HasValue());
|
|
const auto& expected = baseline.Value().state.AcceptedFrames()[1];
|
|
for (bool flush : {false, true}) {
|
|
for (bool mutate_mass : {false, true}) {
|
|
SCOPED_TRACE(flush);
|
|
SCOPED_TRACE(mutate_mass);
|
|
auto domain = fesa::Domain::Create(OneFree(false, .01, 3));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver real_mass;
|
|
fesa::MklPardisoTangentSolver real_tangent;
|
|
EnvironmentMass changed_mass(flush);
|
|
EnvironmentTangent changed_tangent(flush);
|
|
auto status = fesa::Status::Ok();
|
|
{
|
|
RestoreArithmeticEnvironment restore;
|
|
status = mutate_mass ? fesa::SolidHhtDynamicProcedureExecutor(
|
|
parallel, changed_mass, real_tangent)
|
|
.ExecuteInto(domain.Value(), state.Value())
|
|
: fesa::SolidHhtDynamicProcedureExecutor(
|
|
parallel, real_mass, changed_tangent)
|
|
.ExecuteInto(domain.Value(), state.Value());
|
|
EXPECT_FALSE(fesa::math_internal::ArithmeticEnvironmentSupported());
|
|
}
|
|
ASSERT_FALSE(status.IsOk());
|
|
ASSERT_FALSE(status.Diagnostics().empty());
|
|
EXPECT_EQ(status.Diagnostics().front().code,
|
|
"solid-dynamic-arithmetic-environment");
|
|
if (mutate_mass) {
|
|
EXPECT_TRUE(state.Value().AcceptedFrames().empty());
|
|
} else {
|
|
ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U);
|
|
const auto& actual = state.Value().CommittedState();
|
|
EXPECT_EQ(actual.displacement[0], expected.displacement[0]);
|
|
EXPECT_EQ(actual.low.displacement[0], expected.low.displacement[0]);
|
|
EXPECT_EQ(actual.acceleration[0], expected.acceleration[0]);
|
|
EXPECT_EQ(actual.low.acceleration[0], expected.low.acceleration[0]);
|
|
EXPECT_EQ(actual.external_work, expected.external_work);
|
|
EXPECT_EQ(actual.low.external_work, expected.low.external_work);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics, AllConstrainedRetainsPhysicalReactionAndEightPoints) {
|
|
auto d = solid_hht_test::Brick();
|
|
d.dynamic_step->boundaries = {{"ALL", 1, 3, 0, {}}};
|
|
d.dynamic_step->loads = {{"N0", 2, 7, {}}};
|
|
auto domain = fesa::Domain::Create(std::move(d));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
fesa::MklPardisoTangentSolver tangent;
|
|
auto result = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.Execute(domain.Value());
|
|
ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus());
|
|
ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 21U);
|
|
for (const auto& f : result.Value().state.AcceptedFrames()) {
|
|
EXPECT_EQ(f.displacement.Size(), 24U);
|
|
EXPECT_EQ(f.solid_results.size(), 8U);
|
|
EXPECT_DOUBLE_EQ(f.reaction[1], -7);
|
|
EXPECT_DOUBLE_EQ(f.acceleration[1], 0);
|
|
EXPECT_DOUBLE_EQ(f.external_work, 0);
|
|
EXPECT_EQ(f.iteration_count, 0U);
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
ActualEndpointsUseDifferentCurvesKnotsInteriorsAndHolds) {
|
|
auto d = solid_hht_test::Brick(.125, 8);
|
|
d.dynamic_step->boundaries = {{"ALL", 1, 3, 0, {}}};
|
|
d.amplitudes = {{"A", "A", {}, {{0, 0}, {.375, 1}, {.75, 2}}},
|
|
{"B", "B", {}, {{0, 3}, {.25, -1}, {.5, 1}}}};
|
|
d.dynamic_step->loads = {{"N0", 1, 2, {}, 0, "A", {}},
|
|
{"N0", 1, 3, {}, 1, "B", {}},
|
|
{"N1", 3, 7, {}}};
|
|
auto result = RunDynamics(std::move(d));
|
|
ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus());
|
|
for (const auto& f : result.Value().state.AcceptedFrames()) {
|
|
const double a = f.time <= .375 ? f.time / .375
|
|
: f.time <= .75 ? 1 + (f.time - .375) / .375
|
|
: 2;
|
|
const double b = f.time <= .25 ? 3 - 16 * f.time
|
|
: f.time <= .5 ? -1 + 8 * (f.time - .25)
|
|
: 1;
|
|
EXPECT_NEAR(f.external_force[0], 2 * a + 3 * b, 1e-14);
|
|
EXPECT_DOUBLE_EQ(f.reaction[0], -f.external_force[0]);
|
|
EXPECT_DOUBLE_EQ(f.external_force[5], 7);
|
|
EXPECT_DOUBLE_EQ(f.time, static_cast<double>(f.increment_index) * .125);
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics, RampStartsAtZeroAndAuxiliaryAccelerationIsDistinct) {
|
|
auto d = OneFree();
|
|
d.amplitudes = {{"A", "A", {}, {{0, 0}, {.1, 1}, {.2, 1}}}};
|
|
d.dynamic_step->loads[0].amplitude_index = 0;
|
|
d.dynamic_step->loads[0].amplitude_name_raw = "A";
|
|
auto result = RunDynamics(std::move(d));
|
|
ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus());
|
|
const auto& frames = result.Value().state.AcceptedFrames();
|
|
EXPECT_DOUBLE_EQ(frames.front().acceleration[0], 0);
|
|
EXPECT_GT(std::abs(frames[1].acceleration[0] -
|
|
frames[1].algorithmic_acceleration[0]),
|
|
1e-5);
|
|
for (std::size_t n = 1; n < frames.size(); ++n) {
|
|
const auto& f = frames[n];
|
|
const auto& old = frames[n - 1];
|
|
EXPECT_NEAR(f.algorithmic_acceleration[0],
|
|
.95 * f.acceleration[0] + .05 * old.acceleration[0], 1e-14);
|
|
EXPECT_LT(f.iterations.back().update_metric, 1e-12);
|
|
EXPECT_LE(f.iterations.back().residual_metric, 1e-10);
|
|
for (std::size_t i = 1; i < 24; ++i) {
|
|
EXPECT_DOUBLE_EQ(f.displacement[i], 0);
|
|
EXPECT_DOUBLE_EQ(f.velocity[i], 0);
|
|
EXPECT_DOUBLE_EQ(f.acceleration[i], 0);
|
|
EXPECT_DOUBLE_EQ(f.algorithmic_acceleration[i], 0);
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
OneFreeTranslationHasAnalyticLumpedInertialSupportReactions) {
|
|
auto result = RunDynamics(OneFree());
|
|
ASSERT_TRUE(result.HasValue()) << solid_hht_test::Message(result.GetStatus());
|
|
const auto d = solid_hht_test::Brick();
|
|
constexpr double bulk = 8.0;
|
|
constexpr double mu = 4.8;
|
|
const double roundoff = 512 * std::numeric_limits<double>::epsilon();
|
|
const auto& frames = result.Value().state.AcceptedFrames();
|
|
for (std::size_t n = 0; n < frames.size(); ++n) {
|
|
const auto& f = frames[n];
|
|
for (unsigned b = 0; b < 8; ++b) {
|
|
std::array<double, 3> product{}, sign{};
|
|
double scalar_mass = b == 0 ? .25 : 0.0;
|
|
for (unsigned c = 0; c < 3; ++c) {
|
|
const bool same = d.nodes[b].coordinates[c] == 0;
|
|
product[c] = same ? 1.0 / 3 : 1.0 / 6;
|
|
sign[c] = same ? -1 : 1;
|
|
}
|
|
for (unsigned c = 0; c < 3; ++c) {
|
|
double stiffness;
|
|
double stiffness_term_scale;
|
|
if (c == 0) {
|
|
stiffness = (4 * mu / 3) * (-sign[0]) * product[1] * product[2] +
|
|
mu * (-sign[1]) * product[0] * product[2] +
|
|
mu * (-sign[2]) * product[0] * product[1];
|
|
stiffness_term_scale = (4 * mu / 3) * product[1] * product[2] +
|
|
mu * product[0] * product[2] +
|
|
mu * product[0] * product[1];
|
|
} else {
|
|
const unsigned other = c == 1 ? 2 : 1;
|
|
stiffness =
|
|
-.25 * ((-2 * mu / 3) * sign[c] + mu * sign[0]) * product[other];
|
|
stiffness_term_scale = .25 * (2 * mu / 3 + mu) * product[other];
|
|
}
|
|
stiffness -= bulk * sign[c] / 16;
|
|
stiffness_term_scale += bulk / 16;
|
|
const double inertia = c == 0 ? scalar_mass * f.acceleration[0] : 0;
|
|
const double internal = stiffness * f.displacement[0];
|
|
const double applied = b == 0 && c == 0 ? .01 : 0;
|
|
const double reaction =
|
|
b == 0 && c == 0 ? 0 : inertia + internal - applied;
|
|
const double scale =
|
|
std::abs(inertia) +
|
|
stiffness_term_scale * std::abs(f.displacement[0]) +
|
|
std::abs(applied);
|
|
EXPECT_NEAR(f.reaction[3 * b + c], reaction, roundoff * scale);
|
|
}
|
|
}
|
|
if (n != 0) {
|
|
const auto& old = frames[n - 1];
|
|
const double mass = .25;
|
|
const double weighted =
|
|
mass * f.algorithmic_acceleration[0] +
|
|
.95 * (f.internal_force[0] - f.external_force[0]) +
|
|
.05 * (old.internal_force[0] - old.external_force[0]);
|
|
const double endpoint_identity =
|
|
.95 * f.residual[0] + .05 * old.residual[0];
|
|
const double scale = std::abs(mass * f.algorithmic_acceleration[0]) +
|
|
.95 * (std::abs(f.internal_force[0]) +
|
|
std::abs(f.external_force[0])) +
|
|
.05 * (std::abs(old.internal_force[0]) +
|
|
std::abs(old.external_force[0]));
|
|
EXPECT_NEAR(weighted, endpoint_identity, roundoff * scale);
|
|
EXPECT_NEAR(f.iterations.back().weighted_residual_numerator,
|
|
std::abs(weighted), roundoff * scale);
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
BackendNanNonconvergenceAndInversionPreserveAcceptedInitial) {
|
|
for (const auto mode :
|
|
{ControlledTangent::kFail, ControlledTangent::kNan,
|
|
ControlledTangent::kZero, ControlledTangent::kInvert}) {
|
|
SCOPED_TRACE(mode);
|
|
auto d = OneFree(true);
|
|
// Moving node0 past the opposite x corner creates a genuinely folded cell.
|
|
auto domain = fesa::Domain::Create(std::move(d));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
ASSERT_TRUE(state.HasValue());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
ControlledTangent tangent(mode);
|
|
auto status =
|
|
fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteInto(domain.Value(), state.Value());
|
|
EXPECT_FALSE(status.IsOk());
|
|
ASSERT_EQ(state.Value().AcceptedFrames().size(), 1U);
|
|
EXPECT_EQ(state.Value().CommittedState().increment_index, 0U);
|
|
EXPECT_DOUBLE_EQ(state.Value().CommittedState().displacement[0], 0);
|
|
EXPECT_NEAR(state.Value().CommittedState().acceleration[0],
|
|
.01 / .25, 1e-14);
|
|
ASSERT_FALSE(state.Value().IterationRecords().empty());
|
|
EXPECT_EQ(state.Value().IterationRecords().back().stage, "failed");
|
|
EXPECT_FALSE(state.Value().IterationRecords().back().failure_code.empty());
|
|
if (mode == ControlledTangent::kZero) {
|
|
const auto& failed = state.Value().IterationRecords().back();
|
|
EXPECT_EQ(failed.evidence.iteration_index, 100U);
|
|
EXPECT_EQ(failed.failure_code, "solid-dynamic-nonconvergence");
|
|
// Paired physical evidence replaces the retired Float64 shadow sum.
|
|
EXPECT_FALSE(failed.failure_force_sum.has_value());
|
|
const double expected_inertia =
|
|
.01 * (-(.5 - .275625) / .275625 - .05) / .95;
|
|
EXPECT_NEAR(failed.evidence.worst_force_terms[1], expected_inertia,
|
|
32 * std::numeric_limits<double>::epsilon() *
|
|
std::abs(expected_inertia));
|
|
EXPECT_DOUBLE_EQ(failed.evidence.worst_force_terms[2], 0);
|
|
EXPECT_DOUBLE_EQ(failed.evidence.worst_force_terms[3], .01);
|
|
EXPECT_FALSE(failed.evidence.residual_pass);
|
|
EXPECT_TRUE(failed.evidence.update_pass);
|
|
EXPECT_TRUE(failed.evidence.correction_pass);
|
|
for (double value : failed.evidence.worst_force_terms)
|
|
EXPECT_TRUE(std::isfinite(value));
|
|
for (std::size_t i = 0;
|
|
i < state.Value().CommittedState().displacement.Size(); ++i)
|
|
EXPECT_DOUBLE_EQ(state.Value().CommittedState().low.displacement[i], 0);
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
LateBackendFailurePreservesActualCommittedMotionAndHistory) {
|
|
auto domain = fesa::Domain::Create(OneFree());
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
ASSERT_TRUE(state.HasValue());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
ControlledTangent tangent(ControlledTangent::kFail, 2);
|
|
auto status = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteInto(domain.Value(), state.Value());
|
|
EXPECT_FALSE(status.IsOk());
|
|
ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U);
|
|
EXPECT_EQ(state.Value().CommittedState().increment_index, 1U);
|
|
EXPECT_GT(state.Value().CommittedState().displacement[0], 0);
|
|
const auto committed = state.Value().CommittedState();
|
|
auto bad = committed;
|
|
bad.frame_index = 2;
|
|
bad.increment_index = 2;
|
|
bad.time = .02;
|
|
bad.solid_results[0].cauchy_stress[0] =
|
|
std::numeric_limits<double>::infinity();
|
|
EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk());
|
|
ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U);
|
|
EXPECT_DOUBLE_EQ(state.Value().CommittedState().displacement[0],
|
|
committed.displacement[0]);
|
|
EXPECT_DOUBLE_EQ(state.Value().CommittedState().strain_energy,
|
|
committed.strain_energy);
|
|
EXPECT_EQ(state.Value().CommittedState().solid_results[0].cauchy_stress,
|
|
committed.solid_results[0].cauchy_stress);
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
NonfiniteWeightedEvidenceCannotContaminateAcceptedHistory) {
|
|
auto domain = fesa::Domain::Create(OneFree());
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
ControlledTangent tangent(ControlledTangent::kFail, 2);
|
|
EXPECT_FALSE(fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteInto(domain.Value(), state.Value())
|
|
.IsOk());
|
|
ASSERT_EQ(state.Value().AcceptedFrames().size(), 2U);
|
|
const auto committed = state.Value().CommittedState();
|
|
for (double invalid : {std::numeric_limits<double>::quiet_NaN(),
|
|
std::numeric_limits<double>::infinity()}) {
|
|
for (unsigned field = 0; field < 8; ++field) {
|
|
auto bad = committed;
|
|
bad.frame_index = bad.increment_index = 2;
|
|
bad.time = .02;
|
|
auto& evidence = bad.iterations.back();
|
|
if (field == 0) evidence.weighted_residual_numerator = invalid;
|
|
if (field == 1) evidence.recombined_residual_numerator = invalid;
|
|
if (field == 2) evidence.old_residual_numerator = invalid;
|
|
if (field == 3) evidence.weighted_equivalence_metric = invalid;
|
|
if (field >= 4) evidence.worst_weighted_terms[field - 4] = invalid;
|
|
EXPECT_FALSE(state.Value().CommitAcceptedFrame(bad).IsOk());
|
|
EXPECT_EQ(state.Value().AcceptedFrames().size(), 2U);
|
|
EXPECT_EQ(state.Value().CommittedState().solid_results[0].cauchy_stress,
|
|
committed.solid_results[0].cauchy_stress);
|
|
EXPECT_DOUBLE_EQ(state.Value().CommittedState().displacement[0],
|
|
committed.displacement[0]);
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
AlphaZeroLinearUnforcedTailConservesEnergyAndTrapezoidWork) {
|
|
auto d = OneFree(false, .01, 100);
|
|
d.amplitudes = {{"A", "A", {}, {{0, 0}, {.1, 1}, {.2, 0}, {1, 0}}}};
|
|
d.dynamic_step->loads[0].amplitude_index = 0;
|
|
auto domain = fesa::Domain::Create(std::move(d));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
fesa::MklPardisoTangentSolver tangent;
|
|
auto status =
|
|
fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteIntoForVerification(domain.Value(), state.Value(), 0);
|
|
ASSERT_TRUE(status.IsOk()) << solid_hht_test::Message(status);
|
|
const auto& frames = state.Value().AcceptedFrames();
|
|
const double tail = frames[20].kinetic_energy + frames[20].strain_energy;
|
|
ASSERT_GT(tail, 1e-7);
|
|
double work = 0;
|
|
for (std::size_t n = 1; n < frames.size(); ++n) {
|
|
work += .5 *
|
|
(frames[n - 1].external_force[0] + frames[n].external_force[0]) *
|
|
(frames[n].displacement[0] - frames[n - 1].displacement[0]);
|
|
EXPECT_NEAR(frames[n].external_work, work, 1e-16);
|
|
EXPECT_NEAR(frames[n].energy_balance_error, 0, 2e-12 * tail);
|
|
if (n >= 20)
|
|
EXPECT_NEAR(frames[n].kinetic_energy + frames[n].strain_energy, tail,
|
|
2e-12 * tail);
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics, InitialMassBackendFailuresPreserveEmptyHistory) {
|
|
for (bool factorization : {false, true}) {
|
|
auto domain = fesa::Domain::Create(OneFree());
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
ASSERT_TRUE(state.HasValue());
|
|
fesa::SerialParallelFor parallel;
|
|
FailedMass mass(factorization);
|
|
fesa::MklPardisoTangentSolver tangent;
|
|
const auto status =
|
|
fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteInto(domain.Value(), state.Value());
|
|
EXPECT_FALSE(status.IsOk());
|
|
EXPECT_TRUE(state.Value().AcceptedFrames().empty());
|
|
EXPECT_FALSE(state.Value().BeginTrialFromCommitted().HasValue());
|
|
const auto& records = state.Value().IterationRecords();
|
|
ASSERT_EQ(records.size(), factorization ? 1U : 2U);
|
|
if (!factorization) {
|
|
EXPECT_EQ(records.front().stage, "initial-mass");
|
|
EXPECT_EQ(records.front().outcome, "residual-refinement");
|
|
EXPECT_EQ(records.front().evidence.iteration_index, 0U);
|
|
EXPECT_FALSE(records.front().evidence.residual_pass);
|
|
}
|
|
EXPECT_EQ(records.back().stage, "failed");
|
|
EXPECT_EQ(records.back().failure_code, "injected-backend");
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics, AllFreeNoForceIsExactlyAtRestForLinearAndTl) {
|
|
for (bool nonlinear : {false, true}) {
|
|
auto result = RunDynamics(solid_hht_test::Brick(.01, 3, nonlinear));
|
|
ASSERT_TRUE(result.HasValue())
|
|
<< solid_hht_test::Message(result.GetStatus());
|
|
for (const auto& frame : result.Value().state.AcceptedFrames()) {
|
|
EXPECT_EQ(frame.iteration_count, 0U);
|
|
EXPECT_DOUBLE_EQ(frame.kinetic_energy, 0);
|
|
EXPECT_DOUBLE_EQ(frame.strain_energy, 0);
|
|
for (std::size_t i = 0; i < frame.displacement.Size(); ++i) {
|
|
EXPECT_DOUBLE_EQ(frame.displacement[i], 0);
|
|
EXPECT_DOUBLE_EQ(frame.velocity[i], 0);
|
|
EXPECT_DOUBLE_EQ(frame.acceleration[i], 0);
|
|
EXPECT_DOUBLE_EQ(frame.algorithmic_acceleration[i], 0);
|
|
EXPECT_DOUBLE_EQ(frame.residual[i], 0);
|
|
EXPECT_DOUBLE_EQ(frame.reaction[i], 0);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void CheckOriginalSourceDeck(const std::string& name) {
|
|
auto root = std::filesystem::path{__FILE__}.parent_path();
|
|
for (unsigned i = 0; i < 3; ++i) root = root.parent_path();
|
|
{
|
|
auto syntax =
|
|
fesa::AbaqusInputReader{}.Read(root / "reference/hexa" / name);
|
|
ASSERT_TRUE(syntax.HasValue());
|
|
auto domain = fesa::AbaqusDomainMapper{}.Map(syntax.Value());
|
|
ASSERT_TRUE(domain.HasValue())
|
|
<< solid_hht_test::Message(domain.GetStatus());
|
|
fesa::TbbParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
fesa::MklPardisoTangentSolver tangent;
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
ASSERT_TRUE(state.HasValue());
|
|
auto status =
|
|
fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteInto(domain.Value(), state.Value());
|
|
if (!status.IsOk()) {
|
|
const auto& r = state.Value().IterationRecords().back();
|
|
std::cout << name << " accepted=" << state.Value().AcceptedFrames().size()
|
|
<< " increment=" << r.increment_index
|
|
<< " iteration=" << r.evidence.iteration_index
|
|
<< " metric=" << r.evidence.residual_metric
|
|
<< " update=" << r.evidence.update_metric
|
|
<< " worst=" << r.evidence.worst_full_dof << '\n';
|
|
std::cout.precision(17);
|
|
std::cout << "Worst terms r,Ma,fint,F,S: ";
|
|
for (const double term : r.evidence.worst_force_terms)
|
|
std::cout << term << ' ';
|
|
std::cout << '\n';
|
|
std::cout << "Original,recombined,old residual,equivalence scale: ";
|
|
for (double term : r.evidence.worst_weighted_terms)
|
|
std::cout << term << ' ';
|
|
std::cout << " equivalence_metric="
|
|
<< r.evidence.weighted_equivalence_metric << '\n';
|
|
for (const auto& coordinate : r.failure_coordinates) {
|
|
std::cout
|
|
<< "Coordinate element=" << coordinate.element_index
|
|
<< " row=" << coordinate.local_row
|
|
<< " u=" << coordinate.displacement
|
|
<< " previous_u=" << coordinate.previous_displacement
|
|
<< " fresh_du=" << coordinate.fresh_correction
|
|
<< " node0_u=" << coordinate.reference_displacement
|
|
<< " previous_node0_u="
|
|
<< coordinate.previous_reference_displacement
|
|
<< " relative_high=" << coordinate.relative_high
|
|
<< " previous_relative_high=" << coordinate.previous_relative_high
|
|
<< " u_changed="
|
|
<< (coordinate.displacement != coordinate.previous_displacement)
|
|
<< " relative_changed="
|
|
<< (coordinate.relative_high != coordinate.previous_relative_high)
|
|
<< '\n';
|
|
}
|
|
if (r.failure_force_sum) {
|
|
const auto& sums = *r.failure_force_sum;
|
|
std::cout << "Failure force sums "
|
|
"internal,compensated,abs_local,inertia,compensated,abs_"
|
|
"products: "
|
|
<< sums.assembled_internal << ' ' << sums.compensated_internal
|
|
<< ' ' << sums.absolute_local_internal_sum << ' '
|
|
<< sums.assembled_inertia << ' ' << sums.compensated_inertia
|
|
<< ' ' << sums.absolute_inertia_product_sum << '\n';
|
|
}
|
|
double shadow_delta = 0, local_scale = 0;
|
|
for (const auto& dot : r.failure_linear_dot) {
|
|
std::cout << "Local dot element=" << dot.element_index
|
|
<< " row=" << dot.local_row
|
|
<< " sum_abs_products=" << dot.absolute_product_sum
|
|
<< " kernel_force=" << dot.kernel_force
|
|
<< " rounded_relative_compensated_force="
|
|
<< dot.compensated_force
|
|
<< " delta=" << dot.compensated_force - dot.kernel_force
|
|
<< '\n';
|
|
shadow_delta += dot.compensated_force - dot.kernel_force;
|
|
local_scale += dot.absolute_product_sum;
|
|
}
|
|
std::cout << "Assembled rounded-relative shadow delta=" << shadow_delta
|
|
<< " local_product_scale=" << local_scale << '\n';
|
|
for (const auto& rec : state.Value().IterationRecords()) {
|
|
if (rec.newton_repeat &&
|
|
(rec.newton_repeat->repeat || rec.evidence.iteration_index < 6 ||
|
|
rec.evidence.iteration_index > 94)) {
|
|
const auto& probe = *rec.newton_repeat;
|
|
std::cout << "Source exact repeat " << rec.increment_index << ' '
|
|
<< rec.evidence.iteration_index
|
|
<< " history_entries=" << probe.history_entry_count
|
|
<< " matched_entry=" << probe.matched_entry_iteration
|
|
<< " u_changed=" << probe.changed_translations
|
|
<< " holding_mask_changed=" << probe.changed_holding_mask
|
|
<< " correction_changed="
|
|
<< probe.changed_correction_metrics << '\n';
|
|
}
|
|
if (rec.evidence.iteration_index < 6 ||
|
|
rec.evidence.iteration_index > 94)
|
|
std::cout << "Source trace " << rec.increment_index << ' '
|
|
<< rec.evidence.iteration_index << ' ' << rec.outcome
|
|
<< " ratio=" << rec.evidence.residual_metric
|
|
<< " correction=" << rec.evidence.correction_metrics[0]
|
|
<< '\n';
|
|
}
|
|
}
|
|
ASSERT_TRUE(status.IsOk())
|
|
<< name << ' ' << solid_hht_test::Message(status);
|
|
ASSERT_EQ(state.Value().AcceptedFrames().size(), 101U);
|
|
EXPECT_EQ(dofs.Value().FullDofCount(), 1326U);
|
|
const auto& final = state.Value().CommittedState();
|
|
EXPECT_DOUBLE_EQ(final.time, 10);
|
|
EXPECT_EQ(final.solid_results.size(), 2400U);
|
|
std::cout << name
|
|
<< " accepted=101 final_u_norm=" << final.displacement.Norm()
|
|
<< " final_force_metric="
|
|
<< final.iterations.back().residual_metric
|
|
<< " final_update_metric="
|
|
<< final.iterations.back().update_metric << '\n';
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtSourceSmoke, ReadOnlyOriginalLinearAcceptsAll101Endpoints) {
|
|
CheckOriginalSourceDeck("hexa_dynlinear.inp");
|
|
}
|
|
TEST(SolidHhtSourceSmoke, ReadOnlyOriginalNonlinearAcceptsAll101Endpoints) {
|
|
CheckOriginalSourceDeck("hexa_dynnonlinear.inp");
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
WeightedIdentityRetainsNonzeroOldResidualAndEffectiveTangent) {
|
|
const fesa::SparsePattern pattern{{0, 2, 4}, {0, 1, 0, 1}};
|
|
auto mass = fesa::SparseMatrix::FromCoo(
|
|
2, 2,
|
|
{{0, 0, 2, 0, 0}, {0, 1, .3, 0, 1}, {1, 0, .3, 0, 2}, {1, 1, 1.5, 0, 3}},
|
|
pattern);
|
|
ASSERT_TRUE(mass.HasValue());
|
|
constexpr double alpha = -.05, beta = .275625, h = .2;
|
|
fesa::SolidDynamicAcceptedFrameState old;
|
|
old.acceleration = fesa::Vector(2);
|
|
old.acceleration[0] = .4;
|
|
old.acceleration[1] = -.7;
|
|
old.internal_force = fesa::Vector(2);
|
|
old.internal_force[0] = 1;
|
|
old.internal_force[1] = -2;
|
|
old.external_force = fesa::Vector(2);
|
|
old.external_force[0] = .1;
|
|
old.external_force[1] = .2;
|
|
old.residual = mass.Value().Multiply(old.acceleration);
|
|
for (unsigned i = 0; i < 2; ++i)
|
|
old.residual[i] += old.internal_force[i] - old.external_force[i];
|
|
ASSERT_NE(old.residual[0], 0);
|
|
ASSERT_NE(old.residual[1], 0);
|
|
const std::array<double, 2> u{.12, -.03}, predictor{.09, -.01},
|
|
direction{.7, -.4};
|
|
for (double cubic : {0.0, 11.0}) {
|
|
const auto evaluate = [&](double epsilon) {
|
|
fesa::SolidDynamicAcceptedFrameState next;
|
|
next.acceleration = fesa::Vector(2);
|
|
next.algorithmic_acceleration = fesa::Vector(2);
|
|
next.internal_force = fesa::Vector(2);
|
|
next.external_force = fesa::Vector(2);
|
|
next.external_force[0] = .3;
|
|
next.external_force[1] = -.8;
|
|
const double x = u[0] + epsilon * direction[0];
|
|
const double y = u[1] + epsilon * direction[1];
|
|
for (unsigned i = 0; i < 2; ++i) {
|
|
const double q = i == 0 ? x : y;
|
|
next.algorithmic_acceleration[i] = (q - predictor[i]) / (beta * h * h);
|
|
next.acceleration[i] =
|
|
(next.algorithmic_acceleration[i] + alpha * old.acceleration[i]) /
|
|
(1 + alpha);
|
|
}
|
|
next.internal_force[0] = 7 * x + 2 * y + cubic * x * x * x;
|
|
next.internal_force[1] = 2 * x + 5 * y + cubic * y * y * y;
|
|
next.residual = mass.Value().Multiply(next.acceleration);
|
|
const auto MA = mass.Value().Multiply(next.algorithmic_acceleration);
|
|
std::array<fesa::detail::SolidHhtForceIdentity, 2> identities{};
|
|
for (unsigned i = 0; i < 2; ++i) {
|
|
next.residual[i] += next.internal_force[i] - next.external_force[i];
|
|
EXPECT_TRUE(fesa::detail::EvaluateSolidHhtForceIdentity(
|
|
mass.Value(), i, MA[i], next, old, alpha, &identities[i]));
|
|
EXPECT_LE(identities[i].metric, 1e-12);
|
|
EXPECT_NEAR(
|
|
identities[i].original, identities[i].solved,
|
|
512 * std::numeric_limits<double>::epsilon() * identities[i].scale);
|
|
EXPECT_GT(
|
|
std::abs(identities[i].solved - (1 + alpha) * next.residual[i]),
|
|
.01);
|
|
}
|
|
return identities;
|
|
};
|
|
constexpr double epsilon = 1e-6;
|
|
const auto left = evaluate(-epsilon), right = evaluate(epsilon);
|
|
const std::array<double, 2> tangent_direction{
|
|
(2 * direction[0] + .3 * direction[1]) / (beta * h * h) +
|
|
(1 + alpha) * ((7 + 3 * cubic * u[0] * u[0]) * direction[0] +
|
|
2 * direction[1]),
|
|
(.3 * direction[0] + 1.5 * direction[1]) / (beta * h * h) +
|
|
(1 + alpha) * (2 * direction[0] +
|
|
(5 + 3 * cubic * u[1] * u[1]) * direction[1])};
|
|
for (unsigned i = 0; i < 2; ++i)
|
|
EXPECT_NEAR((right[i].solved - left[i].solved) / (2 * epsilon),
|
|
tangent_direction[i], 1e-9 * std::abs(tangent_direction[i]));
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
WeightedIdentityRejectsNonfiniteAndOverflowBeforeNormReduction) {
|
|
const fesa::SparsePattern pattern{{0, 1}, {0}};
|
|
auto mass = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 2, 0, 0}}, pattern);
|
|
ASSERT_TRUE(mass.HasValue());
|
|
fesa::SolidDynamicAcceptedFrameState old, next;
|
|
for (auto* frame : {&old, &next})
|
|
for (auto* v :
|
|
{&frame->acceleration, &frame->algorithmic_acceleration,
|
|
&frame->internal_force, &frame->external_force, &frame->residual})
|
|
*v = fesa::Vector(1);
|
|
fesa::detail::SolidHhtForceIdentity identity;
|
|
EXPECT_TRUE(fesa::detail::EvaluateSolidHhtForceIdentity(
|
|
mass.Value(), 0, 0, next, old, -.05, &identity));
|
|
EXPECT_EQ(identity.scale, 0);
|
|
EXPECT_EQ(identity.metric, 0);
|
|
for (double invalid : {std::numeric_limits<double>::quiet_NaN(),
|
|
std::numeric_limits<double>::infinity()}) {
|
|
EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity(
|
|
mass.Value(), 0, invalid, next, old, -.05, &identity));
|
|
for (unsigned field = 0; field < 8; ++field) {
|
|
auto bad = next;
|
|
fesa::Vector* vectors[]{&bad.acceleration, &bad.algorithmic_acceleration,
|
|
&bad.internal_force, &bad.external_force,
|
|
&bad.residual};
|
|
if (field < 5) (*vectors[field])[0] = invalid;
|
|
auto previous = old;
|
|
if (field == 5) previous.acceleration[0] = invalid;
|
|
if (field == 6) previous.residual[0] = invalid;
|
|
if (field == 7) previous.external_force[0] = invalid;
|
|
EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity(
|
|
mass.Value(), 0, 0, bad, previous, -.05, &identity));
|
|
}
|
|
}
|
|
const double huge = std::numeric_limits<double>::max();
|
|
next.internal_force[0] = next.external_force[0] = huge;
|
|
EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity(
|
|
mass.Value(), 0, 0, next, old, -.05, &identity));
|
|
next.internal_force[0] = next.external_force[0] = 0;
|
|
next.algorithmic_acceleration[0] = huge;
|
|
EXPECT_FALSE(fesa::detail::EvaluateSolidHhtForceIdentity(
|
|
mass.Value(), 0, 0, next, old, -.05, &identity));
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
FullFreeNewtonRecoversAllComponentsAfterCoupledTrialPerturbation) {
|
|
class CoupledKick final : public fesa::TangentSolver {
|
|
public:
|
|
fesa::Status FactorizeFresh(const fesa::SparseMatrix& matrix) override {
|
|
rows_ = matrix.Rows();
|
|
full_size_all_calls = full_size_all_calls && rows_ == 24;
|
|
if (kicked) solved_after_kick = true;
|
|
return real_.FactorizeFresh(matrix);
|
|
}
|
|
fesa::Status SolveCorrection(const fesa::Vector& rhs,
|
|
fesa::Vector& du) const override {
|
|
auto status = real_.SolveCorrection(rhs, du);
|
|
if (status.IsOk() && !kicked) {
|
|
// Deliberately perturb one transverse correction; the next full
|
|
// physical evaluation must recover every coupled original component.
|
|
du[1] += 1e-6;
|
|
kicked = true;
|
|
}
|
|
return status;
|
|
}
|
|
mutable bool kicked{false};
|
|
bool full_size_all_calls{true};
|
|
bool solved_after_kick{false};
|
|
|
|
private:
|
|
std::size_t rows_{0};
|
|
fesa::MklPardisoTangentSolver real_;
|
|
} tangent;
|
|
auto d = solid_hht_test::Brick(.01, 1);
|
|
d.dynamic_step->loads = {{"ALL", 1, .25, {}}};
|
|
auto domain = fesa::Domain::Create(std::move(d));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
auto status = fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteInto(domain.Value(), state.Value());
|
|
EXPECT_TRUE(tangent.kicked);
|
|
EXPECT_TRUE(tangent.full_size_all_calls);
|
|
EXPECT_TRUE(tangent.solved_after_kick);
|
|
ASSERT_TRUE(status.IsOk()) << solid_hht_test::Message(status);
|
|
const double roundoff = 512 * std::numeric_limits<double>::epsilon();
|
|
for (const auto& frame : state.Value().AcceptedFrames()) {
|
|
EXPECT_LE(frame.iterations.back().residual_metric, 1e-10);
|
|
EXPECT_LE(frame.iterations.back().update_metric, 1e-12);
|
|
for (unsigned n = 0; n < 8; ++n)
|
|
EXPECT_NEAR(frame.displacement[3 * n], .5 * frame.time * frame.time,
|
|
roundoff * .5 * frame.time * frame.time);
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtFloat64History,
|
|
RepresentableNeighborUsesHeldCouplingAndOriginalFullForceMetric) {
|
|
const fesa::SparsePattern pattern{{0, 2, 4}, {0, 1, 0, 1}};
|
|
auto mass = fesa::SparseMatrix::FromCoo(2, 2,
|
|
{{0, 0, 1000, 0, 0},
|
|
{0, 1, 0, 0, 1},
|
|
{1, 0, 0, 0, 2},
|
|
{1, 1, 1000, 0, 3}},
|
|
pattern);
|
|
auto tangent = fesa::SparseMatrix::FromCoo(2, 2,
|
|
{{0, 0, 1e6, 0, 0},
|
|
{0, 1, 1e6, 0, 1},
|
|
{1, 0, 1e6, 0, 2},
|
|
{1, 1, 2e6, 0, 3}},
|
|
pattern);
|
|
ASSERT_TRUE(mass.HasValue());
|
|
ASSERT_TRUE(tangent.HasValue());
|
|
fesa::SolidDynamicAcceptedFrameState frame;
|
|
for (auto* v :
|
|
{&frame.displacement, &frame.inertia_force, &frame.internal_force,
|
|
&frame.external_force, &frame.residual, &frame.residual_scale})
|
|
*v = fesa::Vector(2);
|
|
frame.displacement[0] = 1.5;
|
|
frame.displacement[1] = .5;
|
|
frame.inertia_force[0] = -.5;
|
|
frame.internal_force[0] = .500000000102;
|
|
frame.inertia_force[1] = -50;
|
|
frame.internal_force[1] = 50;
|
|
for (unsigned i = 0; i < 2; ++i) {
|
|
frame.residual[i] = frame.inertia_force[i] + frame.internal_force[i];
|
|
frame.residual_scale[i] =
|
|
std::abs(frame.inertia_force[i]) + std::abs(frame.internal_force[i]);
|
|
}
|
|
const double initial = std::abs(frame.residual[0]) / frame.residual_scale[0];
|
|
ASSERT_GT(initial, 1e-10);
|
|
ASSERT_EQ(frame.residual[1], 0); // Equilibrated coupled coordinate is held.
|
|
fesa::detail::SolidHhtRepresentableNeighbor selected;
|
|
ASSERT_TRUE(fesa::detail::SelectSolidHhtRepresentableNeighbor(
|
|
mass.Value(), tangent.Value(), {0, 1}, frame, 1, -.05, .275625,
|
|
&selected));
|
|
ASSERT_TRUE(selected.found);
|
|
EXPECT_EQ(selected.full_dof, 1U);
|
|
EXPECT_EQ(selected.displacement,
|
|
std::nextafter(.5, -std::numeric_limits<double>::infinity()));
|
|
EXPECT_EQ(selected.delta, selected.displacement - frame.displacement[1]);
|
|
EXPECT_LT(selected.predicted_metric, initial);
|
|
EXPECT_LE(selected.predicted_metric, 1e-10);
|
|
// Independent separate inertia/internal-force update checks the original
|
|
// physical component score, including the previously equilibrated row.
|
|
double exact_metric = 0;
|
|
for (unsigned i = 0; i < 2; ++i) {
|
|
const double ma = frame.inertia_force[i] +
|
|
(i == 1 ? 1000 : 0) * selected.delta / (.95 * .275625);
|
|
const double f =
|
|
frame.internal_force[i] + (i == 0 ? 1e6 : 2e6) * selected.delta;
|
|
exact_metric = (std::max)(exact_metric,
|
|
std::abs(ma + f) / (std::abs(ma) + std::abs(f)));
|
|
}
|
|
EXPECT_EQ(selected.predicted_metric, exact_metric);
|
|
auto alternative_tangent = fesa::SparseMatrix::FromCoo(2, 2,
|
|
{{0, 0, 4e5, 0, 0},
|
|
{0, 1, 1e6, 0, 1},
|
|
{1, 0, 1e6, 0, 2},
|
|
{1, 1, 2e6, 0, 3}},
|
|
pattern);
|
|
ASSERT_TRUE(alternative_tangent.HasValue());
|
|
std::vector<fesa::detail::SolidHhtRepresentableNeighbor> ranked;
|
|
ASSERT_TRUE(fesa::detail::SelectSolidHhtRepresentableNeighbor(
|
|
mass.Value(), alternative_tangent.Value(), {0, 1}, frame, 1, -.05,
|
|
.275625, &selected, &ranked));
|
|
ASSERT_EQ(ranked.size(), 2U);
|
|
EXPECT_EQ(ranked[0].full_dof, 0U);
|
|
EXPECT_EQ(ranked[1].full_dof, 1U);
|
|
EXPECT_LT(ranked[0].predicted_metric, ranked[1].predicted_metric);
|
|
EXPECT_LT(ranked[1].predicted_metric, initial);
|
|
// Independent actual-score fixture models a rounded force evaluation whose
|
|
// best linear prediction worsens the physical ratio, while another improves.
|
|
// Neither prediction is accepted as an actual force result.
|
|
const std::array<double, 2> actual_scores{initial * 1.1, initial * .8};
|
|
std::vector<std::size_t> tried;
|
|
auto accepted =
|
|
fesa::detail::TrySolidHhtRankedNeighbors(ranked, [&](const auto& trial) {
|
|
tried.push_back(trial.full_dof);
|
|
return fesa::Result<bool>::Success(actual_scores[trial.full_dof] <
|
|
initial);
|
|
});
|
|
ASSERT_TRUE(accepted.HasValue());
|
|
EXPECT_TRUE(accepted.Value());
|
|
EXPECT_EQ(tried, (std::vector<std::size_t>{0, 1}));
|
|
tried.clear();
|
|
auto invalid =
|
|
fesa::detail::TrySolidHhtRankedNeighbors(ranked, [&](const auto& trial) {
|
|
tried.push_back(trial.full_dof);
|
|
return fesa::Result<bool>::Failure(
|
|
fesa::Status::Failure(fesa::FailureCategory::kSolver, {}));
|
|
});
|
|
EXPECT_FALSE(invalid.HasValue());
|
|
EXPECT_EQ(tried, (std::vector<std::size_t>{0}));
|
|
}
|
|
|
|
TEST(SolidHhtDynamics,
|
|
RepresentableNeighborKeepsExactZeroAndRejectsInvalidArithmetic) {
|
|
const fesa::SparsePattern pattern{{0, 1}, {0}};
|
|
auto mass = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 1, 0, 0}}, pattern);
|
|
auto tangent = fesa::SparseMatrix::FromCoo(1, 1, {{0, 0, 2, 0, 0}}, pattern);
|
|
ASSERT_TRUE(mass.HasValue());
|
|
ASSERT_TRUE(tangent.HasValue());
|
|
fesa::SolidDynamicAcceptedFrameState frame;
|
|
for (auto* v :
|
|
{&frame.displacement, &frame.inertia_force, &frame.internal_force,
|
|
&frame.external_force, &frame.residual, &frame.residual_scale})
|
|
*v = fesa::Vector(1);
|
|
fesa::detail::SolidHhtRepresentableNeighbor selected;
|
|
EXPECT_TRUE(fesa::detail::SelectSolidHhtRepresentableNeighbor(
|
|
mass.Value(), tangent.Value(), {0}, frame, .1, -.05, .275625, &selected));
|
|
EXPECT_FALSE(selected.found);
|
|
for (double invalid : {std::numeric_limits<double>::quiet_NaN(),
|
|
std::numeric_limits<double>::infinity()}) {
|
|
auto bad = frame;
|
|
bad.internal_force[0] = invalid;
|
|
EXPECT_FALSE(fesa::detail::SelectSolidHhtRepresentableNeighbor(
|
|
mass.Value(), tangent.Value(), {0}, bad, .1, -.05, .275625, &selected));
|
|
}
|
|
EXPECT_FALSE(fesa::detail::SelectSolidHhtRepresentableNeighbor(
|
|
mass.Value(), tangent.Value(), {0}, frame, 0, -.05, .275625, &selected));
|
|
}
|
|
|
|
TEST(SolidHhtFloat64History,
|
|
ActualNewtonStallAtLowHalfOfTwoCycleTriesFallback) {
|
|
// The current physical force ratio improved from the previous iterate, but
|
|
// the next complete Newton trial returns to the larger ratio. Historical
|
|
// progress must not suppress actual evaluation or representable fallback.
|
|
constexpr double previous_metric = 1.6e-10;
|
|
constexpr double current_metric = 1.4e-10;
|
|
constexpr double actual_neighbor_metric = 1.2e-10;
|
|
ASSERT_LT(current_metric, previous_metric);
|
|
fesa::Vector current(2), previous(2);
|
|
current[0] = .5;
|
|
previous[0] = 1;
|
|
const std::vector<bool> mask(2, false);
|
|
const std::array<double, 4> metrics{1e-20, 2e-20, 3e-20, 4e-20};
|
|
fesa::detail::SolidHhtNewtonHistory history;
|
|
ASSERT_TRUE(history.Capture(previous, mask, metrics, 0, false));
|
|
ASSERT_TRUE(history.Capture(current, mask, metrics, 1, false));
|
|
std::vector<std::string> evaluated;
|
|
auto selected = fesa::detail::TrySolidHhtActualStep(
|
|
false, {1e-20, 2e-20, 3e-20, 4e-20},
|
|
[&] {
|
|
evaluated.push_back("complete-newton");
|
|
fesa::SolidDynamicNewtonRepeatEvidence repeat;
|
|
EXPECT_TRUE(history.Match(previous, mask, metrics, false, &repeat));
|
|
return fesa::Result<bool>::Success(!repeat.repeat);
|
|
},
|
|
[&] {
|
|
evaluated.push_back("actual-neighbor");
|
|
return fesa::Result<bool>::Success(actual_neighbor_metric <
|
|
current_metric);
|
|
});
|
|
ASSERT_TRUE(selected.HasValue());
|
|
EXPECT_TRUE(selected.Value());
|
|
EXPECT_EQ(evaluated,
|
|
(std::vector<std::string>{"complete-newton", "actual-neighbor"}));
|
|
evaluated.clear();
|
|
auto productive = fesa::detail::TrySolidHhtActualStep(
|
|
false, {1e-20, 2e-20, 3e-20, 4e-20},
|
|
[&] {
|
|
evaluated.push_back("complete-newton");
|
|
return fesa::Result<bool>::Success(true);
|
|
},
|
|
[&] {
|
|
evaluated.push_back("actual-neighbor");
|
|
return fesa::Result<bool>::Success(true);
|
|
});
|
|
ASSERT_TRUE(productive.HasValue());
|
|
EXPECT_FALSE(productive.Value()); // Keep complete Newton common application.
|
|
EXPECT_EQ(evaluated, (std::vector<std::string>{"complete-newton"}));
|
|
}
|
|
|
|
TEST(SolidHhtFloat64History,
|
|
ExactNewtonLoopHistoryDetectsThreeStateRecurrence) {
|
|
// Three distinct entry states form a closed deterministic Newton cycle.
|
|
// The prospective state is neither the current nor its immediate predecessor.
|
|
fesa::detail::SolidHhtNewtonHistory history;
|
|
fesa::Vector first(3), second(3), third(3);
|
|
first[0] = 1;
|
|
second[0] = 2;
|
|
third[0] = 3;
|
|
const std::vector<bool> first_mask{true, false}, second_mask{false, true};
|
|
const std::array<double, 4> first_metrics{1e-20, 2e-20, 3e-20, 4e-20};
|
|
const std::array<double, 4> second_metrics{2e-20, 3e-20, 4e-20, 5e-20};
|
|
ASSERT_TRUE(history.Capture(first, first_mask, first_metrics, 4, false));
|
|
ASSERT_TRUE(history.Capture(second, second_mask, second_metrics, 5, false));
|
|
ASSERT_TRUE(history.Capture(third, first_mask, second_metrics, 6, false));
|
|
EXPECT_EQ(history.Size(), 3U);
|
|
fesa::SolidDynamicNewtonRepeatEvidence e;
|
|
ASSERT_TRUE(history.Match(first, first_mask, first_metrics, false, &e));
|
|
EXPECT_TRUE(e.repeat);
|
|
EXPECT_EQ(e.matched_entry_iteration, 4U);
|
|
EXPECT_EQ(e.history_entry_count, 3U);
|
|
EXPECT_EQ(e.changed_translations, 0U);
|
|
EXPECT_EQ(e.changed_holding_mask, 0U);
|
|
EXPECT_EQ(e.changed_correction_metrics, 0U);
|
|
ASSERT_TRUE(history.Match(first, second_mask, first_metrics, false, &e));
|
|
EXPECT_FALSE(e.repeat);
|
|
auto different_metrics = first_metrics;
|
|
different_metrics[3] = std::nextafter(first_metrics[3], 1.0);
|
|
ASSERT_TRUE(history.Match(first, first_mask, different_metrics, false, &e));
|
|
EXPECT_TRUE(
|
|
e.repeat); // Both finite/PASS controls produce the same transition.
|
|
EXPECT_EQ(e.changed_correction_metrics, 1U);
|
|
auto failing_metrics = first_metrics;
|
|
failing_metrics[2] = 1.1e-10;
|
|
ASSERT_TRUE(history.Match(first, first_mask, failing_metrics, false, &e));
|
|
EXPECT_FALSE(e.repeat);
|
|
first[2] = std::numeric_limits<double>::denorm_min();
|
|
ASSERT_TRUE(history.Match(first, first_mask, first_metrics, false, &e));
|
|
EXPECT_FALSE(e.repeat);
|
|
first[2] = 0;
|
|
ASSERT_TRUE(history.Match(first, first_mask, first_metrics, true, &e));
|
|
EXPECT_FALSE(e.repeat); // Actual physical equilibrium keeps ordinary Newton.
|
|
EXPECT_FALSE(
|
|
history.Match(fesa::Vector(2), first_mask, first_metrics, false, &e));
|
|
EXPECT_FALSE(history.Match(first, {false}, first_metrics, false, &e));
|
|
for (double bad : {std::numeric_limits<double>::quiet_NaN(),
|
|
std::numeric_limits<double>::infinity()}) {
|
|
first[2] = bad;
|
|
EXPECT_FALSE(history.Match(first, first_mask, first_metrics, false, &e));
|
|
EXPECT_FALSE(history.Capture(first, first_mask, first_metrics, 7, false));
|
|
first[2] = 0;
|
|
auto bad_metrics = first_metrics;
|
|
bad_metrics[0] = bad;
|
|
EXPECT_FALSE(history.Match(first, first_mask, bad_metrics, false, &e));
|
|
EXPECT_FALSE(history.Capture(first, first_mask, bad_metrics, 7, false));
|
|
}
|
|
}
|
|
|
|
TEST(SolidHhtFloat64History,
|
|
EquivalentFailedNewtonEntrySurvivesRepresentableMove) {
|
|
fesa::detail::SolidHhtNewtonHistory history;
|
|
fesa::Vector failed(3);
|
|
failed[0] = .5;
|
|
const std::vector<bool> mask{false, true};
|
|
const std::array<double, 4> correction{1e-20, 2e-20, 3e-20, 4e-20};
|
|
ASSERT_TRUE(history.Capture(failed, mask, correction, 23, true));
|
|
EXPECT_EQ(history.Size(), 0U);
|
|
auto large = correction;
|
|
large[0] = 1.1e-10;
|
|
ASSERT_TRUE(history.Capture(failed, mask, large, 23, false));
|
|
EXPECT_EQ(history.Size(), 0U);
|
|
ASSERT_TRUE(history.Capture(failed, mask, correction, 23, false));
|
|
fesa::Vector actual = failed;
|
|
std::vector<std::string> applied;
|
|
auto first = fesa::detail::TrySolidHhtActualStep(
|
|
false, correction,
|
|
[&] {
|
|
fesa::SolidDynamicNewtonRepeatEvidence e;
|
|
EXPECT_TRUE(history.Match(failed, mask, correction, false, &e));
|
|
return fesa::Result<bool>::Success(!e.repeat);
|
|
},
|
|
[&] {
|
|
actual[0] = std::nextafter(actual[0], 0.0);
|
|
applied.push_back("accepted-ulp");
|
|
return fesa::Result<bool>::Success(true);
|
|
});
|
|
ASSERT_TRUE(first.HasValue());
|
|
ASSERT_TRUE(first.Value());
|
|
EXPECT_NE(actual[0], failed[0]);
|
|
auto new_metrics = correction;
|
|
new_metrics[0] = std::nextafter(new_metrics[0], 1.0);
|
|
// The ordinary Newton proposal returns to the exact previously failed state.
|
|
auto second = fesa::detail::TrySolidHhtActualStep(
|
|
false, new_metrics,
|
|
[&] {
|
|
applied.push_back("returning-newton-evaluated");
|
|
fesa::SolidDynamicNewtonRepeatEvidence e;
|
|
EXPECT_TRUE(history.Match(failed, mask, new_metrics, false, &e));
|
|
EXPECT_TRUE(e.repeat);
|
|
EXPECT_EQ(e.matched_entry_iteration, 23U);
|
|
return fesa::Result<bool>::Success(!e.repeat);
|
|
},
|
|
[&] {
|
|
actual[0] = std::nextafter(actual[0], 0.0);
|
|
applied.push_back("another-accepted-ulp");
|
|
return fesa::Result<bool>::Success(true);
|
|
});
|
|
ASSERT_TRUE(second.HasValue());
|
|
EXPECT_TRUE(second.Value());
|
|
EXPECT_EQ(applied, (std::vector<std::string>{"accepted-ulp",
|
|
"returning-newton-evaluated",
|
|
"another-accepted-ulp"}));
|
|
auto fresh = failed;
|
|
fresh[2] = std::numeric_limits<double>::denorm_min();
|
|
fesa::SolidDynamicNewtonRepeatEvidence e;
|
|
ASSERT_TRUE(history.Match(fresh, mask, new_metrics, false, &e));
|
|
EXPECT_FALSE(e.repeat);
|
|
ASSERT_TRUE(history.Match(failed, {true, false}, new_metrics, false, &e));
|
|
EXPECT_FALSE(e.repeat);
|
|
ASSERT_TRUE(history.Match(failed, mask, new_metrics, true, &e));
|
|
EXPECT_FALSE(e.repeat);
|
|
}
|
|
|
|
TEST(SolidHhtDynamics, UnsupportedPairedProductsPreserveEmptyHistory) {
|
|
// Preserve the original extreme finite model. Its coefficient products
|
|
// lie outside the approved error-free paired multiplication range.
|
|
const double q = std::ldexp(1.0, -960);
|
|
const double delta = std::ldexp(1.0, -1010);
|
|
const double force1 = 4 * q, force2 = 8 * delta;
|
|
auto d = OneFree(false, 1, 2);
|
|
d.materials = {{"M", 1e-300, 0, {}}};
|
|
d.densities = {{0, 27, {}}};
|
|
d.amplitudes = {{"A", "A", {}, {{0, 0}, {1, force1}, {2, force2}}}};
|
|
d.dynamic_step->loads = {{"N0", 1, 1, {}, 0, "A", {}}};
|
|
auto domain = fesa::Domain::Create(std::move(d));
|
|
ASSERT_TRUE(domain.HasValue());
|
|
auto model = fesa::AnalysisModel::Create(domain.Value());
|
|
ASSERT_TRUE(model.HasValue());
|
|
auto dofs = fesa::DofManager::Create(model.Value());
|
|
ASSERT_TRUE(dofs.HasValue());
|
|
auto state =
|
|
fesa::SolidDynamicAnalysisState::Create(domain.Value(), dofs.Value());
|
|
ASSERT_TRUE(state.HasValue());
|
|
fesa::SerialParallelFor parallel;
|
|
fesa::MklPardisoSolver mass;
|
|
fesa::MklPardisoTangentSolver tangent;
|
|
const auto status =
|
|
fesa::SolidHhtDynamicProcedureExecutor(parallel, mass, tangent)
|
|
.ExecuteIntoForVerification(domain.Value(), state.Value(), 0);
|
|
ASSERT_FALSE(status.IsOk());
|
|
ASSERT_FALSE(status.Diagnostics().empty());
|
|
EXPECT_NE(status.Diagnostics().front().message.find("supported"),
|
|
std::string::npos);
|
|
EXPECT_TRUE(state.Value().AcceptedFrames().empty());
|
|
EXPECT_FALSE(state.Value().BeginTrialFromCommitted().HasValue());
|
|
}
|
|
|
|
TEST(SolidHhtDynamics, FailureShadowRecoversCancellationAndRejectsOverflow) {
|
|
double sum = 0, correction = 0, absolute_sum = 0;
|
|
for (double term : {1e16, 1.0, -1e16})
|
|
ASSERT_TRUE(fesa::detail::AccumulateSolidHhtForceShadow(
|
|
term, 0, &sum, &correction, &absolute_sum));
|
|
EXPECT_EQ(sum, 0);
|
|
EXPECT_EQ(sum + correction, 1); // Independent exactly representable sum.
|
|
sum = 0;
|
|
correction = 0;
|
|
absolute_sum = 0;
|
|
ASSERT_TRUE(fesa::detail::AccumulateSolidHhtForceShadow(
|
|
1e308, 0, &sum, &correction, &absolute_sum));
|
|
EXPECT_FALSE(fesa::detail::AccumulateSolidHhtForceShadow(
|
|
-1e308, 0, &sum, &correction, &absolute_sum));
|
|
EXPECT_EQ(sum, 1e308); // No invalid partial shadow is returned or committed.
|
|
EXPECT_EQ(absolute_sum, 1e308);
|
|
EXPECT_FALSE(fesa::detail::AccumulateSolidHhtForceShadow(
|
|
0, std::numeric_limits<double>::quiet_NaN(), &sum, &correction,
|
|
&absolute_sum));
|
|
}
|