Switch legacy build to GCC and OpenMPI

.idea/vcs.xml

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

AMSS_NCKU_Program.py

@@ -126,6 +126,11 @@ 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(                                                                                   )
@@ -307,7 +312,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
     
     import generate_TwoPuncture_input
     
-    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input(numerical_grid.puncture_data)
+    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
     
     print(                                                                                              )
     print( " The input parfile for the TwoPunctureABE executable has been generated. " )
@@ -349,7 +354,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
 
 import renew_puncture_parameter
     
-renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory, numerical_grid.puncture_data)
+renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory)
 
 
 ## Generated AMSS-NCKU input filename

AMSS_NCKU_Verify_ASC26.py

@@ -9,11 +9,6 @@ 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
@@ -28,10 +23,6 @@ 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:
@@ -70,132 +61,6 @@ 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.
@@ -319,45 +184,18 @@ def print_constraint_results(results, threshold=2.0):
     return passed
 
 
-def print_figure_results(results, threshold_percent=0.001):
+def print_summary(rms_passed, constraint_passed):
-    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 and figure_passed
+    all_passed = rms_passed and constraint_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}")
@@ -374,8 +212,6 @@ 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")
@@ -394,8 +230,6 @@ 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)
@@ -409,13 +243,7 @@ def main():
     constraint_results = analyze_constraint_violation(constraint_data)
     constraint_passed = print_constraint_results(constraint_results)
 
-    try:
+    all_passed = print_summary(rms_passed, constraint_passed)
-        figure_results = compare_required_figures(reference_figure_dir, target_figure_dir)
-        figure_passed = print_figure_results(figure_results)
-    except (FileNotFoundError, RuntimeError) as exc:
-        figure_passed = print_figure_error(str(exc))
-
-    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

@@ -1,87 +0,0 @@
-
-
-#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

AMSS_NCKU_source/FFT.f90

@@ -0,0 +1,92 @@
+
+
+#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
+
+SUBROUTINE four1(dataa,nn,isign)
+implicit none
+INTEGER::isign,nn
+double precision,dimension(2*nn)::dataa
+INTEGER::i,istep,j,m,mmax,n
+double precision::tempi,tempr
+DOUBLE PRECISION::theta,wi,wpi,wpr,wr,wtemp
+n=2*nn
+j=1
+do i=1,n,2
+  if(j.gt.i)then
+     tempr=dataa(j)
+     tempi=dataa(j+1)
+     dataa(j)=dataa(i)
+     dataa(j+1)=dataa(i+1)
+     dataa(i)=tempr
+     dataa(i+1)=tempi
+  endif
+  m=nn
+1 if ((m.ge.2).and.(j.gt.m)) then
+  j=j-m
+  m=m/2
+goto 1
+  endif
+j=j+m
+enddo
+mmax=2
+2  if (n.gt.mmax) then
+     istep=2*mmax
+     theta=6.28318530717959d0/(isign*mmax)
+     wpr=-2.d0*sin(0.5d0*theta)**2
+     wpi=sin(theta)
+     wr=1.d0
+     wi=0.d0
+     do m=1,mmax,2
+       do i=m,n,istep
+         j=i+mmax
+         tempr=sngl(wr)*dataa(j)-sngl(wi)*dataa(j+1)
+         tempi=sngl(wr)*dataa(j+1)+sngl(wi)*dataa(j)
+         dataa(j)=dataa(i)-tempr
+         dataa(j+1)=dataa(i+1)-tempi
+         dataa(i)=dataa(i)+tempr
+         dataa(i+1)=dataa(i+1)+tempi
+       enddo
+          wtemp=wr
+          wr=wr*wpr-wi*wpi+wr
+          wi=wi*wpr+wtemp*wpi+wi
+     enddo
+mmax=istep
+goto 2
+endif
+return
+END SUBROUTINE four1

AMSS_NCKU_source/Surface_Integral/gaussj.C

@@ -1,106 +0,0 @@
-
-#ifdef newc
-#include <iostream>
-#include <iomanip>
-#include <fstream>
-#include <cstdlib>
-#include <cstring>
-#include <cmath>
-using namespace std;
-#else
-#include <iostream.h>
-#include <iomanip.h>
-#include <fstream.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <math.h>
-#endif
-
-// Intel oneMKL LAPACK interface
-#include <mkl_lapacke.h>
-/* Linear equation solution using Intel oneMKL LAPACK.
-a[0..n-1][0..n-1] is the input matrix. b[0..n-1] is input
-containing the right-hand side vectors. On output a is
-replaced by its matrix inverse, and b is replaced by the
-corresponding set of solution vectors.
-
-Mathematical equivalence:
-  Solves: A * x = b  =>  x = A^(-1) * b
-  Original Gauss-Jordan and LAPACK dgesv/dgetri produce identical results
-  within numerical precision. */
-
-int gaussj(double *a, double *b, int n)
-{
-  // Allocate pivot array and workspace
-  lapack_int *ipiv = new lapack_int[n];
-  lapack_int info;
-
-  // Make a copy of matrix a for solving (dgesv modifies it to LU form)
-  double *a_copy = new double[n * n];
-  for (int i = 0; i < n * n; i++) {
-    a_copy[i] = a[i];
-  }
-
-  // Step 1: Solve linear system A*x = b using LU decomposition
-  // LAPACKE_dgesv uses column-major by default, but we use row-major
-  info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, a_copy, n, ipiv, b, 1);
-
-  if (info != 0) {
-    cout << "gaussj: Singular Matrix (dgesv info=" << info << ")" << endl;
-    delete[] ipiv;
-    delete[] a_copy;
-    return 1;
-  }
-
-  // Step 2: Compute matrix inverse A^(-1) using LU factorization
-  // First do LU factorization of original matrix a
-  info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, a, n, ipiv);
-
-  if (info != 0) {
-    cout << "gaussj: Singular Matrix (dgetrf info=" << info << ")" << endl;
-    delete[] ipiv;
-    delete[] a_copy;
-    return 1;
-  }
-
-  // Then compute inverse from LU factorization
-  info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, a, n, ipiv);
-
-  if (info != 0) {
-    cout << "gaussj: Singular Matrix (dgetri info=" << info << ")" << endl;
-    delete[] ipiv;
-    delete[] a_copy;
-    return 1;
-  }
-
-  delete[] ipiv;
-  delete[] a_copy;
-
-  return 0;
-}
-// for check usage
-/*
-int main()
-{
-  double *A,*b;
-  A=new double[9];
-  b=new double[3];
-
-  A[0]=0.5; A[1]=1.0/3; A[2]=1;
-  A[3]=1;   A[4]=5.0/3; A[5]=3;
-  A[6]=2;   A[7]=4.0/3; A[8]=5;
-
-  b[0]=1; b[1]=3; b[2]=2;
-
-  cout<<"initial data:"<<endl;
-  for(int i=0;i<3;i++) cout<<A[i*3]<<" "<<A[i*3+1]<<" "<<A[i*3+2]<<" "<<b[i]<<endl;
-
-  gaussj(A, b, 3);
-
-  cout<<"final data:"<<endl;
-  for(int i=0;i<3;i++) cout<<A[i*3]<<" "<<A[i*3+1]<<" "<<A[i*3+2]<<" "<<b[i]<<endl;
-
-  delete[] A; delete[] b;
-}
-*/

AMSS_NCKU_source/TwoPunctures.C

@@ -27,7 +27,21 @@ using namespace std;
 #endif
 
 #include "TwoPunctures.h"
-#include <mkl_cblas.h>
+
+extern "C" {
+double cblas_ddot(const int, const double *, const int, const double *, const int);
+double cblas_dnrm2(const int, const double *, const int);
+void cblas_dgemm(const int, const int, const int,
+                 const int, const int, const int,
+                 const double, const double *, const int,
+                 const double *, const int, const double,
+                 double *, const int);
+}
+
+enum {
+  CblasRowMajor = 101,
+  CblasNoTrans = 111
+};
 
 TwoPunctures::TwoPunctures(double mp, double mm, double b,
                            double P_plusx, double P_plusy, double P_plusz,