From 0d24f1503c7f2ec9b88dd777024a850a9beffbd5 Mon Sep 17 00:00:00 2001 From: CGH0S7 Date: Sat, 17 Jan 2026 00:37:30 +0800 Subject: [PATCH] Add accuracy verification script for GW150914 simulation - Verify RMS error < 1% (black hole trajectory vs. post-Newtonian theory) - Verify ADM constraint violation < 2 (Grid Level 0) - Return exit code 0 on pass, 1 on fail Co-Authored-By: Claude Opus 4.5 --- verify_accuracy.py | 216 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 216 insertions(+) create mode 100644 verify_accuracy.py diff --git a/verify_accuracy.py b/verify_accuracy.py new file mode 100644 index 0000000..2b99f2b --- /dev/null +++ b/verify_accuracy.py @@ -0,0 +1,216 @@ +#!/usr/bin/env python3 +""" +AMSS-NCKU GW150914 Simulation Accuracy Verification Script + +Verification Requirements: +1. RMS error < 1% (Black hole trajectory vs. post-Newtonian theory) +2. ADM constraint violation < 2 (Grid Level 0) + +Usage: python3 verify_accuracy.py [output_dir] +Default: output_dir = GW150914/AMSS_NCKU_output +""" + +import numpy as np +import sys +import os + + +def load_bh_trajectory(filepath): + """Load black hole trajectory data""" + data = np.loadtxt(filepath) + return { + 'time': data[:, 0], + 'x1': data[:, 1], 'y1': data[:, 2], 'z1': data[:, 3], + 'x2': data[:, 4], 'y2': data[:, 5], 'z2': data[:, 6] + } + + +def load_constraint_data(filepath): + """Load constraint violation data""" + data = [] + with open(filepath, 'r') as f: + for line in f: + if line.startswith('#'): + continue + parts = line.split() + if len(parts) >= 8: + data.append([float(x) for x in parts[:8]]) + return np.array(data) + + +def calculate_rms_error(bh_data, M1=0.5538461539, M2=0.4461538472): + """ + Calculate RMS error + Compare numerical orbit with post-Newtonian (PN) theory prediction + """ + time = bh_data['time'] + x1, y1, z1 = bh_data['x1'], bh_data['y1'], bh_data['z1'] + x2, y2, z2 = bh_data['x2'], bh_data['y2'], bh_data['z2'] + + # Calculate separation distance + r_num = np.sqrt((x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2) + + # Mass parameters + M_total = M1 + M2 + eta = M1 * M2 / M_total**2 + + # Use inspiral phase data (r > 2M) + inspiral_mask = r_num > 2.0 + t_insp = time[inspiral_mask] + r_insp = r_num[inspiral_mask] + + if len(t_insp) < 10: + return None, None, "Insufficient data points" + + # Post-Newtonian theory prediction + r0 = r_insp[0] + t0 = t_insp[0] + t_coal_PN = (5.0/256.0) * (r0**4) / (eta * M_total**3) + + # Calculate PN predicted separation distance + tau = 1.0 - (t_insp - t0) / t_coal_PN + tau = np.maximum(tau, 1e-10) # Avoid negative values + r_PN = r0 * np.power(tau, 0.25) + + # Use only valid data (t < 0.95 * t_coal) + valid_mask = (t_insp - t0) < t_coal_PN * 0.95 + r_num_valid = r_insp[valid_mask] + r_PN_valid = r_PN[valid_mask] + + # Calculate RMS error + residual = r_num_valid - r_PN_valid + rms_abs = np.sqrt(np.mean(residual**2)) + rms_rel = rms_abs / np.mean(r_num_valid) * 100 + + return rms_abs, rms_rel, None + + +def analyze_constraint_violation(constraint_data, n_levels=9): + """ + Analyze ADM constraint violation + Return maximum constraint violation for Grid Level 0 + """ + # Extract Grid Level 0 data (first entry for each time step) + level0_data = constraint_data[::n_levels] + + # Calculate maximum absolute value for each constraint + results = { + 'Ham': np.max(np.abs(level0_data[:, 1])), + 'Px': np.max(np.abs(level0_data[:, 2])), + 'Py': np.max(np.abs(level0_data[:, 3])), + 'Pz': np.max(np.abs(level0_data[:, 4])), + 'Gx': np.max(np.abs(level0_data[:, 5])), + 'Gy': np.max(np.abs(level0_data[:, 6])), + 'Gz': np.max(np.abs(level0_data[:, 7])) + } + + results['max_violation'] = max(results.values()) + return results + + +def print_header(): + """Print report header""" + print("=" * 60) + print("AMSS-NCKU GW150914 Simulation Accuracy Verification Report") + print("=" * 60) + + +def print_rms_results(rms_abs, rms_rel, error, threshold=1.0): + """Print RMS error results""" + print("\n1. RMS Error Analysis (Black Hole Trajectory)") + print("-" * 40) + + if error: + print(f" Error: {error}") + return False + + print(f" RMS absolute error: {rms_abs:.4f} M") + print(f" RMS relative error: {rms_rel:.4f}%") + print(f" Requirement: < {threshold}%") + + passed = rms_rel < threshold + status = "PASS" if passed else "FAIL" + print(f" Status: {status}") + + return passed + + +def print_constraint_results(results, threshold=2.0): + """Print constraint violation results""" + print("\n2. ADM Constraint Violation Analysis (Grid Level 0)") + print("-" * 40) + + for name in ['Ham', 'Px', 'Py', 'Pz', 'Gx', 'Gy', 'Gz']: + print(f" Max |{name}|: {results[name]:.6f}") + + print(f"\n Maximum constraint violation: {results['max_violation']:.6f}") + print(f" Requirement: < {threshold}") + + passed = results['max_violation'] < threshold + status = "PASS" if passed else "FAIL" + print(f" Status: {status}") + + return passed + + +def print_summary(rms_passed, constraint_passed): + """Print summary""" + print("\n" + "=" * 60) + print("Verification Summary") + print("=" * 60) + + all_passed = rms_passed and constraint_passed + + print(f" RMS error check: {'PASS' if rms_passed else 'FAIL'}") + print(f" Constraint violation check: {'PASS' if constraint_passed else 'FAIL'}") + print(f"\n Overall result: {'All checks passed' if all_passed else 'Some checks failed'}") + + return all_passed + + +def main(): + # Determine output directory + if len(sys.argv) > 1: + output_dir = sys.argv[1] + else: + script_dir = os.path.dirname(os.path.abspath(__file__)) + output_dir = os.path.join(script_dir, "GW150914/AMSS_NCKU_output") + + # Data file paths + bh_file = os.path.join(output_dir, "bssn_BH.dat") + constraint_file = os.path.join(output_dir, "bssn_constraint.dat") + + # Check if files exist + if not os.path.exists(bh_file): + print(f"Error: Black hole trajectory file not found: {bh_file}") + sys.exit(1) + + if not os.path.exists(constraint_file): + print(f"Error: Constraint data file not found: {constraint_file}") + sys.exit(1) + + # Print header + print_header() + print(f"\nData directory: {output_dir}") + + # Load data + bh_data = load_bh_trajectory(bh_file) + constraint_data = load_constraint_data(constraint_file) + + # Calculate RMS error + rms_abs, rms_rel, error = calculate_rms_error(bh_data) + rms_passed = print_rms_results(rms_abs, rms_rel, error) + + # Analyze constraint violation + constraint_results = analyze_constraint_violation(constraint_data) + constraint_passed = print_constraint_results(constraint_results) + + # Print summary + all_passed = print_summary(rms_passed, constraint_passed) + + # Return exit code + sys.exit(0 if all_passed else 1) + + +if __name__ == "__main__": + main()