tboudreaux/MeanField_SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(solver/normalization): Added normalization system

Browse Source 
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.
This commit is contained in:
Emily Boudreaux
2026-11-09 13:30:20 -05:00
parent 9da779417b
commit 24d4c20f82
21 changed files with 1467 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
src/lib/serif/solver/normalization/diagonal.cpp Normal file
View File

@@ -0,0 +1,89 @@
#include "serif/solver/normalization/diagonal.hpp"
#include <functional>
#include <stdexcept>
#include <utility>
#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<const double> factors, std::span<const double> 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<SERiFErrorCode::NORMALIZATION_VECTOR_SIZE_MISMATCH, std::invalid_argument>(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<const double> factors, const std::span<const double> input, std::span<double> 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<const double*> 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<SERiFErrorCode::NORMALIZATION_BUFFER_OVERLAP, std::invalid_argument>(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<double> state, std::vector<double> residual) : m_state_factors(std::move(state)), m_residual_factors(std::move(residual)) {}
void DiagonalNormalization::normalize_state(const std::span<const double> physical, const std::span<double> normalized) const {
detail::apply_normalization(m_state_factors, physical, normalized, false);
}
void DiagonalNormalization::denormalize_state(const std::span<const double> normalized, const std::span<double> physical) const {
detail::apply_normalization(m_state_factors, normalized, physical, true);
}
void DiagonalNormalization::normalize_residual(const std::span<const double> physical, const std::span<double> normalized) const {
detail::apply_normalization(m_residual_factors, physical, normalized, false);
}
void DiagonalNormalization::denormalize_residual(const std::span<const double> normalized, const std::span<double> physical) const {
detail::apply_normalization(m_residual_factors, normalized, physical, true);
}
double DiagonalNormalization::state_norm_squared(const std::span<const double> physical) const {
return normalization_norm_squared(m_state_factors, physical);
}
double DiagonalNormalization::residual_norm_squared(const std::span<const double> physical) const {
return normalization_norm_squared(m_residual_factors, physical);
}
std::expected<DiagonalNormalization, std::string> make_diagonal_normalization(std::vector<double> state, std::vector<double> residual) {
for (const std::span<const double> factors : {std::span<const double>{state}, std::span<const double>{residual}})
{
const auto valid = detail::validate_factors(factors);
if (!valid)
{
return std::unexpected(valid.error());
}
}
return DiagonalNormalization{std::move(state), std::move(residual)};
}
}

46
src/lib/serif/solver/normalization/dimensionless.cpp Normal file
View File

@@ -0,0 +1,46 @@
#include "serif/solver/normalization/dimensionless.hpp"
#include <cmath>
#include <initializer_list>
#include "serif/utils/error/errors.hpp"
namespace serif::solver::normalization::detail {
std::expected<double, std::string> 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<SERiFErrorCode::NORMALIZATION_ERROR>(
"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<SERiFErrorCode::INVALID_ARGUMENT>(
"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<long double>(reference_radius) * std::sqrt(static_cast<long double>(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<long double>(reference_radius) : 1.0L / reference_radius;
break;
case ReferenceGeometry::algebraic:
break;
default:
return std::unexpected(format_error_message<SERiFErrorCode::INVALID_ARGUMENT>(
"Dimensionless normalization received an unknown reference geometry. Select volume, surface, or algebraic coordinates."));
}
const double result = static_cast<double>(scale);
if (!std::isfinite(result) || result <= 0.0) {
return std::unexpected(format_error_message<SERiFErrorCode::NORMALIZATION_ERROR>(
"The dimensionless normalization scale is not representable as a finite positive double. Change the reference units."));
}
return result;
}
}

52
src/lib/serif/solver/normalization/packing.cpp Normal file
View File

@@ -0,0 +1,52 @@
#include "serif/solver/normalization/detail/packing.hpp"
#include <cstddef>
#include <utility>
#include <vector>
#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::vector<double>, std::string> pack_factors(std::span<const std::span<const double>> coordinates)
{
using namespace utils::errors;
std::size_t count = 0;
for (const std::span<const double> factors : coordinates)
{
if (factors.size() > std::vector<double>{}.max_size() - count)
{
return std::unexpected(format_error_message<SERiFErrorCode::INVALID_DOF_COUNT>("The packed normalization layout exceeds the supported vector size. Reduce its coordinate counts."));
}
count += factors.size();
}
std::vector<double> packed;
packed.reserve(count);
for (const std::span<const double> factors : coordinates)
{
packed.insert(packed.end(), factors.begin(), factors.end());
}
return packed;
}
}
std::expected<DiagonalNormalization, std::string> pack_normalization(std::span<const std::span<const double>> states, std::span<const std::span<const double>> 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));
}
}

135
src/lib/serif/solver/normalization/riesz_factors.cpp Normal file
View File

@@ -0,0 +1,135 @@
#include "serif/solver/normalization/point_rows.hpp"
#include <cmath>
#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<void, std::string> validate_scale(double physical_scale)
{
using namespace utils::errors;
if (!std::isfinite(physical_scale) || physical_scale <= 0.0)
{
return std::unexpected(format_error_message<SERiFErrorCode::NORMALIZATION_ERROR>("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<void, std::string> 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<SERiFErrorCode::INVALID_ARGUMENT>("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<double, std::string> 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<SERiFErrorCode::NORMALIZATION_ERROR>("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<SERiFErrorCode::NORMALIZATION_ERROR>("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<double, std::string> 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::vector<double>, std::string> make_riesz_factors(Kind kind, std::span<const double> 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<double> 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::vector<double>, std::string> make_point_replaced_residual_factors(std::span<const double> volume_metric, double volume_scale, std::span<const int> 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<bool> replaced(volume_metric.size(), false);
for (const int row : selected_rows)
{
if (row < 0 || static_cast<std::size_t>(row) >= volume_metric.size())
{
return std::unexpected(format_error_message<SERiFErrorCode::OUT_OF_RANGE_ERROR>("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<SERiFErrorCode::INVALID_ARGUMENT>("Residual row {} was selected more than once for point normalization. Supply disjoint replacement rows.", row));
}
replaced[row] = true;
}
std::vector<double> 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;
}
}

62
src/lib/serif/solver/normalization/scales.cpp Normal file
View File

@@ -0,0 +1,62 @@
#include "serif/solver/normalization/scales.hpp"
#include <cmath>
#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<StellarReferenceScales, std::string> prepare_stellar_reference_scales(
dimensions::QuantityValue<dimensions::Mass> mass,
dimensions::QuantityValue<dimensions::Length> radius,
dimensions::QuantityValue<dimensions::GravitationalConstant> 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<SERiFErrorCode::NORMALIZATION_ERROR>("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<SERiFErrorCode::NORMALIZATION_ERROR>("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<dimensions::MassDensity>{density},
.acceleration = dimensions::QuantityValue<dimensions::Acceleration>{acceleration},
.inverse_time_squared = dimensions::QuantityValue<dimensions::InverseTimeSquared>{inverse_time_squared},
.specific_energy = dimensions::QuantityValue<dimensions::SpecificEnergy>{specific_energy},
.pressure = dimensions::QuantityValue<dimensions::Pressure>{pressure},
.angular_velocity = dimensions::QuantityValue<dimensions::AngularVelocity>{angular_velocity},
.angular_momentum = dimensions::QuantityValue<dimensions::AngularMomentum>{angular_momentum},
.force = dimensions::QuantityValue<dimensions::Force>{force}
};
}
}

49
src/lib/serif/solver/normalization/validation.cpp Normal file
View File

@@ -0,0 +1,49 @@
#include "serif/solver/normalization/detail/validation.hpp"
#include <cmath>
#include <vector>
#include "serif/utils/error/errors.hpp"
namespace serif::solver::normalization::detail
{
std::expected<void, std::string> validate_factors(std::span<const double> 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<SERiFErrorCode::NORMALIZATION_ERROR>("Normalization factor {} is not finite and positive (received {}). Supply finite positive factors for the declared coordinate.", index, factors[index]));
}
}
return {};
}
std::expected<std::size_t, std::string> validate_coordinate_size(std::size_t size, bool global_scalar)
{
using namespace utils::errors;
if (size > std::vector<double>{}.max_size())
{
return std::unexpected(format_error_message<SERiFErrorCode::NORMALIZATION_VECTOR_SIZE_MISMATCH>("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<SERiFErrorCode::NORMALIZATION_VECTOR_SIZE_MISMATCH>("A global scalar normalization coordinate requires one entry, but its layout declares {}. Correct its state or residual layout.", size));
}
return size;
}
std::expected<void, std::string> check_coordinate_size(std::size_t expected, std::size_t received)
{
using namespace utils::errors;
if (expected != received)
{
return std::unexpected(format_error_message<SERiFErrorCode::NORMALIZATION_VECTOR_SIZE_MISMATCH>("Normalization received {} entries for a coordinate with {} entries. Use the metric and ordering of its reference representation.", received, expected));
}
return {};
}
}