#!/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 # ANSI Color Codes class Color: GREEN = '\033[92m' RED = '\033[91m' YELLOW = '\033[93m' BLUE = '\033[94m' BOLD = '\033[1m' RESET = '\033[0m' def get_status_text(passed): if passed: return f"{Color.GREEN}{Color.BOLD}PASS{Color.RESET}" else: return f"{Color.RED}{Color.BOLD}FAIL{Color.RESET}" 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("\n" + Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET) print(Color.BOLD + " AMSS-NCKU GW150914 Simulation Accuracy Verification Report" + Color.RESET) print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET) def print_rms_results(rms_abs, rms_rel, error, threshold=1.0): """Print RMS error results""" print(f"\n{Color.BOLD}1. RMS Error Analysis (Black Hole Trajectory){Color.RESET}") print("-" * 45) if error: print(f" {Color.RED}Error: {error}{Color.RESET}") return False passed = rms_rel < threshold print(f" RMS absolute error: {rms_abs:.4e} M") print(f" RMS relative error: {rms_rel:.4f}%") print(f" Requirement: < {threshold}%") print(f" Status: {get_status_text(passed)}") return passed def print_constraint_results(results, threshold=2.0): """Print constraint violation results""" print(f"\n{Color.BOLD}2. ADM Constraint Violation Analysis (Grid Level 0){Color.RESET}") print("-" * 45) names = ['Ham', 'Px', 'Py', 'Pz', 'Gx', 'Gy', 'Gz'] for i, name in enumerate(names): print(f" Max |{name:3}|: {results[name]:.6f}", end=" ") if (i + 1) % 2 == 0: print() if len(names) % 2 != 0: print() passed = results['max_violation'] < threshold print(f"\n Maximum violation: {results['max_violation']:.6f}") print(f" Requirement: < {threshold}") print(f" Status: {get_status_text(passed)}") return passed def print_summary(rms_passed, constraint_passed): """Print summary""" print("\n" + Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET) print(Color.BOLD + "Verification Summary" + Color.RESET) print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET) all_passed = rms_passed and constraint_passed res_rms = get_status_text(rms_passed) res_con = get_status_text(constraint_passed) print(f" [1] RMS trajectory check: {res_rms}") print(f" [2] ADM constraint check: {res_con}") final_status = f"{Color.GREEN}{Color.BOLD}ALL CHECKS PASSED{Color.RESET}" if all_passed else f"{Color.RED}{Color.BOLD}SOME CHECKS FAILED{Color.RESET}" print(f"\n Overall result: {final_status}") print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET + "\n") 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"{Color.RED}{Color.BOLD}Error:{Color.RESET} Black hole trajectory file not found: {bh_file}") sys.exit(1) if not os.path.exists(constraint_file): print(f"{Color.RED}{Color.BOLD}Error:{Color.RESET} Constraint data file not found: {constraint_file}") sys.exit(1) # Print header print_header() print(f"\n{Color.BOLD}Target Directory:{Color.RESET} {Color.BLUE}{output_dir}{Color.RESET}") # 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()