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

Assess independent log density coordinates and exact vacuum

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 10:35:48 -04:00
parent 6950dc68af
commit 549b1a761d
2 changed files with 616 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

334
experiments/LOG_DENSITY_ASSESSMENT_2026-09-19.md Normal file
View File

@@ -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.

282
experiments/log_density_reference_model.py Normal file
View File

@@ -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()