tboudreaux/MeanField_SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
MeanField_SERiF/tests/surface/boundary.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

545 lines
31 KiB
C++
Raw Permalink Blame History

#include <algorithm>
#include <cmath>
#include <limits>
#include <memory>
#include <stdexcept>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include "serif/backend/mfem/geometry/surface/pressure_boundary.hpp"
#include "serif/tests/surface/boundary_fixture.hpp"
#include "serif/utils/dependencies/computation_cache.hpp"
namespace
{
namespace test = serif::tests::surface;
namespace fixture = test::boundary;
namespace physics = serif::geometry::surface;
namespace boundary = physics::backend::mfem;
namespace dependencies = serif::utils::dependencies;
namespace coefficients = serif::operators::backend::mfem;
using fixture::Vector;
using Status = dependencies::PreparationStatus;
using Linearization = physics::keys::Linearization;
using Response = physics::keys::EOSResponse;
using Catch::Matchers::WithinAbs;
/** @brief Track physical EOS coefficients while excluding observational evaluation counters. */
[[nodiscard]] auto eos_input(const test::GeneralEOS& model)
{
return dependencies::borrow<physics::keys::EquationOfState>(model, dependencies::track_by([](const test::GeneralEOS& eos)
{
return eos.scale;
}));
}
/** @brief The ideal-gas fixture has no mutable physical parameters. */
[[nodiscard]] auto eos_input(const test::GasEOS& model)
{
return dependencies::borrow<physics::keys::EquationOfState>(model, dependencies::track_by([](const test::GasEOS&)
{
return 0;
}));
}
/** @brief Make the prescribed thermodynamic pressure explicit in every candidate. */
[[nodiscard]] physics::IsobaricSurface pressure(const double value)
{
return physics::IsobaricSurface{test::quantities::PressureValue{value}};
}
/** @brief Bind the composition-dependent EOS to three distinct finite-element fields. */
[[nodiscard]] auto composition_state(const fixture::Fixture& mesh, const Vector& enthalpy, const Vector& entropy, const Vector& composition)
{
return coefficients::make_term_coefficients(coefficients::bind(mesh.enthalpy, enthalpy), coefficients::bind(mesh.entropy, entropy), coefficients::bind(mesh.composition, composition));
}
/** @brief Check the same nonlinear gas surface with either carrier-resolved or direct-pressure residuals. */
template <typename Formulation>
void check_gas_rows()
{
fixture::Fixture mesh;
test::GasEOS model;
Vector density = fixture::sequence(mesh.density.size(), 0.6, 0.1);
Vector temperature = fixture::sequence(mesh.temperature.size(), 1.3, 0.2);
const Vector density_direction = fixture::sequence(density.Size(), 0.07, 0.3);
const Vector temperature_direction = fixture::sequence(temperature.Size(), -0.11, 0.4);
const auto state = coefficients::make_term_coefficients(coefficients::bind(mesh.density, density), coefficients::bind(mesh.temperature, temperature));
const auto direction = coefficients::make_term_coefficients(coefficients::bind(mesh.density, density_direction), coefficients::bind(mesh.temperature, temperature_direction));
const auto selected = boundary::make_pressure_boundary<Formulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.8), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<Formulation>(), inputs);
REQUIRE(cache.try_prepare(inputs));
const Vector sampled_temperature = fixture::reference_samples(mesh.density, *selected->selected_rows, mesh.temperature, temperature);
const Vector sampled_direction = fixture::reference_samples(mesh.density, *selected->selected_rows, mesh.temperature, temperature_direction);
const auto& physics = cache.template get<Linearization>();
Vector residual(density.Size()), action(density.Size());
residual = 4.0;
action = -3.0;
boundary::apply_residual_rows(*selected, physics, residual);
boundary::apply_jacobian_rows(*selected, physics, direction, action);
const Vector mask = fixture::boundary_mask(mesh.density);
for (int index = 0; index < selected->selected_rows->size(); ++index)
{
const int row = selected->selected_rows->reduced_dofs()[index];
const double local_temperature = sampled_temperature[index];
const double expected_residual = Formulation::resolves_carrier ? density[row] - 0.8 / local_temperature : density[row] * local_temperature - 0.8;
const double expected_action = Formulation::resolves_carrier ? density_direction[row] + 0.8 * sampled_direction[index] / (local_temperature * local_temperature) : local_temperature * density_direction[row] + density[row] * sampled_direction[index];
CHECK_THAT(residual[row], WithinAbs(expected_residual, 1.0e-13));
CHECK_THAT(action[row], WithinAbs(expected_action, 1.0e-13));
}
for (int row = 0; row < mask.Size(); ++row)
{
if (mask[row] == 0.0)
{
CHECK(residual[row] == 4.0);
CHECK(action[row] == -3.0);
action[row] = 0.0;
}
}
const Vector dual = fixture::sequence(density.Size(), 0.4, 0.8);
Vector density_force(density.Size()), temperature_force(temperature.Size());
density_force = 0.0;
temperature_force = 0.0;
boundary::add_jacobian_transpose_rows<fixture::BoundaryDensity>(*selected, physics, dual, density_force);
boundary::add_jacobian_transpose_rows<test::Temperature>(*selected, physics, dual, temperature_force);
const double forward_work = action * dual;
const double reverse_work = density_direction * density_force + temperature_direction * temperature_force;
CHECK_THAT(forward_work, WithinAbs(reverse_work, 2.0e-13 * std::max(1.0, std::abs(forward_work))));
constexpr double step = 1.0e-7;
const Vector original_density(density), original_temperature(temperature);
Vector plus(density.Size()), minus(density.Size());
plus = 0.0;
minus = 0.0;
density.Add(step, density_direction);
temperature.Add(step, temperature_direction);
REQUIRE(cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.8), state)));
boundary::apply_residual_rows(*selected, cache.template get<Linearization>(), plus);
density = original_density;
temperature = original_temperature;
density.Add(-step, density_direction);
temperature.Add(-step, temperature_direction);
REQUIRE(cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.8), state)));
boundary::apply_residual_rows(*selected, cache.template get<Linearization>(), minus);
for (const int row : selected->selected_rows->reduced_dofs())
{
CHECK_THAT((plus[row] - minus[row]) / (2.0 * step), WithinAbs(action[row], 2.0e-9));
}
}
}
TEST_CASE("General EOS pressure rows sample every field in its own finite-element basis", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GeneralEOS model;
Vector enthalpy = fixture::sequence(mesh.enthalpy.size(), 1.7, 0.1);
Vector entropy = fixture::sequence(mesh.entropy.size(), 0.2, 0.3);
Vector composition = fixture::sequence(mesh.composition.size(), 0.1, 0.5);
const Vector enthalpy_direction = fixture::sequence(enthalpy.Size(), 0.7, 0.2);
const Vector entropy_direction = fixture::sequence(entropy.Size(), -0.2, 0.4);
const Vector composition_direction = fixture::sequence(composition.Size(), 0.05, 0.6);
const auto state = composition_state(mesh, enthalpy, entropy, composition);
const auto direction = composition_state(mesh, enthalpy_direction, entropy_direction, composition_direction);
const auto selected = boundary::make_pressure_boundary<fixture::CompositionFormulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<fixture::CompositionFormulation>(), inputs);
CHECK(model.evaluations == 0);
REQUIRE(cache.try_prepare(inputs));
const auto& accepted = cache.get<Linearization>();
const auto& rows = *selected->selected_rows;
REQUIRE(rows.size() > 0);
REQUIRE(rows.size() < enthalpy.Size());
const Vector mask = fixture::boundary_mask(mesh.enthalpy);
CHECK(rows.size() == static_cast<int>(mask.Sum()));
const Vector sampled_entropy = fixture::reference_samples(mesh.enthalpy, rows, mesh.entropy, entropy);
const Vector sampled_composition = fixture::reference_samples(mesh.enthalpy, rows, mesh.composition, composition);
const Vector sampled_entropy_direction = fixture::reference_samples(mesh.enthalpy, rows, mesh.entropy, entropy_direction);
const Vector sampled_composition_direction = fixture::reference_samples(mesh.enthalpy, rows, mesh.composition, composition_direction);
Vector residual(enthalpy.Size()), action(enthalpy.Size());
residual = 7.0;
action = -4.0;
boundary::apply_residual_rows(*selected, accepted, residual);
boundary::apply_jacobian_rows(*selected, accepted, direction, action);
for (int index = 0; index < rows.size(); ++index)
{
const int row = rows.reduced_dofs()[index];
CHECK(mask[row] == 1.0);
CHECK_THAT(residual[row], WithinAbs(enthalpy[row] - (0.8 + 3.0 * sampled_entropy[index] + 5.0 * sampled_composition[index]), 1.0e-13));
CHECK_THAT(action[row], WithinAbs(enthalpy_direction[row] - 3.0 * sampled_entropy_direction[index] - 5.0 * sampled_composition_direction[index], 1.0e-13));
}
for (int row = 0; row < mask.Size(); ++row)
{
if (mask[row] == 0.0)
{
CHECK(residual[row] == 7.0);
CHECK(action[row] == -4.0);
action[row] = 0.0;
}
}
const Vector dual = fixture::sequence(enthalpy.Size(), 0.3, 0.9);
Vector enthalpy_force(enthalpy.Size()), entropy_force(entropy.Size()), composition_force(composition.Size()), volume_dual;
enthalpy_force = 0.0;
entropy_force = 0.0;
composition_force = 0.0;
boundary::mask_residual_dual(*selected, dual, volume_dual);
boundary::add_jacobian_transpose_rows<test::Enthalpy>(*selected, accepted, dual, enthalpy_force);
boundary::add_jacobian_transpose_rows<test::Entropy>(*selected, accepted, dual, entropy_force);
boundary::add_jacobian_transpose_rows<fixture::QuadraticComposition>(*selected, accepted, dual, composition_force);
const double forward_work = dual * action;
const double reverse_work = enthalpy_direction * enthalpy_force + entropy_direction * entropy_force + composition_direction * composition_force;
CHECK_THAT(forward_work, WithinAbs(reverse_work, 2.0e-13 * std::max(1.0, std::abs(forward_work))));
for (int row = 0; row < dual.Size(); ++row)
{
CHECK(volume_dual[row] == (mask[row] != 0.0 ? 0.0 : dual[row]));
}
Vector aliased_direction(enthalpy_direction);
const auto aliased_state = composition_state(mesh, aliased_direction, entropy_direction, composition_direction);
boundary::apply_jacobian_rows(*selected, accepted, aliased_state, aliased_direction);
Vector aliased_dual(dual);
boundary::add_jacobian_transpose_rows<test::Enthalpy>(*selected, accepted, aliased_dual, aliased_dual);
for (int row = 0; row < enthalpy.Size(); ++row)
{
CHECK(aliased_direction[row] == (mask[row] != 0.0 ? action[row] : enthalpy_direction[row]));
CHECK_THAT(aliased_dual[row], WithinAbs(dual[row] + enthalpy_force[row], 1.0e-15));
}
// Mixed-order FE interpolation adds roundoff before residual subtraction; a larger step retains the original absolute tolerance.
// The local EOS regression separately preserves MeanField's original 1e-7 step.
constexpr double step = 1.0e-6;
const Vector original_enthalpy(enthalpy), original_entropy(entropy), original_composition(composition);
Vector plus(enthalpy.Size()), minus(enthalpy.Size());
plus = 0.0;
minus = 0.0;
enthalpy.Add(step, enthalpy_direction);
entropy.Add(step, entropy_direction);
composition.Add(step, composition_direction);
REQUIRE(cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state)));
boundary::apply_residual_rows(*selected, cache.get<Linearization>(), plus);
enthalpy = original_enthalpy;
entropy = original_entropy;
composition = original_composition;
enthalpy.Add(-step, enthalpy_direction);
entropy.Add(-step, entropy_direction);
composition.Add(-step, composition_direction);
REQUIRE(cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state)));
boundary::apply_residual_rows(*selected, cache.get<Linearization>(), minus);
for (const int row : rows.reduced_dofs())
{
CHECK_THAT((plus[row] - minus[row]) / (2.0 * step), WithinAbs(action[row], 2.0e-9));
}
}
TEST_CASE("A continuous density field enforces isobaric rows without enthalpy", "[surface][unit][boundary]")
{
SECTION("Density from pressure and temperature")
{
check_gas_rows<fixture::DensityFormulation<>>();
}
SECTION("Pressure directly from density and temperature")
{
check_gas_rows<fixture::PressureFormulation<>>();
}
}
TEST_CASE("Uniform EOS inputs scatter their transpose contributions into one global scalar", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GasEOS model;
Vector density = fixture::sequence(mesh.density.size(), 0.6, 0.1);
Vector temperature(1), density_direction = fixture::sequence(density.Size(), 0.1, 0.4), temperature_direction(1);
temperature[0] = 1.4;
temperature_direction[0] = -0.2;
const auto state = coefficients::make_term_coefficients(coefficients::bind(mesh.density, density), coefficients::bind(mesh.uniform_temperature, temperature));
const auto direction = coefficients::make_term_coefficients(coefficients::bind(mesh.density, density_direction), coefficients::bind(mesh.uniform_temperature, temperature_direction));
using Formulation = fixture::DensityFormulation<fixture::UniformTemperature>;
const auto selected = boundary::make_pressure_boundary<Formulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.8), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<Formulation>(), inputs);
REQUIRE(cache.try_prepare(inputs));
const auto& accepted = cache.get<Linearization>();
Vector residual(density.Size()), action(density.Size()), density_force(density.Size()), temperature_force(1);
residual = 0.0;
action = 0.0;
density_force = 0.0;
temperature_force = 0.0;
boundary::apply_residual_rows(*selected, accepted, residual);
boundary::apply_jacobian_rows(*selected, accepted, direction, action);
const Vector dual = fixture::sequence(density.Size(), 0.4, 0.8);
boundary::add_jacobian_transpose_rows<fixture::BoundaryDensity>(*selected, accepted, dual, density_force);
boundary::add_jacobian_transpose_rows<fixture::UniformTemperature>(*selected, accepted, dual, temperature_force);
double expected_temperature_force = 0.0;
for (const int row : selected->selected_rows->reduced_dofs())
{
CHECK_THAT(residual[row], WithinAbs(density[row] - 0.8 / temperature[0], 1.0e-13));
expected_temperature_force += 0.8 / (temperature[0] * temperature[0]) * dual[row];
}
CHECK_THAT(temperature_force[0], WithinAbs(expected_temperature_force, 1.0e-13));
CHECK_THAT(action * dual, WithinAbs(density_direction * density_force + temperature_direction * temperature_force, 2.0e-13));
}
TEST_CASE("Boundary field sampling and EOS work have distinct dependency lifetimes", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GeneralEOS model;
Vector enthalpy = fixture::sequence(mesh.enthalpy.size(), 1.7, 0.1);
Vector entropy = fixture::sequence(mesh.entropy.size(), 0.2, 0.3);
Vector composition = fixture::sequence(mesh.composition.size(), 0.1, 0.5);
const auto state = composition_state(mesh, enthalpy, entropy, composition);
const auto selected = boundary::make_pressure_boundary<fixture::CompositionFormulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<fixture::CompositionFormulation>(), inputs);
CHECK(model.evaluations == 0);
REQUIRE(cache.try_prepare(inputs));
const int evaluations = model.evaluations;
REQUIRE(evaluations > 0);
const auto repeated = cache.try_prepare(inputs);
REQUIRE(repeated);
CHECK(repeated->computed_count() == 0);
enthalpy[selected->selected_rows->reduced_dofs()[0]] += 0.1;
const auto carrier_changed = cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state));
REQUIRE(carrier_changed);
CHECK(carrier_changed->status<boundary::keys::BoundarySamples<test::Enthalpy>>() == Status::computed);
CHECK(carrier_changed->status<boundary::keys::BoundarySamples<test::Entropy>>() == Status::reused);
CHECK(carrier_changed->status<Response>() == Status::reused);
CHECK(carrier_changed->status<Linearization>() == Status::computed);
CHECK(model.evaluations == evaluations);
entropy *= 1.1;
const auto entropy_changed = cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state));
REQUIRE(entropy_changed);
CHECK(entropy_changed->status<boundary::keys::BoundarySamples<test::Entropy>>() == Status::computed);
CHECK(entropy_changed->status<Response>() == Status::computed);
const auto pressure_changed = cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.5), state));
REQUIRE(pressure_changed);
CHECK(pressure_changed->status<Response>() == Status::computed);
CHECK(pressure_changed->status<boundary::keys::BoundarySamples<test::Entropy>>() == Status::reused);
model.scale = 1.2;
const auto eos_changed = cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.5), state));
REQUIRE(eos_changed);
CHECK(eos_changed->status<Response>() == Status::computed);
CHECK(eos_changed->status<boundary::keys::BoundarySamples<fixture::QuadraticComposition>>() == Status::reused);
CHECK(cache.get<Linearization>()[0].field_derivatives.get<test::Entropy>() == -3.0 * model.scale);
}
TEST_CASE("Direct-pressure boundary changes reuse EOS samples when only prescribed pressure changes", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GasEOS model;
Vector density = fixture::sequence(mesh.density.size(), 0.6, 0.1);
Vector temperature = fixture::sequence(mesh.temperature.size(), 1.3, 0.2);
const auto state = coefficients::make_term_coefficients(coefficients::bind(mesh.density, density), coefficients::bind(mesh.temperature, temperature));
using Formulation = fixture::PressureFormulation<>;
const auto selected = boundary::make_pressure_boundary<Formulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.8), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<Formulation>(), inputs);
REQUIRE(cache.try_prepare(inputs));
const int evaluations = model.evaluations;
const double previous = cache.get<Linearization>()[0].residual;
const auto changed = cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.9), state));
REQUIRE(changed);
CHECK(changed->status<Response>() == Status::reused);
CHECK(changed->status<Linearization>() == Status::computed);
CHECK(model.evaluations == evaluations);
CHECK_THAT(cache.get<Linearization>()[0].residual, WithinAbs(previous - 0.1, 1.0e-15));
}
TEST_CASE("A rejected boundary candidate hides old results and permits a valid retry", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GeneralEOS model;
Vector enthalpy = fixture::sequence(mesh.enthalpy.size(), 1.7, 0.1);
Vector entropy = fixture::sequence(mesh.entropy.size(), 0.2, 0.3);
Vector composition = fixture::sequence(mesh.composition.size(), 0.1, 0.5);
const auto state = composition_state(mesh, enthalpy, entropy, composition);
const auto selected = boundary::make_pressure_boundary<fixture::CompositionFormulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<fixture::CompositionFormulation>(), inputs);
REQUIRE(cache.try_prepare(inputs));
const double previous = cache.get<Linearization>()[0].residual;
const Vector original_entropy(entropy);
entropy = std::numeric_limits<double>::quiet_NaN();
const auto rejected = cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state));
REQUIRE_FALSE(rejected);
CHECK_FALSE(cache.ready());
CHECK_THROWS_AS(cache.get<Linearization>(), std::logic_error);
entropy = original_entropy;
REQUIRE(cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state)));
CHECK(cache.get<Linearization>()[0].residual == previous);
}
TEST_CASE("Pressure boundaries reject stale and unrelated coefficient bindings before publishing rows", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GeneralEOS model;
Vector enthalpy = fixture::sequence(mesh.enthalpy.size(), 1.7, 0.1);
Vector entropy = fixture::sequence(mesh.entropy.size(), 0.2, 0.3);
Vector composition = fixture::sequence(mesh.composition.size(), 0.1, 0.5);
const auto state = composition_state(mesh, enthalpy, entropy, composition);
const auto selected = boundary::make_pressure_boundary<fixture::CompositionFormulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<fixture::CompositionFormulation>(), inputs);
REQUIRE(cache.try_prepare(inputs));
SECTION("A same-sized binding cannot silently replace the compiled field")
{
fixture::finite_elements::TermBinding<test::Entropy, fixture::Schema> other_entropy(mesh.linear_space);
const auto other_state = coefficients::make_term_coefficients(coefficients::bind(mesh.enthalpy, enthalpy), coefficients::bind(other_entropy, entropy), coefficients::bind(mesh.composition, composition));
const auto other_inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), other_state);
CHECK_THROWS_AS(cache.try_prepare(other_inputs), std::invalid_argument);
CHECK_FALSE(cache.ready());
REQUIRE(cache.try_prepare(inputs));
}
SECTION("The candidate vector must match its finite-element coefficients")
{
Vector short_entropy(entropy.Size() - 1);
short_entropy = 0.2;
const auto other_state = composition_state(mesh, enthalpy, short_entropy, composition);
CHECK_THROWS_AS(cache.try_prepare(boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), other_state)), std::invalid_argument);
CHECK_FALSE(cache.ready());
REQUIRE(cache.try_prepare(inputs));
}
SECTION("Refinement invalidates traces even if all coefficients are unchanged")
{
mesh.mesh.UniformRefinement();
CHECK_THROWS_AS(cache.try_prepare(inputs), std::logic_error);
CHECK_FALSE(cache.ready());
CHECK_THROWS_AS(cache.get<Linearization>(), std::logic_error);
}
SECTION("Boundary labels invalidate the selected pressure rows without a topology revision")
{
const int previous = mesh.mesh.GetBdrAttribute(0);
mesh.mesh.GetBdrElement(0)->SetAttribute(previous == 1 ? 2 : 1);
mesh.mesh.SetAttributes();
CHECK_THROWS_AS(cache.try_prepare(inputs), std::logic_error);
CHECK_FALSE(cache.ready());
mesh.mesh.GetBdrElement(0)->SetAttribute(previous);
mesh.mesh.SetAttributes();
REQUIRE(cache.try_prepare(inputs));
}
SECTION("Material labels invalidate the supported coefficient binding")
{
mesh.mesh.GetElement(0)->SetAttribute(2);
mesh.mesh.SetAttributes();
CHECK_THROWS_AS(cache.try_prepare(inputs), std::logic_error);
CHECK_FALSE(cache.ready());
mesh.mesh.GetElement(0)->SetAttribute(1);
mesh.mesh.SetAttributes();
REQUIRE(cache.try_prepare(inputs));
}
SECTION("Curved coordinate edits preserve canonical finite-element traces")
{
const double previous = cache.get<Linearization>()[0].residual;
const int evaluations = model.evaluations;
const Vector original_nodes(*mesh.mesh.GetNodes());
::mfem::VectorFunctionCoefficient warp(3, [](const Vector& position, Vector& mapped)
{
mapped = position;
mapped[0] += 0.04 * position[1] * position[1];
mapped[2] += 0.03 * position[0] * position[1];
});
mesh.mesh.Transform(warp);
Vector coordinate_change(*mesh.mesh.GetNodes());
coordinate_change -= original_nodes;
REQUIRE(coordinate_change.Norml2() > 0.0);
const auto reused = cache.try_prepare(inputs);
REQUIRE(reused);
CHECK(reused->computed_count() == 0);
CHECK(model.evaluations == evaluations);
CHECK(cache.get<Linearization>()[0].residual == previous);
}
}
TEST_CASE("Pressure row actions reject malformed buffers before replacing any entries", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GeneralEOS model;
Vector enthalpy = fixture::sequence(mesh.enthalpy.size(), 1.7, 0.1);
Vector entropy = fixture::sequence(mesh.entropy.size(), 0.2, 0.3);
Vector composition = fixture::sequence(mesh.composition.size(), 0.1, 0.5);
const auto state = composition_state(mesh, enthalpy, entropy, composition);
const auto selected = boundary::make_pressure_boundary<fixture::CompositionFormulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(selected, eos_input(model), pressure(0.4), state);
auto cache = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<fixture::CompositionFormulation>(), inputs);
REQUIRE(cache.try_prepare(inputs));
const auto& accepted = cache.get<Linearization>();
Vector result(enthalpy.Size());
result = 3.0;
SECTION("Missing point linearizations cannot index outside the accepted batch")
{
auto incomplete = accepted;
incomplete.pop_back();
CHECK_THROWS_AS(boundary::apply_residual_rows(*selected, incomplete, result), std::invalid_argument);
CHECK_THROWS_AS(boundary::apply_jacobian_rows(*selected, incomplete, state, result), std::invalid_argument);
CHECK_THROWS_AS(boundary::add_jacobian_transpose_rows<test::Enthalpy>(*selected, incomplete, enthalpy, result), std::invalid_argument);
}
SECTION("Short directions cannot index outside an input field")
{
Vector short_entropy(entropy.Size() - 1);
short_entropy = 0.2;
const auto direction = composition_state(mesh, enthalpy, short_entropy, composition);
CHECK_THROWS_AS(boundary::apply_jacobian_rows(*selected, accepted, direction, result), std::invalid_argument);
}
SECTION("A nonfinite direction cannot partially replace the rows")
{
entropy = std::numeric_limits<double>::infinity();
CHECK_THROWS_AS(boundary::apply_jacobian_rows(*selected, accepted, state, result), std::domain_error);
}
SECTION("Short residual output is rejected before writing")
{
result.SetSize(enthalpy.Size() - 1);
result = 3.0;
CHECK_THROWS_AS(boundary::apply_residual_rows(*selected, accepted, result), std::invalid_argument);
CHECK_THROWS_AS(boundary::apply_jacobian_rows(*selected, accepted, state, result), std::invalid_argument);
CHECK_THROWS_AS(boundary::add_jacobian_transpose_rows<test::Enthalpy>(*selected, accepted, enthalpy, result), std::invalid_argument);
}
for (int index = 0; index < result.Size(); ++index)
{
CHECK(result[index] == 3.0);
}
}
TEST_CASE("Accepted boundary snapshots retain their identity until a replacement is published", "[surface][unit][boundary]")
{
fixture::Fixture mesh;
test::GeneralEOS model;
Vector enthalpy = fixture::sequence(mesh.enthalpy.size(), 1.7, 0.1);
Vector entropy = fixture::sequence(mesh.entropy.size(), 0.2, 0.3);
Vector composition = fixture::sequence(mesh.composition.size(), 0.1, 0.5);
const auto state = composition_state(mesh, enthalpy, entropy, composition);
using BoundaryPointer = decltype(boundary::make_pressure_boundary<fixture::CompositionFormulation>(state));
std::weak_ptr<typename BoundaryPointer::element_type> previous_boundary;
auto cache = [&]()
{
const auto original = boundary::make_pressure_boundary<fixture::CompositionFormulation>(state);
previous_boundary = original;
const auto inputs = boundary::make_pressure_boundary_inputs(original, eos_input(model), pressure(0.4), state);
auto prepared = dependencies::make_computation_cache<physics::SurfaceFailure>(boundary::make_pressure_boundary_graph<fixture::CompositionFormulation>(), inputs);
REQUIRE(prepared.try_prepare(inputs));
return prepared;
}();
REQUIRE_FALSE(previous_boundary.expired());
const double previous_residual = cache.get<Linearization>()[0].residual;
const auto replacement = boundary::make_pressure_boundary<fixture::CompositionFormulation>(state);
const auto inputs = boundary::make_pressure_boundary_inputs(replacement, eos_input(model), pressure(0.4), state);
const auto changed = cache.try_prepare(inputs);
REQUIRE(changed);
CHECK(changed->status<boundary::keys::BoundarySamples<test::Enthalpy>>() == Status::computed);
CHECK(changed->status<boundary::keys::BoundarySamples<test::Entropy>>() == Status::computed);
CHECK(changed->status<boundary::keys::BoundarySamples<fixture::QuadraticComposition>>() == Status::computed);
CHECK(changed->status<Response>() == Status::computed);
CHECK(changed->status<Linearization>() == Status::computed);
CHECK(cache.get<Linearization>()[0].residual == previous_residual);
CHECK(previous_boundary.expired());
}
Reference in New Issue View Git Blame Copy Permalink