From a1d7d4e67948d47d9ddd7a38320bf57ede306664 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 16 Sep 2026 14:39:55 -0400 Subject: [PATCH] Keep central density as seed data in the n1 sandbox --- ...ENTRAL_CONSTRAINT_HYPOTHESIS_2026-09-16.md | 38 ++++++++++++++++ sandbox.cpp | 10 ++--- tests/user-api/stellar_equilibrium.cpp | 43 +++++++++++++++++++ 3 files changed, 86 insertions(+), 5 deletions(-) create mode 100644 experiments/CENTRAL_CONSTRAINT_HYPOTHESIS_2026-09-16.md diff --git a/experiments/CENTRAL_CONSTRAINT_HYPOTHESIS_2026-09-16.md b/experiments/CENTRAL_CONSTRAINT_HYPOTHESIS_2026-09-16.md new file mode 100644 index 0000000..32c6622 --- /dev/null +++ b/experiments/CENTRAL_CONSTRAINT_HYPOTHESIS_2026-09-16.md @@ -0,0 +1,38 @@ +# Optional central-density constraint + +The n=1 sandbox now supplies central density only to its Lane–Emden seed. Its +equilibrium model fixes the polytropic EOS, surface pressure, total mass, and +angular momentum. The existing public API already supports this combination; +no operator, normalization, preconditioner, or solver API changes are needed. + +`constraint::FixedCentralDensity` remains an explicit, optional model +specification. For n=3 it can select the scale along the nonrotating homology +family. Omitting it remains permitted and may expose a singular or poorly +conditioned problem. Neither the model nor the solver silently inserts a phase +constraint or regularization based on the polytropic index. Seed-only central +density sets an initial condition, not an additional equilibrium equation. + +The hypothesis is that fixing the spherical central density together with K, M, +and nonzero J can distort the interior of an n=1 equilibrium. The implemented +constraint fixes central enthalpy and adds a scalar source to the hydrostatic +row; a small bordered residual can therefore coexist with a larger physical +residual after that source is removed. This is separate from the late Krylov +failure: the initial no-phase J=.01 control retained that failure while reducing +the unbordered residual, with negligible change in virial balance. + +The next matched control is n=1, J=.05 with and without the constraint, using the +same mesh, seed, normalization, preconditioner, and solver settings. Record formal +solve status separately from the retained state's physical quality. Compare +energy and force virial errors against the required 1e-5 milestone, mass and J, +actual angular velocity, the unbordered residual, and location-resolved central +enthalpy/EOS/Bernoulli errors. The phased baseline's sampled maximum EOS mismatch +is about .934% of the density scale despite virial errors near 3e-7; its cause +requires this control. Higher-spin no-phase cases follow if the matched control +is satisfactory. J or angular velocity in reference units is not itself a +measured fraction of critical rotation. + +The added user-API regression covers the previously missing explicit-seed, +no-phase context construction and solver borrowing/release path. Existing tests +already cover the square fixed-J root and its Jacobian column. This lifecycle +test performs no Newton solve and makes no physical-convergence claim. Retain +the existing phased n=3 tests; unphased n=3 convergence is not a requirement. diff --git a/sandbox.cpp b/sandbox.cpp index fa72b14..437ec51 100644 --- a/sandbox.cpp +++ b/sandbox.cpp @@ -99,9 +99,8 @@ int main( // const double polytropicConstant = 0.3639391; // Use for n=3 test std::println("K = {}", polytropicConstant); - const double centralDensity = + const double seedCentralDensity = std::numbers::pi_v * mass / (4.0 * radius * radius * radius); - //const double centralDensity =12.9351146; // use for n=3 test auto stellarModel = model::StellarModel( eos::Polytrope({.n = 1.0, .K = polytropicConstant}), @@ -111,8 +110,7 @@ int main( .Jtotal = dimensions::AngularMomentumValue{angularMomentum}, .axis = {0.0, 0.0, 1.0}, .center = {0.0, 0.0, 0.0} - }), - constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{centralDensity}}) + }) ); auto discretization = equilibrium::makeStellarDiscretization( @@ -126,8 +124,10 @@ int main( auto preconditioner = preconditioning::makePreconditioner(); auto linearSolver = solver::linear::FGMRES({.restartLength = 40, .printLevel = 1}); const auto contextStart = std::chrono::steady_clock::now(); + // This density initializes the seed; the n=1 equilibrium is free to change it. auto context = solver::makeContext( - std::move(stellarModel), std::move(discretization), std::move(preconditioner), std::move(linearSolver) + std::move(stellarModel), std::move(discretization), std::move(preconditioner), std::move(linearSolver), + seed::LaneEmden({.centralDensity = dimensions::DensityValue{seedCentralDensity}}) ); if (rank == 0) { std::cout << "Solver context setup: " diff --git a/tests/user-api/stellar_equilibrium.cpp b/tests/user-api/stellar_equilibrium.cpp index 7a7a205..0072030 100644 --- a/tests/user-api/stellar_equilibrium.cpp +++ b/tests/user-api/stellar_equilibrium.cpp @@ -116,6 +116,49 @@ TEST_CASE( CHECK(context.hasActiveSolver()); } +TEST_CASE( + "An Explicit Lane Emden Seed Builds A Fixed Angular Momentum Context Without Fixing Central Density", + "[user-api][solver][context][seed][fixed-angular-momentum]" +) { + using namespace mean_field; + + auto finiteElements = makeFiniteElements(); + REQUIRE(finiteElements.okay()); + + constexpr double radius = utils::RADIUS; + const double polytropicConstant = 2.0 * utils::G * radius * radius / std::numbers::pi_v; + const double seedCentralDensity = std::numbers::pi_v * utils::MASS / (4.0 * radius * radius * radius); + auto stellarModel = model::StellarModel( + eos::Polytrope({.n = 1.0, .K = polytropicConstant}), + surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}), + integral::FixedTotalMass({.Mtotal = dimensions::MassValue{utils::MASS}}), + integral::FixedAngularMomentum({.Jtotal = dimensions::AngularMomentumValue{0.05}}) + ); + auto discretization = equilibrium::makeStellarDiscretization( + std::move(finiteElements), normalization::PhysicalRieszDiagonal{dimensions::LengthValue{radius}, utils::G} + ); + auto context = solver::makeContext( + std::move(stellarModel), std::move(discretization), preconditioning::makePreconditioner(), + solver::linear::FGMRES({.restartLength = 20, .printLevel = -1}), + seed::LaneEmden({.centralDensity = dimensions::DensityValue{seedCentralDensity}}) + ); + using Problem = typename std::remove_cvref_t::ProblemType; + STATIC_CHECK(Problem::hasFixedAngularMomentum); + STATIC_CHECK_FALSE(Problem::hasFixedCentralDensity); + + REQUIRE(context.isReady()); + CHECK_FALSE(context.hasActiveSolver()); + { + auto equilibriumSolver = solver::make(context, solver::nonlinear::Newton{}); + CHECK(equilibriumSolver.isReady()); + CHECK(context.hasActiveSolver()); + } + CHECK_FALSE(context.hasActiveSolver()); + REQUIRE(context.isReady()); + auto nextSolver = solver::make(context, solver::nonlinear::Newton{}); + CHECK(nextSolver.isReady()); +} + TEST_CASE( "Fixed Angular Momentum User API Generates Rotation And Its Composable Solver Border", "[user-api][fixed-angular-momentum][type]"