MeanField/experiments/LOG_DENSITY_ASSESSMENT_2026-09-19.md

Independent log density: positivity, vacuum, and closure

Branch: research/log-density-assessment-20260919. Parent: 7ed4d23. Scope: mathematical review and a small uniform one-dimensional profile control. No stellar equilibrium, production operator, field order, normalization, or preconditioner was changed. This is separate from the density-only Q2/Q3/Q4 study.

Assessment

Log density is compatible with density remaining an independent physical unknown. It need not make density a function of enthalpy or restrict the eventual EOS to a barotrope. The important choice is what gets exponentiated:

1. Exponentiating the existing nodal coefficients and then evaluating the usual density polynomial does not guarantee a nonnegative field. It is only a coordinate change on positive coefficient vectors. The measured negative density fields already have positive nodal coefficients.
2. Evaluating rho = rho_scale * exp(u_h(x)) pointwise guarantees strictly positive density in exact arithmetic wherever u_h is finite. Floating-point evaluation also needs checks for underflow and overflow. This is a different approximation family from polynomial density. Every density-dependent operator must use it.
3. A finite polynomial u_h cannot represent an exact zero-density vacuum surface with that raw exponential. A vanishing reference factor can encode the surface, but its asymptotic behavior becomes part of the representation.

The small control below shows that a raw exponential can satisfy all the chosen weak density moments while retaining an incorrect positive surface trace. Thus positivity and a small closure residual do not by themselves settle the vacuum or pointwise EOS question. Log density merits a separate formulation hypothesis; it is not a completed remedy for the current zero-pressure polytropic boundary.

Three distinct representations and their derivatives

Let q denote geometry parameters, h enthalpy coefficients, and z the other unknowns. All variations below are at fixed reference coordinates; physical volume and other mapped geometry contributions must also be differentiated.

Logarithms of polynomial nodal coefficients

Suppose the existing density field is

[ \rho_h(\xi)=\sum_j r_j\varphi_j(\xi),\qquad r_j=\rho_\exp(u_j),\quad \rho_>0. ]

For the existing residual vector F(r,q,h,z), keeping its residual coordinates unchanged gives

[ \widetilde F(u,q,h,z)=F(g(u),q,h,z),\qquad J_u=J_r,\operatorname{diag}(r_j). ]

This retains the original roots whose coefficient vectors are positive. The Gauss-Legendre nodal basis has sign-changing basis functions, so positive r_j does not imply nonnegative rho_h. In particular, taking logarithms of the already-positive coefficients at the measured n=1.5 root and applying this coordinate change would reproduce its negative polynomial exactly. Requiring positive coefficients alone does not remove that root.

A pointwise exponential of an independent FE field

Instead choose

[ u_h=\sum_j u_j\psi_j,\qquad \rho(\xi)=\rho_*\exp(u_h(\xi)). ]

Now density is nonpolynomial, but remains independently parameterized by u. For a fixed positive scale,

[ \delta\rho=\rho,\delta u_h. ]

Keeping the current density-moment tests as an explicit discretization choice would give, for the present barotropic closure,

[ R_i=\int_{\Omega(q)}\varphi_i[\rho_*e^{u_h}-\rho_{\rm EOS}(h_h)],dV, ] [ \frac{\partial R_i}{\partial u_j} =\int\varphi_i\rho\psi_j,dV, \qquad \frac{\partial R_i}{\partial h_k} =-\int\varphi_i\rho'_{\rm EOS}(h_h)\chi_k,dV. ]

With equal log trial and closure test bases, the density block is a density-weighted mass matrix. Positive quadrature weights and full column rank make it positive definite for strictly positive sampled density. Its magnitude can nevertheless become very small near vacuum.

This is not the same residual obtained by minimizing the squared pointwise EOS error in the nonlinear trial family. Differentiating that least-squares objective introduces an additional rho in the test factor. That choice need not retain the constant density moment. The test space, residual and derivatives must therefore be chosen explicitly rather than inferred from the word “log.”

For a general EOS, density can remain independent alongside temperature, internal energy and composition, as appropriate to the chosen physics. For example, an auxiliary pressure closure P_aux - P(rho,T,X) has log-density variation -P_rho * rho * delta u_h. No inversion to rho(h) is required by the log representation itself. The current barotropic Bernoulli equation is a separate modeling choice and does not automatically generalize to varying entropy or composition.

In the current formulation, enthalpy is also an independent trial field. The prepared pressure operator rejects nonfinite or negative enthalpy at its own integration points (libmeanfield/impl/operators/prepared_pressure_force.cpp, lines 788–800). Pointwise-positive density alone therefore does not make a Newton candidate admissible. Any coupled log-density experiment must still test the enthalpy directions and pressure-preparation obstruction; the frozen density controls below establish no improvement in that separate guard.

A vanishing reference factor

Let

[ \rho(\xi)=\rho_{\rm ref}(\xi,q)\exp(u_h(\xi)),\qquad \rho_{\rm ref}>0\text{ in the interior},\quad \rho_{\rm ref}=0\text{ on the intended vacuum surface}. ]

Finite u_h now preserves both interior positivity and the prescribed zero trace. Its directional variation is

[ \delta\rho=\rho,\delta u_h+e^{u_h}\delta\rho_{\rm ref}. ]

The second form is usable at a zero of the reference; division by rho_ref to form a logarithmic derivative is not. A reference fixed in mesh coordinates has no explicit shape derivative at fixed reference coordinates, but mapped volume, physical positions and force transformations still change. A reference defined from physical distance or updated with the accepted shape has additional dependencies that must be represented and differentiated.

For a surface profile rho ~ A s^n and a reference rho_ref ~ B s^m,

[ u\sim\log(A/B)+(n-m)\log s. ]

The amplitude is bounded at the surface only when the leading exponent matches. A reference proportional to s is not a uniformly regular log coordinate for the n=1.5 or n=3 profiles. Choosing rho_ref from the current rho_EOS(h) would also introduce enthalpy dependence; it is not necessary for independent density and would need to be distinguished from a reference chosen independently of the live thermodynamic state. A fixed reference based on a model index is a reasonable polytropic experiment, not yet a general-EOS prescription.

Exact vacuum versus a finite-density outer boundary

For the leading vacuum profile rho=A s^n, raw log density is log(A/rho_scale)+n log(s). It tends to negative infinity, and its derivative is n/s. Moving density DOFs to interior Gauss points avoids evaluating the logarithm at zero, but does not give a finite polynomial the correct limit. Exponentiating that polynomial leaves a positive trace on the boundary.

Imposing an additional exact rho=0 condition on a finite raw log field is incompatible; it can lead to coefficients diverging rather than a finite root. The current weak closure can instead converge with a nonzero density trace, because matching volume moments does not enforce that trace. These are two different failure mechanisms, not a proof that the volume moment equations themselves are overdetermined.

A physical outer boundary with specified positive density avoids this particular singularity. Such a boundary needs its own atmospheric/pressure prescription and EOS interpretation. Replacing a vacuum boundary by a numerical density floor would change the problem, not merely regularize coordinates. General EOS models can have different surface asymptotics, so the choice should follow the boundary physics rather than be hidden inside the linear solver.

Uniform surface-cell control

log_density_reference_model.py compares independent approximations to t^n on [0,1], for n=1,1.5,3. All densities below are normalized to the exact inner-edge value, which is one. Physical curvature, tangential variation, gravity, rotation and stellar equilibrium are absent.

The raw log interpolant is the Q2 interpolant of n log(t) at the actual three Gauss-Legendre nodes. It is exponentiated at evaluation points, not interpolated again as density. The conservative exponential control solves only three small moment equations:

[ \int_0^1\varphi_i(t)e^{u_2(t)}dt =\int_0^1\varphi_i(t)t^n dt,\qquad i=0,1,2. ]

The tests are an L2-orthonormal polynomial basis; the right-hand moments are analytical. A three-coordinate Newton calculation fits these moments using a 64-point rule and verifies them with an independently generated 128-point rule. This is a frozen toy fit, not a new stellar solve or a global existence proof.

n	Raw log interpolation: surface density	Raw log interpolation: relative L2 error	Moment-fitted exponential: surface density	Moment-fitted exponential: relative L2 error
1	0.0615032	0.0669374	0.0673707	0.0578060
1.5	0.0152527	0.0924952	0.0229264	0.0517531
3	0.000232645	0.181788	0.00141091	0.0333959

The moment fits take 3,4,5 tiny Newton updates, respectively; verified relative moment defects are at most 2.38e-14. They are strictly positive everywhere in the cell, preserve mass to the reported arithmetic tolerance, and still have nonzero EOS error and an incorrect positive vacuum trace. This gives a concrete example in which positivity and weak closure convergence coexist with a boundary representation error.

For comparison, the exact Q2 density projections have relative L2 errors 0, 0.00952381, and 0.05, respectively. Thus positivity does not imply that the raw log family is a more accurate approximation at the same degree: the moment-fitted exponential improves this norm for the cubic profile, but worsens it for the n=1 and n=1.5 profiles. All these statements concern the one-dimensional profiles, not the actual stellar fields.

A constant shift of the interpolated log field can enforce the exact cell mass: u_new=u+log(M_target/M_current). It leaves the positive surface trace and nonconstant moment errors. The script records this distinct control; it is not used to alter a stellar state. A matched factor t^n exp(0) reproduces the target and its zero surface exactly by construction. That result demonstrates the role of the factor, not the general adequacy of an a priori reference profile.

The smooth exponential controls have 64/128-point moment differences below 3e-16. For the matched fractional factor t^1.5, the same quadratures differ by about 1.60e-10 in the largest orthonormal moment despite the represented function being exactly the target. That is integration error, not density approximation error. The script records it rather than replacing quadrature values by the known integral silently.

Conditioning and conservation

For rho=A Delta^n t^n, using fixed reference polynomial bases on a cell of physical width Delta gives

[ J_u=A\Delta^{n+1}H_n,\qquad (H_n)_{ij}=\int_0^1 t^n\varphi_i\varphi_j,dt. ]

The ordinary density mass matrix in this orthonormal reference basis is Delta I. Relative to it, the log-coordinate response carries the extra factor A Delta^n. Halving width therefore decreases the response scale by 2^-n. This does not by itself worsen the condition number of the isolated cell matrix, but can increase scale disparities across cells and coupled blocks.

The small weighted matrices for the exact reference profiles have condition numbers approximately 7.87298, 16.8523, and 94.5609 for n=1,1.5,3. These are measurements of 3-by-3 matrices for one-dimensional quadratic trial/test spaces, not estimates of the stellar Jacobian's condition number. For the fitted exponential controls the corresponding numbers are 6.65377, 13.7153, and 74.9786.

The first-order identity delta rho/rho=delta u makes a log step a relative density step. That can be useful across large density contrasts, but large log steps create exponential changes. Underflow can produce zero from finite log coefficients, so finite-value checks alone do not ensure positivity. A raw exponential implementation must also detect or reject underflow and nonpositive reconstructed density in the material interior, including underflow when multiplying by its reference scale or factor. Prescribed zero-density traces from a vanishing reference factor remain intentional. Overflow and weak response near vacuum also need explicit checks and meaningful globalization; silently clamping the exponential would change the equations. A vanishing factor removes the singular log profile when it matches the asymptotics, but retains the density-weighted sensitivity near vacuum.

With the explicit constant test retained, a converged weak exponential closure preserves its target cell mass under the closure quadrature. Exponentiation alone does not ensure this; neither does positivity. Mass measured under a different rule, higher moments, inertia and virial balance still require their own consistency checks. A nonlinear change of trial family also does not automatically preserve a discrete variational identity coupling gravity, pressure and shape. Any integration must retain the true full residual and Jacobian tests and measure the physical consequences.

Architecture implications and recommended boundary of this work

The natural internal abstraction is a typed independent-density representation with a physical-value evaluator and exact directional action. A log coordinate has no density units; its positive reference scale does. It should not be passed through APIs that assume its coefficients are physical density coefficients. The same representation must feed gravity source/force, centrifugal force, mass, angular momentum, diagnostics and exported profiles. Their dependency generation, quadrature, cache invalidation and normalization must agree.

The implementation must preserve density's stellar-only support. Physical density and its directional action must remain identically zero in exterior and vacuum domains; the positive reconstruction applies only in the material domain. An inactive log coefficient stored as zero must not be exponentiated into rho_scale outside the star. Domain support must therefore be part of the representation contract shared by operators, diagnostics and exports. This intentional exterior vacuum is distinct from underflow inside the material and from the intended zero trace of a vanishing reference factor at its boundary.

The current registry represents density as an independent scalar L2 field, and the prepared source and force operators read it as a polynomial grid function. The current field relationship catalog supports independent and differential relationships, not a nonlinear physical-value representation. Therefore a pointwise log field requires an explicit extension; setting a flag or exponentiating the accepted coefficient buffer is insufficient. Relevant local contracts are in field_registry.cppm, field_base.cppm, prepared_gravity_source.cpp, prepared_barotropic_closure.cpp, normalization/stellar_equilibrium.cppm, and preconditioning/material_surface.cppm.

For the currently intended exact-vacuum polytrope, the companion density-only order study supplies the polynomial baseline, including cost and pointwise error. If log density remains attractive, the next bounded question is whether a physically justified vanishing reference can provide a regular independent coordinate without importing a restrictive EOS assumption. For a genuinely finite-density outer boundary, raw pointwise log density is a different, potentially simpler candidate. Neither has been implemented or validated against a coupled star here; no logarithmic formulation, n=3 convergence result, speedup, or virial improvement is claimed.

Primary-source context

The analogy to 1D stellar codes is supported by MESA's current star_utils.f90: get_rho_and_lnd_from_xh reconstructs a zone's density with exp(lnd), while store_rho_in_xh stores log(rho). Its documented T(tau) atmosphere treatment sets surface pressure and temperature at a selected optical depth. An ordinary atmosphere match with positive gas pressure therefore avoids requiring log(0) at an exact vacuum interface; this is the relevant inference for our boundary choice, not a claim about every MESA mode. These zonewise logarithmic variables support investigating independent log density, but do not establish positivity between nodes of our high-order FEM reconstruction or resolve its exact zero-density free surface.

Vijaywargiya and Fu, Two Finite Element Approaches for the Porous Medium Equation That Are Positivity Preserving and Energy Stable (2024), formulate a pointwise exponential log-density method and explicitly discuss the infinite logarithm at zero density. Their compact-support implementation treats that issue through a degree-of-freedom deactivation strategy. Its time discretization, mass lumping and PDE-specific guarantees do not establish a corresponding method for steady stellar equilibrium. This is a relevant precedent and vacuum caution, not a proposed algorithm to transplant. Primary paper.

Fu and Xu, High-order space-time finite element methods for the Poisson-Nernst-Planck equations: Positivity and unconditional energy stability, derive positive, conservative schemes through entropy variables and a specific space-time discretization. The conservation and stability results belong to that complete formulation, not to exponentiating coefficients in an arbitrary existing solver. Author preprint.

The representation identities, vacuum asymptotics and one-dimensional controls above are independent derivations rather than imported numerical evidence.

Reproduction

From this checkout:

/opt/homebrew/anaconda3/envs/general/bin/python experiments/log_density_reference_model.py --self-test
/opt/homebrew/anaconda3/envs/general/bin/python experiments/log_density_reference_model.py --output ../../diagnostics/log_density_2026-09-19/reference_v1

The output directory must be fresh. It contains summary CSV/JSON and sampled profile CSV, with script and artifact SHA256 hashes. The retained invocation log is ../../diagnostics/log_density_reference_v1.log. Four tiny test methods pass: Gauss moments, positive coefficient interpolation still admitting undershoot, moment-fit/vacuum behavior, and the weighted Jacobian directional derivative and small spectral identities. The production code was neither built nor run for this assessment.