From 24d4c20f82e75964b746df49eec1320afa675fd3 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 9 Nov 2026 13:30:20 -0500 Subject: [PATCH] feat(solver/normalization): Added normalization system The dynamical range of different fields SERiF must solve for is massive. This just breaks any reasonable preconditioner / spreads the eigen modes of the system matrix out. This is a normalization system which uses charectaristic scale factors (e.g. the approximate radius estimated using the lane-emden equation) to rescale different fields to a more similar numeric scale. This turns out to be essential for stable convergence. Broadly this implements Riesz normalization along with exact diagnoal normalization. --- .../solver/normalization/computations.hpp | 76 +++++++ .../solver/normalization/configuration.hpp | 95 ++++++++ .../detail/factor_computations.hpp | 213 ++++++++++++++++++ .../solver/normalization/detail/packing.hpp | 13 ++ .../normalization/detail/validation.hpp | 18 ++ .../serif/solver/normalization/diagonal.hpp | 86 +++++++ .../solver/normalization/dimensionless.hpp | 43 ++++ .../serif/solver/normalization/keys.hpp | 37 +++ .../serif/solver/normalization/layout.hpp | 14 ++ .../serif/solver/normalization/metrics.hpp | 85 +++++++ .../solver/normalization/normalization.hpp | 10 + .../normalization/physical_coordinates.hpp | 44 ++++ .../serif/solver/normalization/plan.hpp | 170 ++++++++++++++ .../serif/solver/normalization/point_rows.hpp | 15 ++ .../serif/solver/normalization/scales.hpp | 115 ++++++++++ .../serif/solver/normalization/diagonal.cpp | 89 ++++++++ .../solver/normalization/dimensionless.cpp | 46 ++++ .../serif/solver/normalization/packing.cpp | 52 +++++ .../solver/normalization/riesz_factors.cpp | 135 +++++++++++ src/lib/serif/solver/normalization/scales.cpp | 62 +++++ .../serif/solver/normalization/validation.cpp | 49 ++++ 21 files changed, 1467 insertions(+) create mode 100644 src/include/serif/solver/normalization/computations.hpp create mode 100644 src/include/serif/solver/normalization/configuration.hpp create mode 100644 src/include/serif/solver/normalization/detail/factor_computations.hpp create mode 100644 src/include/serif/solver/normalization/detail/packing.hpp create mode 100644 src/include/serif/solver/normalization/detail/validation.hpp create mode 100644 src/include/serif/solver/normalization/diagonal.hpp create mode 100644 src/include/serif/solver/normalization/dimensionless.hpp create mode 100644 src/include/serif/solver/normalization/keys.hpp create mode 100644 src/include/serif/solver/normalization/layout.hpp create mode 100644 src/include/serif/solver/normalization/metrics.hpp create mode 100644 src/include/serif/solver/normalization/physical_coordinates.hpp create mode 100644 src/include/serif/solver/normalization/point_rows.hpp create mode 100644 src/include/serif/solver/normalization/scales.hpp create mode 100644 src/lib/serif/solver/normalization/diagonal.cpp create mode 100644 src/lib/serif/solver/normalization/dimensionless.cpp create mode 100644 src/lib/serif/solver/normalization/packing.cpp create mode 100644 src/lib/serif/solver/normalization/riesz_factors.cpp create mode 100644 src/lib/serif/solver/normalization/scales.cpp create mode 100644 src/lib/serif/solver/normalization/validation.cpp diff --git a/src/include/serif/solver/normalization/computations.hpp b/src/include/serif/solver/normalization/computations.hpp new file mode 100644 index 00000000..38f26204 --- /dev/null +++ b/src/include/serif/solver/normalization/computations.hpp @@ -0,0 +1,76 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "serif/solver/normalization/detail/factor_computations.hpp" +#include "serif/solver/normalization/detail/packing.hpp" + +namespace serif::solver::normalization +{ + namespace detail + { + /** @brief Test every coordinate factory without instantiating an unusable runtime graph. */ + template + [[nodiscard]] consteval bool can_prepare_coordinates(utils::dependencies::KeyList) + { + return (requires + { + coordinate_computation(utils::dependencies::bound_type_t{}); + } && ...); + } + + /** @brief Expand the supplied layout, looking up methods without imposing plan declaration order. */ + template + [[nodiscard]] auto coordinate_computations(utils::dependencies::KeyList) + { + return std::tuple{coordinate_computation(utils::dependencies::bound_type_t{})...}; + } + + /** @brief Freeze the independently ordered state and residual coordinates into one owning result. */ + template + [[nodiscard]] auto packing_computation(utils::dependencies::KeyList, utils::dependencies::KeyList) + { + using namespace utils::dependencies; + using discretization::blocks::RepresentedValueKind; + + return make_computation..., keys::Factors...>>( + [](const auto&... factors) -> std::expected + { + const std::array, sizeof...(Values) + sizeof...(Residuals)> coordinates{factors...}; + const std::span> blocks{coordinates}; + return pack_normalization(blocks.first(sizeof...(Values)), blocks.subspan(sizeof...(Values))); + }); + } + } + + template + concept NormalizationPreparable = CoordinateLayoutType && + CompleteNormalizationFor::Values, typename std::remove_cvref_t::Residuals> && + detail::can_prepare_coordinates::Values, std::remove_cvref_t, std::remove_cvref_t, discretization::blocks::RepresentedValueKind::value>(typename std::remove_cvref_t::Values{}) && + detail::can_prepare_coordinates::Residuals, std::remove_cvref_t, std::remove_cvref_t, discretization::blocks::RepresentedValueKind::residual>(typename std::remove_cvref_t::Residuals{}); + + /** + * @brief Declare normalization as ordinary dependency computations; no metrics or factors are prepared here. + * @note Supply keys::Layout, keys::Scales, and the plan's metric/selection/prescription sources. Copy the accepted DiagonalNormalization to freeze one solve before preparing later candidates. + */ + template + requires NormalizationPreparable + [[nodiscard]] auto make_normalization_computations(const LayoutT&) + { + using discretization::blocks::RepresentedValueKind; + using Plan = std::remove_cvref_t; + using Layout = std::remove_cvref_t; + using Scales = std::remove_cvref_t; + + return std::tuple_cat( + detail::coordinate_computations(typename Layout::Values{}), + detail::coordinate_computations(typename Layout::Residuals{}), + std::tuple{detail::packing_computation(typename Layout::Values{}, typename Layout::Residuals{})}); + } +} diff --git a/src/include/serif/solver/normalization/configuration.hpp b/src/include/serif/solver/normalization/configuration.hpp new file mode 100644 index 00000000..8ed78f97 --- /dev/null +++ b/src/include/serif/solver/normalization/configuration.hpp @@ -0,0 +1,95 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "serif/dimensions/type_alias.hpp" +#include "serif/solver/normalization/scales.hpp" +#include "serif/utils/dependencies/dependency_context.hpp" +#include "serif/utils/laneemden/laneemden.hpp" + +namespace serif::solver::normalization { + /** @brief Obtain a reference value from the corresponding descriptive model prescription. */ + struct FromModel final {}; + + /** @brief Derive the coordinate rules from the selected fields, equations and boundary formulation. */ + struct AutomaticPlan final {}; + + /** @brief Remove physical and geometric units from the solver's coordinates and norms. */ + struct Dimensionless final {}; + + /** + * @brief Estimate a reference radius from a Lane–Emden shape, mass and central density. + * @note These describe the numerical reference star; they never alter the model's EOS or constraints. + */ + template + requires (std::same_as || std::same_as) && + (std::same_as || std::same_as) + struct LaneEmdenApproximation final { + double n{3.0}; + Mass mass; + CentralDensity central_density; + utils::laneemden::IntegrationOptions integration{}; + }; + + /** + * @brief Declare normalization without evaluating a reference star or allocating coordinate factors. + * @note G uses the model's declared units; its default is the existing G=1 code-unit convention. + * Advanced users may supply a complete plan and ordinary graph sources/computations for custom rules. + */ + template , typename Coordinates = Dimensionless, + typename Plan = AutomaticPlan, typename Sources = utils::dependencies::DependencyContext<>, typename Computations = std::tuple<>> + requires std::same_as && utils::dependencies::IsDependencyContext + struct NormalizationOptions final { + Scheme scheme; + Coordinates coordinates{}; + Plan plan; + dimensions::QuantityValue gravitational_constant{1.0}; + Sources sources; + Computations computations; + }; + + /** @brief Preserve the deduced numerical choices for later preparation with a model and discretization. */ + template + [[nodiscard]] auto make_normalization(NormalizationOptions options) { + return options; + } + + /** @brief Select the usual dimensionless Lane–Emden reference; supply its central density before resolving it against a model. */ + [[nodiscard]] inline auto make_normalization() { + return NormalizationOptions<>{}; + } + + /** @brief Derive the shared stellar scales after both dimensional reference inputs have been resolved. */ + [[nodiscard]] inline std::expected prepare_reference_scales( + const LaneEmdenApproximation& scheme, + dimensions::QuantityValue gravitational_constant) { + const auto solution = utils::laneemden::solve(scheme.n, scheme.integration); + if (!solution) { + return std::unexpected(solution.error()); + } + const auto reference = utils::laneemden::scale_from_mass(*solution, scheme.mass, scheme.central_density, gravitational_constant); + if (!reference) { + return std::unexpected(reference.error()); + } + return prepare_stellar_reference_scales(reference->mass, reference->radius, gravitational_constant); + } + + /** @brief Use explicitly prepared reference scales through the same graph as estimated references. */ + [[nodiscard]] inline std::expected prepare_reference_scales( + const StellarReferenceScales& reference, dimensions::QuantityValue gravitational_constant) { + if (reference.gravitational_constant != gravitational_constant) { + return std::unexpected("The explicit normalization reference and configuration use different gravitational constants. Supply one G in the model's units."); + } + return prepare_stellar_reference_scales(reference.mass, reference.radius, gravitational_constant); + } + + /** @brief A custom stellar reference strategy supplies the same typed scales and ordinary failure value. */ + template + concept ReferenceScheme = requires(const Scheme& scheme, dimensions::QuantityValue gravity) { + { prepare_reference_scales(scheme, gravity) } -> std::same_as>; + }; +} diff --git a/src/include/serif/solver/normalization/detail/factor_computations.hpp b/src/include/serif/solver/normalization/detail/factor_computations.hpp new file mode 100644 index 00000000..e3d1f926 --- /dev/null +++ b/src/include/serif/solver/normalization/detail/factor_computations.hpp @@ -0,0 +1,213 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "serif/solver/normalization/keys.hpp" +#include "serif/solver/normalization/layout.hpp" +#include "serif/solver/normalization/plan.hpp" +#include "serif/solver/normalization/diagonal.hpp" +#include "serif/solver/normalization/detail/validation.hpp" +#include "serif/solver/normalization/point_rows.hpp" +#include "serif/solver/normalization/scales.hpp" +#include "serif/utils/dependencies/computation.hpp" + +namespace serif::solver::normalization::detail +{ + template + concept ProvidesScale = requires(const ScalesT& scales) + { + { characteristic_scale(scales, QuantityT{}) } -> std::same_as>; + }; + + /** @brief Reject a layout that changes the arity of a declared global scalar. */ + template + [[nodiscard]] std::expected coordinate_size(const LayoutT& layout) + { + const std::size_t size = layout.template size(); + using Term = typename decltype(metric_term())::type; + return validate_coordinate_size(size, discretization::blocks::GlobalScalarTerm); + } + + /** @brief Check supplied metrics or custom factors before publishing a coordinate result. */ + template + [[nodiscard]] std::expected check_coordinate_size(const LayoutT& layout, std::size_t received) + { + const auto size = coordinate_size(layout); + if (!size) + { + return std::unexpected(size.error()); + } + return check_coordinate_size(*size, received); + } + + /** @brief Identity preparation requires a layout but no scale or metric sources. */ + template + [[nodiscard]] auto coordinate_computation(IdentityCoordinate) + { + using namespace utils::dependencies; + + return make_computation, InputList>( + [](const LayoutT& layout) -> std::expected, std::string> + { + const auto size = coordinate_size(layout); + if (!size) + { + return std::unexpected(size.error()); + } + return std::vector(*size, 1.0); + }); + } + + /** @brief Algebraic coordinates use a scalar metric without allocating an intermediate vector of ones. */ + template + requires ProvidesScale + [[nodiscard]] auto coordinate_computation(PhysicalRieszCoordinate) + { + using namespace utils::dependencies; + + return make_computation, InputList>( + [](const LayoutT& layout, const ScalesT& scales) -> std::expected, std::string> + { + const auto size = coordinate_size(layout); + if (!size) + { + return std::unexpected(size.error()); + } + const auto factor = make_riesz_factor(Kind, 1.0, characteristic_scale(scales, QuantityT{}).value()); + if (!factor) + { + return std::unexpected(factor.error()); + } + return std::vector(*size, *factor); + }); + } + + /** @brief Both coordinate roles consume the same cached reference metric with different primal/dual formulas. */ + template + requires (!std::same_as && ProvidesScale) + [[nodiscard]] auto coordinate_computation(PhysicalRieszCoordinate) + { + using namespace utils::dependencies; + + return make_computation, InputList>( + [](const LayoutT& layout, const ScalesT& scales, const std::vector& metric) -> std::expected, std::string> + { + const auto valid = check_coordinate_size(layout, metric.size()); + if (!valid) + { + return std::unexpected(valid.error()); + } + return make_riesz_factors(Kind, metric, characteristic_scale(scales, QuantityT{}).value()); + }); + } + + /** @brief Read the replacement's quantity independently of the bulk equation's scale. */ + template + requires (Kind == discretization::blocks::RepresentedValueKind::residual && !std::same_as && ProvidesScale && ProvidesScale) + [[nodiscard]] auto coordinate_computation(WithPointRows, PointQuantityT, SelectionT>) + { + using namespace utils::dependencies; + + return make_computation, InputList>( + [](const LayoutT& layout, const ScalesT& scales, const std::vector& metric, const std::vector& rows) -> std::expected, std::string> + { + const auto valid = check_coordinate_size(layout, metric.size()); + if (!valid) + { + return std::unexpected(valid.error()); + } + return make_point_replaced_residual_factors(metric, characteristic_scale(scales, QuantityT{}).value(), rows, characteristic_scale(scales, PointQuantityT{}).value()); + }); + } + + /** @brief Global scalars have no geometric measure to remove. */ + template + requires ProvidesScale + [[nodiscard]] auto coordinate_computation(DimensionlessRieszCoordinate) { + return coordinate_computation(PhysicalRieszCoordinate{}); + } + + /** @brief Cancel the reference measure while retaining the native finite-element coordinate representation. */ + template + requires (!std::same_as && ProvidesScale && ProvidesScale) + [[nodiscard]] auto coordinate_computation(DimensionlessRieszCoordinate) { + using namespace utils::dependencies; + + return make_computation, InputList>( + [](const LayoutT& layout, const ScalesT& scales, const std::vector& metric) -> std::expected, std::string> { + const auto valid = check_coordinate_size(layout, metric.size()); + if (!valid) { + return std::unexpected(valid.error()); + } + const auto scale = dimensionless_riesz_scale(Kind, Geometry, characteristic_scale(scales, QuantityT{}).value(), characteristic_scale(scales, dimensions::Length{}).value()); + if (!scale) { + return std::unexpected(scale.error()); + } + return make_riesz_factors(Kind, metric, *scale); + }); + } + + /** @brief Remove geometric units only from weak rows, leaving each strong point equation in its own physical units. */ + template + requires (Kind == discretization::blocks::RepresentedValueKind::residual && !std::same_as && ProvidesScale && ProvidesScale && ProvidesScale) + [[nodiscard]] auto coordinate_computation(WithPointRows, PointQuantityT, SelectionT>) { + using namespace utils::dependencies; + + return make_computation, InputList>( + [](const LayoutT& layout, const ScalesT& scales, const std::vector& metric, const std::vector& rows) -> std::expected, std::string> { + const auto valid = check_coordinate_size(layout, metric.size()); + if (!valid) { + return std::unexpected(valid.error()); + } + const auto scale = dimensionless_riesz_scale(Kind, Geometry, characteristic_scale(scales, QuantityT{}).value(), characteristic_scale(scales, dimensions::Length{}).value()); + if (!scale) { + return std::unexpected(scale.error()); + } + return make_point_replaced_residual_factors(metric, *scale, rows, characteristic_scale(scales, PointQuantityT{}).value()); + }); + } + + /** @brief Invoke only the exact declared prescription, supplied as a tracked source rather than a hidden mutable capture. */ + template + requires requires(const std::remove_cvref_t& prescription, std::size_t size) + { + { prescription.prepare_factors(keys::Factors{}, size) } -> std::same_as, std::string>>; + } + [[nodiscard]] auto coordinate_computation(RuntimePreparedCoordinate) + { + using namespace utils::dependencies; + using Prescription = std::remove_cvref_t; + + return make_computation, InputList>>( + [](const LayoutT& layout, const SuppliedPrescription& prescription) -> std::expected, std::string> + requires std::same_as, Prescription> + { + const auto size = coordinate_size(layout); + if (!size) + { + return std::unexpected(size.error()); + } + auto factors = prescription.prepare_factors(keys::Factors{}, *size); + if (!factors) + { + return factors; + } + const auto valid_size = check_coordinate_size(*size, factors->size()); + if (!valid_size) + { + return std::unexpected(valid_size.error()); + } + const auto valid_factors = validate_factors(*factors); + if (!valid_factors) + { + return std::unexpected(valid_factors.error()); + } + return factors; + }); + } +} diff --git a/src/include/serif/solver/normalization/detail/packing.hpp b/src/include/serif/solver/normalization/detail/packing.hpp new file mode 100644 index 00000000..eb923bd9 --- /dev/null +++ b/src/include/serif/solver/normalization/detail/packing.hpp @@ -0,0 +1,13 @@ +#pragma once + +#include +#include +#include + +#include "serif/solver/normalization/diagonal.hpp" + +namespace serif::solver::normalization::detail +{ + /** @brief Freeze independently ordered state and residual blocks into their packed solver coordinates. */ + [[nodiscard]] std::expected pack_normalization(std::span> states, std::span> residuals); +} diff --git a/src/include/serif/solver/normalization/detail/validation.hpp b/src/include/serif/solver/normalization/detail/validation.hpp new file mode 100644 index 00000000..2fd9b8bf --- /dev/null +++ b/src/include/serif/solver/normalization/detail/validation.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include +#include +#include +#include + +namespace serif::solver::normalization::detail +{ + /** @brief Require finite positive factors before accepting custom or packed normalization. */ + [[nodiscard]] std::expected validate_factors(std::span factors); + + /** @brief Check representable coordinate storage and the fixed arity of global scalar terms. */ + [[nodiscard]] std::expected validate_coordinate_size(std::size_t size, bool global_scalar); + + /** @brief Compare metric or factor size with the independently declared coordinate layout. */ + [[nodiscard]] std::expected check_coordinate_size(std::size_t expected, std::size_t received); +} diff --git a/src/include/serif/solver/normalization/diagonal.hpp b/src/include/serif/solver/normalization/diagonal.hpp new file mode 100644 index 00000000..07c53f0d --- /dev/null +++ b/src/include/serif/solver/normalization/diagonal.hpp @@ -0,0 +1,86 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include "serif/discretization/blocks/base.hpp" + +namespace serif::solver::normalization +{ + namespace detail + { + /** @brief Require matching coordinate sizes before any output is modified. */ + void require_normalization_size(std::size_t received, std::size_t expected); + + /** @brief Apply a diagonal in place or to disjoint storage, rejecting partial overlap. */ + void apply_normalization(std::span factors, std::span input, std::span output, bool inverse); + } + + /** @brief Prepare one primal sqrt(metric)/scale or dual 1/(scale sqrt(metric)) factor. */ + [[nodiscard]] std::expected make_riesz_factor(discretization::blocks::RepresentedValueKind kind, double metric, double physical_scale); + + /** @brief Prepare diagonal Riesz factors, reporting which coordinate has invalid data. */ + [[nodiscard]] std::expected, std::string> make_riesz_factors(discretization::blocks::RepresentedValueKind kind, std::span metric, double physical_scale); + + /** @brief Own the frozen, independent state and residual transformations used by one solve. */ + class DiagonalNormalization final + { + public: + /** @brief Read the factors taking physical state coefficients to normalized coordinates. */ + [[nodiscard]] std::span state_factors() const noexcept + { + return m_state_factors; + } + + /** @brief Read the factors taking physical residual coefficients to normalized coordinates. */ + [[nodiscard]] std::span residual_factors() const noexcept + { + return m_residual_factors; + } + + /** @brief Number of state coordinates. */ + [[nodiscard]] std::size_t state_size() const noexcept + { + return m_state_factors.size(); + } + + /** @brief Number of residual coordinates, which need not equal the state size. */ + [[nodiscard]] std::size_t residual_size() const noexcept + { + return m_residual_factors.size(); + } + + /** @brief Convert physical state coefficients to normalized coordinates. */ + void normalize_state(std::span physical, std::span normalized) const; + + /** @brief Recover physical state coefficients from normalized coordinates. */ + void denormalize_state(std::span normalized, std::span physical) const; + + /** @brief Convert physical residual coefficients to normalized coordinates. */ + void normalize_residual(std::span physical, std::span normalized) const; + + /** @brief Recover physical residual coefficients from normalized coordinates. */ + void denormalize_residual(std::span normalized, std::span physical) const; + + /** @brief Squared primal norm of a physical state vector. */ + [[nodiscard]] double state_norm_squared(std::span physical) const; + + /** @brief Squared dual norm of a physical residual vector. */ + [[nodiscard]] double residual_norm_squared(std::span physical) const; + + private: + friend std::expected make_diagonal_normalization(std::vector, std::vector); + + /** @brief Accept factors only after the factory has validated them. */ + DiagonalNormalization(std::vector state, std::vector residual); + + std::vector m_state_factors; + std::vector m_residual_factors; + }; + + /** @brief Validate and freeze independent coordinate factors, including empty or rectangular systems. */ + [[nodiscard]] std::expected make_diagonal_normalization(std::vector state, std::vector residual); +} diff --git a/src/include/serif/solver/normalization/dimensionless.hpp b/src/include/serif/solver/normalization/dimensionless.hpp new file mode 100644 index 00000000..fd86e660 --- /dev/null +++ b/src/include/serif/solver/normalization/dimensionless.hpp @@ -0,0 +1,43 @@ +#pragma once + +#include +#include + +#include "serif/solver/normalization/metrics.hpp" + +namespace serif::solver::normalization { + /** @brief Identify the geometric convention of a reference coordinate metric. */ + enum class ReferenceGeometry { + volume, + surface, + algebraic + }; + + /** @brief Require custom metrics to declare their geometric convention explicitly. */ + template + struct DimensionlessMetricGeometry; + + template + struct DimensionlessMetricGeometry> { + static constexpr ReferenceGeometry value = ReferenceGeometry::volume; + }; + + template <> + struct DimensionlessMetricGeometry { + static constexpr ReferenceGeometry value = ReferenceGeometry::surface; + }; + + template <> + struct DimensionlessMetricGeometry { + static constexpr ReferenceGeometry value = ReferenceGeometry::algebraic; + }; + + namespace detail { + /** + * @brief Include reference geometry in the scale consumed by the existing Riesz factor calculation. + * @note Volume residual quantities are densities before weak integration; surface residual quantities are integrated forces. + */ + [[nodiscard]] std::expected dimensionless_riesz_scale( + discretization::blocks::RepresentedValueKind kind, ReferenceGeometry geometry, double physical_scale, double reference_radius); + } +} diff --git a/src/include/serif/solver/normalization/keys.hpp b/src/include/serif/solver/normalization/keys.hpp new file mode 100644 index 00000000..1e1ebc87 --- /dev/null +++ b/src/include/serif/solver/normalization/keys.hpp @@ -0,0 +1,37 @@ +#pragma once + +#include "serif/discretization/blocks/base.hpp" +#include "serif/discretization/blocks/keys.hpp" + +namespace serif::solver::normalization::keys +{ + /** @brief Normalization consumes the same coordinate layout as every other numerical service. */ + using Layout = discretization::blocks::keys::Layout; + + struct Scales + { + }; + + struct Normalization + { + }; + + /** @brief Distinguish state and residual factors even when they share the same physical term. */ + template + struct Factors + { + using Coordinate = CoordinateT; + static constexpr discretization::blocks::RepresentedValueKind kind = Kind; + }; + + /** @brief Bind the exact custom prescription named by a runtime-prepared coordinate. */ + template + struct Prescription + { + }; + + template + struct ReferenceBinding + { + }; +} diff --git a/src/include/serif/solver/normalization/layout.hpp b/src/include/serif/solver/normalization/layout.hpp new file mode 100644 index 00000000..1b69de6c --- /dev/null +++ b/src/include/serif/solver/normalization/layout.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include "serif/discretization/blocks/system_layout.hpp" + +namespace serif::solver::normalization +{ + /** @brief Compatibility spelling for the shared discretization layout. */ + template + using CoordinateLayout = discretization::blocks::SystemLayout; + + /** @brief Normalization consumes the same coordinate representation as physical operators and solvers. */ + template + concept CoordinateLayoutType = discretization::blocks::SystemLayoutType; +} diff --git a/src/include/serif/solver/normalization/metrics.hpp b/src/include/serif/solver/normalization/metrics.hpp new file mode 100644 index 00000000..c982b7e2 --- /dev/null +++ b/src/include/serif/solver/normalization/metrics.hpp @@ -0,0 +1,85 @@ +#pragma once + +#include +#include + +#include "serif/discretization/blocks/concepts.hpp" +#include "serif/utils/dependencies/key_list_operations.hpp" + +namespace serif::geometry::deformation::keys +{ + struct SurfaceMetric; +} + +namespace serif::solver::normalization +{ + /** @brief A unit Gram metric for algebraic coordinates, still subject to physical scaling. */ + struct UnitMetric final + { + }; + + /** @brief Identify a reference volume metric by its physical term, including its support and FE representation. */ + template + struct ReferenceVolumeMetric final + { + using Coordinates = utils::dependencies::KeyList; + }; + + /** @brief Metric keys advertise the coordinate identities they measure; a same-sized unrelated metric is insufficient. */ + template + struct MetricCoordinates + { + using Type = utils::dependencies::KeyList<>; + }; + + template + struct MetricCoordinates> + { + using Type = typename MetricT::Coordinates; + }; + + /** @brief Permit a Euclidean norm of algebraic surface moments after their measure is made dimensionless. */ + template <> + struct MetricCoordinates + { + using Type = utils::dependencies::KeyList; + }; + + /** @brief Consume deformation's existing metric directly in both native shape-coordinate roles. */ + template <> + struct MetricCoordinates + { + using Type = utils::dependencies::KeyList; + }; + + namespace detail + { + /** @brief Resolve assembled row identity without confusing a point-border row with its sampled field. */ + template + [[nodiscard]] consteval auto metric_term() + { + if constexpr (requires + { + typename CoordinateT::Border; + }) + { + return std::type_identity{}; + } + else if constexpr (requires + { + typename CoordinateT::Test; + }) + { + return std::type_identity{}; + } + else + { + return std::type_identity{}; + } + } + } + + template + concept MetricFor = (std::same_as && discretization::blocks::GlobalScalarTerm())::type>) || + utils::dependencies::ContainsKey())::type, typename MetricCoordinates::Type>; +} diff --git a/src/include/serif/solver/normalization/normalization.hpp b/src/include/serif/solver/normalization/normalization.hpp index e69de29b..a0316ec5 100644 --- a/src/include/serif/solver/normalization/normalization.hpp +++ b/src/include/serif/solver/normalization/normalization.hpp @@ -0,0 +1,10 @@ +#pragma once + +#include "serif/solver/normalization/plan.hpp" +#include "serif/solver/normalization/physical_coordinates.hpp" +#include "serif/solver/normalization/layout.hpp" +#include "serif/solver/normalization/scales.hpp" +#include "serif/solver/normalization/diagonal.hpp" +#include "serif/solver/normalization/computations.hpp" + +#include "serif/solver/normalization/configuration.hpp" diff --git a/src/include/serif/solver/normalization/physical_coordinates.hpp b/src/include/serif/solver/normalization/physical_coordinates.hpp new file mode 100644 index 00000000..03e3cd28 --- /dev/null +++ b/src/include/serif/solver/normalization/physical_coordinates.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include +#include + +#include "serif/solver/normalization/plan.hpp" + +namespace serif::solver::normalization +{ + namespace detail + { + /** @brief Select ordinary FE or algebraic metrics from the existing physical representation. */ + template + [[nodiscard]] consteval auto default_metric() + { + using Term = typename decltype(metric_term())::type; + if constexpr (discretization::blocks::SpatialTerm) + { + return std::type_identity>{}; + } + else if constexpr (discretization::blocks::GlobalScalarTerm) + { + return std::type_identity{}; + } + else + { + return std::type_identity{}; + } + } + } + + /** @brief Describe an unknown from its existing quantity, support, and functional-space metadata. */ + template + requires (!std::same_as())::type, void>) + using PhysicalValueRule = utils::dependencies::Binding())::type>>; + + /** @brief Describe a residual using the equation's declared quantity, which may differ from the unknown it closes. */ + template + requires requires + { + typename EquationT::ResidualQuantity; + } && (!std::same_as())::type, void>) + using PhysicalResidualRule = utils::dependencies::Binding())::type>>; +} diff --git a/src/include/serif/solver/normalization/plan.hpp b/src/include/serif/solver/normalization/plan.hpp index e69de29b..81cda0db 100644 --- a/src/include/serif/solver/normalization/plan.hpp +++ b/src/include/serif/solver/normalization/plan.hpp @@ -0,0 +1,170 @@ +#pragma once + +#include +#include + +#include "serif/dimensions/quantities.hpp" +#include "serif/discretization/blocks/base.hpp" +#include "serif/solver/normalization/metrics.hpp" +#include "serif/solver/normalization/dimensionless.hpp" +#include "serif/utils/dependencies/binding.hpp" + +namespace serif::solver::normalization +{ + struct NormalizationPrescriptionTag + { + }; + + template + concept NormalizationPrescription = std::derived_from, NormalizationPrescriptionTag>; + + /** @brief Leave the declared coordinates unchanged, including their physical units. */ + struct IdentityCoordinate final + { + }; + + /** @brief Measure a physical quantity with the named reference Gram diagonal. */ + template + struct PhysicalRieszCoordinate final + { + using Quantity = QuantityT; + using Metric = MetricKeyT; + }; + + /** + * @brief Measure coordinates using dimensionless physical fields and reference geometry. + * @note The surface convention applies to shape lengths and their integrated force residuals. Custom metrics must select their convention. + */ + template ::value> + struct DimensionlessRieszCoordinate final { + using Quantity = QuantityT; + using Metric = MetricKeyT; + static constexpr ReferenceGeometry geometry = Geometry; + }; + + /** @brief Replace selected volume residual rows by point residuals with their own physical quantity. */ + template + struct WithPointRows final + { + using Selection = SelectionKeyT; + }; + + /** @brief Identify the exact prescription supplying positive factors; this never claims identity scaling. */ + template + struct RuntimePreparedCoordinate final + { + using Prescription = std::remove_cvref_t; + }; + + /** @brief Independent term-to-method and row-to-method maps; their order does not define the packed layout. */ + template + struct NormalizationPlan final + { + using Values = ValueRulesT; + using Residuals = ResidualRulesT; + }; + + namespace detail + { + template + struct RuleKeys; + + template + struct RuleKeys> + { + using Type = utils::dependencies::KeyList; + }; + + /** @brief Compare memberships in both directions without deduplicating ownership declarations. */ + template + [[nodiscard]] consteval bool exact_coordinates(utils::dependencies::KeyList, utils::dependencies::KeyList) + { + return sizeof...(Declared) == sizeof...(Required) && (utils::dependencies::ContainsKey> && ...); + } + + /** @brief Check physics metadata when present; explicit low-level residual declarations remain their own authority. */ + template + [[nodiscard]] consteval bool quantity_matches() + { + if constexpr (Kind == discretization::blocks::RepresentedValueKind::value && requires + { + typename KeyT::Quantity; + }) + { + return std::same_as; + } + else if constexpr (Kind == discretization::blocks::RepresentedValueKind::residual && requires + { + typename KeyT::ResidualQuantity; + }) + { + return std::same_as; + } + else + { + return true; + } + } + + template + struct MethodSupports : std::false_type + { + }; + + template + struct MethodSupports : std::true_type + { + }; + + template + struct MethodSupports, KeyT, Kind> : std::bool_constant && quantity_matches()> + { + }; + + template + struct MethodSupports, PointQuantityT, SelectionT>, KeyT, discretization::blocks::RepresentedValueKind::residual> : MethodSupports, KeyT, discretization::blocks::RepresentedValueKind::residual> + { + }; + + template + struct MethodSupports, KeyT, Kind> : std::bool_constant< + MetricFor && quantity_matches() && + (Geometry == ReferenceGeometry::volume || Geometry == ReferenceGeometry::surface || Geometry == ReferenceGeometry::algebraic) && + (!std::same_as || Geometry == ReferenceGeometry::algebraic)> { + }; + + template + struct MethodSupports, PointQuantityT, SelectionT>, KeyT, discretization::blocks::RepresentedValueKind::residual> : + MethodSupports, KeyT, discretization::blocks::RepresentedValueKind::residual> { + }; + + template + struct MethodSupports, KeyT, Kind> : std::true_type + { + }; + + template + struct ValidRules : std::false_type + { + }; + + template + struct ValidRules, Kind> : std::bool_constant> && (MethodSupports::value && ...)> + { + }; + + template + inline constexpr bool valid_plan = false; + + template + inline constexpr bool valid_plan> = ValidRules::value && ValidRules::value; + } + + template + concept NormalizationPlanType = detail::valid_plan>; + + template + concept CompleteNormalizationFor = NormalizationPlanType && utils::dependencies::IsKeyList && utils::dependencies::IsKeyList && + detail::exact_coordinates(typename detail::RuleKeys::Values>::Type{}, std::remove_cvref_t{}) && + detail::exact_coordinates(typename detail::RuleKeys::Residuals>::Type{}, std::remove_cvref_t{}); +} diff --git a/src/include/serif/solver/normalization/point_rows.hpp b/src/include/serif/solver/normalization/point_rows.hpp new file mode 100644 index 00000000..95683347 --- /dev/null +++ b/src/include/serif/solver/normalization/point_rows.hpp @@ -0,0 +1,15 @@ +#pragma once + +#include +#include +#include +#include + +namespace serif::solver::normalization +{ + /** + * @brief Measure volume and replacement point residuals in their independently specified physical units. + * @note Only replaced rows may have a zero or otherwise unusable volume metric. No bulk factor is evaluated for them. + */ + [[nodiscard]] std::expected, std::string> make_point_replaced_residual_factors(std::span volume_metric, double volume_scale, std::span selected_rows, double point_scale); +} diff --git a/src/include/serif/solver/normalization/scales.hpp b/src/include/serif/solver/normalization/scales.hpp new file mode 100644 index 00000000..7a3fd4c5 --- /dev/null +++ b/src/include/serif/solver/normalization/scales.hpp @@ -0,0 +1,115 @@ +#pragma once + +#include +#include +#include +#include + +#include "serif/dimensions/dimensions.hpp" + +namespace serif::solver::normalization +{ + /** @brief Fixed reference scales used throughout a solve, independent of its candidate state. */ + struct StellarReferenceScales final + { + dimensions::QuantityValue mass; + dimensions::QuantityValue radius; + dimensions::QuantityValue gravitational_constant; + dimensions::QuantityValue density; + dimensions::QuantityValue acceleration; + dimensions::QuantityValue inverse_time_squared; + dimensions::QuantityValue specific_energy; + dimensions::QuantityValue pressure; + dimensions::QuantityValue angular_velocity; + dimensions::QuantityValue angular_momentum; + dimensions::QuantityValue force; + + /** @brief Compare reference values for dependency snapshot invalidation. */ + friend bool operator==(const StellarReferenceScales &, const StellarReferenceScales &) = default; + }; + + /** @brief Derive gravitational stellar scales from a finite, positive reference mass, radius, and G. */ + [[nodiscard]] std::expected prepare_stellar_reference_scales( + dimensions::QuantityValue mass, + dimensions::QuantityValue radius, + dimensions::QuantityValue gravitational_constant = dimensions::QuantityValue{1.0}); + + /** @brief Quantities whose characteristic scales follow from reference mass, radius, and G alone. */ + template + concept StellarScaleQuantity = std::same_as || std::same_as || + std::same_as || std::same_as || + std::same_as || std::same_as || + std::same_as || std::same_as || + std::same_as || + std::same_as || std::same_as || + std::same_as || std::same_as || std::same_as || + std::same_as || std::same_as || + std::same_as || std::same_as || std::same_as; + + /** @brief Select a typed scale; temperature and other unrelated quantities require an explicit prescription. */ + template + [[nodiscard]] dimensions::QuantityValue characteristic_scale(const StellarReferenceScales &scales, Quantity) + { + if constexpr (std::same_as) + { + return dimensions::QuantityValue{1.0}; + } + else if constexpr (std::same_as) + { + return scales.mass; + } + else if constexpr (std::same_as) + { + return scales.radius; + } + else if constexpr (std::same_as) + { + return scales.density; + } + else if constexpr (std::same_as) + { + return scales.acceleration; + } + else if constexpr (std::same_as || std::same_as) + { + return dimensions::QuantityValue{scales.inverse_time_squared.value()}; + } + else if constexpr (std::same_as || std::same_as || + std::same_as || std::same_as) + { + return dimensions::QuantityValue{scales.specific_energy.value()}; + } + else if constexpr (std::same_as) + { + return scales.pressure; + } + else if constexpr (std::same_as) + { + return dimensions::QuantityValue{std::sqrt(scales.specific_energy.value())}; + } + else if constexpr (std::same_as) + { + return dimensions::QuantityValue{scales.density.value() * scales.angular_velocity.value()}; + } + else if constexpr (std::same_as) + { + return dimensions::QuantityValue{scales.density.value() * scales.acceleration.value()}; + } + else if constexpr (std::same_as) + { + return dimensions::QuantityValue{scales.pressure.value() / (scales.radius.value() * scales.radius.value())}; + } + else if constexpr (std::same_as) + { + return scales.angular_velocity; + } + else if constexpr (std::same_as) + { + return scales.angular_momentum; + } + else + { + return scales.force; + } + } +} diff --git a/src/lib/serif/solver/normalization/diagonal.cpp b/src/lib/serif/solver/normalization/diagonal.cpp new file mode 100644 index 00000000..9340fa9d --- /dev/null +++ b/src/lib/serif/solver/normalization/diagonal.cpp @@ -0,0 +1,89 @@ +#include "serif/solver/normalization/diagonal.hpp" + +#include +#include +#include + +#include "serif/solver/normalization/detail/validation.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::solver::normalization +{ + namespace { + /** @brief Measure a physical vector in its prepared diagonal metric. */ + [[nodiscard]] double normalization_norm_squared(std::span factors, std::span physical) { + detail::require_normalization_size(physical.size(), factors.size()); + double norm_squared = 0.0; + for (std::size_t index = 0; index < physical.size(); ++index) + { + const double normalized = factors[index] * physical[index]; + norm_squared += normalized * normalized; + } + return norm_squared; + } + } + + void detail::require_normalization_size(std::size_t received, std::size_t expected) { + using namespace utils::errors; + + serif_ensure(received == expected, "Normalization received a vector of length {} but expected {}. Supply a vector in the prepared coordinate layout.", received, expected); + } + + void detail::apply_normalization(const std::span factors, const std::span input, std::span output, const bool inverse) + { + using namespace utils::errors; + + require_normalization_size(input.size(), factors.size()); + require_normalization_size(output.size(), factors.size()); + if (input.data() != output.data() && !input.empty()) + { + const std::less before; + const bool input_ends_before_output = !before(output.data(), input.data() + input.size()); + const bool output_ends_before_input = !before(input.data(), output.data() + output.size()); + serif_ensure(input_ends_before_output || output_ends_before_input, "Normalization input and output overlap partially. Use the same span for an in-place conversion or provide disjoint storage."); + } + for (std::size_t index = 0; index < input.size(); ++index) + { + output[index] = inverse ? input[index] / factors[index] : factors[index] * input[index]; + } + } + + DiagonalNormalization::DiagonalNormalization(std::vector state, std::vector residual) : m_state_factors(std::move(state)), m_residual_factors(std::move(residual)) {} + + void DiagonalNormalization::normalize_state(const std::span physical, const std::span normalized) const { + detail::apply_normalization(m_state_factors, physical, normalized, false); + } + + void DiagonalNormalization::denormalize_state(const std::span normalized, const std::span physical) const { + detail::apply_normalization(m_state_factors, normalized, physical, true); + } + + void DiagonalNormalization::normalize_residual(const std::span physical, const std::span normalized) const { + + detail::apply_normalization(m_residual_factors, physical, normalized, false); + } + + void DiagonalNormalization::denormalize_residual(const std::span normalized, const std::span physical) const { + detail::apply_normalization(m_residual_factors, normalized, physical, true); + } + + double DiagonalNormalization::state_norm_squared(const std::span physical) const { + return normalization_norm_squared(m_state_factors, physical); + } + + double DiagonalNormalization::residual_norm_squared(const std::span physical) const { + return normalization_norm_squared(m_residual_factors, physical); + } + + std::expected make_diagonal_normalization(std::vector state, std::vector residual) { + for (const std::span factors : {std::span{state}, std::span{residual}}) + { + const auto valid = detail::validate_factors(factors); + if (!valid) + { + return std::unexpected(valid.error()); + } + } + return DiagonalNormalization{std::move(state), std::move(residual)}; + } +} diff --git a/src/lib/serif/solver/normalization/dimensionless.cpp b/src/lib/serif/solver/normalization/dimensionless.cpp new file mode 100644 index 00000000..b5411532 --- /dev/null +++ b/src/lib/serif/solver/normalization/dimensionless.cpp @@ -0,0 +1,46 @@ +#include "serif/solver/normalization/dimensionless.hpp" + +#include +#include + +#include "serif/utils/error/errors.hpp" + +namespace serif::solver::normalization::detail { + std::expected dimensionless_riesz_scale( + discretization::blocks::RepresentedValueKind kind, ReferenceGeometry geometry, double physical_scale, double reference_radius) { + using namespace utils::errors; + using Kind = discretization::blocks::RepresentedValueKind; + + if (!std::isfinite(physical_scale) || physical_scale <= 0.0 || !std::isfinite(reference_radius) || reference_radius <= 0.0) { + return std::unexpected(format_error_message( + "Dimensionless normalization needs finite positive physical and radius scales (received {} and {}). Choose a valid reference star.", physical_scale, reference_radius)); + } + if (kind != Kind::value && kind != Kind::residual) { + return std::unexpected(format_error_message( + "Dimensionless normalization requires a value or residual coordinate kind. Supply a supported RepresentedValueKind.")); + } + + long double scale = physical_scale; + switch (geometry) { + case ReferenceGeometry::volume: + // Raw volume Gram matrices retain the physical basis, including the RT Piola transform. + scale *= static_cast(reference_radius) * std::sqrt(static_cast(reference_radius)); + break; + case ReferenceGeometry::surface: + // Shape values are lengths; their dual coordinates are integrated forces, not force per area. + scale *= kind == Kind::value ? static_cast(reference_radius) : 1.0L / reference_radius; + break; + case ReferenceGeometry::algebraic: + break; + default: + return std::unexpected(format_error_message( + "Dimensionless normalization received an unknown reference geometry. Select volume, surface, or algebraic coordinates.")); + } + const double result = static_cast(scale); + if (!std::isfinite(result) || result <= 0.0) { + return std::unexpected(format_error_message( + "The dimensionless normalization scale is not representable as a finite positive double. Change the reference units.")); + } + return result; + } +} diff --git a/src/lib/serif/solver/normalization/packing.cpp b/src/lib/serif/solver/normalization/packing.cpp new file mode 100644 index 00000000..9e8af40d --- /dev/null +++ b/src/lib/serif/solver/normalization/packing.cpp @@ -0,0 +1,52 @@ +#include "serif/solver/normalization/detail/packing.hpp" + +#include +#include +#include + +#include "serif/utils/error/errors.hpp" + +namespace serif::solver::normalization::detail +{ + namespace + { + /** @brief Preserve the declared block order while checking the total before allocating packed storage. */ + [[nodiscard]] std::expected, std::string> pack_factors(std::span> coordinates) + { + using namespace utils::errors; + + std::size_t count = 0; + for (const std::span factors : coordinates) + { + if (factors.size() > std::vector{}.max_size() - count) + { + return std::unexpected(format_error_message("The packed normalization layout exceeds the supported vector size. Reduce its coordinate counts.")); + } + count += factors.size(); + } + + std::vector packed; + packed.reserve(count); + for (const std::span factors : coordinates) + { + packed.insert(packed.end(), factors.begin(), factors.end()); + } + return packed; + } + } + + std::expected pack_normalization(std::span> states, std::span> residuals) + { + auto state_factors = pack_factors(states); + if (!state_factors) + { + return std::unexpected(state_factors.error()); + } + auto residual_factors = pack_factors(residuals); + if (!residual_factors) + { + return std::unexpected(residual_factors.error()); + } + return make_diagonal_normalization(std::move(*state_factors), std::move(*residual_factors)); + } +} diff --git a/src/lib/serif/solver/normalization/riesz_factors.cpp b/src/lib/serif/solver/normalization/riesz_factors.cpp new file mode 100644 index 00000000..35daed52 --- /dev/null +++ b/src/lib/serif/solver/normalization/riesz_factors.cpp @@ -0,0 +1,135 @@ +#include "serif/solver/normalization/point_rows.hpp" + +#include + +#include "serif/solver/normalization/diagonal.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::solver::normalization +{ + namespace + { + using Kind = discretization::blocks::RepresentedValueKind; + + /** @brief Validate units independently of whether any coordinate uses their reciprocal. */ + [[nodiscard]] std::expected validate_scale(double physical_scale) + { + using namespace utils::errors; + + if (!std::isfinite(physical_scale) || physical_scale <= 0.0) + { + return std::unexpected(format_error_message("A physical normalization scale must be finite and positive (received {}). Choose a valid reference scale.", physical_scale)); + } + return {}; + } + + /** @brief Validate coordinate-wide settings once before scalar or bulk Riesz preparation. */ + [[nodiscard]] std::expected validate_riesz_parameters(Kind kind, double physical_scale) + { + using namespace utils::errors; + + const auto valid_scale = validate_scale(physical_scale); + if (!valid_scale) + { + return valid_scale; + } + if (kind != Kind::value && kind != Kind::residual) + { + return std::unexpected(format_error_message("Riesz factors require a value or residual coordinate kind. Supply a supported RepresentedValueKind.")); + } + return {}; + } + + /** @brief Compute one metric-dependent factor after the coordinate's role and physical scale are validated. */ + [[nodiscard]] std::expected riesz_factor(bool dual, double metric, double physical_scale, std::size_t coordinate) + { + using namespace utils::errors; + + if (!std::isfinite(metric) || metric <= 0.0) + { + return std::unexpected(format_error_message("Normalization metric entry {} must be finite and positive (received {}). Use a valid reference Gram diagonal.", coordinate, metric)); + } + const double root_metric = std::sqrt(metric); + // Inverting the root first avoids overflow of scale * root_metric when its reciprocal is representable. + const double factor = dual ? (1.0 / root_metric) / physical_scale : root_metric / physical_scale; + if (!std::isfinite(factor) || factor <= 0.0) + { + return std::unexpected(format_error_message("Normalization factor {} is not representable as a finite positive double for metric {} and scale {}. Change the reference units.", coordinate, metric, physical_scale)); + } + return factor; + } + } + + std::expected make_riesz_factor(Kind kind, double metric, double physical_scale) + { + const auto valid = validate_riesz_parameters(kind, physical_scale); + if (!valid) + { + return std::unexpected(valid.error()); + } + return riesz_factor(kind == Kind::residual, metric, physical_scale, 0); + } + + std::expected, std::string> make_riesz_factors(Kind kind, std::span metric, double physical_scale) + { + const auto valid = validate_riesz_parameters(kind, physical_scale); + if (!valid) + { + return std::unexpected(valid.error()); + } + const bool dual = kind == Kind::residual; + std::vector factors(metric.size()); + for (std::size_t index = 0; index < metric.size(); ++index) + { + const auto factor = riesz_factor(dual, metric[index], physical_scale, index); + if (!factor) + { + return std::unexpected(factor.error()); + } + factors[index] = *factor; + } + return factors; + } + + std::expected, std::string> make_point_replaced_residual_factors(std::span volume_metric, double volume_scale, std::span selected_rows, double point_scale) + { + using namespace utils::errors; + + for (const double scale : {volume_scale, point_scale}) + { + const auto valid = validate_scale(scale); + if (!valid) + { + return std::unexpected(valid.error()); + } + } + + std::vector replaced(volume_metric.size(), false); + for (const int row : selected_rows) + { + if (row < 0 || static_cast(row) >= volume_metric.size()) + { + return std::unexpected(format_error_message("Replacement normalization row {} lies outside its residual coordinate of size {}. Use the physical equation's reduced row selection.", row, volume_metric.size())); + } + if (replaced[row]) + { + return std::unexpected(format_error_message("Residual row {} was selected more than once for point normalization. Supply disjoint replacement rows.", row)); + } + replaced[row] = true; + } + + std::vector factors(volume_metric.size()); + for (std::size_t row = 0; row < factors.size(); ++row) + { + const double metric = replaced[row] ? 1.0 : volume_metric[row]; + const double scale = replaced[row] ? point_scale : volume_scale; + const auto factor = riesz_factor(true, metric, scale, row); + if (!factor) + { + return std::unexpected(factor.error()); + } + factors[row] = *factor; + } + return factors; + } +} diff --git a/src/lib/serif/solver/normalization/scales.cpp b/src/lib/serif/solver/normalization/scales.cpp new file mode 100644 index 00000000..7b5baf4d --- /dev/null +++ b/src/lib/serif/solver/normalization/scales.cpp @@ -0,0 +1,62 @@ +#include "serif/solver/normalization/scales.hpp" + +#include + +#include "serif/utils/error/errors.hpp" + +namespace serif::solver::normalization +{ + /** @brief Derive gravitational stellar scales from a finite, positive reference mass, radius, and G. */ + std::expected prepare_stellar_reference_scales( + dimensions::QuantityValue mass, + dimensions::QuantityValue radius, + dimensions::QuantityValue gravitational_constant) + { + using namespace utils::errors; + + for (const double value : {mass.value(), radius.value(), gravitational_constant.value()}) + { + if (!std::isfinite(value) || value <= 0.0) + { + return std::unexpected(format_error_message("Stellar normalization requires finite, positive reference mass, radius, and gravitational constant. Correct the reference values.")); + } + } + + const double reference_mass = mass.value(); + const double reference_radius = radius.value(); + const double gravity = gravitational_constant.value(); + const double radius_squared = reference_radius * reference_radius; + const double radius_cubed = radius_squared * reference_radius; + const double density = reference_mass / radius_cubed; + const double acceleration = gravity * reference_mass / radius_squared; + const double inverse_time_squared = gravity * reference_mass / radius_cubed; + const double specific_energy = gravity * reference_mass / reference_radius; + const double pressure = density * specific_energy; + const double angular_velocity = std::sqrt(inverse_time_squared); + const double angular_momentum = reference_mass * std::sqrt(gravity * reference_mass * reference_radius); + const double force = reference_mass * acceleration; + + for (const double value : {density, acceleration, inverse_time_squared, specific_energy, pressure, angular_velocity, angular_momentum, force}) + { + if (!std::isfinite(value) || value <= 0.0) + { + return std::unexpected(format_error_message("A characteristic stellar normalization scale is not finite and positive; choose reference units in which every derived scale is representable.")); + } + } + + return StellarReferenceScales{ + .mass = mass, + .radius = radius, + .gravitational_constant = gravitational_constant, + .density = dimensions::QuantityValue{density}, + .acceleration = dimensions::QuantityValue{acceleration}, + .inverse_time_squared = dimensions::QuantityValue{inverse_time_squared}, + .specific_energy = dimensions::QuantityValue{specific_energy}, + .pressure = dimensions::QuantityValue{pressure}, + .angular_velocity = dimensions::QuantityValue{angular_velocity}, + .angular_momentum = dimensions::QuantityValue{angular_momentum}, + .force = dimensions::QuantityValue{force} + }; + } + +} diff --git a/src/lib/serif/solver/normalization/validation.cpp b/src/lib/serif/solver/normalization/validation.cpp new file mode 100644 index 00000000..699273a0 --- /dev/null +++ b/src/lib/serif/solver/normalization/validation.cpp @@ -0,0 +1,49 @@ +#include "serif/solver/normalization/detail/validation.hpp" + +#include +#include + +#include "serif/utils/error/errors.hpp" + +namespace serif::solver::normalization::detail +{ + std::expected validate_factors(std::span factors) + { + using namespace utils::errors; + + for (std::size_t index = 0; index < factors.size(); ++index) + { + if (!std::isfinite(factors[index]) || factors[index] <= 0.0) + { + return std::unexpected(format_error_message("Normalization factor {} is not finite and positive (received {}). Supply finite positive factors for the declared coordinate.", index, factors[index])); + } + } + return {}; + } + + std::expected validate_coordinate_size(std::size_t size, bool global_scalar) + { + using namespace utils::errors; + + if (size > std::vector{}.max_size()) + { + return std::unexpected(format_error_message("A normalization coordinate requests {} entries, exceeding supported vector storage. Reduce its state or residual layout count.", size)); + } + if (global_scalar && size != 1) + { + return std::unexpected(format_error_message("A global scalar normalization coordinate requires one entry, but its layout declares {}. Correct its state or residual layout.", size)); + } + return size; + } + + std::expected check_coordinate_size(std::size_t expected, std::size_t received) + { + using namespace utils::errors; + + if (expected != received) + { + return std::unexpected(format_error_message("Normalization received {} entries for a coordinate with {} entries. Use the metric and ordering of its reference representation.", received, expected)); + } + return {}; + } +}