From 549b1a761dc4506f5e12ffd3d741cd4756c8417f Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 10:35:48 -0400 Subject: [PATCH] Assess independent log density coordinates and exact vacuum --- .../LOG_DENSITY_ASSESSMENT_2026-09-19.md | 334 ++++++++++++++++++ experiments/log_density_reference_model.py | 282 +++++++++++++++ 2 files changed, 616 insertions(+) create mode 100644 experiments/LOG_DENSITY_ASSESSMENT_2026-09-19.md create mode 100644 experiments/log_density_reference_model.py diff --git a/experiments/LOG_DENSITY_ASSESSMENT_2026-09-19.md b/experiments/LOG_DENSITY_ASSESSMENT_2026-09-19.md new file mode 100644 index 0000000..59ac6b2 --- /dev/null +++ b/experiments/LOG_DENSITY_ASSESSMENT_2026-09-19.md @@ -0,0 +1,334 @@ +# 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 wherever `u_h` is finite. 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. No additional +constraint on positive coefficients removes 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. + +### 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 three-dimensional local polynomial-space measurements, 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, overflow and weak response near +vacuum need explicit finite-value 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 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, continue the separately +authorized density-only order study first. Its cost and pointwise-error results +provide the right baseline. 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 + +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](https://link.springer.com/article/10.1007/s10915-024-02642-x). + +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](https://arxiv.org/abs/2105.01163). + +The representation identities, vacuum asymptotics and one-dimensional controls +above are independent derivations rather than imported numerical evidence. + +## Reproduction + +From this checkout: + +```text +/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. diff --git a/experiments/log_density_reference_model.py b/experiments/log_density_reference_model.py new file mode 100644 index 0000000..73c57a6 --- /dev/null +++ b/experiments/log_density_reference_model.py @@ -0,0 +1,282 @@ +#!/usr/bin/env python3 +"""Small independent log-density controls for rho(t)=t**n on [0,1]. + +No stellar solve or FEM library is used. A three-coordinate exponential moment +fit illustrates a changed density approximation, not a production proposal. +Polynomial extrema are analytic; integrals use independently generated Gauss +rules, with 64/128-point comparison. Only stdlib is required. +""" +import argparse +import csv +from fractions import Fraction +import hashlib +import json +import math +from pathlib import Path +import unittest + + +def dot(a, b): + return math.fsum(x * y for x, y in zip(a, b)) + + +def norm(a): + return math.sqrt(dot(a, a)) + + +def gauss(count): + rule = [] + for j in range(1, count + 1): + z = math.cos(math.pi * (j - .25) / (count + .5)) + for _ in range(40): + p0, p1 = 1., z + for degree in range(2, count + 1): + p0, p1 = p1, ((2 * degree - 1) * z * p1 - (degree - 1) * p0) / degree + derivative = count * (z * p1 - p0) / (z * z - 1) + change = p1 / derivative + z -= change + if abs(change) < 2e-16: + break + rule.append(((1 + z) / 2, 1 / ((1 - z * z) * derivative * derivative))) + return sorted(rule) + + +def polynomial(c, t): + return c[0] + t * (c[1] + t * c[2]) + + +def interpolate(values): + separation = math.sqrt(.6) + slope = (values[2] - values[0]) / separation + curvature = (10 / 3) * (values[2] + values[0] - 2 * values[1]) + return [values[1] - slope / 2 + curvature / 4, slope - curvature, curvature] + + +def project(n): + n = Fraction(str(n)) + moments = [1 / (n + j + 1) for j in range(3)] + return [float(sum(a * b for a, b in zip(row, moments))) + for row in ((9, -36, 30), (-36, 192, -180), (30, -180, 180))] + + +def extrema(c): + points = [0., 1.] + if c[2] and 0 < -c[1] / (2 * c[2]) < 1: + points.append(-c[1] / (2 * c[2])) + return min(polynomial(c, t) for t in points), max(polynomial(c, t) for t in points) + + +def basis(t): + return [1., math.sqrt(3) * (2 * t - 1), math.sqrt(5) * (6 * t * t - 6 * t + 1)] + + +def orthogonal_coefficients(c): + return [c[0] + c[1] / 2 + c[2] / 3, (c[1] + c[2]) / (2 * math.sqrt(3)), c[2] / (6 * math.sqrt(5))] + + +def power_coefficients(a): + return [a[0] - math.sqrt(3) * a[1] + math.sqrt(5) * a[2], + 2 * math.sqrt(3) * a[1] - 6 * math.sqrt(5) * a[2], 6 * math.sqrt(5) * a[2]] + + +def target_moments(n): + return [1 / (n + 1), math.sqrt(3) * (2 / (n + 2) - 1 / (n + 1)), + math.sqrt(5) * (6 / (n + 3) - 6 / (n + 2) + 1 / (n + 1))] + + +def moments_and_jacobian(density, rule): + samples = [(weight * density(t), basis(t)) for t, weight in rule] + moments = [math.fsum(weight * phi[i] for weight, phi in samples) for i in range(3)] + jacobian = [[math.fsum(weight * phi[i] * phi[j] for weight, phi in samples) for j in range(3)] for i in range(3)] + return moments, jacobian + + +def solve(matrix, rhs): + a = [row[:] + [value] for row, value in zip(matrix, rhs)] + for column in range(3): + pivot = max(range(column, 3), key=lambda i: abs(a[i][column])) + a[column], a[pivot] = a[pivot], a[column] + if not abs(a[column][column]) > 1e-16: + raise ArithmeticError("Unresolved tiny moment Jacobian") + for row in range(column + 1, 3): + factor = a[row][column] / a[column][column] + for j in range(column, 4): + a[row][j] -= factor * a[column][j] + result = [0.] * 3 + for row in (2, 1, 0): + result[row] = (a[row][3] - dot(a[row][row + 1:3], result[row + 1:])) / a[row][row] + if norm([dot(row, result) - value for row, value in zip(matrix, rhs)]) > 1e-12 * norm(rhs): + raise ArithmeticError("Tiny solve failed independent defect check") + return result + + +def determinant(a): + return (a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) + - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) + + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])) + + +def eigenvalues(a): + # Eigenvalues of the three-by-three symmetric weighted mass matrix itself; + # no Gram-squared surrogate or full stellar conditioning claim is involved. + mean = sum(a[i][i] for i in range(3)) / 3 + centered = [[a[i][j] - (mean if i == j else 0) for j in range(3)] for i in range(3)] + scale = math.sqrt(sum(value * value for row in centered for value in row) / 6) + if scale == 0: + return [mean] * 3 + angle = math.acos(max(-1., min(1., determinant([[x / scale for x in row] for row in centered]) / 2))) / 3 + values = sorted(mean + 2 * scale * math.cos(angle + j * 2 * math.pi / 3) for j in range(3)) + if abs(sum(values) - 3 * mean) > 1e-12 or abs(math.prod(values) - determinant(a)) > 1e-12: + raise ArithmeticError("Tiny spectral invariants failed") + return values + + +def fit_log_moments(n, initial, rule): + target = target_moments(n) + a = orthogonal_coefficients(initial) + for iteration in range(40): + moments, jacobian = moments_and_jacobian(lambda t: math.exp(dot(a, basis(t))), rule) + residual = [value - wanted for value, wanted in zip(moments, target)] + if norm(residual) < 2e-13 * norm(target): + return power_coefficients(a), iteration, norm(residual) / norm(target) + direction = solve(jacobian, [-value for value in residual]) + old_norm = norm(residual) + for backtrack in range(20): + step = 2. ** -backtrack + trial = [value + step * change for value, change in zip(a, direction)] + try: + trial_moments, _ = moments_and_jacobian(lambda t: math.exp(dot(trial, basis(t))), rule) + new_norm = norm([value - wanted for value, wanted in zip(trial_moments, target)]) + except OverflowError: + continue + if new_norm < (1 - 1e-4 * step) * old_norm: + a = trial + break + else: + raise ArithmeticError("Tiny moment fit failed line search") + raise ArithmeticError("Tiny moment fit did not converge") + + +def controls(n): + nodes = [t for t, _ in gauss(3)] + direct = interpolate([t ** n for t in nodes]) + logs = interpolate([n * math.log(t) for t in nodes]) + log_mass = moments_and_jacobian(lambda t: math.exp(polynomial(logs, t)), gauss(64))[0][0] + shifted = logs[:] + shifted[0] += math.log((1 / (n + 1)) / log_mass) + fitted, iterations, relative = fit_log_moments(n, logs, gauss(64)) + return [("density_q2_interpolation_or_exp_coefficients", "polynomial", direct, None), + ("density_q2_l2_projection", "polynomial", project(n), None), + ("pointwise_exp_q2_log_interpolation", "exp", logs, None), + ("pointwise_exp_q2_log_mass_shift", "exp", shifted, None), + ("pointwise_exp_q2_log_moment_fit", "exp", fitted, + {"iterations": iterations, "relative_moment_defect_64": relative}), + ("matched_factor_t_power_n_exp_zero", "factor", [0., 0., 0.], None)] + + +def evaluate_case(n, name, kind, coefficients, fit): + if kind == "polynomial": + density = lambda t: polynomial(coefficients, t) + minimum, maximum = extrema(coefficients) + elif kind == "exp": + density = lambda t: math.exp(polynomial(coefficients, t)) + lo, hi = extrema(coefficients) + minimum, maximum = math.exp(lo), math.exp(hi) + else: + density = lambda t: t ** n + minimum, maximum = 0., 1. + moments, jacobian = moments_and_jacobian(density, gauss(128)) + coarse_moments, _ = moments_and_jacobian(density, gauss(64)) + residual = [a - b for a, b in zip(moments, target_moments(n))] + error = math.sqrt(math.fsum(w * (density(t) - t ** n) ** 2 for t, w in gauss(128))) + eig = eigenvalues(jacobian) if kind in ("exp", "factor") else None + return dict(n=n, name=name, kind=kind, coefficients=coefficients, + boundary_density=density(0.), minimum_density=minimum, maximum_density=maximum, + cell_mass=moments[0], relative_mass_error=moments[0] * (n + 1) - 1, + relative_l2_error=error * math.sqrt(2 * n + 1), + moment_defect=residual, relative_moment_defect=norm(residual) / norm(target_moments(n)), + moment_quadrature_64_vs_128_maximum=max(abs(a - b) for a, b in zip(moments, coarse_moments)), + log_coordinate_jacobian_eigenvalues=eig, + log_coordinate_jacobian_condition=eig[-1] / eig[0] if eig else None, + fit=fit) + + +class Tests(unittest.TestCase): + def test_gauss(self): + for degree in range(12): + self.assertAlmostEqual(math.fsum(w * t ** degree for t, w in gauss(8)), 1 / (degree + 1), places=14) + + def test_log_coefficients_do_not_fix_polynomial(self): + for n in (1.5, 3.): + nodes = [t for t, _ in gauss(3)] + a = interpolate([math.exp(n * math.log(t)) for t in nodes]) + b = interpolate([t ** n for t in nodes]) + self.assertLess(max(abs(x - y) for x, y in zip(a, b)), 2e-15) + self.assertLess(extrema(a)[0], 0) + + def test_fit_and_vacuum(self): + for n in (1., 1.5, 3.): + results = [evaluate_case(n, *case) for case in controls(n)] + fitted = next(row for row in results if row["name"].endswith("moment_fit")) + self.assertLess(fitted["relative_moment_defect"], 3e-13) + self.assertGreater(fitted["boundary_density"], 0) + self.assertGreater(fitted["relative_l2_error"], 1e-3) + exact = results[-1] + self.assertEqual(exact["boundary_density"], 0) + self.assertEqual(exact["relative_l2_error"], 0) + + def test_jacobian(self): + a, direction, step = [-1., .2, -.1], [.3, -.1, .2], 1e-5 + rule = gauss(64) + _, jacobian = moments_and_jacobian(lambda t: math.exp(dot(a, basis(t))), rule) + plus, _ = moments_and_jacobian(lambda t: math.exp(dot([x + step * y for x, y in zip(a, direction)], basis(t))), rule) + minus, _ = moments_and_jacobian(lambda t: math.exp(dot([x - step * y for x, y in zip(a, direction)], basis(t))), rule) + error = [(x - y) / (2 * step) - dot(row, direction) for x, y, row in zip(plus, minus, jacobian)] + self.assertLess(norm(error), 2e-10) + self.assertEqual(eigenvalues([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), [1., 1., 1.]) + self.assertLess(max(abs(x - y) for x, y in zip(eigenvalues([[2., 1., 0.], [1., 2., 0.], [0., 0., 4.]]), [1., 3., 4.])), 2e-14) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path) + parser.add_argument("--self-test", action="store_true") + args = parser.parse_args() + if args.self_test: + result = unittest.TextTestRunner(verbosity=2).run(unittest.defaultTestLoader.loadTestsFromTestCase(Tests)) + if not result.wasSuccessful(): + raise SystemExit(1) + if args.output is None: + return + cases = [evaluate_case(n, *case) for n in (1., 1.5, 3.) for case in controls(n)] + data = dict(schema="log_density_reference_v1", producer_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), + scope="Uniform 1D frozen t^n profile; no FEM, coupled star solve, virial test, or general-EOS validation.", + eigenvalue_scope="Three-dimensional log-coordinate weighted mass matrix in L2-orthonormal polynomial basis on [0,1].", + scaling="For rho=A*Delta^n*t^n, local reference-basis Jacobian=A*Delta^(n+1)*H_n; relative to ordinary mass the scale is A*Delta^n.", + cases=cases) + if args.output is None: + print(json.dumps(data, indent=2, allow_nan=False)) + return + args.output.mkdir(parents=True, exist_ok=False) + with (args.output / "summary.csv").open("w", newline="") as stream: + writer = csv.DictWriter(stream, fieldnames=list(cases[0])) + writer.writeheader() + for row in cases: + writer.writerow({key: json.dumps(value) if isinstance(value, (list, dict)) else value for key, value in row.items()}) + with (args.output / "profiles.csv").open("w", newline="") as stream: + writer = csv.writer(stream) + writer.writerow(["n", "method", "t", "target", "density"]) + for row in cases: + for j in range(401): + t = j / 400 + value = (t ** row["n"] if row["kind"] == "factor" else + math.exp(polynomial(row["coefficients"], t)) if row["kind"] == "exp" else + polynomial(row["coefficients"], t)) + writer.writerow([row["n"], row["name"], t, t ** row["n"], value]) + data["artifacts_sha256"] = {path.name: hashlib.sha256(path.read_bytes()).hexdigest() for path in sorted(args.output.iterdir())} + (args.output / "summary.json").write_text(json.dumps(data, indent=2, allow_nan=False) + "\n") + print("Wrote", args.output) + + +if __name__ == "__main__": + main()