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.
164 lines
8.0 KiB
C++
164 lines
8.0 KiB
C++
#include <catch2/catch_test_macros.hpp>
|
|
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
|
|
|
#include "serif/tests/equations/test_equations.hpp"
|
|
#include "serif/eos/models/polytropic.hpp"
|
|
#include "serif/geometry/surface/automatic.hpp"
|
|
#include "serif/backend/mfem/geometry/surface/pressure_boundary.hpp"
|
|
#include "serif/tests/surface/boundary_fixture.hpp"
|
|
#include "serif/tests/surface/models.hpp"
|
|
#include "serif/utils/dependencies/computation_cache.hpp"
|
|
|
|
namespace
|
|
{
|
|
using namespace serif;
|
|
using namespace serif::tests::surface;
|
|
namespace backend = physics::backend::mfem;
|
|
|
|
using ContinuousDensity = ScalarTerm<quantities::MassDensity>;
|
|
using Pressure = ScalarTerm<quantities::Pressure>;
|
|
|
|
struct ForwardOnly : GasEOS
|
|
{
|
|
using Relations = eos::relations::RelationCatalog<PressureFromDensityAndTemperature>;
|
|
};
|
|
|
|
struct SymbolicEOS
|
|
{
|
|
};
|
|
}
|
|
|
|
TEST_CASE("Automatic isobaric compilation preserves the polytropic zero-pressure unit derivative", "[surface][automatic]")
|
|
{
|
|
using Fields = dependencies::KeyList<Enthalpy, Density, blocks::barotropic_constant::mass_normalization, blocks::rotation::angular_velocity>;
|
|
using Selection = physics::AutomaticPressureSurface<eos::models::Polytrope, Fields, dependencies::KeyList<Enthalpy>>;
|
|
using Formulation = Selection::Formulation;
|
|
const eos::models::Polytrope equation_of_state{3.0, 0.25};
|
|
physics::SurfacePoint<Formulation> point{};
|
|
point.fields.get<Enthalpy>() = 0.0;
|
|
|
|
const auto result = physics::evaluate_pressure_surface<Formulation>(equation_of_state, physics::Isobaric{quantities::PressureValue{0.0}}, point);
|
|
REQUIRE(result);
|
|
CHECK(result->residual == 0.0);
|
|
CHECK(result->field_derivatives.get<Enthalpy>() == 1.0);
|
|
CHECK(Formulation::resolves_carrier);
|
|
}
|
|
|
|
TEST_CASE("Automatic isobaric compilation uses inverse or direct EOS relations without enthalpy", "[surface][automatic]")
|
|
{
|
|
using Fields = dependencies::KeyList<ContinuousDensity, Temperature>;
|
|
using Rows = dependencies::KeyList<ContinuousDensity>;
|
|
const auto check = []<typename EOS>(const EOS& equation_of_state)
|
|
{
|
|
using Formulation = typename physics::AutomaticPressureSurface<EOS, Fields, Rows>::Formulation;
|
|
physics::SurfacePoint<Formulation> point{};
|
|
point.fields.template get<ContinuousDensity>() = 0.4;
|
|
point.fields.template get<Temperature>() = 2.0;
|
|
const auto result = physics::evaluate_pressure_surface<Formulation>(equation_of_state, physics::Isobaric{quantities::PressureValue{0.8}}, point);
|
|
REQUIRE(result);
|
|
CHECK(result->residual == 0.0);
|
|
|
|
constexpr double step = 1.0e-6;
|
|
auto plus = point;
|
|
auto minus = point;
|
|
plus.fields.template get<ContinuousDensity>() += step * 0.3;
|
|
minus.fields.template get<ContinuousDensity>() -= step * 0.3;
|
|
plus.fields.template get<Temperature>() -= step * 0.2;
|
|
minus.fields.template get<Temperature>() += step * 0.2;
|
|
const auto forward = physics::evaluate_pressure_surface<Formulation>(equation_of_state, physics::Isobaric{quantities::PressureValue{0.8}}, plus);
|
|
const auto backward = physics::evaluate_pressure_surface<Formulation>(equation_of_state, physics::Isobaric{quantities::PressureValue{0.8}}, minus);
|
|
REQUIRE(forward);
|
|
REQUIRE(backward);
|
|
const double derivative = result->field_derivatives.template get<ContinuousDensity>() * 0.3 - result->field_derivatives.template get<Temperature>() * 0.2;
|
|
CHECK_THAT((forward->residual - backward->residual) / (2.0 * step), Catch::Matchers::WithinAbs(derivative, 2.0e-9));
|
|
};
|
|
check(GasEOS{});
|
|
check(ForwardOnly{});
|
|
}
|
|
|
|
TEST_CASE("A pressure unknown needs no EOS relations or EOS graph input", "[surface][automatic]")
|
|
{
|
|
using Fields = dependencies::KeyList<Pressure>;
|
|
using Formulation = physics::AutomaticPressureSurface<SymbolicEOS, Fields, Fields>::Formulation;
|
|
physics::SurfacePoint<Formulation> point{};
|
|
point.fields.get<Pressure>() = 0.2;
|
|
const auto single_point = physics::evaluate_pressure_surface<Formulation>(SymbolicEOS{}, physics::Isobaric{quantities::PressureValue{0.2}}, point);
|
|
REQUIRE(single_point);
|
|
CHECK(single_point->residual == 0.0);
|
|
CHECK(single_point->field_derivatives.get<Pressure>() == 1.0);
|
|
|
|
const auto inputs = [](const double target)
|
|
{
|
|
return dependencies::make_dependency_context(physics::pressure_condition(physics::Isobaric{quantities::PressureValue{target}}), physics::boundary_values<Pressure>({0.0, 0.2, 0.8}));
|
|
};
|
|
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(physics::make_pressure_surface_graph<Formulation>(), inputs(0.2));
|
|
REQUIRE(cache.try_prepare(inputs(0.2)));
|
|
const auto& values = cache.get<physics::keys::Linearization>();
|
|
REQUIRE(values.size() == 3);
|
|
CHECK(values[0].residual == -0.2);
|
|
CHECK(values[1].residual == 0.0);
|
|
CHECK_THAT(values[2].residual, Catch::Matchers::WithinAbs(0.6, 1.0e-15));
|
|
for (const auto& value : values)
|
|
{
|
|
CHECK(value.field_derivatives.get<Pressure>() == 1.0);
|
|
}
|
|
|
|
REQUIRE(cache.try_prepare(inputs(0.8)));
|
|
CHECK(cache.get<physics::keys::Linearization>()[2].residual == 0.0);
|
|
const auto unchanged = cache.try_prepare(inputs(0.8));
|
|
REQUIRE(unchanged);
|
|
CHECK(unchanged->computed_count() == 0);
|
|
}
|
|
|
|
TEST_CASE("Automatic pressure-field boundaries preserve interior rows and exact transpose actions", "[surface][automatic][backend]")
|
|
{
|
|
namespace fixture = serif::tests::surface::boundary;
|
|
namespace coefficients = operators::backend::mfem;
|
|
using Fields = dependencies::KeyList<Pressure>;
|
|
using Formulation = physics::AutomaticPressureSurface<SymbolicEOS, Fields, Fields>::Formulation;
|
|
|
|
fixture::Fixture mesh;
|
|
discretization::backend::mfem::TermBinding<Pressure, fixture::Schema> binding{mesh.linear_space};
|
|
const ::mfem::Vector pressure = fixture::sequence(binding.size(), 0.4, 0.2);
|
|
const ::mfem::Vector direction = fixture::sequence(binding.size(), -0.1, 0.7);
|
|
const ::mfem::Vector dual = fixture::sequence(binding.size(), 0.3, 0.4);
|
|
const auto state = coefficients::make_term_coefficients(coefficients::bind(binding, pressure));
|
|
const auto directions = coefficients::make_term_coefficients(coefficients::bind(binding, direction));
|
|
const auto boundary = backend::make_pressure_boundary<Formulation>(state);
|
|
const SymbolicEOS equation_of_state;
|
|
const auto eos_input = dependencies::borrow<physics::keys::EquationOfState>(equation_of_state, dependencies::track_by([](const SymbolicEOS&)
|
|
{
|
|
return 0;
|
|
}));
|
|
const auto inputs = backend::make_pressure_boundary_inputs(boundary, eos_input, physics::Isobaric{quantities::PressureValue{0.2}}, state);
|
|
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(backend::make_pressure_boundary_graph<Formulation>(), inputs);
|
|
REQUIRE(cache.try_prepare(inputs));
|
|
|
|
const auto& linearization = cache.get<physics::keys::Linearization>();
|
|
::mfem::Vector residual(binding.size()), action(binding.size()), transpose(binding.size());
|
|
residual = 4.0;
|
|
action = 0.0;
|
|
transpose = 0.0;
|
|
backend::apply_residual_rows(*boundary, linearization, residual);
|
|
backend::apply_jacobian_rows(*boundary, linearization, directions, action);
|
|
backend::add_jacobian_transpose_rows<Pressure>(*boundary, linearization, dual, transpose);
|
|
const ::mfem::Vector selected = fixture::boundary_mask(binding);
|
|
|
|
for (int row = 0; row < binding.size(); ++row)
|
|
{
|
|
if (selected[row] != 0.0)
|
|
{
|
|
CHECK_THAT(residual[row], Catch::Matchers::WithinAbs(pressure[row] - 0.2, 1.0e-14));
|
|
CHECK_THAT(action[row], Catch::Matchers::WithinAbs(direction[row], 1.0e-14));
|
|
CHECK_THAT(transpose[row], Catch::Matchers::WithinAbs(dual[row], 1.0e-14));
|
|
}
|
|
else
|
|
{
|
|
CHECK(residual[row] == 4.0);
|
|
CHECK(action[row] == 0.0);
|
|
CHECK(transpose[row] == 0.0);
|
|
}
|
|
}
|
|
CHECK_THAT(action * dual, Catch::Matchers::WithinAbs(direction * transpose, 1.0e-14));
|
|
}
|