Emily Boudreaux
cb20faeb3d
The test suite is extensive. Note that LLMs were used to aid in the development of the test suite. All LLM generated code was manually reviewed.
224 lines
14 KiB
C++
224 lines
14 KiB
C++
#include <array>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <memory>
|
|
|
|
#include <catch2/catch_approx.hpp>
|
|
#include <catch2/catch_test_macros.hpp>
|
|
|
|
#include "serif/solver/continuation/polytropic.hpp"
|
|
#include "serif/solver/continuation/identified_equations.hpp"
|
|
|
|
namespace {
|
|
namespace family = serif::solver::continuation;
|
|
namespace model = serif::model;
|
|
namespace dimensions = serif::dimensions;
|
|
using Polytrope = serif::eos::models::Polytrope;
|
|
|
|
struct Definition {
|
|
double pressure{0.01};
|
|
/** @brief State explicitly which captured physical input belongs to family identity. */
|
|
[[nodiscard]] std::array<double, 1> continuation_identity() const { return {pressure}; }
|
|
};
|
|
|
|
/** @brief Retain all physical prescriptions while changing selected independent inputs. */
|
|
[[nodiscard]] auto target_model(double index = 3.0, double constant = 0.36029966935264141, double momentum = 0.035,
|
|
double mass = 1.0, double pressure = 0.01,
|
|
std::array<double, 3> axis = {1.0, 2.0, 3.0}, std::array<double, 3> center = {0.2, -0.3, 0.4}) {
|
|
return model::make_stellar_model(Polytrope{index, constant}, Definition{pressure},
|
|
model::FixedTotalMass{dimensions::MassValue{mass}},
|
|
model::FixedAngularMomentum{{.angular_momentum = dimensions::AngularMomentumValue{momentum}, .axis = axis, .center = center}});
|
|
}
|
|
|
|
/** @brief Evaluate pressure independently through the EOS relation. */
|
|
[[nodiscard]] double pressure(const Polytrope& eos, double density) {
|
|
return eos.evaluate(serif::eos::relations::PressureFromDensity{}, dimensions::DensityValue{density}).value();
|
|
}
|
|
|
|
/** @brief Compare values using a relative tolerance that does not privilege one system of units. */
|
|
[[nodiscard]] bool close(double actual, double expected) {
|
|
return std::abs(actual - expected) <= 128.0 * std::numeric_limits<double>::epsilon() * std::abs(expected);
|
|
}
|
|
|
|
class TestEOS;
|
|
struct TestPressure {
|
|
using Specification = TestEOS;
|
|
using Value = dimensions::PressureValue;
|
|
static constexpr std::string_view name = "test pressure";
|
|
static constexpr std::string_view stable_id = "test.pressure";
|
|
/** @brief Read the extension's selected physical pressure. */
|
|
[[nodiscard]] static Value read(const TestEOS& value);
|
|
/** @brief Test extension domain validation without exceptions during candidate preparation. */
|
|
[[nodiscard]] static bool valid(Value value) { return std::isfinite(value.value()) && value.value() > 0.0; }
|
|
};
|
|
|
|
class TestEOS {
|
|
public:
|
|
using ContinuationParameters = family::ParameterList<TestPressure>;
|
|
double pressure{3.0};
|
|
bool ignore_updates{false};
|
|
double clone_offset{};
|
|
/** @brief Include all physical values; fault-injection switches are test behavior. */
|
|
[[nodiscard]] std::array<double, 1> continuation_identity() const { return {pressure}; }
|
|
/** @brief Deliberately allow a faulty extension to exercise installed-value readback. */
|
|
template <typename Values>
|
|
[[nodiscard]] TestEOS with_continuation_parameters(const Values& values) const {
|
|
TestEOS result = *this;
|
|
if (!ignore_updates) {
|
|
result.pressure = values.get(TestPressure{}).value();
|
|
}
|
|
return result;
|
|
}
|
|
/** @brief Deliberately allow a faulty clone to exercise exact target identity. */
|
|
[[nodiscard]] TestEOS clone_for_continuation() const {
|
|
TestEOS result = *this;
|
|
result.pressure += clone_offset;
|
|
return result;
|
|
}
|
|
};
|
|
|
|
TestPressure::Value TestPressure::read(const TestEOS& value) { return Value{value.pressure}; }
|
|
|
|
struct MoveOnlyDefinition {
|
|
std::unique_ptr<double> pressure;
|
|
/** @brief Owning captures still identify physical family inputs. */
|
|
[[nodiscard]] std::array<double, 1> continuation_identity() const { return {*pressure}; }
|
|
};
|
|
}
|
|
|
|
TEST_CASE("Physical continuation preserves pressure pivots and exact target prescriptions", "[unit][solver][continuation][family]") {
|
|
const auto target = target_model();
|
|
const auto bound = family::in(family::PolytropicIndex{}, family::AngularMomentum{})
|
|
.with(family::preserve_pressure_at_reference_density(dimensions::DensityValue{0.2})).bind(target);
|
|
CHECK(bound.target_coordinates() == std::array<double, 2>{3.0, 0.035});
|
|
CHECK(bound.starting_coordinates() == std::array<double, 2>{1.0, 0.0});
|
|
CHECK(bound.parameter_descriptors()[0].stable_id == "polytrope.index");
|
|
CHECK(bound.parameter_descriptors()[1].stable_id == "fixed_angular_momentum.target");
|
|
for (double index : {1.0, 1.5, 2.0, 2.999999, 3.0, 3.5}) {
|
|
const auto stage = bound.model_at({index, 0.02});
|
|
CHECK(close(pressure(stage.equation_of_state(), 0.2), pressure(target.equation_of_state(), 0.2)));
|
|
CHECK(stage.template specification<model::FixedTotalMass>().target().value() == 1.0);
|
|
CHECK(stage.equations().pressure == 0.01);
|
|
const auto& rotation = stage.template specification<model::FixedAngularMomentum>();
|
|
CHECK(rotation.target().value() == 0.02);
|
|
CHECK(rotation.axis() == target.template specification<model::FixedAngularMomentum>().axis());
|
|
CHECK(rotation.center() == target.template specification<model::FixedAngularMomentum>().center());
|
|
REQUIRE(bound.try_validate_member(stage));
|
|
}
|
|
const auto endpoint = bound.model_at(bound.target_coordinates());
|
|
CHECK(endpoint.equation_of_state().polytropic_constant() == target.equation_of_state().polytropic_constant());
|
|
CHECK(bound.model_signature(endpoint) == bound.model_signature(target));
|
|
CHECK(bound.coordinates_of(endpoint) == bound.target_coordinates());
|
|
}
|
|
|
|
TEST_CASE("Pressure pivot families are invariant under consistent changes of physical units", "[unit][solver][continuation][family]") {
|
|
constexpr double density = 0.4, density_units = 100.0, pressure_units = 0.01, index = 1.5, constant = 0.73;
|
|
const double converted = pressure_units * constant / std::pow(density_units, 1.0 + 1.0 / index);
|
|
const auto first = family::in(family::PolytropicIndex{}).with(family::preserve_pressure_at_reference_density(dimensions::DensityValue{density})).bind(target_model(index, constant));
|
|
const auto second = family::in(family::PolytropicIndex{}).with(family::preserve_pressure_at_reference_density(dimensions::DensityValue{density * density_units})).bind(target_model(index, converted));
|
|
for (double candidate : {1.0, 2.0, 3.0}) {
|
|
const auto original = first.model_at({candidate}), rescaled = second.model_at({candidate});
|
|
for (double sample : {0.01, 0.4, 3.0}) {
|
|
CHECK(close(pressure(rescaled.equation_of_state(), sample * density_units), pressure_units * pressure(original.equation_of_state(), sample)));
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Continuation family membership checks all unchanged prescriptions and captured equations", "[unit][solver][continuation][family]") {
|
|
const auto bound = family::in(family::AngularMomentum{}).bind(target_model());
|
|
CHECK(bound.try_validate_member(target_model(3.0, 0.36029966935264141, 0.01)));
|
|
CHECK_FALSE(bound.try_validate_member(target_model(2.0)));
|
|
CHECK_FALSE(bound.try_validate_member(target_model(3.0, 0.3)));
|
|
CHECK_FALSE(bound.try_validate_member(target_model(3.0, 0.36029966935264141, 0.01, 2.0)));
|
|
CHECK_FALSE(bound.try_validate_member(target_model(3.0, 0.36029966935264141, 0.01, 1.0, 0.02)));
|
|
CHECK_FALSE(bound.try_validate_member(target_model(3.0, 0.36029966935264141, 0.01, 1.0, 0.01, {0.0, 0.0, 1.0})));
|
|
CHECK_FALSE(bound.try_validate_member(target_model(3.0, 0.36029966935264141, 0.01, 1.0, 0.01, {1.0, 2.0, 3.0}, {})));
|
|
const auto index_family = family::in(family::PolytropicIndex{})
|
|
.with(family::preserve_pressure_at_reference_density(dimensions::DensityValue{0.2})).bind(target_model());
|
|
CHECK_FALSE(index_family.try_validate_member(target_model(1.5, 0.36029966935264141, 0.035)));
|
|
}
|
|
|
|
TEST_CASE("Continuation laws report invalid numerical candidates and reject inconsistent configuration", "[unit][solver][continuation][family]") {
|
|
const double nan = std::numeric_limits<double>::quiet_NaN();
|
|
for (double invalid : {0.0, -1.0, nan, std::numeric_limits<double>::infinity()}) {
|
|
CHECK_THROWS_AS(family::preserve_pressure_at_reference_density(dimensions::DensityValue{invalid}), std::invalid_argument);
|
|
}
|
|
const auto mismatched = family::in(family::PolytropicIndex{}).with(family::law(family::PolytropicConstant{}, family::depends_on(family::PolytropicIndex{}), [](double) { return 2.0; }));
|
|
CHECK_THROWS_AS(mismatched.bind(target_model()), std::invalid_argument);
|
|
const auto broken = family::in(family::PolytropicIndex{}).with(family::law(family::PolytropicConstant{}, family::depends_on(family::PolytropicIndex{}),
|
|
[nan](double index) { return index == 3.0 ? 0.36029966935264141 : nan; })).bind(target_model());
|
|
CHECK_FALSE(broken.try_model_at({2.0}));
|
|
CHECK_FALSE(broken.try_model_at({0.5}));
|
|
CHECK_FALSE(broken.try_model_at({nan}));
|
|
const auto physical = family::in(family::PolytropicIndex{}, family::AngularMomentum{})
|
|
.with(family::preserve_pressure_at_reference_density(dimensions::DensityValue{0.2})).bind(target_model());
|
|
CHECK_FALSE(physical.try_model_at({0.5, 0.01}));
|
|
CHECK_FALSE(physical.try_model_at({2.0, -0.01}));
|
|
CHECK_FALSE(physical.try_model_at({std::numeric_limits<double>::infinity(), 0.01}));
|
|
CHECK_FALSE(physical.try_model_at({2.0, nan}));
|
|
const auto overflowing = family::in(family::PolytropicIndex{}).with(
|
|
family::preserve_pressure_at_reference_density(dimensions::DensityValue{std::numeric_limits<double>::max()}));
|
|
CHECK_THROWS_AS(overflowing.bind(target_model()), std::invalid_argument);
|
|
using Function = double (*)(double);
|
|
CHECK_THROWS_AS(family::law(family::PolytropicConstant{}, family::depends_on(family::PolytropicIndex{}), static_cast<Function>(nullptr)), std::invalid_argument);
|
|
const family::ParameterLaw<family::PolytropicConstant, family::ParameterList<family::PolytropicIndex>, Function> bypassed{nullptr};
|
|
CHECK_THROWS_AS(family::in(family::PolytropicIndex{}).with(bypassed).bind(target_model()), std::invalid_argument);
|
|
using Member = double (dimensions::AngularMomentumValue::*)() const noexcept;
|
|
CHECK_THROWS_AS(family::law(family::PolytropicConstant{}, family::depends_on(family::AngularMomentum{}), static_cast<Member>(nullptr)), std::invalid_argument);
|
|
}
|
|
|
|
TEST_CASE("Continuation owns move-only laws and explicitly reproducible equation definitions", "[unit][solver][continuation][family]") {
|
|
auto definition = family::in(family::PolytropicIndex{}).with(family::law(family::PolytropicConstant{}, family::depends_on(family::PolytropicIndex{}),
|
|
[constant = std::make_unique<double>(0.3)](double index) { return *constant * index / 3.0; }));
|
|
const auto bound = std::move(definition).bind(target_model(3.0, 0.3));
|
|
CHECK(bound.model_at({1.5}).equation_of_state().polytropic_constant() == Catch::Approx(0.15));
|
|
CHECK(bound.model_at({3.0}).equation_of_state().polytropic_constant() == 0.3);
|
|
auto target = model::make_stellar_model(Polytrope{3.0, 0.3}, MoveOnlyDefinition{std::make_unique<double>(0.01)},
|
|
model::FixedAngularMomentum{dimensions::AngularMomentumValue{0.1}});
|
|
const auto movable = family::in(family::AngularMomentum{}).bind(std::move(target), [](const MoveOnlyDefinition& original) {
|
|
return MoveOnlyDefinition{std::make_unique<double>(*original.pressure)};
|
|
});
|
|
auto first = movable.model_at({0.02}), second = movable.model_at({0.04});
|
|
CHECK(first.equations().pressure.get() != second.equations().pressure.get());
|
|
CHECK(movable.try_validate_member(first));
|
|
}
|
|
|
|
TEST_CASE("Custom physical families verify installed values and exact target cloning", "[unit][solver][continuation][family]") {
|
|
const auto target = model::make_stellar_model(TestEOS{}, Definition{});
|
|
const auto custom = family::in(TestPressure{}).bind(target);
|
|
CHECK(custom.starting_coordinates() == std::array<double, 1>{3.0});
|
|
CHECK(custom.target_coordinates() == std::array<double, 1>{3.0});
|
|
const auto stage = custom.model_at({2.0});
|
|
CHECK(stage.equation_of_state().pressure == 2.0);
|
|
CHECK(custom.try_validate_member(stage));
|
|
const auto faulty = model::make_stellar_model(TestEOS{.ignore_updates = true}, Definition{},
|
|
model::FixedAngularMomentum{dimensions::AngularMomentumValue{0.1}});
|
|
const auto ignores = family::in(TestPressure{}).bind(faulty);
|
|
CHECK(ignores.try_model_at(ignores.target_coordinates()));
|
|
CHECK_FALSE(ignores.try_model_at({2.0}));
|
|
CHECK_FALSE(ignores.try_model_at({std::nextafter(3.0, 2.0)}));
|
|
CHECK(ignores.target_model().equation_of_state().pressure == 3.0);
|
|
const auto derived = family::in(family::AngularMomentum{}).with(family::law(TestPressure{}, family::depends_on(family::AngularMomentum{}),
|
|
[](dimensions::AngularMomentumValue momentum) { return dimensions::PressureValue{30.0 * momentum.value()}; })).bind(faulty);
|
|
CHECK(derived.try_model_at(derived.target_coordinates()));
|
|
CHECK_FALSE(derived.try_model_at({0.05}));
|
|
CHECK(derived.target_model().equation_of_state().pressure == 3.0);
|
|
const auto bad_clone = model::make_stellar_model(TestEOS{.clone_offset = std::numeric_limits<double>::epsilon() * 3.0}, Definition{},
|
|
model::FixedAngularMomentum{dimensions::AngularMomentumValue{0.1}});
|
|
CHECK_THROWS_AS(family::in(family::AngularMomentum{}).bind(bad_clone), std::invalid_argument);
|
|
CHECK(bad_clone.equation_of_state().pressure == 3.0);
|
|
}
|
|
|
|
TEST_CASE("Identified equation definitions retain typed captures and support move-only callables", "[unit][solver][continuation][family]") {
|
|
auto equations = family::identified_equations([offset = std::make_unique<double>(2.0)](const auto& model) {
|
|
return model.equation_of_state().polytropic_index() + *offset;
|
|
}, dimensions::PressureValue{0.02}, dimensions::QuantityValue<dimensions::GravitationalConstant>{1.0});
|
|
STATIC_CHECK_FALSE(std::copy_constructible<decltype(equations)>);
|
|
STATIC_CHECK(std::same_as<std::tuple_element_t<0, decltype(equations.identity)>, dimensions::PressureValue>);
|
|
CHECK(equations.continuation_identity() == std::array{0.02, 1.0});
|
|
CHECK(equations(target_model()) == 5.0);
|
|
auto moved = std::move(equations);
|
|
CHECK(moved(target_model(1.5)) == 3.5);
|
|
CHECK(moved.continuation_identity() == std::array{0.02, 1.0});
|
|
}
|