From 1d746b028ef001782158a0a82d81b40eb933b68a Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 13:58:08 -0400 Subject: [PATCH] Add explicit mass-compatible rotating n3 physical-family control --- .../N3_COMPATIBLE_FAMILY_2026-09-19.md | 50 +++++++++++++++++++ ...n3_compatible_family_cases_2026-09-19.json | 39 +++++++++++++++ experiments/plateau_diagnostic.cpp | 15 +++++- 3 files changed, 103 insertions(+), 1 deletion(-) create mode 100644 experiments/N3_COMPATIBLE_FAMILY_2026-09-19.md create mode 100644 experiments/n3_compatible_family_cases_2026-09-19.json diff --git a/experiments/N3_COMPATIBLE_FAMILY_2026-09-19.md b/experiments/N3_COMPATIBLE_FAMILY_2026-09-19.md new file mode 100644 index 0000000..ab32510 --- /dev/null +++ b/experiments/N3_COMPATIBLE_FAMILY_2026-09-19.md @@ -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/`. diff --git a/experiments/n3_compatible_family_cases_2026-09-19.json b/experiments/n3_compatible_family_cases_2026-09-19.json new file mode 100644 index 0000000..8df41a7 --- /dev/null +++ b/experiments/n3_compatible_family_cases_2026-09-19.json @@ -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" + ] + } +] diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 70262a5..64e9ece 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -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 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;