diff --git a/src/include/serif/model/configured_discretization.hpp b/src/include/serif/model/configured_discretization.hpp new file mode 100644 index 00000000..0a23c097 --- /dev/null +++ b/src/include/serif/model/configured_discretization.hpp @@ -0,0 +1,35 @@ +#pragma once + +#include + +#include "serif/model/normalization/reference.hpp" + +namespace serif::model::preconditioning { + /** @brief Evaluate physical equations without preparing an inverse until a recipe is selected. */ + struct None final {}; +} + +namespace serif::model { + /** @brief Select a descriptive star, its existing finite-element resources, and its numerical coordinates. */ + template , typename Preconditioner = preconditioning::None> + struct DiscretizationOptions final { + Model model; + FiniteElements finite_elements; + Normalization normalization; + Preconditioner preconditioner; + }; + + /** + * @brief Resolve model defaults while preserving compile-time field, scheme, and coordinate choices. + * @note Finite-element resources remain explicit: normalization does not choose a mesh or approximation order. + */ + template + requires normalization::ReferenceResolvable + [[nodiscard]] auto make_discretization(DiscretizationOptions options) { + auto reference = normalization::resolve_configuration(std::move(options.normalization), options.model); + return DiscretizationOptions{ + .model = std::move(options.model), .finite_elements = std::move(options.finite_elements), .normalization = std::move(reference), + .preconditioner = std::move(options.preconditioner) + }; + } +} diff --git a/src/include/serif/model/description.hpp b/src/include/serif/model/description.hpp new file mode 100644 index 00000000..d6c80b42 --- /dev/null +++ b/src/include/serif/model/description.hpp @@ -0,0 +1,7 @@ +#pragma once + +#include "serif/model/equilibrium_equations.hpp" +#include "serif/model/specification.hpp" +#include "serif/model/specifications/fixed_angular_momentum.hpp" +#include "serif/model/specifications/fixed_total_mass.hpp" +#include "serif/model/stellar_model.hpp" diff --git a/src/include/serif/model/detail/specification_endpoints.hpp b/src/include/serif/model/detail/specification_endpoints.hpp new file mode 100644 index 00000000..ae0c6401 --- /dev/null +++ b/src/include/serif/model/detail/specification_endpoints.hpp @@ -0,0 +1,111 @@ +#pragma once + +#include "serif/model/specification.hpp" +#include "serif/model/specifications/physical_names.hpp" +#include "serif/utils/dependencies/binding.hpp" +#include "serif/utils/dependencies/key_list_operations.hpp" +#include "serif/utils/misc/concepts/filter.hpp" + +namespace serif::model::detail +{ + template + using GeneratesState = std::bool_constant; + + /** @brief Numerical checks operate only on prescriptions that introduce a coordinate and equation. */ + template + using GeneratedSpecifications = utils::misc::concepts::FilterT; + + /** @brief Recognize generated endpoint identities without interpreting unrelated physical names. */ + template + struct EndpointOwner { using Type = void; }; + template + struct EndpointOwner> { using Type = Specification; }; + template + struct EndpointOwner> { using Type = Specification; }; + + /** @brief Resolve owned aliases and explicit bindings without merging same-unit identities. */ + template + [[nodiscard]] consteval auto resolved_endpoint() + { + using namespace utils::dependencies; + if constexpr (std::same_as) + { + return resolved_endpoint, Bindings>(); + } + else if constexpr (std::same_as) + { + return resolved_endpoint, Bindings>(); + } + else if constexpr (!std::same_as::Type, void> && + BindingsWithBinding::Type>>) + { + using Specification = typename EndpointOwner::Type; + using Enforcement = bound_type_t>; + if constexpr (std::same_as>) + { + return std::type_identity{}; + } + else + { + return std::type_identity{}; + } + } + else if constexpr (BindingsWithBinding) + { + return std::type_identity>{}; + } + else + { + return std::type_identity{}; + } + } + + template + using ResolvedEndpoint = typename decltype(resolved_endpoint())::type; + + /** @brief A coordinate binding cannot supply an absent or non-generating physical prescription. */ + template + [[nodiscard]] consteval bool endpoint_owner_present() + { + if constexpr (requires { typename Endpoint::SpecificationType; }) + { + return Model::template contains_specification && DimensionalScalarSpecification; + } + return true; + } + + template + struct ResolvedEndpoints; + + template + struct ResolvedEndpoints, Bindings> + { + using Keys = utils::dependencies::CombinedKeyList>...>; + + template + static constexpr bool owners_present = (endpoint_owner_present() && ...); + }; + + /** @brief Legacy prescriptions retain their original dependencies unless an equation explicitly registers another enforcement. */ + template >> + struct EnforcementDependencies + { + using Type = typename Specification::ModelDefinition; + }; + template + struct EnforcementDependencies + { + using Type = utils::dependencies::bound_type_t>; + }; + + /** @brief Share resolved physical identities between endpoint and numerical-provider checks. */ + template + struct ResolvedSpecification + { + using Coordinate = ResolvedEndpoint, Bindings>; + using Constraint = ResolvedEndpoint, Bindings>; + using Enforcement = typename EnforcementDependencies::Type; + using Reads = ResolvedEndpoints; + using Changes = ResolvedEndpoints; + }; +} diff --git a/src/include/serif/model/detail/specification_list.hpp b/src/include/serif/model/detail/specification_list.hpp new file mode 100644 index 00000000..38bde844 --- /dev/null +++ b/src/include/serif/model/detail/specification_list.hpp @@ -0,0 +1,97 @@ +#pragma once + +#include +#include + +#include "serif/model/specification.hpp" +#include "serif/utils/misc/concepts/filter.hpp" + +namespace serif::model::detail +{ + /** @brief Compare physical identities in the historical role, coordinate-kind, stable-name order. */ + template + [[nodiscard]] consteval bool specification_precedes() + { + using L = typename Left::ModelDefinition; + using R = typename Right::ModelDefinition; + return std::tuple{L::role, L::generated_state_kind, L::name} < std::tuple{R::role, R::generated_state_kind, R::name}; + } + + template + struct InsertSpecification; + + template + struct InsertSpecification> + { + using Type = utils::types::TypeList; + }; + + template + struct InsertSpecification> + { + using Type = std::conditional_t(), utils::types::TypeList, + utils::misc::concepts::PrependT>::Type>>; + }; + + template + struct CanonicalSpecifications; + + template <> + struct CanonicalSpecifications<> + { + using Type = utils::types::TypeList<>; + }; + + template + struct CanonicalSpecifications + { + using Type = typename InsertSpecification::Type>::Type; + }; + + /** @brief Reject semantic collisions even when two prescriptions have different C++ types or scalar kinds. */ + template + [[nodiscard]] consteval bool unique_specification_names() + { + constexpr std::array descriptions{{ + {Specifications::ModelDefinition::name, Specifications::ModelDefinition::role, Specifications::ModelDefinition::generated_state_kind}... + }}; + for (std::size_t left = 0; left < descriptions.size(); ++left) + { + for (std::size_t right = left + 1; right < descriptions.size(); ++right) + { + if (descriptions[left].role == descriptions[right].role && descriptions[left].name == descriptions[right].name) + { + return false; + } + } + } + return true; + } + + template + struct SelectRole + { + template + using HasRole = std::bool_constant; + using Type = utils::misc::concepts::FilterT, HasRole>; + }; + + template + struct ArgumentsMatch; + + template + struct ArgumentsMatch, Arguments...> + { + template + static constexpr std::size_t count = (std::size_t{0} + ... + std::size_t{std::same_as>}); + static constexpr bool value = sizeof...(Specifications) == sizeof...(Arguments) && ((count == 1) && ...); + }; + + /** @brief Consume the argument tuple once, moving each specification into its canonical slot once. */ + template + [[nodiscard]] std::tuple canonical_specifications(utils::types::TypeList, Arguments&&... arguments) + { + std::tuple...> values{std::forward(arguments)...}; + return {std::get(std::move(values))...}; + } +} // namespace serif::model::detail diff --git a/src/include/serif/model/detail/surface_plan.hpp b/src/include/serif/model/detail/surface_plan.hpp new file mode 100644 index 00000000..4abfb3c5 --- /dev/null +++ b/src/include/serif/model/detail/surface_plan.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include "serif/geometry/surface/automatic.hpp" +#include "serif/model/equilibrium_equations.hpp" + +namespace serif::model::detail +{ + template + struct SurfaceFieldHint + { + using Type = utils::dependencies::KeyList<>; + }; + + /** @brief An equation author may disambiguate surface enforcement without changing the physical condition. */ + template + struct SurfaceFieldHint> + { + using Type = utils::dependencies::KeyList; + }; + + template + struct UniqueSurfaceField; + + template <> + struct UniqueSurfaceField> + { + using Type = void; + }; + + template + struct UniqueSurfaceField> + { + static_assert(sizeof...(Others) == 0, "The selected equations name conflicting IsobaricSurfaceField types. Choose one physical field for the isobaric condition."); + using Type = Field; + }; + + template + struct PreferredSurfaceField; + + template + struct PreferredSurfaceField> : UniqueSurfaceField::Type...>> + { + }; + + /** @brief A disabled or unavailable boundary has no compiled field or row requirements. */ + template + struct SelectedSurface + { + using Formulation = void; + using Coefficients = utils::dependencies::KeyList<>; + using Rows = utils::dependencies::KeyList<>; + static constexpr bool supported = Supported; + }; + + template + struct SelectedSurface + { + using Formulation = typename Selection::Formulation; + using Coefficients = typename Formulation::StateTerms; + using Rows = utils::dependencies::KeyList; + static constexpr bool supported = true; + }; + + template + struct SurfacePlan; + + template + struct SurfacePlan : SelectedSurface + { + static constexpr bool supported = true; + }; + + /** @brief Compile the condition once; planning, graph assembly and validation share the selected formulation. */ + template + struct SurfacePlan : SelectedSurface< + geometry::surface::AutomaticPressureSurface::Type>, + geometry::surface::AutomaticPressureSurfaceCompilable::Type>> + { + }; +} diff --git a/src/include/serif/model/eos_reference.hpp b/src/include/serif/model/eos_reference.hpp new file mode 100644 index 00000000..8665001c --- /dev/null +++ b/src/include/serif/model/eos_reference.hpp @@ -0,0 +1,96 @@ +#pragma once + +#include +#include +#include +#include + +#include "serif/utils/error/errors.hpp" + +namespace serif::model +{ + /** @brief Reuse one immutable EOS in equation objects, optionally retaining shared ownership. */ + template + requires requires { typename EOS::Relations; } + class EosReference final + { + public: + using Relations = typename EOS::Relations; + + /** @brief Retain an owned EOS, including an aliasing handle into an immutable model. */ + explicit EosReference(std::shared_ptr equation_of_state) : m_owner(std::move(equation_of_state)), m_eos(m_owner.get()) + { + using namespace utils::errors; + serif_ensure(m_eos != nullptr, "An EOS reference requires an equation of state; supply a non-null shared owner."); + } + + /** @brief Borrow an EOS whose model owner must outlive all equations and prepared actions using it. */ + explicit EosReference(const EOS& equation_of_state) noexcept : m_eos(&equation_of_state) + { + } + + /** @brief Never retain a pointer into a temporary EOS. */ + EosReference(EOS&&) = delete; + /** @brief Const temporaries have the same borrowing lifetime restriction. */ + EosReference(const EOS&&) = delete; + + /** @brief Forward only the relations actually implemented by the underlying EOS. */ + template + requires requires(const EOS& eos, Relation relation, Values... values) { eos.evaluate(relation, values...); } + [[nodiscard]] decltype(auto) evaluate(Relation relation, Values... values) const noexcept(noexcept(m_eos->evaluate(relation, values...))) + { + return m_eos->evaluate(relation, values...); + } + + /** @brief Preserve exact derivative availability and physical return types. */ + template + requires requires(const EOS& eos, Relation relation, Input input, Values... values) { eos.partial_derivative(relation, input, values...); } + [[nodiscard]] decltype(auto) partial_derivative(Relation relation, Input input, Values... values) const noexcept(noexcept(m_eos->partial_derivative(relation, input, values...))) + { + return m_eos->partial_derivative(relation, input, values...); + } + + /** @brief Forward second partials only when the underlying EOS provides them. */ + template + requires requires(const EOS& eos, Relation relation, FirstInput first, SecondInput second, Values... values) { + eos.second_partial_derivative(relation, first, second, values...); + } + [[nodiscard]] decltype(auto) second_partial_derivative(Relation relation, FirstInput first, SecondInput second, Values... values) const + noexcept(noexcept(m_eos->second_partial_derivative(relation, first, second, values...))) { + return m_eos->second_partial_derivative(relation, first, second, values...); + } + + /** @brief Inspect EOS-specific physical parameters without storing another copy. */ + [[nodiscard]] const EOS& get() const noexcept + { + return *m_eos; + } + + private: + std::shared_ptr m_owner; + const EOS* m_eos; + }; + + /** @brief Share the EOS subobject while retaining the model that owns it. */ + template + [[nodiscard]] auto share_equation_of_state(std::shared_ptr model) + { + using namespace utils::errors; + serif_ensure(model != nullptr, "Sharing an equation of state requires a non-null stellar model owner."); + using EOS = std::remove_cvref_tequation_of_state())>; + const EOS* eos = &model->equation_of_state(); + return EosReference{std::shared_ptr{std::move(model), eos}}; + } + + /** @brief Borrow the model-owned EOS when the prepared graph already retains its immutable model owner. */ + template + [[nodiscard]] EosReference borrow_eos(const EOS& eos) noexcept + { + return EosReference{eos}; + } + + /** @brief Borrowing never extends a temporary EOS lifetime. */ + template + requires (!std::is_lvalue_reference_v) + void borrow_eos(EOS&&) = delete; +} // namespace serif::model diff --git a/src/include/serif/model/equilibrium_equations.hpp b/src/include/serif/model/equilibrium_equations.hpp new file mode 100644 index 00000000..3c9f57ca --- /dev/null +++ b/src/include/serif/model/equilibrium_equations.hpp @@ -0,0 +1,101 @@ +#pragma once + +#include +#include +#include +#include + +namespace serif::model +{ + /** @brief Explicitly omit a shape equation from a declarative formulation. */ + struct NoShapeForce + { + }; + + /** @brief Explicitly omit a surface prescription from a declarative formulation. */ + struct NoSurfaceCondition + { + }; + + /** @brief State the chosen equations explicitly, without selecting physics from the EOS. */ + template + struct EquilibriumEquations + { + using InteriorEquations = std::remove_cvref_t; + using ShapeForce = std::remove_cvref_t; + using SurfaceCondition = std::remove_cvref_t; + + Interior interior_equations; + Shape shape_force{}; + Surface surface_condition{}; + }; + + /** + * @brief Preserve every authored member when a declaration passes through a templated model factory. + * @note Explicit deduction avoids compiler substitution dropping defaulted aggregate members. Deduced + * declarations supply preceding slots: use NoShapeForce before a surface condition when the shape force is + * absent. Explicitly typed aggregates may still skip members. + */ + template + EquilibriumEquations(Interior, Shape = {}, Surface = {}) -> EquilibriumEquations; + + namespace detail + { + template + inline constexpr bool is_equilibrium_equations = false; + + template + inline constexpr bool is_equilibrium_equations> = true; + + /** @brief Evaluate a physical definition or retain its immutable object without requiring a copy. */ + template + [[nodiscard]] decltype(auto) resolve_equation_member(Definition&& definition, const Model& model) + { + if constexpr (std::invocable) + { + return std::invoke(std::forward(definition), model); + } + else if constexpr (std::is_lvalue_reference_v) + { + return std::forward(definition); + } + else + { + return std::remove_cvref_t{std::forward(definition)}; + } + } + + /** @brief Resolve each explicit part, owning temporary results and borrowing unchanged model members. */ + template + [[nodiscard]] auto resolve_equation_members(Definition&& definition, const Model& model) + { + using Interior = decltype(resolve_equation_member(std::forward(definition).interior_equations, model)); + using Shape = decltype(resolve_equation_member(std::forward(definition).shape_force, model)); + using Surface = decltype(resolve_equation_member(std::forward(definition).surface_condition, model)); + return EquilibriumEquations + { + .interior_equations = resolve_equation_member(std::forward(definition).interior_equations, model), + .shape_force = resolve_equation_member(std::forward(definition).shape_force, model), + .surface_condition = resolve_equation_member(std::forward(definition).surface_condition, model) + }; + } + } // namespace detail + + /** @brief Recognize composition without instantiating any numerical equation contract. */ + template + concept EquilibriumEquationDescription = detail::is_equilibrium_equations>; + + /** @brief Resolve an explicit composition at binding; retain the model while the result borrows its members. */ + template + requires requires(const Model& model) { model.equations(); } && + EquilibriumEquationDescription().equations(), std::declval()))> + [[nodiscard]] auto resolve_equilibrium_equations(const Model& model) + { + return detail::resolve_equation_members(detail::resolve_equation_member(model.equations(), model), model); + } + + /** @brief A temporary model cannot retain the equation members borrowed by the resolved description. */ + template + requires (!std::is_lvalue_reference_v) + void resolve_equilibrium_equations(Model&&) = delete; +} // namespace serif::model diff --git a/src/include/serif/model/equilibrium_plan.hpp b/src/include/serif/model/equilibrium_plan.hpp new file mode 100644 index 00000000..11fcc19a --- /dev/null +++ b/src/include/serif/model/equilibrium_plan.hpp @@ -0,0 +1,66 @@ +#pragma once + +#include "serif/model/equilibrium_equations.hpp" +#include "serif/model/plan.hpp" +#include "serif/model/detail/surface_plan.hpp" +#include "serif/operators/concepts.hpp" +#include "serif/equations/surface_potential.hpp" +#include "serif/equations/primal_gravity.hpp" + +namespace serif::model +{ + namespace detail + { + template + struct ShapePlan; + + template <> + struct ShapePlan + { + using Coefficients = utils::dependencies::KeyList<>; + using Values = utils::dependencies::KeyList<>; + using Rows = utils::dependencies::KeyList<>; + }; + + /** @brief The isobar level couples only existing scalar potential, rotation and mass coordinates. */ + template <> + struct ShapePlan + { + using Coefficients = typename OperandTerms::Type; + using Values = utils::dependencies::KeyList; + using Rows = utils::dependencies::KeyList; + }; + + } + + /** + * @brief Collect the explicitly selected physical coefficients and supported constraint coordinates. + * @note This plan verifies representation, not mathematical closure. Equal units never merge distinct + * physical identities, and prescribed coefficients remain available as graph inputs without becoming unknowns. + */ + template , + utils::dependencies::IsKeyList AdditionalFields = utils::dependencies::KeyList<>> + struct EquilibriumPlan final + { + static constexpr bool has_shape = !std::same_as; + static constexpr bool has_surface = !std::same_as; + static constexpr bool has_primal_gravity = [](utils::types::TypeList) { + return (std::same_as || ...); + }(typename Description::InteriorEquations::Equations{}); + using Interior = detail::InteriorPlan; + using Shape = detail::ShapePlan; + using ResidualEquations = typename Interior::EquationsList; + using AvailableFields = utils::dependencies::CombinedKeyList; + using SolvedFields = typename detail::RemoveKeys::Type; + using Residuals = utils::dependencies::CombinedKeyList; + using Surface = detail::SurfacePlan; + using Coefficients = utils::dependencies::CombinedKeyList; + static_assert(!has_shape || !utils::dependencies::ContainsKey, + "Prescribe geometry by selecting NoShapeForce and supplying a displacement field. An active shape-force equation uses independent solved surface coordinates."); + using Values = utils::dependencies::CombinedKeyList::Type, typename Shape::Values>; + // Geometry is supplied separately even when no equation solves for the shape. + using PhysicalInputs = utils::dependencies::CombinedKeyList>; + using InactiveRows = std::conditional_t, utils::dependencies::KeyList>; + using Layout = discretization::blocks::SystemLayout; + }; +} diff --git a/src/include/serif/model/execution_validation.hpp b/src/include/serif/model/execution_validation.hpp new file mode 100644 index 00000000..71124362 --- /dev/null +++ b/src/include/serif/model/execution_validation.hpp @@ -0,0 +1,144 @@ +#pragma once + +#include "serif/model/equilibrium_plan.hpp" +#include "serif/model/validation.hpp" + +namespace serif::model +{ + namespace detail + { + template + struct CartesianCouplings; + + template + struct CartesianCouplings, utils::dependencies::KeyList> + { + template + using ForRow = utils::dependencies::KeyList...>; + using Type = utils::dependencies::CombinedKeyList...>; + }; + + template + struct ExecutionSurface : SelectedSurface + { + }; + + template + struct ExecutionSurface> : Plan::Surface + { + }; + + template + struct ExecutionInputs + { + using PhysicalInputs = typename Plan::Values; + using InactiveRows = utils::dependencies::KeyList<>; + }; + + template + struct ExecutionInputs> + { + using PhysicalInputs = typename Plan::PhysicalInputs; + using InactiveRows = typename Plan::InactiveRows; + }; + + /** @brief Demand an explicitly supplied generated equation and effect before inferring closed-row zeros. */ + template + [[nodiscard]] consteval bool generated_physics_is_present() + { + using namespace utils::dependencies; + using Endpoints = ResolvedSpecification; + return ContainsKey && [](KeyList) + { + return (((ContainsKey && !ContainsKey) || ContainsKey, Actual>) && ...); + }(typename Endpoints::Changes::Keys{}); + } + + /** @brief Check only numerical-provider admission; physical endpoint ownership is checked separately. */ + template + [[nodiscard]] consteval bool generated_physics_are_present(utils::types::TypeList) + { + return (generated_physics_is_present() && ...); + } + + /** @brief Every surface carrier must have an equation row in the selected layout. */ + template + [[nodiscard]] consteval bool equation_rows_exist(utils::dependencies::KeyList) + { + return (utils::dependencies::ContainsKey && ...); + } + } // namespace detail + + /** @brief Recognize current equation backends without instantiating an invalid complete plan. */ + template + concept ExecutableEquilibriumDescription = EquilibriumEquationDescription && ExecutableInteriorEquations && requires + { + typename detail::ShapePlan::Rows; + requires std::same_as || std::same_as; + requires std::same_as || + std::same_as; + }; + + namespace detail + { + template + struct EquilibriumCouplings + { + using Interior = InteriorPlan; + using Shape = ShapePlan; + using GeometryRows = std::conditional_t, utils::dependencies::KeyList<>, + utils::dependencies::CombinedKeyList>; + using Geometry = typename CartesianCouplings>::Type; + using Force = typename CartesianCouplings::Type; + using Type = utils::dependencies::CombinedKeyList, Force, Geometry, + typename CartesianCouplings::Type>; + using AuthoredRows = utils::dependencies::CombinedKeyList; + using RequiredRows = utils::dependencies::CombinedKeyList; + }; + } // namespace detail + + /** @brief Inspect the actual derivative producers admitted by the explicitly supplied physical formulation. */ + template + using ActualEquilibriumCouplings = typename detail::EquilibriumCouplings>::Type; + + /** + * @brief Validate a complete selected formulation without inferring solvability or selecting physics. + * @note Closed equations define zero derivatives for other coordinates, but that fact cannot admit a + * missing custom multiplier force. ExplicitZeroCouplings is an explicit author registration, not a fallback. + */ + template > + [[nodiscard]] consteval bool executable_model_topology() + { + using namespace utils::dependencies; + if constexpr (!StellarModelType || !ExecutableEquilibriumDescription || !IsBindings || !IsKeyList) + { + return false; + } + else if constexpr (!requires { typename Plan::Values; typename Plan::Residuals; }) + { + return false; + } + else if constexpr (!IsKeyList || !IsKeyList) + { + return false; + } + else + { + using Surface = detail::ExecutionSurface; + using Inputs = detail::ExecutionInputs; + using Couplings = detail::EquilibriumCouplings; + if constexpr (!SpecificationEndpointsValid) + { + return false; + } + else + { + using Actual = CombinedKeyList; + return (std::same_as || Surface::supported) && + detail::equation_rows_exist(typename Couplings::RequiredRows{}) && + detail::generated_physics_are_present(detail::GeneratedSpecifications>{}); + } + } + } +} // namespace serif::model diff --git a/src/include/serif/model/finite_elements.hpp b/src/include/serif/model/finite_elements.hpp new file mode 100644 index 00000000..4d78d3ed --- /dev/null +++ b/src/include/serif/model/finite_elements.hpp @@ -0,0 +1,13 @@ +#pragma once + +#include "serif/backend/mfem/model/finite_elements.hpp" + +namespace serif::model { + // Public construction names retain the backend's native resource ownership without forwarding calls. + using backend::mfem::make_stellar_discretization; + using backend::mfem::with_deformation; +} + +namespace serif::geometry::deformation { + using backend::mfem::deformation_discretization; +} diff --git a/src/include/serif/model/keys.hpp b/src/include/serif/model/keys.hpp new file mode 100644 index 00000000..e3fe140e --- /dev/null +++ b/src/include/serif/model/keys.hpp @@ -0,0 +1,9 @@ +#pragma once + +namespace serif::model::keys +{ + /** @brief The complete accepted physical residual and its matching Jacobian. */ + struct PhysicalSystem + { + }; +} diff --git a/src/include/serif/model/mesh_options.hpp b/src/include/serif/model/mesh_options.hpp new file mode 100644 index 00000000..248b647f --- /dev/null +++ b/src/include/serif/model/mesh_options.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include +#include "serif/dimensions/dimensions.hpp" + +namespace serif::model { + /** @brief Add mesh resolution at the stellar interface before creating field spaces. */ + struct StellarSurfaceRefinementOptions { + /// Additional isotropic levels on both sides of the entire stellar surface. + int levels{0}; + /// Additional latitude splits in closed rings around the rotation axis. + int equatorial_levels{0}; + /// Minimum half-width of the resolved equatorial band, in degrees; 0 < angle <= 35. + double equatorial_half_angle_degrees{15.0}; + /// Additional radial halvings of the stellar-side surface shell only; angular surface faces are unchanged. + int radial_levels{0}; + /// Align the mesh's polar direction with this finite nonzero axis through the stellar center. + std::array axis{0.0, 0.0, 1.0}; + }; + + /** + * @brief Resolve a spherical stellar interior and a compactified exterior on a curved reference mesh. + * @note Field approximation orders come from the model's term types. Geometry order controls only + * the mesh; its default resolves the coarse spherical boundary for radial seed projection. + */ + struct StellarMeshOptions { + dimensions::LengthValue radius{1.0}; + double core_radius_fraction{0.25}; + double infinity_radius_fraction{4.0}; + int refinement{0}; + int geometry_order{6}; + StellarSurfaceRefinementOptions surface_refinement{}; + }; +} diff --git a/src/include/serif/model/model.hpp b/src/include/serif/model/model.hpp new file mode 100644 index 00000000..a8aade37 --- /dev/null +++ b/src/include/serif/model/model.hpp @@ -0,0 +1,7 @@ +#pragma once + +#include "serif/model/description.hpp" +#include "serif/model/eos_reference.hpp" +#include "serif/model/equilibrium_plan.hpp" +#include "serif/model/validation.hpp" +#include "serif/model/state.hpp" diff --git a/src/include/serif/model/normalization/plan.hpp b/src/include/serif/model/normalization/plan.hpp new file mode 100644 index 00000000..c42c6794 --- /dev/null +++ b/src/include/serif/model/normalization/plan.hpp @@ -0,0 +1,131 @@ +#pragma once + +#include "serif/model/equilibrium_plan.hpp" +#include "serif/solver/normalization/configuration.hpp" +#include "serif/solver/normalization/physical_coordinates.hpp" + +namespace serif::model::normalization { + /** @brief Reduced equation rows whose volume residual is replaced by the selected surface condition. */ + struct SurfaceRows final {}; + /** @brief Potential rows replaced by the zero trace at compactified infinity. */ + struct InfinityRows final {}; + + namespace detail { + namespace numerical = solver::normalization; + namespace dependencies = utils::dependencies; + namespace blocks = discretization::blocks; + + template + using Coordinate = numerical::DimensionlessRieszCoordinate; + + template > + struct Metric { + using Type = numerical::UnitMetric; + }; + + template + struct Metric { + using Type = numerical::ReferenceVolumeMetric; + }; + + template + struct ValueRule { + using Type = dependencies::Binding::Type>>; + }; + + template + struct ValueRule { + using Type = dependencies::Binding>; + }; + + template + struct SurfaceRule { + using Type = Method; + }; + + template + struct SurfaceRule { + using Formulation = typename Plan::Surface::Formulation; + using Type = std::conditional_t, + numerical::WithPointRows, Method>; + }; + + /** @brief Pointwise infinity values have potential units rather than integrated Poisson units. */ + template + struct InfinityRule { + using Type = std::conditional_t, + numerical::WithPointRows, Method>; + }; + + template + struct ResidualRules; + + template + struct ResidualRules, dependencies::Bindings> { + using Type = dependencies::Bindings; + }; + + template + struct AutomaticRules; + + /** @brief Equation residual quantities remain independent from the unknowns occupying their rows. */ + template + struct AutomaticRules, utils::types::TypeList> { + static_assert((requires { typename Equations::ResidualQuantity; } && ...), + "Automatic normalization requires a residual quantity for each equation; otherwise supply an explicit normalization plan."); + using ValueRules = dependencies::Bindings::Type...>; + using InteriorRules = dependencies::Bindings::Type>, typename Equations::Test, Plan>::Type, + typename Equations::Test, Plan>::Type>...>; + using ShapeRules = std::conditional_t>>, dependencies::Bindings<>>; + using Type = numerical::NormalizationPlan::Type>; + }; + + template + struct SelectPlan { + using Type = Selection; + }; + + template + struct SelectPlan : AutomaticRules {}; + + template + struct MethodMetrics { + using Type = dependencies::KeyList<>; + }; + + template + struct MethodMetrics> { + using Type = dependencies::KeyList; + }; + + template + struct MethodMetrics> { + using Type = dependencies::KeyList; + }; + + template + struct MethodMetrics> : MethodMetrics {}; + + template + struct RuleMetrics; + + template + struct RuleMetrics> { + using Type = dependencies::CombinedKeyList::Type...>; + }; + } + + /** @brief Compile automatic physics-aware rules or retain the user's complete explicit plan. */ + template + using PlanFor = typename detail::SelectPlan::Type; + + /** @brief Assemble each metric once even when several state and residual rules share it. */ + template + using MetricsFor = utils::dependencies::CombinedKeyList::Type, + typename detail::RuleMetrics::Type>; +} diff --git a/src/include/serif/model/normalization/reference.hpp b/src/include/serif/model/normalization/reference.hpp new file mode 100644 index 00000000..e7e94b97 --- /dev/null +++ b/src/include/serif/model/normalization/reference.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include + +#include "serif/model/description.hpp" +#include "serif/solver/normalization/configuration.hpp" + +namespace serif::model::normalization { + namespace detail { + /** @brief Explicit reference values take precedence over model defaults without changing physical prescriptions. */ + template + concept ReferenceAvailable = !std::same_as || Model::template contains_specification; + + /** @brief Resolve a typed reference source only when its model default is available. */ + template + requires ReferenceAvailable + [[nodiscard]] auto reference_value(const Source& source, const Model& model) { + if constexpr (std::same_as) { + return model.template specification().target(); + } else { + return source; + } + } + } + + /** @brief Resolve the Lane–Emden reference from explicit inputs or the model's mass target; no model prescribes central density. */ + template + requires detail::ReferenceAvailable && std::same_as + [[nodiscard]] auto resolve_reference(const solver::normalization::LaneEmdenApproximation& scheme, const Model& model) { + return solver::normalization::LaneEmdenApproximation{ + .n = scheme.n, + .mass = detail::reference_value(scheme.mass, model), + .central_density = scheme.central_density, + .integration = scheme.integration + }; + } + + /** @brief Explicit or custom reference schemes need no model-specific defaults. */ + template + [[nodiscard]] Scheme resolve_reference(const Scheme& scheme, const Model&) { + return scheme; + } + + /** @brief Check model-provided reference requirements before instantiating numerical preparation. */ + template + concept ReferenceResolvable = requires(const Options& options, const Model& model) { + resolve_reference(options.scheme, model); + }; + + /** @brief Freeze resolved reference choices while retaining the deduced coordinate and plan types. */ + template + requires ReferenceResolvable + [[nodiscard]] auto resolve_configuration(Options options, const Model& model) { + return solver::normalization::NormalizationOptions{ + .scheme = resolve_reference(options.scheme, model), + .coordinates = options.coordinates, + .plan = options.plan, + .gravitational_constant = options.gravitational_constant, + .sources = std::move(options.sources), + .computations = std::move(options.computations) + }; + } +} diff --git a/src/include/serif/model/plan.hpp b/src/include/serif/model/plan.hpp new file mode 100644 index 00000000..ebf6a245 --- /dev/null +++ b/src/include/serif/model/plan.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include +#include + +#include "serif/discretization/blocks/system_layout.hpp" +#include "serif/equations/registry.hpp" +#include "serif/operators/plan.hpp" + +namespace serif::model +{ + namespace detail + { + template + struct OperandTerms; + + template + struct OperandTerms> + { + using Type = utils::dependencies::CombinedKeyList...>; + }; + + template + struct InteriorPlan; + + template + struct InteriorPlan> + { + using EquationsList = utils::types::TypeList; + using Rows = utils::dependencies::KeyList; + using Coefficients = utils::dependencies::CombinedKeyList::Columns>::Type...>; + }; + + template + struct RemoveKeys; + + template + struct RemoveKeys, Excluded> + { + using Type = utils::dependencies::CombinedKeyList, utils::dependencies::KeyList<>, utils::dependencies::KeyList>...>; + }; + } + + /** @brief The current interior backend can evaluate every explicitly listed equation. */ + template + concept ExecutableInteriorEquations = requires + { + typename detail::InteriorPlan>::Rows; + }; + + /** + * @brief Derive physical coordinates from explicit equations, independently of numerical solvability. + * @note Prescribed coefficients remain equation inputs but are excluded from solved coordinates. + * Additional coordinates and rows describe explicit shape or constraint contributions. Neither + * equation counting nor this plan claims absence of null spaces or existence of an equilibrium. + */ + template , + utils::dependencies::IsKeyList AdditionalValues = utils::dependencies::KeyList<>, + utils::dependencies::IsKeyList AdditionalRows = utils::dependencies::KeyList<>> + struct ModelPlan final : detail::InteriorPlan> + { + using Interior = detail::InteriorPlan>; + using PrescribedCoefficients = Prescribed; + using Values = utils::dependencies::CombinedKeyList::Type, AdditionalValues>; + using Residuals = utils::dependencies::CombinedKeyList; + using Layout = discretization::blocks::SystemLayout; + }; +} diff --git a/src/include/serif/model/specification.hpp b/src/include/serif/model/specification.hpp new file mode 100644 index 00000000..1a80254a --- /dev/null +++ b/src/include/serif/model/specification.hpp @@ -0,0 +1,155 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "serif/dimensions/quantities.hpp" +#include "serif/utils/types/strings.hpp" +#include "serif/utils/types/type_list.hpp" + +namespace serif::model +{ + /** @brief Physical purpose of a prescription, independently of its numerical implementation. */ + enum class SpecificationRole + { + constitutive_law, boundary_condition, invariant, phase_condition, gauge_choice, rotation_law + }; + + /** @brief Distinguish prescribed physics from the physical meaning of a generated scalar. */ + enum class GeneratedStateKind + { + none, multiplier, physical_coordinate, solver_border + }; + + /** @brief Declare independent target, unknown and residual units and stable descriptive identities. */ + template + struct ScalarConstraint + { + using TargetQuantity = Target; + using CoordinateQuantity = Coordinate; + using ResidualQuantity = Residual; + using TargetValue = dimensions::QuantityValue; + static constexpr std::string_view coordinate_id{CoordinateId.data(), CoordinateId.size()}; + static constexpr std::string_view coordinate_symbol{CoordinateSymbol.data(), CoordinateSymbol.size()}; + static constexpr std::string_view residual_id{ResidualId.data(), ResidualId.size()}; + static constexpr std::string_view residual_symbol{ResidualSymbol.data(), ResidualSymbol.size()}; + }; + + template + using Reads = utils::types::TypeList; + + template + using Changes = utils::types::TypeList; + + /** @brief Declare a prescription once; numerical providers and normalization consume this metadata. */ + template , typename ChangedEquations = Changes<>> + struct SpecificationDefinition + { + using SpecificationType = Owner; + using ScalarDescription = Scalar; + using DependsOn = ReadFields; + using Affects = ChangedEquations; + static constexpr std::string_view name{Name.data(), Name.size()}; + static constexpr SpecificationRole role = Role; + static constexpr GeneratedStateKind generated_state_kind = Kind; + }; + + namespace detail + { + /** @brief Reject unknown enums and mismatched physical roles before inspecting optional scalar metadata. */ + [[nodiscard]] consteval bool compatible_generated_state(SpecificationRole role, GeneratedStateKind kind) + { + switch (role) + { + case SpecificationRole::constitutive_law: + case SpecificationRole::boundary_condition: + case SpecificationRole::rotation_law: + return kind == GeneratedStateKind::none; + case SpecificationRole::invariant: + return kind == GeneratedStateKind::multiplier || kind == GeneratedStateKind::physical_coordinate; + case SpecificationRole::phase_condition: + case SpecificationRole::gauge_choice: + return kind == GeneratedStateKind::solver_border; + } + return false; + } + + template + inline constexpr bool is_type_list = false; + + template + inline constexpr bool is_type_list> = true; + } // namespace detail + + /** @brief Recognize a self-describing prescription without demanding executable physics or optional metadata. */ + template + concept ModelSpecification = requires + { + typename std::remove_cvref_t::ModelDefinition; + requires std::same_as::ModelDefinition::SpecificationType, std::remove_cvref_t>; + { std::remove_cvref_t::ModelDefinition::name } -> std::convertible_to; + requires (!std::remove_cvref_t::ModelDefinition::name.empty()); + requires detail::compatible_generated_state(std::remove_cvref_t::ModelDefinition::role, std::remove_cvref_t::ModelDefinition::generated_state_kind); + requires detail::is_type_list::ModelDefinition::DependsOn>; + requires detail::is_type_list::ModelDefinition::Affects>; + }; + + /** @brief Require exact dimensional agreement only when scalar numerical metadata is needed. */ + template + concept DimensionalScalarSpecification = ModelSpecification && requires(const std::remove_cvref_t& specification) + { + requires (std::remove_cvref_t::ModelDefinition::generated_state_kind != GeneratedStateKind::none); + typename std::remove_cvref_t::ModelDefinition::ScalarDescription::TargetQuantity; + typename std::remove_cvref_t::ModelDefinition::ScalarDescription::CoordinateQuantity; + typename std::remove_cvref_t::ModelDefinition::ScalarDescription::ResidualQuantity; + requires dimensions::PhysicalQuantityType::ModelDefinition::ScalarDescription::CoordinateQuantity>; + requires dimensions::PhysicalQuantityType::ModelDefinition::ScalarDescription::ResidualQuantity>; + { specification.target() } -> std::same_as::ModelDefinition::ScalarDescription::TargetQuantity>>; + }; + + /** @brief Name another prescription, including a forward declaration, without erasing its identity. */ + template + struct GeneratedCoordinateOf + { + using SpecificationType = Specification; + static constexpr std::size_t scalar_arity = 1; + static constexpr bool global_scalar = true; + }; + + /** @brief Identify a generated residual independently from its target and unknown units. */ + template + struct ConstraintOf + { + using SpecificationType = Specification; + static constexpr std::size_t scalar_arity = 1; + static constexpr bool global_scalar = true; + }; + + /** @brief Select an explicit numerical enforcement while preserving the public physical prescription. */ + template + struct EnforcementOf + { + using SpecificationType = Specification; + }; + + /** @brief Query scalar units only once the owning prescription is complete. */ + template + using GeneratedCoordinateQuantity = typename Specification::ModelDefinition::ScalarDescription::CoordinateQuantity; + + /** @brief Residual units remain independent from the target and generated coordinate. */ + template + using ConstraintQuantity = typename Specification::ModelDefinition::ScalarDescription::ResidualQuantity; + + /** @brief Inspect canonical prescription identity without claiming numerical availability. */ + struct SpecificationDescriptor + { + std::string_view name; + SpecificationRole role; + GeneratedStateKind generated_state_kind; + }; +} // namespace serif::model diff --git a/src/include/serif/model/specifications/fixed_angular_momentum.hpp b/src/include/serif/model/specifications/fixed_angular_momentum.hpp new file mode 100644 index 00000000..00f795fb --- /dev/null +++ b/src/include/serif/model/specifications/fixed_angular_momentum.hpp @@ -0,0 +1,58 @@ +#pragma once + +#include + +#include "serif/dimensions/type_alias.hpp" +#include "serif/model/specification.hpp" +#include "serif/model/specifications/physical_names.hpp" + +namespace serif::model +{ + /** @brief Prescribe angular momentum while angular velocity remains an independent solved coordinate. */ + class FixedAngularMomentum final + { + public: + struct Parameters + { + dimensions::AngularMomentumValue angular_momentum; + std::array axis{0.0, 0.0, 1.0}; + std::array center{}; + }; + + using ScalarDescription = ScalarConstraint; + using TargetValue = ScalarDescription::TargetValue; + using ModelDefinition = SpecificationDefinition, + Changes>; + + /** @brief Require finite nonnegative angular momentum about the default z axis. */ + explicit FixedAngularMomentum(TargetValue angular_momentum); + + /** @brief Validate the physical target and normalize a finite nonzero rotation axis. */ + explicit FixedAngularMomentum(Parameters parameters); + + /** @brief Inspect the prescribed total angular momentum. */ + [[nodiscard]] TargetValue target() const noexcept + { + return m_target; + } + + /** @brief Read the normalized rotation axis. */ + [[nodiscard]] const std::array& axis() const noexcept + { + return m_axis; + } + + /** @brief Read the point through which the rotation axis passes. */ + [[nodiscard]] const std::array& center() const noexcept + { + return m_center; + } + + private: + TargetValue m_target; + std::array m_axis; + std::array m_center; + }; +} // namespace serif::model diff --git a/src/include/serif/model/specifications/fixed_total_mass.hpp b/src/include/serif/model/specifications/fixed_total_mass.hpp new file mode 100644 index 00000000..c886a3c9 --- /dev/null +++ b/src/include/serif/model/specifications/fixed_total_mass.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include "serif/dimensions/type_alias.hpp" +#include "serif/model/specification.hpp" +#include "serif/model/specifications/physical_names.hpp" + +namespace serif::model +{ + /** @brief Prescribe total stellar mass while the effective surface-potential level remains an independent solved coordinate. */ + class FixedTotalMass final + { + public: + struct Parameters + { + dimensions::MassValue total_mass; + }; + + // Preserve the serialized coordinate identity while describing its present physical meaning. + using ScalarDescription = ScalarConstraint; + using TargetValue = ScalarDescription::TargetValue; + // Us retains the mass-constraint multiplier role; changing its physical units must not reorder public prescriptions. + using ModelDefinition = SpecificationDefinition, Changes>; + + /** @brief Require a positive finite target mass. */ + explicit FixedTotalMass(TargetValue mass); + + /** @brief Construct from named physical parameters. */ + explicit FixedTotalMass(Parameters parameters); + + /** @brief Inspect the prescribed total mass. */ + [[nodiscard]] TargetValue target() const noexcept + { + return m_target; + } + + private: + TargetValue m_target; + }; +} // namespace serif::model diff --git a/src/include/serif/model/specifications/physical_names.hpp b/src/include/serif/model/specifications/physical_names.hpp new file mode 100644 index 00000000..56c2ec47 --- /dev/null +++ b/src/include/serif/model/specifications/physical_names.hpp @@ -0,0 +1,36 @@ +#pragma once + +namespace serif::model +{ + namespace state + { + struct Density + { + }; + struct SurfaceShape + { + }; + struct GravitationalPotential + { + }; + struct OwnGeneratedCoordinate + { + }; + } // namespace state + + namespace equation + { + struct PoissonBalance + { + }; + struct SurfaceShapeBalance + { + }; + struct HydrostaticBalance + { + }; + struct OwnConstraint + { + }; + } // namespace equation +} // namespace serif::model diff --git a/src/include/serif/model/state.hpp b/src/include/serif/model/state.hpp new file mode 100644 index 00000000..183e7ef7 --- /dev/null +++ b/src/include/serif/model/state.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include +#include +#include + +#include "serif/discretization/blocks/vector_view.hpp" + +namespace serif::model { + /** + * @brief Own one stellar candidate in physical coordinates, including any solved surface shape. + * @note Copying a state makes an independent candidate. Its layout fixes the field identities and + * sizes; normalization changes solver coordinates, never these physical coefficients. Prescribed + * inputs are supplied separately to evaluation because they have no solved coordinate in this layout. + */ + template + struct StellarState final { + using Layout = std::remove_cvref_t; + + /** @brief Allocate one zero-initialized coefficient vector in the declared physical ordering. */ + explicit StellarState(Layout layout) : m_layout(std::move(layout)), + m_values(m_layout.template total_size()) { + } + + /** @brief Adopt existing physical coefficients after checking their complete coordinate count. */ + StellarState(Layout layout, std::vector values) : m_layout(std::move(layout)), m_values(std::move(values)) { + (void)discretization::blocks::StateView(m_layout, std::span{m_values}); + } + + /** @brief Inspect the fixed field and equation ordering associated with this state. */ + [[nodiscard]] const Layout& layout() const & noexcept { + return m_layout; + } + + /** @brief Layout references must not outlive a temporary state. */ + const Layout& layout() const && = delete; + + /** @brief Edit one physical field or scalar without maintaining a second coefficient copy. */ + template requires utils::dependencies::ContainsKey + [[nodiscard]] std::span get() & { + return discretization::blocks::StateView(m_layout, std::span{m_values}).template get(); + } + + /** @brief Inspect one physical field while preserving the constness of its owning state. */ + template requires utils::dependencies::ContainsKey + [[nodiscard]] std::span get() const & { + return discretization::blocks::StateView(m_layout, std::span{m_values}).template get(); + } + + /** @brief A field view cannot outlive a temporary owning state. */ + template + std::span get() const && = delete; + + /** @brief Access packed physical coordinates for an explicitly chosen numerical update. */ + [[nodiscard]] std::span values() & noexcept { + return m_values; + } + + /** @brief Read packed physical coordinates without permitting a change in their count. */ + [[nodiscard]] std::span values() const & noexcept { + return m_values; + } + + /** @brief Do not return a span into a temporary owning state. */ + std::span values() const && = delete; + + private: + Layout m_layout; + std::vector m_values; + }; + + /** @brief Infer a state's layout type from its physical coordinate description. */ + template + StellarState(Layout) -> StellarState>; + + /** @brief Infer a state's layout type when adopting physical coefficients. */ + template + StellarState(Layout, std::vector) -> StellarState>; + +} diff --git a/src/include/serif/model/stellar_model.hpp b/src/include/serif/model/stellar_model.hpp new file mode 100644 index 00000000..400ec6f4 --- /dev/null +++ b/src/include/serif/model/stellar_model.hpp @@ -0,0 +1,186 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "serif/model/detail/specification_list.hpp" +#include "serif/model/equilibrium_equations.hpp" + +namespace serif::model +{ + /** @brief Prescriptions accompany one separately owned EOS; no numerical capability is assumed. */ + template + concept ValidSpecificationPack = (ModelSpecification && ...) && + ((Specifications::ModelDefinition::role != SpecificationRole::constitutive_law) && ...) && detail::unique_specification_names(); + + template + class StellarModel; + + /** @brief Own immutable physics and an explicit formulation, independently of meshes, caches and solvers. */ + template + requires ValidSpecificationPack && std::is_object_v && std::is_object_v + class StellarModel> final + { + public: + using EquationOfState = EquationOfStateT; + using Equations = EquationDefinitionT; + using SpecificationTypes = utils::types::TypeList; + static constexpr std::size_t specification_count = sizeof...(Specifications) + 1; + + template + using SpecificationsForRole = std::conditional_t, typename detail::SelectRole::Type>; + + template + static constexpr std::size_t specification_role_count = std::size_t{Role == SpecificationRole::constitutive_law} + (std::size_t{0} + ... + std::size_t{Specifications::ModelDefinition::role == Role}); + + template + static constexpr bool contains_specification = std::same_as, EquationOfState> || (std::same_as, Specifications> || ...); + + /** @brief Capture each physical input once; constructor order does not change specification identity. */ + template + requires detail::ArgumentsMatch::value && (std::constructible_from, Arguments> && ...) + explicit StellarModel(EquationOfState equation_of_state, Equations equations, Arguments&&... specifications) : + m_equation_of_state(std::move(equation_of_state)), m_equations(std::move(equations)), + m_specifications(detail::canonical_specifications(SpecificationTypes{}, std::forward(specifications)...)) + { + } + + /** @brief Inspect the authoritative EOS without requiring its numerical relations. */ + [[nodiscard]] const EquationOfState& equation_of_state() const noexcept + { + return m_equation_of_state; + } + + /** @brief Inspect the supplied equation composition or its deferred definition. */ + [[nodiscard]] const Equations& equations() const noexcept + { + return m_equations; + } + + /** @brief Access a physical prescription by exact type, never by a coincidentally shared unit. */ + template + requires contains_specification + [[nodiscard]] const std::remove_cvref_t& specification() const noexcept + { + if constexpr (std::same_as, EquationOfState>) + { + return m_equation_of_state; + } + else + { + return std::get>(m_specifications); + } + } + + /** @brief Return immutable references to every prescription of one physical role in canonical order. */ + template + [[nodiscard]] auto specifications_for_role() const noexcept + { + return [this](utils::types::TypeList) + { + return std::tuple{specification()...}; + }(SpecificationsForRole{}); + } + + /** @brief Offer a singular physical accessor only when that role is unambiguous. */ + template + requires (specification_role_count == 1) + [[nodiscard]] const auto& specification_for_role() const noexcept + { + return std::get<0>(specifications_for_role()); + } + + /** @brief Visit the EOS first, then physical prescriptions in canonical order. */ + template + requires std::invocable && (std::invocable && ...) + void visit_specifications(Visitor&& visitor) const + { + std::invoke(visitor, m_equation_of_state); + (std::invoke(visitor, specification()), ...); + } + + /** @brief Rebuild same-type physical values; deferred equations will read the new values at binding. */ + template + requires std::copy_constructible && std::invocable && + std::same_as>, EquationOfState> && + (std::invocable && ...) && + (std::same_as>, Specifications> && ...) + [[nodiscard]] StellarModel transform_specifications(Transform&& transform) const + { + return transform_specifications(std::forward(transform), m_equations); + } + + /** @brief Rebuild move-only equation definitions explicitly instead of imposing a cloning policy. */ + template + requires std::invocable && + std::same_as>, EquationOfState> && + (std::invocable && ...) && + (std::same_as>, Specifications> && ...) + [[nodiscard]] StellarModel transform_specifications(Transform&& transform, Equations equations) const + { + return StellarModel{std::invoke(transform, m_equation_of_state), std::move(equations), std::invoke(transform, specification())...}; + } + + /** @brief Describe canonical non-EOS prescriptions without predicting runtime support. */ + [[nodiscard]] static constexpr std::span specification_descriptors() noexcept + { + static constexpr std::array descriptions{{ + {Specifications::ModelDefinition::name, Specifications::ModelDefinition::role, Specifications::ModelDefinition::generated_state_kind}... + }}; + return descriptions; + } + + private: + EquationOfState m_equation_of_state; + Equations m_equations; + std::tuple m_specifications; + }; + + template + requires ValidSpecificationPack...> + StellarModel(EOS&&, Definition&&, Specifications&&...) -> StellarModel, std::remove_cvref_t, typename detail::CanonicalSpecifications...>::Type>; + + /** @brief Construct a typed physical description with programmer-selected equations. */ + template + requires ValidSpecificationPack...> + [[nodiscard]] auto make_stellar_model(EOS&& eos, Definition&& definition, Specifications&&... specifications) + { + return StellarModel{std::forward(eos), std::forward(definition), std::forward(specifications)...}; + } + + namespace detail + { + template + inline constexpr bool is_stellar_model = false; + + template + inline constexpr bool is_stellar_model, std::void_t::SpecificationTypes>> = true; + } // namespace detail + + /** @brief Safely recognize descriptions, including when presented with an incomplete specialization. */ + template + concept StellarModelType = detail::is_stellar_model>; + + /** @brief Inspect role counts without instantiating invalid model specializations. */ + template + inline constexpr std::size_t specification_role_count = [] + { + if constexpr (StellarModelType) + { + return std::remove_cvref_t::template specification_role_count; + } + else + { + return std::size_t{0}; + } + }(); + + template + concept HasUniqueSpecificationForRole = StellarModelType && specification_role_count == 1; + + template + using SpecificationsForRole = typename std::remove_cvref_t::template SpecificationsForRole; +} // namespace serif::model diff --git a/src/include/serif/model/validation.hpp b/src/include/serif/model/validation.hpp new file mode 100644 index 00000000..f68257a7 --- /dev/null +++ b/src/include/serif/model/validation.hpp @@ -0,0 +1,138 @@ +#pragma once + +#include "serif/model/description.hpp" +#include "serif/model/detail/specification_endpoints.hpp" +#include "serif/model/plan.hpp" +#include "serif/operators/linearized_system_view.hpp" + +namespace serif::model +{ + /** @brief Name a physical derivative using the shared block-source key convention. */ + template + using Derivative = operators::PhysicalBlock, utils::dependencies::KeyList>; + + /** @brief Map physical names and built-in generated identities once; custom identities require explicit bindings. */ + template + using SpecificationBindings = utils::dependencies::Bindings< + utils::dependencies::Binding, + utils::dependencies::Binding, + utils::dependencies::Binding, + utils::dependencies::Binding, + utils::dependencies::Binding, + utils::dependencies::Binding, + utils::dependencies::Binding, discretization::blocks::barotropic_constant::mass_normalization>, + utils::dependencies::Binding, discretization::blocks::barotropic_constant::mass_normalization>, + utils::dependencies::Binding, discretization::blocks::rotation::angular_velocity>, + utils::dependencies::Binding, discretization::blocks::rotation::angular_velocity>, Additional...>; + + using DefaultSpecificationBindings = SpecificationBindings<>; + + namespace detail + { + template + struct EquationCouplings; + + template + struct EquationCouplings> + { + using Type = utils::dependencies::KeyList...>; + }; + + template + struct InteriorCouplings; + + template + struct InteriorCouplings> + { + using Type = utils::dependencies::CombinedKeyList::Columns>::Type...>; + }; + + /** @brief Check physical inputs separately from unknowns and explicitly inactive equations. */ + template + [[nodiscard]] consteval bool valid_specification_endpoints() + { + using namespace utils::dependencies; + using Endpoints = ResolvedSpecification; + if constexpr (!ContainsKey || !ContainsKey || + !Endpoints::Reads::template owners_present || !Endpoints::Changes::template owners_present) + { + return false; + } + return [](KeyList, KeyList) + { + return sizeof...(Read) != 0 && (ContainsKey && ...) && (ContainsKey || ...) && + ((ContainsKey || ContainsKey) && ...); + }(typename Endpoints::Reads::Keys{}, typename Endpoints::Changes::Keys{}); + } + + /** @brief Fold independent physical declarations; no equation counting or rank inference is involved. */ + template + [[nodiscard]] consteval bool valid_specifications(utils::types::TypeList) + { + using namespace utils::dependencies; + constexpr bool declarations = (((!GeneratesState::value && std::same_as> && + std::same_as>) || DimensionalScalarSpecification) && ...); + if constexpr (!declarations) + { + return false; + } + else + { + return [](utils::types::TypeList) + { + using Coordinates = CombinedKeyList::Coordinate>...>; + using Residuals = CombinedKeyList::Constraint>...>; + return key_count_v == sizeof...(Generated) && key_count_v == sizeof...(Generated) && + (valid_specification_endpoints() && ...); + }(GeneratedSpecifications{}); + } + } + + /** @brief Validate declared derivative contracts without inferring zeros from missing providers. */ + template + [[nodiscard]] consteval bool specification_derivatives_available() + { + using namespace utils::dependencies; + using Endpoints = ResolvedSpecification; + return [](KeyList, KeyList) + { + constexpr auto changed_row = []() + { + return ContainsKey, Couplings> && (ContainsKey, Couplings> && ...); + }; + return (ContainsKey, Couplings> && ...) && (changed_row.template operator()() && ...); + }(typename Endpoints::Reads::Keys{}, typename Endpoints::Changes::Keys{}); + } + + /** @brief Check each prescription against independently supplied derivative registrations. */ + template + [[nodiscard]] consteval bool specification_derivatives_available(utils::types::TypeList) + { + return (specification_derivatives_available() && ...); + } + } // namespace detail + + /** @brief Derive declared interior couplings from actual operands rather than a second physical sparsity table. */ + template + using DeclaredEquationCouplings = typename detail::InteriorCouplings>::Type; + + /** @brief Validate physical ownership, available inputs and uniquely owned generated equations. */ + template , typename Bindings = DefaultSpecificationBindings> + concept SpecificationEndpointsValid = StellarModelType && utils::dependencies::IsKeyList && utils::dependencies::IsKeyList && + utils::dependencies::IsKeyList && utils::dependencies::IsKeyList && utils::dependencies::IsBindings && + detail::valid_specifications, Values, Inputs, Rows, InactiveRows, Bindings>(typename std::remove_cvref_t::SpecificationTypes{}); + + /** @brief Validate physical endpoints and every derivative in an explicitly declared contract. */ + template + concept SpecificationTopologyValid = SpecificationEndpointsValid, Bindings> && + utils::dependencies::IsKeyList && detail::specification_derivatives_available(detail::GeneratedSpecifications>{}); + + /** + * @brief Require a registered implementation for every declared generated derivative, including explicit zeros. + * @note ImplementedCouplings must come from admitted providers. Missing methods never imply zero; this checks + * specification contributions only, leaving the equation backend's own execution contract independent. + */ + template + concept SpecificationBackendAvailable = SpecificationTopologyValid && + SpecificationTopologyValid; +} // namespace serif::model diff --git a/src/lib/serif/model/specifications.cpp b/src/lib/serif/model/specifications.cpp new file mode 100644 index 00000000..80d7af35 --- /dev/null +++ b/src/lib/serif/model/specifications.cpp @@ -0,0 +1,60 @@ +#include +#include +#include + +#include "serif/model/specifications/fixed_angular_momentum.hpp" +#include "serif/model/specifications/fixed_total_mass.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::model +{ + /** @brief Validate a dimensional scalar target before storing a physical prescription. */ + static void validate_target(double value, bool zero_allowed, std::string_view name) + { + using namespace utils::errors; + serif_ensure(std::isfinite(value) && (zero_allowed ? value >= 0.0 : value > 0.0), + "The {} target must be finite and {}; received {}. Supply a physical target in the permitted range.", name, zero_allowed ? "nonnegative" : "positive", value); + } + + /** @brief Store a physically admissible total mass. */ + FixedTotalMass::FixedTotalMass(TargetValue mass) : m_target(mass) + { + validate_target(mass.value(), false, "total mass"); + } + + /** @brief Delegate named mass parameters to the dimensional constructor. */ + FixedTotalMass::FixedTotalMass(Parameters parameters) : FixedTotalMass(parameters.total_mass) + { + } + + /** @brief Use the default rotation axis for a dimensional angular momentum target. */ + FixedAngularMomentum::FixedAngularMomentum(TargetValue angular_momentum) : FixedAngularMomentum(Parameters{.angular_momentum = angular_momentum}) + { + } + + /** @brief Preserve arbitrary finite axes and centers without overflowing the normalization calculation. */ + FixedAngularMomentum::FixedAngularMomentum(Parameters parameters) : m_target(parameters.angular_momentum), m_axis(parameters.axis), m_center(parameters.center) + { + using namespace utils::errors; + validate_target(m_target.value(), true, "angular momentum"); + + double maximum_component = 0.0; + for (std::size_t component = 0; component < m_axis.size(); ++component) + { + serif_ensure(std::isfinite(m_axis[component]) && std::isfinite(m_center[component]), + "Rotation axis and center components must be finite; component {} contains axis {} and center {}. Supply a finite rotation axis and center.", component, m_axis[component], m_center[component]); + maximum_component = std::max(maximum_component, std::abs(m_axis[component])); + } + serif_ensure(maximum_component > 0.0, "A rotation axis must be nonzero; supply a direction such as (0, 0, 1)."); + + for (double& component : m_axis) + { + component /= maximum_component; + } + const double norm = std::hypot(m_axis[0], m_axis[1], m_axis[2]); + for (double& component : m_axis) + { + component /= norm; + } + } +} // namespace serif::model