AMSS-NCKU/verify_accuracy.py

#!/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()