[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_Program.py

@@ -126,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(                                                                                   )
@@ -312,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. " )
@@ -354,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

AMSS_NCKU_Verify_ASC26.py

@@ -9,6 +9,11 @@ Verification Requirements:
    - 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
@@ -23,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:
@@ -61,6 +70,132 @@ def load_constraint_data(filepath):
                 data.append([float(x) for x in parts[:8]])
     return np.array(data)
 
+
+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 3D Vector RMS and component-wise RMS (X, Y, Z) independently.
@@ -184,18 +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(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] 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}")
@@ -212,6 +374,8 @@ def main():
 
     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")
 
     bh_file_ref = os.path.join(reference_dir, "bssn_BH.dat")
     bh_file_target = os.path.join(target_dir, "bssn_BH.dat")
@@ -230,6 +394,8 @@ def main():
     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}")
 
     bh_data_ref = load_bh_trajectory(bh_file_ref)
     bh_data_target = load_bh_trajectory(bh_file_target)
@@ -243,7 +409,13 @@ def main():
     constraint_results = analyze_constraint_violation(constraint_data)
     constraint_passed = print_constraint_results(constraint_results)
 
-    all_passed = print_summary(rms_passed, constraint_passed)
+    try:
+        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))
+
+    all_passed = print_summary(rms_passed, constraint_passed, figure_passed)
     sys.exit(0 if all_passed else 1)
 
 if __name__ == "__main__":

AMSS_NCKU_source/AHF_Direct/FFT.f90

@@ -0,0 +1,87 @@
+
+
+#if 0
+program checkFFT
+use dfport
+implicit none
+double precision::x
+integer,parameter::N=256
+double precision,dimension(N*2)::p
+double precision,dimension(N/2)::s
+integer::ncount,j,idum
+character(len=8)::tt
+tt=clock()
+idum=iachar(tt(8:8))-48
+p=0.0
+open(77,file='prime.dat',status='unknown')
+loop1:do ncount=1,N
+   x=ran(idum)
+   p(2*ncount-1)=x
+   write(77,'(f15.3)')x
+enddo loop1
+close(77)
+call four1(p,N,1)
+do j=1,N/2
+  s(j)=p(2*j)*p(2*j)+p(2*j-1)*p(2*j-1)
+enddo
+x=0.0
+do j=1,N/2
+  x=x+s(j)
+enddo
+s=s/x
+open(77,file='power.dat',status='unknown')
+do j=1,N/2
+  write(77,'(2(1x,f15.3))')dble(j-1)/dble(N),s(j)
+enddo
+close(77)
+end program checkFFT
+#endif
+
+!-------------
+! Optimized FFT using Intel oneMKL DFTI
+! Mathematical equivalence: Standard DFT definition
+!   Forward (isign=1):  X[k] = sum_{n=0}^{N-1} x[n] * exp(-2*pi*i*k*n/N)
+!   Backward (isign=-1): X[k] = sum_{n=0}^{N-1} x[n] * exp(+2*pi*i*k*n/N)
+! Input/Output: dataa is interleaved complex array [Re(0),Im(0),Re(1),Im(1),...]
+!-------------
+SUBROUTINE four1(dataa,nn,isign)
+use MKL_DFTI
+implicit none
+INTEGER, intent(in) :: isign, nn
+DOUBLE PRECISION, dimension(2*nn), intent(inout) :: dataa
+
+type(DFTI_DESCRIPTOR), pointer :: desc
+integer :: status
+
+! Create DFTI descriptor for 1D complex-to-complex transform
+status = DftiCreateDescriptor(desc, DFTI_DOUBLE, DFTI_COMPLEX, 1, nn)
+if (status /= 0) return
+
+! Set input/output storage as interleaved complex (default)
+status = DftiSetValue(desc, DFTI_PLACEMENT, DFTI_INPLACE)
+if (status /= 0) then
+   status = DftiFreeDescriptor(desc)
+   return
+endif
+
+! Commit the descriptor
+status = DftiCommitDescriptor(desc)
+if (status /= 0) then
+   status = DftiFreeDescriptor(desc)
+   return
+endif
+
+! Execute FFT based on direction
+if (isign == 1) then
+   ! Forward FFT: exp(-2*pi*i*k*n/N)
+   status = DftiComputeForward(desc, dataa)
+else
+   ! Backward FFT: exp(+2*pi*i*k*n/N)
+   status = DftiComputeBackward(desc, dataa)
+endif
+
+! Free descriptor
+status = DftiFreeDescriptor(desc)
+
+return
+END SUBROUTINE four1