根据组委会邮件要求更新检测脚本，增加对3D向量和三个分量分别检测RMS小于1.0%

AMSS_NCKU_Verify_ASC26.py

@@ -1,9 +1,13 @@
 #!/usr/bin/env python3
 """
-AMSS-NCKU GW150914 Simulation Regression Test Script
+AMSS-NCKU GW150914 Simulation Regression Test Script (Comprehensive Version)
 
 Verification Requirements:
-1. XY-plane trajectory RMS error < 1% (Optimized vs. baseline, max of BH1 and BH2)
+1. RMS errors < 1% for:
+   - 3D Vector Total RMS
+   - X Component RMS
+   - Y Component RMS
+   - Z Component RMS
 2. ADM constraint violation < 2 (Grid Level 0)
 
 RMS Calculation Method:
@@ -57,79 +61,62 @@ def load_constraint_data(filepath):
                 data.append([float(x) for x in parts[:8]])
     return np.array(data)
 
-
+def calculate_all_rms_errors(bh_data_ref, bh_data_target):
-def calculate_rms_error(bh_data_ref, bh_data_target):
     """
-    Calculate trajectory-based RMS error on the XY plane between baseline and optimized simulations.
+    Calculate 3D Vector RMS and component-wise RMS (X, Y, Z) independently.
-
+    Uses r = sqrt(x^2 + y^2) as the denominator for all error normalizations.
-    This function computes the RMS error independently for BH1 and BH2 trajectories,
+    Returns the maximum error between BH1 and BH2 for each category.
-    then returns the maximum of the two as the final RMS error metric.
-
-    For each black hole, the RMS is calculated as:
-        RMS = sqrt( (1/M) * sum( (Δr_i / r_i^max)^2 ) ) × 100%
-
-    where:
-        Δr_i = sqrt((x_ref,i - x_new,i)^2 + (y_ref,i - y_new,i)^2)
-        r_i^max = max(sqrt(x_ref,i^2 + y_ref,i^2), sqrt(x_new,i^2 + y_new,i^2))
-
-    Args:
-        bh_data_ref: Reference (baseline) trajectory data
-        bh_data_target: Target (optimized) trajectory data
-
-    Returns:
-        rms_value: Final RMS error as a percentage (max of BH1 and BH2)
-        error: Error message if any
     """
-    # Align data: truncate to the length of the shorter dataset
     M = min(len(bh_data_ref['time']), len(bh_data_target['time']))
 
     if M < 10:
         return None, "Insufficient data points for comparison"
 
-    # Extract XY coordinates for both black holes
+    results = {}
-    x1_ref = bh_data_ref['x1'][:M]
-    y1_ref = bh_data_ref['y1'][:M]
-    x2_ref = bh_data_ref['x2'][:M]
-    y2_ref = bh_data_ref['y2'][:M]
 
-    x1_new = bh_data_target['x1'][:M]
+    for bh in ['1', '2']:
-    y1_new = bh_data_target['y1'][:M]
+        x_r, y_r, z_r = bh_data_ref[f'x{bh}'][:M], bh_data_ref[f'y{bh}'][:M], bh_data_ref[f'z{bh}'][:M]
-    x2_new = bh_data_target['x2'][:M]
+        x_n, y_n, z_n = bh_data_target[f'x{bh}'][:M], bh_data_target[f'y{bh}'][:M], bh_data_target[f'z{bh}'][:M]
-    y2_new = bh_data_target['y2'][:M]
 
-    # Calculate RMS for BH1
+        # 核心修改：根据组委会的邮件指示，分母统一使用 r = sqrt(x^2 + y^2)
-    delta_r1 = np.sqrt((x1_ref - x1_new)**2 + (y1_ref - y1_new)**2)
+        r_ref = np.sqrt(x_r**2 + y_r**2)
-    r1_ref = np.sqrt(x1_ref**2 + y1_ref**2)
+        r_new = np.sqrt(x_n**2 + y_n**2)
-    r1_new = np.sqrt(x1_new**2 + y1_new**2)
+        denom_max = np.maximum(r_ref, r_new)
-    r1_max = np.maximum(r1_ref, r1_new)
 
-    # Calculate RMS for BH2
+        valid = denom_max > 1e-15
-    delta_r2 = np.sqrt((x2_ref - x2_new)**2 + (y2_ref - y2_new)**2)
+        if np.sum(valid) < 10:
-    r2_ref = np.sqrt(x2_ref**2 + y2_ref**2)
+            results[f'BH{bh}'] = { '3D_Vector': 0.0, 'X_Component': 0.0, 'Y_Component': 0.0, 'Z_Component': 0.0 }
-    r2_new = np.sqrt(x2_new**2 + y2_new**2)
+            continue
-    r2_max = np.maximum(r2_ref, r2_new)
 
-    # Avoid division by zero for BH1
+        def calc_rms(delta):
-    valid_mask1 = r1_max > 1e-15
+            # 将对应分量的偏差除以统一的轨道半径分母 denom_max
-    if np.sum(valid_mask1) < 10:
+            return np.sqrt(np.mean((delta[valid] / denom_max[valid])**2)) * 100
-        return None, "Insufficient valid data points for BH1"
 
-    terms1 = (delta_r1[valid_mask1] / r1_max[valid_mask1])**2
+        # 1. Total 3D Vector RMS
-    rms_bh1 = np.sqrt(np.mean(terms1)) * 100
+        delta_vec = np.sqrt((x_r - x_n)**2 + (y_r - y_n)**2 + (z_r - z_n)**2)
+        rms_3d = calc_rms(delta_vec)
 
-    # Avoid division by zero for BH2
+        # 2. Component-wise RMS (分离计算各轴，但共用半径分母)
-    valid_mask2 = r2_max > 1e-15
+        rms_x = calc_rms(np.abs(x_r - x_n))
-    if np.sum(valid_mask2) < 10:
+        rms_y = calc_rms(np.abs(y_r - y_n))
-        return None, "Insufficient valid data points for BH2"
+        rms_z = calc_rms(np.abs(z_r - z_n))
 
-    terms2 = (delta_r2[valid_mask2] / r2_max[valid_mask2])**2
+        results[f'BH{bh}'] = {
-    rms_bh2 = np.sqrt(np.mean(terms2)) * 100
+            '3D_Vector': rms_3d,
+            'X_Component': rms_x,
+            'Y_Component': rms_y,
+            'Z_Component': rms_z
+        }
 
-    # Final RMS is the maximum of BH1 and BH2
+    # 获取 BH1 和 BH2 中的最大误差
-    rms_final = max(rms_bh1, rms_bh2)
+    max_rms = {
-
+        '3D_Vector': max(results['BH1']['3D_Vector'], results['BH2']['3D_Vector']),
-    return rms_final, None
+        'X_Component': max(results['BH1']['X_Component'], results['BH2']['X_Component']),
+        'Y_Component': max(results['BH1']['Y_Component'], results['BH2']['Y_Component']),
+        'Z_Component': max(results['BH1']['Z_Component'], results['BH2']['Z_Component'])
+    }
 
+    return max_rms, None
 
 def analyze_constraint_violation(constraint_data, n_levels=9):
     """
@@ -155,34 +142,32 @@ def analyze_constraint_violation(constraint_data, n_levels=9):
 
 
 def print_header():
-    """Print report header"""
     print("\n" + Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
-    print(Color.BOLD + "   AMSS-NCKU GW150914 Simulation Regression Test Report" + Color.RESET)
+    print(Color.BOLD + "   AMSS-NCKU GW150914 Comprehensive Regression Test" + Color.RESET)
     print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
 
-
+def print_rms_results(rms_dict, error, threshold=1.0):
-def print_rms_results(rms_rel, error, threshold=1.0):
+    print(f"\n{Color.BOLD}1. RMS Error Analysis (Maximums of BH1 & BH2){Color.RESET}")
-    """Print RMS error results"""
+    print("-" * 65)
-    print(f"\n{Color.BOLD}1. RMS Error Analysis (Baseline vs Optimized){Color.RESET}")
-    print("-" * 45)
 
     if error:
         print(f"   {Color.RED}Error: {error}{Color.RESET}")
         return False
 
-    passed = rms_rel < threshold
+    all_passed = True
+    print(f"   Requirement: < {threshold}%\n")
 
-    print(f"   RMS relative error: {rms_rel:.4f}%")
+    for key, val in rms_dict.items():
-    print(f"   Requirement:        < {threshold}%")
+        passed = val < threshold
-    print(f"   Status:             {get_status_text(passed)}")
+        all_passed = all_passed and passed
-
+        status = get_status_text(passed)
-    return passed
+        print(f"   {key:15}: {val:8.4f}%   |   Status: {status}")
 
+    return all_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)
+    print("-" * 65)
 
     names = ['Ham', 'Px', 'Py', 'Pz', 'Gx', 'Gy', 'Gz']
     for i, name in enumerate(names):
@@ -200,7 +185,6 @@ def print_constraint_results(results, threshold=2.0):
 
 
 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)
@@ -210,7 +194,7 @@ def print_summary(rms_passed, 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"   [1] Comprehensive RMS 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}"
@@ -219,61 +203,48 @@ def print_summary(rms_passed, constraint_passed):
 
     return all_passed
 
-
 def main():
-    # Determine target (optimized) output directory
     if len(sys.argv) > 1:
         target_dir = sys.argv[1]
     else:
         script_dir = os.path.dirname(os.path.abspath(__file__))
         target_dir = os.path.join(script_dir, "GW150914/AMSS_NCKU_output")
 
-    # Determine reference (baseline) directory
     script_dir = os.path.dirname(os.path.abspath(__file__))
     reference_dir = os.path.join(script_dir, "GW150914-origin/AMSS_NCKU_output")
 
-    # Data file paths
     bh_file_ref = os.path.join(reference_dir, "bssn_BH.dat")
     bh_file_target = os.path.join(target_dir, "bssn_BH.dat")
     constraint_file = os.path.join(target_dir, "bssn_constraint.dat")
 
-    # Check if files exist
     if not os.path.exists(bh_file_ref):
         print(f"{Color.RED}{Color.BOLD}Error:{Color.RESET} Baseline trajectory file not found: {bh_file_ref}")
         sys.exit(1)
-
     if not os.path.exists(bh_file_target):
         print(f"{Color.RED}{Color.BOLD}Error:{Color.RESET} Target trajectory file not found: {bh_file_target}")
         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}Reference (Baseline):{Color.RESET} {Color.BLUE}{reference_dir}{Color.RESET}")
     print(f"{Color.BOLD}Target (Optimized):  {Color.RESET} {Color.BLUE}{target_dir}{Color.RESET}")
 
-    # Load data
     bh_data_ref = load_bh_trajectory(bh_file_ref)
     bh_data_target = load_bh_trajectory(bh_file_target)
     constraint_data = load_constraint_data(constraint_file)
 
-    # Calculate RMS error
+    # Output modified RMS results
-    rms_rel, error = calculate_rms_error(bh_data_ref, bh_data_target)
+    rms_dict, error = calculate_all_rms_errors(bh_data_ref, bh_data_target)
-    rms_passed = print_rms_results(rms_rel, error)
+    rms_passed = print_rms_results(rms_dict, error)
 
-    # Analyze constraint violation
+    # Output constraint results
     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()