删除diff_new.f90中冗余部分，方便后续工作

对fmisc.f90的polint修改
2026-02-08 00:54:23 +08:00 · 2026-02-07 01:56:44 +08:00
18 changed files with 1386 additions and 5091 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,3 @@
 __pycache__
 GW150914
-GW150914*
-.codex
-docs/
+GW150914-origin
--- a/AMSS_NCKU_ABEtest.py
+++ b/AMSS_NCKU_ABEtest.py
@@ -0,0 +1,445 @@
+
+##################################################################
+##
+## 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
+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()
+
+##################################################################
+
+
--- a/AMSS_NCKU_Program.py
+++ b/AMSS_NCKU_Program.py
@@ -9,19 +9,9 @@
 ##################################################################


-##################################################################
-
-## 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.
-if __name__ != '__main__':
-    import sys as _sys
-    _sys.exit(0)
-
-
-##################################################################
-
-## Print program introduction
+##################################################################
+
+## Print program introduction

 import print_information

@@ -432,36 +422,31 @@ print( " Plotting the txt and binary results data from the AMSS-NCKU simulation
 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_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot,   (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_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_tasks.append( ( plot_xiaoqu.generate_gravitational_wave_psi4_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_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_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 )
+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 " )
--- a/AMSS_NCKU_Verify_ASC26.py
+++ b/AMSS_NCKU_Verify_ASC26.py
@@ -0,0 +1,280 @@
+#!/usr/bin/env python3
+"""
+AMSS-NCKU GW150914 Simulation Regression Test Script
+
+Verification Requirements:
+1. XY-plane trajectory RMS error < 1% (Optimized vs. baseline, max of BH1 and BH2)
+2. ADM constraint violation < 2 (Grid Level 0)
+
+RMS Calculation Method:
+- Computes trajectory deviation on the XY plane independently for BH1 and BH2
+- For each black hole: RMS = sqrt((1/M) * sum((Δr_i / r_i^max)^2)) × 100%
+- Final RMS = max(RMS_BH1, RMS_BH2)
+
+Usage: python3 AMSS_NCKU_Verify_ASC26.py [output_dir]
+Default: output_dir = GW150914/AMSS_NCKU_output
+Reference: GW150914-origin (baseline simulation)
+"""
+
+import numpy as np
+import sys
+import os
+
+# ANSI Color Codes
+class Color:
+    GREEN = '\033[92m'
+    RED = '\033[91m'
+    YELLOW = '\033[93m'
+    BLUE = '\033[94m'
+    BOLD = '\033[1m'
+    RESET = '\033[0m'
+
+def get_status_text(passed):
+    if passed:
+        return f"{Color.GREEN}{Color.BOLD}PASS{Color.RESET}"
+    else:
+        return f"{Color.RED}{Color.BOLD}FAIL{Color.RESET}"
+
+def load_bh_trajectory(filepath):
+    """Load black hole trajectory data"""
+    data = np.loadtxt(filepath)
+    return {
+        'time': data[:, 0],
+        'x1': data[:, 1], 'y1': data[:, 2], 'z1': data[:, 3],
+        'x2': data[:, 4], 'y2': data[:, 5], 'z2': data[:, 6]
+    }
+
+
+def load_constraint_data(filepath):
+    """Load constraint violation data"""
+    data = []
+    with open(filepath, 'r') as f:
+        for line in f:
+            if line.startswith('#'):
+                continue
+            parts = line.split()
+            if len(parts) >= 8:
+                data.append([float(x) for x in parts[:8]])
+    return np.array(data)
+
+
+def calculate_rms_error(bh_data_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
+    """
+    # 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]
+
+    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]
+
+    # 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)
+
+    # 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)
+
+    # 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"
+
+    terms1 = (delta_r1[valid_mask1] / r1_max[valid_mask1])**2
+    rms_bh1 = np.sqrt(np.mean(terms1)) * 100
+
+    # 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"
+
+    terms2 = (delta_r2[valid_mask2] / r2_max[valid_mask2])**2
+    rms_bh2 = np.sqrt(np.mean(terms2)) * 100
+
+    # Final RMS is the maximum of BH1 and BH2
+    rms_final = max(rms_bh1, rms_bh2)
+
+    return rms_final, None
+
+
+def analyze_constraint_violation(constraint_data, n_levels=9):
+    """
+    Analyze ADM constraint violation
+    Return maximum constraint violation for Grid Level 0
+    """
+    # Extract Grid Level 0 data (first entry for each time step)
+    level0_data = constraint_data[::n_levels]
+
+    # Calculate maximum absolute value for each constraint
+    results = {
+        'Ham': np.max(np.abs(level0_data[:, 1])),
+        'Px': np.max(np.abs(level0_data[:, 2])),
+        'Py': np.max(np.abs(level0_data[:, 3])),
+        'Pz': np.max(np.abs(level0_data[:, 4])),
+        'Gx': np.max(np.abs(level0_data[:, 5])),
+        'Gy': np.max(np.abs(level0_data[:, 6])),
+        'Gz': np.max(np.abs(level0_data[:, 7]))
+    }
+
+    results['max_violation'] = max(results.values())
+    return results
+
+
+def print_header():
+    """Print report header"""
+    print("\n" + Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
+    print(Color.BOLD + "   AMSS-NCKU GW150914 Simulation Regression Test Report" + 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)
+
+    if error:
+        print(f"   {Color.RED}Error: {error}{Color.RESET}")
+        return False
+
+    passed = rms_rel < threshold
+
+    print(f"   RMS relative error: {rms_rel:.4f}%")
+    print(f"   Requirement:        < {threshold}%")
+    print(f"   Status:             {get_status_text(passed)}")
+
+    return passed
+
+
+def print_constraint_results(results, threshold=2.0):
+    """Print constraint violation results"""
+    print(f"\n{Color.BOLD}2. ADM Constraint Violation Analysis (Grid Level 0){Color.RESET}")
+    print("-" * 45)
+
+    names = ['Ham', 'Px', 'Py', 'Pz', 'Gx', 'Gy', 'Gz']
+    for i, name in enumerate(names):
+        print(f"   Max |{name:3}|: {results[name]:.6f}", end="   ")
+        if (i + 1) % 2 == 0: print()
+    if len(names) % 2 != 0: print()
+
+    passed = results['max_violation'] < threshold
+    
+    print(f"\n   Maximum violation:  {results['max_violation']:.6f}")
+    print(f"   Requirement:        < {threshold}")
+    print(f"   Status:             {get_status_text(passed)}")
+
+    return passed
+
+
+def print_summary(rms_passed, constraint_passed):
+    """Print summary"""
+    print("\n" + Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
+    print(Color.BOLD + "Verification Summary" + Color.RESET)
+    print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
+
+    all_passed = rms_passed and constraint_passed
+    
+    res_rms = get_status_text(rms_passed)
+    res_con = get_status_text(constraint_passed)
+
+    print(f"   [1] RMS trajectory check:         {res_rms}")
+    print(f"   [2] ADM constraint check:         {res_con}")
+    
+    final_status = f"{Color.GREEN}{Color.BOLD}ALL CHECKS PASSED{Color.RESET}" if all_passed else f"{Color.RED}{Color.BOLD}SOME CHECKS FAILED{Color.RESET}"
+    print(f"\n   Overall result: {final_status}")
+    print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET + "\n")
+
+    return all_passed
+
+
+def main():
+    # Determine 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)
+
+    # Analyze constraint violation
+    constraint_results = analyze_constraint_violation(constraint_data)
+    constraint_passed = print_constraint_results(constraint_results)
+
+    # Print summary
+    all_passed = print_summary(rms_passed, constraint_passed)
+
+    # Return exit code
+    sys.exit(0 if all_passed else 1)
+
+
+if __name__ == "__main__":
+    main()
+
--- a/AMSS_NCKU_source/TwoPunctures.C
+++ b/AMSS_NCKU_source/TwoPunctures.C
--- a/AMSS_NCKU_source/TwoPunctures.h
+++ b/AMSS_NCKU_source/TwoPunctures.h
@@ -1,8 +1,7 @@
+
 #ifndef TWO_PUNCTURES_H
 #define TWO_PUNCTURES_H

-#include <omp.h>
-
 #define StencilSize 19
 #define N_PlaneRelax 1
 #define NRELAX 200
@@ -33,7 +32,7 @@ private:
       int npoints_A, npoints_B, npoints_phi;

       double target_M_plus, target_M_minus;
-
+       
       double admMass;

       double adm_tol;
@@ -43,18 +42,6 @@ private:

       int ntotal;

-       // ===== Precomputed spectral derivative matrices =====
-       double *D1_A, *D2_A;
-       double *D1_B, *D2_B;
-       double *DF1_phi, *DF2_phi;
-
-       // ===== Pre-allocated workspace for LineRelax (per-thread) =====
-       int max_threads;
-       double **ws_diag_be, **ws_e_be, **ws_f_be, **ws_b_be, **ws_x_be;
-       double **ws_l_be, **ws_u_be, **ws_d_be, **ws_y_be;
-       double **ws_diag_al, **ws_e_al, **ws_f_al, **ws_b_al, **ws_x_al;
-       double **ws_l_al, **ws_u_al, **ws_d_al, **ws_y_al;
-
       struct parameters
       {
              int nvar, n1, n2, n3;
@@ -71,28 +58,6 @@ public:
                    int Newtonmaxit);
       ~TwoPunctures();

-       // 02/07: New/modified methods
-       void allocate_workspace();
-       void free_workspace();
-       void precompute_derivative_matrices();
-       void build_cheb_deriv_matrices(int n, double *D1, double *D2);
-       void build_fourier_deriv_matrices(int N, double *DF1, double *DF2);
-       void Derivatives_AB3_MatMul(int nvar, int n1, int n2, int n3, derivs v);
-       void ThomasAlgorithm_ws(int N, double *b, double *a, double *c, double *x, double *q,
-                                double *l, double *u_ws, double *d, double *y);
-       void LineRelax_be_omp(double *dv,
-                         int const i, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols, int **cols,
-                         double **JFD, int tid);
-       void LineRelax_al_omp(double *dv,
-                         int const j, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols,
-                         int **cols, double **JFD, int tid);
-       void relax_omp(double *dv, int const nvar, int const n1, int const n2, int const n3,
-                  double const *rhs, int const *ncols, int **cols, double **JFD);
-
       void Solve();
       void set_initial_guess(derivs v);
       int index(int i, int j, int k, int l, int a, int b, int c, int d);
@@ -151,11 +116,23 @@ public:
       double BY_KKofxyz(double x, double y, double z);
       void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix);
       void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u);
+       void relax(double *dv, int const nvar, int const n1, int const n2, int const n3,
+                  double const *rhs, int const *ncols, int **cols, double **JFD);
+       void LineRelax_be(double *dv,
+                         int const i, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols, int **cols,
+                         double **JFD);
       void JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
                         int n3, derivs dv, derivs u, double *values);
       void LinEquations(double A, double B, double X, double R,
                         double x, double r, double phi,
                         double y, double z, derivs dU, derivs U, double *values);
+       void LineRelax_al(double *dv,
+                         int const j, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols,
+                         int **cols, double **JFD);
       void ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q);
       void Save(char *fname);
       // provided by Vasileios Paschalidis (vpaschal@illinois.edu)
@@ -164,4 +141,4 @@ public:
       void SpecCoef(parameters par, int ivar, double *v, double *cf);
 };

-#endif /* TWO_PUNCTURES_H */
+#endif /* TWO_PUNCTURES_H */
--- a/AMSS_NCKU_source/diff_new.f90
+++ b/AMSS_NCKU_source/diff_new.f90
--- a/AMSS_NCKU_source/enforce_algebra.f90
+++ b/AMSS_NCKU_source/enforce_algebra.f90
@@ -17,62 +17,50 @@
  real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Axx,Axy,Axz
  real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Ayy,Ayz,Azz

-!~~~~~~~> Local variable:
-
-  integer :: i,j,k
-  real*8 :: lgxx,lgyy,lgzz,ldetg
-  real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
-  real*8 :: ltrA,lscale
-  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
-
-!~~~~~~>
-
-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-
-    lgxx = dxx(i,j,k) + ONE
-    lgyy = dyy(i,j,k) + ONE
-    lgzz = dzz(i,j,k) + ONE
-
-    ldetg =  lgxx * lgyy * lgzz &
-           + gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) &
-           + gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) &
-           - gxz(i,j,k) * lgyy * gxz(i,j,k) &
-           - gxy(i,j,k) * gxy(i,j,k) * lgzz &
-           - lgxx * gyz(i,j,k) * gyz(i,j,k)
-
-    lgupxx =   ( lgyy * lgzz - gyz(i,j,k) * gyz(i,j,k) ) / ldetg
-    lgupxy = - ( gxy(i,j,k) * lgzz - gyz(i,j,k) * gxz(i,j,k) ) / ldetg
-    lgupxz =   ( gxy(i,j,k) * gyz(i,j,k) - lgyy * gxz(i,j,k) ) / ldetg
-    lgupyy =   ( lgxx * lgzz - gxz(i,j,k) * gxz(i,j,k) ) / ldetg
-    lgupyz = - ( lgxx * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / ldetg
-    lgupzz =   ( lgxx * lgyy - gxy(i,j,k) * gxy(i,j,k) ) / ldetg
-
-    ltrA =         lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
-                 + lgupzz * Azz(i,j,k) &
-         + TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
-                 + lgupyz * Ayz(i,j,k))
-
-    Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
-    Axy(i,j,k) = Axy(i,j,k) - F1o3 * gxy(i,j,k) * ltrA
-    Axz(i,j,k) = Axz(i,j,k) - F1o3 * gxz(i,j,k) * ltrA
-    Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
-    Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * gyz(i,j,k) * ltrA
-    Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
-
-    lscale = ONE / ( ldetg ** F1o3 )
-
-    dxx(i,j,k) = lgxx * lscale - ONE
-    gxy(i,j,k) = gxy(i,j,k) * lscale
-    gxz(i,j,k) = gxz(i,j,k) * lscale
-    dyy(i,j,k) = lgyy * lscale - ONE
-    gyz(i,j,k) = gyz(i,j,k) * lscale
-    dzz(i,j,k) = lgzz * lscale - ONE
-
-  enddo
-  enddo
-  enddo
+!~~~~~~~> Local variable:
+  
+  real*8, dimension(ex(1),ex(2),ex(3)) :: trA,detg
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
+  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
+
+!~~~~~~>
+
+  gxx = dxx + ONE
+  gyy = dyy + ONE
+  gzz = dzz + ONE
+
+  detg =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
+          gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
+  gupxx =   ( gyy * gzz - gyz * gyz ) / detg
+  gupxy = - ( gxy * gzz - gyz * gxz ) / detg
+  gupxz =   ( gxy * gyz - gyy * gxz ) / detg
+  gupyy =   ( gxx * gzz - gxz * gxz ) / detg
+  gupyz = - ( gxx * gyz - gxy * gxz ) / detg
+  gupzz =   ( gxx * gyy - gxy * gxy ) / detg
+
+  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
+       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
+
+  Axx = Axx - F1o3 * gxx * trA
+  Axy = Axy - F1o3 * gxy * trA
+  Axz = Axz - F1o3 * gxz * trA
+  Ayy = Ayy - F1o3 * gyy * trA
+  Ayz = Ayz - F1o3 * gyz * trA
+  Azz = Azz - F1o3 * gzz * trA
+
+  detg = ONE / ( detg ** F1o3 ) 
+  
+  gxx = gxx * detg
+  gxy = gxy * detg
+  gxz = gxz * detg
+  gyy = gyy * detg
+  gyz = gyz * detg
+  gzz = gzz * detg
+
+  dxx = gxx - ONE
+  dyy = gyy - ONE
+  dzz = gzz - ONE

  return

@@ -93,72 +81,52 @@
  real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Axx,Axy,Axz
  real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Ayy,Ayz,Azz

-!~~~~~~~> Local variable:
-
-  integer :: i,j,k
-  real*8 :: lgxx,lgyy,lgzz,lscale
-  real*8 :: lgxy,lgxz,lgyz
-  real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
-  real*8 :: ltrA
-  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
-
-!~~~~~~>
-
-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-
-! for g: normalize determinant first
-    lgxx = dxx(i,j,k) + ONE
-    lgyy = dyy(i,j,k) + ONE
-    lgzz = dzz(i,j,k) + ONE
-    lgxy = gxy(i,j,k)
-    lgxz = gxz(i,j,k)
-    lgyz = gyz(i,j,k)
-
-    lscale =  lgxx * lgyy * lgzz + lgxy * lgyz * lgxz &
-            + lgxz * lgxy * lgyz - lgxz * lgyy * lgxz &
-            - lgxy * lgxy * lgzz - lgxx * lgyz * lgyz
-
-    lscale = ONE / ( lscale ** F1o3 )
-
-    lgxx = lgxx * lscale
-    lgxy = lgxy * lscale
-    lgxz = lgxz * lscale
-    lgyy = lgyy * lscale
-    lgyz = lgyz * lscale
-    lgzz = lgzz * lscale
-
-    dxx(i,j,k) = lgxx - ONE
-    gxy(i,j,k) = lgxy
-    gxz(i,j,k) = lgxz
-    dyy(i,j,k) = lgyy - ONE
-    gyz(i,j,k) = lgyz
-    dzz(i,j,k) = lgzz - ONE
-
-! for A: trace-free using normalized metric (det=1, no division needed)
-    lgupxx =   ( lgyy * lgzz - lgyz * lgyz )
-    lgupxy = - ( lgxy * lgzz - lgyz * lgxz )
-    lgupxz =   ( lgxy * lgyz - lgyy * lgxz )
-    lgupyy =   ( lgxx * lgzz - lgxz * lgxz )
-    lgupyz = - ( lgxx * lgyz - lgxy * lgxz )
-    lgupzz =   ( lgxx * lgyy - lgxy * lgxy )
-
-    ltrA =         lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
-                 + lgupzz * Azz(i,j,k) &
-         + TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
-                 + lgupyz * Ayz(i,j,k))
-
-    Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
-    Axy(i,j,k) = Axy(i,j,k) - F1o3 * lgxy * ltrA
-    Axz(i,j,k) = Axz(i,j,k) - F1o3 * lgxz * ltrA
-    Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
-    Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * lgyz * ltrA
-    Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
-
-  enddo
-  enddo
-  enddo
+!~~~~~~~> Local variable:
+  
+  real*8, dimension(ex(1),ex(2),ex(3)) :: trA
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
+  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
+
+!~~~~~~>
+
+  gxx = dxx + ONE
+  gyy = dyy + ONE
+  gzz = dzz + ONE
+! for g
+  gupzz =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
+           gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
+
+  gupzz = ONE / ( gupzz ** F1o3 ) 
+  
+  gxx = gxx * gupzz
+  gxy = gxy * gupzz
+  gxz = gxz * gupzz
+  gyy = gyy * gupzz
+  gyz = gyz * gupzz
+  gzz = gzz * gupzz
+
+  dxx = gxx - ONE
+  dyy = gyy - ONE
+  dzz = gzz - ONE
+! for A  
+
+  gupxx =   ( gyy * gzz - gyz * gyz )
+  gupxy = - ( gxy * gzz - gyz * gxz )
+  gupxz =   ( gxy * gyz - gyy * gxz )
+  gupyy =   ( gxx * gzz - gxz * gxz )
+  gupyz = - ( gxx * gyz - gxy * gxz )
+  gupzz =   ( gxx * gyy - gxy * gxy )
+
+  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
+       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
+
+  Axx = Axx - F1o3 * gxx * trA
+  Axy = Axy - F1o3 * gxy * trA
+  Axz = Axz - F1o3 * gxz * trA
+  Ayy = Ayy - F1o3 * gyy * trA
+  Ayz = Ayz - F1o3 * gyz * trA
+  Azz = Azz - F1o3 * gzz * trA

  return

--- a/AMSS_NCKU_source/fmisc.f90
+++ b/AMSS_NCKU_source/fmisc.f90
@@ -324,7 +324,8 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
   enddo
@@ -349,7 +350,8 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1)
@@ -377,7 +379,8 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1)
@@ -883,20 +886,17 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)

  integer::i

-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
-   do i=0,ord-1
-      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
-   enddo
-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
-   do i=0,ord-1
-      funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
-   enddo
-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
-   do i=0,ord-1
-      funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
-   enddo
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+   do i=0,ord-1
+      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
+   enddo
+   do i=0,ord-1
+      funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
+   enddo
+   do i=0,ord-1
+      funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
+   enddo

 end subroutine symmetry_bd

@@ -912,7 +912,8 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
@@ -940,7 +941,8 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
@@ -1111,353 +1113,142 @@ end subroutine d2dump
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! common code for cell and vertex
-!------------------------------------------------------------------------------
-! Lagrangian polynomial interpolation
-!------------------------------------------------------------------------------
-#ifndef POLINT6_USE_BARYCENTRIC
-#define POLINT6_USE_BARYCENTRIC 1
-#endif
-
-!DIR$ ATTRIBUTES FORCEINLINE :: polint6_neville
-  subroutine polint6_neville(xa, ya, x, y, dy)
-  implicit none
-
-  real*8, dimension(6), intent(in) :: xa, ya
-  real*8, intent(in) :: x
-  real*8, intent(out) :: y, dy
-
-  integer :: i, m, ns, n_m
-  real*8, dimension(6) :: c, d, ho
-  real*8 :: dif, dift, hp, h, den_val
-
-  c = ya
-  d = ya
-  ho = xa - x
-
-  ns = 1
-  dif = abs(x - xa(1))
-
-  do i = 2, 6
-    dift = abs(x - xa(i))
-    if (dift < dif) then
-      ns = i
-      dif = dift
-    end if
-  end do
-
-  y = ya(ns)
-  ns = ns - 1
-
-  do m = 1, 5
-    n_m = 6 - m
-    do i = 1, n_m
-      hp = ho(i)
-      h  = ho(i+m)
-      den_val = hp - h
-
-      if (den_val == 0.0d0) then
-        write(*,*) 'failure in polint for point',x
-        write(*,*) 'with input points: ',xa
-        stop
-      end if
-
-      den_val = (c(i+1) - d(i)) / den_val
-
-      d(i) = h * den_val
-      c(i) = hp * den_val
-    end do
-
-    if (2 * ns < n_m) then
-      dy = c(ns + 1)
-    else
-      dy = d(ns)
-      ns = ns - 1
-    end if
-    y = y + dy
-  end do
-
-  return
-  end subroutine polint6_neville
-
-!DIR$ ATTRIBUTES FORCEINLINE :: polint6_barycentric
-  subroutine polint6_barycentric(xa, ya, x, y, dy)
-  implicit none
-
-  real*8, dimension(6), intent(in) :: xa, ya
-  real*8, intent(in) :: x
-  real*8, intent(out) :: y, dy
-
-  integer :: i, j
-  logical :: is_uniform
-  real*8, dimension(6) :: lambda
-  real*8 :: dx, den_i, term, num, den, step, tol
-  real*8, parameter :: c_uniform(6) = (/ -1.d0, 5.d0, -10.d0, 10.d0, -5.d0, 1.d0 /)
-
-  do i = 1, 6
-    if (x == xa(i)) then
-      y = ya(i)
-      dy = 0.d0
-      return
-    end if
-  end do
-
-  step = xa(2) - xa(1)
-  is_uniform = (step /= 0.d0)
-  if (is_uniform) then
-    tol = 64.d0 * epsilon(1.d0) * max(1.d0, abs(step))
-    do i = 3, 6
-      if (abs((xa(i) - xa(i-1)) - step) > tol) then
-        is_uniform = .false.
-        exit
-      end if
-    end do
-  end if
-
-  if (is_uniform) then
-    num = 0.d0
-    den = 0.d0
-    do i = 1, 6
-      term = c_uniform(i) / (x - xa(i))
-      num = num + term * ya(i)
-      den = den + term
-    end do
-    y = num / den
-    dy = 0.d0
-    return
-  end if
-
-  do i = 1, 6
-    den_i = 1.d0
-    do j = 1, 6
-      if (j /= i) then
-        dx = xa(i) - xa(j)
-        if (dx == 0.0d0) then
-          write(*,*) 'failure in polint for point',x
-          write(*,*) 'with input points: ',xa
-          stop
-        end if
-        den_i = den_i * dx
-      end if
-    end do
-    lambda(i) = 1.d0 / den_i
-  end do
-
-  num = 0.d0
-  den = 0.d0
-  do i = 1, 6
-    term = lambda(i) / (x - xa(i))
-    num = num + term * ya(i)
-    den = den + term
-  end do
-
-  y = num / den
-  dy = 0.d0
-
-  return
-  end subroutine polint6_barycentric
-
-!DIR$ ATTRIBUTES FORCEINLINE :: polint
-  subroutine polint(xa, ya, x, y, dy, ordn)
-  implicit none
-
-  integer, intent(in) :: ordn
-  real*8, dimension(ordn), intent(in) :: xa, ya
-  real*8, intent(in) :: x
-  real*8, intent(out) :: y, dy
-
-  integer :: i, m, ns, n_m
-  real*8, dimension(ordn) :: c, d, ho
-  real*8 :: dif, dift, hp, h, den_val
-
-  if (ordn == 6) then
-#if POLINT6_USE_BARYCENTRIC
-    call polint6_barycentric(xa, ya, x, y, dy)
-#else
-    call polint6_neville(xa, ya, x, y, dy)
-#endif
-    return
-  end if
-
-  c = ya
-  d = ya
-  ho = xa - x
-
-  ns = 1
-  dif = abs(x - xa(1))
-
-  do i = 2, ordn
-    dift = abs(x - xa(i))
-    if (dift < dif) then
-      ns = i
-      dif = dift
-    end if
-  end do
-
-  y = ya(ns)
-  ns = ns - 1
-
-  do m = 1, ordn - 1
-    n_m = ordn - m
-    do i = 1, n_m
-      hp = ho(i)
-      h  = ho(i+m)
-      den_val = hp - h
-
-      if (den_val == 0.0d0) then
-        write(*,*) 'failure in polint for point',x
-        write(*,*) 'with input points: ',xa
-        stop
-      end if
-
-      den_val = (c(i+1) - d(i)) / den_val
-
-      d(i) = h * den_val
-      c(i) = hp * den_val
-    end do
-
-    if (2 * ns < n_m) then
-      dy = c(ns + 1)
-    else
-      dy = d(ns)
-      ns = ns - 1
-    end if
-    y = y + dy
-  end do
-
-  return
-  end subroutine polint
-!------------------------------------------------------------------------------
-! Compute Lagrange interpolation basis weights for one target point.
-!------------------------------------------------------------------------------
-!DIR$ ATTRIBUTES FORCEINLINE :: polint_lagrange_weights
-  subroutine polint_lagrange_weights(xa, x, w, ordn)
-  implicit none
-
-  integer, intent(in) :: ordn
-  real*8, dimension(1:ordn), intent(in) :: xa
-  real*8, intent(in) :: x
-  real*8, dimension(1:ordn), intent(out) :: w
-
-  integer :: i, j
-  real*8 :: num, den, dx
-
-  do i = 1, ordn
-    num = 1.d0
-    den = 1.d0
-    do j = 1, ordn
-      if (j /= i) then
-        dx = xa(i) - xa(j)
-        if (dx == 0.0d0) then
-          write(*,*) 'failure in polint for point',x
-          write(*,*) 'with input points: ',xa
-          stop
-        end if
-        num = num * (x - xa(j))
-        den = den * dx
-      end if
-    end do
-    w(i) = num / den
-  end do
-
-  return
-  end subroutine polint_lagrange_weights
-!------------------------------------------------------------------------------
-!
-! interpolation in 2 dimensions, follow yx order
-!
-!------------------------------------------------------------------------------
-  subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
-  implicit none
-
-  integer,intent(in) :: ordn
-  real*8, dimension(1:ordn), intent(in) :: x1a,x2a
-  real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
-  real*8, intent(in) :: x1,x2
-  real*8, intent(out) :: y,dy
-
-#ifdef POLINT_LEGACY_ORDER
-  integer  :: i,m
-  real*8, dimension(ordn) :: ymtmp
-  real*8, dimension(ordn) :: yntmp
-
-  m=size(x1a)
-  do i=1,m
-    yntmp=ya(i,:)
-    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-  end do
-  call polint(x1a,ymtmp,x1,y,dy,ordn)
-#else
-  integer  :: j
-  real*8, dimension(ordn) :: ymtmp
-  real*8 :: dy_temp
-
-  do j=1,ordn
-    call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
-  end do
-  call polint(x2a, ymtmp, x2, y, dy, ordn)
-#endif
-
-  return
-  end subroutine polin2
-!------------------------------------------------------------------------------
-!
-! interpolation in 3 dimensions, follow zyx order
-!
-!------------------------------------------------------------------------------
-  subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
-  implicit none
-
-  integer,intent(in) :: ordn
-  real*8, dimension(1:ordn), intent(in) :: x1a,x2a,x3a
-  real*8, dimension(1:ordn,1:ordn,1:ordn), intent(in) :: ya
-  real*8, intent(in) :: x1,x2,x3
-  real*8, intent(out) :: y,dy
-
-#ifdef POLINT_LEGACY_ORDER
-  integer  :: i,j,m,n
-  real*8, dimension(ordn,ordn) :: yatmp
-  real*8, dimension(ordn) :: ymtmp
-  real*8, dimension(ordn) :: yntmp
-  real*8, dimension(ordn) :: yqtmp
-
-  m=size(x1a)
-  n=size(x2a)
-  do i=1,m
-   do j=1,n
-    yqtmp=ya(i,j,:)
-    call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
-   end do
-    yntmp=yatmp(i,:)
-    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-  end do
-  call polint(x1a,ymtmp,x1,y,dy,ordn)
-#else
-  integer  :: i, j, k
-  real*8, dimension(ordn) :: w1, w2
-  real*8, dimension(ordn) :: ymtmp
-  real*8 :: yx_sum, x_sum
-
-  call polint_lagrange_weights(x1a, x1, w1, ordn)
-  call polint_lagrange_weights(x2a, x2, w2, ordn)
-
-  do k = 1, ordn
-    yx_sum = 0.d0
-    do j = 1, ordn
-      x_sum = 0.d0
-      do i = 1, ordn
-        x_sum = x_sum + w1(i) * ya(i,j,k)
-      end do
-      yx_sum = yx_sum + w2(j) * x_sum
-    end do
-    ymtmp(k) = yx_sum
-  end do
-
-  call polint(x3a, ymtmp, x3, y, dy, ordn)
-#endif
-
-  return
-  end subroutine polin3
+! common code for cell and vertex
+!------------------------------------------------------------------------------
+! Lagrangian polynomial interpolation
+!------------------------------------------------------------------------------
+ subroutine polint(xa, ya, x, y, dy, ordn)
+  implicit none
+
+  integer, intent(in) :: ordn
+  real*8, dimension(ordn), intent(in) :: xa, ya
+  real*8, intent(in) :: x
+  real*8, intent(out) :: y, dy
+
+  integer :: i, m, ns, n_m
+  real*8, dimension(ordn) :: c, d, ho
+  real*8 :: dif, dift, hp, h, den_val
+
+  ! Initialization
+  c = ya
+  d = ya
+  ho = xa - x
+  
+  ns = 1
+  dif = abs(x - xa(1))
+  
+  ! Find the index of the closest table entry
+  do i = 2, ordn
+    dift = abs(x - xa(i))
+    if (dift < dif) then
+      ns = i
+      dif = dift
+    end if
+  end do
+
+  y = ya(ns)
+  ns = ns - 1
+  
+  ! Main Neville's algorithm loop
+  do m = 1, ordn - 1
+    n_m = ordn - m
+    do i = 1, n_m
+      hp = ho(i)
+      h  = ho(i+m)
+      den_val = hp - h
+      
+      ! Check for division by zero locally
+      if (den_val == 0.0d0) then
+        write(*,*) 'failure in polint for point',x
+        write(*,*) 'with input points: ',xa
+        stop
+      end if
+      
+      ! Reuse den_val to avoid redundant divisions
+      den_val = (c(i+1) - d(i)) / den_val
+      
+      ! Update c and d in place
+      d(i) = h * den_val
+      c(i) = hp * den_val
+    end do
+
+    ! Decide which path (up or down the tableau) to take
+    if (2 * ns < n_m) then
+      dy = c(ns + 1)
+    else
+      dy = d(ns)
+      ns = ns - 1
+    end if
+    y = y + dy
+  end do
+
+  return
+  end subroutine polint
+!------------------------------------------------------------------------------
+!
+! interpolation in 2 dimensions, follow yx order
+!
+!------------------------------------------------------------------------------
+subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
+    implicit none
+    integer,intent(in) :: ordn
+    real*8, dimension(ordn), intent(in) :: x1a,x2a
+    real*8, dimension(ordn,ordn), intent(in) :: ya
+    real*8, intent(in) :: x1,x2
+    real*8, intent(out) :: y,dy
+
+    integer  :: j
+    real*8, dimension(ordn) :: ymtmp
+    real*8 :: dy_temp ! Local variable to prevent overwriting result
+
+    ! Optimized sequence: Loop over columns (j) 
+    ! ya(:,j) is a contiguous memory block in Fortran
+    do j=1,ordn
+      call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
+    end do
+
+    ! Final interpolation on the results
+    call polint(x2a, ymtmp, x2, y, dy, ordn)
+
+    return
+  end subroutine polin2
+!------------------------------------------------------------------------------
+!
+! interpolation in 3 dimensions, follow zyx order
+!
+!------------------------------------------------------------------------------
+subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
+    implicit none
+    integer,intent(in) :: ordn
+    real*8, dimension(ordn), intent(in) :: x1a,x2a,x3a
+    real*8, dimension(ordn,ordn,ordn), intent(in) :: ya
+    real*8, intent(in) :: x1,x2,x3
+    real*8, intent(out) :: y,dy
+
+    integer  :: j, k
+    real*8, dimension(ordn,ordn) :: yatmp
+    real*8, dimension(ordn) :: ymtmp
+    real*8 :: dy_temp
+
+    ! Sequence change: Process the contiguous first dimension (x1) first.
+    ! We loop through the 'slow' planes (j, k) to extract 'fast' columns.
+    do k=1,ordn
+      do j=1,ordn
+        ! ya(:,j,k) is contiguous; much faster than ya(i,j,:)
+        call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
+      end do
+    end do
+
+    ! Now process the second dimension
+    do k=1,ordn
+      call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
+    end do
+
+    ! Final dimension
+    call polint(x3a, ymtmp, x3, y, dy, ordn)
+
+    return
+  end subroutine polin3
 !--------------------------------------------------------------------------------------
 ! calculate L2norm  
  subroutine l2normhelper(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
@@ -1476,9 +1267,7 @@ end subroutine d2dump
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
-  integer::i,j,k,n_elements
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
+  integer::i,j,k

  dX = X(2) - X(1)
  dY = Y(2) - Y(1)
@@ -1502,91 +1291,15 @@ if(dabs(X(1)-xmin) < dX) imin = 1
 if(dabs(Y(1)-ymin) < dY) jmin = 1
 if(dabs(Z(1)-zmin) < dZ) kmin = 1

-  n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
-  allocate(f_flat(n_elements))
-  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
-  f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
-  deallocate(f_flat)
+f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))

 f_out = f_out*dX*dY*dZ

  return

  end subroutine l2normhelper
-!--------------------------------------------------------------------------------------
-  subroutine l2normhelper7(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
-                           f1,f2,f3,f4,f5,f6,f7,f_out,gw)
-
-  implicit none
-!~~~~~~> Input parameters:
-  integer,intent(in ):: ex(1:3)
-  real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)),xmin,ymin,zmin,xmax,ymax,zmax
-  integer,intent(in)::gw
-  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f1,f2,f3,f4,f5,f6,f7
-  real*8, intent(out) :: f_out(7)
-!~~~~~~> Other variables:
-
-  real*8            :: dX, dY, dZ
-  integer::imin,jmin,kmin
-  integer::imax,jmax,kmax
-  integer::i,j,k
-  real*8 :: s1,s2,s3,s4,s5,s6,s7
-
-  dX = X(2) - X(1)
-  dY = Y(2) - Y(1)
-  dZ = Z(2) - Z(1)
-
-   imin = gw+1
-   jmin = gw+1
-   kmin = gw+1
-
-   imax = ex(1) - gw
-   jmax = ex(2) - gw
-   kmax = ex(3) - gw
-
-if(dabs(X(ex(1))-xmax) < dX) imax = ex(1)
-if(dabs(Y(ex(2))-ymax) < dY) jmax = ex(2)
-if(dabs(Z(ex(3))-zmax) < dZ) kmax = ex(3)
-if(dabs(X(1)-xmin) < dX) imin = 1
-if(dabs(Y(1)-ymin) < dY) jmin = 1
-if(dabs(Z(1)-zmin) < dZ) kmin = 1
-
-  s1 = 0.d0
-  s2 = 0.d0
-  s3 = 0.d0
-  s4 = 0.d0
-  s5 = 0.d0
-  s6 = 0.d0
-  s7 = 0.d0
-
-  do k=kmin,kmax
-    do j=jmin,jmax
-!DIR$ SIMD REDUCTION(+:s1,s2,s3,s4,s5,s6,s7)
-      do i=imin,imax
-        s1 = s1 + f1(i,j,k)*f1(i,j,k)
-        s2 = s2 + f2(i,j,k)*f2(i,j,k)
-        s3 = s3 + f3(i,j,k)*f3(i,j,k)
-        s4 = s4 + f4(i,j,k)*f4(i,j,k)
-        s5 = s5 + f5(i,j,k)*f5(i,j,k)
-        s6 = s6 + f6(i,j,k)*f6(i,j,k)
-        s7 = s7 + f7(i,j,k)*f7(i,j,k)
-      enddo
-    enddo
-  enddo
-
-  f_out(1) = s1*dX*dY*dZ
-  f_out(2) = s2*dX*dY*dZ
-  f_out(3) = s3*dX*dY*dZ
-  f_out(4) = s4*dX*dY*dZ
-  f_out(5) = s5*dX*dY*dZ
-  f_out(6) = s6*dX*dY*dZ
-  f_out(7) = s7*dX*dY*dZ
-
-  return
-
-  end subroutine l2normhelper7
-!--------------------------------------------------------------------------------------
-! calculate L2norm especially for shell Blocks
+!--------------------------------------------------------------------------------------
+! calculate L2norm especially for shell Blocks
  subroutine l2normhelper_sh(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
                          f,f_out,gw,ogw,Symmetry)

@@ -1603,9 +1316,7 @@ if(dabs(Z(1)-zmin) < dZ) kmin = 1
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
-  integer::i,j,k,n_elements
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
+  integer::i,j,k

  real*8 :: PIo4

@@ -1668,11 +1379,7 @@ if(Symmetry==2)then
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif

-  n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
-  allocate(f_flat(n_elements))
-  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
-  f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
-  deallocate(f_flat)
+f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))

 f_out = f_out*dX*dY*dZ

@@ -1699,9 +1406,7 @@ f_out = f_out*dX*dY*dZ
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
-  integer::i,j,k
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
+  integer::i,j,k

  real*8 :: PIo4

@@ -1764,11 +1469,11 @@ if(Symmetry==2)then
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif

-Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
-  allocate(f_flat(Nout))
-  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [Nout])
-  f_out = DDOT(Nout, f_flat, 1, f_flat, 1)
-  deallocate(f_flat)
+f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
+
+f_out = f_out
+
+Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)

  return

@@ -1869,12 +1574,9 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
 !       ^
 ! f=3/8*f_1 + 3/4*f_2 - 1/8*f_3

-  real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
-  integer :: i,j,k
-
-  do concurrent (k=1:ext(3), j=1:ext(2), i=1:ext(1))
-    fout(i,j,k) = C1*f1(i,j,k)+C2*f2(i,j,k)+C3*f3(i,j,k)
-  end do
+  real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
+
+  fout = C1*f1+C2*f2+C3*f3

  return

@@ -1968,8 +1670,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
  real*8, dimension(ORDN,ORDN,ORDN) :: ya
  real*8, dimension(ORDN,ORDN) :: tmp2
  real*8, dimension(ORDN) :: tmp1
-  real*8, dimension(3) :: SoAh
-  real*8, external :: DDOT
+  real*8, dimension(3) :: SoAh

 ! +1 because c++ gives 0 for first point
  cxB = inds+1  
@@ -2015,7 +1716,10 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
    tmp1 = tmp1 + coef(ORDN+m)*tmp2(:,m)
  enddo

-  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*tmp1(m)
+  enddo

  return

@@ -2044,8 +1748,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
  integer,dimension(2) :: cxB,cxT
  real*8, dimension(ORDN,ORDN) :: ya
  real*8, dimension(ORDN) :: tmp1
-  real*8, dimension(2) :: SoAh
-  real*8, external :: DDOT
+  real*8, dimension(2) :: SoAh

 ! +1 because c++ gives 0 for first point
  cxB = inds(1:2)+1  
@@ -2080,7 +1783,10 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
    tmp1 = tmp1 + coef(ORDN+m)*ya(:,m)
  enddo

-  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*tmp1(m)
+  enddo

  return

@@ -2106,12 +1812,11 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
 !~~~~~~> Other parameters:

  real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,ex(3)) :: fh
-  integer :: m
-  integer :: cxB,cxT
-  real*8, dimension(ORDN) :: ya
-  real*8 :: SoAh
-  integer,dimension(3) :: inds
-  real*8, external :: DDOT
+  integer :: m
+  integer :: cxB,cxT
+  real*8, dimension(ORDN) :: ya
+  real*8 :: SoAh
+  integer,dimension(3) :: inds

 ! +1 because c++ gives 0 for first point
  inds = indsi + 1
@@ -2172,7 +1877,10 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
          write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
  endif

-  f_int = DDOT(ORDN, coef, 1, ya, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*ya(m)
+  enddo

  return

@@ -2408,28 +2116,20 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
  implicit none
  integer,intent(in) :: N

-  real*8 :: gont
-
-  integer :: i
-  real*8, parameter, dimension(0:20) :: fact_table = [ &
-    1.d0, 1.d0, 2.d0, 6.d0, 24.d0, 120.d0, 720.d0, 5040.d0, 40320.d0, &
-    362880.d0, 3628800.d0, 39916800.d0, 479001600.d0, 6227020800.d0, &
-    87178291200.d0, 1307674368000.d0, 20922789888000.d0, &
-    355687428096000.d0, 6402373705728000.d0, 121645100408832000.d0, &
-    2432902008176640000.d0 ]
+  real*8 :: gont
+
+  integer :: i

 ! sanity check
-  if(N < 0)then
-     write(*,*) "ffact: error input for factorial"
-     gont = 1.d0
-     return
-  endif
-
-  if(N <= 20)then
-     gont = fact_table(N)
-  else
-     gont = exp(log_gamma(dble(N+1)))
-  endif
+  if(N < 0)then
+     write(*,*) "ffact: error input for factorial"
+     return
+  endif
+
+  gont = 1.d0
+  do i=1,N
+     gont = gont*i
+  enddo

  return

@@ -2563,3 +2263,4 @@ subroutine find_maximum(ext,X,Y,Z,fun,val,pos,llb,uub)
  return

 end subroutine
+
--- a/AMSS_NCKU_source/fmisc.h
+++ b/AMSS_NCKU_source/fmisc.h
@@ -12,10 +12,9 @@
 #define f_global_interpind global_interpind
 #define f_global_interpind2d global_interpind2d
 #define f_global_interpind1d global_interpind1d
-#define f_l2normhelper l2normhelper
-#define f_l2normhelper7 l2normhelper7
-#define f_l2normhelper_sh l2normhelper_sh
-#define f_l2normhelper_sh_rms l2normhelper_sh_rms
+#define f_l2normhelper l2normhelper
+#define f_l2normhelper_sh l2normhelper_sh
+#define f_l2normhelper_sh_rms l2normhelper_sh_rms
 #define f_average average
 #define f_average3 average3
 #define f_average2 average2
@@ -42,10 +41,9 @@
 #define f_global_interpind GLOBAL_INTERPIND
 #define f_global_interpind2d GLOBAL_INTERPIND2D
 #define f_global_interpind1d GLOBAL_INTERPIND1D
-#define f_l2normhelper L2NORMHELPER
-#define f_l2normhelper7 L2NORMHELPER7
-#define f_l2normhelper_sh L2NORMHELPER_SH
-#define f_l2normhelper_sh_rms L2NORMHELPER_SH_RMS
+#define f_l2normhelper L2NORMHELPER
+#define f_l2normhelper_sh L2NORMHELPER_SH
+#define f_l2normhelper_sh_rms L2NORMHELPER_SH_RMS
 #define f_average AVERAGE
 #define f_average3 AVERAGE3
 #define f_average2 AVERAGE2
@@ -72,10 +70,9 @@
 #define f_global_interpind global_interpind_
 #define f_global_interpind2d global_interpind2d_
 #define f_global_interpind1d global_interpind1d_
-#define f_l2normhelper l2normhelper_
-#define f_l2normhelper7 l2normhelper7_
-#define f_l2normhelper_sh l2normhelper_sh_
-#define f_l2normhelper_sh_rms l2normhelper_sh_rms_
+#define f_l2normhelper l2normhelper_
+#define f_l2normhelper_sh l2normhelper_sh_
+#define f_l2normhelper_sh_rms l2normhelper_sh_rms_
 #define f_average average_
 #define f_average3 average3_
 #define f_average2 average2_
@@ -159,30 +156,21 @@ extern "C"
 							  int *, double *, int &, int &);
 }

-extern "C"
-{
-	void f_l2normhelper(int *, double *, double *, double *,
-						double &, double &, double &,
-						double &, double &, double &,
-						double *, double &, int &);
-}
-
-extern "C"
-{
-	void f_l2normhelper7(int *, double *, double *, double *,
-						 double &, double &, double &,
-						 double &, double &, double &,
-						 double *, double *, double *, double *,
-						 double *, double *, double *, double *, int &);
-}
-
-extern "C"
-{
-	void f_l2normhelper_sh(int *, double *, double *, double *,
-						   double &, double &, double &,
-						   double &, double &, double &,
-						   double *, double &, int &, int &, int &);
-}
+extern "C"
+{
+	void f_l2normhelper(int *, double *, double *, double *,
+						double &, double &, double &,
+						double &, double &, double &,
+						double *, double &, int &);
+}
+
+extern "C"
+{
+	void f_l2normhelper_sh(int *, double *, double *, double *,
+						   double &, double &, double &,
+						   double &, double &, double &,
+						   double *, double &, int &, int &, int &);
+}

 extern "C"
 {
--- a/AMSS_NCKU_source/macrodef.h
+++ b/AMSS_NCKU_source/macrodef.h
@@ -2,7 +2,7 @@
 #ifndef MICRODEF_H
 #define MICRODEF_H

-#include "macrodef.fh"
+#include "microdef.fh"

 // application parameters

--- a/AMSS_NCKU_source/makefile
+++ b/AMSS_NCKU_source/makefile
@@ -1,25 +1,11 @@
-
-
-include makefile.inc
-
-## polint(ordn=6) kernel selector:
-##   1 (default): barycentric fast path
-##   0          : fallback to Neville path
-POLINT6_USE_BARY ?= 1
-POLINT6_FLAG = -DPOLINT6_USE_BARYCENTRIC=$(POLINT6_USE_BARY)
-
-ARCH_OPT = -march=x86-64-v4
-CXXAPPFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
-              -Dfortran3 -Dnewc -I${MKLROOT}/include
-f90appflags = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
-              -align array64byte -fpp -I${MKLROOT}/include $(POLINT6_FLAG)
-TP_OPTFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
-              -Dfortran3 -Dnewc -I${MKLROOT}/include
-
-.SUFFIXES: .o .f90 .C .for .cu
-
-.f90.o:
-	$(f90) $(f90appflags) -c $< -o $@
+
+
+include makefile.inc
+
+.SUFFIXES: .o .f90 .C .for .cu
+
+.f90.o:
+	$(f90) $(f90appflags) -c $< -o $@

 .C.o:
 	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
@@ -27,14 +13,8 @@ TP_OPTFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
 .for.o:
 	$(f77) -c $< -o $@

-.cu.o:
-	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)
-
-TwoPunctures.o: TwoPunctures.C
-	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
-
-TwoPunctureABE.o: TwoPunctureABE.C
-	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
+.cu.o:
+	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)

 # Input files
 C++FILES = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o\
@@ -115,8 +95,8 @@ ABE: $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS)
 ABEGPU: $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
 	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)

-TwoPunctureABE: $(TwoPunctureFILES)
-	$(CLINKER) $(TP_OPTFLAGS) -qopenmp -o $@ $(TwoPunctureFILES) $(LDLIBS)
+TwoPunctureABE: $(TwoPunctureFILES)
+	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(TwoPunctureFILES) $(LDLIBS)

 clean:
 	rm *.o ABE ABEGPU TwoPunctureABE make.log -f
--- a/AMSS_NCKU_source/makefile.inc
+++ b/AMSS_NCKU_source/makefile.inc
@@ -1,34 +1,33 @@
-
 ## GCC version (commented out)
 ## filein  = -I/usr/include -I/usr/lib/x86_64-linux-gnu/mpich/include -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
 ## filein  = -I/usr/include/ -I/usr/include/openmpi-x86_64/ -I/usr/lib/x86_64-linux-gnu/openmpi/include/ -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
 ## LDLIBS  = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran -lmpi -lgfortran

-## Intel oneAPI version with oneMKL
+## Intel oneAPI version with oneMKL (Optimized for performance)
 filein  = -I/usr/include/ -I${MKLROOT}/include

-## Use sequential oneMKL to avoid introducing extra OpenMP behavior into ABE.
-LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl -liomp5
-
-## Optional Intel oneTBB allocator, kept aligned with main's build environment.
-USE_TBBMALLOC ?= 1
-TBBMALLOC_SO ?= /home/intel/oneapi/2025.3/lib/libtbbmalloc.so
-ifneq ($(wildcard $(TBBMALLOC_SO)),)
-TBBMALLOC_LIBS = -Wl,--no-as-needed $(TBBMALLOC_SO) -Wl,--as-needed
-else
-TBBMALLOC_LIBS = -Wl,--no-as-needed -ltbbmalloc -Wl,--as-needed
-endif
-ifeq ($(USE_TBBMALLOC),1)
-LDLIBS := $(TBBMALLOC_LIBS) $(LDLIBS)
-endif
+## Using sequential MKL (OpenMP disabled for better single-threaded performance)
+## Added -lifcore for Intel Fortran runtime and -limf for Intel math library
+LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl

+## Aggressive optimization flags:
+## -O3: Maximum optimization
+## -xHost: Optimize for the host CPU architecture (Intel/AMD compatible)
+## -fp-model fast=2: Aggressive floating-point optimizations
+## -fma: Enable fused multiply-add instructions
+## Note: OpenMP has been disabled (-qopenmp removed) due to performance issues
+CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma \
+               -Dfortran3 -Dnewc -I${MKLROOT}/include
+f90appflags  = -O3 -xHost -fp-model fast=2 -fma \
+               -fpp -I${MKLROOT}/include
 f90          = ifx
 f77          = ifx
 CXX          = icpx
 CC           = icx
-CLINKER      = mpiicpx
+CLINKER      = mpiicpx 

 Cu = nvcc
 CUDA_LIB_PATH = -L/usr/lib/cuda/lib64 -I/usr/include -I/usr/lib/cuda/include
 #CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -arch compute_13 -code compute_13,sm_13 -Dfortran3 -Dnewc
 CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -Dfortran3 -Dnewc
+
--- a/makefile_and_run.py
+++ b/makefile_and_run.py
@@ -28,7 +28,7 @@ def makefile_ABE():

    ## Build command
    if (input_data.GPU_Calculation == "no"):
-        makefile_command  = "make -j96" + " ABE"
+        makefile_command  = "make -j4" + " ABE"
    elif (input_data.GPU_Calculation == "yes"):
        makefile_command  = "make -j4" + " ABEGPU"
    else:
--- a/parallel_plot_helper.py
+++ b/parallel_plot_helper.py
@@ -1,12 +0,0 @@
-import multiprocessing
-
-
-def run_plot_task(task):
-    func, args = task
-    return func(*args)
-
-
-def run_plot_tasks_parallel(plot_tasks):
-    ctx = multiprocessing.get_context('fork')
-    with ctx.Pool() as pool:
-        pool.map(run_plot_task, plot_tasks)
--- a/plot_GW_strain_amplitude_xiaoqu.py
+++ b/plot_GW_strain_amplitude_xiaoqu.py
@@ -8,13 +8,11 @@
 ##
 #################################################

-import numpy                               ## numpy for array operations
-import scipy                               ## scipy for interpolation and signal processing
-import math
-import matplotlib
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
-import matplotlib.pyplot    as     plt     ## matplotlib for plotting
-import os                                  ## os for system/file operations
+import numpy                               ## numpy for array operations
+import scipy                               ## scipy for interpolation and signal processing
+import math
+import matplotlib.pyplot    as     plt     ## matplotlib for plotting
+import os                                  ## os for system/file operations

 import AMSS_NCKU_Input as input_data

--- a/plot_binary_data.py
+++ b/plot_binary_data.py
@@ -6,22 +6,17 @@
 ## Author: Xiaoqu
 ## Dates: 2024/10/01 --- 2025/09/14
 ##
-#################################################
-
-## Restrict OpenMP to one thread per process so that parallel
-## subprocess plotting does not multiply BLAS thread counts.
-import os
-os.environ.setdefault("OMP_NUM_THREADS", "1")
-
-import numpy
-import scipy
-import matplotlib
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
-import matplotlib.pyplot    as     plt
-from   matplotlib.colors    import LogNorm
-from   mpl_toolkits.mplot3d import Axes3D
-## import torch
-import AMSS_NCKU_Input      as input_data
+#################################################
+
+import numpy
+import scipy
+import matplotlib.pyplot    as     plt
+from   matplotlib.colors    import LogNorm
+from   mpl_toolkits.mplot3d import Axes3D
+## import torch
+import AMSS_NCKU_Input      as input_data
+
+import os


 #########################################################################################
@@ -97,9 +92,9 @@ def plot_binary_data( filename, binary_outdir, figure_outdir ):



-####################################################################################
-
-# Plot a single binary dataset (2D slices and 3D surface)
+####################################################################################
+
+# Plot a single binary dataset (2D slices and 3D surface)

 def get_data_xy( Rmin, Rmax, n, data0, time, figure_title, figure_outdir ):

@@ -193,15 +188,7 @@ def get_data_xy( Rmin, Rmax, n, data0, time, figure_title, figure_outdir ):
    plt.savefig( os.path.join(figure_surfaceplot_outdir, figure_title + " time = " + str(time) + " surface_plot.pdf") )   # save figure
    plt.close()

-    return
-
-####################################################################################
-
-## Allow standalone subprocess execution for parallel binary-data plotting.
-if __name__ == '__main__':
-    import sys
-    if len(sys.argv) != 4:
-        print(f"Usage: {sys.argv[0]} <filename> <binary_outdir> <figure_outdir>")
-        sys.exit(1)
-    plot_binary_data(sys.argv[1], sys.argv[2], sys.argv[3])
+    return
+
+####################################################################################

--- a/plot_xiaoqu.py
+++ b/plot_xiaoqu.py
@@ -6,20 +6,15 @@
 ## 2024/10/01 --- 2025/09/14
 ##
 #################################################
-
-import numpy                               ## numpy for array operations
-import matplotlib
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
-import matplotlib.pyplot    as     plt     ## matplotlib for plotting
-from   mpl_toolkits.mplot3d import Axes3D  ## needed for 3D plots
-import glob
-import os                                  ## operating system utilities
-
-import plot_binary_data
-import AMSS_NCKU_Input as input_data
-import subprocess
-import sys
-import multiprocessing
+
+import numpy                               ## numpy for array operations
+import matplotlib.pyplot    as     plt     ## matplotlib for plotting
+from   mpl_toolkits.mplot3d import Axes3D  ## needed for 3D plots
+import glob
+import os                                  ## operating system utilities
+
+import plot_binary_data
+import AMSS_NCKU_Input as input_data

 # plt.rcParams['text.usetex'] = True  ## enable LaTeX fonts in plots

@@ -55,37 +50,13 @@ def generate_binary_data_plot( binary_outdir, figure_outdir ):
        file_list.append(x)
        print(x)

-    ## Plot each file in parallel using subprocesses.
-    ## Each subprocess starts with BLAS thread limits in plot_binary_data.py.
-    script = os.path.join( os.path.dirname(__file__), "plot_binary_data.py" )
-    max_workers = min( multiprocessing.cpu_count(), len(file_list) ) if file_list else 0
-
-    running = []
-    failed  = []
-    for filename in file_list:
-        print(filename)
-        proc = subprocess.Popen(
-            [sys.executable, script, filename, binary_outdir, figure_outdir],
-        )
-        running.append( (proc, filename) )
-        if len(running) >= max_workers:
-            p, fn = running.pop(0)
-            p.wait()
-            if p.returncode != 0:
-                failed.append(fn)
-
-    for p, fn in running:
-        p.wait()
-        if p.returncode != 0:
-            failed.append(fn)
-
-    if failed:
-        print( " WARNING: the following binary data plots failed:" )
-        for fn in failed:
-            print( "   ", fn )
-
-    print(                        )
-    print( " Binary Data Plot Has been Finished " )
+    ## Plot each file in the list
+    for filename in file_list:
+        print(filename)
+        plot_binary_data.plot_binary_data(filename, binary_outdir, figure_outdir)
+
+    print(                        )
+    print( " Binary Data Plot Has been Finished " )
    print(                                        )

    return
Author	SHA1	Message	Date
jaunatisblue	3f7e20f702	删除diff_new.f90中冗余部分，方便后续工作	2026-02-08 00:54:23 +08:00
jaunatisblue	673dd20722	对fmisc.f90的polint修改	2026-02-07 01:56:44 +08:00