/* * Focused regression for the explicit DP5(4) adaptive geodesic core and its * RayPool integration. * * The production stepper (src/geodesic.c), the production RHS, the analytic * Minkowski and Schwarzschild backends, the observer builder and the * production RayPool scheduler are exercised directly. The single-trace * geodesic core, `geodesic_advance_past_ray`, and the batch pool's activation, * OMP bulk advance and resume-state handling are covered here. * * Both analytic providers define spacetime_create_default, so they are * textually included with a renamed provider symbol; COMMON_SOURCES excludes * them from the link line. */ #define spacetime_create_default spacetime_create_default_minkowski_adaptive #include "../src/spacetime_minkowski.c" #undef spacetime_create_default #define spacetime_create_default spacetime_create_default_schwarzschild_adaptive #include "../src/spacetime_schwarzschild.c" #undef spacetime_create_default #include "asymptotic.h" #include "asymptotic_schwarzschild.h" #include "geodesic.h" #include "observer.h" #include "ray.h" #include "spacetime.h" #include #include #include #include #include static int failures = 0; #define CHECK(condition, message) \ do { \ if (!(condition)) { \ fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \ ++failures; \ } \ } while (0) static double dot3(const double a[3], const double b[3]) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } static int invert3(const double g[3][3], double inv[3][3]) { const double d = g[0][0] * (g[1][1] * g[2][2] - g[1][2] * g[2][1]) - g[0][1] * (g[1][0] * g[2][2] - g[1][2] * g[2][0]) + g[0][2] * (g[1][0] * g[2][1] - g[1][1] * g[2][0]); if (!isfinite(d) || fabs(d) < 1e-300) return -1; inv[0][0] = (g[1][1] * g[2][2] - g[1][2] * g[2][1]) / d; inv[0][1] = (g[0][2] * g[2][1] - g[0][1] * g[2][2]) / d; inv[0][2] = (g[0][1] * g[1][2] - g[0][2] * g[1][1]) / d; inv[1][0] = (g[1][2] * g[2][0] - g[1][0] * g[2][2]) / d; inv[1][1] = (g[0][0] * g[2][2] - g[0][2] * g[2][0]) / d; inv[1][2] = (g[0][2] * g[1][0] - g[0][0] * g[1][2]) / d; inv[2][0] = (g[1][0] * g[2][1] - g[1][1] * g[2][0]) / d; inv[2][1] = (g[0][1] * g[2][0] - g[0][0] * g[2][1]) / d; inv[2][2] = (g[0][0] * g[1][1] - g[0][1] * g[1][0]) / d; return 0; } static double null_residual(const MetricData *m, const double Pi[3]) { double inv[3][3]; if (invert3(m->gamma, inv)) return NAN; double value = 0.0; for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) value += inv[i][j] * Pi[i] * Pi[j]; return value; } /* Stationary KS Killing energy reference (oracle cross-check only). */ static double killing_energy(const MetricData *m, const double Pi[3], double log_alpha_p0) { return exp(log_alpha_p0) * (m->alpha - dot3(m->beta, Pi)); } static RayEndpoint blank_endpoint(void) { RayEndpoint out; memset(&out, 0, sizeof out); out.magnification = 1.0; out.end_id = SPACETIME_END_NONE; out.outcome = RAY_OUTCOME_INCOMPLETE; out.reason = RAY_REASON_NONE; out.threshold_value = NAN; out.stop_coordinate_time = NAN; return out; } static GeodesicTraceConfig dp_config(double tol, double initial_step, double lookback) { GeodesicTraceConfig c; memset(&c, 0, sizeof c); c.stepper = GEODESIC_STEPPER_DP54; c.coordinate_time_step = initial_step; c.max_steps = 1000000u; c.threshold = (ThresholdPolicy){.kind = THRESHOLD_DISABLED, .value = 0.0, .policy_version = 0}; c.atol_x = tol; c.atol_Pi = tol; c.atol_L = tol; c.rtol = tol; c.min_step = 1e-12; c.max_step = 1e6; c.consecutive_rejection_limit = 1000u; c.max_lookback_time = lookback; return c; } static ObserverState flat_observer_at(const double p[3]) { ObserverState o; memset(&o, 0, sizeof o); o.coordinate_position[0] = p[0]; o.coordinate_position[1] = p[1]; o.coordinate_position[2] = p[2]; o.tetrad[0][0] = 1.0; o.tetrad[1][1] = 1.0; o.tetrad[2][2] = 1.0; o.tetrad[3][3] = 1.0; return o; } static MetricData flat_metric(void) { MetricData m; memset(&m, 0, sizeof m); m.alpha = 1.0; m.gamma[0][0] = m.gamma[1][1] = m.gamma[2][2] = 1.0; return m; } /* ------------------------------------------------------------------ */ /* Synthetic fixture: flat / linear-potential curved metric with a */ /* configurable spatial domain, invalid-metric and unavailable-time */ /* injection windows. Declares a single fixed Minkowski end. */ /* ------------------------------------------------------------------ */ typedef struct { int has_domain; /* enable the [domain_lo, domain_hi] wall check */ double domain_lo, domain_hi; SpacetimePointStatus bad_status; /* injection status (OK disables) */ double bad_lo, bad_hi; double radius; /* escape worldtube radius */ } FixtureContext; static SpacetimePointStatus fixture_eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric) { const FixtureContext *c = source->context; if (c->bad_status != SPACETIME_POINT_OK && t >= c->bad_lo && t <= c->bad_hi) return c->bad_status; if (c->has_domain && x[0] >= c->domain_lo && x[0] <= c->domain_hi) return SPACETIME_POINT_OUT_OF_DOMAIN; *metric = flat_metric(); return SPACETIME_POINT_OK; } static SpacetimeRayStatus fixture_classify(const SpacetimeSource *source, double t, const double x[3]) { (void)source; (void)t; (void)x; return SPACETIME_RAY_ACTIVE; } static size_t fixture_end_count(const SpacetimeSource *source) { (void)source; return 1; } static int fixture_end(const SpacetimeSource *source, size_t index, SpacetimeAsymptoticEnd *out) { (void)source; if (index != 0) return -1; *out = (SpacetimeAsymptoticEnd){ .end_id = 0, .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, .mass = 0.0, .frame_origin = {0.0, 0.0, 0.0}, .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; return 0; } static int fixture_worldtube(const SpacetimeSource *source, SpacetimeEndId end_id, double t, SpacetimeEscapeWorldtubeSample *out) { const FixtureContext *c = source->context; (void)t; if (end_id != 0) return -1; *out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0}, .velocity = {0.0, 0.0, 0.0}, .radius = c->radius, .radius_rate = 0.0, .velocity_constant = 1, .valid = 1}; return 0; } static void fixture_destroy(SpacetimeSource *source) { source->context = NULL; source->ops = NULL; } static const SpacetimeOps fixture_ops = { .eval = fixture_eval, .classify = fixture_classify, .asymptotic_end_count = fixture_end_count, .asymptotic_end = fixture_end, .escape_worldtube_sample = fixture_worldtube, .destroy = fixture_destroy, }; /* ------------------------------------------------------------------ */ /* 1. Minkowski finite interval: DP accuracy, translation, slabs, count */ /* ------------------------------------------------------------------ */ static void test_minkowski_finite_interval(void) { SpacetimeSource source = {0}; CHECK(spacetime_create_minkowski(&source, 1.0e9) == 0, "create minkowski"); const double position[3] = {3.0, -4.0, 5.0}; const double direction[3] = {0.36, 0.48, 0.8}; const ObserverState observer = flat_observer_at(position); MetricData metric; CHECK(spacetime_eval(&source, 0.0, position, &metric) == SPACETIME_POINT_OK, "minkowski metric"); GeodesicRayState state; CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, &state) == 0, "minkowski init"); const GeodesicTraceConfig config = dp_config(1e-9, 0.5, 1.0e6); state.next_step = config.coordinate_time_step; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -3.0, &slab) == 0, "mink slab"); RayEndpoint endpoint = blank_endpoint(); const GeodesicAdvanceResult result = geodesic_advance_past_ray(slab, &state, -2.0, &config, &endpoint); CHECK(result == GEODESIC_ADVANCE_ACTIVE, "finite interval stays active"); /* Flat RHS: x(t) = x0 + (t - t0) * Pi and L constant. */ CHECK(fabs(state.coordinate_time + 2.0) < 1e-12, "reached interval end"); for (int i = 0; i < 3; ++i) { const double expected = position[i] - 2.0 * state.Pi[i]; CHECK(fabs(state.x[i] - expected) < 1e-9, "minkowski translation exact"); } CHECK(fabs(state.log_alpha_p0) < 1e-12, "flat L stays zero"); CHECK(state.rhs_evaluations == 7u * ((unsigned long)state.steps + (unsigned long)state.rejected_steps), "actual RHS count equals full DP trials"); CHECK(state.steps > 0, "accepted at least one step"); spacetime_free_slab(slab); /* Translation independence: the same boost shifted by a constant vector * produces the same state shifted by that vector. */ { const double shift[3] = {100.0, -200.0, 300.0}; double shifted[3]; for (int i = 0; i < 3; ++i) shifted[i] = position[i] + shift[i]; const ObserverState observer2 = flat_observer_at(shifted); MetricData metric2; CHECK(spacetime_eval(&source, 0.0, shifted, &metric2) == SPACETIME_POINT_OK, "shifted metric"); GeodesicRayState state2; CHECK(geodesic_initialize_past_ray_metric(&metric2, &observer2, direction, &state2) == 0, "shifted init"); state2.next_step = config.coordinate_time_step; MetricSlab *slab2 = NULL; CHECK(spacetime_load_slab(&source, 0.0, -3.0, &slab2) == 0, "shifted slab"); RayEndpoint endpoint2 = blank_endpoint(); (void)geodesic_advance_past_ray(slab2, &state2, -2.0, &config, &endpoint2); for (int i = 0; i < 3; ++i) CHECK(fabs((state2.x[i] - shift[i]) - state.x[i]) < 1e-9, "translation invariance"); spacetime_free_slab(slab2); } /* Cross-slab continuity: two partial slabs must match one long slab. */ { GeodesicRayState split; CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, &split) == 0, "split init"); split.next_step = config.coordinate_time_step; MetricSlab *first = NULL, *second = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.5, &first) == 0, "first slab"); RayEndpoint e1 = blank_endpoint(); const GeodesicAdvanceResult r1 = geodesic_advance_past_ray(first, &split, -1.0, &config, &e1); spacetime_free_slab(first); CHECK(r1 == GEODESIC_ADVANCE_ACTIVE, "first partial slab active"); CHECK(spacetime_load_slab(&source, split.coordinate_time, -3.0, &second) == 0, "second slab"); RayEndpoint e2 = blank_endpoint(); const GeodesicAdvanceResult r2 = geodesic_advance_past_ray(second, &split, -2.0, &config, &e2); spacetime_free_slab(second); CHECK(r2 == GEODESIC_ADVANCE_ACTIVE, "second partial slab active"); for (int i = 0; i < 3; ++i) CHECK(fabs(split.x[i] - state.x[i]) < 1e-7, "cross-slab continuity x"); CHECK(fabs(split.log_alpha_p0 - state.log_alpha_p0) < 1e-9, "cross-slab continuity L"); } /* A tiny slab-left interval must still accept a boundary-limited step even * below min_step, without permanently depressing the next suggestion. */ { GeodesicRayState tiny; CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, &tiny) == 0, "tiny init"); GeodesicTraceConfig tiny_config = dp_config(1e-9, 1.0, 1.0e6); tiny_config.min_step = 0.25; tiny.next_step = tiny_config.coordinate_time_step; MetricSlab *s = NULL; CHECK(spacetime_load_slab(&source, 0.0, -0.1, &s) == 0, "tiny slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(s, &tiny, -1e-6, &tiny_config, &out); spacetime_free_slab(s); CHECK(r == GEODESIC_ADVANCE_ACTIVE, "tiny boundary step active"); CHECK(fabs(tiny.coordinate_time + 1e-6) < 1e-18, "tiny boundary reached"); CHECK(tiny.next_step == tiny_config.coordinate_time_step, "tiny boundary kept the original proposal"); } /* Invalid DP configuration must be rejected explicitly, not defaulted. */ { GeodesicRayState bad_state; CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, &bad_state) == 0, "bad config init"); MetricSlab *s = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &s) == 0, "bad slab"); GeodesicTraceConfig bad = dp_config(1e-9, 0.5, 1.0e6); bad.atol_x = 0.0; RayEndpoint out = blank_endpoint(); CHECK(geodesic_advance_past_ray(s, &bad_state, -0.5, &bad, &out) == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_INVALID_STEPPER_CONFIG, "zero atol rejected"); bad = dp_config(1e-9, 0.5, 1.0e6); bad.max_lookback_time = 0.0; CHECK(geodesic_advance_past_ray(s, &bad_state, -0.5, &bad, &out) == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_INVALID_STEPPER_CONFIG, "zero lookback rejected"); bad = dp_config(1e-9, 0.5, 1.0e6); bad.min_step = 2.0; bad.max_step = 1.0; CHECK(geodesic_advance_past_ray(s, &bad_state, -0.5, &bad, &out) == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_INVALID_STEPPER_CONFIG, "inverted step bounds rejected"); bad = dp_config(1e-9, 0.5, 1.0e6); bad.consecutive_rejection_limit = 0; CHECK(geodesic_advance_past_ray(s, &bad_state, -0.5, &bad, &out) == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_INVALID_STEPPER_CONFIG, "zero reject limit rejected"); spacetime_free_slab(s); } spacetime_destroy(&source); } /* ------------------------------------------------------------------ */ /* 2. Schwarzschild radial branches, interior observer, residuals */ /* ------------------------------------------------------------------ */ static void radial_static_observer(const MetricData *metric, double r0, ObserverState *out) { *out = (ObserverState){.coordinate_time = 0.0, .coordinate_position = {r0, 0.0, 0.0}}; const double alpha = metric->alpha; out->tetrad[0][0] = 1.0 / alpha; for (int i = 0; i < 3; ++i) out->tetrad[0][i + 1] = -metric->beta[i] / alpha; out->tetrad[1][1] = 1.0 / sqrt(metric->gamma[0][0]); out->tetrad[2][2] = 1.0; out->tetrad[3][3] = 1.0; } static int trace_radial_dp(const SpacetimeSource *source, const ObserverState *o, double direction, double duration, const GeodesicTraceConfig *config, GeodesicRayState *state) { MetricData metric; if (spacetime_eval(source, o->coordinate_time, o->coordinate_position, &metric) != SPACETIME_POINT_OK) return -1; const double n[3] = {direction, 0.0, 0.0}; if (geodesic_initialize_past_ray_metric(&metric, o, n, state)) return -1; state->next_step = config->coordinate_time_step; MetricSlab *slab = NULL; if (spacetime_load_slab(source, 0.0, -duration - 1.0, &slab)) return -1; RayEndpoint endpoint = blank_endpoint(); const GeodesicAdvanceResult result = geodesic_advance_past_ray(slab, state, -duration, config, &endpoint); spacetime_free_slab(slab); return result == GEODESIC_ADVANCE_FAILED ? -1 : 0; } static void test_schwarzschild_radial(void) { SpacetimeSource source = {0}; CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); const double r0 = 10.0; MetricData metric; CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &metric) == SPACETIME_POINT_OK, "r0 metric"); ObserverState observer; radial_static_observer(&metric, r0, &observer); const GeodesicTraceConfig config = dp_config(1e-7, 0.5, 1.0e6); GeodesicRayState outward; CHECK(trace_radial_dp(&source, &observer, 1.0, 2.0, &config, &outward) == 0, "outward trace"); const double s_out = -outward.coordinate_time; CHECK(fabs(outward.x[0] - r0 - s_out) < 1e-4, "outward closed-form invariant"); MetricData out_metric; CHECK(spacetime_eval(&source, outward.coordinate_time, outward.x, &out_metric) == SPACETIME_POINT_OK, "outward end metric"); CHECK(fabs(null_residual(&out_metric, outward.Pi) - 1.0) < 1e-6, "outward null constraint"); GeodesicRayState inward; CHECK(trace_radial_dp(&source, &observer, -1.0, 2.0, &config, &inward) == 0, "inward trace"); const double s_in = -inward.coordinate_time; const double c0 = (r0 - 2.0) + 4.0 * log(r0 - 2.0); const double c1 = (inward.x[0] - 2.0) + 4.0 * log(inward.x[0] - 2.0) + s_in; CHECK(inward.x[0] > 2.0, "inward stays outside horizon"); CHECK(inward.x[0] < r0, "inward decreases r"); CHECK(fabs(c1 - c0) < 1e-4, "inward closed-form invariant"); MetricData in_metric; CHECK(spacetime_eval(&source, inward.coordinate_time, inward.x, &in_metric) == SPACETIME_POINT_OK, "inward end metric"); CHECK(fabs(null_residual(&in_metric, inward.Pi) - 1.0) < 1e-6, "inward null constraint"); /* Killing energy conservation along a radial ray. */ { GeodesicRayState start; MetricData m0 = metric; CHECK(geodesic_initialize_past_ray_metric( &m0, &observer, (double[]){1.0, 0.0, 0.0}, &start) == 0, "killing init"); GeodesicRayState end; CHECK(trace_radial_dp(&source, &observer, 1.0, 2.0, &config, &end) == 0, "killing trace"); MetricData m1; CHECK(spacetime_eval(&source, end.coordinate_time, end.x, &m1) == SPACETIME_POINT_OK, "killing end metric"); const double ek_start = killing_energy(&m0, start.Pi, start.log_alpha_p0); const double ek_end = killing_energy(&m1, end.Pi, end.log_alpha_p0); CHECK(fabs(ek_end / ek_start - 1.0) < 1e-5, "Killing energy conserved along DP ray"); } /* Interior r = 1.5 free-fall-from-rest-at-infinity observer: legal * timelike observer, null constraint along the integrated ray. */ { const double r = 1.5; const double f = 1.0 - 2.0 / r; const double ur = -sqrt(2.0 / r); const double uks = 1.0 / f + (2.0 / (r - 2.0)) * ur; const double vx = ur / uks; ObserverCamera camera = {.position = {r, 0.0, 0.0}, .velocity = {vx, 0.0, 0.0}, .look_ra_deg = 0.0, .look_dec_deg = 0.0, .roll_deg = 0.0}; MetricData m; ObserverState inner; CHECK(spacetime_eval(&source, 0.0, camera.position, &m) == SPACETIME_POINT_OK && observer_from_coordinate_camera(&m, &camera, &inner, NULL) == OBSERVER_BUILD_OK, "interior free-fall observer"); GeodesicRayState state; CHECK(trace_radial_dp(&source, &inner, 1.0, 1.0, &config, &state) == 0, "interior trace"); MetricData end_metric; CHECK(spacetime_eval(&source, state.coordinate_time, state.x, &end_metric) == SPACETIME_POINT_OK, "interior end metric"); CHECK(fabs(null_residual(&end_metric, state.Pi) - 1.0) < 1e-5, "interior null constraint"); } /* Tolerance tightening must improve (or preserve) the radial residual. */ { const GeodesicTraceConfig loose = dp_config(1e-5, 0.5, 1.0e6); const GeodesicTraceConfig tight = dp_config(1e-9, 0.5, 1.0e6); GeodesicRayState sl, st; CHECK(trace_radial_dp(&source, &observer, 1.0, 2.0, &loose, &sl) == 0 && trace_radial_dp(&source, &observer, 1.0, 2.0, &tight, &st) == 0, "tolerance scan traces"); const double slow = fabs(sl.x[0] - r0 + sl.coordinate_time); const double stight = fabs(st.x[0] - r0 + st.coordinate_time); CHECK(stight <= slow + 1e-9, "tighter tolerance no worse residual"); CHECK(stight < 1e-6, "tight residual small"); } spacetime_destroy(&source); } /* ------------------------------------------------------------------ */ /* 3. Over-large initial step: rejection without accepted-state damage */ /* ------------------------------------------------------------------ */ static void test_over_large_initial_step(void) { SpacetimeSource source = {0}; CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); const double r0 = 10.0; MetricData metric; CHECK(spacetime_eval(&source, 0.0, (double[]){r0, 0.0, 0.0}, &metric) == SPACETIME_POINT_OK, "metric"); ObserverState observer; radial_static_observer(&metric, r0, &observer); const GeodesicTraceConfig small = dp_config(1e-9, 0.05, 1.0e6); GeodesicRayState reference; CHECK(trace_radial_dp(&source, &observer, -1.0, 2.0, &small, &reference) == 0, "reference small-step trace"); GeodesicTraceConfig big = dp_config(1e-9, 4.0, 1.0e6); GeodesicRayState large; CHECK(trace_radial_dp(&source, &observer, -1.0, 2.0, &big, &large) == 0, "large initial step trace"); CHECK(large.rejected_steps > 0, "large initial step was rejected"); /* Comparing final states is the available evidence that rejected trials did * not pollute the accepted trajectory; it is not a direct internal proof. */ CHECK(fabs(large.coordinate_time - reference.coordinate_time) < 1e-12, "large-step trace reaches the same time"); for (int i = 0; i < 3; ++i) CHECK(fabs(large.x[i] - reference.x[i]) < 1e-4, "large-step result matches small-step reference"); CHECK(fabs(large.log_alpha_p0 - reference.log_alpha_p0) < 1e-5, "large-step L matches reference"); spacetime_destroy(&source); } /* ------------------------------------------------------------------ */ /* 4. Synthetic fixture: retryable stages, fatal points, hmin/maxreject */ /* ------------------------------------------------------------------ */ /* Build a manual flat state at x = (x0, 0, 0) moving in +x in the past * (Pi = -xhat), used by the fixture tests. */ static void fixture_state(double x0, GeodesicRayState *state) { memset(state, 0, sizeof *state); state->x[0] = x0; state->Pi[0] = -1.0; } static void test_fixture_stage_retry(void) { /* Step-size-dependent retryable stage failure: a thin forbidden wall lies * across the outward path. The initial trial places a stage inside the * wall (OUT_OF_DOMAIN), is rejected and shrunk; a smaller/reshaped step then * jumps over the wall and the integration completes. */ FixtureContext context; memset(&context, 0, sizeof context); context.has_domain = 1; context.domain_lo = 0.9; context.domain_hi = 0.95; context.radius = 1.0e9; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(0.0, &state); state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-8, 1.15, 1.0e6); state.next_step = config.coordinate_time_step; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -3.0, &slab) == 0, "wall slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult result = geodesic_advance_past_ray(slab, &state, -2.0, &config, &out); spacetime_free_slab(slab); CHECK(result == GEODESIC_ADVANCE_ACTIVE, "wall overshoot eventually active"); CHECK(fabs(state.coordinate_time + 2.0) < 1e-9, "wall reached target"); CHECK(state.rejected_steps > 0, "wall overshoot caused a rejection"); CHECK(fabs(state.x[0] - 2.0) < 1e-9, "wall integration stayed on the ray"); } static void test_fixture_fatal_and_bounds(void) { /* The initial accepted state itself is invalid: direct specific failure, * one RHS call, no shrink. */ { FixtureContext context; memset(&context, 0, sizeof context); context.bad_status = SPACETIME_POINT_INVALID_METRIC; context.bad_lo = -1.0e30; context.bad_hi = 1.0e30; context.radius = 1.0e9; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(0.0, &state); GeodesicTraceConfig config = dp_config(1e-8, 0.5, 1.0e6); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "invalid slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -0.5, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_INVALID_METRIC, "initial invalid metric is a direct specific failure"); CHECK(state.rhs_evaluations == 1u, "initial failure used exactly one RHS"); CHECK(state.rejected_steps == 0u, "initial failure did not reject-retry"); } /* A later stage hits TIME_UNAVAILABLE: direct specific failure, no shrink. */ { FixtureContext context; memset(&context, 0, sizeof context); context.bad_status = SPACETIME_POINT_TIME_UNAVAILABLE; context.bad_lo = -0.3; context.bad_hi = -0.01; context.radius = 1.0e9; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(0.0, &state); GeodesicTraceConfig config = dp_config(1e-8, 0.5, 1.0e6); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "unavail slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -0.5, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_TIME_RANGE_EXHAUSTED, "stage TIME_UNAVAILABLE is a direct specific failure"); CHECK(state.rhs_evaluations == 2u, "stage failure was not retried with repeated shrink"); CHECK(state.rejected_steps == 0u, "stage TIME_UNAVAILABLE not a reject"); } /* A later stage hits INTERNAL_ERROR: direct protocol failure. */ { FixtureContext context; memset(&context, 0, sizeof context); context.bad_status = SPACETIME_POINT_INTERNAL_ERROR; context.bad_lo = -0.3; context.bad_hi = -0.01; context.radius = 1.0e9; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(0.0, &state); GeodesicTraceConfig config = dp_config(1e-8, 0.5, 1.0e6); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "internal slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -0.5, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_METRIC_INTERNAL_ERROR, "stage INTERNAL_ERROR is a direct protocol failure"); CHECK(state.rhs_evaluations == 2u, "internal error not retried"); } /* Retryable stage OUT_OF_DOMAIN that never clears: hmin and maxreject. */ { FixtureContext context; memset(&context, 0, sizeof context); context.bad_status = SPACETIME_POINT_OUT_OF_DOMAIN; context.bad_lo = -1.0e30; context.bad_hi = -0.001; /* stage 1 (t = -h/5) is always in the window */ context.radius = 1.0e9; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(0.0, &state); GeodesicTraceConfig hmin = dp_config(1e-8, 1.0, 1.0e6); hmin.min_step = 0.1; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -10.0, &slab) == 0, "hmin slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -10.0, &hmin, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_MIN_STEP_REACHED, "step below hmin is an integration error"); CHECK(state.rejected_steps >= 1u, "hmin path rejected before failing"); FixtureContext context2; memset(&context2, 0, sizeof context2); context2.bad_status = SPACETIME_POINT_OUT_OF_DOMAIN; context2.bad_lo = -1.0e30; context2.bad_hi = -0.001; context2.radius = 1.0e9; SpacetimeSource source2 = {.ops = &fixture_ops, .context = &context2}; GeodesicRayState state2; fixture_state(0.0, &state2); GeodesicTraceConfig maxreject = dp_config(1e-8, 1.0, 1.0e6); maxreject.min_step = 1e-12; maxreject.consecutive_rejection_limit = 2; MetricSlab *slab2 = NULL; CHECK(spacetime_load_slab(&source2, 0.0, -10.0, &slab2) == 0, "maxreject slab"); RayEndpoint out2 = blank_endpoint(); const GeodesicAdvanceResult r2 = geodesic_advance_past_ray(slab2, &state2, -10.0, &maxreject, &out2); spacetime_free_slab(slab2); CHECK(r2 == GEODESIC_ADVANCE_FAILED && out2.reason == RAY_REASON_OUT_OF_DOMAIN, "consecutive rejection limit reports the specific stage reason"); CHECK(state2.rejected_steps == 2u, "maxreject counted exactly two rejects"); CHECK(state2.coordinate_time == 0.0, "rejected trials did not advance the accepted time"); } } /* ------------------------------------------------------------------ */ /* 5. Budget, explicit lookback, history exhaustion, resume */ /* ------------------------------------------------------------------ */ static GeodesicRayState continuation_from(const RayEndpoint *e) { GeodesicRayState s; memset(&s, 0, sizeof s); s.coordinate_time = e->stop_coordinate_time; for (int i = 0; i < 3; ++i) { s.x[i] = e->final_x[i]; s.Pi[i] = e->final_Pi[i]; } s.log_alpha_p0 = e->final_log_alpha_p0; s.log_alpha_p0_0 = e->final_log_alpha_p0_0; s.steps = e->accepted_steps; s.integration_start_time = e->integration_start_time; s.next_step = e->next_step; s.rejected_steps = e->rejected_steps; s.rhs_evaluations = e->rhs_evaluations; s.previous_rejected = e->previous_rejected; return s; } static void test_budget_lookback_and_resume(void) { SpacetimeSource source = {0}; CHECK(spacetime_create_minkowski(&source, 10.0) == 0, "create minkowski"); const ObserverState observer = flat_observer_at((double[]){0.0, 0.0, 0.0}); const double direction[3] = {1.0, 0.0, 0.0}; MetricData metric; CHECK(spacetime_eval(&source, 0.0, (double[]){0.0, 0.0, 0.0}, &metric) == SPACETIME_POINT_OK, "metric"); /* Accepted-step budget exhaustion -> UNRESOLVED/BUDGET_EXHAUSTED. */ { GeodesicRayState state; CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, &state) == 0, "budget init"); GeodesicTraceConfig config = dp_config(1e-9, 0.5, 1.0e6); config.max_steps = 3; /* Keep three accepted steps well inside the r = 10 escape sphere. */ config.max_step = 1.0; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -50.0, &slab) == 0, "budget slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -50.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_UNRESOLVED && out.reason == RAY_REASON_BUDGET_EXHAUSTED, "max_steps yields UNRESOLVED/BUDGET_EXHAUSTED"); CHECK(state.steps == 3u, "budget consumed exactly max_steps"); } /* Explicit lookback budget exhaustion -> UNRESOLVED, not history. */ { GeodesicRayState state; CHECK(geodesic_initialize_past_ray_metric(&metric, &observer, direction, &state) == 0, "lookback init"); GeodesicTraceConfig config = dp_config(1e-9, 0.5, 0.25); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -50.0, &slab) == 0, "lookback slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -50.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_UNRESOLVED && out.reason == RAY_REASON_BUDGET_EXHAUSTED, "explicit lookback yields UNRESOLVED/BUDGET_EXHAUSTED"); CHECK(fabs(state.coordinate_time + 0.25) < 1e-9, "lookback stopped at the explicit budget"); } /* A source slab data time hole is TIME_RANGE_EXHAUSTED, not UNRESOLVED. */ { FixtureContext context; memset(&context, 0, sizeof context); context.bad_status = SPACETIME_POINT_TIME_UNAVAILABLE; context.bad_lo = -1.0e30; context.bad_hi = -0.3; context.radius = 1.0e9; SpacetimeSource fs = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(0.0, &state); GeodesicTraceConfig config = dp_config(1e-9, 0.25, 1.0e6); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&fs, 0.0, -1.0, &slab) == 0, "hole slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.outcome == RAY_OUTCOME_INCOMPLETE && out.reason == RAY_REASON_TIME_RANGE_EXHAUSTED, "slab data hole is a specific INCOMPLETE reason, not UNRESOLVED"); } /* Resume preserves L0, start time, counters and next step; an explicitly * larger lookback budget continues the trace to escape without replay. */ { GeodesicTraceConfig first = dp_config(1e-9, 0.5, 1.0); const RayEndpoint partial = geodesic_trace_past(&source, &observer, direction, &first); CHECK(partial.outcome == RAY_OUTCOME_UNRESOLVED && partial.reason == RAY_REASON_BUDGET_EXHAUSTED, "first lookback trace is unresolved"); const double l0 = partial.final_log_alpha_p0_0; const double start = partial.integration_start_time; const unsigned long rhs_before = partial.rhs_evaluations; GeodesicRayState continuation = continuation_from(&partial); GeodesicTraceConfig retry = dp_config(1e-9, 0.5, 100.0); const RayEndpoint resumed = geodesic_trace_past_from_state(&source, &continuation, &retry); CHECK(resumed.outcome == RAY_OUTCOME_ESCAPED, "resumed trace escapes"); CHECK(resumed.final_log_alpha_p0_0 == l0, "resume preserves L0"); CHECK(resumed.integration_start_time == start, "resume preserves integration start time"); CHECK(resumed.rhs_evaluations > rhs_before, "resume accumulates RHS cost without replaying from scratch"); CHECK(resumed.accepted_steps >= partial.accepted_steps, "resume keeps accepted step count"); } spacetime_destroy(&source); } /* ------------------------------------------------------------------ */ /* 6. Escape event localization: flat analytic + strong-field compare */ /* ------------------------------------------------------------------ */ static void test_flat_escape_localization(void) { FixtureContext context; memset(&context, 0, sizeof context); context.has_domain = 0; context.radius = 5.0; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; /* Start inside the worldtube and integrate outward so the crossing is * produced by the advance loop itself, not by the pre-route. */ GeodesicRayState state; fixture_state(2.0, &state); GeodesicTraceConfig config = dp_config(1e-10, 0.5, 1.0e6); state.next_step = config.coordinate_time_step; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -20.0, &slab) == 0, "escape slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -20.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED && out.end_id == 0, "flat DP escape detected by the advance loop"); CHECK(fabs(out.n_infinity[0] - 1.0) < 1e-6 && fabs(out.n_infinity[1]) < 1e-6 && fabs(out.n_infinity[2]) < 1e-6, "flat DP escape direction"); const double expected_crossing = -(context.radius - 2.0); CHECK(fabs(out.stop_coordinate_time - expected_crossing) < 1e-4, "flat DP crossing time matches analytic root"); } /* Event localization with a translated time origin. At t0 = 1e12 the local * coordinate-time ULP is ~1.2e-4, so the crossing cannot be resolved below * that; the old relative epsilon (= 1.0 at 1e12) skipped the subintegration * and fabricated an off-trajectory endpoint. The endpoint must stay on the * flat ray x(t) = 2 + (t0 - t) to a few ULP. */ static void check_flat_crossing_at(double t0, double position_tol, double time_tol) { FixtureContext context; memset(&context, 0, sizeof context); context.radius = 5.0; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(2.0, &state); state.coordinate_time = t0; state.integration_start_time = t0; GeodesicTraceConfig config = dp_config(1e-10, 0.5, 100.0); config.max_step = 0.5; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, t0, t0 - 20.0, &slab) == 0, "translated slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, t0 - 20.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "translated flat escape"); CHECK(fabs(out.final_x[0] - 5.0) <= position_tol, "translated crossing position matches analytic root"); CHECK(fabs((t0 - out.stop_coordinate_time) - 3.0) <= time_tol, "translated crossing time matches analytic root"); CHECK(fabs(out.final_x[0] - 2.0 - (t0 - out.stop_coordinate_time)) <= 2.0 * position_tol, "translated endpoint lies on the flat ray"); } static void test_event_time_translation(void) { check_flat_crossing_at(0.0, 1e-6, 1e-6); check_flat_crossing_at(1e12, 2e-3, 2e-3); } /* A TIME_UNAVAILABLE window that the accepted main trial does not hit but a * bisection subintegration does must surface as INCOMPLETE / * TIME_RANGE_EXHAUSTED (the metric-integration reason), not as a worldtube * protocol error, with the last trusted state restored and the failed * subintegration RHS calls accumulated. */ static void test_event_subintegration_time_hole(void) { FixtureContext context; memset(&context, 0, sizeof context); context.radius = 2.7; context.bad_status = SPACETIME_POINT_TIME_UNAVAILABLE; /* The dense escape bracket re-localizes on the accepted trajectory; the * bisection reaches t = -0.595 here, so the injected window must cover it. * Only the subintegration sees the hole, not the accepted main trial. */ context.bad_lo = -0.596; context.bad_hi = -0.594; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(2.0, &state); GeodesicTraceConfig config = dp_config(1e-10, 1.0, 100.0); config.max_step = 1.0; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -2.0, &slab) == 0, "event hole slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -2.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.outcome == RAY_OUTCOME_INCOMPLETE, "event subintegration hole fails as incomplete"); CHECK(out.reason == RAY_REASON_TIME_RANGE_EXHAUSTED, "event subintegration propagates TIME_RANGE_EXHAUSTED"); CHECK(out.reason != RAY_REASON_PROTOCOL_ERROR, "event integration reason is not a worldtube protocol error"); CHECK(state.coordinate_time == 0.0 && out.stop_coordinate_time == 0.0, "event hole restores the last trusted time"); CHECK(out.final_x[0] == 2.0, "event hole restores the last trusted position"); CHECK(state.rhs_evaluations > 7u, "event hole RHS includes the failed subintegration calls"); } static void test_schwarzschild_dp_escape_matches_finish(void) { SpacetimeSource source = {0}; CHECK(spacetime_create_schwarzschild_ks(&source, 1.0, 256.0) == 0, "create schwarzschild"); SpacetimeAsymptoticEnd end; CHECK(spacetime_asymptotic_end(&source, 0, &end) == 0, "end"); SchwarzschildCanonical canonical = {.end_id = 0, .t = 0.0, .rho = 256.0, .rhat = {0.8, 0.6, 0.0}, .Lhat = {0.0, 0.0, 1.0}, .beta = 5.0, .energy = 1.0, .radial_sign = 1}; double x[3], Pi[3], log_alpha_p0; CHECK(asymptotic_schwarzschild_state_from_canonical( &end, &canonical, x, Pi, &log_alpha_p0) == 0, "canonical state"); double n_analytic[3], freq_analytic; CHECK(asymptotic_schwarzschild_finish(&end, &canonical, n_analytic, &freq_analytic) == 0, "analytic finish"); GeodesicRayState state; memset(&state, 0, sizeof state); state.x[0] = x[0]; state.x[1] = x[1]; state.x[2] = x[2]; state.Pi[0] = Pi[0]; state.Pi[1] = Pi[1]; state.Pi[2] = Pi[2]; state.log_alpha_p0 = log_alpha_p0; state.log_alpha_p0_0 = log_alpha_p0; GeodesicTraceConfig config = dp_config(1e-9, 0.5, 1.0e6); state.next_step = config.coordinate_time_step; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0e6, &slab) == 0, "far slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0e6, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "strong-field DP integration escapes"); const double angle = acos(fmax( -1.0, fmin(1.0, dot3(n_analytic, out.n_infinity)))); CHECK(angle < 1e-3, "DP escape direction matches analytic finish"); CHECK(fabs(out.frequency_ratio / freq_analytic - 1.0) < 1e-3, "DP escape frequency matches analytic finish"); spacetime_destroy(&source); } /* ------------------------------------------------------------------ */ /* 7. Legacy RK4 selection is unchanged (zero-initialized config) */ /* ------------------------------------------------------------------ */ static void test_rk4_compatibility(void) { SpacetimeSource source = {0}; CHECK(spacetime_create_minkowski(&source, 10.0) == 0, "create minkowski"); const ObserverState observer = flat_observer_at((double[]){0.0, 0.0, 0.0}); GeodesicTraceConfig config; memset(&config, 0, sizeof config); config.coordinate_time_step = 0.25; config.max_steps = 100; config.threshold.kind = THRESHOLD_DISABLED; CHECK(config.stepper == GEODESIC_STEPPER_RK4, "zero-initialized config selects RK4"); const RayEndpoint out = geodesic_trace_past(&source, &observer, (double[]){1.0, 0.0, 0.0}, &config); CHECK(out.outcome == RAY_OUTCOME_ESCAPED, "RK4 escape still works"); CHECK(fabs(out.frequency_ratio - 1.0) < 1e-12, "RK4 flat frequency"); spacetime_destroy(&source); } /* ------------------------------------------------------------------ */ /* 8. RayPool batch scheduler: OMP determinism, single-trace agreement, */ /* activation/control state and cross-slab reject accounting */ /* ------------------------------------------------------------------ */ /* Sweep a pool from `top` downward until every ray terminates or the bottom * bound is reached, activating pending rays in each slab. */ static int pool_sweep(RayPool *pool, const SpacetimeSource *source, const GeodesicTraceConfig *config, double top, double bottom, double slab_duration) { ray_pool_preroute(pool, source); double slab_hi = top; while (ray_pool_has_live(pool) && slab_hi > bottom) { const double slab_lo = fmax(slab_hi - slab_duration, bottom); MetricSlab *slab = NULL; if (spacetime_load_slab(source, slab_hi, slab_lo, &slab)) return -1; ray_pool_activate_in_time_range(pool, slab); ray_pool_advance_active(pool, slab, config); spacetime_free_slab(slab); slab_hi = slab_lo; } return 0; } static void test_ray_pool_batch_matches_single(void) { FixtureContext context; memset(&context, 0, sizeof context); context.radius = 50.0; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; const ObserverState observer = flat_observer_at((double[]){2.0, 0.0, 0.0}); const double camera_directions[4][3] = {{1.0, 0.0, 0.0}, {-1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; const GeodesicTraceConfig config = dp_config(1e-9, 0.5, 1.0e6); RayEndpoint by_threads[2][4]; for (int threads = 1; threads <= 2; ++threads) { omp_set_num_threads(threads); RayPool pool; CHECK(ray_pool_init(&pool, 4) == 0, "batch pool init"); for (int i = 0; i < 4; ++i) CHECK(ray_pool_append(&pool, &observer, camera_directions[i], 0, (size_t)i) == 0, "batch pool append"); ray_pool_preroute(&pool, &source); for (int i = 0; i < 4; ++i) { CHECK(pool.continuation[i] == 0, "new pool ray is not a continuation"); CHECK(pool.observer[i] == &observer, "preroute keeps the observer"); CHECK(pool.integration_start_time[i] == pool.activate_t[i], "adaptive window starts at the activation time"); CHECK(pool.next_step[i] == 0.0, "first trial step is deferred to the trace config"); CHECK(pool.rejected_steps[i] == 0 && pool.rhs_evaluations[i] == 0, "new pool ray starts with zero adaptive cost"); } CHECK(pool_sweep(&pool, &source, &config, 0.0, -1000.0, 10.0) == 0, "batch pool sweep"); for (int i = 0; i < 4; ++i) { CHECK(pool.status[i] == RAY_POOL_TERMINATED, "batch ray terminated"); CHECK(pool.endpoint[i].outcome == RAY_OUTCOME_ESCAPED, "batch ray escapes"); CHECK(pool.steps[i] > 0, "batch ray accepted steps"); by_threads[threads - 1][i] = pool.endpoint[i]; } ray_pool_destroy(&pool); } omp_set_num_threads(1); for (int i = 0; i < 4; ++i) { for (int axis = 0; axis < 3; ++axis) CHECK(fabs(by_threads[0][i].n_infinity[axis] - by_threads[1][i].n_infinity[axis]) < 1e-12, "OMP1/OMP2 escape direction identical"); CHECK(fabs(by_threads[0][i].stop_coordinate_time - by_threads[1][i].stop_coordinate_time) < 1e-12, "OMP1/OMP2 crossing time identical"); } for (int i = 0; i < 4; ++i) { const RayEndpoint single = geodesic_trace_past(&source, &observer, camera_directions[i], &config); CHECK(single.outcome == RAY_OUTCOME_ESCAPED, "single trace escapes"); for (int axis = 0; axis < 3; ++axis) CHECK(fabs(single.n_infinity[axis] - by_threads[0][i].n_infinity[axis]) < 1e-3, "pool matches single-trace escape direction"); CHECK(fabs(single.stop_coordinate_time - by_threads[0][i].stop_coordinate_time) < 1e-2, "pool matches single-trace crossing time"); } } static void test_ray_pool_continuation_state(void) { FixtureContext context; memset(&context, 0, sizeof context); context.radius = 1.0e9; /* no escape inside the traced window */ SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; const ObserverState observer = flat_observer_at((double[]){0.0, 0.0, 0.0}); GeodesicTraceConfig first = dp_config(1e-9, 0.5, 1.0e6); first.max_steps = 3; first.max_step = 1.0; const RayEndpoint partial = geodesic_trace_past( &source, &observer, (double[]){-1.0, 0.0, 0.0}, &first); CHECK(partial.outcome == RAY_OUTCOME_UNRESOLVED && partial.reason == RAY_REASON_BUDGET_EXHAUSTED, "partial trace is budget-unresolved"); GeodesicRayState state = continuation_from(&partial); const double l0 = state.log_alpha_p0_0; const double start = state.integration_start_time; const double t_before = state.coordinate_time; const unsigned long rhs_before = state.rhs_evaluations; RayPool pool; CHECK(ray_pool_init(&pool, 1) == 0, "continuation pool init"); CHECK(ray_pool_append_continuation_state(&pool, 0, 0, &state, 6, partial.lookback_limit) == 0, "append continuation state"); CHECK(pool.continuation[0] == 1, "continuation is flagged"); CHECK(pool.observer[0] == NULL, "continuation carries no camera observer"); CHECK(pool.integration_start_time[0] == start && pool.log_alpha_p0_0[0] == l0 && pool.rejected_steps[0] == state.rejected_steps && pool.rhs_evaluations[0] == rhs_before && pool.next_step[0] == state.next_step, "continuation copies the full adaptive state"); ray_pool_preroute(&pool, &source); CHECK(pool.status[0] == RAY_POOL_PENDING && pool.observer[0] == NULL, "continuation skips preroute"); CHECK(pool.integration_start_time[0] == start, "preroute does not touch continuation control state"); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -10.0, &slab) == 0, "continuation slab"); ray_pool_activate_in_time_range(&pool, slab); CHECK(pool.status[0] == RAY_POOL_ACTIVE && pool.steps[0] == state.steps, "activation preserves the continuation step count"); ray_pool_advance_active(&pool, slab, &first); spacetime_free_slab(slab); CHECK(pool.steps[0] == 6u, "continuation reached the new step budget"); CHECK(pool.endpoint[0].outcome == RAY_OUTCOME_UNRESOLVED, "continuation is still budget-unresolved"); CHECK(pool.t[0] < t_before, "continuation moved further into the past"); CHECK(pool.integration_start_time[0] == start && pool.log_alpha_p0_0[0] == l0, "continuation preserves window start and L0"); CHECK(pool.rhs_evaluations[0] > rhs_before, "continuation accumulates RHS without replay"); CHECK(pool.rhs_evaluations[0] == 7ul * ((unsigned long)pool.steps[0] + (unsigned long)pool.rejected_steps[0]), "pool RHS count matches full DP trials across the resume"); ray_pool_destroy(&pool); } static void test_ray_pool_cross_slab_reject(void) { /* A thin forbidden wall lies across the outward path: the first trial is * rejected and shrunk. Splitting the sweep across slabs verifies that the * accumulated rejection/RHS cost and the accepted state survive both the * slab boundary and the pool gather/scatter. */ FixtureContext context; memset(&context, 0, sizeof context); context.has_domain = 1; context.domain_lo = 0.9; context.domain_hi = 0.95; context.radius = 1.0e9; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; const ObserverState observer = flat_observer_at((double[]){0.0, 0.0, 0.0}); RayPool pool; CHECK(ray_pool_init(&pool, 1) == 0, "reject pool init"); /* A flat camera direction n gives Pi = -n, so n = +xhat sends the past ray * toward increasing x, across the forbidden wall. */ CHECK(ray_pool_append(&pool, &observer, (double[]){1.0, 0.0, 0.0}, 0, 0) == 0, "reject pool append"); const GeodesicTraceConfig config = dp_config(1e-8, 1.15, 1.0e6); ray_pool_preroute(&pool, &source); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.5, &slab) == 0, "reject slab 1"); ray_pool_activate_in_time_range(&pool, slab); ray_pool_advance_active(&pool, slab, &config); spacetime_free_slab(slab); const unsigned int rejected_first = pool.rejected_steps[0]; const unsigned long rhs_first = pool.rhs_evaluations[0]; CHECK(rejected_first > 0, "wall overshoot rejected in the first slab"); CHECK(rhs_first > 0, "first slab spent RHS calls"); CHECK(spacetime_load_slab(&source, -1.5, -2.0, &slab) == 0, "reject slab 2"); ray_pool_activate_in_time_range(&pool, slab); ray_pool_advance_active(&pool, slab, &config); spacetime_free_slab(slab); CHECK(fabs(pool.t[0] + 2.0) < 1e-9, "reject pool reached the target time"); CHECK(fabs(pool.x0[0] - 2.0) < 1e-9, "reject pool stayed on the ray"); CHECK(pool.steps[0] > 0, "reject run accepted steps"); CHECK(pool.rejected_steps[0] >= rejected_first, "rejection count survives the slab boundary"); CHECK(pool.rhs_evaluations[0] > rhs_first, "RHS cost accumulates across slabs"); /* Every accepted trial costs 7 RHS; a rejected trial costs at least the * initial evaluation plus the failing stage. */ CHECK(pool.rhs_evaluations[0] >= 7ul * (unsigned long)pool.steps[0] + 2ul * (unsigned long)pool.rejected_steps[0], "cross-slab reject/RHS accounting is consistent"); ray_pool_destroy(&pool); } /* ------------------------------------------------------------------ */ /* 9. Dense event layer: polynomial roots, multi-end order, segments, */ /* endpoint/tangent roots, and threshold-versus-escape ordering. */ /* ------------------------------------------------------------------ */ /* Focused hook into the production dense-event polynomial root finder. */ extern int geodesic_event_polynomial_exit(const double *coeff, int degree, double *theta, double *lo, double *hi); extern int geodesic_event_worldtube_dense_probe(const MetricSlab *slab, const GeodesicTraceConfig *config, SpacetimeEndId end_id, GeodesicRayState *state, double left, double theta, double *dense_value, double *actual_value); static void test_event_polynomial_roots(void) { double theta, lo, hi; /* inside -> outside -> inside: the earliest exit at the FIRST root. The * production root isolator must not return the later re-entry. */ { const double c[3] = {-0.1875, 1.0, -1.0}; /* -theta^2+theta-0.1875 */ CHECK(geodesic_event_polynomial_exit(c, 2, &theta, &lo, &hi) == 1, "in-out-in polynomial has an exit"); CHECK(fabs(theta - 0.25) < 1e-9, "earliest exit is the first root"); CHECK(lo <= 1e-12 && hi > 0.25 && hi < 0.75, "exit bracket stops before the re-entry"); } /* Tangent (even multiplicity, no sign change) is not an escape. */ { const double c[3] = {-0.25, 1.0, -1.0}; /* -(theta-0.5)^2 */ CHECK(geodesic_event_polynomial_exit(c, 2, &theta, &lo, &hi) == 0, "tangent touch is not an escape"); } /* A root exactly at the step endpoint is left for the following step. */ { const double c[2] = {-1.0, 1.0}; /* theta - 1 */ CHECK(geodesic_event_polynomial_exit(c, 1, &theta, &lo, &hi) == 0, "root at theta == 1 is not an in-step exit"); } /* A boundary start that immediately moves outside is an exit at theta 0. */ { const double c[2] = {0.0, 1.0}; /* theta */ CHECK(geodesic_event_polynomial_exit(c, 1, &theta, &lo, &hi) == 1, "boundary start moving outside is an exit"); CHECK(theta <= 1e-12, "boundary exit is at theta 0"); } /* outside -> inside (entry) is not an escape. */ { const double c[2] = {1.0, -2.0}; /* 1 - 2 theta */ CHECK(geodesic_event_polynomial_exit(c, 1, &theta, &lo, &hi) == 0, "entry is not an escape"); } } /* Flat worldtube fixture with constant-velocity ends, optional motion * segments (piecewise-constant velocity/radius_rate) and several ends. */ #define EVENT_ENDS 2 typedef struct { size_t end_count; SpacetimeEndId id[EVENT_ENDS]; double center[EVENT_ENDS][3]; double velocity[EVENT_ENDS][3]; double radius_a[EVENT_ENDS]; double radius_rate_a[EVENT_ENDS]; double radius_b[EVENT_ENDS]; double radius_rate_b[EVENT_ENDS]; double segment_t[EVENT_ENDS]; double t_ref[EVENT_ENDS]; /* center reference time for the moving end */ int has_segment[EVENT_ENDS]; int velocity_constant[EVENT_ENDS]; } EventContext; static SpacetimePointStatus event_eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric) { (void)source; (void)t; (void)x; *metric = flat_metric(); return SPACETIME_POINT_OK; } static SpacetimeRayStatus event_classify(const SpacetimeSource *source, double t, const double x[3]) { (void)source; (void)t; (void)x; return SPACETIME_RAY_ACTIVE; } static size_t event_end_count(const SpacetimeSource *source) { return ((const EventContext *)source->context)->end_count; } static int event_end(const SpacetimeSource *source, size_t index, SpacetimeAsymptoticEnd *out) { const EventContext *c = source->context; if (index >= c->end_count) return -1; *out = (SpacetimeAsymptoticEnd){ .end_id = c->id[index], .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, .mass = 0.0, .frame_origin = {0.0, 0.0, 0.0}, .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; return 0; } static int event_worldtube(const SpacetimeSource *source, SpacetimeEndId end_id, double t, SpacetimeEscapeWorldtubeSample *out) { const EventContext *c = source->context; for (size_t i = 0; i < c->end_count; ++i) { if (c->id[i] != end_id) continue; double radius = c->radius_a[i]; double radius_rate = c->radius_rate_a[i]; if (c->has_segment[i] && t < c->segment_t[i]) { radius = c->radius_b[i] + c->radius_rate_b[i] * (t - c->segment_t[i]); radius_rate = c->radius_rate_b[i]; } const double dt = t - c->t_ref[i]; *out = (SpacetimeEscapeWorldtubeSample){ .center = {c->center[i][0] + c->velocity[i][0] * dt, c->center[i][1] + c->velocity[i][1] * dt, c->center[i][2] + c->velocity[i][2] * dt}, .velocity = {c->velocity[i][0], c->velocity[i][1], c->velocity[i][2]}, .radius = radius, .radius_rate = radius_rate, .velocity_constant = c->velocity_constant[i], .valid = 1}; return 0; } return -1; } static double event_next_segment(const SpacetimeSource *source, SpacetimeEndId end_id, double t) { const EventContext *c = source->context; for (size_t i = 0; i < c->end_count; ++i) if (c->id[i] == end_id && c->has_segment[i] && t > c->segment_t[i]) return c->segment_t[i]; return NAN; } static void event_destroy(SpacetimeSource *source) { source->context = NULL; source->ops = NULL; } static const SpacetimeOps event_ops = { .eval = event_eval, .classify = event_classify, .asymptotic_end_count = event_end_count, .asymptotic_end = event_end, .escape_worldtube_sample = event_worldtube, .escape_worldtube_next_segment = event_next_segment, .destroy = event_destroy, }; /* Two ends at the same center: end id 0 is the smaller (nearer) sphere, so it * is exited first no matter which descriptor position it occupies. */ static void test_event_two_ends_order(void) { for (int swap = 0; swap < 2; ++swap) { EventContext c; memset(&c, 0, sizeof c); c.end_count = 2; c.velocity_constant[0] = c.velocity_constant[1] = 1; const SpacetimeEndId near_id = 3, far_id = 7; const double near_radius = 3.0, far_radius = 5.0; if (!swap) { c.id[0] = near_id; c.radius_a[0] = near_radius; c.id[1] = far_id; c.radius_a[1] = far_radius; } else { c.id[0] = far_id; c.radius_a[0] = far_radius; c.id[1] = near_id; c.radius_a[1] = near_radius; } SpacetimeSource source = {.ops = &event_ops, .context = &c}; GeodesicRayState state; fixture_state(2.0, &state); GeodesicTraceConfig config = dp_config(1e-10, 10.0, 1.0e6); config.max_step = 10.0; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -20.0, &slab) == 0, "two-end slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -20.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "two-end escape terminates"); CHECK(out.end_id == near_id, "earlier end wins regardless of descriptor order"); CHECK(fabs(out.stop_coordinate_time + 1.0) < 1e-6, "two-end crossing is the nearest surface"); } } /* Piecewise-constant worldtube: the step must be clamped to the segment * boundary or the constant-velocity polynomial extrapolates the wrong radius * and predicts the wrong crossing. */ static void test_event_segmented_worldtube(void) { EventContext c; memset(&c, 0, sizeof c); c.end_count = 1; c.id[0] = 0; c.velocity_constant[0] = 1; c.radius_a[0] = 5.0; c.radius_rate_a[0] = 0.0; c.radius_b[0] = 5.0; c.radius_rate_b[0] = 1.0; c.segment_t[0] = -2.0; c.has_segment[0] = 1; SpacetimeSource source = {.ops = &event_ops, .context = &c}; GeodesicRayState state; fixture_state(2.0, &state); GeodesicTraceConfig config = dp_config(1e-10, 4.0, 1.0e6); config.max_step = 4.0; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -20.0, &slab) == 0, "segment slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -20.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "segmented worldtube escape"); CHECK(fabs(out.stop_coordinate_time + 2.5) < 1e-6, "segmented crossing uses the second segment"); } /* Arbitrary accelerated worldtubes are explicitly unsupported rather than * routed with a coarse sign search that can miss a narrow first entry. */ static void test_event_nonconstant_unsupported(void) { EventContext c; memset(&c, 0, sizeof c); c.end_count = 1; c.id[0] = 0; c.radius_a[0] = 5.0; c.velocity_constant[0] = 0; /* arbitrary acceleration */ SpacetimeSource source = {.ops = &event_ops, .context = &c}; /* Pre-route from outside must report UNSUPPORTED, not a fabricated miss. */ const ObserverState outside = flat_observer_at((double[]){10.0, 0.0, 0.0}); AsymptoticRoute route; CHECK(asymptotic_route_camera(&source, &outside, (double[]){-1.0, 0.0, 0.0}, &route) == ASYMPTOTIC_UNSUPPORTED, "nonconstant preroute is unsupported"); /* An inside state that reaches the advance loop is rejected there too. */ GeodesicRayState state; fixture_state(1.0, &state); GeodesicTraceConfig config = dp_config(1e-9, 0.5, 1.0e6); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -5.0, &slab) == 0, "nonconstant slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -5.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.outcome == RAY_OUTCOME_INCOMPLETE && out.reason == RAY_REASON_UNSUPPORTED, "nonconstant advance is an explicit unsupported failure"); } /* Synthetic curved metric with alpha = 1 - k x. The Eulerian energy * L = ln(alpha p^0) grows along the past-directed ray, so a threshold crossing * and a worldtube escape can be placed in the same accepted step. */ typedef struct { double k; double radius; int has_bad; double bad_lo, bad_hi; } AlphaContext; static SpacetimePointStatus alpha_eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric) { const AlphaContext *c = source->context; if (c->has_bad && t >= c->bad_lo && t <= c->bad_hi) return SPACETIME_POINT_OUT_OF_DOMAIN; const double alpha = 1.0 - c->k * x[0]; if (!(alpha > 0.0)) return SPACETIME_POINT_INVALID_METRIC; *metric = flat_metric(); metric->alpha = alpha; metric->d_alpha[0] = -c->k; return SPACETIME_POINT_OK; } static size_t alpha_end_count(const SpacetimeSource *source) { (void)source; return 1; } static int alpha_end(const SpacetimeSource *source, size_t index, SpacetimeAsymptoticEnd *out) { (void)source; if (index != 0) return -1; *out = (SpacetimeAsymptoticEnd){ .end_id = 0, .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, .mass = 0.0, .frame_origin = {0.0, 0.0, 0.0}, .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; return 0; } static int alpha_worldtube(const SpacetimeSource *source, SpacetimeEndId end_id, double t, SpacetimeEscapeWorldtubeSample *out) { const AlphaContext *c = source->context; (void)t; if (end_id != 0) return -1; *out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0}, .velocity = {0.0, 0.0, 0.0}, .radius = c->radius, .radius_rate = 0.0, .velocity_constant = 1, .valid = 1}; return 0; } static const SpacetimeOps alpha_ops = { .eval = alpha_eval, .classify = event_classify, .asymptotic_end_count = alpha_end_count, .asymptotic_end = alpha_end, .escape_worldtube_sample = alpha_worldtube, .destroy = event_destroy, }; static void alpha_state(double x0, double k, GeodesicRayState *state) { memset(state, 0, sizeof *state); state->x[0] = x0; state->Pi[0] = -1.0; state->log_alpha_p0 = log(1.0 - k * x0); state->log_alpha_p0_0 = state->log_alpha_p0; } /* Threshold and escape in the same accepted step: the earlier event wins. */ static void test_event_threshold_order(void) { struct { double radius; RayOutcome outcome; const char *label; } cases[2] = {{0.5, RAY_OUTCOME_DARK, "threshold first"}, {0.3, RAY_OUTCOME_ESCAPED, "escape first"}}; for (int i = 0; i < 2; ++i) { AlphaContext c; memset(&c, 0, sizeof c); c.k = 1.0; c.radius = cases[i].radius; SpacetimeSource source = {.ops = &alpha_ops, .context = &c}; GeodesicRayState state; alpha_state(0.2, c.k, &state); GeodesicTraceConfig config = dp_config(1e-4, 1.0, 1.0e6); config.max_step = 1.0; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = 0.35, .policy_version = 1}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -3.0, &slab) == 0, "threshold slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -3.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == cases[i].outcome, cases[i].label); if (out.outcome == RAY_OUTCOME_DARK) { CHECK(out.accepted_steps == 1u, "threshold and escape were compared within one accepted step"); CHECK(out.reason == RAY_REASON_REDSHIFT_LIMIT, "dark reason"); CHECK(out.threshold_value >= 0.35 && out.threshold_value < 0.35 + 0.05, "dark records the actual threshold value"); } } } /* A trial whose stages cross the threshold but is rejected must not produce a * dark terminal. Here the long first trial hits a forbidden window, is * rejected and shrunk to a step that stays below the threshold; with a * one-step budget the ray ends UNRESOLVED, not DARK. */ static void test_event_rejected_stage_threshold(void) { AlphaContext c; memset(&c, 0, sizeof c); c.k = 1.0; c.radius = 10.0; c.has_bad = 1; c.bad_lo = -0.81; c.bad_hi = -0.79; SpacetimeSource source = {.ops = &alpha_ops, .context = &c}; GeodesicRayState state; alpha_state(0.2, c.k, &state); GeodesicTraceConfig config = dp_config(1e-9, 1.0, 1.0e6); config.max_step = 1.0; config.max_steps = 1; /* Threshold crossing near s = 0.8, well beyond the shrunk accepted step. */ config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = 0.8, .policy_version = 1}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -2.0, &slab) == 0, "reject threshold slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -2.0, &config, &out); spacetime_free_slab(slab); CHECK(state.rejected_steps > 0u, "long trial was rejected"); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_UNRESOLVED, "rejected trial does not create DARK"); CHECK(out.reason == RAY_REASON_BUDGET_EXHAUSTED, "unresolved reason is budget exhaustion"); } /* ------------------------------------------------------------------ */ /* 10. RK4 real RHS accounting and explicit unknown-stepper rejection */ /* ------------------------------------------------------------------ */ static void test_rk4_rhs_cost(void) { /* Four fixed RK4 steps over a finite interval cost exactly four RHS each. */ { SpacetimeSource source = {0}; CHECK(spacetime_create_minkowski(&source, 1.0e9) == 0, "create minkowski"); GeodesicRayState state; fixture_state(0.0, &state); GeodesicTraceConfig config; memset(&config, 0, sizeof config); config.stepper = GEODESIC_STEPPER_RK4; config.coordinate_time_step = 0.25; config.max_steps = 100; config.threshold.kind = THRESHOLD_DISABLED; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "rk4 slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_ACTIVE, "rk4 finite interval is active"); CHECK(state.steps == 4u, "rk4 took four accepted steps"); CHECK(state.rhs_evaluations == 16u, "rk4 four steps cost sixteen RHS"); CHECK(state.rejected_steps == 0u, "rk4 never rejects"); spacetime_destroy(&source); } /* A stage failure costs the real number of evaluations performed. */ { FixtureContext context; memset(&context, 0, sizeof context); context.bad_status = SPACETIME_POINT_OUT_OF_DOMAIN; context.bad_lo = -1.0e30; context.bad_hi = -0.01; /* stage 1 at t = -0.05 fails, the initial is OK */ context.radius = 1.0e9; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(0.0, &state); GeodesicTraceConfig config; memset(&config, 0, sizeof config); config.stepper = GEODESIC_STEPPER_RK4; config.coordinate_time_step = 0.1; config.max_steps = 10; config.threshold.kind = THRESHOLD_DISABLED; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "rk4 failure slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.outcome == RAY_OUTCOME_INCOMPLETE && out.reason == RAY_REASON_OUT_OF_DOMAIN, "rk4 failed stage reports the point status"); CHECK(state.rhs_evaluations == 2u, "rk4 failed stage-1 cost is the real two RHS"); CHECK(state.rejected_steps == 0u, "rk4 failure is not a rejection"); } /* Escape localization adds its own real RHS cost on top of the steps. */ { FixtureContext context; memset(&context, 0, sizeof context); context.radius = 5.0; SpacetimeSource source = {.ops = &fixture_ops, .context = &context}; GeodesicRayState state; fixture_state(2.0, &state); GeodesicTraceConfig config; memset(&config, 0, sizeof config); config.stepper = GEODESIC_STEPPER_RK4; config.coordinate_time_step = 0.5; config.max_steps = 100; config.threshold.kind = THRESHOLD_DISABLED; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -20.0, &slab) == 0, "rk4 escape slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -20.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "rk4 escape still terminates"); CHECK(out.accepted_steps > 0u && out.rhs_evaluations > 4ul * out.accepted_steps, "rk4 escape localization adds real RHS cost"); } } static void test_unknown_stepper_rejected(void) { SpacetimeSource source = {0}; CHECK(spacetime_create_minkowski(&source, 10.0) == 0, "create minkowski"); const ObserverState observer = flat_observer_at((double[]){0.0, 0.0, 0.0}); GeodesicTraceConfig config; memset(&config, 0, sizeof config); config.stepper = (GeodesicStepper)99; config.coordinate_time_step = 0.25; config.max_steps = 10; config.threshold.kind = THRESHOLD_DISABLED; const RayEndpoint traced = geodesic_trace_past( &source, &observer, (double[]){1.0, 0.0, 0.0}, &config); CHECK(traced.outcome == RAY_OUTCOME_INCOMPLETE && traced.reason == RAY_REASON_UNKNOWN_STEPPER, "trace_past rejects an unknown stepper"); GeodesicRayState state; fixture_state(0.0, &state); const RayEndpoint resumed = geodesic_trace_past_from_state(&source, &state, &config); CHECK(resumed.outcome == RAY_OUTCOME_INCOMPLETE && resumed.reason == RAY_REASON_UNKNOWN_STEPPER, "trace_past_from_state rejects an unknown stepper"); MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "unknown slab"); RayEndpoint advanced = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &advanced); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && advanced.outcome == RAY_OUTCOME_INCOMPLETE && advanced.reason == RAY_REASON_UNKNOWN_STEPPER, "advance rejects an unknown stepper"); spacetime_destroy(&source); } /* ------------------------------------------------------------------ */ /* 11. Production advance on a curved synthetic metric: real single- */ /* step inside->outside->inside, tangent, and endpoint roots. */ /* */ /* gamma = I, alpha = 1, K = 0, beta(t) = (3 - 8 (t0 - t), 0, 0) with */ /* no spatial beta derivative. This is a flat metric in accelerating */ /* translated coordinates: Pi and L are constant, the past spatial path */ /* is x(s) = 4 s (1 - s) for s = t0 - t in [0, 1]. It is a geometric */ /* numerical fixture for the production event layer, not a physical */ /* outer-region model; the terminal Minkowski conversion only reads the */ /* legal local metric at the crossing. */ /* ------------------------------------------------------------------ */ typedef struct { double t0; double radius; } CurvedContext; static SpacetimePointStatus curved_eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric) { const CurvedContext *c = source->context; (void)x; *metric = flat_metric(); metric->beta[0] = 3.0 - 8.0 * (c->t0 - t); return SPACETIME_POINT_OK; } static size_t curved_end_count(const SpacetimeSource *source) { (void)source; return 1; } static int curved_end(const SpacetimeSource *source, size_t index, SpacetimeAsymptoticEnd *out) { (void)source; if (index != 0) return -1; *out = (SpacetimeAsymptoticEnd){ .end_id = 0, .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, .mass = 0.0, .frame_origin = {0.0, 0.0, 0.0}, .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; return 0; } static int curved_worldtube(const SpacetimeSource *source, SpacetimeEndId end_id, double t, SpacetimeEscapeWorldtubeSample *out) { const CurvedContext *c = source->context; (void)t; if (end_id != 0) return -1; *out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0}, .velocity = {0.0, 0.0, 0.0}, .radius = c->radius, .radius_rate = 0.0, .velocity_constant = 1, .valid = 1}; return 0; } static const SpacetimeOps curved_ops = { .eval = curved_eval, .classify = event_classify, .asymptotic_end_count = curved_end_count, .asymptotic_end = curved_end, .escape_worldtube_sample = curved_worldtube, .destroy = event_destroy, }; /* Null past state: Pi = (-1, 0, 0) has gamma^{ij} Pi_i Pi_j = 1. */ static void curved_state(GeodesicRayState *state) { memset(state, 0, sizeof *state); state->Pi[0] = -1.0; } /* The earliest exit of the inside->outside->inside path is s = 0.1464466. */ static double curved_first_exit_s(void) { return (1.0 - sqrt(0.5)) / 2.0; } static void test_event_curved_in_out_in(void) { CurvedContext c = {.t0 = 0.0, .radius = 0.5}; SpacetimeSource source = {.ops = &curved_ops, .context = &c}; GeodesicRayState state; curved_state(&state); const double expected_s = curved_first_exit_s(); GeodesicTraceConfig config = dp_config(1e-10, 1.0, 1.0e6); config.max_step = 1.0; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "curved slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "curved in-out-in escapes at the first exit"); CHECK(out.end_id == 0, "curved escape end id"); CHECK(fabs(out.final_x[0] - 0.5) < 1e-6, "curved earliest exit position"); CHECK(fabs(out.stop_coordinate_time + expected_s) < 1e-6, "curved earliest exit time"); MetricData metric; CHECK(spacetime_eval(&source, 0.0, state.x, &metric) == SPACETIME_POINT_OK && fabs(null_residual(&metric, state.Pi) - 1.0) < 1e-12, "curved fixture state is null"); } static void test_event_curved_tangent(void) { /* The path peaks at x = 1, exactly touching the R = 1 sphere at s = 0.5. * The tangent touch must not terminate the ray. */ CurvedContext c = {.t0 = 0.0, .radius = 1.0}; SpacetimeSource source = {.ops = &curved_ops, .context = &c}; GeodesicRayState state; curved_state(&state); GeodesicTraceConfig config = dp_config(1e-10, 1.0, 1.0e6); config.max_step = 1.0; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "tangent slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_ACTIVE, "tangent touch does not escape"); CHECK(state.coordinate_time == -1.0, "tangent trace reached the slab boundary"); } static void test_event_curved_endpoint_root(void) { /* Truncate the step exactly at the first exit so F == 0 lands on the step * endpoint. The in-step enumeration must not accept it; the following step * starts at the boundary root and escapes reliably. */ const double expected_s = curved_first_exit_s(); CurvedContext c = {.t0 = 0.0, .radius = 0.5}; SpacetimeSource source = {.ops = &curved_ops, .context = &c}; GeodesicRayState state; curved_state(&state); GeodesicTraceConfig config = dp_config(1e-12, expected_s, 1.0e6); config.max_step = expected_s; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "endpoint slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "endpoint root escapes on a later step"); CHECK(out.accepted_steps == 2u, "endpoint root is detected on the following boundary step"); CHECK(fabs(out.stop_coordinate_time + expected_s) < 1e-6, "endpoint root crossing time"); CHECK(fabs(out.final_x[0] - 0.5) < 1e-6, "endpoint root position"); } /* LOG_P0 is not polynomial in the dense output, so under DP54 it is still only * checked on trusted accepted states; it is not the CLI default and needs no * dense candidate. */ static void test_event_log_p0_accepted_state(void) { AlphaContext c; memset(&c, 0, sizeof c); c.k = 1.0; c.radius = 10.0; SpacetimeSource source = {.ops = &alpha_ops, .context = &c}; GeodesicRayState state; alpha_state(0.2, c.k, &state); GeodesicTraceConfig config = dp_config(1e-6, 0.25, 1.0e6); config.max_step = 0.25; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_P0, .value = 0.2, .policy_version = 1}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "log_p0 slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_DARK, "LOG_P0 truncates on an accepted state"); CHECK(out.reason == RAY_REASON_REDSHIFT_LIMIT, "LOG_P0 dark reason"); CHECK(out.threshold_value >= 0.2, "LOG_P0 records the reached monitored value"); } /* ------------------------------------------------------------------ */ /* 12. Review-fix regressions: representable step times, absorbed */ /* endpoint roots, and actual-time escape/threshold ordering. */ /* ------------------------------------------------------------------ */ /* A nominal step must commit the representable interval it actually * integrated: x must equal the elapsed coordinate time to local ULP even when * t + h rounds badly. */ static void test_event_step_time_precision(void) { const double origins[3] = {0.0, 1.0e12, 1.0e14}; for (int i = 0; i < 3; ++i) { FixtureContext c; memset(&c, 0, sizeof c); c.radius = 1.0e6; SpacetimeSource source = {.ops = &fixture_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.coordinate_time = origins[i]; state.integration_start_time = origins[i]; state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-9, 0.1, 10.0); config.max_step = 0.1; config.max_steps = 10; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, origins[i], origins[i] - 2.0, &slab) == 0, "time-precision slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray( slab, &state, origins[i] - 2.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_UNRESOLVED, "ten-step budget is unresolved"); CHECK(state.steps == 10u, "ten accepted steps"); const double elapsed = origins[i] - state.coordinate_time; CHECK(fabs(state.x[0] - elapsed) < 1e-10, "committed position matches elapsed coordinate time"); MetricData metric; CHECK(spacetime_eval(&source, state.coordinate_time, state.x, &metric) == SPACETIME_POINT_OK && fabs(null_residual(&metric, state.Pi) - 1.0) < 1e-12, "time-precision fixture stays null"); } /* At a huge coordinate-time origin a nominal step below the local ULP * cannot define a representable nonzero interval: report the exact * TIME_STEP_UNREPRESENTABLE cause instead of committing a zero-length step. */ { const double origin = 1.0e17; FixtureContext c; memset(&c, 0, sizeof c); c.radius = 1.0e30; SpacetimeSource source = {.ops = &fixture_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.coordinate_time = origin; state.integration_start_time = origin; state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-9, 1e-12, 1.0e30); config.min_step = 1e-12; config.max_step = 1.0; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, origin, origin - 1.0e30, &slab) == 0, "unrepresentable slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, origin - 1.0e30, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_TIME_STEP_UNREPRESENTABLE, "sub-ULP step is unrepresentable"); } } /* A crossing whose dense root is absorbed at theta == 1 (F_after strictly * positive) must still escape via the actual-trajectory localizer. */ static void test_event_absorbed_endpoint_exit(void) { double radius = 1.0; for (int j = 0; j < 3; ++j) radius = nextafter(radius, 0.0); FixtureContext c; memset(&c, 0, sizeof c); c.radius = radius; SpacetimeSource source = {.ops = &fixture_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-9, 1.0, 10.0); config.max_step = 1.0; config.max_steps = 4; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -5.0, &slab) == 0, "absorbed-root slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -5.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "absorbed endpoint root still escapes"); CHECK(fabs(out.final_x[0] - radius) <= 4.0 * DBL_EPSILON, "absorbed-root crossing stays on the surface"); } /* A start within the local geometric ULP outside the boundary is a boundary * state and must escape when it moves outside. */ static void test_event_boundary_roundoff_exit(void) { FixtureContext c; memset(&c, 0, sizeof c); c.radius = 1.0; SpacetimeSource source = {.ops = &fixture_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.x[0] = nextafter(1.0, 2.0); state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-9, 0.5, 10.0); config.max_step = 0.5; config.max_steps = 4; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -5.0, &slab) == 0, "roundoff slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -5.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "boundary roundoff start escapes outward"); } /* A single-end state clearly outside its worldtube is reported with the exact * OUTSIDE_WORLDTUBE cause, not an unbounded search that can only end as budget * exhaustion, and the last trusted accepted state plus the real cost counters * are preserved (the failed event step is never committed). */ static void test_event_clearly_outside_protocol(void) { FixtureContext c; memset(&c, 0, sizeof c); c.radius = 1.0; SpacetimeSource source = {.ops = &fixture_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.x[0] = 2.0; state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-9, 0.5, 10.0); config.max_step = 0.5; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -5.0, &slab) == 0, "outside slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -5.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.outcome == RAY_OUTCOME_INCOMPLETE && out.reason == RAY_REASON_OUTSIDE_WORLDTUBE, "clearly-outside single end is an explicit failure"); CHECK(out.stop_coordinate_time == 0.0 && out.final_x[0] == 2.0 && out.final_x[1] == 0.0 && out.final_x[2] == 0.0 && out.accepted_steps == 0u, "outside-worldtube failure preserves the last trusted state"); CHECK(state.coordinate_time == 0.0 && state.x[0] == 2.0 && state.steps == 0u && state.rhs_evaluations > 0ul, "outside-worldtube failure keeps the trusted state and real cost"); } /* Escape and threshold are both localized on the actual trajectory and ordered * by real coordinate time; the ODE tolerance is not a time tie. */ static void test_event_threshold_actual_order(void) { const double exact = log(0.8 / 0.5); const double tols[2] = {1e-4, 1e-9}; for (int t = 0; t < 2; ++t) { /* Below the escape energy -> DARK first, even at loose tolerance. */ { AlphaContext c; memset(&c, 0, sizeof c); c.k = 1.0; c.radius = 0.5; SpacetimeSource source = {.ops = &alpha_ops, .context = &c}; GeodesicRayState state; alpha_state(0.2, c.k, &state); GeodesicTraceConfig config = dp_config(tols[t], 1.0, 10.0); config.max_step = 1.0; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = exact - 2.0e-5, .policy_version = 3}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -3.0, &slab) == 0, "order dark slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -3.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_DARK && out.reason == RAY_REASON_REDSHIFT_LIMIT, "sub-escape threshold is DARK"); CHECK(out.threshold_value >= exact - 2.0e-5, "DARK records a reached threshold"); CHECK(out.stop_coordinate_time > -exact, "DARK stop precedes the escape crossing"); } /* Above the escape energy -> ESCAPED. */ { AlphaContext c; memset(&c, 0, sizeof c); c.k = 1.0; c.radius = 0.5; SpacetimeSource source = {.ops = &alpha_ops, .context = &c}; GeodesicRayState state; alpha_state(0.2, c.k, &state); GeodesicTraceConfig config = dp_config(tols[t], 1.0, 10.0); config.max_step = 1.0; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = exact + 1.0e-5, .policy_version = 3}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -3.0, &slab) == 0, "order escape slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -3.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "super-escape threshold escapes"); } } } /* ------------------------------------------------------------------ */ /* 13. Non-monotonic L pulse: a mid-trace dark truncation must fire */ /* even when the accepted endpoint returns below the threshold. */ /* */ /* Time-dependent gamma metric: alpha = 1, beta = 0, gamma_ij = */ /* a(t)^2 delta_ij, K_ij = -a(t)^2 H(t) delta_ij with a = exp(-F) and */ /* H = d ln a/dt = F'(s), s = t0 - t. Spatial derivatives are zero. */ /* With Pi_i = a n_i the past path has dL/ds = F'(s), so L(s) = F(s): */ /* the metric's gamma time dependence and K are consistent (a partial */ /* alpha(t) with K = 0 would not produce this evolution). */ /* ------------------------------------------------------------------ */ typedef struct { double A; double t0; double radius; int kind; /* 0 = two-pulse 16 A s(1-s)(s-.5)^2, 1 = single 4 A s(1-s) */ size_t end_count; } CosmoContext; static double cosmo_amp(double s, double A, int kind) { if (kind == 0) return 16.0 * A * s * (1.0 - s) * (s - 0.5) * (s - 0.5); if (kind == 2) return A * sin(6.283185307179586476925286766559 * s); return 4.0 * A * s * (1.0 - s); } static double cosmo_amp_prime(double s, double A, int kind) { if (kind == 0) return 16.0 * A * (-4.0 * s * s * s + 6.0 * s * s - 2.5 * s + 0.25); if (kind == 2) return A * 6.283185307179586476925286766559 * cos(6.283185307179586476925286766559 * s); return 4.0 * A * (1.0 - 2.0 * s); } static SpacetimePointStatus cosmo_eval(const SpacetimeSource *source, double t, const double x[3], MetricData *metric) { const CosmoContext *c = source->context; (void)x; const double s = c->t0 - t; const double F = cosmo_amp(s, c->A, c->kind); const double H = cosmo_amp_prime(s, c->A, c->kind); const double a = exp(-F); const double g = a * a; *metric = flat_metric(); metric->alpha = 1.0; metric->gamma[0][0] = metric->gamma[1][1] = metric->gamma[2][2] = g; metric->K[0][0] = metric->K[1][1] = metric->K[2][2] = -g * H; return SPACETIME_POINT_OK; } static size_t cosmo_end_count(const SpacetimeSource *source) { return ((const CosmoContext *)source->context)->end_count; } static int cosmo_end(const SpacetimeSource *source, size_t index, SpacetimeAsymptoticEnd *out) { const CosmoContext *c = source->context; if (index >= c->end_count) return -1; *out = (SpacetimeAsymptoticEnd){ .end_id = 0, .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI, .mass = 0.0, .frame_origin = {0.0, 0.0, 0.0}, .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; return 0; } static int cosmo_worldtube(const SpacetimeSource *source, SpacetimeEndId end_id, double t, SpacetimeEscapeWorldtubeSample *out) { const CosmoContext *c = source->context; (void)t; if (end_id != 0) return -1; *out = (SpacetimeEscapeWorldtubeSample){.center = {0.0, 0.0, 0.0}, .velocity = {0.0, 0.0, 0.0}, .radius = c->radius, .radius_rate = 0.0, .velocity_constant = 1, .valid = 1}; return 0; } static const SpacetimeOps cosmo_ops = { .eval = cosmo_eval, .classify = event_classify, .asymptotic_end_count = cosmo_end_count, .asymptotic_end = cosmo_end, .escape_worldtube_sample = cosmo_worldtube, .destroy = event_destroy, }; /* Analytic first upcrossing of F = threshold for each pulse shape. */ static double cosmo_first_crossing(double A, double threshold, int kind) { if (kind == 0) { const double disc = 0.0625 - threshold / (4.0 * A); return 0.5 - sqrt((0.25 + sqrt(disc)) / 2.0); } return (1.0 - sqrt(1.0 - threshold / A)) / 2.0; } static void cosmo_run(double A, double threshold, int kind, size_t end_count, double tol, int expect_dark) { CosmoContext c; memset(&c, 0, sizeof c); c.A = A; c.t0 = 0.0; c.radius = 100.0; c.kind = kind; c.end_count = end_count; SpacetimeSource source = {.ops = &cosmo_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.Pi[0] = -1.0; /* null: gamma^{ij} Pi_i Pi_j = |Pi|^2/a^2 = 1 */ state.log_alpha_p0 = 0.0; state.log_alpha_p0_0 = 0.0; GeodesicTraceConfig config = dp_config(tol, 1.0, 10.0); config.max_step = 1.0; config.max_steps = 100; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = threshold, .policy_version = 1}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "cosmo slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED, "cosmo terminates"); if (!expect_dark) { CHECK(out.outcome != RAY_OUTCOME_DARK, "cosmo expected non-dark"); return; } CHECK(out.outcome == RAY_OUTCOME_DARK && out.reason == RAY_REASON_REDSHIFT_LIMIT, "mid-pulse dark truncation fires"); CHECK(out.threshold_value >= threshold, "dark records a reached threshold"); const double s_ref = cosmo_first_crossing(A, threshold, kind); CHECK(fabs(-out.stop_coordinate_time - s_ref) < 1e-5, "dark stop is the first upcrossing"); CHECK(-out.stop_coordinate_time < 0.1464466, "dark stop is before the first peak"); } static void test_event_nonmonotonic_threshold(void) { const double A = 0.08; const double T = 0.01; /* Both accepted-step tolerances must find the same first crossing even * though the accepted endpoint returns to L = 0 < T. */ cosmo_run(A, T, 0, 1, 1e-3, 1); cosmo_run(A, T, 0, 1, 1e-9, 1); /* The same energy policy applies to a backend that declares no ends. */ cosmo_run(A, T, 0, 0, 1e-3, 1); /* A single non-monotonic pulse (one upcrossing) behaves the same way. */ cosmo_run(A, T, 1, 1, 1e-3, 1); cosmo_run(A, T, 1, 1, 1e-9, 1); } /* Large-origin moving worldtube: the production dense polynomial must match * the actual F at the same theta, so the event timing cannot drift with the * rounding of the sampled half-step. */ static void test_event_worldtube_polynomial_time_origin(void) { const double t0 = 1.0e14; EventContext c; memset(&c, 0, sizeof c); c.end_count = 1; c.id[0] = 0; c.center[0][0] = 0.0; c.velocity[0][0] = -0.99; c.t_ref[0] = t0; c.radius_a[0] = 1.0; c.velocity_constant[0] = 1; SpacetimeSource source = {.ops = &event_ops, .context = &c}; GeodesicTraceConfig config = dp_config(1e-9, 0.05, 1.0e6); config.max_step = 0.05; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, t0, t0 - 1.0, &slab) == 0, "probe slab"); /* Only theta values whose target t_before + h*theta is exactly representable * compare cleanly; 1/3 lands one ULP inside the step for h = 3 ULP. */ const double thetas[3] = {0.0, 1.0 / 3.0, 1.0}; for (int k = 0; k < 3; ++k) { GeodesicRayState state; memset(&state, 0, sizeof state); state.coordinate_time = t0; state.integration_start_time = t0; state.x[0] = 0.5; state.Pi[0] = -1.0; double dense_value = NAN, actual_value = NAN; CHECK(geodesic_event_worldtube_dense_probe(slab, &config, 0, &state, t0 - 1.0, thetas[k], &dense_value, &actual_value) == 1, "worldtube dense probe runs"); const double scale = fmax(1.0, fabs(actual_value)); CHECK(fabs(dense_value - actual_value) < 1e-10 * scale, "dense worldtube F matches actual F at a large time origin"); } spacetime_free_slab(slab); } /* Explicitly prove the loose-tolerance pulse step is accepted as one step * with both endpoints below the threshold (so the test covers the miss). */ static void test_event_nonmonotonic_accepts_one_step(void) { const double A = 0.08; const double T = 0.01; CosmoContext c; memset(&c, 0, sizeof c); c.A = A; c.radius = 100.0; c.kind = 0; c.end_count = 1; SpacetimeSource source = {.ops = &cosmo_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-3, 1.0, 10.0); config.max_step = 1.0; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = T, .policy_version = 1}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "cosmo one-step slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_DARK, "loose-tolerance pulse is dark"); CHECK(out.accepted_steps == 1u, "loose tolerance accepts the whole pulse step in one step"); /* Both endpoints of the accepted step are genuinely below the threshold, so * the dark stop can only come from the mid-step pulse. */ CHECK(cosmo_amp(0.0, A, 0) < T && cosmo_amp(1.0, A, 0) < T, "pulse endpoints are below the threshold"); CHECK(out.threshold_value >= T, "pulse dark reached the threshold"); } /* The threshold is a closed (>=) event: a crossing exactly at the step/slab * endpoint must be accepted through the closed-end candidate, not sent into a * permanent shrink. */ static void test_event_threshold_closed_endpoint(void) { AlphaContext c; memset(&c, 0, sizeof c); c.k = 1.0; c.radius = 10.0; SpacetimeSource source = {.ops = &alpha_ops, .context = &c}; /* Phase 1: with the threshold disabled, find the exactly accepted endpoint * monitored value at s = 1 (L - L0 = s along this ray). */ GeodesicRayState probe; alpha_state(0.2, c.k, &probe); GeodesicTraceConfig pconfig = dp_config(1e-4, 1.0, 10.0); pconfig.max_step = 1.0; pconfig.threshold.kind = THRESHOLD_DISABLED; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "closed endpoint slab"); RayEndpoint out = blank_endpoint(); (void)geodesic_advance_past_ray(slab, &probe, -1.0, &pconfig, &out); const double endpoint_value = probe.log_alpha_p0 - probe.log_alpha_p0_0; CHECK(fabs(probe.coordinate_time + 1.0) < 1e-9, "endpoint probe reached the slab boundary"); CHECK(fabs(endpoint_value - 1.0) < 1e-9, "endpoint value is near one"); /* Phase 2: threshold set to that endpoint value: a closed (>=) event at the * step/slab endpoint must be a DARK event, not a permanent shrink. */ GeodesicRayState state; alpha_state(0.2, c.k, &state); GeodesicTraceConfig config = dp_config(1e-4, 1.0, 10.0); config.max_step = 1.0; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = endpoint_value, .policy_version = 1}; (void)geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(out.outcome == RAY_OUTCOME_DARK && out.reason == RAY_REASON_REDSHIFT_LIMIT, "endpoint threshold is a dark event"); CHECK(fabs(out.stop_coordinate_time + 1.0) < 1e-6, "endpoint threshold is localized at the endpoint"); CHECK(out.threshold_value >= endpoint_value, "endpoint threshold reached"); } /* When an event candidate cannot be confirmed on the actual trajectory and * the bounded shrink is exhausted, this is an INCOMPLETE integration failure, * never a rewritten DARK with a restored (steps = 0) state. */ static void test_event_threshold_confirmation_exhausted(void) { CosmoContext c; memset(&c, 0, sizeof c); c.A = 0.1; c.radius = 100.0; c.kind = 2; /* non-polynomial L pulse; the coarse dense quartic cannot be confirmed at this step/tolerance */ c.end_count = 1; SpacetimeSource source = {.ops = &cosmo_ops, .context = &c}; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -1.0, &slab) == 0, "confirmation slab"); /* Phase 1: with the threshold disabled the same coarse step is accepted, so * the phase-2 failure is an event-confirmation failure, not a step rejection. */ GeodesicRayState probe; memset(&probe, 0, sizeof probe); probe.Pi[0] = -1.0; GeodesicTraceConfig pconfig = dp_config(1.0, 1.0, 10.0); pconfig.max_step = 1.0; pconfig.max_steps = 100; pconfig.consecutive_rejection_limit = 1; pconfig.atol_L = 1.0e6; pconfig.rtol = 1.0e6; pconfig.threshold.kind = THRESHOLD_DISABLED; RayEndpoint probe_out = blank_endpoint(); (void)geodesic_advance_past_ray(slab, &probe, -1.0, &pconfig, &probe_out); CHECK(probe.steps >= 1u, "the coarse step itself is accepted"); GeodesicRayState state; memset(&state, 0, sizeof state); state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1.0, 1.0, 10.0); config.max_step = 1.0; config.max_steps = 100; config.consecutive_rejection_limit = 1; config.atol_L = 1.0e6; config.rtol = 1.0e6; config.threshold = (ThresholdPolicy){.kind = THRESHOLD_LOG_ENERGY_GROWTH, .value = 0.13, .policy_version = 1}; RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -1.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.outcome == RAY_OUTCOME_INCOMPLETE && out.reason == RAY_REASON_THRESHOLD_EVENT_UNCONFIRMED, "unconfirmable threshold is an incomplete integration failure"); CHECK(out.outcome != RAY_OUTCOME_DARK, "unconfirmable threshold is not rewritten into DARK"); CHECK(state.steps == 0u, "last trusted state carries no fabricated accepted step"); CHECK(state.rejected_steps >= 1u, "rejection cost is real"); CHECK(out.rhs_evaluations > 0ul, "actual RHS cost is recorded"); } /* A pure error-estimate rejection that exhausts the consecutive-rejection quota * reports the exact REJECTION_LIMIT cause. With a one-rejection quota the * limit is reached before any minimum-step check, so this is not MIN_STEP. */ static void test_rejection_limit_reason(void) { AlphaContext c; memset(&c, 0, sizeof c); c.k = 1.0; c.radius = 10.0; SpacetimeSource source = {.ops = &alpha_ops, .context = &c}; GeodesicRayState state; alpha_state(0.2, c.k, &state); GeodesicTraceConfig config = dp_config(1e-15, 1.0, 10.0); config.max_step = 1.0; config.min_step = 1.0; config.consecutive_rejection_limit = 1; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -5.0, &slab) == 0, "reject-limit slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -5.0, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_FAILED && out.reason == RAY_REASON_REJECTION_LIMIT, "error-estimate rejection quota reports REJECTION_LIMIT"); CHECK(state.coordinate_time == 0.0 && state.rejected_steps >= 1u, "rejected trial did not advance the accepted state"); } /* Uniform unit scaling: with every length/time quantity scaled by 1e-9 the * same flat escape must give event time/scale ~ 3 and x/scale ~ 5, i.e. no * absolute coordinate-time floor breaks small units. */ static void test_event_unit_scaling(void) { const double scale = 1.0e-9; FixtureContext c; memset(&c, 0, sizeof c); c.radius = 5.0 * scale; SpacetimeSource source = {.ops = &fixture_ops, .context = &c}; GeodesicRayState state; memset(&state, 0, sizeof state); state.x[0] = 2.0 * scale; state.Pi[0] = -1.0; GeodesicTraceConfig config = dp_config(1e-10, 0.5, 100.0); config.coordinate_time_step *= scale; config.min_step *= scale; config.max_step *= scale; config.atol_x *= scale; config.max_lookback_time *= scale; MetricSlab *slab = NULL; CHECK(spacetime_load_slab(&source, 0.0, -20.0 * scale, &slab) == 0, "scaled slab"); RayEndpoint out = blank_endpoint(); const GeodesicAdvanceResult r = geodesic_advance_past_ray(slab, &state, -20.0 * scale, &config, &out); spacetime_free_slab(slab); CHECK(r == GEODESIC_ADVANCE_TERMINATED && out.outcome == RAY_OUTCOME_ESCAPED, "scaled flat escape"); CHECK(fabs((0.0 - out.stop_coordinate_time) / scale - 3.0) < 1e-6, "scaled event time / scale is 3"); CHECK(fabs(out.final_x[0] / scale - 5.0) < 1e-6, "scaled crossing position / scale is 5"); } int main(void) { test_minkowski_finite_interval(); test_schwarzschild_radial(); test_over_large_initial_step(); test_fixture_stage_retry(); test_fixture_fatal_and_bounds(); test_budget_lookback_and_resume(); test_flat_escape_localization(); test_event_time_translation(); test_event_subintegration_time_hole(); test_schwarzschild_dp_escape_matches_finish(); test_rk4_compatibility(); test_ray_pool_batch_matches_single(); test_ray_pool_continuation_state(); test_ray_pool_cross_slab_reject(); test_event_polynomial_roots(); test_event_two_ends_order(); test_event_segmented_worldtube(); test_event_nonconstant_unsupported(); test_event_threshold_order(); test_event_rejected_stage_threshold(); test_rk4_rhs_cost(); test_unknown_stepper_rejected(); test_event_curved_in_out_in(); test_event_curved_tangent(); test_event_curved_endpoint_root(); test_event_log_p0_accepted_state(); test_event_step_time_precision(); test_event_absorbed_endpoint_exit(); test_event_boundary_roundoff_exit(); test_event_clearly_outside_protocol(); test_event_threshold_actual_order(); test_event_nonmonotonic_threshold(); test_event_nonmonotonic_accepts_one_step(); test_event_worldtube_polynomial_time_origin(); test_event_threshold_closed_endpoint(); test_event_threshold_confirmation_exhausted(); test_rejection_limit_reason(); test_event_unit_scaling(); if (failures != 0) { fprintf(stderr, "geodesic adaptive regression: %d failure(s)\n", failures); return 1; } puts("geodesic adaptive regression passed"); return 0; }