From 30af52cceadb87902d6359b40549d8205970818d Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 9 Nov 2026 13:13:53 -0500 Subject: [PATCH] feat(geometry/surface): System for prescripting surface condition and building computational node for surface condition The surface needs some defined condition (e.g. Isobaric) This module defines how the surface is defined. Again, much like the deformation system this does not actually implement much of the code to deform the surface to that condition; rather, it implements code to allow the prescription to be described robustly at compile time --- .../serif/geometry/surface/automatic.hpp | 83 +++++++++++ .../serif/geometry/surface/compilation.hpp | 108 ++++++++++++++ .../serif/geometry/surface/computations.hpp | 108 ++++++++++++++ .../surface/detail/automatic_relations.hpp | 87 +++++++++++ .../serif/geometry/surface/detail/batches.hpp | 114 +++++++++++++++ .../geometry/surface/detail/evaluation.hpp | 135 ++++++++++++++++++ .../geometry/surface/detail/relations.hpp | 97 +++++++++++++ .../serif/geometry/surface/evaluation.hpp | 26 ++++ .../serif/geometry/surface/isobaric.hpp | 40 ++++++ .../serif/geometry/surface/linearization.hpp | 41 ++++++ .../serif/geometry/surface/samples.hpp | 70 +++++++++ .../serif/geometry/surface/surface.hpp | 17 +++ src/lib/serif/geometry/surface/isobaric.cpp | 31 ++++ 13 files changed, 957 insertions(+) create mode 100644 src/include/serif/geometry/surface/automatic.hpp create mode 100644 src/include/serif/geometry/surface/compilation.hpp create mode 100644 src/include/serif/geometry/surface/computations.hpp create mode 100644 src/include/serif/geometry/surface/detail/automatic_relations.hpp create mode 100644 src/include/serif/geometry/surface/detail/batches.hpp create mode 100644 src/include/serif/geometry/surface/detail/evaluation.hpp create mode 100644 src/include/serif/geometry/surface/detail/relations.hpp create mode 100644 src/include/serif/geometry/surface/evaluation.hpp create mode 100644 src/include/serif/geometry/surface/linearization.hpp create mode 100644 src/include/serif/geometry/surface/samples.hpp create mode 100644 src/include/serif/geometry/surface/surface.hpp diff --git a/src/include/serif/geometry/surface/automatic.hpp b/src/include/serif/geometry/surface/automatic.hpp new file mode 100644 index 00000000..bdef120b --- /dev/null +++ b/src/include/serif/geometry/surface/automatic.hpp @@ -0,0 +1,83 @@ +#pragma once + +#include "serif/geometry/surface/detail/automatic_relations.hpp" + +namespace serif::geometry::surface { + enum class PressureSurfaceSelection { + supported, + unavailable, + ambiguous + }; + + namespace detail { + /** @brief Strong boundary enforcement requires an actual solved, scalar H1 equation row. */ + template + concept EligibleSurfaceRow = discretization::blocks::SpatialTerm && discretization::blocks::HasResidualBlock && + utils::dependencies::ContainsKey && utils::dependencies::ContainsKey && + (std::same_as || std::same_as); + + /** @brief Search one precedence level, without asking the EOS for a relation when pressure is already an unknown. */ + template + [[nodiscard]] consteval auto field_candidates() { + if constexpr (!EligibleSurfaceRow) { + return std::tuple<>{}; + } + + else if constexpr (Priority == 0) { + if constexpr (std::same_as) { + using Bindings = utils::dependencies::Bindings>; + return std::tuple>>{}; + } + else { + return std::tuple<>{}; + } + } + else if constexpr (eos::models::EOSModel) { + return relation_candidates(typename EOS::Relations{}); + } + else { + return std::tuple<>{}; + } + } + + /** @brief Prefer a pressure unknown, then an explicit inverse, then a direct pressure residual. */ + template + [[nodiscard]] consteval auto select_candidates(utils::dependencies::KeyList fields) { + using FieldsList = utils::dependencies::KeyList; + using Candidates = decltype(std::tuple_cat(field_candidates()...)); + + if constexpr (std::tuple_size_v != 0 || Priority == 2) { + return Candidates{}; + } + else { + return select_candidates(fields); + } + } + + template + struct SelectedPressureSurface { + static constexpr PressureSurfaceSelection status = std::tuple_size_v == 0 ? PressureSurfaceSelection::unavailable : PressureSurfaceSelection::ambiguous; + }; + + template + struct SelectedPressureSurface> { + static constexpr PressureSurfaceSelection status = PressureSurfaceSelection::supported; + using Formulation = typename Candidate::type; + }; + } + + /** + * @brief Compile an isobaric condition from available physical fields, actual equation rows, and the EOS catalog. + * @note This searches declared algebraic relations, never performs inversion, and rejects equal-priority choices. + * Global thermodynamic inputs may opt in with thermodynamic_state = true; units alone do not establish that role. + * An optional preferred field narrows the row search without choosing an EOS relation for the model user. + */ + template + struct AutomaticPressureSurface : detail::SelectedPressureSurface, ActualRows, SolvedTerms, PreferredField>(detail::thermodynamic_fields(AvailableTerms{})))> { + }; + + /** @brief Whether the available equations and numerical backend select one implementable isobaric residual. */ + template + concept AutomaticPressureSurfaceCompilable = requires { typename AutomaticPressureSurface::Formulation; }; +} diff --git a/src/include/serif/geometry/surface/compilation.hpp b/src/include/serif/geometry/surface/compilation.hpp new file mode 100644 index 00000000..2e5d3f25 --- /dev/null +++ b/src/include/serif/geometry/surface/compilation.hpp @@ -0,0 +1,108 @@ +#pragma once + +#include "serif/geometry/surface/detail/relations.hpp" +#include "serif/geometry/surface/isobaric.hpp" + +namespace serif::geometry::surface +{ + /** + * @brief Whether the available state selects exactly one pressure-to-carrier EOS relation. + * @tparam Carrier Optional physical quantity whose equation enforces the boundary condition. + * @note An ambiguous catalog never selects its first entry. Specifying a carrier narrows the same search. + */ + template + concept PressureSurfaceCompilable = eos::models::EOSModel && SurfaceStateBindings && + (detail::MatchingPressureRelations, std::remove_cvref_t, Carrier, typename std::remove_cvref_t::Relations>::count == 1); + + template + requires PressureSurfaceCompilable + using pressure_relation_t = std::tuple_element_t<0, typename detail::MatchingPressureRelations, std::remove_cvref_t, Carrier, typename std::remove_cvref_t::Relations>::Matches>; + + namespace detail + { + /** @brief Common physical metadata; numerical values and preparation belong to the dependency context and cache. */ + template + struct PressureSurfaceDescription + { + using EquationOfState = std::remove_cvref_t; + using StateBindings = std::remove_cvref_t; + using Relation = RelationT; + using PhysicalQuantity = dimensions::Pressure; + using ResidualQuantity = std::conditional_t; + using CarrierQuantity = RowQuantity; + using CarrierTerm = utils::dependencies::bound_type_t; + using Test = CarrierTerm; + using Support = discretization::domain::StellarSurfaceBoundary; + using EOSStateQuantities = std::conditional_t::NonPressureInputs, typename RelationQuantities::AllInputs>; + using StateQuantities = utils::dependencies::CombinedKeyList, EOSStateQuantities>; + using StateTerms = BoundTerms; + using EOSStateTerms = BoundTerms; + using Operands = typename SurfaceOperands::Type; + + static constexpr bool resolves_carrier = ResolveCarrier; + static constexpr bool uses_eos = true; + }; + } + + /** @brief The boundary residual is carrier minus the EOS value at prescribed pressure; no carrier quantity is privileged. */ + template + requires PressureSurfaceCompilable + struct PressureSurfaceFormulation final : detail::PressureSurfaceDescription::OutputQuantity, pressure_relation_t, true> + { + }; + + /** @brief A selected explicit inverse relation; choosing it requires no numerical EOS inversion. */ + template + requires eos::relations::CatalogContainingRelation && (detail::matches_pressure_relation()) + struct InversePressureSurfaceFormulation final : detail::PressureSurfaceDescription + { + }; + + /** @brief A direct pressure residual requires an explicitly chosen pressure relation and an input quantity whose equation row it replaces. */ + template + concept DirectPressureSurfaceCompilable = eos::models::EOSModel && SurfaceStateBindings && eos::relations::ValidThermodynamicRelation && + eos::relations::CatalogContainingRelation::Relations, Relation> && + std::same_as && eos::relations::RelationContainingInput && + SurfaceResidualRow && + utils::dependencies::SatisfiesKeyList::AllInputs> && + detail::supports_partials(typename detail::RelationQuantities::AllInputs{}); + + /** + * @brief Enforce P_EOS(state) - P_surface = 0 without inventing an inverse EOS relation. + * @note The explicit row and relation make the choice of residual, units, and Newton scaling visible. + */ + template + requires DirectPressureSurfaceCompilable + struct DirectPressureSurfaceFormulation final : detail::PressureSurfaceDescription + { + }; + + /** @brief An existing pressure unknown enforces P - P_surface directly, without evaluating the EOS. */ + template + requires SurfaceResidualRow + struct PressureFieldSurfaceFormulation final : detail::PressureSurfaceDescription, true> + { + static constexpr bool uses_eos = false; + }; + + namespace detail + { + template + inline constexpr bool is_pressure_surface_formulation = false; + + template + inline constexpr bool is_pressure_surface_formulation> = true; + + template + inline constexpr bool is_pressure_surface_formulation> = true; + + template + inline constexpr bool is_pressure_surface_formulation> = true; + + template + inline constexpr bool is_pressure_surface_formulation> = true; + } + + template + concept PressureSurfaceFormulationType = detail::is_pressure_surface_formulation; +} diff --git a/src/include/serif/geometry/surface/computations.hpp b/src/include/serif/geometry/surface/computations.hpp new file mode 100644 index 00000000..4d432a2f --- /dev/null +++ b/src/include/serif/geometry/surface/computations.hpp @@ -0,0 +1,108 @@ +#pragma once + +#include +#include +#include + +#include "serif/geometry/surface/detail/batches.hpp" +#include "serif/utils/dependencies/computation_graph.hpp" + +namespace serif::geometry::surface +{ + namespace detail + { + /** @brief Declare the exact EOS inputs from the formulation, so derivative metadata and cache dependencies cannot diverge. */ + template typename SampleKey, typename... Terms> + [[nodiscard]] auto make_surface_eos_computation(utils::dependencies::KeyList) + { + using namespace utils::dependencies; + using EOS = typename Formulation::EquationOfState; + + if constexpr (!Formulation::uses_eos) + { + // Reuse the carrier residual stage with one constant target; pressure unknowns need no EOS dependency. + return make_computation>( + [](const IsobaricSurface& condition) + { + return std::vector>{{.value = condition.target_pressure()}}; + }); + } + else if constexpr (Formulation::resolves_carrier) + { + return make_computation...>>( + [](const EOS& equation_of_state, const IsobaricSurface& condition, const BoundaryValues&... values) + { + return evaluate_response_samples(equation_of_state, condition.target_pressure(), values...); + }); + } + else + { + return make_computation...>>( + [](const EOS& equation_of_state, const BoundaryValues&... values) + { + return evaluate_response_samples(equation_of_state, dimensions::PressureValue{0.0}, values...); + }); + } + } + + /** @brief Keep the inexpensive boundary residual separate from the reusable EOS response and its partials. */ + template typename SampleKey> + [[nodiscard]] auto make_surface_linearization_computation() + { + using namespace utils::dependencies; + using Responses = std::vector>; + + if constexpr (Formulation::resolves_carrier) + { + using Carrier = typename Formulation::CarrierTerm; + return make_computation>>( + [](const Responses& responses, const BoundaryValues& carrier) + { + return linearize_carrier_samples(responses, carrier); + }); + } + else + { + return make_computation>( + [](const Responses& responses, const IsobaricSurface& condition) + { + return linearize_pressure_samples(condition, responses); + }); + } + } + } + + /** @brief Expose the existing EOS and residual stages for composition with backend sampling or other physics. */ + template typename SampleKey = std::type_identity_t> + [[nodiscard]] auto make_pressure_surface_computations() + { + return std::tuple{ + detail::make_surface_eos_computation(typename Formulation::EOSStateTerms{}), + detail::make_surface_linearization_computation()}; + } + + /** + * @brief Describe the isobaric calculation as an ordinary dependency graph, without evaluating the EOS. + * @tparam Formulation A carrier-resolved or explicitly selected direct-pressure surface formulation. + * @return Two computations: the EOS response and the boundary residual with analytical derivatives. + * @note Supply pressure_condition and boundary_values, and bind keys::EquationOfState with a tracking policy when the formulation uses the EOS. + * A pressure-field formulation needs no EOS source. Use make_computation_cache and its normal try_prepare/get interface; + * the cache owns the complete preparation lifecycle. + */ + template + [[nodiscard]] auto make_pressure_surface_graph() + { + return std::apply([](auto... computations) + { + return utils::dependencies::make_computation_graph(std::move(computations)...); + }, make_pressure_surface_computations()); + } + + /** @brief Infer a unique carrier-resolved formulation directly from an EOS and the available physical state. */ + template + requires PressureSurfaceCompilable + [[nodiscard]] auto make_pressure_surface_graph() + { + return make_pressure_surface_graph>(); + } +} diff --git a/src/include/serif/geometry/surface/detail/automatic_relations.hpp b/src/include/serif/geometry/surface/detail/automatic_relations.hpp new file mode 100644 index 00000000..a50d7e5c --- /dev/null +++ b/src/include/serif/geometry/surface/detail/automatic_relations.hpp @@ -0,0 +1,87 @@ +#pragma once + +#include "serif/geometry/surface/compilation.hpp" +#include "serif/discretization/blocks/traits.hpp" + +namespace serif::geometry::surface::detail { + /** @brief Globals require an explicit physical role: matching units alone would mistake C or a multiplier for enthalpy. */ + template + concept ThermodynamicBoundaryField = discretization::blocks::ScalarTraceTerm && + discretization::blocks::HasQuantity && dimensions::ThermodynamicQuantityType && + (discretization::blocks::SpatialTerm || requires { requires Term::thermodynamic_state; }); + + /** @brief Collect physical fields, retaining distinct terms with equal quantities so ambiguity remains visible. */ + template + [[nodiscard]] consteval auto thermodynamic_fields(utils::dependencies::KeyList) { + return utils::dependencies::CombinedKeyList, utils::dependencies::KeyList, utils::dependencies::KeyList<>>...>{}; + } + + /** @brief Return a unique field for a quantity; the chosen residual carrier supplies its own quantity unambiguously. */ + template + [[nodiscard]] consteval auto quantity_field(utils::dependencies::KeyList) { + using Matches = decltype(std::tuple_cat(std::declval, std::tuple, std::tuple<>>>()...)); + if constexpr (std::same_as) { + return std::type_identity{}; + } + else if constexpr (std::tuple_size_v == 1) { + return std::type_identity>{}; + } + else { + return std::type_identity{}; + } + } + + /** @brief Bind only the quantities read by this candidate, ignoring unrelated fields and quantities. */ + template + [[nodiscard]] consteval auto relation_bindings(utils::dependencies::KeyList) { + using namespace utils::dependencies; + + if constexpr ((!std::same_as(Fields{}))::type> && ...)) { + return std::type_identity(Fields{}))::type>...>>{}; + } + else { + return std::type_identity{}; + } + } + + /** @brief Test a declared inverse and every required analytic partial before publishing its formulation. */ + template + [[nodiscard]] consteval auto inverse_candidate() { + using namespace utils::dependencies; + using Quantities = CombinedKeyList, typename surface::detail::RelationQuantities::NonPressureInputs>; + using Bindings = typename decltype(relation_bindings(Quantities{}))::type; + + if constexpr (SurfaceStateBindings && surface::detail::matches_pressure_relation()) { + return std::tuple>>{}; + } + else { + return std::tuple<>{}; + } + } + + /** @brief A direct pressure relation may replace only a row belonging to one of its solved inputs. */ + template + [[nodiscard]] consteval auto direct_candidate() { + using namespace utils::dependencies; + using Quantities = CombinedKeyList, typename surface::detail::RelationQuantities::AllInputs>; + using Bindings = typename decltype(relation_bindings(Quantities{}))::type; + + if constexpr (DirectPressureSurfaceCompilable) { + return std::tuple>>{}; + } + else { + return std::tuple<>{}; + } + } + + /** @brief Keep every admissible relation; catalog order never selects a physical implementation. */ + template + [[nodiscard]] consteval auto relation_candidates(eos::relations::RelationCatalog) { + if constexpr (Inverse) { + return std::tuple_cat(inverse_candidate()...); + } + else { + return std::tuple_cat(direct_candidate()...); + } + } +} diff --git a/src/include/serif/geometry/surface/detail/batches.hpp b/src/include/serif/geometry/surface/detail/batches.hpp new file mode 100644 index 00000000..44f78813 --- /dev/null +++ b/src/include/serif/geometry/surface/detail/batches.hpp @@ -0,0 +1,114 @@ +#pragma once + +#include + +#include "serif/geometry/surface/detail/evaluation.hpp" +#include "serif/geometry/surface/samples.hpp" + +namespace serif::geometry::surface::detail +{ + /** @brief All EOS columns must describe the same ordered, nonempty set of boundary points. */ + template + [[nodiscard]] std::expected surface_sample_count(const BoundaryValues&... columns) + { + if constexpr (sizeof...(Terms) == 0) + { + // A pressure-only relation has one constant response, broadcast over the carrier samples later. + return 1; + } + else + { + const std::array counts{columns.values.size()...}; + if (counts.front() == 0) + { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::empty_surface, .message = "A surface calculation requires at least one boundary sample."}); + } + for (const std::size_t count : counts) + { + if (count != counts.front()) + { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::inconsistent_sample_count, .message = "Surface terms must be sampled at the same ordered boundary points; their sample counts differ."}); + } + } + return counts.front(); + } + } + + /** @brief Assemble the conventional equation point from typed columns; no mesh assumptions enter the physical calculation. */ + template + [[nodiscard]] SurfacePoint surface_point(const std::size_t sample, const BoundaryValues&... columns) + { + SurfacePoint point{}; + ((point.fields.template get() = columns.values[sample]), ...); + return point; + } + + /** @brief Collect one ordered batch, attaching the first failure's sample index before any result is published. */ + template + [[nodiscard]] std::expected, SurfaceFailure> collect_surface_samples(std::size_t count, Evaluate evaluate) + { + std::vector result; + result.reserve(count); + for (std::size_t sample = 0; sample < count; ++sample) + { + auto value = evaluate(sample); + if (!value) + { + SurfaceFailure failure = value.error(); + failure.sample = sample; + return std::unexpected(failure); + } + result.push_back(std::move(*value)); + } + return result; + } + + /** @brief Evaluate the EOS once per state sample, or once for a genuinely pressure-only relation. */ + template + [[nodiscard]] std::expected>, SurfaceFailure> evaluate_response_samples(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const BoundaryValues&... columns) + { + const auto count = surface_sample_count(columns...); + if (!count) + { + return std::unexpected(count.error()); + } + return collect_surface_samples>(*count, [&](std::size_t sample) + { + return evaluate_surface_response(equation_of_state, pressure, surface_point(sample, columns...)); + }); + } + + /** @brief Broadcast only a pressure-only relation; every state-dependent carrier sample needs its own response. */ + template + [[nodiscard]] std::expected>, SurfaceFailure> linearize_carrier_samples(const std::vector>& responses, const BoundaryValues& carrier) + { + constexpr bool constant_response = std::same_as>; + const auto count = surface_sample_count(carrier); + if (!count) + { + return std::unexpected(count.error()); + } + if (responses.size() != (constant_response ? 1 : *count)) + { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::inconsistent_sample_count, .message = "The carrier and EOS inputs must refer to the same boundary points."}); + } + return collect_surface_samples>(*count, [&](std::size_t sample) -> std::expected, SurfaceFailure> + { + if (!std::isfinite(carrier.values[sample])) + { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_state, .message = "A boundary carrier value is non-finite."}); + } + return checked_linearization(dimensions::PressureValue{0.0}, surface_point(sample, carrier), responses[constant_response ? 0 : sample]); + }); + } + + /** @brief Changing prescribed pressure reuses all values and partials of a direct-pressure EOS relation. */ + template + [[nodiscard]] std::expected>, SurfaceFailure> linearize_pressure_samples(const IsobaricSurface& condition, const std::vector>& responses) + { + return collect_surface_samples>(responses.size(), [&](std::size_t sample) + { + return checked_linearization(condition.target_pressure(), SurfacePoint{}, responses[sample]); + }); + } +} // namespace serif::geometry::surface::detail diff --git a/src/include/serif/geometry/surface/detail/evaluation.hpp b/src/include/serif/geometry/surface/detail/evaluation.hpp new file mode 100644 index 00000000..5bb2eaca --- /dev/null +++ b/src/include/serif/geometry/surface/detail/evaluation.hpp @@ -0,0 +1,135 @@ +#pragma once + +#include +#include + +#include "serif/eos/evaluation.hpp" +#include "serif/geometry/surface/linearization.hpp" + +namespace serif::geometry::surface::detail +{ + /** @brief Read one thermodynamic role from its bound term, preserving the EOS argument's physical type. */ + template + [[nodiscard]] dimensions::QuantityValue state_value(const SurfacePoint& point) + { + using Term = utils::dependencies::bound_type_t; + return dimensions::QuantityValue{point.fields.template get()}; + } + + /** @brief Only the carrier-resolved formulation substitutes prescribed pressure for an EOS input. */ + template + [[nodiscard]] dimensions::QuantityValue relation_input(const dimensions::PressureValue pressure, const SurfacePoint& point) + { + if constexpr (Formulation::resolves_carrier && std::same_as) + { + return pressure; + } + else + { + return state_value(point); + } + } + + /** @brief Validate precisely the state read by this calculation, leaving unrelated state outside its dependencies. */ + template + [[nodiscard]] bool finite_state(const SurfacePoint& point, utils::dependencies::KeyList) + { + return (std::isfinite(state_value(point).value()) && ...); + } + + /** @brief Add each quantity's chain-rule contribution, including distinct roles bound to the same unknown. */ + template + void accumulate_partial(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const SurfacePoint& point, SurfaceEOSResponse& response, dimensions::QuantityList) + { + using Output = typename Formulation::Relation::OutputQuantity; + using Term = utils::dependencies::bound_type_t; + + const dimensions::PartialDerivative derivative = eos::partial_derivative(equation_of_state, relation_input(pressure, point)...); + response.field_derivatives.template get() += derivative.value(); + } + + /** @brief Evaluate one selected relation and all required partials without any hidden EOS inversion. */ + template + [[nodiscard]] SurfaceEOSResponse evaluate_relation(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const SurfacePoint& point, dimensions::QuantityList inputs, utils::dependencies::KeyList) + { + using Output = typename Formulation::Relation::OutputQuantity; + + SurfaceEOSResponse response{}; + response.value = eos::evaluate(equation_of_state, relation_input(pressure, point)...); + (accumulate_partial(equation_of_state, pressure, point, response, inputs), ...); + return response; + } + + /** + * @brief Reject non-finite thermodynamic inputs, EOS outputs, and partials as values. + * @note Existing EOS implementations can still throw their documented domain errors; these propagate to the caller. + */ + template + [[nodiscard]] std::expected, SurfaceFailure> evaluate_surface_response(const typename Formulation::EquationOfState& equation_of_state, const dimensions::PressureValue pressure, const SurfacePoint& point) + { + if constexpr (!Formulation::uses_eos) + { + // A pressure unknown uses the prescribed pressure itself as its target, without an EOS relation. + return SurfaceEOSResponse{.value = pressure}; + } + else + { + if (!finite_state(point, typename Formulation::EOSStateQuantities{})) + { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_state, .message = "The surface EOS inputs must be finite; supply an admissible boundary state."}); + } + + const SurfaceEOSResponse response = evaluate_relation(equation_of_state, pressure, point, typename Formulation::Relation::InputQuantities{}, typename Formulation::EOSStateQuantities{}); + bool finite = std::isfinite(response.value.value()); + for (const double derivative : response.field_derivatives.values) + { + finite = finite && std::isfinite(derivative); + } + if (!finite) + { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_eos_response, .message = "The surface EOS returned a non-finite value or partial derivative; check its thermodynamic domain."}); + } + return response; + } + } + + /** @brief Form the declared residual from an already evaluated EOS response, keeping the carrier identity explicit. */ + template + [[nodiscard]] SurfaceLinearization linearize_surface(const dimensions::PressureValue pressure, const SurfacePoint& point, const SurfaceEOSResponse& response) + { + SurfaceLinearization result{}; + result.field_derivatives = response.field_derivatives; + + if constexpr (Formulation::resolves_carrier) + { + result.residual = state_value(point).value() - response.value.value(); + for (double& derivative : result.field_derivatives.values) + { + derivative = -derivative; + } + result.field_derivatives.template get() += 1.0; + } + else + { + result.residual = response.value.value() - pressure.value(); + } + return result; + } + + /** @brief Catch overflow in residual subtraction or derivative accumulation before publishing a linearization. */ + template + [[nodiscard]] std::expected, SurfaceFailure> checked_linearization(const dimensions::PressureValue pressure, const SurfacePoint& point, const SurfaceEOSResponse& response) + { + const SurfaceLinearization result = linearize_surface(pressure, point, response); + bool finite = std::isfinite(result.residual); + for (const double derivative : result.field_derivatives.values) + { + finite = finite && std::isfinite(derivative); + } + if (!finite) + { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_linearization, .message = "The surface residual or derivative is non-finite; reduce the trial state before using this linearization."}); + } + return result; + } +} diff --git a/src/include/serif/geometry/surface/detail/relations.hpp b/src/include/serif/geometry/surface/detail/relations.hpp new file mode 100644 index 00000000..f378d7d1 --- /dev/null +++ b/src/include/serif/geometry/surface/detail/relations.hpp @@ -0,0 +1,97 @@ +#pragma once + +#include +#include + +#include "serif/discretization/forms/value.hpp" +#include "serif/eos/models/concepts.hpp" +#include "serif/utils/dependencies/key_list_operations.hpp" + +namespace serif::geometry::surface { + /** @brief Scalar thermodynamic data may come from any declared scalar space; boundary enforcement is a separate backend capability. */ + template + concept SurfaceScalarTerm = discretization::blocks::Term && requires { + typename discretization::forms::OperandValue>; + requires std::same_as>, double>; + }; + + namespace detail { + template + struct AreSurfaceBindings : std::false_type { + }; + + template + struct AreSurfaceBindings> : std::bool_constant< + sizeof...(Entries) != 0 && (dimensions::ThermodynamicQuantityType && ...) && + (SurfaceScalarTerm && ...)> { + }; + } + + /** @brief Reuse dependency bindings, retaining distinct quantity roles even when they refer to one term. */ + template + concept SurfaceStateBindings = utils::dependencies::IsBindings && detail::AreSurfaceBindings>::value; + + /** @brief The quantity enforcing the boundary must have an actual equation row, not just sampled values. */ + template + concept SurfaceResidualRow = utils::dependencies::BindingsWithBinding && + discretization::blocks::HasResidualBlock>; + + namespace detail { + template + struct RelationQuantities; + + template + struct RelationQuantities> { + using AllInputs = utils::dependencies::KeyList; + using NonPressureInputs = utils::dependencies::CombinedKeyList< + std::conditional_t, utils::dependencies::KeyList<>, utils::dependencies::KeyList>...>; + }; + + /** @brief Check only the derivatives of unknown state; prescribed pressure has no Newton column. */ + template + [[nodiscard]] consteval bool supports_partials(utils::dependencies::KeyList) { + return (eos::models::EOSSupportingPartialDerivative && ...); + } + + /** @brief A resolved carrier must be supplied by one pressure-input relation, with all remaining inputs and partials available. */ + template + [[nodiscard]] consteval bool matches_pressure_relation() { + using Output = typename Relation::OutputQuantity; + using Inputs = typename RelationQuantities::NonPressureInputs; + + return !std::same_as && + (std::same_as || std::same_as) && + eos::relations::RelationContainingInput && + SurfaceResidualRow && + utils::dependencies::SatisfiesKeyList && supports_partials(Inputs{}); + } + + template + struct MatchingPressureRelations; + + template + struct MatchingPressureRelations> { + using Matches = decltype(std::tuple_cat(std::declval(), std::tuple, std::tuple<>>>()...)); + static constexpr std::size_t count = std::tuple_size_v; + }; + + template + struct SurfaceOperands; + + template + struct SurfaceOperands> { + using Type = utils::types::TypeList...>; + }; + + template + struct UniqueTerms; + + template + struct UniqueTerms> { + using Type = utils::dependencies::CombinedKeyList...>; + }; + + template + using BoundTerms = typename UniqueTerms>::Type; + } +} diff --git a/src/include/serif/geometry/surface/evaluation.hpp b/src/include/serif/geometry/surface/evaluation.hpp new file mode 100644 index 00000000..805b6058 --- /dev/null +++ b/src/include/serif/geometry/surface/evaluation.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include "serif/geometry/surface/detail/evaluation.hpp" + +namespace serif::geometry::surface { + /** + * @brief Evaluate the chosen isobaric residual and its analytical derivatives at one boundary point. + * @param equation_of_state The EOS declared by Formulation. + * @param condition The prescribed surface pressure, independent of the choice of carrier. + * @param point Values of the formulation's bound scalar terms, already sampled at a common boundary point. + * @return A conventional equation linearization, or a numerical rejection. No preparation or borrowed state is retained. + * @note Domain exceptions from an existing throwing EOS propagate; no exception is used for surface validation. + */ + template + [[nodiscard]] std::expected, SurfaceFailure> evaluate_pressure_surface(const typename Formulation::EquationOfState& equation_of_state, const IsobaricSurface& condition, const SurfacePoint& point) { + if (!detail::finite_state(point, typename Formulation::StateQuantities{})) { + return std::unexpected(SurfaceFailure{.reason = SurfaceFailureReason::non_finite_state, .message = "The surface state must be finite; supply finite values for every required thermodynamic quantity."}); + } + + const auto response = detail::evaluate_surface_response(equation_of_state, condition.target_pressure(), point); + if (!response) { + return std::unexpected(response.error()); + } + return detail::checked_linearization(condition.target_pressure(), point, *response); + } +} diff --git a/src/include/serif/geometry/surface/isobaric.hpp b/src/include/serif/geometry/surface/isobaric.hpp index e69de29b..8b902a64 100644 --- a/src/include/serif/geometry/surface/isobaric.hpp +++ b/src/include/serif/geometry/surface/isobaric.hpp @@ -0,0 +1,40 @@ +#pragma once + +#include +#include + +#include "serif/dimensions/dimensions.hpp" +#include "serif/dimensions/type_alias.hpp" + +namespace serif::geometry::surface { + struct PressureSurfaceDescriptor final { + double target_pressure; + }; + + class IsobaricSurface final { + public: + using PhysicalQuantity = dimensions::Pressure; + using TargetValue = dimensions::PressureValue; + + public: + struct Parameters final { + dimensions::PressureValue PSurf; + }; + + public: + explicit IsobaricSurface(const Parameters parameters); + explicit IsobaricSurface (const TargetValue targetPressure); + + [[nodiscard]] TargetValue target_pressure() const noexcept; + [[nodiscard]] std::array continuation_identity() const noexcept; + [[nodiscard]] PressureSurfaceDescriptor descriptor() const noexcept; + + private: + TargetValue m_targetPressure; + }; + + template + concept IsobaricSurfaceType = std::same_as, IsobaricSurface>; + + using Isobaric = IsobaricSurface; +} \ No newline at end of file diff --git a/src/include/serif/geometry/surface/linearization.hpp b/src/include/serif/geometry/surface/linearization.hpp new file mode 100644 index 00000000..db6d6c8c --- /dev/null +++ b/src/include/serif/geometry/surface/linearization.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include +#include +#include + +#include "serif/equations/point.hpp" +#include "serif/geometry/surface/compilation.hpp" + +namespace serif::geometry::surface { + enum class SurfaceFailureReason { + non_finite_state, + non_finite_eos_response, + non_finite_linearization, + inconsistent_sample_count, + empty_surface + }; + + /** @brief A rejected thermodynamic candidate, with a sample index when evaluating a boundary batch. */ + struct SurfaceFailure { + SurfaceFailureReason reason; + std::string_view message; + std::size_t sample = std::numeric_limits::max(); + }; + + template + using SurfacePoint = equations::EquationPoint; + + template + using SurfaceLinearization = equations::PointLinearization; + + /** + * @brief An EOS value and all of its local derivatives, evaluated together at one state. + * @note Derivatives sharing a term are accumulated; the later residual calculation supplies the sign and carrier identity. + */ + template + struct SurfaceEOSResponse { + dimensions::QuantityValue value{0.0}; + equations::OperandValues field_derivatives{}; + }; +} diff --git a/src/include/serif/geometry/surface/samples.hpp b/src/include/serif/geometry/surface/samples.hpp new file mode 100644 index 00000000..fc1ab3ee --- /dev/null +++ b/src/include/serif/geometry/surface/samples.hpp @@ -0,0 +1,70 @@ +#pragma once + +#include +#include + +#include "serif/geometry/surface/linearization.hpp" +#include "serif/utils/dependencies/resource.hpp" + +namespace serif::geometry::surface { + namespace keys { + struct EquationOfState { + }; + struct PressureCondition { + }; + struct EOSResponse { + }; + struct Linearization { + }; + } + + /** + * @brief Values of one scalar term at an ordered collection of boundary points. + * @note Every term in a preparation must use the same point ordering. Sampling and finite-element trace choices belong to the backend. + */ + template + struct BoundaryValues { + std::vector values; + + /// Compare the actual sampled values, including their count, without a lossy checksum. + [[nodiscard]] bool operator==(const BoundaryValues&) const = default; + }; + + namespace detail { + struct BoundaryValueTracking : utils::dependencies::ValueSnapshot { + /// Numerical samples have value semantics; rebuilding an equal context need not reevaluate the physics. + template + [[nodiscard]] const void* identity(const BoundaryValues&) const { + return nullptr; + } + }; + + struct PressureConditionTracking { + /// Only prescribed pressure can affect the isobaric calculation. + [[nodiscard]] dimensions::PressureValue capture(const IsobaricSurface& condition) const { + return condition.target_pressure(); + } + + /// Compare pressure exactly, retaining the dependency engine's normal change semantics. + [[nodiscard]] bool matches(const dimensions::PressureValue captured, const IsobaricSurface& condition) const { + return captured == condition.target_pressure(); + } + + /// An equivalent pressure condition may be supplied by a newly constructed context. + [[nodiscard]] const void* identity(const IsobaricSurface&) const { + return nullptr; + } + }; + } + + /** @brief Supply owned boundary samples to an ordinary dependency context, keyed by their physical term. */ + template + [[nodiscard]] auto boundary_values(std::vector values) { + return utils::dependencies::own(BoundaryValues{.values = std::move(values)}, detail::BoundaryValueTracking{}); + } + + /** @brief Supply the isobaric condition as an owned dependency source, with pressure tracked by value. */ + [[nodiscard]] inline auto pressure_condition(IsobaricSurface condition) { + return utils::dependencies::own(std::move(condition), detail::PressureConditionTracking{}); + } +} diff --git a/src/include/serif/geometry/surface/surface.hpp b/src/include/serif/geometry/surface/surface.hpp new file mode 100644 index 00000000..ebb1b5e5 --- /dev/null +++ b/src/include/serif/geometry/surface/surface.hpp @@ -0,0 +1,17 @@ +#pragma once + +/** + * @file surface.hpp + * @brief Isobaric conditions, EOS-selected formulations, and reusable boundary calculations. + * + * IsobaricSurface describes the prescribed pressure. PressureSurfaceFormulation selects one compatible EOS relation and its carrier; + * DirectPressureSurfaceFormulation explicitly chooses a forward pressure residual and its equation row. Neither selects a finite-element trace. + * + * Use evaluate_pressure_surface for one sampled point, or make_pressure_surface_graph for an ordered batch. The latter returns an ordinary + * dependency graph: the caller supplies the EOS with its tracking policy, the pressure condition, and typed boundary values. Results use the + * same PointLinearization representation as other equations. Finite-element row enforcement remains a separate numerical capability. + */ + +#include "serif/geometry/surface/compilation.hpp" +#include "serif/geometry/surface/computations.hpp" +#include "serif/geometry/surface/evaluation.hpp" diff --git a/src/lib/serif/geometry/surface/isobaric.cpp b/src/lib/serif/geometry/surface/isobaric.cpp index e69de29b..fed8ac4a 100644 --- a/src/lib/serif/geometry/surface/isobaric.cpp +++ b/src/lib/serif/geometry/surface/isobaric.cpp @@ -0,0 +1,31 @@ +#include +#include +#include + +#include "serif/geometry/surface/isobaric.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::geometry::surface { + using namespace utils::errors; + + IsobaricSurface::IsobaricSurface(const Parameters parameters) : IsobaricSurface(parameters.PSurf) {} + + IsobaricSurface::IsobaricSurface(const TargetValue targetPressure) : m_targetPressure(targetPressure) { + serif_ensure(std::isfinite(targetPressure.value()), "The target surface pressure must be finite, got {}", targetPressure.value()); + serif_ensure(targetPressure.value() >= 0, "The target surface pressure must be non-negative, got {}", targetPressure.value()); + } + + IsobaricSurface::TargetValue IsobaricSurface::target_pressure() const noexcept { + return m_targetPressure; + } + + std::array IsobaricSurface::continuation_identity() const noexcept { + return {m_targetPressure.value()}; + } + + PressureSurfaceDescriptor IsobaricSurface::descriptor() const noexcept { + return PressureSurfaceDescriptor{ + .target_pressure = m_targetPressure.value() + }; + } +}