diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 7a9cba3..cf383ff 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -29,7 +29,7 @@ repos: files: *python_files - repo: https://github.com/pre-commit/mirrors-mypy - rev: v1.9.0 + rev: v1.19.1 hooks: - id: mypy args: [--ignore-missing-imports, --no-error-summary] diff --git a/CHANGELOG.md b/CHANGELOG.md index ed9c8c5..c068f72 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,19 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/). +## [1.4.3] - 2026-07-06 + +### Added + +- R1 light-current helpers for package-native gravitational-lensing probes: + `planar_r1_light_packet()`, `r1_vacuum_subtracted_potential()`, + `r1_light_acceleration()`, and `r1_light_step()`. +- Noether-current readouts: `noether_spatial_current()` and + `positive_noether_current()`. +- 19-point-stencil-consistent FFT Poisson/equilibrium helpers: + `poisson_solve_fft_19pt()`, `equilibrate_chi_19pt()`, and + `equilibrate_from_fields_19pt()`. + ## [1.4.1] - 2026-04-24 ### Fixed diff --git a/lfm/__init__.py b/lfm/__init__.py index 0451e02..582e1cd 100644 --- a/lfm/__init__.py +++ b/lfm/__init__.py @@ -15,7 +15,7 @@ print(sim.metrics()) """ -__version__ = "1.4.1" +__version__ = "1.4.5" from lfm.analysis import ( angular_momentum_density, @@ -60,9 +60,14 @@ measure_force, measure_separation, metric_perturbation, + metric_refractive_index, momentum_density, + noether_spatial_current, + op05_spherical_chi_deflection, phase_coherence, + phase_current_energy_density, phase_field, + positive_noether_current, power_spectrum, precession_rate, project_field_onto_modes, @@ -71,6 +76,7 @@ rotation_curve, rotation_curve_fit, schwarzschild_chi, + schwarzschild_radius_si, smoothed_color_variance, sparc_load, spinor_center_of_energy, @@ -99,7 +105,9 @@ AGE_UNIVERSE_GYR, ALPHA_EM, ALPHA_S, + ARCSEC_PER_RADIAN, BETA_0, + C_SI, CHI0, D_ST, DT_DEFAULT, @@ -107,6 +115,7 @@ E_AMPLITUDE_BY_GRID, EPSILON_CC, EPSILON_W, + G_SI, KAPPA, KAPPA_C, KAPPA_STRING, @@ -125,6 +134,8 @@ SA_GAMMA, SA_L, SIN2_THETA_W, + SOLAR_MASS_KG, + SOLAR_RADIUS_M, TOTAL_RADIUS_LOWER_BOUND_PLANCK, Z2_COORD, D, @@ -156,13 +167,20 @@ disk_positions, disk_velocities, equilibrate_chi, + equilibrate_chi_19pt, equilibrate_from_fields, + equilibrate_from_fields_19pt, gaussian_soliton, gaussian_spinor, grid_positions, initialize_disk, place_solitons, + planar_r1_light_packet, poisson_solve_fft, + poisson_solve_fft_19pt, + r1_light_acceleration, + r1_light_step, + r1_vacuum_subtracted_potential, seed_noise, sparse_positions, spherical_phase_source, @@ -305,6 +323,11 @@ "CHI0", "D", "D_ST", + "G_SI", + "C_SI", + "SOLAR_MASS_KG", + "SOLAR_RADIUS_M", + "ARCSEC_PER_RADIAN", "KAPPA", "KAPPA_C", "KAPPA_STRING", @@ -362,13 +385,20 @@ "place_solitons", "wave_kick", "poisson_solve_fft", + "poisson_solve_fft_19pt", "equilibrate_chi", + "equilibrate_chi_19pt", "equilibrate_from_fields", + "equilibrate_from_fields_19pt", "seed_noise", "uniform_chi", "tetrahedral_positions", "sparse_positions", "spherical_phase_source", + "planar_r1_light_packet", + "r1_vacuum_subtracted_potential", + "r1_light_acceleration", + "r1_light_step", "grid_positions", "disk_positions", "disk_velocities", @@ -454,11 +484,17 @@ "time_dilation_factor", "gravitational_potential", "schwarzschild_chi", + "schwarzschild_radius_si", + "metric_refractive_index", + "op05_spherical_chi_deflection", "find_apparent_horizon", "horizon_mass", # Phase (EM / charge) "phase_field", "charge_density", + "noether_spatial_current", + "positive_noether_current", + "phase_current_energy_density", "phase_coherence", "coulomb_interaction_energy", # Angular momentum diff --git a/lfm/analysis/__init__.py b/lfm/analysis/__init__.py index dd9e93c..262138b 100644 --- a/lfm/analysis/__init__.py +++ b/lfm/analysis/__init__.py @@ -63,7 +63,10 @@ gravitational_potential, horizon_mass, metric_perturbation, + metric_refractive_index, + op05_spherical_chi_deflection, schwarzschild_chi, + schwarzschild_radius_si, time_dilation_factor, ) from lfm.analysis.metrics import compute_metrics @@ -83,8 +86,11 @@ from lfm.analysis.phase import ( charge_density, coulomb_interaction_energy, + noether_spatial_current, phase_coherence, + phase_current_energy_density, phase_field, + positive_noether_current, ) from lfm.analysis.ringdown import ( fit_ringdown_series, @@ -170,11 +176,17 @@ "time_dilation_factor", "gravitational_potential", "schwarzschild_chi", + "schwarzschild_radius_si", + "metric_refractive_index", + "op05_spherical_chi_deflection", "find_apparent_horizon", "horizon_mass", # phase (EM / charge) "phase_field", "charge_density", + "noether_spatial_current", + "positive_noether_current", + "phase_current_energy_density", "phase_coherence", "coulomb_interaction_energy", # ringdown extraction diff --git a/lfm/analysis/metric.py b/lfm/analysis/metric.py index 411f0c3..ef7081a 100644 --- a/lfm/analysis/metric.py +++ b/lfm/analysis/metric.py @@ -13,7 +13,7 @@ import numpy as np -from lfm.constants import CHI0 +from lfm.constants import ARCSEC_PER_RADIAN, C_SI, CHI0, G_SI if TYPE_CHECKING: from numpy.typing import NDArray @@ -152,6 +152,109 @@ def schwarzschild_chi( return chi.astype(np.float32) +def schwarzschild_radius_si( + mass_kg: float, + gravitational_constant: float = G_SI, + c_si: float = C_SI, +) -> float: + """Return the Schwarzschild radius in meters for an SI mass input.""" + if mass_kg <= 0.0: + raise ValueError("mass_kg must be positive") + if gravitational_constant <= 0.0: + raise ValueError("gravitational_constant must be positive") + if c_si <= 0.0: + raise ValueError("c_si must be positive") + return float(2.0 * gravitational_constant * mass_kg / (c_si * c_si)) + + +def metric_refractive_index( + chi: NDArray, + chi0: float = CHI0, + ppn_gamma: float = 1.0, +) -> NDArray: + """Return the LFM geometric-optics refractive index from the chi metric. + + The GOV-01 metric map gives g00 = -(chi/chi0)^2. In the weak-field + optical limit, the PPN spatial metric contribution gives + n = (chi0 / chi) ** (1 + gamma). The canonical LFM weak-GR closure has + gamma = 1, so n = (chi0 / chi) ** 2. + """ + if chi0 <= 0.0: + raise ValueError("chi0 must be positive") + if ppn_gamma < 0.0: + raise ValueError("ppn_gamma must be non-negative") + chi_f = np.asarray(chi, dtype=np.float64) + if np.any(chi_f <= 0.0): + raise ValueError("chi must be positive for metric refractive index") + return np.power(chi0 / chi_f, 1.0 + ppn_gamma) + + +def op05_spherical_chi_deflection( + mass_kg: float, + impact_parameter_m: float, + x_extent_multiplier: float = 500.0, + sample_count: int = 20001, + ppn_gamma: float = 1.0, +) -> dict[str, object]: + """Integrate OP-05 for a spherical GR-16 chi profile. + + This uses the canonical LFM chain: + + - GR-16: chi/chi0 = sqrt(1 - r_s / r) + - GR-24: gamma = 1 for the weak-field spatial metric + - OP-05: dtheta/dx = (1/n) * partial_y n + + The returned comparator is not used by the integration; it is the + closed-form weak-field value 2*r_s/b for checking the numerical result. + """ + if impact_parameter_m <= 0.0: + raise ValueError("impact_parameter_m must be positive") + if x_extent_multiplier <= 1.0: + raise ValueError("x_extent_multiplier must be greater than 1") + if sample_count < 101: + raise ValueError("sample_count must be at least 101") + if sample_count % 2 == 0: + sample_count += 1 + + rs_m = schwarzschild_radius_si(mass_kg) + x_extent_m = float(x_extent_multiplier) * impact_parameter_m + x_m = np.linspace(-x_extent_m, x_extent_m, sample_count, dtype=np.float64) + y0_m = np.full_like(x_m, impact_parameter_m) + radius_m = np.sqrt(x_m * x_m + y0_m * y0_m) + exponent = 0.5 * (1.0 + ppn_gamma) + safe = np.maximum(1.0 - rs_m / radius_m, 1.0e-15) + n_eff = np.power(safe, -exponent) + + dn_dr = -exponent * rs_m / (radius_m * radius_m) * np.power(safe, -exponent - 1.0) + dn_dy = dn_dr * y0_m / radius_m + dtheta_dx = dn_dy / n_eff + angle_rad = abs(float(np.trapezoid(dtheta_dx, x_m))) + comparator_rad = 2.0 * rs_m / impact_parameter_m + + increments = 0.5 * (dtheta_dx[1:] + dtheta_dx[:-1]) * np.diff(x_m) + theta_rad = np.concatenate([[0.0], np.cumsum(increments)]) + y_increments = 0.5 * (theta_rad[1:] + theta_rad[:-1]) * np.diff(x_m) + y_m = impact_parameter_m + np.concatenate([[0.0], np.cumsum(y_increments)]) + + return { + "mass_kg": float(mass_kg), + "impact_parameter_m": float(impact_parameter_m), + "schwarzschild_radius_m": float(rs_m), + "x_extent_multiplier": float(x_extent_multiplier), + "sample_count": int(sample_count), + "ppn_gamma": float(ppn_gamma), + "x_over_b": (x_m / impact_parameter_m).tolist(), + "y_over_b": (y_m / impact_parameter_m).tolist(), + "theta_arcsec": (np.abs(theta_rad) * ARCSEC_PER_RADIAN).tolist(), + "n_eff": n_eff.tolist(), + "dtheta_dx": dtheta_dx.tolist(), + "recovered_angle_radians": float(angle_rad), + "recovered_angle_arcsec": float(angle_rad * ARCSEC_PER_RADIAN), + "canonical_comparator_arcsec": float(comparator_rad * ARCSEC_PER_RADIAN), + "comparator_relative_error": float((angle_rad - comparator_rad) / comparator_rad), + } + + # --------------------------------------------------------------------------- # Apparent horizon detection (v16 black-hole analysis) # --------------------------------------------------------------------------- diff --git a/lfm/analysis/phase.py b/lfm/analysis/phase.py index 3345337..5265d52 100644 --- a/lfm/analysis/phase.py +++ b/lfm/analysis/phase.py @@ -19,6 +19,14 @@ from numpy.typing import NDArray +def _runtime_float_dtype(*arrays: object) -> type[np.float32] | type[np.float64]: + for arr in arrays: + dtype = getattr(arr, "dtype", None) + if dtype is not None and np.dtype(dtype) == np.dtype(np.float64): + return np.float64 + return np.float32 + + def phase_field( psi_r: NDArray, psi_i: NDArray, @@ -76,6 +84,84 @@ def charge_density( return psi_r * dpsi_i_dt - psi_i * dpsi_r_dt +def noether_spatial_current( + psi_r: NDArray, + psi_i: NDArray, + axis: int = 0, +) -> NDArray: + """Compute spatial Noether current along one lattice axis. + + j_axis = Im(conj(Psi) * d_axis Psi). A centered finite difference is + used on the periodic lattice. + """ + if axis not in (0, 1, 2): + raise ValueError("axis must be 0, 1, or 2") + dpsi_r = 0.5 * (np.roll(psi_r, -1, axis=axis) - np.roll(psi_r, 1, axis=axis)) + dpsi_i = 0.5 * (np.roll(psi_i, -1, axis=axis) - np.roll(psi_i, 1, axis=axis)) + out_dtype = _runtime_float_dtype(psi_r, psi_i) + return (psi_r * dpsi_i - psi_i * dpsi_r).astype(out_dtype) + + +def positive_noether_current( + psi_r: NDArray, + psi_i: NDArray, + axis: int = 0, +) -> NDArray: + """Return only the positive outgoing part of spatial Noether current.""" + out_dtype = _runtime_float_dtype(psi_r, psi_i) + return np.maximum(noether_spatial_current(psi_r, psi_i, axis=axis), out_dtype(0.0)).astype( + out_dtype + ) + + +def phase_current_energy_density( + psi_r: NDArray, + psi_i: NDArray, + psi_r_prev: NDArray, + psi_i_prev: NDArray, + dt: float, + c_speed: float = 1.0, + amplitude_floor: float = 1.0e-30, +) -> NDArray: + """Return the phase-current stress-energy component of a complex wave. + + A pure phase photon can carry energy while ``|Psi|^2`` remains nearly + constant. This observable extracts that missing source from the U(1) + Noether current: + + rho_phase = 0.5 * (j_0^2 + c^2 |j_space|^2) / |Psi|^2 + + where ``j_0 = Im(conj(Psi) d_t Psi)`` and + ``j_i = Im(conj(Psi) d_i Psi)``. For ``Psi = A exp(i theta)``, this is + ``0.5 * A^2 * (theta_t^2 + c^2 |grad theta|^2)``. It is invariant under + global phase rotations and vanishes for a static uniform phase. + """ + if dt <= 0.0: + raise ValueError("dt must be positive") + if c_speed < 0.0: + raise ValueError("c_speed must be non-negative") + out_dtype = _runtime_float_dtype(psi_r, psi_i, psi_r_prev, psi_i_prev) + + psi_r_f = psi_r.astype(out_dtype, copy=False) + psi_i_f = psi_i.astype(out_dtype, copy=False) + psi_r_prev_f = psi_r_prev.astype(out_dtype, copy=False) + psi_i_prev_f = psi_i_prev.astype(out_dtype, copy=False) + + dpsi_r_dt = (psi_r_f - psi_r_prev_f) / out_dtype(dt) + dpsi_i_dt = (psi_i_f - psi_i_prev_f) / out_dtype(dt) + j0 = psi_r_f * dpsi_i_dt - psi_i_f * dpsi_r_dt + + jx = noether_spatial_current(psi_r_f, psi_i_f, axis=0) + jy = noether_spatial_current(psi_r_f, psi_i_f, axis=1) + jz = noether_spatial_current(psi_r_f, psi_i_f, axis=2) + + amp_sq = psi_r_f * psi_r_f + psi_i_f * psi_i_f + amp_safe = np.maximum(amp_sq, out_dtype(amplitude_floor)) + c2 = out_dtype(c_speed * c_speed) + energy = 0.5 * (j0 * j0 + c2 * (jx * jx + jy * jy + jz * jz)) / amp_safe + return energy.astype(out_dtype) + + def phase_coherence( psi_r: NDArray, psi_i: NDArray, diff --git a/lfm/constants.py b/lfm/constants.py index e5e5b27..797d2bb 100644 --- a/lfm/constants.py +++ b/lfm/constants.py @@ -33,6 +33,21 @@ CHI0: float = float(3**D - 2**D) """Background χ value = 19.0. From 3D discrete Laplacian: 1 center + 6 face + 12 edge modes.""" +G_SI: float = 6.67430e-11 +"""CODATA Newtonian constant in m^3 kg^-1 s^-2. Observed SI bridge input.""" + +C_SI: float = 299792458.0 +"""Speed of light in m/s. Exact SI bridge input.""" + +SOLAR_MASS_KG: float = 1.98847e30 +"""Nominal solar mass in kg. Observed input for solar-limb validation.""" + +SOLAR_RADIUS_M: float = 6.957e8 +"""Nominal solar radius in m. Observed impact parameter for solar-limb validation.""" + +ARCSEC_PER_RADIAN: float = 206264.80624709636 +"""Arcseconds per radian.""" + N_COLORS: int = 3 """Number of color components (Ψₐ, a = 1,2,3).""" diff --git a/lfm/core/__init__.py b/lfm/core/__init__.py index 590c677..377cd27 100644 --- a/lfm/core/__init__.py +++ b/lfm/core/__init__.py @@ -19,5 +19,6 @@ """ from lfm.core.evolver import Evolver +from lfm.core.stencils import eigenvalue_19pt, laplacian_19pt -__all__ = ["Evolver"] +__all__ = ["Evolver", "eigenvalue_19pt", "laplacian_19pt"] diff --git a/lfm/core/stencils.py b/lfm/core/stencils.py index 7fdff0d..6c1715a 100644 --- a/lfm/core/stencils.py +++ b/lfm/core/stencils.py @@ -71,6 +71,34 @@ def laplacian_19pt(field: NDArray[np.floating]) -> NDArray[np.floating]: return STENCIL_FACE_WEIGHT * faces + STENCIL_EDGE_WEIGHT * edges + STENCIL_CENTER_WEIGHT * field +def eigenvalue_19pt( + kx: NDArray[np.floating], + ky: NDArray[np.floating], + kz: NDArray[np.floating], +) -> NDArray[np.floating]: + """Return the spectral eigenvalue of the 19-point stencil. + + The result matches :func:`laplacian_19pt` exactly on a periodic grid + with dx = 1. It is useful for FFT Poisson solves whose equilibrium + must be consistent with the same lattice operator used for evolution. + """ + face = ( + (2.0 * np.cos(kx) - 2.0) / 3.0 + + (2.0 * np.cos(ky) - 2.0) / 3.0 + + (2.0 * np.cos(kz) - 2.0) / 3.0 + ) + edge = ( + np.cos(kx + ky) + + np.cos(kx - ky) + + np.cos(kx + kz) + + np.cos(kx - kz) + + np.cos(ky + kz) + + np.cos(ky - kz) + - 6.0 + ) / 3.0 + return face + edge + + def laplacian_7pt(field: NDArray[np.floating]) -> NDArray[np.floating]: """Compute standard 7-point Laplacian on a 3D periodic grid. diff --git a/lfm/experiment/collision.py b/lfm/experiment/collision.py index 1c1b51a..ce194d2 100644 --- a/lfm/experiment/collision.py +++ b/lfm/experiment/collision.py @@ -489,7 +489,7 @@ def _build_collision_sim( np.clip(chi_template, 0.01, None, out=chi_template) dchi_template = chi_template - np.float32(geo.chi0) # Use chi_min as approximate eigenvalue (wave frequency inside well) - eigenvalue = float(max(chi_template.min(), 1.0)) + eigenvalue = max(float(chi_template.min()), 1.0) if verbose: print( f" Poisson-only ready: chi_min={chi_template.min():.4f} " diff --git a/lfm/fields/__init__.py b/lfm/fields/__init__.py index b8b82b7..2841fe2 100644 --- a/lfm/fields/__init__.py +++ b/lfm/fields/__init__.py @@ -26,10 +26,19 @@ from lfm.fields.boosted import boosted_soliton from lfm.fields.equilibrium import ( equilibrate_chi, + equilibrate_chi_19pt, equilibrate_from_fields, + equilibrate_from_fields_19pt, poisson_solve_fft, + poisson_solve_fft_19pt, +) +from lfm.fields.light import ( + planar_r1_light_packet, + r1_light_acceleration, + r1_light_step, + r1_vacuum_subtracted_potential, + spherical_phase_source, ) -from lfm.fields.light import spherical_phase_source from lfm.fields.random import seed_noise, uniform_chi from lfm.fields.soliton import gaussian_soliton, place_solitons, wave_kick from lfm.fields.spinor import ( @@ -45,8 +54,11 @@ "wave_kick", "boosted_soliton", "poisson_solve_fft", + "poisson_solve_fft_19pt", "equilibrate_chi", + "equilibrate_chi_19pt", "equilibrate_from_fields", + "equilibrate_from_fields_19pt", "seed_noise", "uniform_chi", "tetrahedral_positions", @@ -62,4 +74,8 @@ "apply_rotation_z", # light "spherical_phase_source", + "planar_r1_light_packet", + "r1_vacuum_subtracted_potential", + "r1_light_acceleration", + "r1_light_step", ] diff --git a/lfm/fields/equilibrium.py b/lfm/fields/equilibrium.py index 53aa4aa..11952f3 100644 --- a/lfm/fields/equilibrium.py +++ b/lfm/fields/equilibrium.py @@ -18,6 +18,7 @@ import numpy as np from lfm.constants import CHI0, KAPPA +from lfm.core.stencils import eigenvalue_19pt if TYPE_CHECKING: from numpy.typing import NDArray @@ -59,6 +60,31 @@ def poisson_solve_fft( return np.fft.irfftn(phi_hat, s=(N, N, N), axes=(0, 1, 2)).astype(np.float32) +def poisson_solve_fft_19pt( + source: NDArray[np.floating], + N: int | None = None, + dx: float = 1.0, +) -> NDArray[np.floating]: + """Solve L19(phi) = source with the exact 19-point stencil symbol.""" + if source.ndim != 3: + raise ValueError("source must have shape (N, N, N)") + N = int(N or source.shape[0]) + if source.shape != (N, N, N): + raise ValueError("source shape must match N") + out_dtype = np.float64 if np.dtype(source.dtype) == np.dtype(np.float64) else np.float32 + + src_hat = np.fft.rfftn(source.astype(np.float64)) + kx = np.fft.fftfreq(N) * 2.0 * np.pi + ky = np.fft.fftfreq(N) * 2.0 * np.pi + kz = np.fft.rfftfreq(N) * 2.0 * np.pi + KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing="ij") + lam = eigenvalue_19pt(KX, KY, KZ) / (dx * dx) + lam[0, 0, 0] = 1.0 + phi_hat = src_hat / lam + phi_hat[0, 0, 0] = 0.0 + return np.fft.irfftn(phi_hat, s=(N, N, N), axes=(0, 1, 2)).astype(out_dtype) + + def equilibrate_chi( psi_sq: NDArray[np.floating], chi0: float = CHI0, @@ -99,6 +125,23 @@ def equilibrate_chi( return chi +def equilibrate_chi_19pt( + psi_sq: NDArray[np.floating], + chi0: float = CHI0, + kappa: float = KAPPA, + e0_sq: float = 0.0, + boundary_mask: NDArray[np.bool_] | None = None, +) -> NDArray[np.floating]: + """Compute chi equilibrium with a 19-point-consistent Poisson solve.""" + N = psi_sq.shape[0] + rhs = kappa * (psi_sq - e0_sq) + delta_chi = poisson_solve_fft_19pt(rhs, N) + chi = (chi0 + delta_chi).astype(delta_chi.dtype, copy=False) + if boundary_mask is not None: + chi[boundary_mask] = chi0 + return chi + + def equilibrate_from_fields( psi_r: NDArray[np.float32], psi_i: NDArray[np.float32] | None = None, @@ -137,10 +180,33 @@ def equilibrate_from_fields( psi_sq = psi_sq + psi_i**2 elif psi_r.ndim == 4: # Multi-color: (n_colors, N, N, N) - psi_sq = np.sum(psi_r**2, axis=0) + psi_sq = np.sum(psi_r**2, axis=0).astype(np.float32) if psi_i is not None: - psi_sq = psi_sq + np.sum(psi_i**2, axis=0) + psi_sq = (psi_sq + np.sum(psi_i**2, axis=0)).astype(np.float32) else: raise ValueError(f"Unexpected psi_r shape: {psi_r.shape}") return equilibrate_chi(psi_sq, chi0, kappa, e0_sq, boundary_mask) + + +def equilibrate_from_fields_19pt( + psi_r: NDArray[np.float32], + psi_i: NDArray[np.float32] | None = None, + chi0: float = CHI0, + kappa: float = KAPPA, + e0_sq: float = 0.0, + boundary_mask: NDArray[np.bool_] | None = None, +) -> NDArray[np.float32]: + """Compute chi from fields using the 19-point-consistent Poisson solve.""" + if psi_r.ndim == 3: + psi_sq = psi_r**2 + if psi_i is not None: + psi_sq = psi_sq + psi_i**2 + elif psi_r.ndim == 4: + psi_sq = np.sum(psi_r**2, axis=0).astype(np.float32) + if psi_i is not None: + psi_sq = (psi_sq + np.sum(psi_i**2, axis=0)).astype(np.float32) + else: + raise ValueError(f"Unexpected psi_r shape: {psi_r.shape}") + + return equilibrate_chi_19pt(psi_sq, chi0, kappa, e0_sq, boundary_mask) diff --git a/lfm/fields/light.py b/lfm/fields/light.py index fffd736..1850871 100644 --- a/lfm/fields/light.py +++ b/lfm/fields/light.py @@ -51,10 +51,22 @@ import numpy as np +from lfm.constants import CHI0 +from lfm.core.stencils import laplacian_19pt + if TYPE_CHECKING: from numpy.typing import NDArray +def _runtime_float_dtype(*arrays: object) -> type[np.float32] | type[np.float64]: + """Use float64 only when a caller explicitly supplies float64 arrays.""" + for arr in arrays: + dtype = getattr(arr, "dtype", None) + if dtype is not None and np.dtype(dtype) == np.dtype(np.float64): + return np.float64 + return np.float32 + + def spherical_phase_source( N: int, center: tuple[float, float, float], @@ -141,3 +153,128 @@ def _shell(r_centre: float) -> NDArray[np.float64]: (shell_tm1 * cos_p).astype(np.float32), (shell_tm1 * sin_p).astype(np.float32), ) + + +def planar_r1_light_packet( + N: int, + center: tuple[float, float, float], + sigma: tuple[float, float, float], + carrier_k: float, + amplitude: float = 0.20, + dt: float = 0.32, + c_speed: float = 1.0, + axis: int = 0, + dtype: type[np.float32] | type[np.float64] = np.float32, +) -> tuple[NDArray[np.floating], NDArray[np.floating], NDArray[np.floating], NDArray[np.floating]]: + """Return a localized R1 complex packet moving in the +axis direction.""" + if axis not in (0, 1, 2): + raise ValueError("axis must be 0, 1, or 2") + out_dtype = np.dtype(dtype).type + if out_dtype not in (np.float32, np.float64): + raise ValueError("dtype must be np.float32 or np.float64") + + coords = np.arange(N, dtype=out_dtype) + grids = np.meshgrid(coords, coords, coords, indexing="ij") + center_arr = np.asarray(center, dtype=out_dtype) + sigma_arr = np.asarray(sigma, dtype=out_dtype) + if np.any(sigma_arr <= 0.0): + raise ValueError("sigma values must be positive") + + def _packet( + packet_center: NDArray[np.floating], + ) -> tuple[NDArray[np.floating], NDArray[np.floating]]: + radius_sq = np.zeros((N, N, N), dtype=out_dtype) + for idx, grid in enumerate(grids): + radius_sq += ((grid - packet_center[idx]) / sigma_arr[idx]) ** 2 + envelope = amplitude * np.exp(-0.5 * radius_sq) + phase = carrier_k * (grids[axis] - packet_center[axis]) + return ( + (envelope * np.cos(phase)).astype(out_dtype), + (envelope * np.sin(phase)).astype(out_dtype), + ) + + prev_center = center_arr.copy() + prev_center[axis] -= c_speed * dt + psi_r, psi_i = _packet(center_arr) + psi_r_prev, psi_i_prev = _packet(prev_center) + return psi_r, psi_i, psi_r_prev, psi_i_prev + + +def r1_vacuum_subtracted_potential( + chi: NDArray[np.floating], + chi0: float = CHI0, +) -> NDArray[np.floating]: + """Return chi^2 - chi0^2 for the massless R1 light perturbation.""" + out_dtype = _runtime_float_dtype(chi) + chi_f = chi.astype(out_dtype, copy=False) + return (chi_f * chi_f - out_dtype(chi0 * chi0)).astype(out_dtype) + + +def r1_light_acceleration( + psi_r: NDArray[np.floating], + psi_i: NDArray[np.floating], + chi: NDArray[np.floating] | None = None, + chi0: float = CHI0, + c_speed: float = 1.0, + vacuum_subtracted: bool = True, +) -> tuple[NDArray[np.floating], NDArray[np.floating]]: + """Compute the R1 light acceleration for one leapfrog update. + + With chi=None this is the flat U(1) massless phase/current sector. + With chi provided and vacuum_subtracted=True the potential is + chi^2 - chi0^2, so uniform vacuum remains massless while nonuniform + chi affects the full complex R1 field. + """ + out_dtype = _runtime_float_dtype(psi_r, psi_i, chi) + acc_r = (c_speed * c_speed * laplacian_19pt(psi_r)).astype(out_dtype) + acc_i = (c_speed * c_speed * laplacian_19pt(psi_i)).astype(out_dtype) + if chi is None: + return acc_r, acc_i + + chi_f = chi.astype(out_dtype, copy=False) + if vacuum_subtracted: + potential = r1_vacuum_subtracted_potential(chi_f, chi0) + else: + potential = (chi_f * chi_f).astype(out_dtype) + return ( + (acc_r - potential * psi_r).astype(out_dtype), + (acc_i - potential * psi_i).astype(out_dtype), + ) + + +def r1_light_step( + psi_r: NDArray[np.floating], + psi_i: NDArray[np.floating], + psi_r_prev: NDArray[np.floating], + psi_i_prev: NDArray[np.floating], + dt: float, + chi: NDArray[np.floating] | None = None, + chi0: float = CHI0, + c_speed: float = 1.0, + vacuum_subtracted: bool = True, + sponge: NDArray[np.floating] | None = None, +) -> tuple[NDArray[np.floating], NDArray[np.floating], NDArray[np.floating], NDArray[np.floating]]: + """Advance one leapfrog step for a massless R1 light packet.""" + out_dtype = _runtime_float_dtype(psi_r, psi_i, psi_r_prev, psi_i_prev, chi) + acc_r, acc_i = r1_light_acceleration( + psi_r, + psi_i, + chi=chi, + chi0=chi0, + c_speed=c_speed, + vacuum_subtracted=vacuum_subtracted, + ) + dt2 = out_dtype(dt * dt) + psi_r_next = (2.0 * psi_r - psi_r_prev + dt2 * acc_r).astype(out_dtype) + psi_i_next = (2.0 * psi_i - psi_i_prev + dt2 * acc_i).astype(out_dtype) + psi_r_prev_next = psi_r.astype(out_dtype, copy=True) + psi_i_prev_next = psi_i.astype(out_dtype, copy=True) + + if sponge is not None: + sponge_f = sponge.astype(out_dtype, copy=False) + psi_r_next *= sponge_f + psi_i_next *= sponge_f + psi_r_prev_next *= sponge_f + psi_i_prev_next *= sponge_f + + return psi_r_next, psi_i_next, psi_r_prev_next, psi_i_prev_next diff --git a/pyproject.toml b/pyproject.toml index 69856fb..f96d8c4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "lfm-physics" -version = "1.4.1" +version = "1.4.5" description = "Lattice Field Medium physics simulation library" readme = "README.md" license = "MIT" diff --git a/tests/test_fields.py b/tests/test_fields.py index 23edb33..19a2c4a 100644 --- a/tests/test_fields.py +++ b/tests/test_fields.py @@ -4,13 +4,20 @@ import pytest from lfm.constants import CHI0 +from lfm.core.stencils import laplacian_19pt from lfm.fields import ( equilibrate_chi, + equilibrate_chi_19pt, equilibrate_from_fields, gaussian_soliton, grid_positions, place_solitons, + planar_r1_light_packet, poisson_solve_fft, + poisson_solve_fft_19pt, + r1_light_acceleration, + r1_light_step, + r1_vacuum_subtracted_potential, seed_noise, sparse_positions, tetrahedral_positions, @@ -141,6 +148,25 @@ def test_roundtrip_spectral(self): laplacian = np.fft.irfftn(lap_hat, s=(N, N, N), axes=(0, 1, 2)) assert np.allclose(laplacian, source, atol=1e-3) + def test_roundtrip_19pt(self): + N = 16 + rng = np.random.default_rng(123) + source = rng.standard_normal((N, N, N)).astype(np.float32) + source -= np.mean(source) + phi = poisson_solve_fft_19pt(source, N) + assert phi.dtype == np.float32 + lap = laplacian_19pt(phi) + np.testing.assert_allclose(lap, source, atol=1e-4) + + def test_19pt_preserves_float64_for_small_sources(self): + N = 16 + source = np.zeros((N, N, N), dtype=np.float64) + source[N // 2, N // 2, N // 2] = 1e-9 + source -= np.mean(source) + phi = poisson_solve_fft_19pt(source, N) + assert phi.dtype == np.float64 + assert np.max(np.abs(phi)) > 0.0 + class TestEquilibrateChi: def test_zero_energy_gives_chi0(self): @@ -168,6 +194,92 @@ def test_boundary_mask(self): assert np.allclose(chi[0], CHI0) assert np.allclose(chi[-1], CHI0) + def test_19pt_energy_creates_well(self): + N = 24 + psi_sq = np.zeros((N, N, N), dtype=np.float32) + c = N // 2 + psi_sq[c - 1 : c + 2, c - 1 : c + 2, c - 1 : c + 2] = 20.0 + chi = equilibrate_chi_19pt(psi_sq) + assert chi[c, c, c] < CHI0 + + def test_19pt_preserves_float64_chi_perturbation(self): + N = 24 + psi_sq = np.zeros((N, N, N), dtype=np.float64) + c = N // 2 + psi_sq[c, c, c] = 1e-7 + chi = equilibrate_chi_19pt(psi_sq) + assert chi.dtype == np.float64 + assert np.min(chi) < CHI0 + + +class TestR1Light: + def test_vacuum_subtracted_potential_zero_in_uniform_vacuum(self): + chi = np.full((8, 8, 8), CHI0, dtype=np.float32) + pot = r1_vacuum_subtracted_potential(chi) + assert pot.dtype == np.float32 + np.testing.assert_allclose(pot, 0.0) + + def test_vacuum_subtracted_potential_preserves_float64_delta(self): + chi = np.full((8, 8, 8), CHI0, dtype=np.float64) + chi[4, 4, 4] -= 1e-9 + pot = r1_vacuum_subtracted_potential(chi) + assert pot.dtype == np.float64 + assert pot[4, 4, 4] < 0.0 + + def test_uniform_vacuum_acceleration_matches_flat_laplacian(self): + rng = np.random.default_rng(5) + psi_r = rng.standard_normal((8, 8, 8)).astype(np.float32) + psi_i = rng.standard_normal((8, 8, 8)).astype(np.float32) + chi = np.full((8, 8, 8), CHI0, dtype=np.float32) + acc_r, acc_i = r1_light_acceleration(psi_r, psi_i, chi=chi) + np.testing.assert_allclose(acc_r, laplacian_19pt(psi_r), atol=1e-6) + np.testing.assert_allclose(acc_i, laplacian_19pt(psi_i), atol=1e-6) + + def test_planar_packet_shapes_and_forward_current(self): + pr, pi, prp, pip = planar_r1_light_packet( + 16, + center=(5.0, 8.0, 8.0), + sigma=(2.0, 3.0, 3.0), + carrier_k=0.4, + ) + assert pr.shape == (16, 16, 16) + assert pi.dtype == np.float32 + assert prp.shape == pr.shape + assert pip.shape == pi.shape + dpr = 0.5 * (np.roll(pr, -1, axis=0) - np.roll(pr, 1, axis=0)) + dpi = 0.5 * (np.roll(pi, -1, axis=0) - np.roll(pi, 1, axis=0)) + current = np.sum(pr * dpi - pi * dpr) + assert current > 0.0 + + def test_planar_packet_float64_and_step_preserve_dtype(self): + pr, pi, prp, pip = planar_r1_light_packet( + 12, + center=(4.0, 6.0, 6.0), + sigma=(2.0, 2.5, 2.5), + carrier_k=0.3, + dtype=np.float64, + ) + chi = np.full((12, 12, 12), CHI0, dtype=np.float64) + nr, ni, npr, npi = r1_light_step(pr, pi, prp, pip, dt=0.1, chi=chi) + assert pr.dtype == np.float64 + assert nr.dtype == np.float64 + assert ni.dtype == np.float64 + assert npr.dtype == np.float64 + assert npi.dtype == np.float64 + + def test_r1_light_step_shapes(self): + pr, pi, prp, pip = planar_r1_light_packet( + 12, + center=(4.0, 6.0, 6.0), + sigma=(2.0, 2.5, 2.5), + carrier_k=0.3, + ) + nr, ni, npr, npi = r1_light_step(pr, pi, prp, pip, dt=0.1) + assert nr.shape == pr.shape + assert ni.shape == pi.shape + assert npr.shape == pr.shape + assert npi.shape == pi.shape + class TestEquilibrateFromFields: def test_real_field(self): @@ -295,5 +407,8 @@ def test_fields_importable_from_lfm(self): assert hasattr(lfm, "gaussian_soliton") assert hasattr(lfm, "equilibrate_chi") + assert hasattr(lfm, "equilibrate_chi_19pt") + assert hasattr(lfm, "planar_r1_light_packet") + assert hasattr(lfm, "r1_light_step") assert hasattr(lfm, "seed_noise") assert hasattr(lfm, "tetrahedral_positions") diff --git a/tests/test_physics_modules.py b/tests/test_physics_modules.py index 656eb7c..cfbd2e9 100644 --- a/tests/test_physics_modules.py +++ b/tests/test_physics_modules.py @@ -12,16 +12,23 @@ effective_metric_00, gravitational_potential, metric_perturbation, + metric_refractive_index, + op05_spherical_chi_deflection, schwarzschild_chi, + schwarzschild_radius_si, time_dilation_factor, ) from lfm.analysis.phase import ( charge_density, coulomb_interaction_energy, + noether_spatial_current, phase_coherence, + phase_current_energy_density, phase_field, + positive_noether_current, ) from lfm.config import SimulationConfig +from lfm.constants import SOLAR_MASS_KG, SOLAR_RADIUS_M from lfm.fields.boosted import boosted_soliton from lfm.sweep import sweep_2d @@ -77,6 +84,42 @@ def test_schwarzschild_chi_boundary(self): # At center → 0 (inside horizon) assert chi[16, 16, 16] == pytest.approx(0.0, abs=0.1) + def test_metric_refractive_index_increases_in_chi_well(self): + chi = np.array([19.0, 18.99], dtype=np.float64) + n_eff = metric_refractive_index(chi) + assert n_eff[0] == pytest.approx(1.0) + assert n_eff[1] > 1.0 + + def test_op05_solar_limb_deflection_recovers_arcsecond_scale(self): + result = op05_spherical_chi_deflection( + SOLAR_MASS_KG, + SOLAR_RADIUS_M, + x_extent_multiplier=500.0, + sample_count=20001, + ) + assert result["schwarzschild_radius_m"] == pytest.approx( + schwarzschild_radius_si(SOLAR_MASS_KG), + rel=1e-14, + ) + assert result["recovered_angle_arcsec"] == pytest.approx(1.751243281, rel=1e-3) + assert abs(result["comparator_relative_error"]) < 1e-4 + + def test_op05_deflection_is_linear_in_weak_lens_mass(self): + half = op05_spherical_chi_deflection( + 0.5 * SOLAR_MASS_KG, + SOLAR_RADIUS_M, + x_extent_multiplier=250.0, + sample_count=5001, + ) + full = op05_spherical_chi_deflection( + SOLAR_MASS_KG, + SOLAR_RADIUS_M, + x_extent_multiplier=250.0, + sample_count=5001, + ) + ratio = half["recovered_angle_arcsec"] / full["recovered_angle_arcsec"] + assert ratio == pytest.approx(0.5, rel=1e-5) + # ── Phase ──────────────────────────────────────────────────────────── @@ -146,6 +189,124 @@ def test_coulomb_opposite_phase(self): e_int = coulomb_interaction_energy(psi, zero, -psi, zero) assert e_int < 0 + def test_noether_spatial_current_plane_wave(self): + N = 16 + k = 2.0 * np.pi / N + x = np.arange(N, dtype=np.float32) + phase = k * x[:, None, None] + psi_r = np.broadcast_to(np.cos(phase), (N, N, N)).astype(np.float32) + psi_i = np.broadcast_to(np.sin(phase), (N, N, N)).astype(np.float32) + jx = noether_spatial_current(psi_r, psi_i, axis=0) + np.testing.assert_allclose(jx, np.sin(k), rtol=1e-5, atol=1e-6) + + def test_noether_spatial_current_preserves_float64(self): + N = 16 + k = 2.0 * np.pi / N + x = np.arange(N, dtype=np.float64) + phase = k * x[:, None, None] + psi_r = np.broadcast_to(np.cos(phase), (N, N, N)).astype(np.float64) + psi_i = np.broadcast_to(np.sin(phase), (N, N, N)).astype(np.float64) + jx = positive_noether_current(psi_r, psi_i, axis=0) + assert jx.dtype == np.float64 + np.testing.assert_allclose(jx, np.sin(k), rtol=1e-12, atol=1e-12) + + def test_positive_noether_current_clips_reverse_wave(self): + N = 16 + k = 2.0 * np.pi / N + x = np.arange(N, dtype=np.float32) + phase = -k * x[:, None, None] + psi_r = np.broadcast_to(np.cos(phase), (N, N, N)).astype(np.float32) + psi_i = np.broadcast_to(np.sin(phase), (N, N, N)).astype(np.float32) + jx = positive_noether_current(psi_r, psi_i, axis=0) + np.testing.assert_allclose(jx, 0.0, atol=1e-6) + + def test_phase_current_energy_static_uniform_zero(self): + psi_r = np.ones((8, 8, 8), dtype=np.float64) + psi_i = np.zeros((8, 8, 8), dtype=np.float64) + rho = phase_current_energy_density(psi_r, psi_i, psi_r, psi_i, dt=0.02) + assert rho.dtype == np.float64 + np.testing.assert_allclose(rho, 0.0, atol=1e-14) + + def test_phase_current_energy_global_phase_invariant(self): + N = 16 + dt = 0.02 + k = 2.0 * np.pi / N + omega = 0.35 + x = np.arange(N, dtype=np.float64) + phase = k * x[:, None, None] + psi_r = np.broadcast_to(np.cos(phase), (N, N, N)).astype(np.float64) + psi_i = np.broadcast_to(np.sin(phase), (N, N, N)).astype(np.float64) + prev_phase = phase - omega * dt + psi_r_prev = np.broadcast_to(np.cos(prev_phase), (N, N, N)).astype(np.float64) + psi_i_prev = np.broadcast_to(np.sin(prev_phase), (N, N, N)).astype(np.float64) + + rho = phase_current_energy_density(psi_r, psi_i, psi_r_prev, psi_i_prev, dt=dt) + + phi0 = 1.234 + psi_r_rot = np.cos(phi0) * psi_r - np.sin(phi0) * psi_i + psi_i_rot = np.sin(phi0) * psi_r + np.cos(phi0) * psi_i + psi_r_prev_rot = np.cos(phi0) * psi_r_prev - np.sin(phi0) * psi_i_prev + psi_i_prev_rot = np.sin(phi0) * psi_r_prev + np.cos(phi0) * psi_i_prev + rho_rot = phase_current_energy_density( + psi_r_rot, + psi_i_rot, + psi_r_prev_rot, + psi_i_prev_rot, + dt=dt, + ) + np.testing.assert_allclose(rho_rot, rho, rtol=1e-12, atol=1e-12) + + def test_phase_current_energy_sees_constant_amplitude_plane_wave(self): + N = 32 + dt = 0.01 + amp = 0.4 + k = 2.0 * np.pi / N + omega = k + x = np.arange(N, dtype=np.float64) + phase = k * x[:, None, None] + prev_phase = phase - omega * dt + psi_r = np.broadcast_to(amp * np.cos(phase), (N, N, N)).astype(np.float64) + psi_i = np.broadcast_to(amp * np.sin(phase), (N, N, N)).astype(np.float64) + psi_r_prev = np.broadcast_to(amp * np.cos(prev_phase), (N, N, N)).astype(np.float64) + psi_i_prev = np.broadcast_to(amp * np.sin(prev_phase), (N, N, N)).astype(np.float64) + + amp_sq = psi_r * psi_r + psi_i * psi_i + rho = phase_current_energy_density(psi_r, psi_i, psi_r_prev, psi_i_prev, dt=dt) + + np.testing.assert_allclose(amp_sq, amp * amp, rtol=1e-14, atol=1e-14) + assert float(np.mean(rho)) > 0.0 + assert float(np.std(rho)) < 1e-12 + + def test_phase_current_energy_scales_with_carrier(self): + N = 64 + dt = 0.005 + amp = 0.25 + x = np.arange(N, dtype=np.float64) + + def mean_rho(mode: int) -> float: + k = 2.0 * np.pi * mode / N + phase = k * x[:, None, None] + prev_phase = phase - k * dt + psi_r = np.broadcast_to(amp * np.cos(phase), (N, N, N)).astype(np.float64) + psi_i = np.broadcast_to(amp * np.sin(phase), (N, N, N)).astype(np.float64) + psi_r_prev = np.broadcast_to(amp * np.cos(prev_phase), (N, N, N)).astype(np.float64) + psi_i_prev = np.broadcast_to(amp * np.sin(prev_phase), (N, N, N)).astype(np.float64) + return float( + np.mean( + phase_current_energy_density( + psi_r, + psi_i, + psi_r_prev, + psi_i_prev, + dt=dt, + ) + ) + ) + + rho_1 = mean_rho(1) + rho_2 = mean_rho(2) + assert rho_2 / rho_1 == pytest.approx(4.0, rel=0.08) + # ── Angular Momentum ───────────────────────────────────────────────── diff --git a/tests/test_stencils.py b/tests/test_stencils.py index 21b430e..884cbdb 100644 --- a/tests/test_stencils.py +++ b/tests/test_stencils.py @@ -2,7 +2,7 @@ import numpy as np -from lfm.core.stencils import laplacian_7pt, laplacian_19pt +from lfm.core.stencils import eigenvalue_19pt, laplacian_7pt, laplacian_19pt class TestLaplacian19pt: @@ -55,6 +55,19 @@ def test_output_shape(self): lap = laplacian_19pt(field) assert lap.shape == field.shape + def test_eigenvalue_matches_single_fourier_mode(self): + N = 16 + coords = np.arange(N, dtype=np.float64) + X, Y, Z = np.meshgrid(coords, coords, coords, indexing="ij") + kx_i, ky_i, kz_i = 2, 1, 0 + kx = 2.0 * np.pi * kx_i / N + ky = 2.0 * np.pi * ky_i / N + kz = 2.0 * np.pi * kz_i / N + field = np.cos(kx * X + ky * Y + kz * Z) + lap = laplacian_19pt(field) + lam = float(eigenvalue_19pt(np.array(kx), np.array(ky), np.array(kz))) + np.testing.assert_allclose(lap, lam * field, atol=1e-12) + class TestLaplacian7pt: def test_constant_field_zero(self): diff --git a/tests/validation/test_string_tension.py b/tests/validation/test_string_tension.py index 9f022c7..a3f43c9 100644 --- a/tests/validation/test_string_tension.py +++ b/tests/validation/test_string_tension.py @@ -12,6 +12,7 @@ from __future__ import annotations import numpy as np +import pytest from lfm import Simulation from lfm.config_presets import full_physics @@ -64,6 +65,7 @@ def _make_pair(sep: int) -> Simulation: class TestStringTension: """Tube energy between colored quarks should grow with separation.""" + @pytest.mark.timeout(300) def test_string_tension(self) -> None: """Tube energy should increase monotonically and yield positive σ.