tboudreaux/MeanField_SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
MeanField_SERiF/tests/solver/continuation_family.cpp
Emily Boudreaux cb20faeb3d test(tests): added test suite 
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.
2026-11-09 13:53:17 -05:00

224 lines
14 KiB
C++
Raw Permalink Blame History

#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});
}
Reference in New Issue View Git Blame Copy Permalink