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.
133 lines
7.7 KiB
C++
133 lines
7.7 KiB
C++
#include <limits>
|
|
#include <type_traits>
|
|
|
|
#include <catch2/catch_test_macros.hpp>
|
|
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
|
|
|
#include "serif/eos/models/polytropic.hpp"
|
|
#include "serif/tests/surface/models.hpp"
|
|
|
|
namespace {
|
|
using namespace serif;
|
|
using namespace serif::tests::surface;
|
|
}
|
|
|
|
TEST_CASE("Isobaric surface prescribes pressure without selecting a carrier", "[surface][unit]") {
|
|
// MeanField: Constant Pressure Surface Prescribes Only A Pressure Quantity.
|
|
const physics::IsobaricSurface condition{quantities::PressureValue{0.03125}};
|
|
CHECK(condition.target_pressure() == quantities::PressureValue{0.03125});
|
|
CHECK(condition.descriptor().target_pressure == 0.03125);
|
|
CHECK(condition.continuation_identity() == std::array<double, 1>{0.03125});
|
|
CHECK(physics::IsobaricSurface{physics::IsobaricSurface::Parameters{.PSurf = quantities::PressureValue{0.4}}}.target_pressure() == quantities::PressureValue{0.4});
|
|
|
|
for (const double invalid : {-0.1, std::numeric_limits<double>::infinity(), -std::numeric_limits<double>::infinity(), std::numeric_limits<double>::quiet_NaN()}) {
|
|
CHECK_THROWS_AS(physics::IsobaricSurface{quantities::PressureValue{invalid}}, std::invalid_argument);
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Polytropic surface resolves pressure through its declared enthalpy relation", "[surface][unit]") {
|
|
// MeanField: Polytropic EOS Resolves Constant Surface Pressure Through Its Enthalpy Relation.
|
|
using Formulation = physics::PressureSurfaceFormulation<eos::models::Polytrope, BarotropicState>;
|
|
const eos::models::Polytrope equation_of_state{3.0, 0.25};
|
|
const physics::IsobaricSurface condition{quantities::PressureValue{0.03125}};
|
|
const double required = eos::evaluate<quantities::SpecificEnthalpy>(equation_of_state, condition.target_pressure()).value();
|
|
|
|
physics::SurfacePoint<Formulation> point{};
|
|
point.fields.get<Enthalpy>() = required;
|
|
const auto result = physics::evaluate_pressure_surface<Formulation>(equation_of_state, condition, point);
|
|
REQUIRE(result);
|
|
CHECK(result->residual == 0.0);
|
|
CHECK(result->field_derivatives.get<Enthalpy>() * -0.19 == -0.19);
|
|
CHECK(result->position_derivative == std::array<double, 3>{});
|
|
}
|
|
|
|
TEST_CASE("General pressure surface uses the local composition", "[surface][unit]") {
|
|
// MeanField: General EOS Resolves Constant Surface Pressure With Local Composition.
|
|
const GeneralEOS equation_of_state;
|
|
const physics::IsobaricSurface condition{quantities::PressureValue{0.4}};
|
|
const auto first = general_point(2.0 * 0.4 + 3.0 * 0.2 + 5.0 * 0.1, 0.2, 0.1);
|
|
const auto second = general_point(2.0 * 0.4 + 3.0 * 0.3 + 5.0 * 0.1, 0.3, 0.1);
|
|
REQUIRE(first.fields.get<Enthalpy>() != second.fields.get<Enthalpy>());
|
|
|
|
const auto first_result = physics::evaluate_pressure_surface<GeneralFormulation>(equation_of_state, condition, first);
|
|
const auto second_result = physics::evaluate_pressure_surface<GeneralFormulation>(equation_of_state, condition, second);
|
|
REQUIRE(first_result);
|
|
REQUIRE(second_result);
|
|
CHECK(first_result->residual == 0.0);
|
|
CHECK(second_result->residual == 0.0);
|
|
}
|
|
|
|
TEST_CASE("General pressure surface Jacobian retains every MeanField state dependency", "[surface][unit]") {
|
|
// MeanField: preserve its finite-difference step, direction, and 2e-9 absolute gate.
|
|
const GeneralEOS equation_of_state;
|
|
const physics::IsobaricSurface condition{quantities::PressureValue{0.4}};
|
|
constexpr double step = 1.0e-7;
|
|
const auto point = general_point(1.7, 0.2, 0.1);
|
|
const auto forward = general_point(1.7 + step * 0.7, 0.2 - step * 0.2, 0.1 + step * 0.05);
|
|
const auto backward = general_point(1.7 - step * 0.7, 0.2 + step * 0.2, 0.1 - step * 0.05);
|
|
|
|
const auto result = physics::evaluate_pressure_surface<GeneralFormulation>(equation_of_state, condition, point);
|
|
const auto plus = physics::evaluate_pressure_surface<GeneralFormulation>(equation_of_state, condition, forward);
|
|
const auto minus = physics::evaluate_pressure_surface<GeneralFormulation>(equation_of_state, condition, backward);
|
|
REQUIRE(result);
|
|
REQUIRE(plus);
|
|
REQUIRE(minus);
|
|
const double action = result->field_derivatives.get<Enthalpy>() * 0.7 + result->field_derivatives.get<Entropy>() * -0.2 + result->field_derivatives.get<Composition>() * 0.05;
|
|
CHECK(action == 0.7 - 3.0 * -0.2 - 5.0 * 0.05);
|
|
CHECK_THAT((plus->residual - minus->residual) / (2.0 * step), Catch::Matchers::WithinAbs(action, 2.0e-9));
|
|
}
|
|
|
|
TEST_CASE("Isobaric surface supports density without an enthalpy quantity or EOS relation", "[surface][unit]") {
|
|
const GasEOS equation_of_state;
|
|
for (const double pressure : {0.0, 0.4, 2.0}) {
|
|
for (const double temperature : {0.5, 1.0, 3.0}) {
|
|
CAPTURE(pressure, temperature);
|
|
const physics::IsobaricSurface condition{quantities::PressureValue{pressure}};
|
|
const auto point = gas_point(pressure / temperature, temperature);
|
|
const auto inverse = physics::evaluate_pressure_surface<GasFormulation>(equation_of_state, condition, point);
|
|
const auto direct = physics::evaluate_pressure_surface<DirectGasFormulation>(equation_of_state, condition, point);
|
|
REQUIRE(inverse);
|
|
REQUIRE(direct);
|
|
CHECK(inverse->residual == 0.0);
|
|
CHECK_THAT(direct->residual, Catch::Matchers::WithinAbs(0.0, 1.0e-15));
|
|
CHECK(inverse->field_derivatives.get<Density>() == 1.0);
|
|
CHECK(inverse->field_derivatives.get<Temperature>() == pressure / (temperature * temperature));
|
|
CHECK(direct->field_derivatives.get<Density>() == temperature);
|
|
CHECK(direct->field_derivatives.get<Temperature>() == pressure / temperature);
|
|
}
|
|
}
|
|
}
|
|
|
|
TEST_CASE("Nonlinear density and direct pressure residuals have independent finite-difference checks", "[surface][unit]") {
|
|
const GasEOS equation_of_state;
|
|
const physics::IsobaricSurface condition{quantities::PressureValue{0.8}};
|
|
constexpr double step = 1.0e-6;
|
|
const auto check = [&]<typename Formulation>() {
|
|
const auto result = physics::evaluate_pressure_surface<Formulation>(equation_of_state, condition, gas_point(0.6, 1.3));
|
|
const auto plus = physics::evaluate_pressure_surface<Formulation>(equation_of_state, condition, gas_point(0.6 + step * 0.31, 1.3 - step * 0.27));
|
|
const auto minus = physics::evaluate_pressure_surface<Formulation>(equation_of_state, condition, gas_point(0.6 - step * 0.31, 1.3 + step * 0.27));
|
|
REQUIRE(result);
|
|
REQUIRE(plus);
|
|
REQUIRE(minus);
|
|
const double action = result->field_derivatives.template get<Density>() * 0.31 - result->field_derivatives.template get<Temperature>() * 0.27;
|
|
CHECK_THAT((plus->residual - minus->residual) / (2.0 * step), Catch::Matchers::WithinAbs(action, 2.0e-9));
|
|
};
|
|
check.template operator()<GasFormulation>();
|
|
check.template operator()<DirectGasFormulation>();
|
|
}
|
|
|
|
TEST_CASE("Shared thermodynamic targets accumulate their Jacobian contributions", "[surface][unit]") {
|
|
using SharedState = dependencies::Bindings<
|
|
dependencies::Binding<quantities::SpecificEnthalpy, Enthalpy>,
|
|
dependencies::Binding<quantities::Entropy, Enthalpy>,
|
|
dependencies::Binding<ElectronFraction, Enthalpy>
|
|
>;
|
|
using Formulation = physics::PressureSurfaceFormulation<GeneralEOS, SharedState>;
|
|
physics::SurfacePoint<Formulation> point{};
|
|
point.fields.get<Enthalpy>() = 0.2;
|
|
const auto result = physics::evaluate_pressure_surface<Formulation>(GeneralEOS{}, physics::IsobaricSurface{quantities::PressureValue{0.4}}, point);
|
|
REQUIRE(result);
|
|
CHECK_THAT(result->residual, Catch::Matchers::WithinAbs(0.2 - (2.0 * 0.4 + 8.0 * 0.2), 1.0e-15));
|
|
CHECK(result->field_derivatives.get<Enthalpy>() == 1.0 - 3.0 - 5.0);
|
|
}
|