From 318b5254ccd8816f607f7870f066ee0489ec5495 Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Fri, 27 Feb 2026 17:38:21 +0800 Subject: [PATCH] =?UTF-8?q?=E6=A0=B9=E6=8D=AE=E7=BB=84=E5=A7=94=E4=BC=9A?= =?UTF-8?q?=E9=82=AE=E4=BB=B6=E8=A6=81=E6=B1=82=E6=9B=B4=E6=96=B0=E6=A3=80?= =?UTF-8?q?=E6=B5=8B=E8=84=9A=E6=9C=AC=EF=BC=8C=E5=A2=9E=E5=8A=A0=E5=AF=B9?= =?UTF-8?q?3D=E5=90=91=E9=87=8F=E5=92=8C=E4=B8=89=E4=B8=AA=E5=88=86?= =?UTF-8?q?=E9=87=8F=E5=88=86=E5=88=AB=E6=A3=80=E6=B5=8BRMS=E5=B0=8F?= =?UTF-8?q?=E4=BA=8E1.0%?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- AMSS_NCKU_Verify_ASC26.py | 157 ++++++++++++++++---------------------- 1 file changed, 64 insertions(+), 93 deletions(-) diff --git a/AMSS_NCKU_Verify_ASC26.py b/AMSS_NCKU_Verify_ASC26.py index ed386e7..6ea76d5 100644 --- a/AMSS_NCKU_Verify_ASC26.py +++ b/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_rms_error(bh_data_ref, bh_data_target): +def calculate_all_rms_errors(bh_data_ref, bh_data_target): """ - Calculate trajectory-based RMS error on the XY plane between baseline and optimized simulations. - - This function computes the RMS error independently for BH1 and BH2 trajectories, - 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 + 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. + Returns the maximum error between BH1 and BH2 for each category. """ - # 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 - 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] + results = {} - x1_new = bh_data_target['x1'][:M] - y1_new = bh_data_target['y1'][:M] - x2_new = bh_data_target['x2'][:M] - y2_new = bh_data_target['y2'][:M] + for bh in ['1', '2']: + 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] + 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] - # Calculate RMS for BH1 - delta_r1 = np.sqrt((x1_ref - x1_new)**2 + (y1_ref - y1_new)**2) - r1_ref = np.sqrt(x1_ref**2 + y1_ref**2) - r1_new = np.sqrt(x1_new**2 + y1_new**2) - r1_max = np.maximum(r1_ref, r1_new) + # 核心修改：根据组委会的邮件指示，分母统一使用 r = sqrt(x^2 + y^2) + r_ref = np.sqrt(x_r**2 + y_r**2) + r_new = np.sqrt(x_n**2 + y_n**2) + denom_max = np.maximum(r_ref, r_new) - # Calculate RMS for BH2 - delta_r2 = np.sqrt((x2_ref - x2_new)**2 + (y2_ref - y2_new)**2) - r2_ref = np.sqrt(x2_ref**2 + y2_ref**2) - r2_new = np.sqrt(x2_new**2 + y2_new**2) - r2_max = np.maximum(r2_ref, r2_new) + valid = denom_max > 1e-15 + if np.sum(valid) < 10: + results[f'BH{bh}'] = { '3D_Vector': 0.0, 'X_Component': 0.0, 'Y_Component': 0.0, 'Z_Component': 0.0 } + continue - # Avoid division by zero for BH1 - valid_mask1 = r1_max > 1e-15 - if np.sum(valid_mask1) < 10: - return None, "Insufficient valid data points for BH1" + def calc_rms(delta): + # 将对应分量的偏差除以统一的轨道半径分母 denom_max + return np.sqrt(np.mean((delta[valid] / denom_max[valid])**2)) * 100 - terms1 = (delta_r1[valid_mask1] / r1_max[valid_mask1])**2 - rms_bh1 = np.sqrt(np.mean(terms1)) * 100 + # 1. Total 3D Vector RMS + 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 - valid_mask2 = r2_max > 1e-15 - if np.sum(valid_mask2) < 10: - return None, "Insufficient valid data points for BH2" + # 2. Component-wise RMS (分离计算各轴，但共用半径分母) + rms_x = calc_rms(np.abs(x_r - x_n)) + rms_y = calc_rms(np.abs(y_r - y_n)) + rms_z = calc_rms(np.abs(z_r - z_n)) - terms2 = (delta_r2[valid_mask2] / r2_max[valid_mask2])**2 - rms_bh2 = np.sqrt(np.mean(terms2)) * 100 + results[f'BH{bh}'] = { + '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 - rms_final = max(rms_bh1, rms_bh2) - - return rms_final, None + # 获取 BH1 和 BH2 中的最大误差 + max_rms = { + '3D_Vector': max(results['BH1']['3D_Vector'], results['BH2']['3D_Vector']), + '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_rel, error, threshold=1.0): - """Print RMS error results""" - print(f"\n{Color.BOLD}1. RMS Error Analysis (Baseline vs Optimized){Color.RESET}") - print("-" * 45) +def print_rms_results(rms_dict, error, threshold=1.0): + print(f"\n{Color.BOLD}1. RMS Error Analysis (Maximums of BH1 & BH2){Color.RESET}") + print("-" * 65) 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}%") - print(f" Requirement: < {threshold}%") - print(f" Status: {get_status_text(passed)}") - - return passed + for key, val in rms_dict.items(): + passed = val < threshold + all_passed = all_passed and passed + status = get_status_text(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 - rms_rel, error = calculate_rms_error(bh_data_ref, bh_data_target) - rms_passed = print_rms_results(rms_rel, error) + # Output modified RMS results + rms_dict, error = calculate_all_rms_errors(bh_data_ref, bh_data_target) + 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()