wyj
0a46a7095b
Add error-controlled DP5(4) integration and trusted first-crossing localization, including non-monotonic energy thresholds and representable-time stepping. Preserve adaptive state and independent step/time retry grants across RayPool, refinement and movie scheduling. Expose numerical controls, record actual persistent-sample costs, and add v3 lens-map provenance with legacy v2 RK4 import. Use DP54 by default and select an 8M Schwarzschild maximum step from bounded scans and a two-run 4K comparison. Retain the conservative minimum-step guard and document critical-ray and backend capability limits. Archive self-contained benchmark inputs and raw output; keep fixed RK4 HDR references explicit. Validation: make -B -j4 BUILD_TYPE=Debug test passed; explicit RK4 HDR references have zero differences. Bounded convergence checks, benchmark reproduction, Release build and focused reviews passed. No numerical-relativity backend is added.
342 lines
12 KiB
C
342 lines
12 KiB
C
/*
|
|
* Standalone critical-classification reference check.
|
|
*
|
|
* Three static Schwarzschild cameras (M=1, escape sphere 256):
|
|
* r30 look 180 (inward), r100 look 180 (inward), r2.1 look 0 (outward),
|
|
* each aimed at theta_crit exactly and at theta_crit +- 1e-7, where
|
|
* theta_crit = asin(3 sqrt(3) sqrt(1 - 2/r) / r)
|
|
* is the local tetrad angle from the camera forward direction.
|
|
*
|
|
* For every direction four DP54 configs are run (initial step 0.1,
|
|
* min_step 1e-14, max_steps 65536, max_lookback_time 6553.6, reject limit 32,
|
|
* camera-relative dark threshold 8):
|
|
*
|
|
* ref_tol1e-12_h0.25 tol 1e-12, max_step 0.25
|
|
* ref_tol1e-13_h0.125 tol 1e-13, max_step 0.125
|
|
* tol1e-11_h2 tol 1e-11, max_step 2
|
|
* tol1e-11_h8 tol 1e-11, max_step 8
|
|
*
|
|
* Total 3 cameras * 3 directions * 4 configs = 36 endpoint traces.
|
|
*
|
|
* Raw per-run rows go to raw/critical_ref.csv; stdout prints, for all nine
|
|
* directions, ref_vs_tighter / h2_vs_h8 / ref_vs_h2 / ref_vs_h8 comparisons
|
|
* (outcome/reason/end, sky angle, relative g, dark stop-time and threshold
|
|
* shifts) and a NEIGHBOR_REFERENCE_CHECK that verifies the two reference
|
|
* tiers agree at the +-1e-7 neighbours. The exactly-critical direction is a
|
|
* separatrix: a class flip or a large direction/stop-time difference there is
|
|
* reported as EXACT_CRITICAL_REFERENCE_NOT_CONVERGED and is diagnostic only.
|
|
*
|
|
* Exit status: 0 when the neighbours' reference tiers agree and no neighbour
|
|
* returned an error/unresolved endpoint or a zero-step trace; 1 when a
|
|
* neighbour reference check fails; 2 on setup error. The exact-critical
|
|
* separatrix is never an assertion.
|
|
*
|
|
* Link line (the Schwarzschild backend is textually included):
|
|
* cc -std=c11 -O2 -Isrc critical_ref_check.c geodesic.c asymptotic.c \
|
|
* asymptotic_schwarzschild.c spacetime_common.c observer.c -lm
|
|
*/
|
|
#define _POSIX_C_SOURCE 200809L
|
|
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
#include "../../src/spacetime_schwarzschild.c"
|
|
|
|
#include "geodesic.h"
|
|
#include "observer.h"
|
|
#include "spacetime.h"
|
|
|
|
#define PI 3.14159265358979323846
|
|
|
|
typedef struct {
|
|
const char *name;
|
|
double tol;
|
|
double max_step;
|
|
} CfgDef;
|
|
|
|
static const CfgDef cfgs[4] = {
|
|
{"ref_tol1e-12_h0.25", 1e-12, 0.25},
|
|
{"ref_tol1e-13_h0.125", 1e-13, 0.125},
|
|
{"tol1e-11_h2", 1e-11, 2.0},
|
|
{"tol1e-11_h8", 1e-11, 8.0},
|
|
};
|
|
|
|
typedef struct {
|
|
const char *name;
|
|
double r;
|
|
double look_ra_deg;
|
|
} CamDef;
|
|
|
|
static const CamDef cams[3] = {
|
|
{"r30_inward", 30.0, 180.0},
|
|
{"r100_inward", 100.0, 180.0},
|
|
{"r2p1_outward", 2.1, 0.0},
|
|
};
|
|
|
|
typedef struct {
|
|
int outcome, reason;
|
|
unsigned int end_id;
|
|
double stop_t, g, thr;
|
|
double L, L0;
|
|
double n[3];
|
|
unsigned int steps, rejected;
|
|
unsigned long rhs;
|
|
} Res;
|
|
|
|
static double b_dot3(const double a[3], const double b[3]) {
|
|
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
|
|
}
|
|
|
|
static double b_ang_delta(const double a[3], const double b[3]) {
|
|
const double na = sqrt(b_dot3(a, a)), nb = sqrt(b_dot3(b, b));
|
|
if (!(na > 0.0) || !(nb > 0.0))
|
|
return NAN;
|
|
const double cross[3] = {a[1] * b[2] - a[2] * b[1],
|
|
a[2] * b[0] - a[0] * b[2],
|
|
a[0] * b[1] - a[1] * b[0]};
|
|
return atan2(sqrt(b_dot3(cross, cross)) / (na * nb),
|
|
b_dot3(a, b) / (na * nb));
|
|
}
|
|
|
|
static double critical_angle(double r) {
|
|
const double b = 3.0 * sqrt(3.0) * sqrt(1.0 - 2.0 / r) / r;
|
|
return (b < 1.0) ? asin(b) : PI / 2.0;
|
|
}
|
|
|
|
static void dir_from_theta(double th, double n[3]) {
|
|
n[0] = cos(th);
|
|
n[1] = sin(th);
|
|
n[2] = 0.0;
|
|
const double nn = sqrt(b_dot3(n, n));
|
|
for (int i = 0; i < 3; ++i)
|
|
n[i] /= nn;
|
|
}
|
|
|
|
static GeodesicTraceConfig make_cfg(const CfgDef *d) {
|
|
GeodesicTraceConfig c;
|
|
memset(&c, 0, sizeof c);
|
|
c.coordinate_time_step = 0.1;
|
|
c.max_steps = 65536;
|
|
c.threshold.kind = THRESHOLD_LOG_ENERGY_GROWTH;
|
|
c.threshold.value = 8.0;
|
|
c.threshold.policy_version = 3;
|
|
c.stepper = GEODESIC_STEPPER_DP54;
|
|
c.atol_x = c.atol_Pi = c.atol_L = c.rtol = d->tol;
|
|
c.min_step = 1e-14;
|
|
c.max_step = d->max_step;
|
|
c.consecutive_rejection_limit = 32;
|
|
c.max_lookback_time = 6553.6;
|
|
return c;
|
|
}
|
|
|
|
static int build_observer(const CamDef *cam, const SpacetimeSource *s,
|
|
ObserverState *o) {
|
|
const double pos[3] = {cam->r, 0.0, 0.0};
|
|
MetricData m;
|
|
if (spacetime_eval(s, 0.0, pos, &m) != SPACETIME_POINT_OK)
|
|
return -1;
|
|
ObserverCamera camera;
|
|
memset(&camera, 0, sizeof camera);
|
|
camera.position[0] = pos[0];
|
|
camera.position[1] = pos[1];
|
|
camera.position[2] = pos[2];
|
|
camera.look_ra_deg = cam->look_ra_deg;
|
|
camera.look_dec_deg = 0.0;
|
|
return observer_from_coordinate_camera(&m, &camera, o, NULL) ==
|
|
OBSERVER_BUILD_OK
|
|
? 0
|
|
: -1;
|
|
}
|
|
|
|
static Res run_one(const SpacetimeSource *s, const ObserverState *o,
|
|
const double dir[3], const GeodesicTraceConfig *cfg) {
|
|
Res r;
|
|
memset(&r, 0, sizeof r);
|
|
const RayEndpoint e = geodesic_trace_past(s, o, dir, cfg);
|
|
r.outcome = e.outcome;
|
|
r.reason = e.reason;
|
|
r.end_id = e.end_id;
|
|
r.stop_t = e.stop_coordinate_time;
|
|
r.g = e.frequency_ratio;
|
|
r.thr = e.threshold_value;
|
|
r.L = e.final_log_alpha_p0;
|
|
r.L0 = e.final_log_alpha_p0_0;
|
|
for (int i = 0; i < 3; ++i)
|
|
r.n[i] = e.n_infinity[i];
|
|
r.steps = e.accepted_steps;
|
|
r.rejected = e.rejected_steps;
|
|
r.rhs = e.rhs_evaluations;
|
|
return r;
|
|
}
|
|
|
|
static const char *outcome_name(int o) {
|
|
switch (o) {
|
|
case RAY_OUTCOME_ESCAPED:
|
|
return "ESC";
|
|
case RAY_OUTCOME_DARK:
|
|
return "DARK";
|
|
case RAY_OUTCOME_UNRESOLVED:
|
|
return "UNRES";
|
|
default:
|
|
return "INC";
|
|
}
|
|
}
|
|
|
|
static void print_pair(const char *cam, int di, double theta, const char *pn,
|
|
const Res *a, const Res *b) {
|
|
const int same = a->outcome == b->outcome && a->reason == b->reason &&
|
|
a->end_id == b->end_id;
|
|
const int esc = a->outcome == RAY_OUTCOME_ESCAPED &&
|
|
b->outcome == RAY_OUTCOME_ESCAPED;
|
|
const double dn = esc ? b_ang_delta(a->n, b->n) : NAN;
|
|
const double dg =
|
|
(esc && a->g > 0.0 && b->g > 0.0) ? fabs(a->g / b->g - 1.0) : NAN;
|
|
const double dstop =
|
|
(isfinite(a->stop_t) && isfinite(b->stop_t)) ? a->stop_t - b->stop_t
|
|
: NAN;
|
|
const double dthr = (isfinite(a->thr) && isfinite(b->thr))
|
|
? a->thr - b->thr
|
|
: NAN;
|
|
printf("CAPAIR camera=%s dir=%d theta=%.17g pair=%s class=%s/%s same=%d "
|
|
"dn_ang=%.17g dgrel=%.17g dstop_t=%.17g dthr=%.17g end=%u/%u\n",
|
|
cam, di, theta, pn, outcome_name(a->outcome), outcome_name(b->outcome),
|
|
same, dn, dg, dstop, dthr, a->end_id, b->end_id);
|
|
}
|
|
|
|
/* A two-reference-tier agreement check for one neighbour direction. Returns
|
|
* 1 when the pair is a valid, agreeing ESC/DARK reference, 0 otherwise. */
|
|
static int neighbour_reference_check(const char *cam, int di, double theta,
|
|
const Res *a, const Res *b) {
|
|
const int class_same = a->outcome == b->outcome && a->reason == b->reason &&
|
|
a->end_id == b->end_id;
|
|
const int clean = (a->outcome == RAY_OUTCOME_ESCAPED ||
|
|
a->outcome == RAY_OUTCOME_DARK) &&
|
|
(b->outcome == RAY_OUTCOME_ESCAPED ||
|
|
b->outcome == RAY_OUTCOME_DARK);
|
|
const int stepped = a->steps > 0 && b->steps > 0;
|
|
const int esc = a->outcome == RAY_OUTCOME_ESCAPED &&
|
|
b->outcome == RAY_OUTCOME_ESCAPED;
|
|
const double dn = esc ? b_ang_delta(a->n, b->n) : NAN;
|
|
const double dg =
|
|
(esc && a->g > 0.0 && b->g > 0.0) ? fabs(a->g / b->g - 1.0) : NAN;
|
|
const double dstop =
|
|
(isfinite(a->stop_t) && isfinite(b->stop_t)) ? a->stop_t - b->stop_t
|
|
: NAN;
|
|
printf("NEIGHBOR_REFERENCE_CHECK camera=%s dir=%d theta=%.17g "
|
|
"ref_class0=%s/%d/%u ref_class1=%s/%d/%u same=%d clean=%d stepped=%d "
|
|
"dn_ang=%.17g dgrel=%.17g dstop_t=%.17g\n",
|
|
cam, di, theta, outcome_name(a->outcome), a->reason, a->end_id,
|
|
outcome_name(b->outcome), b->reason, b->end_id, class_same, clean,
|
|
stepped, dn, dg, dstop);
|
|
return class_same && clean && stepped;
|
|
}
|
|
|
|
int main(void) {
|
|
SpacetimeSource source;
|
|
if (spacetime_create_schwarzschild_ks(&source, 1.0, 256.0)) {
|
|
fprintf(stderr, "FATAL: cannot create Schwarzschild source\n");
|
|
return 2;
|
|
}
|
|
|
|
FILE *f = fopen("raw/critical_ref.csv", "w");
|
|
if (!f) {
|
|
fprintf(stderr, "FATAL: cannot open raw/critical_ref.csv\n");
|
|
return 2;
|
|
}
|
|
fprintf(f, "camera,dir,theta,cfg,outcome,reason,end_id,stop_t,steps,"
|
|
"rejected,rhs,nx,ny,nz,g,thr,L,L0\n");
|
|
|
|
Res res[3][3][4];
|
|
double thetas[3][3];
|
|
for (int ci = 0; ci < 3; ++ci) {
|
|
ObserverState observer;
|
|
if (build_observer(&cams[ci], &source, &observer)) {
|
|
fprintf(stderr, "FATAL: cannot build observer %s\n", cams[ci].name);
|
|
return 2;
|
|
}
|
|
const double tc = critical_angle(cams[ci].r);
|
|
const double th[3] = {tc - 1e-7, tc, tc + 1e-7};
|
|
for (int di = 0; di < 3; ++di) {
|
|
double dir[3];
|
|
dir_from_theta(th[di], dir);
|
|
thetas[ci][di] = th[di];
|
|
for (int fi = 0; fi < 4; ++fi) {
|
|
const GeodesicTraceConfig cfg = make_cfg(&cfgs[fi]);
|
|
const Res r = run_one(&source, &observer, dir, &cfg);
|
|
res[ci][di][fi] = r;
|
|
fprintf(f, "%s,%d,%.17g,%s,%d,%d,%u,%.17g,%u,%u,%lu,%.17g,%.17g,"
|
|
"%.17g,%.17g,%.17g,%.17g,%.17g\n",
|
|
cams[ci].name, di, th[di], cfgs[fi].name, r.outcome, r.reason,
|
|
r.end_id, r.stop_t, r.steps, r.rejected, r.rhs, r.n[0], r.n[1],
|
|
r.n[2], r.g, r.thr, r.L, r.L0);
|
|
printf("RAW camera=%s dir=%d theta=%.17g cfg=%s outcome=%d reason=%d "
|
|
"end=%u stop_t=%.17g steps=%u rejected=%u rhs=%lu g=%.17g "
|
|
"thr=%.17g L=%.17g L0=%.17g nx=%.17g ny=%.17g nz=%.17g\n",
|
|
cams[ci].name, di, th[di], cfgs[fi].name, r.outcome, r.reason,
|
|
r.end_id, r.stop_t, r.steps, r.rejected, r.rhs, r.g, r.thr, r.L,
|
|
r.L0, r.n[0], r.n[1], r.n[2]);
|
|
fflush(stdout);
|
|
}
|
|
}
|
|
}
|
|
fclose(f);
|
|
|
|
/* All nine directions: ref_vs_tighter, h2_vs_h8, ref_vs_h2, ref_vs_h8. */
|
|
for (int ci = 0; ci < 3; ++ci) {
|
|
for (int di = 0; di < 3; ++di) {
|
|
print_pair(cams[ci].name, di, thetas[ci][di], "ref_vs_tighter",
|
|
&res[ci][di][0], &res[ci][di][1]);
|
|
print_pair(cams[ci].name, di, thetas[ci][di], "h2_vs_h8",
|
|
&res[ci][di][2], &res[ci][di][3]);
|
|
print_pair(cams[ci].name, di, thetas[ci][di], "ref_vs_h2",
|
|
&res[ci][di][0], &res[ci][di][2]);
|
|
print_pair(cams[ci].name, di, thetas[ci][di], "ref_vs_h8",
|
|
&res[ci][di][0], &res[ci][di][3]);
|
|
}
|
|
}
|
|
|
|
/* Exact-critical: diagnostic only, never asserted. */
|
|
for (int ci = 0; ci < 3; ++ci) {
|
|
const Res *a = &res[ci][1][0];
|
|
const Res *b = &res[ci][1][1];
|
|
const int same = a->outcome == b->outcome && a->reason == b->reason &&
|
|
a->end_id == b->end_id;
|
|
const int esc = a->outcome == RAY_OUTCOME_ESCAPED &&
|
|
b->outcome == RAY_OUTCOME_ESCAPED;
|
|
const int dark = a->outcome == RAY_OUTCOME_DARK &&
|
|
b->outcome == RAY_OUTCOME_DARK;
|
|
const double dn = esc ? b_ang_delta(a->n, b->n) : NAN;
|
|
const double dstop =
|
|
(isfinite(a->stop_t) && isfinite(b->stop_t)) ? a->stop_t - b->stop_t
|
|
: NAN;
|
|
const int nonconverged =
|
|
!same || (esc && !(dn <= 1e-6)) || (dark && fabs(dstop) > 1e-4);
|
|
printf("EXACT_CRITICAL_DIRECTION camera=%s theta=%.17g ref0=%s/%d/%u "
|
|
"ref1=%s/%d/%u same=%d dn_ang=%.17g dstop_t=%.17g\n",
|
|
cams[ci].name, thetas[ci][1], outcome_name(a->outcome), a->reason,
|
|
a->end_id, outcome_name(b->outcome), b->reason, b->end_id, same, dn,
|
|
dstop);
|
|
if (nonconverged)
|
|
printf("EXACT_CRITICAL_REFERENCE_NOT_CONVERGED camera=%s dn_ang=%.17g "
|
|
"dstop_t=%.17g class_same=%d\n",
|
|
cams[ci].name, dn, dstop, same);
|
|
}
|
|
|
|
/* Neighbour reference check: must pass; this is the assertion-valid path. */
|
|
int fail = 0;
|
|
for (int ci = 0; ci < 3; ++ci) {
|
|
for (int di = 0; di < 3; di += 2) {
|
|
if (!neighbour_reference_check(cams[ci].name, di, thetas[ci][di],
|
|
&res[ci][di][0], &res[ci][di][1]))
|
|
fail = 1;
|
|
}
|
|
}
|
|
|
|
spacetime_destroy(&source);
|
|
printf("CRITICAL_REFERENCE_CHECK status=%d\n", fail ? 1 : 0);
|
|
printf("CRITICAL_REF_DONE\n");
|
|
return fail ? 1 : 0;
|
|
}
|