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

Add explicit mass-compatible rotating n3 physical-family control

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 13:58:08 -04:00
parent 6d423cba46
commit 1d746b028e
3 changed files with 103 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

50
experiments/N3_COMPATIBLE_FAMILY_2026-09-19.md Normal file
View File

@@ -0,0 +1,50 @@
# Explicit rotating n=3 physical-family control
Branch `research/n3-compatible-family-20260919`; parent `6d423cb`.
This changes experimental physical inputs, not solver equations or production
defaults. It must not be presented as a preconditioner speedup comparison.
The driver normally chooses K for a nonrotating mass-one Lane–Emden star.
The sharp gravitational-energy inequality and the steady virial identity
require, for n=3 and a nonnegative exact-EOS density,
beta = T/|W| <= [1 - (M_LE(K)/M)^(2/3)]/2.
At exactly M=M_LE(K), nonzero rotation is impossible. An independent DOP853
Lane–Emden integration with five tolerance/step/start controls gives the
actual driver's K=0.36393905995216302 a Lane–Emden mass
0.9999999918534539. Thus the numerical K is slightly on the permissible
side, but beta cannot exceed approximately 2.71552e-9. The control-to-control
mass spread is 1.95e-14; this is a numerical sensitivity observation, not a
rigorous interval enclosure. The intended moderate rotation cannot be reached
while retaining that essentially critical mass/K combination.
The optional experiment flag `--n3-K-factor f` multiplies K by a positive
factor, with index required to be exactly three when f differs from one.
The default is one. The analytic seed central density is multiplied by
f^(3/2), keeping its radius equal to the unchanged reference radius. Its
mass becomes approximately f^(3/2), while the target mass remains one.
This is a deliberately off-constraint initial guess, not a density rescaling
during Newton. The public EOS and independent density formulation are unchanged.
The seed's Bernoulli scalar continues to use the target mass; its discrepancy
is part of the initial nonlinear residual and is not hidden by this control.
The bounded first case uses f=.99, J=.035, no central border, the existing
gravity-Schur preconditioner and gravity-curvature trial path, 10 updates,
eta=.03 and absolute nonlinear target 1e-12. The necessary beta bound is then
0.00500000269. This does not guarantee existence, convergence, or identify
a fraction of critical rotation. J already breaks the nonrotating n=3
homology, so this case does not add an extra central-density phase condition.
The energy inequality's assumptions are not automatically satisfied by the
discrete independent rho/P(h) fields, approximate gravity, or a central
artificial forcing. This compatibility argument selects a sensible continuum
target; it is not a new density limiter or a substitute for numerical tests.
Primary source: Cheng, Cheng and Lin, equations (1.14)–(1.15) and Theorem 3.1,
[Expanding solutions near unstable Lane–Emden stars](https://arxiv.org/html/2305.10218#S3.SS1).
The displayed beta bound is our consequence of their inequality and
2T+W+3 integral(P)=0. Their G=1 convention matches this diagnostic.
Independent integration evidence is retained in
`../../diagnostics/n3_mass_compatibility_2026-09-19/`.

39
experiments/n3_compatible_family_cases_2026-09-19.json Normal file
View File

@@ -0,0 +1,39 @@
[
{
"id": "n3_J0p035_K0p99",
"arguments": [
"--index",
"3",
"--angular-momentum",
".035",
"--n3-K-factor",
".99",
"--no-central",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--max-newton",
"10",
"--restart",
"160",
"--max-linear",
"300",
"--linear-tol",
".03",
"--rtol",
"0",
"--atol",
"1e-12",
"--model-slices",
"--slice-subdivisions",
"4",
"--capture-iterations",
"0,3,7",
"--strategy",
"curvature",
"--preconditioner-variant",
"gravity-schur"
]
}
]

15
experiments/plateau_diagnostic.cpp
View File

@@ -45,6 +45,7 @@ namespace {
using Clock = std::chrono::steady_clock;
struct Options {
double n = 1., angularMomentum = 0., atol = 0., rtol = 1e-8, linearTolerance = .03;
double n3KFactor = 1.;
int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40;
bool probes = false, noCentral = false;
bool strategyProbe = false, positivityProbe = false, positivityCertificate = false;
@@ -89,6 +90,8 @@ Options Parse(int argc, char **argv) {
};
if (a == "--index")
o.n = std::stod(v());
else if (a == "--n3-K-factor")
o.n3KFactor = std::stod(v());
else if (a == "--angular-momentum")
o.angularMomentum = std::stod(v());
else if (a == "--output")
@@ -238,6 +241,8 @@ Options Parse(int argc, char **argv) {
}
if (o.output.empty() || o.n < 1 || o.n >= 5 || o.maxNewton < 0 || o.maxLinear < 1 || o.restart < 1)
throw std::invalid_argument("Invalid options");
if (!std::isfinite(o.n3KFactor) || o.n3KFactor <= 0. || (o.n3KFactor != 1. && o.n != 3.))
throw std::invalid_argument("--n3-K-factor requires a positive factor and index exactly three");
if ((o.continuation || o.allowUnconvergedRestart) && o.restartSnapshot.empty())
throw std::invalid_argument("Continuation/restart permission requires --restart-snapshot");
if (!o.restartSnapshot.empty() && (o.densityDiagnostic || !o.probeState.empty() ||
@@ -767,6 +772,8 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
<< "\ncontext_gravity_schur_relaxation=" << (contextUsesGravitySchur ? o.schurRelaxation : 0.0)
<< "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation
<< "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n';
metadata << "n3_K_factor=" << o.n3KFactor << "\nseed_analytic_mass_estimate="
<< std::pow(o.n3KFactor, 1.5) << '\n';
metadata << "restart_snapshot=" << o.restartSnapshot.string() << "\ncontinuation=" << o.continuation
<< "\nallow_unconverged_restart=" << o.allowUnconvergedRestart
<< "\nbenchmark_trials=" << o.benchmarkTrials << '\n';
@@ -1393,7 +1400,13 @@ template <class Model> int Setup(Model model, fem::FEM f, const Options &o, cons
deformation::NodalRadialSurface{std::move(center)}, o, scales);
}
int Run(const Options &o) {
const auto scales = Derive(o.n);
auto scales = Derive(o.n);
// At n=3, rotation requires mass above the nonrotating mass at fixed K.
// This explicit physical-family control leaves M=1 while changing K.
// Rescale only the analytic seed's central density to retain R_reference=1;
// its mass is intentionally off-constraint and all fields remain unknowns.
scales.K *= o.n3KFactor;
scales.rho *= std::pow(o.n3KFactor, 1.5);
std::cout << "n=" << o.n << " K=" << scales.K << " rho_c=" << scales.rho << " J=" << o.angularMomentum << std::endl;
utils::Args args;
args.mesh_file = o.mesh;