Compare commits
1 Commits
cjy-oneapi
...
baseline
| Author | SHA1 | Date | |
|---|---|---|---|
|
79af79d471 |
4
.gitignore
vendored
4
.gitignore
vendored
@@ -1,6 +1,2 @@
|
||||
__pycache__
|
||||
GW150914
|
||||
GW150914-origin
|
||||
docs
|
||||
*.tmp
|
||||
|
||||
|
||||
@@ -1,279 +0,0 @@
|
||||
#!/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()
|
||||
@@ -37,51 +37,57 @@ 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
|
||||
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
|
||||
|
||||
! 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
|
||||
|
||||
@@ -27,7 +27,6 @@ using namespace std;
|
||||
#endif
|
||||
|
||||
#include "TwoPunctures.h"
|
||||
#include <mkl_cblas.h>
|
||||
|
||||
TwoPunctures::TwoPunctures(double mp, double mm, double b,
|
||||
double P_plusx, double P_plusy, double P_plusz,
|
||||
@@ -892,17 +891,25 @@ double TwoPunctures::norm1(double *v, int n)
|
||||
/* -------------------------------------------------------------------------*/
|
||||
double TwoPunctures::norm2(double *v, int n)
|
||||
{
|
||||
// Optimized with oneMKL BLAS DNRM2
|
||||
// Computes: sqrt(sum(v[i]^2))
|
||||
return cblas_dnrm2(n, v, 1);
|
||||
int i;
|
||||
double result = 0;
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
result += v[i] * v[i];
|
||||
|
||||
return sqrt(result);
|
||||
}
|
||||
|
||||
/* -------------------------------------------------------------------------*/
|
||||
double TwoPunctures::scalarproduct(double *v, double *w, int n)
|
||||
{
|
||||
// Optimized with oneMKL BLAS DDOT
|
||||
// Computes: sum(v[i] * w[i])
|
||||
return cblas_ddot(n, v, 1, w, 1);
|
||||
int i;
|
||||
double result = 0;
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
result += v[i] * w[i];
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
/* -------------------------------------------------------------------------*/
|
||||
|
||||
@@ -1276,9 +1276,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)
|
||||
@@ -1302,12 +1300,7 @@ if(dabs(X(1)-xmin) < dX) imin = 1
|
||||
if(dabs(Y(1)-ymin) < dY) jmin = 1
|
||||
if(dabs(Z(1)-zmin) < dZ) kmin = 1
|
||||
|
||||
! Optimized with oneMKL BLAS DDOT for dot product
|
||||
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
|
||||
|
||||
@@ -1332,9 +1325,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,n_elements
|
||||
real*8, dimension(:), allocatable :: f_flat
|
||||
real*8, external :: DDOT
|
||||
integer::i,j,k
|
||||
|
||||
real*8 :: PIo4
|
||||
|
||||
@@ -1397,12 +1388,7 @@ if(Symmetry==2)then
|
||||
if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
|
||||
endif
|
||||
|
||||
! Optimized with oneMKL BLAS DDOT for dot product
|
||||
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
|
||||
|
||||
@@ -1430,8 +1416,6 @@ f_out = f_out*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
|
||||
|
||||
real*8 :: PIo4
|
||||
|
||||
@@ -1494,12 +1478,11 @@ if(Symmetry==2)then
|
||||
if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
|
||||
endif
|
||||
|
||||
! Optimized with oneMKL BLAS DDOT for dot product
|
||||
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)
|
||||
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)
|
||||
|
||||
return
|
||||
|
||||
@@ -1697,7 +1680,6 @@ deallocate(f_flat)
|
||||
real*8, dimension(ORDN,ORDN) :: tmp2
|
||||
real*8, dimension(ORDN) :: tmp1
|
||||
real*8, dimension(3) :: SoAh
|
||||
real*8, external :: DDOT
|
||||
|
||||
! +1 because c++ gives 0 for first point
|
||||
cxB = inds+1
|
||||
@@ -1733,21 +1715,20 @@ deallocate(f_flat)
|
||||
ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),cxB(3):cxT(3))
|
||||
endif
|
||||
|
||||
! Optimized with BLAS operations for better performance
|
||||
! First dimension: z-direction weighted sum
|
||||
tmp2=0
|
||||
do m=1,ORDN
|
||||
tmp2 = tmp2 + coef(2*ORDN+m)*ya(:,:,m)
|
||||
enddo
|
||||
|
||||
! Second dimension: y-direction weighted sum
|
||||
tmp1=0
|
||||
do m=1,ORDN
|
||||
tmp1 = tmp1 + coef(ORDN+m)*tmp2(:,m)
|
||||
enddo
|
||||
|
||||
! Third dimension: x-direction weighted sum using BLAS DDOT
|
||||
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
|
||||
|
||||
@@ -1777,7 +1758,6 @@ deallocate(f_flat)
|
||||
real*8, dimension(ORDN,ORDN) :: ya
|
||||
real*8, dimension(ORDN) :: tmp1
|
||||
real*8, dimension(2) :: SoAh
|
||||
real*8, external :: DDOT
|
||||
|
||||
! +1 because c++ gives 0 for first point
|
||||
cxB = inds(1:2)+1
|
||||
@@ -1807,14 +1787,15 @@ deallocate(f_flat)
|
||||
ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),inds(3))
|
||||
endif
|
||||
|
||||
! Optimized with BLAS operations
|
||||
tmp1=0
|
||||
do m=1,ORDN
|
||||
tmp1 = tmp1 + coef(ORDN+m)*ya(:,m)
|
||||
enddo
|
||||
|
||||
! Use BLAS DDOT for final weighted sum
|
||||
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
|
||||
|
||||
@@ -1845,7 +1826,6 @@ deallocate(f_flat)
|
||||
real*8, dimension(ORDN) :: ya
|
||||
real*8 :: SoAh
|
||||
integer,dimension(3) :: inds
|
||||
real*8, external :: DDOT
|
||||
|
||||
! +1 because c++ gives 0 for first point
|
||||
inds = indsi + 1
|
||||
@@ -1906,8 +1886,10 @@ deallocate(f_flat)
|
||||
write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
|
||||
endif
|
||||
|
||||
! Optimized with BLAS DDOT for weighted sum
|
||||
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
|
||||
|
||||
@@ -2139,38 +2121,24 @@ deallocate(f_flat)
|
||||
|
||||
end function fWigner_d_function
|
||||
!----------------------------------
|
||||
! Optimized factorial function using lookup table for small N
|
||||
! and log-gamma for large N to avoid overflow
|
||||
function ffact(N) result(gont)
|
||||
implicit none
|
||||
integer,intent(in) :: N
|
||||
|
||||
real*8 :: gont
|
||||
integer :: i
|
||||
|
||||
! Lookup table for factorials 0! to 20! (precomputed)
|
||||
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 ]
|
||||
integer :: i
|
||||
|
||||
! sanity check
|
||||
if(N < 0)then
|
||||
write(*,*) "ffact: error input for factorial"
|
||||
gont = 1.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! Use lookup table for small N (fast path)
|
||||
if(N <= 20)then
|
||||
gont = fact_table(N)
|
||||
else
|
||||
! Use log-gamma function for large N: N! = exp(log_gamma(N+1))
|
||||
! This avoids overflow and is computed efficiently
|
||||
gont = exp(log_gamma(dble(N+1)))
|
||||
endif
|
||||
gont = 1.d0
|
||||
do i=1,N
|
||||
gont = gont*i
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
|
||||
@@ -16,66 +16,115 @@ using namespace std;
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#endif
|
||||
|
||||
// Intel oneMKL LAPACK interface
|
||||
#include <mkl_lapacke.h>
|
||||
/* Linear equation solution using Intel oneMKL LAPACK.
|
||||
/* Linear equation solution by Gauss-Jordan elimination.
|
||||
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. */
|
||||
corresponding set of solution vectors */
|
||||
|
||||
int gaussj(double *a, double *b, int n)
|
||||
{
|
||||
// Allocate pivot array and workspace
|
||||
lapack_int *ipiv = new lapack_int[n];
|
||||
lapack_int info;
|
||||
double swap;
|
||||
|
||||
// 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];
|
||||
int *indxc, *indxr, *ipiv;
|
||||
indxc = new int[n];
|
||||
indxr = new int[n];
|
||||
ipiv = new int[n];
|
||||
|
||||
int i, icol, irow, j, k, l, ll;
|
||||
double big, dum, pivinv, temp;
|
||||
|
||||
for (j = 0; j < n; j++)
|
||||
ipiv[j] = 0;
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
big = 0.0;
|
||||
for (j = 0; j < n; j++)
|
||||
if (ipiv[j] != 1)
|
||||
for (k = 0; k < n; k++)
|
||||
{
|
||||
if (ipiv[k] == 0)
|
||||
{
|
||||
if (fabs(a[j * n + k]) >= big)
|
||||
{
|
||||
big = fabs(a[j * n + k]);
|
||||
irow = j;
|
||||
icol = k;
|
||||
}
|
||||
}
|
||||
else if (ipiv[k] > 1)
|
||||
{
|
||||
cout << "gaussj: Singular Matrix-1" << endl;
|
||||
for (int ii = 0; ii < n; ii++)
|
||||
{
|
||||
for (int jj = 0; jj < n; jj++)
|
||||
cout << a[ii * n + jj] << " ";
|
||||
cout << endl;
|
||||
}
|
||||
return 1; // error return
|
||||
}
|
||||
}
|
||||
|
||||
ipiv[icol] = ipiv[icol] + 1;
|
||||
if (irow != icol)
|
||||
{
|
||||
for (l = 0; l < n; l++)
|
||||
{
|
||||
swap = a[irow * n + l];
|
||||
a[irow * n + l] = a[icol * n + l];
|
||||
a[icol * n + l] = swap;
|
||||
}
|
||||
|
||||
swap = b[irow];
|
||||
b[irow] = b[icol];
|
||||
b[icol] = swap;
|
||||
}
|
||||
|
||||
indxr[i] = irow;
|
||||
indxc[i] = icol;
|
||||
|
||||
if (a[icol * n + icol] == 0.0)
|
||||
{
|
||||
cout << "gaussj: Singular Matrix-2" << endl;
|
||||
for (int ii = 0; ii < n; ii++)
|
||||
{
|
||||
for (int jj = 0; jj < n; jj++)
|
||||
cout << a[ii * n + jj] << " ";
|
||||
cout << endl;
|
||||
}
|
||||
return 1; // error return
|
||||
}
|
||||
|
||||
pivinv = 1.0 / a[icol * n + icol];
|
||||
a[icol * n + icol] = 1.0;
|
||||
for (l = 0; l < n; l++)
|
||||
a[icol * n + l] *= pivinv;
|
||||
b[icol] *= pivinv;
|
||||
for (ll = 0; ll < n; ll++)
|
||||
if (ll != icol)
|
||||
{
|
||||
dum = a[ll * n + icol];
|
||||
a[ll * n + icol] = 0.0;
|
||||
for (l = 0; l < n; l++)
|
||||
a[ll * n + l] -= a[icol * n + l] * dum;
|
||||
b[ll] -= b[icol] * dum;
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
for (l = n - 1; l >= 0; l--)
|
||||
{
|
||||
if (indxr[l] != indxc[l])
|
||||
for (k = 0; k < n; k++)
|
||||
{
|
||||
swap = a[k * n + indxr[l]];
|
||||
a[k * n + indxr[l]] = a[k * n + indxc[l]];
|
||||
a[k * n + indxc[l]] = swap;
|
||||
}
|
||||
}
|
||||
|
||||
delete[] indxc;
|
||||
delete[] indxr;
|
||||
delete[] ipiv;
|
||||
delete[] a_copy;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
@@ -512,10 +512,11 @@
|
||||
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
||||
DIMENSION V(N),W(N)
|
||||
! SUBROUTINE TO COMPUTE DOUBLE PRECISION VECTOR DOT PRODUCT.
|
||||
! Optimized using Intel oneMKL BLAS ddot
|
||||
! Mathematical equivalence: DGVV = sum_{i=1}^{N} V(i)*W(i)
|
||||
|
||||
DOUBLE PRECISION, EXTERNAL :: DDOT
|
||||
DGVV = DDOT(N, V, 1, W, 1)
|
||||
SUM = 0.0D0
|
||||
DO 10 I = 1,N
|
||||
SUM = SUM + V(I)*W(I)
|
||||
10 CONTINUE
|
||||
DGVV = SUM
|
||||
RETURN
|
||||
END
|
||||
|
||||
@@ -2,7 +2,7 @@
|
||||
#ifndef MICRODEF_H
|
||||
#define MICRODEF_H
|
||||
|
||||
#include "macrodef.fh"
|
||||
#include "microdef.fh"
|
||||
|
||||
// application parameters
|
||||
|
||||
|
||||
@@ -1,30 +1,19 @@
|
||||
## 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 (Optimized for performance)
|
||||
filein = -I/usr/include/ -I${MKLROOT}/include
|
||||
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/
|
||||
|
||||
## 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
|
||||
## LDLIBS = -L/usr/lib/x86_64-linux-gnu -lmpich -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran
|
||||
LDLIBS = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran -lmpi -lgfortran
|
||||
|
||||
## 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
|
||||
CXXAPPFLAGS = -O0 -Wno-deprecated -Dfortran3 -Dnewc
|
||||
#f90appflags = -O0 -fpp
|
||||
f90appflags = -O0 -x f95-cpp-input
|
||||
f90 = gfortran
|
||||
f77 = gfortran
|
||||
CXX = g++
|
||||
CC = gcc
|
||||
CLINKER = mpic++
|
||||
|
||||
Cu = nvcc
|
||||
CUDA_LIB_PATH = -L/usr/lib/cuda/lib64 -I/usr/include -I/usr/lib/cuda/include
|
||||
|
||||
@@ -11,17 +11,6 @@
|
||||
import AMSS_NCKU_Input as input_data
|
||||
import subprocess
|
||||
|
||||
## CPU core binding configuration using taskset
|
||||
## taskset ensures all child processes inherit the CPU affinity mask
|
||||
## This forces make and all compiler processes to use only nohz_full cores (4-55, 60-111)
|
||||
## Format: taskset -c 4-55,60-111 ensures processes only run on these cores
|
||||
NUMACTL_CPU_BIND = "taskset -c 0-111"
|
||||
|
||||
## Build parallelism configuration
|
||||
## Use nohz_full cores (4-55, 60-111) for compilation: 52 + 52 = 104 cores
|
||||
## Set make -j to utilize available cores for faster builds
|
||||
BUILD_JOBS = 104
|
||||
|
||||
|
||||
##################################################################
|
||||
|
||||
@@ -37,11 +26,11 @@ def makefile_ABE():
|
||||
print( " Compiling the AMSS-NCKU executable file ABE/ABEGPU " )
|
||||
print( )
|
||||
|
||||
## Build command with CPU binding to nohz_full cores
|
||||
## Build command
|
||||
if (input_data.GPU_Calculation == "no"):
|
||||
makefile_command = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} ABE"
|
||||
makefile_command = "make -j96" + " ABE"
|
||||
elif (input_data.GPU_Calculation == "yes"):
|
||||
makefile_command = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} ABEGPU"
|
||||
makefile_command = "make -j4" + " ABEGPU"
|
||||
else:
|
||||
print( " CPU/GPU numerical calculation setting is wrong " )
|
||||
print( )
|
||||
@@ -78,8 +67,8 @@ def makefile_TwoPunctureABE():
|
||||
print( " Compiling the AMSS-NCKU executable file TwoPunctureABE " )
|
||||
print( )
|
||||
|
||||
## Build command with CPU binding to nohz_full cores
|
||||
makefile_command = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} TwoPunctureABE"
|
||||
## Build command
|
||||
makefile_command = "make" + " TwoPunctureABE"
|
||||
|
||||
## Execute the command with subprocess.Popen and stream output
|
||||
makefile_process = subprocess.Popen(makefile_command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)
|
||||
@@ -116,10 +105,10 @@ def run_ABE():
|
||||
## Define the command to run; cast other values to strings as needed
|
||||
|
||||
if (input_data.GPU_Calculation == "no"):
|
||||
mpi_command = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
|
||||
mpi_command = "mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
|
||||
mpi_command_outfile = "ABE_out.log"
|
||||
elif (input_data.GPU_Calculation == "yes"):
|
||||
mpi_command = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + " ./ABEGPU"
|
||||
mpi_command = "mpirun -np " + str(input_data.MPI_processes) + " ./ABEGPU"
|
||||
mpi_command_outfile = "ABEGPU_out.log"
|
||||
|
||||
## Execute the MPI command and stream output
|
||||
@@ -158,7 +147,7 @@ def run_TwoPunctureABE():
|
||||
print( )
|
||||
|
||||
## Define the command to run
|
||||
TwoPuncture_command = NUMACTL_CPU_BIND + " ./TwoPunctureABE"
|
||||
TwoPuncture_command = "./TwoPunctureABE"
|
||||
TwoPuncture_command_outfile = "TwoPunctureABE_out.log"
|
||||
|
||||
## Execute the command with subprocess.Popen and stream output
|
||||
|
||||
Reference in New Issue
Block a user