[TEST]UPSTREAM: Pick some source changes from 48080d0a97

* Sync new folder structure

.idea/vcs.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="VcsDirectoryMappings">
+    <mapping directory="" vcs="Git" />
+  </component>
+</project>

AMSS_NCKU_ABEtest.py

@@ -1,447 +0,0 @@
-
-##################################################################
-##
-## AMSS-NCKU ABE Test Program (Skip TwoPuncture if data exists)
-## Modified from AMSS_NCKU_Program.py
-## Author: Xiaoqu
-## Modified: 2026/02/01
-##
-##################################################################
-
-
-##################################################################
-
-## Print program introduction
-
-import print_information
-
-print_information.print_program_introduction()
-
-##################################################################
-
-import AMSS_NCKU_Input as input_data
-
-##################################################################
-
-## Create directories to store program run data
-
-import os
-import shutil
-import sys
-import time
-
-## Set the output directory according to the input file
-File_directory = os.path.join(input_data.File_directory)
-
-## Check if output directory exists and if TwoPuncture data is available
-#skip_twopuncture = False
-skip_twopuncture = True
-output_directory = os.path.join(File_directory, "AMSS_NCKU_output")
-binary_results_directory = os.path.join(output_directory, input_data.Output_directory)
-
-if os.path.exists(File_directory):
-    print( " Output directory already exists." )
-    print()
-    '''
-    # Check if TwoPuncture initial data files exist
-    if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture"):
-        twopuncture_output = os.path.join(output_directory, "TwoPunctureABE")
-        input_par = os.path.join(output_directory, "input.par")
-
-        if os.path.exists(twopuncture_output) and os.path.exists(input_par):
-            print( " Found existing TwoPuncture initial data." )
-            print( " Do you want to skip TwoPuncture phase and reuse existing data?" )
-            print( " Input 'skip' to skip TwoPuncture and start ABE directly" )
-            print( " Input 'regenerate' to regenerate everything from scratch" )
-            print()
-            
-            while True:
-                try:
-                    inputvalue = input()
-                    if ( inputvalue == "skip" ):
-                        print( " Skipping TwoPuncture phase, will reuse existing initial data." )
-                        print()
-                        skip_twopuncture = True
-                        break
-                    elif ( inputvalue == "regenerate" ):
-                        print( " Regenerating everything from scratch." )
-                        print()
-                        skip_twopuncture = False
-                        break
-                    else:
-                        print( " Please input 'skip' or 'regenerate'." )
-                except ValueError:
-                    print( " Please input 'skip' or 'regenerate'." )
-            
-        else:
-            print( " TwoPuncture initial data not found, will regenerate everything." )
-            print()
-'''
-    # If not skipping, remove and recreate directory
-    if not skip_twopuncture:
-        shutil.rmtree(File_directory, ignore_errors=True)
-        os.mkdir(File_directory)
-        os.mkdir(output_directory)
-        os.mkdir(binary_results_directory)
-        figure_directory = os.path.join(File_directory, "figure")
-        os.mkdir(figure_directory)
-        shutil.copy("AMSS_NCKU_Input.py", File_directory)
-        print( " Output directory has been regenerated." )
-        print()
-else:
-    # Create fresh directory structure
-    os.mkdir(File_directory)
-    shutil.copy("AMSS_NCKU_Input.py", File_directory)
-    os.mkdir(output_directory)
-    os.mkdir(binary_results_directory)
-    figure_directory = os.path.join(File_directory, "figure")
-    os.mkdir(figure_directory)
-    print( " Output directory has been generated." )
-    print()
-
-# Ensure figure directory exists
-figure_directory = os.path.join(File_directory, "figure")
-if not os.path.exists(figure_directory):
-    os.mkdir(figure_directory)
-
-##################################################################
-
-## Output related parameter information
-
-import setup
-
-## Print and save input parameter information
-setup.print_input_data( File_directory )
-
-if not skip_twopuncture:
-    setup.generate_AMSSNCKU_input()
-
-setup.print_puncture_information()
-
-
-##################################################################
-
-## Generate AMSS-NCKU program input files based on the configured parameters
-
-if not skip_twopuncture:
-    print()
-    print( " Generating the AMSS-NCKU input parfile for the ABE executable." )
-    print()
-
-    ## Generate cgh-related input files from the grid information
-
-    import numerical_grid
-
-    numerical_grid.append_AMSSNCKU_cgh_input()
-
-    print()
-    print( " The input parfile for AMSS-NCKU C++ executable file ABE has been generated." )
-    print( " However, the input relevant to TwoPuncture need to be appended later." )
-    print()
-
-
-##################################################################
-
-## Plot the initial grid configuration
-
-if not skip_twopuncture:
-    print()
-    print( " Schematically plot the numerical grid structure." )
-    print()
-
-    import numerical_grid
-    numerical_grid.plot_initial_grid()
-
-
-##################################################################
-
-## Generate AMSS-NCKU macro files according to the numerical scheme and parameters
-
-if not skip_twopuncture:
-    print()
-    print( " Automatically generating the macro file for AMSS-NCKU C++ executable file ABE " )
-    print( " (Based on the finite-difference numerical scheme) " )
-    print()
-
-    import generate_macrodef
-
-    generate_macrodef.generate_macrodef_h()
-    print( " AMSS-NCKU macro file macrodef.h has been generated. " )
-
-    generate_macrodef.generate_macrodef_fh()
-    print( " AMSS-NCKU macro file macrodef.fh has been generated. " )
-
-
-##################################################################
-
-# Compile the AMSS-NCKU program according to user requirements
-# NOTE: ABE compilation is always performed, even when skipping TwoPuncture
-
-print()
-print( " Preparing to compile and run the AMSS-NCKU code as requested " )
-print( " Compiling the AMSS-NCKU code based on the generated macro files " )
-print()
-
-AMSS_NCKU_source_path = "AMSS_NCKU_source"
-AMSS_NCKU_source_copy = os.path.join(File_directory, "AMSS_NCKU_source_copy")
-
-## If AMSS_NCKU source folder is missing, create it and prompt the user
-if not os.path.exists(AMSS_NCKU_source_path):
-    os.makedirs(AMSS_NCKU_source_path)
-    print( " The AMSS-NCKU source files are incomplete; copy all source files into ./AMSS_NCKU_source. " )
-    print( " Press Enter to continue. " )
-    inputvalue = input()
-
-# Copy AMSS-NCKU source files to prepare for compilation
-# If skipping TwoPuncture and source_copy already exists, remove it first
-if skip_twopuncture and os.path.exists(AMSS_NCKU_source_copy):
-    shutil.rmtree(AMSS_NCKU_source_copy)
-
-shutil.copytree(AMSS_NCKU_source_path, AMSS_NCKU_source_copy)
-
-# Copy the generated macro files into the AMSS_NCKU source folder
-if not skip_twopuncture:
-    macrodef_h_path  = os.path.join(File_directory, "macrodef.h")
-    macrodef_fh_path = os.path.join(File_directory, "macrodef.fh")
-else:
-    # When skipping TwoPuncture, use existing macro files from previous run
-    macrodef_h_path  = os.path.join(File_directory, "macrodef.h")
-    macrodef_fh_path = os.path.join(File_directory, "macrodef.fh")
-
-shutil.copy2(macrodef_h_path,  AMSS_NCKU_source_copy)
-shutil.copy2(macrodef_fh_path, AMSS_NCKU_source_copy)
-
-# Compile related programs
-import makefile_and_run
-
-## Change working directory to the target source copy
-os.chdir(AMSS_NCKU_source_copy)
-
-## Build the main AMSS-NCKU executable (ABE or ABEGPU)
-makefile_and_run.makefile_ABE()
-
-## If the initial-data method is Ansorg-TwoPuncture, build the TwoPunctureABE executable
-## Only build TwoPunctureABE if not skipping TwoPuncture phase
-if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and not skip_twopuncture:
-    makefile_and_run.makefile_TwoPunctureABE()
-
-## Change current working directory back up two levels
-os.chdir('..')
-os.chdir('..')
-
-print()
-
-##################################################################
-
-## Copy the AMSS-NCKU executable (ABE/ABEGPU) to the run directory
-
-if (input_data.GPU_Calculation == "no"):
-    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE")
-elif (input_data.GPU_Calculation == "yes"):
-    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABEGPU")
-
-if not os.path.exists( ABE_file ):
-    print()
-    print( " Lack of AMSS-NCKU executable file ABE/ABEGPU; recompile AMSS_NCKU_source manually. " )
-    print( " When recompilation is finished, press Enter to continue. " )
-    inputvalue = input()
-
-## Copy the executable ABE (or ABEGPU) into the run directory
-shutil.copy2(ABE_file, output_directory)
-
-## If the initial-data method is TwoPuncture, copy the TwoPunctureABE executable to the run directory
-## Only copy TwoPunctureABE if not skipping TwoPuncture phase
-if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and not skip_twopuncture:
-    TwoPuncture_file = os.path.join(AMSS_NCKU_source_copy, "TwoPunctureABE")
-
-    if not os.path.exists( TwoPuncture_file ):
-        print()
-        print( " Lack of AMSS-NCKU executable file TwoPunctureABE; recompile TwoPunctureABE in AMSS_NCKU_source. " )
-        print( " When recompilation is finished, press Enter to continue. " )
-        inputvalue = input()
-
-    ## Copy the TwoPunctureABE executable into the run directory
-    shutil.copy2(TwoPuncture_file, output_directory)
-
-##################################################################
-
-## If the initial-data method is TwoPuncture, generate the TwoPuncture input files
-
-if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and not skip_twopuncture:
-
-    print()
-    print( " Initial data is chosen as Ansorg-TwoPuncture" )
-    print()
-    
-    print()
-    print( " Automatically generating the input parfile for the TwoPunctureABE executable " )
-    print()
-    
-    import generate_TwoPuncture_input
-    
-    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
-    
-    print()
-    print( " The input parfile for the TwoPunctureABE executable has been generated. " )
-    print()
-    
-    ## Generated AMSS-NCKU TwoPuncture input filename
-    AMSS_NCKU_TwoPuncture_inputfile      = 'AMSS-NCKU-TwoPuncture.input'
-    AMSS_NCKU_TwoPuncture_inputfile_path = os.path.join( File_directory, AMSS_NCKU_TwoPuncture_inputfile )
- 
-    ## Copy and rename the file
-    shutil.copy2( AMSS_NCKU_TwoPuncture_inputfile_path, os.path.join(output_directory, 'TwoPunctureinput.par') )
-    
-    ## Run TwoPuncture to generate initial-data files
-    
-    start_time = time.time()  # Record start time
-
-    print()
-    print()
-    
-    ## Change to the output (run) directory
-    os.chdir(output_directory)
-
-    ## Run the TwoPuncture executable
-    import makefile_and_run
-    makefile_and_run.run_TwoPunctureABE()
-    
-    ## Change current working directory back up two levels
-    os.chdir('..')
-    os.chdir('..')
-
-elif (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and skip_twopuncture:
-    print()
-    print( " Skipping TwoPuncture execution, using existing initial data." )
-    print()
-    start_time = time.time()  # Record start time for ABE only
-else:
-    start_time = time.time()  # Record start time
-    
-##################################################################
-    
-## Update puncture data based on TwoPuncture run results
-
-if not skip_twopuncture:
-    import renew_puncture_parameter
-    renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory)
-
-    ## Generated AMSS-NCKU input filename
-    AMSS_NCKU_inputfile      = 'AMSS-NCKU.input'
-    AMSS_NCKU_inputfile_path = os.path.join(File_directory, AMSS_NCKU_inputfile)
- 
-    ## Copy and rename the file
-    shutil.copy2( AMSS_NCKU_inputfile_path, os.path.join(output_directory, 'input.par') )
-
-    print()
-    print( " Successfully copy all AMSS-NCKU input parfile to target dictionary. " )
-    print()
-else:
-    print()
-    print( " Using existing input.par file from previous run." )
-    print()
-
-##################################################################
-
-## Launch the AMSS-NCKU program
-
-print()
-print()
-
-## Change to the run directory
-os.chdir( output_directory )
-
-import makefile_and_run
-makefile_and_run.run_ABE()
-
-## Change current working directory back up two levels
-os.chdir('..')
-os.chdir('..')
-
-end_time = time.time()
-elapsed_time = end_time - start_time
-
-##################################################################
-
-## Copy some basic input and log files out to facilitate debugging
-
-## Path to the file that stores calculation settings
-AMSS_NCKU_error_file_path = os.path.join(binary_results_directory, "setting.par")
-## Copy and rename the file for easier inspection
-shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directory, "AMSSNCKU_setting_parameter") )
-
-## Path to the error log file
-AMSS_NCKU_error_file_path = os.path.join(binary_results_directory, "Error.log")
-## Copy and rename the error log
-shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directory, "Error.log") )
-
-## Primary program outputs
-AMSS_NCKU_BH_data         = os.path.join(binary_results_directory, "bssn_BH.dat"        )
-AMSS_NCKU_ADM_data        = os.path.join(binary_results_directory, "bssn_ADMQs.dat"     )
-AMSS_NCKU_psi4_data       = os.path.join(binary_results_directory, "bssn_psi4.dat"      )
-AMSS_NCKU_constraint_data = os.path.join(binary_results_directory, "bssn_constraint.dat")
-## copy and rename the file
-shutil.copy( AMSS_NCKU_BH_data,         os.path.join(output_directory, "bssn_BH.dat"        ) )
-shutil.copy( AMSS_NCKU_ADM_data,        os.path.join(output_directory, "bssn_ADMQs.dat"     ) )
-shutil.copy( AMSS_NCKU_psi4_data,       os.path.join(output_directory, "bssn_psi4.dat"      ) )
-shutil.copy( AMSS_NCKU_constraint_data, os.path.join(output_directory, "bssn_constraint.dat") )
-
-## Additional program outputs
-if (input_data.Equation_Class == "BSSN-EM"):
-    AMSS_NCKU_phi1_data = os.path.join(binary_results_directory, "bssn_phi1.dat" )
-    AMSS_NCKU_phi2_data = os.path.join(binary_results_directory, "bssn_phi2.dat" )
-    shutil.copy( AMSS_NCKU_phi1_data, os.path.join(output_directory, "bssn_phi1.dat" ) )
-    shutil.copy( AMSS_NCKU_phi2_data, os.path.join(output_directory, "bssn_phi2.dat" ) )
-elif (input_data.Equation_Class == "BSSN-EScalar"):
-    AMSS_NCKU_maxs_data = os.path.join(binary_results_directory, "bssn_maxs.dat" )
-    shutil.copy( AMSS_NCKU_maxs_data, os.path.join(output_directory, "bssn_maxs.dat" ) )
-
-##################################################################
-
-## Plot the AMSS-NCKU program results
-
-print()
-print( " Plotting the txt and binary results data from the AMSS-NCKU simulation " )
-print()
-
-
-import plot_xiaoqu
-import plot_GW_strain_amplitude_xiaoqu
-
-## Plot black hole trajectory
-plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directory, figure_directory )
-plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directory, figure_directory )
-
-## Plot black hole separation vs. time
-plot_xiaoqu.generate_puncture_distence_plot( binary_results_directory, figure_directory )
-
-## Plot gravitational waveforms (psi4 and strain amplitude)
-for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_gravitational_wave_psi4_plot( binary_results_directory, figure_directory, i )
-    plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot( binary_results_directory, figure_directory, i )
-
-## Plot ADM mass evolution
-for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_ADMmass_plot( binary_results_directory, figure_directory, i )
-
-## Plot Hamiltonian constraint violation over time
-for i in range(input_data.grid_level):
-    plot_xiaoqu.generate_constraint_check_plot( binary_results_directory, figure_directory, i )
-
-## Plot stored binary data
-plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )
-
-print()
-print( f" This Program Cost = {elapsed_time} Seconds " )
-print()
-
-
-##################################################################
-
-print()
-print( " The AMSS-NCKU-Python simulation is successfully finished, thanks for using !!! " )
-print()
-
-##################################################################
-
-

AMSS_NCKU_Program.py

@@ -8,6 +8,14 @@
 ##
 ##################################################################
 
+## Guard against re-execution by multiprocessing child processes.
+## Without this, using 'spawn' or 'forkserver' context would cause every
+## worker to re-run the entire script, spawning exponentially more
+## workers (fork bomb).
+if __name__ != '__main__':
+    import sys as _sys
+    _sys.exit(0)
+
 
 ##################################################################
 
@@ -118,11 +126,6 @@ setup.generate_AMSSNCKU_input()
 #inputvalue = input()  ## Wait for user input (press Enter) to proceed
 #print()
 
-setup.print_puncture_information()
-
-
-##################################################################
-
 ## Generate AMSS-NCKU program input files based on the configured parameters
 
 print(                                                                                   )
@@ -262,6 +265,12 @@ if not os.path.exists( ABE_file ):
 ## Copy the executable ABE (or ABEGPU) into the run directory
 shutil.copy2(ABE_file, output_directory)
 
+## Copy interp load balance profile if present (for optimize pass)
+interp_lb_profile = os.path.join(AMSS_NCKU_source_copy, "interp_lb_profile.bin")
+if os.path.exists(interp_lb_profile):
+    shutil.copy2(interp_lb_profile, output_directory)
+    print( " Copied interp_lb_profile.bin to run directory " )
+
 ###########################
 
 ## If the initial-data method is TwoPuncture, copy the TwoPunctureABE executable to the run directory
@@ -298,7 +307,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
     
     import generate_TwoPuncture_input
     
-    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
+    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input(numerical_grid.puncture_data)
     
     print(                                                                                              )
     print( " The input parfile for the TwoPunctureABE executable has been generated. " )
@@ -340,7 +349,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
 
 import renew_puncture_parameter
     
-renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory)
+renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory, numerical_grid.puncture_data)
 
 
 ## Generated AMSS-NCKU input filename
@@ -424,26 +433,31 @@ print(
 
 import plot_xiaoqu
 import plot_GW_strain_amplitude_xiaoqu
+from parallel_plot_helper import run_plot_tasks_parallel
+
+plot_tasks = []
 
 ## Plot black hole trajectory
-plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot,   (binary_results_directory, figure_directory) ) )
-plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot3D, (binary_results_directory, figure_directory) ) )
 
 ## Plot black hole separation vs. time
-plot_xiaoqu.generate_puncture_distence_plot( binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_distence_plot, (binary_results_directory, figure_directory) ) )
 
 ## Plot gravitational waveforms (psi4 and strain amplitude)
 for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_gravitational_wave_psi4_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_gravitational_wave_psi4_plot, (binary_results_directory, figure_directory, i) ) )
-    plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot, (binary_results_directory, figure_directory, i) ) )
 
 ## Plot ADM mass evolution
 for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_ADMmass_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_ADMmass_plot, (binary_results_directory, figure_directory, i) ) )
 
 ## Plot Hamiltonian constraint violation over time
 for i in range(input_data.grid_level):
-    plot_xiaoqu.generate_constraint_check_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_constraint_check_plot, (binary_results_directory, figure_directory, i) ) )
+
+run_plot_tasks_parallel(plot_tasks)
 
 ## Plot stored binary data
 plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )

AMSS_NCKU_Verify_ASC26.py

@@ -1,10 +1,19 @@
 #!/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)
+3. The following figure PDFs must match GW150914-origin exactly after rasterization:
+   - ADM_Constraint_Grid_Level_0.pdf
+   - BH_Trajectory_21_XY.pdf
+   - BH_Trajectory_XY.pdf
+   The script also reports the percentage of differing pixels for each figure.
 
 RMS Calculation Method:
 - Computes trajectory deviation on the XY plane independently for BH1 and BH2
@@ -19,6 +28,10 @@ Reference: GW150914-origin (baseline simulation)
 import numpy as np
 import sys
 import os
+import shutil
+import subprocess
+import tempfile
+from PIL import Image
 
 # ANSI Color Codes
 class Color:
@@ -58,78 +71,187 @@ def load_constraint_data(filepath):
     return np.array(data)
 
 
-def calculate_rms_error(bh_data_ref, bh_data_target):
+def resolve_figure_dir(path):
+    """Resolve the sibling figure directory from an output or figure path."""
+    normalized = os.path.normpath(path)
+    if os.path.basename(normalized) == "figure":
+        return normalized
+    return os.path.join(os.path.dirname(normalized), "figure")
+
+
+def render_pdf_to_images(pdf_path, dpi=150):
+    """Render a PDF to RGB images using Ghostscript."""
+    gs_path = shutil.which("gs")
+    if gs_path is None:
+        raise RuntimeError("Ghostscript executable 'gs' was not found in PATH")
+
+    with tempfile.TemporaryDirectory(prefix="amss_verify_pdf_") as temp_dir:
+        output_pattern = os.path.join(temp_dir, "page-%03d.ppm")
+        cmd = [
+            gs_path,
+            "-q",
+            "-dSAFER",
+            "-dBATCH",
+            "-dNOPAUSE",
+            "-sDEVICE=ppmraw",
+            f"-r{dpi}",
+            f"-o{output_pattern}",
+            pdf_path
+        ]
+
+        try:
+            subprocess.run(cmd, check=True, stdout=subprocess.DEVNULL, stderr=subprocess.PIPE, text=True)
+        except subprocess.CalledProcessError as exc:
+            message = exc.stderr.strip() or str(exc)
+            raise RuntimeError(f"Failed to render PDF '{pdf_path}': {message}") from exc
+
+        ppm_files = sorted(
+            os.path.join(temp_dir, filename)
+            for filename in os.listdir(temp_dir)
+            if filename.endswith(".ppm")
+        )
+
+        if not ppm_files:
+            raise RuntimeError(f"No rendered pages were produced for '{pdf_path}'")
+
+        images = []
+        for ppm_file in ppm_files:
+            with Image.open(ppm_file) as img:
+                images.append(np.array(img.convert("RGB"), dtype=np.uint8))
+
+        return images
+
+
+def compare_rendered_pages(ref_img, target_img):
+    """Return (different_pixels, total_pixels) for two rendered RGB pages."""
+    ref_h, ref_w = ref_img.shape[:2]
+    tgt_h, tgt_w = target_img.shape[:2]
+    total_pixels = max(ref_h, tgt_h) * max(ref_w, tgt_w)
+
+    if ref_h == tgt_h and ref_w == tgt_w:
+        different_pixels = int(np.count_nonzero(np.any(ref_img != target_img, axis=2)))
+        return different_pixels, total_pixels
+
+    diff_mask = np.ones((max(ref_h, tgt_h), max(ref_w, tgt_w)), dtype=bool)
+    overlap_h = min(ref_h, tgt_h)
+    overlap_w = min(ref_w, tgt_w)
+    overlap_diff = np.any(ref_img[:overlap_h, :overlap_w] != target_img[:overlap_h, :overlap_w], axis=2)
+    diff_mask[:overlap_h, :overlap_w] = overlap_diff
+    different_pixels = int(np.count_nonzero(diff_mask))
+    return different_pixels, total_pixels
+
+
+def compare_pdf_images(ref_pdf, target_pdf, dpi=150, threshold_percent=0.001):
+    """Compare two PDFs by rasterizing them and counting differing pixels."""
+    ref_pages = render_pdf_to_images(ref_pdf, dpi=dpi)
+    target_pages = render_pdf_to_images(target_pdf, dpi=dpi)
+
+    total_pixels = 0
+    different_pixels = 0
+    max_pages = max(len(ref_pages), len(target_pages))
+
+    for page_idx in range(max_pages):
+        if page_idx < len(ref_pages) and page_idx < len(target_pages):
+            page_diff, page_total = compare_rendered_pages(ref_pages[page_idx], target_pages[page_idx])
+        else:
+            existing_page = ref_pages[page_idx] if page_idx < len(ref_pages) else target_pages[page_idx]
+            page_total = existing_page.shape[0] * existing_page.shape[1]
+            page_diff = page_total
+
+        total_pixels += page_total
+        different_pixels += page_diff
+
+    diff_percent = (different_pixels / total_pixels * 100.0) if total_pixels else 0.0
+    return {
+        "different_pixels": different_pixels,
+        "total_pixels": total_pixels,
+        "diff_percent": diff_percent,
+        "pages_ref": len(ref_pages),
+        "pages_target": len(target_pages),
+        "passed": diff_percent < threshold_percent
+    }
+
+
+def compare_required_figures(reference_figure_dir, target_figure_dir):
+    """Compare the required GW150914 figure PDFs."""
+    figure_names = [
+        "ADM_Constraint_Grid_Level_0.pdf",
+        "BH_Trajectory_21_XY.pdf",
+        "BH_Trajectory_XY.pdf"
+    ]
+
+    results = []
+    for figure_name in figure_names:
+        ref_pdf = os.path.join(reference_figure_dir, figure_name)
+        target_pdf = os.path.join(target_figure_dir, figure_name)
+
+        if not os.path.exists(ref_pdf):
+            raise FileNotFoundError(f"Reference figure not found: {ref_pdf}")
+        if not os.path.exists(target_pdf):
+            raise FileNotFoundError(f"Target figure not found: {target_pdf}")
+
+        comparison = compare_pdf_images(ref_pdf, target_pdf)
+        comparison["name"] = figure_name
+        results.append(comparison)
+
+    return results
+
+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.
+    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 +277,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):
@@ -199,19 +319,45 @@ def print_constraint_results(results, threshold=2.0):
     return passed
 
 
-def print_summary(rms_passed, constraint_passed):
+def print_figure_results(results, threshold_percent=0.001):
-    """Print summary"""
+    print(f"\n{Color.BOLD}3. Figure Pixel Comparison (PDF Rasterization){Color.RESET}")
+    print("-" * 65)
+    print(f"   Requirement: < {threshold_percent:.3f}% differing pixels\n")
+
+    all_passed = True
+    for result in results:
+        passed = result["passed"]
+        all_passed = all_passed and passed
+        status = get_status_text(passed)
+        print(f"   {result['name']:32}: {result['diff_percent']:10.6f}%   |   Status: {status}")
+
+        if result["pages_ref"] != result["pages_target"]:
+            print(f"   {'':32}  pages(ref/target): {result['pages_ref']}/{result['pages_target']}")
+
+    return all_passed
+
+
+def print_figure_error(error_message):
+    print(f"\n{Color.BOLD}3. Figure Pixel Comparison (PDF Rasterization){Color.RESET}")
+    print("-" * 65)
+    print(f"   {Color.RED}Error: {error_message}{Color.RESET}")
+    return False
+
+
+def print_summary(rms_passed, constraint_passed, figure_passed):
     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
+    all_passed = rms_passed and constraint_passed and figure_passed
     
     res_rms = get_status_text(rms_passed)
     res_con = get_status_text(constraint_passed)
+    res_fig = get_status_text(figure_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}")
+    print(f"   [3] Figure pixel comparison:      {res_fig}")
     
     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}")
@@ -219,61 +365,58 @@ 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")
+    target_figure_dir = resolve_figure_dir(target_dir)
+    reference_figure_dir = os.path.join(script_dir, "GW150914-origin/figure")
 
-    # 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}")
+    print(f"{Color.BOLD}Reference Figures:   {Color.RESET} {Color.BLUE}{reference_figure_dir}{Color.RESET}")
+    print(f"{Color.BOLD}Target Figures:      {Color.RESET} {Color.BLUE}{target_figure_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
+    try:
-    all_passed = print_summary(rms_passed, constraint_passed)
+        figure_results = compare_required_figures(reference_figure_dir, target_figure_dir)
+        figure_passed = print_figure_results(figure_results)
+    except (FileNotFoundError, RuntimeError) as exc:
+        figure_passed = print_figure_error(str(exc))
 
-    # Return exit code
+    all_passed = print_summary(rms_passed, constraint_passed, figure_passed)
     sys.exit(0 if all_passed else 1)
 
-
 if __name__ == "__main__":
     main()

AMSS_NCKU_source/BSSN/bssn_class.h

@@ -48,6 +48,7 @@ public:
        double StartTime, TotalTime;
        double AnasTime, DumpTime, d2DumpTime, CheckTime;
        double LastAnas, LastConsOut;
+       int *ConstraintRefreshLevels;
        double Courant;
        double numepss, numepsb, numepsh;
        int Symmetry;
@@ -126,8 +127,15 @@ public:
        MyList<var> *OldStateList, *DumpList;
        MyList<var> *ConstraintList;
 
+       Parallel::SyncCache *sync_cache_pre;  // per-level cache for predictor sync
+       Parallel::SyncCache *sync_cache_cor;  // per-level cache for corrector sync
+       Parallel::SyncCache *sync_cache_rp_coarse;  // RestrictProlong sync on PatL[lev-1]
+       Parallel::SyncCache *sync_cache_rp_fine;    // RestrictProlong sync on PatL[lev]
+       Parallel::SyncCache *sync_cache_restrict;   // cached Restrict in RestrictProlong
+       Parallel::SyncCache *sync_cache_outbd;      // cached OutBdLow2Hi in RestrictProlong
+
        monitor *ErrorMonitor, *Psi4Monitor, *BHMonitor, *MAPMonitor;
-       monitor *ConVMonitor;
+       monitor *ConVMonitor, *TimingMonitor;
        surface_integral *Waveshell;
        checkpoint *CheckPoint;
 

AMSS_NCKU_source/BSSN/bssn_rhs.h

@@ -32,6 +32,19 @@
 #define f_compute_rhs_Z4c_ss compute_rhs_z4c_ss_
 #define f_compute_constraint_fr compute_constraint_fr_
 #endif
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+        void f_bssn_rhs_kernel_timing_reset();
+        int f_bssn_rhs_kernel_timing_bucket_count();
+        const double *f_bssn_rhs_kernel_timing_local_seconds();
+        const char *f_bssn_rhs_kernel_timing_label(int);
+#ifdef __cplusplus
+}
+#endif
+
 extern "C"
 {
         int f_compute_rhs_bssn(int *, double &, double *, double *, double *,                                                      // ex,T,X,Y,Z

AMSS_NCKU_source/BSSN/lopsided_c.C

@@ -0,0 +1,255 @@
+#include "tool.h"
+/*
+ * 你需要提供 symmetry_bd 的 C 版本（或 Fortran 绑到 C 的接口）。
+ * Fortran: call symmetry_bd(3,ex,f,fh,SoA)
+ *
+ * 约定：
+ *   nghost = 3
+ *   ex[3]  = {ex1,ex2,ex3}
+ *   f      = 原始网格 (ex1*ex2*ex3)
+ *   fh     = 扩展网格 ((ex1+3)*(ex2+3)*(ex3+3))，对应 Fortran 的 (-2:ex1, ...)
+ *   SoA[3] = 输入参数
+ */
+void lopsided(const int ex[3],
+              const double *X, const double *Y, const double *Z,
+              const double *f, double *f_rhs,
+              const double *Sfx, const double *Sfy, const double *Sfz,
+              int Symmetry, const double SoA[3])
+{
+    const double ZEO = 0.0, ONE = 1.0, F3 = 3.0;
+    const double TWO = 2.0, F6 = 6.0, F18 = 18.0;
+    const double F12 = 12.0, F10 = 10.0, EIT = 8.0;
+
+    const int NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2;
+    (void)OCTANT; // 这里和 Fortran 一样只是定义了不用也没关系
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    // 对应 Fortran: dX = X(2)-X(1)  （Fortran 1-based）
+    // C: X[1]-X[0]
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+
+    const double d12dx = ONE / F12 / dX;
+    const double d12dy = ONE / F12 / dY;
+    const double d12dz = ONE / F12 / dZ;
+
+    // Fortran 里算了 d2dx/d2dy/d2dz 但本 subroutine 里没用到（保持一致也算出来）
+    const double d2dx  = ONE / TWO / dX;
+    const double d2dy  = ONE / TWO / dY;
+    const double d2dz  = ONE / TWO / dZ;
+    (void)d2dx; (void)d2dy; (void)d2dz;
+
+    // Fortran:
+    // imax = ex(1); jmax = ex(2); kmax = ex(3)
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    // Fortran:
+    // imin=jmin=kmin=1; 若满足对称条件则设为 -2
+    int iminF = 1, jminF = 1, kminF = 1;
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -2;
+
+    // 分配 fh：大小 (ex1+3)*(ex2+3)*(ex3+3)
+    const size_t nx = (size_t)ex1 + 3;
+    const size_t ny = (size_t)ex2 + 3;
+    const size_t nz = (size_t)ex3 + 3;
+    const size_t fh_size = nx * ny * nz;
+
+    double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return; // 内存不足：直接返回（你也可以改成 abort/报错）
+
+    // Fortran: call symmetry_bd(3,ex,f,fh,SoA)
+    symmetry_bd(3, ex, f, fh, SoA);
+
+    /*
+     * Fortran 主循环：
+     * do k=1,ex(3)-1
+     * do j=1,ex(2)-1
+     * do i=1,ex(1)-1
+     *
+     * 转成 C 0-based：
+     * k0 = 0..ex3-2, j0 = 0..ex2-2, i0 = 0..ex1-2
+     *
+     * 并且 Fortran 里的 i/j/k 在 fh 访问时，仍然是 Fortran 索引值：
+     * iF=i0+1, jF=j0+1, kF=k0+1
+     */
+    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
+                const int iF = i0 + 1;
+
+                const size_t p = idx_ex(i0, j0, k0, ex);
+
+                // ---------------- x direction ----------------
+                const double sfx = Sfx[p];
+                if (sfx > ZEO) {
+                    // Fortran: if(i+3 <= imax)
+                    // iF+3 <= ex1  <=> i0+4 <= ex1 <=> i0 <= ex1-4
+                    if (i0 <= ex1 - 4) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                    // elseif(i+2 <= imax)  <=> i0 <= ex1-3
+                    else if (i0 <= ex1 - 3) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    }
+                    // elseif(i+1 <= imax)  <=> i0 <= ex1-2（循环里总成立）
+                    else if (i0 <= ex1 - 2) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                } else if (sfx < ZEO) {
+                    // Fortran: if(i-3 >= imin)
+                    // (iF-3) >= iminF  <=> (i0-2) >= iminF
+                    if ((i0 - 2) >= iminF) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                    // elseif(i-2 >= imin) <=> (i0-1) >= iminF
+                    else if ((i0 - 1) >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    }
+                    // elseif(i-1 >= imin) <=> i0 >= iminF
+                    else if (i0 >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                }
+
+                // ---------------- y direction ----------------
+                const double sfy = Sfy[p];
+                if (sfy > ZEO) {
+                    // jF+3 <= ex2 <=> j0+4 <= ex2 <=> j0 <= ex2-4
+                    if (j0 <= ex2 - 4) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    } else if (j0 <= ex2 - 3) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 <= ex2 - 2) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    }
+                } else if (sfy < ZEO) {
+                    if ((j0 - 2) >= jminF) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    } else if ((j0 - 1) >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    }
+                }
+
+                // ---------------- z direction ----------------
+                const double sfz = Sfz[p];
+                if (sfz > ZEO) {
+                    if (k0 <= ex3 - 4) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    } else if (k0 <= ex3 - 3) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 <= ex3 - 2) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    }
+                } else if (sfz < ZEO) {
+                    if ((k0 - 2) >= kminF) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    } else if ((k0 - 1) >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    }
+                }
+            }
+        }
+    }
+    free(fh);
+}
+
+
+
+
+

AMSS_NCKU_source/BSSN/lopsided_kodis_c.C

@@ -0,0 +1,248 @@
+#include "tool.h"
+
+/*
+ * Combined advection (lopsided) + KO dissipation (kodis).
+ * Uses one shared symmetry_bd buffer per call.
+ */
+void lopsided_kodis(const int ex[3],
+                    const double *X, const double *Y, const double *Z,
+                    const double *f, double *f_rhs,
+                    const double *Sfx, const double *Sfy, const double *Sfz,
+                    int Symmetry, const double SoA[3], double eps)
+{
+    const double ZEO = 0.0, ONE = 1.0, F3 = 3.0;
+    const double F6 = 6.0, F18 = 18.0;
+    const double F12 = 12.0, F10 = 10.0, EIT = 8.0;
+    const double SIX = 6.0, FIT = 15.0, TWT = 20.0;
+    const double cof = 64.0; // 2^6
+
+    const int NO_SYMM = 0, EQ_SYMM = 1;
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+
+    const double d12dx = ONE / F12 / dX;
+    const double d12dy = ONE / F12 / dY;
+    const double d12dz = ONE / F12 / dZ;
+
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    int iminF = 1, jminF = 1, kminF = 1;
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -2;
+
+    // fh for Fortran-style domain (-2:ex1,-2:ex2,-2:ex3)
+    const size_t nx = (size_t)ex1 + 3;
+    const size_t ny = (size_t)ex2 + 3;
+    const size_t nz = (size_t)ex3 + 3;
+    const size_t fh_size = nx * ny * nz;
+
+    double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return;
+
+    symmetry_bd(3, ex, f, fh, SoA);
+
+    // Advection (same stencil logic as lopsided_c.C)
+    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
+                const int iF = i0 + 1;
+                const size_t p = idx_ex(i0, j0, k0, ex);
+
+                const double sfx = Sfx[p];
+                if (sfx > ZEO) {
+                    if (i0 <= ex1 - 4) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    } else if (i0 <= ex1 - 3) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    } else if (i0 <= ex1 - 2) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                } else if (sfx < ZEO) {
+                    if ((i0 - 2) >= iminF) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    } else if ((i0 - 1) >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    } else if (i0 >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                }
+
+                const double sfy = Sfy[p];
+                if (sfy > ZEO) {
+                    if (j0 <= ex2 - 4) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    } else if (j0 <= ex2 - 3) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 <= ex2 - 2) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    }
+                } else if (sfy < ZEO) {
+                    if ((j0 - 2) >= jminF) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    } else if ((j0 - 1) >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    }
+                }
+
+                const double sfz = Sfz[p];
+                if (sfz > ZEO) {
+                    if (k0 <= ex3 - 4) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    } else if (k0 <= ex3 - 3) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 <= ex3 - 2) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    }
+                } else if (sfz < ZEO) {
+                    if ((k0 - 2) >= kminF) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    } else if ((k0 - 1) >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    }
+                }
+            }
+        }
+    }
+
+    // KO dissipation (same domain restriction as kodiss_c.C)
+    if (eps > ZEO) {
+        const int i0_lo = (iminF + 2 > 0) ? iminF + 2 : 0;
+        const int j0_lo = (jminF + 2 > 0) ? jminF + 2 : 0;
+        const int k0_lo = (kminF + 2 > 0) ? kminF + 2 : 0;
+        const int i0_hi = imaxF - 4; // inclusive
+        const int j0_hi = jmaxF - 4;
+        const int k0_hi = kmaxF - 4;
+
+        if (!(i0_lo > i0_hi || j0_lo > j0_hi || k0_lo > k0_hi)) {
+            for (int k0 = k0_lo; k0 <= k0_hi; ++k0) {
+                const int kF = k0 + 1;
+                for (int j0 = j0_lo; j0 <= j0_hi; ++j0) {
+                    const int jF = j0 + 1;
+                    for (int i0 = i0_lo; i0 <= i0_hi; ++i0) {
+                        const int iF = i0 + 1;
+                        const size_t p = idx_ex(i0, j0, k0, ex);
+
+                        const double Dx_term =
+                            ((fh[idx_fh_F(iF - 3, jF, kF, ex)] + fh[idx_fh_F(iF + 3, jF, kF, ex)]) -
+                             SIX * (fh[idx_fh_F(iF - 2, jF, kF, ex)] + fh[idx_fh_F(iF + 2, jF, kF, ex)]) +
+                             FIT * (fh[idx_fh_F(iF - 1, jF, kF, ex)] + fh[idx_fh_F(iF + 1, jF, kF, ex)]) -
+                             TWT *  fh[idx_fh_F(iF,     jF, kF, ex)]) / dX;
+
+                        const double Dy_term =
+                            ((fh[idx_fh_F(iF, jF - 3, kF, ex)] + fh[idx_fh_F(iF, jF + 3, kF, ex)]) -
+                             SIX * (fh[idx_fh_F(iF, jF - 2, kF, ex)] + fh[idx_fh_F(iF, jF + 2, kF, ex)]) +
+                             FIT * (fh[idx_fh_F(iF, jF - 1, kF, ex)] + fh[idx_fh_F(iF, jF + 1, kF, ex)]) -
+                             TWT *  fh[idx_fh_F(iF, jF,     kF, ex)]) / dY;
+
+                        const double Dz_term =
+                            ((fh[idx_fh_F(iF, jF, kF - 3, ex)] + fh[idx_fh_F(iF, jF, kF + 3, ex)]) -
+                             SIX * (fh[idx_fh_F(iF, jF, kF - 2, ex)] + fh[idx_fh_F(iF, jF, kF + 2, ex)]) +
+                             FIT * (fh[idx_fh_F(iF, jF, kF - 1, ex)] + fh[idx_fh_F(iF, jF, kF + 1, ex)]) -
+                             TWT *  fh[idx_fh_F(iF, jF, kF,     ex)]) / dZ;
+
+                        f_rhs[p] += (eps / cof) * (Dx_term + Dy_term + Dz_term);
+                    }
+                }
+            }
+        }
+    }
+
+    free(fh);
+}

AMSS_NCKU_source/BSSN/prolongrestrict_cell.f90

@@ -1934,18 +1934,35 @@
 ! when if=1 -> ic=0, this is different to vertex center grid 
   real*8, dimension(-2:extc(1),-2:extc(2),-2:extc(3))   :: funcc
   integer,dimension(3) :: cxI
-  integer :: i,j,k,ii,jj,kk
+  integer :: i,j,k,ii,jj,kk,px,py,pz
   real*8, dimension(6,6) :: tmp2
   real*8, dimension(6) :: tmp1
+  integer, dimension(extf(1)) :: cix
+  integer, dimension(extf(2)) :: ciy
+  integer, dimension(extf(3)) :: ciz
+  integer, dimension(extf(1)) :: pix
+  integer, dimension(extf(2)) :: piy
+  integer, dimension(extf(3)) :: piz
 
   real*8, parameter :: C1=7.7d1/8.192d3,C2=-6.93d2/8.192d3,C3=3.465d3/4.096d3
   real*8, parameter :: C6=6.3d1/8.192d3,C5=-4.95d2/8.192d3,C4=1.155d3/4.096d3
+  real*8, dimension(6,2), parameter :: WC = reshape((/&
+      C1,C2,C3,C4,C5,C6,&
+      C6,C5,C4,C3,C2,C1/), (/6,2/))
 
   integer::imini,imaxi,jmini,jmaxi,kmini,kmaxi
   integer::imino,imaxo,jmino,jmaxo,kmino,kmaxo
+  integer::maxcx,maxcy,maxcz
 
   real*8,dimension(3) :: CD,FD
-  
+  real*8 :: tmp_yz(extc(1), 6)      ! 存储整条 X 线上 6 个 Y 轴偏置的 Z 向插值结果
+  real*8 :: tmp_xyz_line(-2:extc(1))   ! 包含 X 向 6 点模板访问所需下界
+  real*8 :: v1, v2, v3, v4, v5, v6
+  integer :: ic, jc, kc, ix_offset,ix,iy,iz,jc_min,jc_max,ic_min,ic_max,kc_min,kc_max
+  integer :: i_lo, i_hi, j_lo, j_hi, k_lo, k_hi
+  logical :: need_full_symmetry
+  real*8 :: res_line
+  real*8 :: tmp_z_slab(-2:extc(1), -2:extc(2))  ! 包含 Y/X 向模板访问所需下界
   if(wei.ne.3)then
      write(*,*)"prolongrestrict.f90::prolong3: this routine only surport 3 dimension"
      write(*,*)"dim = ",wei
@@ -2020,145 +2037,140 @@
           return
   endif
 
-  call symmetry_bd(3,extc,func,funcc,SoA)
+  do i = imino,imaxo
-     
+     ii = i + lbf(1) - 1
-!~~~~~~> prolongation start...
+     cix(i) = ii/2 - lbc(1) + 1
+     if(ii/2*2 == ii)then
+        pix(i) = 1
+     else
+        pix(i) = 2
+     endif
+  enddo
+  do j = jmino,jmaxo
+     jj = j + lbf(2) - 1
+     ciy(j) = jj/2 - lbc(2) + 1
+     if(jj/2*2 == jj)then
+        piy(j) = 1
+     else
+        piy(j) = 2
+     endif
+  enddo
   do k = kmino,kmaxo
-   do j = jmino,jmaxo
+     kk = k + lbf(3) - 1
-    do i = imino,imaxo
+     ciz(k) = kk/2 - lbc(3) + 1
-       cxI(1) = i
+     if(kk/2*2 == kk)then
-       cxI(2) = j
+        piz(k) = 1
-       cxI(3) = k
+     else
-! change to coarse level reference
+        piz(k) = 2
-!|---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*---| 
+     endif
-!|=======x===============x===============x===============x=======|
-       cxI = (cxI+lbf-1)/2
-! change to array index      
-       cxI = cxI - lbc + 1
-
-       if(any(cxI+3 > extc)) write(*,*)"error in prolong"
-       ii=i+lbf(1)-1
-       jj=j+lbf(2)-1
-       kk=k+lbf(3)-1
-#if 0
-       if(ii/2*2==ii)then
-         if(jj/2*2==jj)then
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)=  C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           endif
-         else
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)=  C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           endif
-         endif
-       else
-         if(jj/2*2==jj)then
-           if(kk/2*2==kk)then               
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)=  C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           endif
-         else
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)=  C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           endif
-         endif
-       endif
-#else 
-       if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-       else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-       endif
-
-       if(jj/2*2==jj)then
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-       else
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-       endif
-
-       if(ii/2*2==ii)then
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-       else
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-       endif
-#endif
-    enddo
-   enddo
   enddo
 
+  ic_min = minval(cix(imino:imaxo))
+  ic_max = maxval(cix(imino:imaxo))
+  jc_min = minval(ciy(jmino:jmaxo))
+  jc_max = maxval(ciy(jmino:jmaxo))
+  kc_min = minval(ciz(kmino:kmaxo))
+  kc_max = maxval(ciz(kmino:kmaxo))
+
+  maxcx = ic_max
+  maxcy = jc_max
+  maxcz = kc_max
+  if(maxcx+3 > extc(1) .or. maxcy+3 > extc(2) .or. maxcz+3 > extc(3))then
+     write(*,*)"error in prolong"
+     return
+  endif
+
+  i_lo = ic_min - 2
+  i_hi = ic_max + 3
+  j_lo = jc_min - 2
+  j_hi = jc_max + 3
+  k_lo = kc_min - 2
+  k_hi = kc_max + 3
+  need_full_symmetry = (i_lo < 1) .or. (j_lo < 1) .or. (k_lo < 1)
+  if(need_full_symmetry)then
+     call symmetry_bd(3,extc,func,funcc,SoA)
+  else
+     funcc(i_lo:i_hi,j_lo:j_hi,k_lo:k_hi) = func(i_lo:i_hi,j_lo:j_hi,k_lo:k_hi)
+  endif
+
+     ! 对每个 k（pz, kc 固定）预计算 Z 向插值的 2D 切片
+
+do k = kmino, kmaxo
+    pz = piz(k); kc = ciz(k)
+    ! --- Pass 1: Z 方向，只算一次 ---
+    do iy = jc_min-2, jc_max+3   ! 仅需的 iy 范围（对应 jc-2:jc+3）
+        do ii = ic_min-2, ic_max+3  ! 仅需的 ii 范围（对应 cix-2:cix+3）
+            tmp_z_slab(ii, iy) = sum(WC(:,pz) * funcc(ii, iy, kc-2:kc+3))
+        end do
+    end do
+
+    do j = jmino, jmaxo
+        py = piy(j); jc = ciy(j)
+        ! --- Pass 2: Y 方向 ---
+        do ii = ic_min-2, ic_max+3
+            tmp_xyz_line(ii) = sum(WC(:,py) * tmp_z_slab(ii, jc-2:jc+3))
+        end do
+        ! --- Pass 3: X 方向 ---
+        do i = imino, imaxo
+            funf(i,j,k) = sum(WC(:,pix(i)) * tmp_xyz_line(cix(i)-2:cix(i)+3))
+        end do
+    end do
+end do
+
+!~~~~~~> prolongation start...
+#if 0
+ do k = kmino, kmaxo
+     pz = piz(k)
+     kc = ciz(k)
+
+     do j = jmino, jmaxo
+        py = piy(j)
+        jc = ciy(j)
+
+! --- 步骤 1 & 2 融合：分段处理 X 轴，提升 Cache 命中率 ---
+        ! 我们将 ii 循环逻辑重组，减少对 funcc 的跨行重复访问
+        do ii = 1, extc(1)
+           ! 1. 先做 Z 方向的 6 条线插值（针对当前的 ii 和当前的 6 个 iy）
+           ! 我们直接在这里把 Y 方向的加权也做了，省去 tmp_yz 数组
+           ! 这样 funcc 的数据读进来后立即完成所有维度的贡献，不再写回内存
+           
+           res_line = 0.0d0
+           do jj = 1, 6
+              iy = jc - 3 + jj
+              ! 这一行代码是核心：一次性完成 Z 插值并加上 Y 的权重
+              ! 编译器会把 WC(jj, py) 存在寄存器里
+              res_line = res_line + WC(jj, py) * ( &
+                         WC(1, pz) * funcc(ii, iy, kc-2) + &
+                         WC(2, pz) * funcc(ii, iy, kc-1) + &
+                         WC(3, pz) * funcc(ii, iy, kc  ) + &
+                         WC(4, pz) * funcc(ii, iy, kc+1) + &
+                         WC(5, pz) * funcc(ii, iy, kc+2) + &
+                         WC(6, pz) * funcc(ii, iy, kc+3) )
+           end do
+           tmp_xyz_line(ii) = res_line
+        end do
+
+
+
+
+        ! 3. 【降维：X 向】最后在最内层只处理 X 方向的 6 点加权
+        ! 此时每个点的计算量从原来的 200+ 次乘法降到了仅 6 次
+        do i = imino, imaxo
+           px = pix(i)
+           ic = cix(i)
+           
+           ! 直接从预计算好的 line 中读取连续的 6 个点
+           ! ic-2 到 ic+3 对应原始 6 点算子
+           funf(i,j,k) = WC(1,px)*tmp_xyz_line(ic-2) + &
+                         WC(2,px)*tmp_xyz_line(ic-1) + &
+                         WC(3,px)*tmp_xyz_line(ic  ) + &
+                         WC(4,px)*tmp_xyz_line(ic+1) + &
+                         WC(5,px)*tmp_xyz_line(ic+2) + &
+                         WC(6,px)*tmp_xyz_line(ic+3)
+        end do
+     end do
+  end do
+#endif
   return
 
   end subroutine prolong3
@@ -2358,6 +2370,13 @@
 
   real*8,dimension(3) :: CD,FD
 
+  real*8 :: tmp_xz_plane(-1:extf(1), 6)
+  real*8 :: tmp_x_line(-1:extf(1))
+  integer :: fi, fj, fk, ii, jj, kk
+  integer :: fi_min, fi_max, ii_lo, ii_hi
+  integer :: fj_min, fj_max, fk_min, fk_max, jj_lo, jj_hi, kk_lo, kk_hi
+  logical :: need_full_symmetry
+
   if(wei.ne.3)then
      write(*,*)"prolongrestrict.f90::restrict3: this routine only surport 3 dimension"
      write(*,*)"dim = ",wei
@@ -2436,9 +2455,86 @@
           stop
   endif
 
-  call symmetry_bd(2,extf,funf,funff,SoA)
+  ! 仅计算 X 向最终写回所需的窗口：
+  ! func(i,j,k) 只访问 tmp_x_line(fi-2:fi+3)
+  fi_min = 2*(imino + lbc(1) - 1) - 1 - lbf(1) + 1
+  fi_max = 2*(imaxo + lbc(1) - 1) - 1 - lbf(1) + 1
+  fj_min = 2*(jmino + lbc(2) - 1) - 1 - lbf(2) + 1
+  fj_max = 2*(jmaxo + lbc(2) - 1) - 1 - lbf(2) + 1
+  fk_min = 2*(kmino + lbc(3) - 1) - 1 - lbf(3) + 1
+  fk_max = 2*(kmaxo + lbc(3) - 1) - 1 - lbf(3) + 1
+  ii_lo = fi_min - 2
+  ii_hi = fi_max + 3
+  jj_lo = fj_min - 2
+  jj_hi = fj_max + 3
+  kk_lo = fk_min - 2
+  kk_hi = fk_max + 3
+  if(ii_lo < -1 .or. ii_hi > extf(1) .or. &
+     jj_lo < -1 .or. jj_hi > extf(2) .or. &
+     kk_lo < -1 .or. kk_hi > extf(3))then
+      write(*,*)"restrict3: invalid stencil window"
+      write(*,*)"ii=",ii_lo,ii_hi," jj=",jj_lo,jj_hi," kk=",kk_lo,kk_hi
+      write(*,*)"extf=",extf
+      stop
+  endif
+  need_full_symmetry = (ii_lo < 1) .or. (jj_lo < 1) .or. (kk_lo < 1)
+  if(need_full_symmetry)then
+      call symmetry_bd(2,extf,funf,funff,SoA)
+  else
+      funff(ii_lo:ii_hi,jj_lo:jj_hi,kk_lo:kk_hi) = funf(ii_lo:ii_hi,jj_lo:jj_hi,kk_lo:kk_hi)
+  endif
 
 !~~~~~~> restriction start...
+do k = kmino, kmaxo
+    fk = 2*(k + lbc(3) - 1) - 1 - lbf(3) + 1
+
+    do j = jmino, jmaxo
+        fj = 2*(j + lbc(2) - 1) - 1 - lbf(2) + 1
+        
+        ! 优化点 1: 显式展开 Z 方向计算，减少循环开销
+        ! 确保 ii 循环是最内层且连续访问
+        !DIR$ VECTOR ALWAYS
+        do ii = ii_lo, ii_hi
+            ! 预计算当前 j 对应的 6 行在 Z 方向的压缩结果
+            ! 这里直接硬编码 jj 的偏移，彻底消除一层循环
+            tmp_xz_plane(ii, 1) = C1*(funff(ii,fj-2,fk-2)+funff(ii,fj-2,fk+3)) + &
+                                  C2*(funff(ii,fj-2,fk-1)+funff(ii,fj-2,fk+2)) + &
+                                  C3*(funff(ii,fj-2,fk  )+funff(ii,fj-2,fk+1))
+            tmp_xz_plane(ii, 2) = C1*(funff(ii,fj-1,fk-2)+funff(ii,fj-1,fk+3)) + &
+                                  C2*(funff(ii,fj-1,fk-1)+funff(ii,fj-1,fk+2)) + &
+                                  C3*(funff(ii,fj-1,fk  )+funff(ii,fj-1,fk+1))
+            tmp_xz_plane(ii, 3) = C1*(funff(ii,fj  ,fk-2)+funff(ii,fj  ,fk+3)) + &
+                                  C2*(funff(ii,fj  ,fk-1)+funff(ii,fj  ,fk+2)) + &
+                                  C3*(funff(ii,fj  ,fk  )+funff(ii,fj  ,fk+1))
+            tmp_xz_plane(ii, 4) = C1*(funff(ii,fj+1,fk-2)+funff(ii,fj+1,fk+3)) + &
+                                  C2*(funff(ii,fj+1,fk-1)+funff(ii,fj+1,fk+2)) + &
+                                  C3*(funff(ii,fj+1,fk  )+funff(ii,fj+1,fk+1))
+            tmp_xz_plane(ii, 5) = C1*(funff(ii,fj+2,fk-2)+funff(ii,fj+2,fk+3)) + &
+                                  C2*(funff(ii,fj+2,fk-1)+funff(ii,fj+2,fk+2)) + &
+                                  C3*(funff(ii,fj+2,fk  )+funff(ii,fj+2,fk+1))
+            tmp_xz_plane(ii, 6) = C1*(funff(ii,fj+3,fk-2)+funff(ii,fj+3,fk+3)) + &
+                                  C2*(funff(ii,fj+3,fk-1)+funff(ii,fj+3,fk+2)) + &
+                                  C3*(funff(ii,fj+3,fk  )+funff(ii,fj+3,fk+1))
+        end do
+
+        ! 优化点 2: 同样向量化 Y 方向压缩
+        !DIR$ VECTOR ALWAYS
+        do ii = ii_lo, ii_hi
+            tmp_x_line(ii) = C1*(tmp_xz_plane(ii, 1) + tmp_xz_plane(ii, 6)) + &
+                            C2*(tmp_xz_plane(ii, 2) + tmp_xz_plane(ii, 5)) + &
+                            C3*(tmp_xz_plane(ii, 3) + tmp_xz_plane(ii, 4))
+        end do
+
+        ! 优化点 3: 最终写入，利用已经缓存在 tmp_x_line 的数据
+        do i = imino, imaxo
+            fi = 2*(i + lbc(1) - 1) - 1 - lbf(1) + 1
+            func(i, j, k) = C1*(tmp_x_line(fi-2) + tmp_x_line(fi+3)) + &
+                            C2*(tmp_x_line(fi-1) + tmp_x_line(fi+2)) + &
+                            C3*(tmp_x_line(fi  ) + tmp_x_line(fi+1))
+        end do
+    end do
+end do
+#if 0
   do k = kmino,kmaxo
    do j = jmino,jmaxo
     do i = imino,imaxo
@@ -2462,7 +2558,7 @@
     enddo
    enddo
   enddo
-  
+#endif
   return
 
   end subroutine restrict3