wyj/GR-raytracing
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
0a46a7095b5950ecac29aab0e3029b088375fc3c
GR-raytracing/benchmarks/adaptive_step_bounds_2026-10-05/a_public_endpoints.c
T
wyj 0a46a7095b Feat: Complete adaptive geodesic tracing with DP54 
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.
2026-10-05 20:27:42 -04:00

868 lines
31 KiB
C
Raw Blame History

/*
* Experiment A: DP54 step-bound scan through the *public* production
* endpoint geodesic_trace_past().
*
* Sub-commands:
* schwarzschild : r30/r100 static inward, r2.1 static outward, r1.5 free-fall
* minkowski : moving observer, escape sphere 64, analytic flat check
* alcubierre : comoving bubble-center camera, vs .3/.9, sigma 1/10/100
*
* Each sub-command computes one DP54 tol=1e-12 reference per ray, then scans
* - upper scan : max_step over a list, min_step fixed
* - min scan : min_step over a list, max_step fixed
* and compares every result against the reference (class, sky angle, g,
* dark threshold margin / stop time, cost). Raw per-ray CSV and stdout
* summaries are produced. No production source is modified.
*
* Link line (do NOT link geodesic.c twice; the backends are textually included
* below, so only the common sources are linked):
* cc -std=c11 -O2 -Isrc a_public_endpoints.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 <time.h>
/* Textually include all three analytic backends with renamed default
* constructors so one binary can exercise every provider (the files each
* define spacetime_create_default, which would collide at link time). */
#define spacetime_create_default spacetime_create_default_minkowski_local
#include "../../src/spacetime_minkowski.c"
#undef spacetime_create_default
#define spacetime_create_default spacetime_create_default_schwarzschild_local
#include "../../src/spacetime_schwarzschild.c"
#undef spacetime_create_default
#define spacetime_create_default spacetime_create_default_alcubierre_local
#include "../../src/spacetime_alcubierre.c"
#undef spacetime_create_default
#include "geodesic.h"
#include "observer.h"
#include "spacetime.h"
#define PI 3.14159265358979323846
#define MAX_CASES 16
#define MAX_DIRS 16
#define RAY_CAP 2500
static long g_rays = 0;
/* ------------------------------------------------------------------ */
/* Small helpers */
/* ------------------------------------------------------------------ */
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 double ang_delta(const double a[3], const double b[3]) {
const double na = sqrt(dot3(a, a)), nb = sqrt(dot3(b, b));
if (!(na > 0.0) || !(nb > 0.0))
return NAN;
const double ca = dot3(a, b) / (na * nb);
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]};
const double sn = sqrt(dot3(cross, cross)) / (na * nb);
return atan2(sn, ca);
}
static double now_s(void) {
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC, &ts);
return (double)ts.tv_sec + 1e-9 * (double)ts.tv_nsec;
}
static int is_dark(int outcome) { return outcome == RAY_OUTCOME_DARK; }
static int is_escaped(int outcome) { return outcome == RAY_OUTCOME_ESCAPED; }
/* Critical local angle (static observer, M=1) for the Schwarzschild monopole.
* theta_c = asin(3 sqrt(3) sqrt(1-2/r) / r). */
static double critical_angle(double r) {
const double b = 3.0 * sqrt(3.0) * sqrt(1.0 - 2.0 / r) / r;
if (!(b < 1.0))
return PI / 2.0;
return asin(b);
}
/* ------------------------------------------------------------------ */
/* Case / result / comparison structures */
/* ------------------------------------------------------------------ */
typedef struct {
char name[32];
int prov; /* 0 schwarzschild, 1 minkowski, 2 alcubierre */
double mass, escape_radius;
double alc_vs, alc_radius, alc_sigma;
double look_ra_deg, look_dec_deg;
double position[3], velocity[3];
double initial_step, lookback;
unsigned int max_steps;
unsigned int ref_max_steps; /* budget for the tol=1e-12 reference only */
double max_step_for_min_scan; /* fixed upper during the min scan */
int n_dirs;
double dirs[MAX_DIRS][3];
double theta[MAX_DIRS]; /* NAN for non-angle direction sets */
} Case;
typedef struct {
int outcome, reason, end_id;
double stop_t, g, thr;
double n[3];
unsigned int steps, rejected;
unsigned long rhs;
double wall;
} Res;
typedef struct {
int has_ref;
int class_match;
double dn_ang, dlogg, dgrel, dstopT, dmargin;
} Cmp;
typedef struct {
long n, class_mismatch;
double max_dn_ang, max_dlogg, max_dgrel, max_dstopT, max_dmargin;
unsigned long sum_rhs, max_rhs;
unsigned int max_rejected;
double sum_wall;
long escaped, dark, unresolved, incomplete;
} Agg;
static void agg_init(Agg *a) {
memset(a, 0, sizeof *a);
a->max_dn_ang = a->max_dlogg = a->max_dgrel = a->max_dstopT =
a->max_dmargin = -1.0;
}
static void bump(double *m, double v) {
if (!isfinite(v))
return;
if (v > *m)
*m = v;
}
static void agg_add(Agg *a, const Res *r, const Cmp *c) {
++a->n;
if (c->has_ref && !c->class_match)
++a->class_mismatch;
if (c->has_ref) {
bump(&a->max_dn_ang, c->dn_ang);
bump(&a->max_dlogg, c->dlogg);
bump(&a->max_dgrel, c->dgrel);
bump(&a->max_dstopT, c->dstopT);
bump(&a->max_dmargin, c->dmargin);
}
a->sum_rhs += r->rhs;
if (r->rhs > a->max_rhs)
a->max_rhs = r->rhs;
if (r->rejected > a->max_rejected)
a->max_rejected = r->rejected;
a->sum_wall += r->wall;
if (is_dark(r->outcome))
++a->dark;
else if (is_escaped(r->outcome))
++a->escaped;
else if (r->outcome == RAY_OUTCOME_UNRESOLVED)
++a->unresolved;
else
++a->incomplete;
}
/* ------------------------------------------------------------------ */
/* Config builders */
/* ------------------------------------------------------------------ */
static GeodesicTraceConfig make_dp(double initial, double min_step,
double max_step, double tol,
unsigned int max_steps, double lookback) {
GeodesicTraceConfig c;
memset(&c, 0, sizeof c);
c.coordinate_time_step = initial;
c.max_steps = max_steps;
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 = tol;
c.atol_Pi = tol;
c.atol_L = tol;
c.rtol = tol;
c.min_step = min_step;
c.max_step = max_step;
c.consecutive_rejection_limit = 32;
c.max_lookback_time = lookback;
return c;
}
static GeodesicTraceConfig make_rk4(double step, unsigned int max_steps) {
GeodesicTraceConfig c;
memset(&c, 0, sizeof c);
c.coordinate_time_step = step;
c.max_steps = max_steps;
c.threshold.kind = THRESHOLD_LOG_ENERGY_GROWTH;
c.threshold.value = 8.0;
c.threshold.policy_version = 3;
c.stepper = GEODESIC_STEPPER_RK4;
return c;
}
/* ------------------------------------------------------------------ */
/* Trace + compare */
/* ------------------------------------------------------------------ */
static Res run_ray(const SpacetimeSource *source, const ObserverState *observer,
const double dir[3], const GeodesicTraceConfig *config) {
Res r;
memset(&r, 0, sizeof r);
++g_rays;
if (g_rays > RAY_CAP) {
fprintf(stderr, "FATAL: ray cap %d exceeded\n", RAY_CAP);
exit(3);
}
const double t0 = now_s();
const RayEndpoint e = geodesic_trace_past(source, observer, dir, config);
r.wall = now_s() - t0;
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;
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 void compare(const Res *a, const Res *ref, Cmp *c) {
memset(c, 0, sizeof *c);
c->has_ref = 1;
c->class_match = (a->outcome == ref->outcome && a->reason == ref->reason);
if (is_escaped(a->outcome) && is_escaped(ref->outcome)) {
c->dn_ang = ang_delta(a->n, ref->n);
if (a->g > 0.0 && ref->g > 0.0 && isfinite(a->g) && isfinite(ref->g)) {
c->dlogg = fabs(log(a->g) - log(ref->g));
c->dgrel = fabs(a->g / ref->g - 1.0);
} else {
c->dlogg = c->dgrel = NAN;
}
} else {
c->dn_ang = c->dlogg = c->dgrel = NAN;
}
c->dstopT = (isfinite(a->stop_t) && isfinite(ref->stop_t))
? fabs(a->stop_t - ref->stop_t)
: NAN;
c->dmargin = (is_dark(a->outcome) && is_dark(ref->outcome) &&
isfinite(a->thr) && isfinite(ref->thr))
? fabs(a->thr - ref->thr)
: NAN;
}
static void write_row(FILE *f, const char *phase, const Case *c, int di,
double upper, double min_step, double tol, const Res *r,
const Cmp *cmp) {
fprintf(f,
"%s,%s,%d,%.17g,%.6g,%.6g,%.6g,%d,%d,%u,%.17g,%u,%u,%lu,%.6g,"
"%.17g,%.17g,%.17g,%.17g,%.17g",
phase, c->name, di, c->theta[di], upper, min_step, tol, r->outcome,
r->reason, r->end_id, r->stop_t, r->steps, r->rejected, r->rhs,
r->wall, r->n[0], r->n[1], r->n[2], r->g, r->thr);
if (cmp->has_ref)
fprintf(f, ",%d,%.17g,%.17g,%.17g,%.17g,%.17g", cmp->class_match,
cmp->dn_ang, cmp->dlogg, cmp->dgrel, cmp->dstopT, cmp->dmargin);
else
fprintf(f, ",,nan,nan,nan,nan,nan");
fprintf(f, "\n");
}
static const char *ROW_HEADER =
"phase,case,dir,theta,upper,min_step,tol,outcome,reason,end_id,stop_t,"
"steps,rejected,rhs,wall_s,nx,ny,nz,g,thr,class_match,dn_ang,dlogg,dgrel,"
"dstopT,dmargin\n";
static void write_agg(FILE *f, const char *phase, const Case *c, double upper,
double min_step, double tol, const Agg *a) {
fprintf(f,
"%s,%s,%.6g,%.6g,%.6g,%ld,%ld,%.6g,%.6g,%.6g,%.6g,%.6g,%lu,%lu,%u,"
"%.4f,%ld,%ld,%ld,%ld\n",
phase, c->name, upper, min_step, tol, a->n, a->class_mismatch,
a->max_dn_ang, a->max_dlogg, a->max_dgrel, a->max_dstopT,
a->max_dmargin, a->sum_rhs, a->max_rhs, a->max_rejected, a->sum_wall,
a->escaped, a->dark, a->unresolved, a->incomplete);
}
static const char *AGG_HEADER =
"phase,case,upper,min_step,tol,n,class_mismatch,max_dn_ang,max_dlogg,"
"max_dgrel,max_dstopT,max_dmargin,sum_rhs,max_rhs,max_rejected,sum_wall,"
"escaped,dark,unresolved,incomplete\n";
/* ------------------------------------------------------------------ */
/* Case builders */
/* ------------------------------------------------------------------ */
static void dir_from_theta(double th, double n[3]) {
if (th == 0.0) {
n[0] = 1.0;
n[1] = n[2] = 0.0;
} else if (th == PI) {
n[0] = -1.0;
n[1] = n[2] = 0.0;
} else {
n[0] = cos(th);
n[1] = sin(th);
n[2] = 0.0;
}
const double nn = sqrt(dot3(n, n));
for (int i = 0; i < 3; ++i)
n[i] /= nn;
}
static void fill_static_dirs(Case *c, double r) {
const double tc = critical_angle(r);
const double th[MAX_DIRS] = {0.0, PI, PI / 2.0, tc,
tc - 1e-3, tc + 1e-3, tc - 1e-5, tc + 1e-5,
tc - 1e-7, tc + 1e-7, tc + 0.05, tc - 0.05};
c->n_dirs = 12;
for (int i = 0; i < c->n_dirs; ++i) {
c->theta[i] = th[i];
dir_from_theta(th[i], c->dirs[i]);
}
}
static void case_static(Case *c, const char *name, double r, double look_ra) {
memset(c, 0, sizeof *c);
snprintf(c->name, sizeof c->name, "%s", name);
c->prov = 0;
c->mass = 1.0;
c->escape_radius = 256.0;
c->look_ra_deg = look_ra;
c->look_dec_deg = 0.0;
c->position[0] = r;
c->position[1] = c->position[2] = 0.0;
c->velocity[0] = c->velocity[1] = c->velocity[2] = 0.0;
c->initial_step = 0.1;
c->lookback = 6553.6;
c->max_steps = 65536;
c->ref_max_steps = 65536;
c->max_step_for_min_scan = 2.0;
fill_static_dirs(c, r);
}
static void case_freefall(Case *c, const char *name, double r) {
memset(c, 0, sizeof *c);
snprintf(c->name, sizeof c->name, "%s", name);
c->prov = 0;
c->mass = 1.0;
c->escape_radius = 256.0;
c->look_ra_deg = 180.0;
c->look_dec_deg = 0.0;
c->position[0] = r;
c->position[1] = c->position[2] = 0.0;
/* Free-fall from rest at infinity, coordinate velocity dx/dt (the formula
* used by local/p2d_validation/observer_endpoints.c). */
const double y = sqrt(2.0 / r);
c->velocity[0] = -y * (1.0 + y) / (1.0 + y + y * y);
c->velocity[1] = c->velocity[2] = 0.0;
c->initial_step = 0.1;
c->lookback = 6553.6;
c->max_steps = 65536;
c->ref_max_steps = 65536;
c->max_step_for_min_scan = 2.0;
{
const double th[4] = {0.0, PI / 2.0, PI, 3.0 * PI / 4.0};
c->n_dirs = 4;
for (int i = 0; i < 4; ++i) {
c->theta[i] = th[i];
dir_from_theta(th[i], c->dirs[i]);
}
}
}
static void case_minkowski(Case *c) {
memset(c, 0, sizeof *c);
snprintf(c->name, sizeof c->name, "mink_moving");
c->prov = 1;
c->escape_radius = 64.0;
c->look_ra_deg = 0.0;
c->look_dec_deg = 0.0;
c->position[0] = 10.0;
c->position[1] = 20.0;
c->position[2] = -15.0;
c->velocity[0] = 0.3;
c->velocity[1] = 0.2;
c->velocity[2] = 0.1;
c->initial_step = 1.0;
c->lookback = 2048.0;
c->max_steps = 2048;
c->ref_max_steps = 4096;
c->max_step_for_min_scan = 16.0;
{
const double d[6][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},
{0.6, 0.8, 0.0}, {-0.6, 0.8, 0.0}};
c->n_dirs = 6;
for (int i = 0; i < 6; ++i) {
c->theta[i] = NAN;
for (int k = 0; k < 3; ++k)
c->dirs[i][k] = d[i][k];
}
}
}
static void case_alcubierre(Case *c, const char *name, double vs, double sigma) {
memset(c, 0, sizeof *c);
snprintf(c->name, sizeof c->name, "%s", name);
c->prov = 2;
c->alc_vs = vs;
c->alc_radius = 1.0;
c->alc_sigma = sigma;
c->look_ra_deg = 0.0;
c->look_dec_deg = 0.0;
c->position[0] = c->position[1] = c->position[2] = 0.0;
c->velocity[0] = vs; /* comoving with the bubble center */
c->velocity[1] = c->velocity[2] = 0.0;
c->initial_step = fmin(0.1, 0.05 / sigma);
c->max_steps = 100000;
c->ref_max_steps = 200000;
{
const double esc = spacetime_alcubierre_escape_radius(1.0, sigma);
c->lookback = 1.25 * 4.0 * esc / (1.0 - fabs(vs));
}
c->max_step_for_min_scan = c->initial_step * 8.0;
{
const double d[6][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},
{0.6, 0.8, 0.0}, {-0.6, 0.8, 0.0}};
c->n_dirs = 6;
for (int i = 0; i < 6; ++i) {
c->theta[i] = NAN;
for (int k = 0; k < 3; ++k)
c->dirs[i][k] = d[i][k];
}
}
}
/* ------------------------------------------------------------------ */
/* Source/observer construction */
/* ------------------------------------------------------------------ */
static int build_source(const Case *c, SpacetimeSource *s) {
if (c->prov == 0)
return spacetime_create_schwarzschild_ks(s, c->mass, c->escape_radius);
if (c->prov == 1)
return spacetime_create_minkowski(s, c->escape_radius);
return spacetime_create_alcubierre(s, c->alc_vs, c->alc_radius,
c->alc_sigma);
}
static int build_observer(const Case *c, const SpacetimeSource *s,
ObserverState *o) {
MetricData m;
if (spacetime_eval(s, 0.0, c->position, &m) != SPACETIME_POINT_OK)
return -1;
ObserverCamera cam;
memset(&cam, 0, sizeof cam);
cam.coordinate_time = 0.0;
for (int i = 0; i < 3; ++i) {
cam.position[i] = c->position[i];
cam.velocity[i] = c->velocity[i];
}
cam.look_ra_deg = c->look_ra_deg;
cam.look_dec_deg = c->look_dec_deg;
cam.roll_deg = 0.0;
return observer_from_coordinate_camera(&m, &cam, o, NULL) ==
OBSERVER_BUILD_OK
? 0
: -1;
}
/* ------------------------------------------------------------------ */
/* Generic scan driver */
/* ------------------------------------------------------------------ */
static void run_grid(const Case *cases, int ncases, const double *uppers,
int n_up, const double *floors, int n_fl, double tol,
const char *phase, const char *rows_path,
const char *agg_path, Res refs[][MAX_DIRS]) {
FILE *rf = fopen(rows_path, "w");
FILE *af = fopen(agg_path, "w");
if (!rf || !af) {
fprintf(stderr, "FATAL: cannot open %s / %s\n", rows_path, agg_path);
exit(2);
}
fprintf(rf, "%s", ROW_HEADER);
fprintf(af, "%s", AGG_HEADER);
for (int ci = 0; ci < ncases; ++ci) {
const Case *c = &cases[ci];
SpacetimeSource source;
ObserverState observer;
if (build_source(c, &source) || build_observer(c, &source, &observer)) {
fprintf(stderr, "FATAL: cannot build case %s\n", c->name);
exit(2);
}
for (int ui = 0; ui < n_up; ++ui) {
for (int fi = 0; fi < n_fl; ++fi) {
Agg agg;
agg_init(&agg);
for (int di = 0; di < c->n_dirs; ++di) {
const GeodesicTraceConfig cfg =
make_dp(c->initial_step, floors[fi], uppers[ui], tol, c->max_steps,
c->lookback);
const Res r = run_ray(&source, &observer, c->dirs[di], &cfg);
Cmp cmp;
compare(&r, &refs[ci][di], &cmp);
write_row(rf, phase, c, di, uppers[ui], floors[fi], tol, &r, &cmp);
agg_add(&agg, &r, &cmp);
}
write_agg(af, phase, c, uppers[ui], floors[fi], tol, &agg);
printf("SUMMARY %s %s upper=%.6g min=%.6g n=%ld mismatch=%ld "
"max_dn_ang=%.6g max_dgrel=%.6g max_dstopT=%.6g sum_rhs=%lu "
"esc=%ld dark=%ld unres=%ld inc=%ld wall=%.3fs\n",
phase, c->name, uppers[ui], floors[fi], agg.n,
agg.class_mismatch, agg.max_dn_ang, agg.max_dgrel,
agg.max_dstopT, agg.sum_rhs, agg.escaped, agg.dark,
agg.unresolved, agg.incomplete, agg.sum_wall);
fflush(stdout);
}
}
spacetime_destroy(&source);
}
fclose(rf);
fclose(af);
}
/* Compute DP tol=1e-12 reference for every direction of every case.
* upper_ref <= 0 means "use the case's own initial step" as the max_step. */
static void compute_refs(const Case *cases, int ncases, double upper_ref,
const char *path, Res refs[][MAX_DIRS]) {
FILE *f = fopen(path, "w");
if (!f) {
fprintf(stderr, "FATAL: cannot open %s\n", path);
exit(2);
}
fprintf(f, "%s", ROW_HEADER);
for (int ci = 0; ci < ncases; ++ci) {
const Case *c = &cases[ci];
SpacetimeSource source;
ObserverState observer;
if (build_source(c, &source) || build_observer(c, &source, &observer)) {
fprintf(stderr, "FATAL: cannot build case %s\n", c->name);
exit(2);
}
const double upper = (upper_ref > 0.0) ? upper_ref : c->initial_step;
for (int di = 0; di < c->n_dirs; ++di) {
const GeodesicTraceConfig cfg = make_dp(c->initial_step, 1e-12, upper,
1e-12, c->ref_max_steps,
c->lookback);
refs[ci][di] = run_ray(&source, &observer, c->dirs[di], &cfg);
Cmp none;
memset(&none, 0, sizeof none);
write_row(f, "ref", c, di, upper, 1e-12, 1e-12, &refs[ci][di], &none);
printf("REF %s dir=%d theta=%.6g outcome=%d reason=%d stop_t=%.9g "
"steps=%u rejected=%u rhs=%lu g=%.12g thr=%.12g\n",
c->name, di, c->theta[di], refs[ci][di].outcome,
refs[ci][di].reason, refs[ci][di].stop_t, refs[ci][di].steps,
refs[ci][di].rejected, refs[ci][di].rhs, refs[ci][di].g,
refs[ci][di].thr);
fflush(stdout);
}
spacetime_destroy(&source);
}
fclose(f);
}
/* ------------------------------------------------------------------ */
/* Sub-command: Schwarzschild */
/* ------------------------------------------------------------------ */
static void cmd_schwarzschild(void) {
static Case cases[MAX_CASES];
static Res refs[MAX_CASES][MAX_DIRS];
memset(refs, 0, sizeof refs);
int nc = 0;
case_static(&cases[nc++], "r30_inward", 30.0, 180.0);
case_static(&cases[nc++], "r100_inward", 100.0, 180.0);
case_static(&cases[nc++], "r2p1_outward", 2.1, 0.0);
case_freefall(&cases[nc++], "r1p5_freefall", 1.5);
compute_refs(cases, nc, 0.25, "raw/a_sch_reference.csv", refs);
static const double uppers[9] = {0.25, 0.5, 1.0, 2.0, 4.0,
8.0, 16.0, 32.0, 64.0};
static const double floors[10] = {1e-1, 1e-2, 1e-3, 1e-4, 1e-5,
1e-6, 1e-8, 1e-10, 1e-12, 1e-14};
run_grid(cases, nc, uppers, 9, (double[]){1e-12}, 1, 1e-9, "upper",
"raw/a_sch_upper.csv", "raw/a_sch_upper_summary.csv", refs);
run_grid(cases, nc, (double[]){2.0}, 1, floors, 10, 1e-9, "min",
"raw/a_sch_min.csv", "raw/a_sch_min_summary.csv", refs);
/* Sensitive-ray RK4 h-halving reference (extra independent check).
* (case index, dir index): r30 and r2.1, theta_c +- 1e-5 / +- 1e-7. */
const int sens[8][2] = {{0, 8}, {0, 9}, {0, 6}, {0, 7},
{2, 8}, {2, 9}, {2, 6}, {2, 7}};
FILE *f = fopen("raw/a_sch_rk4_sensitive.csv", "w");
if (!f) {
fprintf(stderr, "FATAL: cannot open rk4 csv\n");
exit(2);
}
fprintf(f, "%s", ROW_HEADER);
for (int k = 0; k < 8; ++k) {
const Case *c = &cases[sens[k][0]];
const int di = sens[k][1];
SpacetimeSource source;
ObserverState observer;
if (build_source(c, &source) || build_observer(c, &source, &observer)) {
fprintf(stderr, "FATAL: cannot build case %s\n", c->name);
exit(2);
}
GeodesicTraceConfig rk01 = make_rk4(0.01, 262144);
GeodesicTraceConfig rk005 = make_rk4(0.005, 262144);
Res r01 = run_ray(&source, &observer, c->dirs[di], &rk01);
Res r005 = run_ray(&source, &observer, c->dirs[di], &rk005);
Cmp hh, vsdp;
compare(&r005, &r01, &hh);
compare(&r005, &refs[sens[k][0]][di], &vsdp);
Cmp none;
memset(&none, 0, sizeof none);
write_row(f, "rk4_0.01", c, di, 0.01, 0.0, 0.0, &r01, &none);
write_row(f, "rk4_0.005", c, di, 0.005, 0.0, 0.0, &r005, &none);
printf("RK4 %s dir=%d theta=%.6g o01=%d/%d o005=%d/%d "
"hhalve_dn=%.6g hhalve_dgrel=%.6g vsdp_dn=%.6g vsdp_dgrel=%.6g "
"vsdp_class=%d\n",
c->name, di, c->theta[di], r01.outcome, r01.reason, r005.outcome,
r005.reason, hh.dn_ang, hh.dgrel, vsdp.dn_ang, vsdp.dgrel,
vsdp.class_match);
spacetime_destroy(&source);
}
fclose(f);
}
/* ------------------------------------------------------------------ */
/* Sub-command: Minkowski (analytic flat verification) */
/* ------------------------------------------------------------------ */
/* Analytic flat past escape from an observer at (t0,x0) with tetrad-frame
* direction n: coordinate photon velocity u = k_vec/k^0, crossing radius R at
* s = t-t0 < 0, sky direction n_inf = -u, frequency ratio g = 1/k^0. */
static int mink_analytic(const Case *c, const ObserverState *o,
const double n[3], double R, double *t_cross,
double *xc, double *ninf, double *g) {
double k[4];
for (int mu = 0; mu < 4; ++mu)
k[mu] = o->tetrad[0][mu];
for (int a = 0; a < 3; ++a)
for (int mu = 0; mu < 4; ++mu)
k[mu] -= n[a] * o->tetrad[a + 1][mu];
if (!(k[0] > 0.0))
return -1;
double u[3];
for (int i = 0; i < 3; ++i)
u[i] = k[i + 1] / k[0];
const double uu = dot3(u, u);
const double B = 2.0 * dot3(c->position, u);
const double C = dot3(c->position, c->position) - R * R;
const double disc = B * B - 4.0 * uu * C;
if (disc < 0.0)
return -1;
const double sq = sqrt(disc);
const double s1 = (-B + sq) / (2.0 * uu);
const double s2 = (-B - sq) / (2.0 * uu);
const double s = (s1 < 0.0) ? fmin(s1, s2) : s2;
if (!(s < 0.0))
return -1;
*t_cross = o->coordinate_time + s;
for (int i = 0; i < 3; ++i)
xc[i] = c->position[i] + u[i] * s;
const double nu = sqrt(uu);
for (int i = 0; i < 3; ++i)
ninf[i] = -u[i] / nu;
*g = 1.0 / k[0];
return 0;
}
static void cmd_minkowski(void) {
Case c;
case_minkowski(&c);
static Res refs[MAX_CASES][MAX_DIRS];
memset(refs, 0, sizeof refs);
SpacetimeSource source;
ObserverState observer;
if (build_source(&c, &source) || build_observer(&c, &source, &observer)) {
fprintf(stderr, "FATAL: cannot build minkowski case\n");
exit(2);
}
FILE *mf = fopen("raw/a_mink_analytic.csv", "w");
if (!mf) {
fprintf(stderr, "FATAL: cannot open analytic csv\n");
exit(2);
}
fprintf(mf,
"case,dir,upper,min_step,tol,outcome,stop_t,x_err,t_err,dn_ang,"
"dgrel,g_analytic,stop_t_analytic\n");
/* Reference: tol=1e-12, max_step = initial (1.0). */
const GeodesicTraceConfig refcfg =
make_dp(c.initial_step, 1e-12, c.initial_step, 1e-12, c.max_steps,
c.lookback);
for (int di = 0; di < c.n_dirs; ++di)
refs[0][di] = run_ray(&source, &observer, c.dirs[di], &refcfg);
static const double uppers[5] = {1.0, 4.0, 16.0, 64.0, 256.0};
static const double floors[3] = {1e-1, 1e-6, 1e-12};
for (int ui = 0; ui < 5; ++ui) {
for (int fi = 0; fi < 3; ++fi) {
Agg agg;
agg_init(&agg);
for (int di = 0; di < c.n_dirs; ++di) {
const GeodesicTraceConfig cfg =
make_dp(c.initial_step, floors[fi], uppers[ui], 1e-9, c.max_steps,
c.lookback);
const Res r = run_ray(&source, &observer, c.dirs[di], &cfg);
Cmp cmp;
compare(&r, &refs[0][di], &cmp);
agg_add(&agg, &r, &cmp);
double t_an, x_an[3], n_an[3], g_an;
const int ok = mink_analytic(&c, &observer, c.dirs[di], c.escape_radius,
&t_an, x_an, n_an, &g_an);
double x_err = NAN, t_err = NAN, dn_an = NAN, dg_an = NAN;
if (ok == 0 && is_escaped(r.outcome)) {
/* Crossing position is checked through the direct trace in
* a_mink_analytic_x.csv; here compare sky direction, g and time. */
t_err = fabs(r.stop_t - t_an);
dn_an = ang_delta(r.n, n_an);
if (r.g > 0.0 && g_an > 0.0)
dg_an = fabs(r.g / g_an - 1.0);
}
fprintf(mf,
"%s,%d,%.6g,%.6g,%.6g,%d,%.17g,%.17g,%.17g,%.17g,%.17g,%.17g,"
"%.17g\n",
c.name, di, uppers[ui], floors[fi], 1e-9, r.outcome, r.stop_t,
x_err, t_err, dn_an, dg_an, g_an, t_an);
}
printf("SUMMARY mink upper=%.6g min=%.6g n=%ld mismatch=%ld "
"max_dn_ang=%.6g max_dgrel=%.6g sum_rhs=%lu esc=%ld dark=%ld "
"unres=%ld inc=%ld\n",
uppers[ui], floors[fi], agg.n, agg.class_mismatch,
agg.max_dn_ang, agg.max_dgrel, agg.sum_rhs, agg.escaped, agg.dark,
agg.unresolved, agg.incomplete);
fflush(stdout);
}
}
/* Direct endpoint check for the crossing position (x) on a few configs. */
FILE *xf = fopen("raw/a_mink_analytic_x.csv", "w");
if (!xf) {
fprintf(stderr, "FATAL: cannot open analytic x csv\n");
exit(2);
}
fprintf(xf, "case,dir,upper,min_step,tol,outcome,x_err_x,x_err_y,x_err_z,"
"t_err\n");
for (int ui = 0; ui < 5; ++ui) {
for (int fi = 0; fi < 3; ++fi) {
for (int di = 0; di < c.n_dirs; ++di) {
const GeodesicTraceConfig cfg =
make_dp(c.initial_step, floors[fi], uppers[ui], 1e-9, c.max_steps,
c.lookback);
const RayEndpoint e =
geodesic_trace_past(&source, &observer, c.dirs[di], &cfg);
++g_rays;
double t_an, x_an[3], n_an[3], g_an;
const int ok = mink_analytic(&c, &observer, c.dirs[di], c.escape_radius,
&t_an, x_an, n_an, &g_an);
double ex = NAN, ey = NAN, ez = NAN, te = NAN;
if (ok == 0 && e.outcome == RAY_OUTCOME_ESCAPED) {
ex = fabs(e.final_x[0] - x_an[0]);
ey = fabs(e.final_x[1] - x_an[1]);
ez = fabs(e.final_x[2] - x_an[2]);
te = fabs(e.stop_coordinate_time - t_an);
}
fprintf(xf, "%s,%d,%.6g,%.6g,%.6g,%d,%.17g,%.17g,%.17g,%.17g\n",
c.name, di, uppers[ui], floors[fi], 1e-9, e.outcome, ex, ey, ez,
te);
}
}
}
fclose(xf);
fclose(mf);
spacetime_destroy(&source);
}
/* ------------------------------------------------------------------ */
/* Sub-command: Alcubierre */
/* ------------------------------------------------------------------ */
static void cmd_alcubierre(void) {
static Case cases[MAX_CASES];
static Res refs[MAX_CASES][MAX_DIRS];
memset(refs, 0, sizeof refs);
int nc = 0;
case_alcubierre(&cases[nc++], "alc_v3_s1", 0.3, 1.0);
case_alcubierre(&cases[nc++], "alc_v3_s10", 0.3, 10.0);
case_alcubierre(&cases[nc++], "alc_v9_s1", 0.9, 1.0);
case_alcubierre(&cases[nc++], "alc_v9_s10", 0.9, 10.0);
case_alcubierre(&cases[nc++], "alc_v9_s100", 0.9, 100.0);
compute_refs(cases, nc, -1.0, "raw/a_alc_reference.csv", refs);
/* Standard grid: upper = initial*{1,4,8,16,32,64}, floor {1e-4,1e-8,1e-12}. */
{
static const double floors[3] = {1e-4, 1e-8, 1e-12};
Case sub[4];
for (int i = 0; i < 4; ++i)
sub[i] = cases[i];
/* run_grid uses one upper array for all cases, so use per-case initial by
* calling run_grid four times. */
for (int i = 0; i < 4; ++i) {
char rp[128], ap[128];
snprintf(rp, sizeof rp, "raw/a_alc_%.31s_upper.csv", sub[i].name);
snprintf(ap, sizeof ap, "raw/a_alc_%.31s_upper_summary.csv", sub[i].name);
const double base = sub[i].initial_step;
const double up[6] = {base, 4 * base, 8 * base,
16 * base, 32 * base, 64 * base};
run_grid(&sub[i], 1, up, 6, floors, 3, 1e-9, "upper", rp, ap, &refs[i]);
}
}
/* Focused lower-bound stress: vs=.9, sigma=100, upper = initial*8 = 0.004,
* floors {1e-4,1e-5,1e-6,1e-8,1e-12,1e-14}. */
{
const Case *c = &cases[4];
const double up[1] = {c->initial_step * 8.0};
const double floors[6] = {1e-4, 1e-5, 1e-6, 1e-8, 1e-12, 1e-14};
run_grid(c, 1, up, 1, floors, 6, 1e-9, "min_stress",
"raw/a_alc_v9_s100_min.csv", "raw/a_alc_v9_s100_min_summary.csv",
&refs[4]);
}
}
/* ------------------------------------------------------------------ */
int main(int argc, char **argv) {
if (argc < 2) {
fprintf(stderr,
"usage: %s schwarzschild|minkowski|alcubierre\n",
argv[0]);
return 1;
}
if (!strcmp(argv[1], "schwarzschild"))
cmd_schwarzschild();
else if (!strcmp(argv[1], "minkowski"))
cmd_minkowski();
else if (!strcmp(argv[1], "alcubierre"))
cmd_alcubierre();
else {
fprintf(stderr, "unknown sub-command %s\n", argv[1]);
return 1;
}
printf("TOTAL_RAYS %ld\n", g_rays);
return 0;
}
Reference in new issue View Git Blame Copy Permalink