From 384b59393c4f5d8db1929dfea6ee7f998463faa7 Mon Sep 17 00:00:00 2001 From: xjoepec Date: Sat, 22 Aug 2026 16:12:08 -0400 Subject: [PATCH] Add central-force potential family with conservation certificates Generalize the potential layer from the two hardcoded laws (Harmonic, Coulomb) into an extensible central-force family, without touching the integrator, the Landau-Lifshitz reduction, or the certificate discipline. Potentials (src/blueberry_circus/potentials.py): - CentralPotential base derives the Cartesian force and an analytic force_jacobian (needed by the RR term) from three radial functions U, U', U'' of the softened radius s = sqrt(r^2 + softening^2). - PowerLaw, Yukawa, Morse, AnharmonicOscillator (isotropic Duffing), LennardJones. PowerLaw reduces exactly to Harmonic (p=2) and Coulomb (p=-1), which the tests pin as reduction oracles. Declarative ops (src/blueberry_circus/program.py): - Matching Operations for each new law. Potential collection in compile() now keys on a _PotentialOp marker instead of a hardcoded (Harmonic, Coulomb) tuple, so new families are picked up automatically. Conservation certificates (src/blueberry_circus/conservation.py): - energy_conservation_certificate / angular_momentum_conservation_certificate encode fractional invariant drift against the canonical residual_le_tol rule; an under-resolved or tampered trajectory re-derives FAIL. This cross-checks that each new potential's analytic force is consistent with its potential energy, purely from the integrated dynamics. Tests (tests/test_central_potentials.py, 33 cases): finite-difference checks of force = -grad U and of the Jacobian, Jacobian symmetry, PowerLaw->Harmonic/Coulomb reductions, energy + |L| conservation on closed orbits (all families to <~1e-11), a tight-tolerance FAIL guard, RR stability, and declarative-program integration. Full suite: 167 passed, 2 expected xfails, no regressions. Also adds docs/central_forces.md and examples/central_forces.py. --- docs/central_forces.md | 85 +++++++++++++ examples/central_forces.py | 63 ++++++++++ src/blueberry_circus/__init__.py | 18 ++- src/blueberry_circus/conservation.py | 73 +++++++++++ src/blueberry_circus/potentials.py | 181 +++++++++++++++++++++++++++ src/blueberry_circus/program.py | 78 +++++++++++- tests/test_central_potentials.py | 175 ++++++++++++++++++++++++++ 7 files changed, 663 insertions(+), 10 deletions(-) create mode 100644 docs/central_forces.md create mode 100644 examples/central_forces.py create mode 100644 src/blueberry_circus/conservation.py create mode 100644 tests/test_central_potentials.py diff --git a/docs/central_forces.md b/docs/central_forces.md new file mode 100644 index 0000000..c16a4c6 --- /dev/null +++ b/docs/central_forces.md @@ -0,0 +1,85 @@ +# Central-force potentials + +BlueberryCircus began with two hardcoded binding laws — `Harmonic` and +`Coulomb`. This module generalizes the potential layer into an extensible +**central-force family** while keeping the engine, the radiation-reaction +reduction, and the certificate discipline untouched. + +## The abstraction + +A potential is any object exposing `potential(x)`, `force(x)`, and +`force_jacobian(x)`. The Jacobian is not optional: the Landau–Lifshitz +radiation-reaction reduction needs `dF/dt = J · v`. To make new laws cheap and +correct, `CentralPotential` derives all three Cartesian quantities from three +*radial* functions of the (optionally softened) radius `s = √(r² + softening²)`: + +| supply | meaning | +|---|---| +| `_U(s)` | potential energy `U(r)` | +| `_dU(s)` | `U'(r)` | +| `_d2U(s)` | `U''(r)` | + +With `F(x) = g(s) x`, `g(s) = −U'(s)/s`, the analytic Jacobian is + +``` +J = g(s) I + (g'(s)/s) x xᵀ , g'(s) = (U'(s) − s U''(s)) / s² . +``` + +Every subclass is a handful of lines and inherits an exact `force_jacobian`. + +## The families + +| class | `U(r)` | notes | +|---|---|---| +| `PowerLaw(coeff, p)` | `coeff · rᵖ` | continuous exponent; `p=2` → harmonic, `p=−1` → Coulomb | +| `Yukawa(g, lam)` | `−(g/r) e^{−r/λ}` | screened Coulomb (`λ→∞` → Coulomb) | +| `Morse(De, a, re)` | `De(1−e^{−a(r−re)})² − De` | bounded well, finite dissociation energy | +| `AnharmonicOscillator(omega0, beta)` | `½ m ω₀² r² + ¼ β r⁴` | isotropic Duffing (`β=0` → harmonic) | +| `LennardJones(eps, sigma)` | `4ε[(σ/r)¹² − (σ/r)⁶]` | stiff repulsive core; use `softening` or keep the orbit off the core | + +Singular laws take a Plummer `softening` length that regularizes the `r→0` core, +identical in spirit to `Coulomb.softening`. + +## Declarative use + +The new operations plug into the existing `op | q[i]` program model. Potential +collection in `Program.compile()` is now keyed on the `_PotentialOp` marker, so +future potentials are picked up automatically (no more hardcoded tuple). + +```python +import blueberry_circus as bc + +prog = bc.Program(n_particles=1) +with prog.context as q: + bc.Yukawa(g=0.8, lam=2.0, softening=1e-3) | q[0] + bc.ZPF(band=(0.3, 3.0), n_modes=200, seed=0) | q[0] + bc.RadiationReaction("landau_lifshitz") | q[0] +result = bc.Engine(dt=0.02, t_max=600).run(prog, x0=[0, 0, 0], v0=[0, 0, 0]) +``` + +## Conservation certificates + +`blueberry_circus.conservation` adds two re-checkable certificates for closed +conservative systems (`rr="none"`, no field): + +- `energy_conservation_certificate(traj, potential, particle, tol)` +- `angular_momentum_conservation_certificate(traj, particle, tol)` + +Both encode the fractional peak-to-peak drift of the invariant against the +canonical `residual_le_tol` rule, so an under-resolved integration or a tampered +trajectory re-derives `FAIL`. For a central force `|L|` conservation is the +signature of centrality, and energy conservation cross-checks that a new +potential's analytic `force` is consistent with its `potential` energy — purely +from the integrated dynamics. See `examples/central_forces.py`; all six families +conserve both invariants to ≲10⁻¹¹ at `dt = 5·10⁻³`. + +## Correctness + +`tests/test_central_potentials.py` pins the family three independent ways: + +1. **Finite differences** — analytic `force` and `force_jacobian` match central + differences of `potential` and `force`; the Jacobian is symmetric (curl-free). +2. **Reduction oracles** — `PowerLaw` reproduces `Harmonic` (`p=2`) and `Coulomb` + (`p=−1`) to floating-point closeness. +3. **Conservation** — closed orbits pass the energy and angular-momentum + certificates, and a real drift `FAIL`s an impossibly tight tolerance. diff --git a/examples/central_forces.py b/examples/central_forces.py new file mode 100644 index 0000000..4b282c4 --- /dev/null +++ b/examples/central_forces.py @@ -0,0 +1,63 @@ +"""Central-force potential family: declarative programs + conservation certificates. + +Runs a near-circular orbit under each new potential with radiation reaction off +(a closed conservative system) and prints the re-checkable energy and +angular-momentum conservation certificates. Because the ZPF background is a pure +function of its seed, swapping the potential while holding the seed fixed gives a +*matched-seed counterfactual*: identical driving noise, one edited force law. + + python examples/central_forces.py +""" +import numpy as np + +import blueberry_circus as bc +from blueberry_circus.dynamics import Particle, integrate +from blueberry_circus.conservation import ( + energy_conservation_certificate, angular_momentum_conservation_certificate, +) + +U = bc.Units.scaled(gamma_over_omega0=0.05, omega0=1.0) +P = Particle(U.charge, U.mass) + +POTENTIALS = { + "PowerLaw(p=-1, Kepler)": bc.potentials.PowerLaw(coeff=-0.7, p=-1.0, softening=1e-3), + "PowerLaw(p=2.5)": bc.potentials.PowerLaw(coeff=0.4, p=2.5), + "Yukawa": bc.potentials.Yukawa(g=0.8, lam=3.0, softening=1e-3), + "Morse": bc.potentials.Morse(De=3.0, a=0.8, re=1.0), + "Anharmonic(Duffing)": bc.potentials.AnharmonicOscillator(omega0=1.0, beta=0.3), + "LennardJones": bc.potentials.LennardJones(eps=1.0, sigma=0.6), +} + + +def orbit(pot, r0=1.0, tmax=120.0, dt=0.005): + Fmag = np.linalg.norm(pot.force([r0, 0, 0])) + v = np.sqrt(Fmag * r0 / P.mass) + t = np.arange(0.0, tmax, dt) + return integrate(field=None, potential=pot, particle=P, t_grid=t, + x0=[r0, 0, 0], v0=[0, v, 0], rr="none", units=U, dipole=False) + + +def main(): + print(f"{'potential':24s} {'E drift':>12s} {'|L| drift':>12s} verdicts") + for name, pot in POTENTIALS.items(): + tr = orbit(pot) + e = energy_conservation_certificate(tr, pot, P, tol=1e-6) + l = angular_momentum_conservation_certificate(tr, P, tol=1e-8) + print(f"{name:24s} {e.residual:12.2e} {l.residual:12.2e} " + f"[E:{e.recheck()}] [L:{l.recheck()}]") + + # matched-seed counterfactual teaser: same ZPF seed, two force laws. + print("\nmatched-seed counterfactual (same noise, edited law):") + for law in (bc.Harmonic(omega0=1.0), bc.Anharmonic(omega0=1.0, beta=0.4)): + prog = bc.Program(n_particles=1, units=U) + with prog.context as q: + law | q[0] + bc.ZPF(band=(0.3, 3.0), n_modes=64, seed=7, + mode="one_dimensional", axis=0) | q[0] + bc.RadiationReaction("landau_lifshitz") | q[0] + res = bc.Engine(dt=0.02, t_max=60.0).run(prog, x0=[0.1, 0, 0], v0=[0, 0, 0]) + print(f" {type(law).__name__:12s} ={res.observables['position_variance'][0]:.4e}") + + +if __name__ == "__main__": + main() diff --git a/src/blueberry_circus/__init__.py b/src/blueberry_circus/__init__.py index bfbfef5..bae2cf8 100644 --- a/src/blueberry_circus/__init__.py +++ b/src/blueberry_circus/__init__.py @@ -13,14 +13,17 @@ ALPHA, A0, radiation_reaction_time, setterfield_rescale) from . import (spectrum, oracles, observables, potentials, symplectic, - rectification, tournament) + rectification, tournament, conservation) from .spectrum import rho, spectral_density_Ex, mode_density, mode_energy from .zpf import ZPFBackground from .dynamics import Particle, Trajectory, integrate from .certify import Certificate, RULES, PASS, FAIL, NULL, audit_overclaim, \ save_bundle, load_bundle -from .program import (Program, Harmonic, Coulomb, ZPF, RadiationReaction, +from .program import (Program, Harmonic, Coulomb, PowerLaw, Yukawa, Morse, + Anharmonic, LennardJones, ZPF, RadiationReaction, Operation) +from .conservation import (energy_conservation_certificate, + angular_momentum_conservation_certificate) from .engine import Engine, Result from .tournament import (OrbitState, EnergyLedger, TournamentConfig, HypothesisResult) @@ -31,11 +34,14 @@ "Units", "SI", "BOHR", "EPS0", "HBAR", "C", "E_CHARGE", "M_E", "K_E", "ALPHA", "A0", "radiation_reaction_time", "setterfield_rescale", "spectrum", "oracles", "observables", "potentials", "symplectic", - "rectification", "tournament", "rho", "spectral_density_Ex", "mode_density", - "mode_energy", + "rectification", "tournament", "conservation", "rho", "spectral_density_Ex", + "mode_density", "mode_energy", "ZPFBackground", "Particle", "Trajectory", "integrate", "Certificate", "RULES", "PASS", "FAIL", "NULL", "audit_overclaim", "save_bundle", - "load_bundle", "Program", "Harmonic", "Coulomb", "ZPF", "RadiationReaction", - "Operation", "Engine", "Result", "OrbitState", "EnergyLedger", + "load_bundle", "Program", "Harmonic", "Coulomb", "PowerLaw", "Yukawa", + "Morse", "Anharmonic", "LennardJones", "ZPF", "RadiationReaction", + "Operation", "energy_conservation_certificate", + "angular_momentum_conservation_certificate", + "Engine", "Result", "OrbitState", "EnergyLedger", "TournamentConfig", "HypothesisResult", "__version__", ] diff --git a/src/blueberry_circus/conservation.py b/src/blueberry_circus/conservation.py new file mode 100644 index 0000000..0570f30 --- /dev/null +++ b/src/blueberry_circus/conservation.py @@ -0,0 +1,73 @@ +"""Conservation-law certificates for closed, conservative systems. + +A closed central-force system integrated with ``rr="none"`` and no external +field must conserve mechanical energy ``E = 1/2 m v^2 + U(r)`` and, because the +force is central, the angular-momentum magnitude ``|L| = |m x cross v|``. These +are re-checkable :class:`~blueberry_circus.certify.Certificate` claims in exactly +the repo's assurance style: the residual is the fractional drift of the quantity +over the window, encoded against the canonical ``residual_le_tol`` rule, so a +tampered trajectory (or an under-resolved integration) re-derives ``FAIL``. + +They are the correctness backbone for the central-force potential family: any new +conservative potential can self-certify that its analytic force is consistent +with its potential energy, purely from the integrated dynamics. +""" +from __future__ import annotations + +import numpy as np + +from .certify import Certificate, finalize +from .observables import total_energy, angular_momentum + +_METHOD = "blueberry_circus closed-system RK4 conservation check" + + +def _fractional_drift(series) -> float: + """(max - min) / |mean| -- the fractional peak-to-peak drift of a conserved + quantity. Returns a finite over-threshold sentinel if the series is not all + finite, so the certificate stays canonicalizable and re-derives FAIL.""" + series = np.asarray(series, float) + if not np.all(np.isfinite(series)): + return np.inf + mean = series.mean() + denom = abs(mean) if mean != 0.0 else 1.0 + return float((series.max() - series.min()) / denom) + + +def _residual_certificate(kind, claim, value, residual, tolerance, provenance): + # Keep the hash surface clean: a non-finite drift is recorded as a finite + # FAIL sentinel and the value is dropped to None (mirrors rel_error_certificate). + if not np.isfinite(residual): + residual = abs(float(tolerance)) * 2.0 + 1.0 + value = None + cert = Certificate( + kind=kind, claim=claim, method=_METHOD, rule="residual_le_tol", + value=(float(value) if value is not None else None), + residual=float(residual), tolerance=float(tolerance), + provenance=provenance) + return finalize(cert) + + +def energy_conservation_certificate(traj, potential, particle, tol=1e-6): + """Certify that mechanical energy drifts by less than ``tol`` (fractional).""" + E = total_energy(traj, potential, particle) + return _residual_certificate( + kind="energy_conservation", + claim="mechanical energy is conserved along the trajectory", + value=(np.mean(E) if np.all(np.isfinite(E)) else None), + residual=_fractional_drift(E), tolerance=tol, + provenance={"n_steps": int(len(traj.t)), + "quantity": "E_peak_to_peak_over_mean"}) + + +def angular_momentum_conservation_certificate(traj, particle, tol=1e-6): + """Certify that ``|L|`` drifts by less than ``tol`` (fractional) -- the + signature of a central force.""" + Lmag = np.linalg.norm(angular_momentum(traj, particle), axis=1) + return _residual_certificate( + kind="angular_momentum_conservation", + claim="angular-momentum magnitude is conserved (central force)", + value=(np.mean(Lmag) if np.all(np.isfinite(Lmag)) else None), + residual=_fractional_drift(Lmag), tolerance=tol, + provenance={"n_steps": int(len(traj.t)), + "quantity": "absL_peak_to_peak_over_mean"}) diff --git a/src/blueberry_circus/potentials.py b/src/blueberry_circus/potentials.py index b0156ff..54410f1 100644 --- a/src/blueberry_circus/potentials.py +++ b/src/blueberry_circus/potentials.py @@ -67,3 +67,184 @@ def force_jacobian(self, x): rs = self._rs(x) I = np.eye(len(x)) return -self._coef * (I / rs**3 - 3.0 * np.outer(x, x) / rs**5) + + +# --------------------------------------------------------------------------- +# Central-force family +# --------------------------------------------------------------------------- +class CentralPotential: + """Base for isotropic central potentials ``U = U(r)``. + + A subclass supplies only the three *radial* functions ``_U(s)``, + ``_dU(s) = dU/ds`` and ``_d2U(s) = d^2U/ds^2`` of the (optionally softened) + radius ``s = sqrt(r^2 + softening^2)``. The Cartesian force and its Jacobian + then follow from the chain rule, so every central potential inherits an + analytic ``force_jacobian`` -- the quantity the Landau--Lifshitz + radiation-reaction reduction needs (``dF/dt = J . v``). + + With ``F(x) = g(s) x`` and ``g(s) = -U'(s)/s`` the Jacobian is:: + + J = g(s) I + (g'(s)/s) x x^T , g'(s) = (U'(s) - s U''(s)) / s^2 . + + The ``softening`` length regularizes the ``r -> 0`` core for singular laws + (Coulomb-like power laws, Yukawa, Lennard-Jones) exactly as in + :class:`Coulomb`; leave it at 0 for regular wells (harmonic, Morse, Duffing). + """ + + softening: float = 0.0 + + def _s(self, x): + x = np.asarray(x, float) + return np.sqrt(np.dot(x, x) + self.softening**2) + + # radial functions -- implemented by subclasses + def _U(self, s): # pragma: no cover - abstract + raise NotImplementedError + def _dU(self, s): # pragma: no cover - abstract + raise NotImplementedError + def _d2U(self, s): # pragma: no cover - abstract + raise NotImplementedError + + # Cartesian interface (shared by the integrator) + def potential(self, x): + return float(self._U(self._s(x))) + + def force(self, x): + x = np.asarray(x, float) + s = self._s(x) + return -(self._dU(s) / s) * x + + def force_jacobian(self, x): + x = np.asarray(x, float) + s = self._s(x) + g = -self._dU(s) / s + gprime = (self._dU(s) - s * self._d2U(s)) / s**2 + return g * np.eye(len(x)) + (gprime / s) * np.outer(x, x) + + +@dataclass +class PowerLaw(CentralPotential): + """Central power law ``U = coeff * r**p``. + + Generalizes the two hardcoded laws: ``p=2, coeff=1/2 m omega0^2`` reproduces + :class:`Harmonic`; ``p=-1, coeff=-Z k_e q^2`` reproduces :class:`Coulomb`. + Continuous ``p`` gives the ``F ~ r^(p-1)`` force families used as the + inverse-dynamics generalization axis. + """ + coeff: float + p: float + mass: float = 1.0 + softening: float = 0.0 + + def _U(self, s): + return self.coeff * s**self.p + + def _dU(self, s): + return self.coeff * self.p * s**(self.p - 1.0) + + def _d2U(self, s): + return self.coeff * self.p * (self.p - 1.0) * s**(self.p - 2.0) + + +@dataclass +class Yukawa(CentralPotential): + """Screened Coulomb (Yukawa) ``U = -(g / r) exp(-r / lam)``. + + Attractive for ``g > 0``; ``lam`` is the screening length (``lam -> inf`` + recovers Coulomb). Set ``softening`` below the orbit scale to regularize the + ``r -> 0`` core. + """ + g: float = 1.0 + lam: float = 1.0 + mass: float = 1.0 + softening: float = 0.0 + + def _U(self, s): + return -self.g * np.exp(-s / self.lam) / s + + def _dU(self, s): + e = np.exp(-s / self.lam) + return self.g * e * (1.0 / (self.lam * s) + 1.0 / s**2) + + def _d2U(self, s): + e = np.exp(-s / self.lam) + L = self.lam + return self.g * e * (-1.0 / (L**2 * s) - 2.0 / (L * s**2) - 2.0 / s**3) + + +@dataclass +class Morse(CentralPotential): + """Radial Morse well ``U = De (1 - e^{-a(r - re)})^2 - De``. + + Minimum ``U(re) = -De``; ``U(inf) = 0``. A bounded, anharmonic binding well + with a soft outer wall -- a qualitatively different force family from the + power laws (finite dissociation energy). + """ + De: float = 1.0 + a: float = 1.0 + re: float = 1.0 + mass: float = 1.0 + softening: float = 0.0 + + def _w(self, s): + return np.exp(-self.a * (s - self.re)) + + def _U(self, s): + w = self._w(s) + return self.De * (1.0 - w)**2 - self.De + + def _dU(self, s): + w = self._w(s) + return 2.0 * self.De * self.a * w * (1.0 - w) + + def _d2U(self, s): + w = self._w(s) + return -2.0 * self.De * self.a**2 * w * (1.0 - 2.0 * w) + + +@dataclass +class AnharmonicOscillator(CentralPotential): + """Isotropic Duffing well ``U = 1/2 m omega0^2 r^2 + 1/4 beta r^4``. + + ``beta > 0`` hardening, ``beta < 0`` softening; ``beta = 0`` recovers + :class:`Harmonic`. The canonical nonlinear-oscillator testbed. + """ + omega0: float = 1.0 + beta: float = 0.0 + mass: float = 1.0 + softening: float = 0.0 + + def _U(self, s): + return 0.5 * self.mass * self.omega0**2 * s**2 + 0.25 * self.beta * s**4 + + def _dU(self, s): + return self.mass * self.omega0**2 * s + self.beta * s**3 + + def _d2U(self, s): + return self.mass * self.omega0**2 + 3.0 * self.beta * s**2 + + +@dataclass +class LennardJones(CentralPotential): + """12-6 Lennard-Jones ``U = 4 eps [ (sigma/r)^12 - (sigma/r)^6 ]``. + + Stiff repulsive core plus a shallow attractive well (minimum at + ``r = 2^(1/6) sigma``). The steep ``r^-12`` wall makes this the stress test + for the integrator; keep the orbit away from the core or add ``softening``. + """ + eps: float = 1.0 + sigma: float = 1.0 + mass: float = 1.0 + softening: float = 0.0 + + def _U(self, s): + u = self.sigma / s + return 4.0 * self.eps * (u**12 - u**6) + + def _dU(self, s): + u = self.sigma / s + return -(4.0 * self.eps / s) * (12.0 * u**12 - 6.0 * u**6) + + def _d2U(self, s): + u = self.sigma / s + return (4.0 * self.eps / s**2) * (156.0 * u**12 - 42.0 * u**6) diff --git a/src/blueberry_circus/program.py b/src/blueberry_circus/program.py index 7393ec8..2078c6b 100644 --- a/src/blueberry_circus/program.py +++ b/src/blueberry_circus/program.py @@ -32,8 +32,17 @@ def __or__(self, regref: "RegRef"): return self +class _PotentialOp(Operation): + """Marker for ops that compile to a single binding potential. + + ``compile()`` collects potential ops by this base rather than a hardcoded + ``(Harmonic, Coulomb)`` tuple, so new potential families are picked up + automatically. Every subclass must provide ``build_potential(units, mass)``. + """ + + @dataclass -class Harmonic(Operation): +class Harmonic(_PotentialOp): omega0: float def build_potential(self, units: Units, mass: float): @@ -41,7 +50,7 @@ def build_potential(self, units: Units, mass: float): @dataclass -class Coulomb(Operation): +class Coulomb(_PotentialOp): Z: float = 1.0 softening: float = 0.0 @@ -50,6 +59,67 @@ def build_potential(self, units: Units, mass: float): mass=mass) +@dataclass +class PowerLaw(_PotentialOp): + """Central power law ``U = coeff * r**p`` (see :class:`potentials.PowerLaw`).""" + coeff: float + p: float + softening: float = 0.0 + + def build_potential(self, units: Units, mass: float): + return _pot.PowerLaw(coeff=self.coeff, p=self.p, mass=mass, + softening=self.softening) + + +@dataclass +class Yukawa(_PotentialOp): + """Screened Coulomb ``U = -(g/r) exp(-r/lam)``.""" + g: float = 1.0 + lam: float = 1.0 + softening: float = 0.0 + + def build_potential(self, units: Units, mass: float): + return _pot.Yukawa(g=self.g, lam=self.lam, mass=mass, + softening=self.softening) + + +@dataclass +class Morse(_PotentialOp): + """Radial Morse well ``U = De (1 - e^{-a(r-re)})^2 - De``.""" + De: float = 1.0 + a: float = 1.0 + re: float = 1.0 + softening: float = 0.0 + + def build_potential(self, units: Units, mass: float): + return _pot.Morse(De=self.De, a=self.a, re=self.re, mass=mass, + softening=self.softening) + + +@dataclass +class Anharmonic(_PotentialOp): + """Isotropic Duffing well ``U = 1/2 m w0^2 r^2 + 1/4 beta r^4``.""" + omega0: float = 1.0 + beta: float = 0.0 + softening: float = 0.0 + + def build_potential(self, units: Units, mass: float): + return _pot.AnharmonicOscillator(omega0=self.omega0, beta=self.beta, + mass=mass, softening=self.softening) + + +@dataclass +class LennardJones(_PotentialOp): + """12-6 Lennard-Jones ``U = 4 eps [(sigma/r)^12 - (sigma/r)^6]``.""" + eps: float = 1.0 + sigma: float = 1.0 + softening: float = 0.0 + + def build_potential(self, units: Units, mass: float): + return _pot.LennardJones(eps=self.eps, sigma=self.sigma, mass=mass, + softening=self.softening) + + @dataclass class ZPF(Operation): band: tuple @@ -139,7 +209,7 @@ def compile(self, index: int = 0) -> "CompiledProgram": passes.append("validate_single_particle") ops = self.ops_for(index) - pots = [op for op in ops if isinstance(op, (Harmonic, Coulomb))] + pots = [op for op in ops if isinstance(op, _PotentialOp)] if len(pots) != 1: raise ValueError(f"collect_potential: exactly one potential required " f"for particle {index}, got {len(pots)}") @@ -160,7 +230,7 @@ def compile(self, index: int = 0) -> "CompiledProgram": rr = reactions[0].build() if reactions else "none" passes.append("resolve_reaction") - handled = (Harmonic, Coulomb, ZPF, RadiationReaction) + handled = (_PotentialOp, ZPF, RadiationReaction) unknown = [op for op in ops if not isinstance(op, handled)] if unknown: raise TypeError(f"compile: unhandled operation(s) {unknown!r}") diff --git a/tests/test_central_potentials.py b/tests/test_central_potentials.py new file mode 100644 index 0000000..c0c5e60 --- /dev/null +++ b/tests/test_central_potentials.py @@ -0,0 +1,175 @@ +"""Tests for the central-force potential family. + +Correctness is pinned three independent ways: + +* analytic ``force`` and ``force_jacobian`` agree with central finite differences + of ``potential`` and ``force`` (so the RR Jacobian is trustworthy); +* the general :class:`PowerLaw` reproduces the hardcoded :class:`Harmonic` and + :class:`Coulomb` bit-closely (a reduction oracle); +* closed conservative orbits self-certify energy and angular-momentum + conservation via :mod:`blueberry_circus.conservation`. +""" +import numpy as np +import pytest + +import blueberry_circus as bc +from blueberry_circus.constants import Units +from blueberry_circus.dynamics import Particle, integrate +from blueberry_circus.potentials import ( + Harmonic, Coulomb, PowerLaw, Yukawa, Morse, AnharmonicOscillator, + LennardJones, +) +from blueberry_circus.conservation import ( + energy_conservation_certificate, angular_momentum_conservation_certificate, +) +from blueberry_circus.certify import PASS + +U = Units.scaled(gamma_over_omega0=0.05, omega0=1.0) +P = Particle(U.charge, U.mass) + +# Evaluation points kept away from the r->0 core (softening handles the rest). +PTS = [np.array([0.9, 0.0, 0.0]), + np.array([0.7, -0.4, 0.5]), + np.array([-0.6, 0.8, -0.3]), + np.array([1.3, 0.2, -0.9])] + +# One representative instance of each family (softening where the core is steep). +INSTANCES = [ + PowerLaw(coeff=0.5 * 1.3 * 0.9**2, p=2.0, mass=1.3), # harmonic-like + PowerLaw(coeff=-0.7, p=-1.0, softening=1e-3), # Coulomb-like + PowerLaw(coeff=0.4, p=2.5, mass=1.0), # fractional exponent + Yukawa(g=0.8, lam=1.5, softening=1e-3), + Morse(De=1.2, a=1.1, re=1.0), + AnharmonicOscillator(omega0=0.9, beta=0.5, mass=1.1), + LennardJones(eps=1.0, sigma=0.6), # well near r=0.67 +] + + +def _num_grad_U(pot, x, h=1e-6): + g = np.zeros(3) + for i in range(3): + d = np.zeros(3); d[i] = h + g[i] = (pot.potential(x + d) - pot.potential(x - d)) / (2 * h) + return g + + +def _num_jac_F(pot, x, h=1e-6): + J = np.zeros((3, 3)) + for j in range(3): + d = np.zeros(3); d[j] = h + J[:, j] = (pot.force(x + d) - pot.force(x - d)) / (2 * h) + return J + + +@pytest.mark.parametrize("pot", INSTANCES) +def test_force_is_minus_grad_potential(pot): + for x in PTS: + assert np.allclose(pot.force(x), -_num_grad_U(pot, x), + rtol=1e-4, atol=1e-6) + + +@pytest.mark.parametrize("pot", INSTANCES) +def test_force_jacobian_matches_finite_difference(pot): + for x in PTS: + assert np.allclose(pot.force_jacobian(x), _num_jac_F(pot, x), + rtol=1e-4, atol=1e-5) + + +@pytest.mark.parametrize("pot", INSTANCES) +def test_jacobian_is_symmetric(pot): + # A conservative (curl-free) force has a symmetric Jacobian (Hessian of -U). + for x in PTS: + J = pot.force_jacobian(x) + assert np.allclose(J, J.T, atol=1e-10) + + +def test_powerlaw_reduces_to_harmonic(): + m, w = 1.3, 0.9 + pl = PowerLaw(coeff=0.5 * m * w**2, p=2.0, mass=m) + ha = Harmonic(w, mass=m) + for x in PTS: + assert np.isclose(pl.potential(x), ha.potential(x)) + assert np.allclose(pl.force(x), ha.force(x)) + assert np.allclose(pl.force_jacobian(x), ha.force_jacobian(x)) + + +def test_powerlaw_reduces_to_coulomb(): + coul = Coulomb(Z=1.0, units=U, charge=U.charge, mass=U.mass) + pl = PowerLaw(coeff=-coul._coef, p=-1.0) + for x in PTS: + assert np.isclose(pl.potential(x), coul.potential(x)) + assert np.allclose(pl.force(x), coul.force(x)) + assert np.allclose(pl.force_jacobian(x), coul.force_jacobian(x)) + + +def _circular_orbit(pot, r0=1.0, tmax=120.0, dt=0.005): + """Integrate a near-circular orbit in the plane for a central potential.""" + Fmag = np.linalg.norm(pot.force([r0, 0, 0])) + v = np.sqrt(Fmag * r0 / P.mass) # circular-orbit speed + t = np.arange(0.0, tmax, dt) + return integrate(field=None, potential=pot, particle=P, t_grid=t, + x0=[r0, 0, 0], v0=[0, v, 0], rr="none", units=U, + dipole=False) + + +@pytest.mark.parametrize("pot", [ + PowerLaw(coeff=-0.7, p=-1.0, softening=1e-3), + Yukawa(g=0.8, lam=3.0, softening=1e-3), + AnharmonicOscillator(omega0=1.0, beta=0.3), + Morse(De=3.0, a=0.8, re=1.0), +]) +def test_closed_orbit_conserves_energy_and_angular_momentum(pot): + tr = _circular_orbit(pot) + assert np.all(np.isfinite(tr.x)) + e_cert = energy_conservation_certificate(tr, pot, P, tol=1e-6) + l_cert = angular_momentum_conservation_certificate(tr, P, tol=1e-8) + assert e_cert.recheck() == PASS, f"energy drift residual={e_cert.residual:.2e}" + assert l_cert.recheck() == PASS, f"|L| drift residual={l_cert.residual:.2e}" + + +def test_conservation_certificate_fails_on_impossible_tolerance(): + # A real (small but nonzero) drift must FAIL an absurdly tight tolerance: + # guards against a certificate that trivially passes everything. + tr = _circular_orbit(AnharmonicOscillator(omega0=1.0, beta=0.3)) + cert = energy_conservation_certificate(tr, AnharmonicOscillator( + omega0=1.0, beta=0.3), P, tol=1e-18) + assert cert.recheck() != PASS + + +@pytest.mark.parametrize("pot", [ + Morse(De=2.0, a=1.0, re=1.0), + AnharmonicOscillator(omega0=1.0, beta=0.5), + Yukawa(g=1.0, lam=2.0, softening=1e-2), +]) +def test_radiation_reaction_stays_finite(pot): + # The analytic Jacobian feeds the Landau-Lifshitz term; integration must not + # blow up for the new families. + Ud = Units.scaled(gamma_over_omega0=0.02, omega0=1.0) + Pd = Particle(Ud.charge, Ud.mass) + t = np.arange(0.0, 60.0, 0.01) + tr = integrate(field=None, potential=pot, particle=Pd, t_grid=t, + x0=[1.0, 0, 0], v0=[0, 0.5, 0], rr="landau_lifshitz", + units=Ud, dipole=False) + assert np.all(np.isfinite(tr.x)) and np.all(np.isfinite(tr.v)) + + +def test_declarative_program_runs_new_potential(): + prog = bc.Program(n_particles=1, units=U) + with prog.context as q: + bc.Yukawa(g=0.8, lam=2.0, softening=1e-3) | q[0] + bc.ZPF(band=(0.3, 3.0), n_modes=64, seed=1, + mode="one_dimensional", axis=0) | q[0] + bc.RadiationReaction("landau_lifshitz") | q[0] + res = bc.Engine(backend="numpy", dt=0.02, t_max=40.0).run( + prog, x0=[1.0, 0, 0], v0=[0, 0.3, 0]) + assert res.observables["trajectory_finite"] is True + assert res.certificates[0].recheck() == PASS + + +def test_program_still_requires_exactly_one_potential(): + prog = bc.Program(n_particles=1, units=U) + with prog.context as q: + bc.Harmonic(omega0=1.0) | q[0] + bc.Morse(De=1.0, a=1.0, re=1.0) | q[0] + with pytest.raises(ValueError): + prog.compile(0)