Emily Boudreaux
75cc638739
This commit uses global pre allocated work space to dramatically reduce memory usage and allocation time
222 lines
13 KiB
C++
222 lines
13 KiB
C++
#pragma once
|
|
|
|
#include "polytrope_analytic_reference.hpp"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <limits>
|
|
#include <numbers>
|
|
#include <string>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
namespace experiment::polytrope_validation {
|
|
struct AnalyticSelfCheck final {
|
|
std::string name;
|
|
double observed{};
|
|
double expected{};
|
|
double scale{1.0};
|
|
double tolerance{};
|
|
bool passed{};
|
|
|
|
[[nodiscard]] double AbsoluteError() const { return std::abs(observed - expected); }
|
|
[[nodiscard]] double ScaledError() const { return AbsoluteError() / scale; }
|
|
};
|
|
|
|
struct AnalyticSelfCheckReport final {
|
|
std::vector<AnalyticSelfCheck> checks;
|
|
|
|
[[nodiscard]] bool Passed() const {
|
|
return std::all_of(checks.begin(), checks.end(), [](const AnalyticSelfCheck &check) {
|
|
return check.passed;
|
|
});
|
|
}
|
|
};
|
|
|
|
namespace analytic_detail {
|
|
struct RadialIntegrals final {
|
|
long double mass{};
|
|
long double pressure{};
|
|
long double potentialEnergy{};
|
|
long double gradientEnergy{};
|
|
long double fieldEnergy{};
|
|
long double momentOfInertia{};
|
|
};
|
|
|
|
// Independent physical radial integration: no mesh, projection, seed,
|
|
// or production quadrature implementation is involved in these checks.
|
|
inline RadialIntegrals IntegrateReference(const N1Reference &reference) {
|
|
constexpr int intervals = 8192;
|
|
constexpr long double pi = std::numbers::pi_v<long double>;
|
|
const long double spacing = static_cast<long double>(reference.radius) / intervals;
|
|
RadialIntegrals result;
|
|
for (int index = 0; index <= intervals; ++index) {
|
|
const double radius = reference.radius * (static_cast<double>(index) / intervals);
|
|
const auto values = reference.AtRadius(radius);
|
|
const long double r = radius;
|
|
const long double density = values.density;
|
|
const long double gradient = values.radialPotentialGradient;
|
|
const long double weight = (index == 0 || index == intervals) ? 1.0L
|
|
: ((index % 2 == 0) ? 2.0L : 4.0L);
|
|
const long double volumeWeight = weight * 4.0L * pi * r * r;
|
|
result.mass += volumeWeight * density;
|
|
result.pressure += volumeWeight * values.pressure;
|
|
result.potentialEnergy += 0.5L * volumeWeight * density * values.potential;
|
|
result.gradientEnergy -= volumeWeight * density * r * gradient;
|
|
result.fieldEnergy -= volumeWeight * gradient * gradient
|
|
/ (8.0L * pi * reference.gravitationalConstant);
|
|
result.momentOfInertia += (2.0L / 3.0L) * volumeWeight * density * r * r;
|
|
}
|
|
const long double factor = spacing / 3.0L;
|
|
result.mass *= factor;
|
|
result.pressure *= factor;
|
|
result.potentialEnergy *= factor;
|
|
result.gradientEnergy *= factor;
|
|
result.fieldEnergy *= factor;
|
|
result.momentOfInertia *= factor;
|
|
// The gravitational field outside the star is not zero. Its
|
|
// analytic contribution is essential to the field-energy identity.
|
|
result.fieldEnergy -= static_cast<long double>(reference.gravitationalConstant)
|
|
* reference.mass * reference.mass / (2.0L * reference.radius);
|
|
return result;
|
|
}
|
|
} // namespace analytic_detail
|
|
|
|
inline AnalyticSelfCheckReport RunAnalyticSelfChecks() {
|
|
AnalyticSelfCheckReport report;
|
|
auto check = [&](std::string name, const double observed, const double expected,
|
|
const double scale, const double tolerance) {
|
|
const bool passed = std::isfinite(observed) && std::isfinite(expected)
|
|
&& std::isfinite(scale) && scale > 0.0
|
|
&& std::abs(observed - expected) <= tolerance * scale;
|
|
report.checks.push_back({std::move(name), observed, expected, scale, tolerance, passed});
|
|
};
|
|
const std::array<std::pair<std::string, N1Reference>, 2> references{{
|
|
{"unit", {}},
|
|
{"nonunit", {.gravitationalConstant = 2.3, .mass = 3.7, .radius = 1.9}}
|
|
}};
|
|
for (const auto &[label, reference] : references) {
|
|
reference.Validate();
|
|
const double densityScale = reference.CentralDensity();
|
|
const double enthalpyScale = reference.CentralEnthalpy();
|
|
const double gradientScale = enthalpyScale / reference.radius;
|
|
const double energyScale = enthalpyScale * reference.mass;
|
|
const double inertiaScale = reference.mass * reference.radius * reference.radius;
|
|
const auto origin = reference.AtRadius(0.0);
|
|
check(label + ".origin.density", origin.density, densityScale, densityScale, 0.0);
|
|
check(label + ".origin.enthalpy", origin.enthalpy, enthalpyScale, enthalpyScale, 0.0);
|
|
check(label + ".origin.potential", origin.potential, -2.0 * enthalpyScale, enthalpyScale, 0.0);
|
|
check(label + ".origin.enclosed_mass", origin.enclosedMass, 0.0, reference.mass, 0.0);
|
|
check(label + ".origin.gradient", origin.radialPotentialGradient, 0.0, gradientScale, 0.0);
|
|
|
|
constexpr double smallFraction = 1.0e-8;
|
|
const auto nearOrigin = reference.AtRadius(reference.radius * smallFraction);
|
|
constexpr double centralSlope = std::numbers::pi * std::numbers::pi / 3.0;
|
|
check(label + ".origin.mass_cubic_coefficient",
|
|
nearOrigin.enclosedMass / (reference.mass * smallFraction * smallFraction * smallFraction),
|
|
centralSlope, centralSlope, 5.0e-15);
|
|
check(label + ".origin.gradient_linear_coefficient",
|
|
nearOrigin.radialPotentialGradient / (gradientScale * smallFraction),
|
|
centralSlope, centralSlope, 5.0e-15);
|
|
constexpr double tinyFraction = 1.0e-200;
|
|
const auto tinyRadius = reference.AtRadius(reference.radius * tinyFraction);
|
|
check(label + ".origin.gradient_without_mass_underflow_division",
|
|
tinyRadius.radialPotentialGradient / (gradientScale * tinyFraction),
|
|
centralSlope, centralSlope, 5.0e-15);
|
|
|
|
const auto surface = reference.AtRadius(reference.radius);
|
|
check(label + ".surface.density", surface.density, 0.0, densityScale, 0.0);
|
|
check(label + ".surface.enthalpy", surface.enthalpy, 0.0, enthalpyScale, 0.0);
|
|
check(label + ".surface.pressure", surface.pressure, 0.0, enthalpyScale * densityScale, 0.0);
|
|
check(label + ".surface.potential", surface.potential, -enthalpyScale, enthalpyScale, 0.0);
|
|
check(label + ".surface.enclosed_mass", surface.enclosedMass, reference.mass, reference.mass, 0.0);
|
|
check(label + ".surface.gradient", surface.radialPotentialGradient, gradientScale, gradientScale, 0.0);
|
|
const double innerRadius = std::nextafter(reference.radius, 0.0);
|
|
const double outerRadius = std::nextafter(reference.radius, std::numeric_limits<double>::infinity());
|
|
const auto justInside = reference.AtRadius(innerRadius);
|
|
const auto justOutside = reference.AtRadius(outerRadius);
|
|
check(label + ".surface.potential_join", justInside.potential, justOutside.potential, enthalpyScale, 2.0e-15);
|
|
check(label + ".surface.gradient_join", justInside.radialPotentialGradient,
|
|
justOutside.radialPotentialGradient, gradientScale, 3.0e-15);
|
|
check(label + ".surface.theta_linear_coefficient",
|
|
reference.DimensionlessTheta(innerRadius) / ((reference.radius - innerRadius) / reference.radius),
|
|
1.0, 1.0, 2.0e-15);
|
|
check(label + ".surface.positive_density_inside", justInside.density > 0.0 ? 1.0 : 0.0, 1.0, 1.0, 0.0);
|
|
|
|
const auto exterior = reference.AtRadius(2.0 * reference.radius);
|
|
check(label + ".exterior.vacuum_density", exterior.density, 0.0, densityScale, 0.0);
|
|
check(label + ".exterior.point_mass_potential", exterior.potential, -0.5 * enthalpyScale, enthalpyScale, 0.0);
|
|
check(label + ".exterior.point_mass_gradient", exterior.radialPotentialGradient, 0.25 * gradientScale, gradientScale, 0.0);
|
|
check(label + ".normalization.does_not_clip_negative_theta",
|
|
reference.NormalizedPotential(exterior.potential), -0.5, 1.0, 0.0);
|
|
check(label + ".normalization.does_not_clip_positive_potential",
|
|
reference.NormalizedPotential(enthalpyScale), -2.0, 1.0, 0.0);
|
|
|
|
int radialIndex = 0;
|
|
for (const double fraction : {0.0, 0.1, 0.25, 0.5, 0.75, 0.99, 1.0}) {
|
|
const auto values = reference.AtRadius(reference.radius * fraction);
|
|
const std::string prefix = label + ".radial_" + std::to_string(radialIndex++);
|
|
check(prefix + ".enthalpy_eos", values.enthalpy,
|
|
2.0 * reference.PolytropicConstant() * values.density, enthalpyScale, 2.0e-15);
|
|
check(prefix + ".pressure_eos", values.pressure,
|
|
reference.PolytropicConstant() * values.density * values.density,
|
|
enthalpyScale * densityScale, 2.0e-15);
|
|
check(prefix + ".hydrostatic_constant", values.enthalpy + values.potential,
|
|
-enthalpyScale, enthalpyScale, 2.0e-15);
|
|
check(prefix + ".normalized_potential", reference.NormalizedPotential(values.potential),
|
|
reference.DimensionlessTheta(reference.radius * fraction), 1.0, 2.0e-15);
|
|
}
|
|
|
|
radialIndex = 0;
|
|
for (const double fraction : {0.1, 0.25, 0.5, 0.75, 0.9}) {
|
|
const double radius = reference.radius * fraction;
|
|
const double spacing = 1.0e-4 * reference.radius;
|
|
const auto minusTwo = reference.AtRadius(radius - 2.0 * spacing);
|
|
const auto minusOne = reference.AtRadius(radius - spacing);
|
|
const auto plusOne = reference.AtRadius(radius + spacing);
|
|
const auto plusTwo = reference.AtRadius(radius + 2.0 * spacing);
|
|
const double enthalpyDerivative = (minusTwo.enthalpy - 8.0 * minusOne.enthalpy
|
|
+ 8.0 * plusOne.enthalpy - plusTwo.enthalpy) / (12.0 * spacing);
|
|
const double massDerivative = (minusTwo.enclosedMass - 8.0 * minusOne.enclosedMass
|
|
+ 8.0 * plusOne.enclosedMass - plusTwo.enclosedMass) / (12.0 * spacing);
|
|
const auto values = reference.AtRadius(radius);
|
|
const std::string prefix = label + ".derivative_" + std::to_string(radialIndex++);
|
|
check(prefix + ".hydrostatic_balance", enthalpyDerivative + values.radialPotentialGradient,
|
|
0.0, gradientScale, 5.0e-11);
|
|
check(prefix + ".enclosed_mass", massDerivative,
|
|
4.0 * std::numbers::pi * radius * radius * values.density,
|
|
reference.mass / reference.radius, 5.0e-11);
|
|
}
|
|
|
|
const auto integrals = analytic_detail::IntegrateReference(reference);
|
|
check(label + ".integral.mass", static_cast<double>(integrals.mass), reference.mass, reference.mass, 2.0e-12);
|
|
check(label + ".integral.pressure", static_cast<double>(integrals.pressure),
|
|
reference.PressureIntegral(), energyScale, 2.0e-12);
|
|
check(label + ".integral.binding_from_potential", static_cast<double>(integrals.potentialEnergy),
|
|
reference.BindingEnergy(), energyScale, 2.0e-12);
|
|
check(label + ".integral.binding_from_gradient", static_cast<double>(integrals.gradientEnergy),
|
|
reference.BindingEnergy(), energyScale, 2.0e-12);
|
|
check(label + ".integral.binding_from_field_with_exterior", static_cast<double>(integrals.fieldEnergy),
|
|
reference.BindingEnergy(), energyScale, 2.0e-12);
|
|
check(label + ".integral.binding_potential_vs_gradient", static_cast<double>(integrals.potentialEnergy),
|
|
static_cast<double>(integrals.gradientEnergy), energyScale, 2.0e-12);
|
|
check(label + ".integral.scalar_virial_nonrotating",
|
|
static_cast<double>(integrals.gradientEnergy + 3.0L * integrals.pressure),
|
|
0.0, energyScale, 2.0e-12);
|
|
check(label + ".integral.axial_moment_of_inertia", static_cast<double>(integrals.momentOfInertia),
|
|
reference.MomentOfInertia(), inertiaScale, 2.0e-12);
|
|
}
|
|
|
|
bool rejectedNegativeRadius = false;
|
|
try { (void)N1Reference{}.AtRadius(-1.0); }
|
|
catch (const std::invalid_argument &) { rejectedNegativeRadius = true; }
|
|
check("contract.negative_radius_rejected", rejectedNegativeRadius ? 1.0 : 0.0, 1.0, 1.0, 0.0);
|
|
bool rejectedInvalidScale = false;
|
|
try { N1Reference{.gravitationalConstant = -1.0}.Validate(); }
|
|
catch (const std::invalid_argument &) { rejectedInvalidScale = true; }
|
|
check("contract.invalid_reference_rejected", rejectedInvalidScale ? 1.0 : 0.0, 1.0, 1.0, 0.0);
|
|
return report;
|
|
}
|
|
} // namespace experiment::polytrope_validation
|