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
|
||||
|
||||
! Commit the descriptor
|
||||
status = DftiCommitDescriptor(desc)
|
||||
if (status /= 0) then
|
||||
status = DftiFreeDescriptor(desc)
|
||||
return
|
||||
m=nn
|
||||
1 if ((m.ge.2).and.(j.gt.m)) then
|
||||
j=j-m
|
||||
m=m/2
|
||||
goto 1
|
||||
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)
|
||||
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
|
||||
|
||||
! Free descriptor
|
||||
status = DftiFreeDescriptor(desc)
|
||||
|
||||
return
|
||||
END SUBROUTINE four1
|
||||
|
||||
@@ -5,7 +5,6 @@
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <string>
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <fstream>
|
||||
@@ -28,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,
|
||||
@@ -61,110 +59,13 @@ TwoPunctures::TwoPunctures(double mp, double mm, double b,
|
||||
F = dvector(0, ntotal - 1);
|
||||
allocate_derivs(&u, ntotal);
|
||||
allocate_derivs(&v, ntotal);
|
||||
|
||||
// Allocate workspace buffers for hot-path allocation elimination
|
||||
int N = maximum3(n1, n2, n3);
|
||||
int maxn = maximum2(n1, n2);
|
||||
|
||||
// LineRelax_be workspace (sized for n2)
|
||||
ws_diag_be = new double[n2];
|
||||
ws_e_be = new double[n2 - 1];
|
||||
ws_f_be = new double[n2 - 1];
|
||||
ws_b_be = new double[n2];
|
||||
ws_x_be = new double[n2];
|
||||
|
||||
// LineRelax_al workspace (sized for n1)
|
||||
ws_diag_al = new double[n1];
|
||||
ws_e_al = new double[n1 - 1];
|
||||
ws_f_al = new double[n1 - 1];
|
||||
ws_b_al = new double[n1];
|
||||
ws_x_al = new double[n1];
|
||||
|
||||
// ThomasAlgorithm workspace (sized for max(n1,n2))
|
||||
ws_thomas_y = new double[maxn];
|
||||
|
||||
// JFD_times_dv workspace (sized for nvar)
|
||||
ws_jfd_values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&ws_jfd_dU, nvar);
|
||||
allocate_derivs(&ws_jfd_U, nvar);
|
||||
|
||||
// chebft_Zeros workspace (sized for N+1)
|
||||
ws_cheb_c = dvector(0, N);
|
||||
|
||||
// fourft workspace (sized for N/2+1 each)
|
||||
ws_four_a = dvector(0, N / 2);
|
||||
ws_four_b = dvector(0, N / 2);
|
||||
|
||||
// Derivatives_AB3 workspace
|
||||
ws_deriv_p = dvector(0, N);
|
||||
ws_deriv_dp = dvector(0, N);
|
||||
ws_deriv_d2p = dvector(0, N);
|
||||
ws_deriv_q = dvector(0, N);
|
||||
ws_deriv_dq = dvector(0, N);
|
||||
ws_deriv_r = dvector(0, N);
|
||||
ws_deriv_dr = dvector(0, N);
|
||||
ws_deriv_indx = ivector(0, N);
|
||||
|
||||
// F_of_v workspace
|
||||
ws_fov_sources = new double[n1 * n2 * n3];
|
||||
ws_fov_values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&ws_fov_U, nvar);
|
||||
|
||||
// J_times_dv workspace
|
||||
ws_jtdv_values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&ws_jtdv_dU, nvar);
|
||||
allocate_derivs(&ws_jtdv_U, nvar);
|
||||
}
|
||||
|
||||
TwoPunctures::~TwoPunctures()
|
||||
{
|
||||
int const nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
|
||||
int N = maximum3(n1, n2, n3);
|
||||
|
||||
free_dvector(F, 0, ntotal - 1);
|
||||
free_derivs(&u, ntotal);
|
||||
free_derivs(&v, ntotal);
|
||||
|
||||
// Free workspace buffers
|
||||
delete[] ws_diag_be;
|
||||
delete[] ws_e_be;
|
||||
delete[] ws_f_be;
|
||||
delete[] ws_b_be;
|
||||
delete[] ws_x_be;
|
||||
|
||||
delete[] ws_diag_al;
|
||||
delete[] ws_e_al;
|
||||
delete[] ws_f_al;
|
||||
delete[] ws_b_al;
|
||||
delete[] ws_x_al;
|
||||
|
||||
delete[] ws_thomas_y;
|
||||
|
||||
free_dvector(ws_jfd_values, 0, nvar - 1);
|
||||
free_derivs(&ws_jfd_dU, nvar);
|
||||
free_derivs(&ws_jfd_U, nvar);
|
||||
|
||||
free_dvector(ws_cheb_c, 0, N);
|
||||
|
||||
free_dvector(ws_four_a, 0, N / 2);
|
||||
free_dvector(ws_four_b, 0, N / 2);
|
||||
|
||||
free_dvector(ws_deriv_p, 0, N);
|
||||
free_dvector(ws_deriv_dp, 0, N);
|
||||
free_dvector(ws_deriv_d2p, 0, N);
|
||||
free_dvector(ws_deriv_q, 0, N);
|
||||
free_dvector(ws_deriv_dq, 0, N);
|
||||
free_dvector(ws_deriv_r, 0, N);
|
||||
free_dvector(ws_deriv_dr, 0, N);
|
||||
free_ivector(ws_deriv_indx, 0, N);
|
||||
|
||||
delete[] ws_fov_sources;
|
||||
free_dvector(ws_fov_values, 0, nvar - 1);
|
||||
free_derivs(&ws_fov_U, nvar);
|
||||
|
||||
free_dvector(ws_jtdv_values, 0, nvar - 1);
|
||||
free_derivs(&ws_jtdv_dU, nvar);
|
||||
free_derivs(&ws_jtdv_U, nvar);
|
||||
}
|
||||
|
||||
void TwoPunctures::Solve()
|
||||
@@ -753,7 +654,7 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
|
||||
int k, j, isignum;
|
||||
double fac, sum, Pion, *c;
|
||||
|
||||
c = ws_cheb_c;
|
||||
c = dvector(0, n);
|
||||
Pion = Pi / n;
|
||||
if (inv == 0)
|
||||
{
|
||||
@@ -784,6 +685,7 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
|
||||
}
|
||||
for (j = 0; j < n; j++)
|
||||
u[j] = c[j];
|
||||
free_dvector(c, 0, n);
|
||||
}
|
||||
|
||||
/* --------------------------------------------------------------------------*/
|
||||
@@ -871,8 +773,8 @@ void TwoPunctures::fourft(double *u, int N, int inv)
|
||||
double x, x1, fac, Pi_fac, *a, *b;
|
||||
|
||||
M = N / 2;
|
||||
a = ws_four_a;
|
||||
b = ws_four_b - 1; /* offset to match dvector(1,M) indexing */
|
||||
a = dvector(0, M);
|
||||
b = dvector(1, M); /* Actually: b=vector(1,M-1) but this is problematic if M=1*/
|
||||
fac = 1. / M;
|
||||
Pi_fac = Pi * fac;
|
||||
if (inv == 0)
|
||||
@@ -921,6 +823,8 @@ void TwoPunctures::fourft(double *u, int N, int inv)
|
||||
iy = -iy;
|
||||
}
|
||||
}
|
||||
free_dvector(a, 0, M);
|
||||
free_dvector(b, 1, M);
|
||||
}
|
||||
|
||||
/* -----------------------------------------*/
|
||||
@@ -987,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;
|
||||
}
|
||||
|
||||
/* -------------------------------------------------------------------------*/
|
||||
@@ -1213,14 +1125,14 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
|
||||
double *p, *dp, *d2p, *q, *dq, *r, *dr;
|
||||
|
||||
N = maximum3(n1, n2, n3);
|
||||
p = ws_deriv_p;
|
||||
dp = ws_deriv_dp;
|
||||
d2p = ws_deriv_d2p;
|
||||
q = ws_deriv_q;
|
||||
dq = ws_deriv_dq;
|
||||
r = ws_deriv_r;
|
||||
dr = ws_deriv_dr;
|
||||
indx = ws_deriv_indx;
|
||||
p = dvector(0, N);
|
||||
dp = dvector(0, N);
|
||||
d2p = dvector(0, N);
|
||||
q = dvector(0, N);
|
||||
dq = dvector(0, N);
|
||||
r = dvector(0, N);
|
||||
dr = dvector(0, N);
|
||||
indx = ivector(0, N);
|
||||
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
{
|
||||
@@ -1303,6 +1215,14 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
|
||||
}
|
||||
}
|
||||
}
|
||||
free_dvector(p, 0, N);
|
||||
free_dvector(dp, 0, N);
|
||||
free_dvector(d2p, 0, N);
|
||||
free_dvector(q, 0, N);
|
||||
free_dvector(dq, 0, N);
|
||||
free_dvector(r, 0, N);
|
||||
free_dvector(dr, 0, N);
|
||||
free_ivector(indx, 0, N);
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
void TwoPunctures::Newton(int const nvar, int const n1, int const n2, int const n3,
|
||||
@@ -1371,11 +1291,10 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
|
||||
derivs U;
|
||||
double *sources;
|
||||
|
||||
values = ws_fov_values;
|
||||
U = ws_fov_U;
|
||||
values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&U, nvar);
|
||||
|
||||
sources = ws_fov_sources;
|
||||
memset(sources, 0, n1 * n2 * n3 * sizeof(double));
|
||||
sources = (double *)calloc(n1 * n2 * n3, sizeof(double));
|
||||
if (0)
|
||||
{
|
||||
double *s_x, *s_y, *s_z;
|
||||
@@ -1530,6 +1449,9 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
|
||||
{
|
||||
fclose(debugfile);
|
||||
}
|
||||
free(sources);
|
||||
free_dvector(values, 0, nvar - 1);
|
||||
free_derivs(&U, nvar);
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
double TwoPunctures::norm_inf(double const *F, int const ntotal)
|
||||
@@ -1935,12 +1857,11 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl
|
||||
|
||||
Derivatives_AB3(nvar, n1, n2, n3, dv);
|
||||
|
||||
values = ws_jtdv_values;
|
||||
dU = ws_jtdv_dU;
|
||||
U = ws_jtdv_U;
|
||||
|
||||
for (i = 0; i < n1; i++)
|
||||
{
|
||||
values = dvector(0, nvar - 1);
|
||||
allocate_derivs(&dU, nvar);
|
||||
allocate_derivs(&U, nvar);
|
||||
for (j = 0; j < n2; j++)
|
||||
{
|
||||
for (k = 0; k < n3; k++)
|
||||
@@ -1994,6 +1915,9 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl
|
||||
}
|
||||
}
|
||||
}
|
||||
free_dvector(values, 0, nvar - 1);
|
||||
free_derivs(&dU, nvar);
|
||||
free_derivs(&U, nvar);
|
||||
}
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
@@ -2040,11 +1964,17 @@ void TwoPunctures::LineRelax_be(double *dv,
|
||||
{
|
||||
int j, m, Ic, Ip, Im, col, ivar;
|
||||
|
||||
double *diag = ws_diag_be;
|
||||
double *e = ws_e_be; /* above diagonal */
|
||||
double *f = ws_f_be; /* below diagonal */
|
||||
double *b = ws_b_be; /* rhs */
|
||||
double *x = ws_x_be; /* solution vector */
|
||||
double *diag = new double[n2];
|
||||
double *e = new double[n2 - 1]; /* above diagonal */
|
||||
double *f = new double[n2 - 1]; /* below diagonal */
|
||||
double *b = new double[n2]; /* rhs */
|
||||
double *x = new double[n2]; /* solution vector */
|
||||
|
||||
// gsl_vector *diag = gsl_vector_alloc(n2);
|
||||
// gsl_vector *e = gsl_vector_alloc(n2-1); /* above diagonal */
|
||||
// gsl_vector *f = gsl_vector_alloc(n2-1); /* below diagonal */
|
||||
// gsl_vector *b = gsl_vector_alloc(n2); /* rhs */
|
||||
// gsl_vector *x = gsl_vector_alloc(n2); /* solution vector */
|
||||
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
{
|
||||
@@ -2054,35 +1984,62 @@ void TwoPunctures::LineRelax_be(double *dv,
|
||||
}
|
||||
diag[n2 - 1] = 0;
|
||||
|
||||
// gsl_vector_set_zero(diag);
|
||||
// gsl_vector_set_zero(e);
|
||||
// gsl_vector_set_zero(f);
|
||||
for (j = 0; j < n2; j++)
|
||||
{
|
||||
Ip = Index(ivar, i, j + 1, k, nvar, n1, n2, n3);
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
Im = Index(ivar, i, j - 1, k, nvar, n1, n2, n3);
|
||||
b[j] = rhs[Ic];
|
||||
// gsl_vector_set(b,j,rhs[Ic]);
|
||||
for (m = 0; m < ncols[Ic]; m++)
|
||||
{
|
||||
col = cols[Ic][m];
|
||||
if (col != Ip && col != Ic && col != Im)
|
||||
b[j] -= JFD[Ic][m] * dv[col];
|
||||
// *gsl_vector_ptr(b, j) -= JFD[Ic][m] * dv[col];
|
||||
else
|
||||
{
|
||||
if (col == Im && j > 0)
|
||||
f[j - 1] = JFD[Ic][m];
|
||||
// gsl_vector_set(f,j-1,JFD[Ic][m]);
|
||||
if (col == Ic)
|
||||
diag[j] = JFD[Ic][m];
|
||||
// gsl_vector_set(diag,j,JFD[Ic][m]);
|
||||
if (col == Ip && j < n2 - 1)
|
||||
e[j] = JFD[Ic][m];
|
||||
// gsl_vector_set(e,j,JFD[Ic][m]);
|
||||
}
|
||||
}
|
||||
}
|
||||
// A x = b
|
||||
// A = ( d_0 e_0 0 0 )
|
||||
// ( f_0 d_1 e_1 0 )
|
||||
// ( 0 f_1 d_2 e_2 )
|
||||
// ( 0 0 f_2 d_3 )
|
||||
//
|
||||
ThomasAlgorithm(n2, f, diag, e, x, b);
|
||||
// gsl_linalg_solve_tridiag(diag, e, f, b, x);
|
||||
for (j = 0; j < n2; j++)
|
||||
{
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
dv[Ic] = x[j];
|
||||
// dv[Ic] = gsl_vector_get(x, j);
|
||||
}
|
||||
}
|
||||
|
||||
delete[] diag;
|
||||
delete[] e;
|
||||
delete[] f;
|
||||
delete[] b;
|
||||
delete[] x;
|
||||
// gsl_vector_free(diag);
|
||||
// gsl_vector_free(e);
|
||||
// gsl_vector_free(f);
|
||||
// gsl_vector_free(b);
|
||||
// gsl_vector_free(x);
|
||||
}
|
||||
/* --------------------------------------------------------------------------*/
|
||||
void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
|
||||
@@ -2099,8 +2056,8 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
|
||||
ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
|
||||
derivs dU, U;
|
||||
|
||||
dU = ws_jfd_dU;
|
||||
U = ws_jfd_U;
|
||||
allocate_derivs(&dU, nvar);
|
||||
allocate_derivs(&U, nvar);
|
||||
|
||||
if (k < 0)
|
||||
k = k + n3;
|
||||
@@ -2218,6 +2175,9 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
|
||||
LinEquations(A, B, X, R, x, r, phi, y, z, dU, U, values);
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
values[ivar] *= FAC;
|
||||
|
||||
free_derivs(&dU, nvar);
|
||||
free_derivs(&U, nvar);
|
||||
}
|
||||
#undef FAC
|
||||
/*-----------------------------------------------------------*/
|
||||
@@ -2249,11 +2209,17 @@ void TwoPunctures::LineRelax_al(double *dv,
|
||||
{
|
||||
int i, m, Ic, Ip, Im, col, ivar;
|
||||
|
||||
double *diag = ws_diag_al;
|
||||
double *e = ws_e_al; /* above diagonal */
|
||||
double *f = ws_f_al; /* below diagonal */
|
||||
double *b = ws_b_al; /* rhs */
|
||||
double *x = ws_x_al; /* solution vector */
|
||||
double *diag = new double[n1];
|
||||
double *e = new double[n1 - 1]; /* above diagonal */
|
||||
double *f = new double[n1 - 1]; /* below diagonal */
|
||||
double *b = new double[n1]; /* rhs */
|
||||
double *x = new double[n1]; /* solution vector */
|
||||
|
||||
// gsl_vector *diag = gsl_vector_alloc(n1);
|
||||
// gsl_vector *e = gsl_vector_alloc(n1-1); /* above diagonal */
|
||||
// gsl_vector *f = gsl_vector_alloc(n1-1); /* below diagonal */
|
||||
// gsl_vector *b = gsl_vector_alloc(n1); /* rhs */
|
||||
// gsl_vector *x = gsl_vector_alloc(n1); /* solution vector */
|
||||
|
||||
for (ivar = 0; ivar < nvar; ivar++)
|
||||
{
|
||||
@@ -2263,35 +2229,57 @@ void TwoPunctures::LineRelax_al(double *dv,
|
||||
}
|
||||
diag[n1 - 1] = 0;
|
||||
|
||||
// gsl_vector_set_zero(diag);
|
||||
// gsl_vector_set_zero(e);
|
||||
// gsl_vector_set_zero(f);
|
||||
for (i = 0; i < n1; i++)
|
||||
{
|
||||
Ip = Index(ivar, i + 1, j, k, nvar, n1, n2, n3);
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
Im = Index(ivar, i - 1, j, k, nvar, n1, n2, n3);
|
||||
b[i] = rhs[Ic];
|
||||
// gsl_vector_set(b,i,rhs[Ic]);
|
||||
for (m = 0; m < ncols[Ic]; m++)
|
||||
{
|
||||
col = cols[Ic][m];
|
||||
if (col != Ip && col != Ic && col != Im)
|
||||
b[i] -= JFD[Ic][m] * dv[col];
|
||||
// *gsl_vector_ptr(b, i) -= JFD[Ic][m] * dv[col];
|
||||
else
|
||||
{
|
||||
if (col == Im && i > 0)
|
||||
f[i - 1] = JFD[Ic][m];
|
||||
// gsl_vector_set(f,i-1,JFD[Ic][m]);
|
||||
if (col == Ic)
|
||||
diag[i] = JFD[Ic][m];
|
||||
// gsl_vector_set(diag,i,JFD[Ic][m]);
|
||||
if (col == Ip && i < n1 - 1)
|
||||
e[i] = JFD[Ic][m];
|
||||
// gsl_vector_set(e,i,JFD[Ic][m]);
|
||||
}
|
||||
}
|
||||
}
|
||||
ThomasAlgorithm(n1, f, diag, e, x, b);
|
||||
// gsl_linalg_solve_tridiag(diag, e, f, b, x);
|
||||
for (i = 0; i < n1; i++)
|
||||
{
|
||||
Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
|
||||
dv[Ic] = x[i];
|
||||
// dv[Ic] = gsl_vector_get(x, i);
|
||||
}
|
||||
}
|
||||
|
||||
delete[] diag;
|
||||
delete[] e;
|
||||
delete[] f;
|
||||
delete[] b;
|
||||
delete[] x;
|
||||
|
||||
// gsl_vector_free(diag);
|
||||
// gsl_vector_free(e);
|
||||
// gsl_vector_free(f);
|
||||
// gsl_vector_free(b);
|
||||
// gsl_vector_free(x);
|
||||
}
|
||||
/* -------------------------------------------------------------------------*/
|
||||
// a[N], b[N-1], c[N-1], x[N], q[N]
|
||||
@@ -2303,29 +2291,44 @@ void TwoPunctures::LineRelax_al(double *dv,
|
||||
//"Parallel Scientific Computing in C++ and MPI" P361
|
||||
void TwoPunctures::ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q)
|
||||
{
|
||||
// In-place Thomas algorithm: uses a[] as d workspace, b[] as l workspace.
|
||||
// c[] is already u (above-diagonal). ws_thomas_y is pre-allocated workspace.
|
||||
int i;
|
||||
double *y = ws_thomas_y;
|
||||
double *l, *u, *d, *y;
|
||||
l = new double[N - 1];
|
||||
u = new double[N - 1];
|
||||
d = new double[N];
|
||||
y = new double[N];
|
||||
|
||||
/* LU Decomposition */
|
||||
d[0] = a[0];
|
||||
u[0] = c[0];
|
||||
|
||||
/* LU Decomposition (in-place: a becomes d, b becomes l) */
|
||||
for (i = 0; i < N - 2; i++)
|
||||
{
|
||||
b[i] = b[i] / a[i];
|
||||
a[i + 1] = a[i + 1] - b[i] * c[i];
|
||||
l[i] = b[i] / d[i];
|
||||
d[i + 1] = a[i + 1] - l[i] * u[i];
|
||||
u[i + 1] = c[i + 1];
|
||||
}
|
||||
b[N - 2] = b[N - 2] / a[N - 2];
|
||||
a[N - 1] = a[N - 1] - b[N - 2] * c[N - 2];
|
||||
|
||||
l[N - 2] = b[N - 2] / d[N - 2];
|
||||
d[N - 1] = a[N - 1] - l[N - 2] * u[N - 2];
|
||||
|
||||
/* Forward Substitution [L][y] = [q] */
|
||||
y[0] = q[0];
|
||||
for (i = 1; i < N; i++)
|
||||
y[i] = q[i] - b[i - 1] * y[i - 1];
|
||||
y[i] = q[i] - l[i - 1] * y[i - 1];
|
||||
|
||||
/* Backward Substitution [U][x] = [y] */
|
||||
x[N - 1] = y[N - 1] / a[N - 1];
|
||||
x[N - 1] = y[N - 1] / d[N - 1];
|
||||
|
||||
for (i = N - 2; i >= 0; i--)
|
||||
x[i] = (y[i] - c[i] * x[i + 1]) / a[i];
|
||||
x[i] = (y[i] - u[i] * x[i + 1]) / d[i];
|
||||
|
||||
delete[] l;
|
||||
delete[] u;
|
||||
delete[] d;
|
||||
delete[] y;
|
||||
|
||||
return;
|
||||
}
|
||||
// --------------------------------------------------------------------------*/
|
||||
// Calculates the value of v at an arbitrary position (x,y,z) if the spectral coefficients are know*/*/
|
||||
|
||||
@@ -42,33 +42,6 @@ private:
|
||||
|
||||
int ntotal;
|
||||
|
||||
// Pre-allocated workspace buffers for hot-path allocation elimination
|
||||
// LineRelax_be workspace (sized for n2)
|
||||
double *ws_diag_be, *ws_e_be, *ws_f_be, *ws_b_be, *ws_x_be;
|
||||
// LineRelax_al workspace (sized for n1)
|
||||
double *ws_diag_al, *ws_e_al, *ws_f_al, *ws_b_al, *ws_x_al;
|
||||
// ThomasAlgorithm workspace (sized for max(n1,n2))
|
||||
double *ws_thomas_y;
|
||||
// JFD_times_dv workspace (sized for nvar)
|
||||
double *ws_jfd_values;
|
||||
derivs ws_jfd_dU, ws_jfd_U;
|
||||
// chebft_Zeros workspace (sized for max(n1,n2,n3)+1)
|
||||
double *ws_cheb_c;
|
||||
// fourft workspace (sized for max(n1,n2,n3)/2+1 each)
|
||||
double *ws_four_a, *ws_four_b;
|
||||
// Derivatives_AB3 workspace
|
||||
double *ws_deriv_p, *ws_deriv_dp, *ws_deriv_d2p;
|
||||
double *ws_deriv_q, *ws_deriv_dq;
|
||||
double *ws_deriv_r, *ws_deriv_dr;
|
||||
int *ws_deriv_indx;
|
||||
// F_of_v workspace
|
||||
double *ws_fov_sources;
|
||||
double *ws_fov_values;
|
||||
derivs ws_fov_U;
|
||||
// J_times_dv workspace
|
||||
double *ws_jtdv_values;
|
||||
derivs ws_jtdv_dU, ws_jtdv_U;
|
||||
|
||||
struct parameters
|
||||
{
|
||||
int nvar, n1, n2, n3;
|
||||
|
||||
@@ -106,8 +106,7 @@
|
||||
call getpbh(BHN,Porg,Mass)
|
||||
#endif
|
||||
|
||||
!!! sanity check (disabled in production builds for performance)
|
||||
#ifdef DEBUG
|
||||
!!! sanity check
|
||||
dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) &
|
||||
+sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) &
|
||||
+sum(Gamx)+sum(Gamy)+sum(Gamz) &
|
||||
@@ -137,7 +136,6 @@
|
||||
gont = 1
|
||||
return
|
||||
endif
|
||||
#endif
|
||||
|
||||
PI = dacos(-ONE)
|
||||
|
||||
@@ -168,8 +166,6 @@
|
||||
call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev)
|
||||
call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
|
||||
|
||||
!$OMP PARALLEL
|
||||
!$OMP WORKSHARE
|
||||
gxx_rhs = - TWO * alpn1 * Axx - F2o3 * gxx * div_beta + &
|
||||
TWO *( gxx * betaxx + gxy * betayx + gxz * betazx)
|
||||
|
||||
@@ -203,8 +199,6 @@
|
||||
gupyy = ( gxx * gzz - gxz * gxz ) / gupzz
|
||||
gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
|
||||
gupzz = ( gxx * gyy - gxy * gxy ) / gupzz
|
||||
!$OMP END WORKSHARE
|
||||
!$OMP END PARALLEL
|
||||
|
||||
if(co == 0)then
|
||||
! Gam^i_Res = Gam^i + gup^ij_,j
|
||||
@@ -238,8 +232,6 @@
|
||||
endif
|
||||
|
||||
! second kind of connection
|
||||
!$OMP PARALLEL
|
||||
!$OMP WORKSHARE
|
||||
Gamxxx =HALF*( gupxx*gxxx + gupxy*(TWO*gxyx - gxxy ) + gupxz*(TWO*gxzx - gxxz ))
|
||||
Gamyxx =HALF*( gupxy*gxxx + gupyy*(TWO*gxyx - gxxy ) + gupyz*(TWO*gxzx - gxxz ))
|
||||
Gamzxx =HALF*( gupxz*gxxx + gupyz*(TWO*gxyx - gxxy ) + gupzz*(TWO*gxzx - gxxz ))
|
||||
@@ -288,8 +280,6 @@
|
||||
(gupxy * gupyz + gupyy * gupxz)* Axy + &
|
||||
(gupxy * gupzz + gupyz * gupxz)* Axz + &
|
||||
(gupyy * gupzz + gupyz * gupyz)* Ayz
|
||||
!$OMP END WORKSHARE
|
||||
!$OMP END PARALLEL
|
||||
|
||||
! Right hand side for Gam^i without shift terms...
|
||||
call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
|
||||
@@ -344,8 +334,6 @@
|
||||
call fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev)
|
||||
call fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev)
|
||||
|
||||
!$OMP PARALLEL
|
||||
!$OMP WORKSHARE
|
||||
Gamx_rhs = Gamx_rhs + F2o3 * Gamxa * div_beta - &
|
||||
Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz + &
|
||||
F1o3 * (gupxx * fxx + gupxy * fxy + gupxz * fxz ) + &
|
||||
@@ -385,8 +373,6 @@
|
||||
gyyz = gxz * Gamxyy + gyz * Gamyyy + gzz * Gamzyy
|
||||
gyzz = gxz * Gamxyz + gyz * Gamyyz + gzz * Gamzyz
|
||||
gzzz = gxz * Gamxzz + gyz * Gamyzz + gzz * Gamzzz
|
||||
!$OMP END WORKSHARE
|
||||
!$OMP END PARALLEL
|
||||
|
||||
!compute Ricci tensor for tilted metric
|
||||
call fdderivs(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev)
|
||||
@@ -413,8 +399,6 @@
|
||||
Ryz = gupxx * fxx + gupyy * fyy + gupzz * fzz + &
|
||||
( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO
|
||||
|
||||
!$OMP PARALLEL
|
||||
!$OMP WORKSHARE
|
||||
Rxx = - HALF * Rxx + &
|
||||
gxx * Gamxx+ gxy * Gamyx + gxz * Gamzx + &
|
||||
Gamxa * gxxx + Gamya * gxyx + Gamza * gxzx + &
|
||||
@@ -615,13 +599,9 @@
|
||||
Gamxyz * gxzz + Gamyyz * gyzz + Gamzyz * gzzz + &
|
||||
Gamxzz * gxzy + Gamyzz * gyzy + Gamzzz * gzzy + &
|
||||
Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz )
|
||||
!$OMP END WORKSHARE
|
||||
!$OMP END PARALLEL
|
||||
!covariant second derivative of chi respect to tilted metric
|
||||
call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
|
||||
|
||||
!$OMP PARALLEL
|
||||
!$OMP WORKSHARE
|
||||
fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz
|
||||
fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz
|
||||
fxz = fxz - Gamxxz * chix - Gamyxz * chiy - Gamzxz * chiz
|
||||
@@ -644,15 +624,11 @@
|
||||
Rxy = Rxy + (fxy - chix*chiy/chin1/TWO + gxy * f)/chin1/TWO
|
||||
Rxz = Rxz + (fxz - chix*chiz/chin1/TWO + gxz * f)/chin1/TWO
|
||||
Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO
|
||||
!$OMP END WORKSHARE
|
||||
!$OMP END PARALLEL
|
||||
|
||||
! covariant second derivatives of the lapse respect to physical metric
|
||||
call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
|
||||
SYM,SYM,SYM,symmetry,Lev)
|
||||
|
||||
!$OMP PARALLEL
|
||||
!$OMP WORKSHARE
|
||||
gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1
|
||||
gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1
|
||||
gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1
|
||||
@@ -813,8 +789,6 @@
|
||||
!!!! gauge variable part
|
||||
|
||||
Lap_rhs = -TWO*alpn1*trK
|
||||
!$OMP END WORKSHARE
|
||||
!$OMP END PARALLEL
|
||||
#if (GAUGE == 0)
|
||||
betax_rhs = FF*dtSfx
|
||||
betay_rhs = FF*dtSfy
|
||||
|
||||
@@ -997,10 +997,10 @@
|
||||
fy = ZEO
|
||||
fz = ZEO
|
||||
|
||||
#if 0
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
#if 0
|
||||
! x direction
|
||||
if(i+2 <= imax .and. i-2 >= imin)then
|
||||
!
|
||||
@@ -1040,13 +1040,9 @@
|
||||
|
||||
! set kmax and kmin 0
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#elif 0
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
! x direction
|
||||
if(i+2 <= imax .and. i-2 >= imin)then
|
||||
!
|
||||
! f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
|
||||
! fx(i) = ---------------------------------------------
|
||||
@@ -1083,32 +1079,8 @@
|
||||
|
||||
! set kmax and kmin 0
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#else
|
||||
! for bam comparison — split into branch-free interior + serial boundary
|
||||
! Interior: all stencil points guaranteed in-bounds, no branches needed
|
||||
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
|
||||
do k=max(3,1),min(ex(3)-1,kmax-2)
|
||||
do j=max(3,1),min(ex(2)-1,jmax-2)
|
||||
!DIR$ IVDEP
|
||||
do i=max(3,1),min(ex(1)-1,imax-2)
|
||||
fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
|
||||
fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
|
||||
fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
! Boundary shell: original branching logic for points near edges
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
if(i >= 3 .and. i <= imax-2 .and. &
|
||||
j >= 3 .and. j <= jmax-2 .and. &
|
||||
k >= 3 .and. k <= kmax-2) cycle
|
||||
! for bam comparison
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. &
|
||||
j+2 <= jmax .and. j-2 >= jmin .and. &
|
||||
k+2 <= kmax .and. k-2 >= kmin) then
|
||||
@@ -1122,10 +1094,10 @@
|
||||
fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
|
||||
fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
@@ -1429,10 +1401,10 @@
|
||||
fxz = ZEO
|
||||
fyz = ZEO
|
||||
|
||||
#if 0
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
#if 0
|
||||
!~~~~~~ fxx
|
||||
if(i+2 <= imax .and. i-2 >= imin)then
|
||||
!
|
||||
@@ -1510,47 +1482,8 @@
|
||||
elseif(j+1 <= jmax .and. j-1 >= jmin .and. k+1 <= kmax .and. k-1 >= kmin)then
|
||||
fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#else
|
||||
! for bam comparison — split into branch-free interior + serial boundary
|
||||
! Interior: all stencil points guaranteed in-bounds, no branches needed
|
||||
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
|
||||
do k=max(3,1),min(ex(3)-1,kmax-2)
|
||||
do j=max(3,1),min(ex(2)-1,jmax-2)
|
||||
!DIR$ IVDEP
|
||||
do i=max(3,1),min(ex(1)-1,imax-2)
|
||||
fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) &
|
||||
-fh(i+2,j,k)+F16*fh(i+1,j,k) )
|
||||
fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) &
|
||||
-fh(i,j+2,k)+F16*fh(i,j+1,k) )
|
||||
fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) &
|
||||
-fh(i,j,k+2)+F16*fh(i,j,k+1) )
|
||||
fxy(i,j,k) = Fdxdy*( (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k)) &
|
||||
-F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k)) &
|
||||
+F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k)) &
|
||||
- (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k)))
|
||||
fxz(i,j,k) = Fdxdz*( (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2)) &
|
||||
-F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1)) &
|
||||
+F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1)) &
|
||||
- (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2)))
|
||||
fyz(i,j,k) = Fdydz*( (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2)) &
|
||||
-F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1)) &
|
||||
+F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1)) &
|
||||
- (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
! Boundary shell: original branching logic for points near edges
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
if(i >= 3 .and. i <= imax-2 .and. &
|
||||
j >= 3 .and. j <= jmax-2 .and. &
|
||||
k >= 3 .and. k <= kmax-2) cycle
|
||||
! for bam comparison
|
||||
if(i+2 <= imax .and. i-2 >= imin .and. &
|
||||
j+2 <= jmax .and. j-2 >= jmin .and. &
|
||||
k+2 <= kmax .and. k-2 >= kmin) then
|
||||
@@ -1585,10 +1518,10 @@
|
||||
fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1))
|
||||
fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
return
|
||||
|
||||
|
||||
@@ -19,60 +19,48 @@
|
||||
|
||||
!~~~~~~~> 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, 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
|
||||
|
||||
!~~~~~~>
|
||||
|
||||
do k=1,ex(3)
|
||||
do j=1,ex(2)
|
||||
do i=1,ex(1)
|
||||
gxx = dxx + ONE
|
||||
gyy = dyy + ONE
|
||||
gzz = dzz + ONE
|
||||
|
||||
lgxx = dxx(i,j,k) + ONE
|
||||
lgyy = dyy(i,j,k) + ONE
|
||||
lgzz = dzz(i,j,k) + 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
|
||||
|
||||
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)
|
||||
trA = gupxx * Axx + gupyy * Ayy + gupzz * Azz &
|
||||
+ TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
|
||||
|
||||
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
|
||||
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
|
||||
|
||||
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))
|
||||
detg = ONE / ( detg ** F1o3 )
|
||||
|
||||
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
|
||||
gxx = gxx * detg
|
||||
gxy = gxy * detg
|
||||
gxz = gxz * detg
|
||||
gyy = gyy * detg
|
||||
gyz = gyz * detg
|
||||
gzz = gzz * detg
|
||||
|
||||
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
|
||||
dxx = gxx - ONE
|
||||
dyy = gyy - ONE
|
||||
dzz = gzz - ONE
|
||||
|
||||
return
|
||||
|
||||
@@ -95,70 +83,50 @@
|
||||
|
||||
!~~~~~~~> 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, 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
|
||||
|
||||
!~~~~~~>
|
||||
|
||||
do k=1,ex(3)
|
||||
do j=1,ex(2)
|
||||
do i=1,ex(1)
|
||||
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
|
||||
|
||||
! 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)
|
||||
gupzz = ONE / ( gupzz ** F1o3 )
|
||||
|
||||
lscale = lgxx * lgyy * lgzz + lgxy * lgyz * lgxz &
|
||||
+ lgxz * lgxy * lgyz - lgxz * lgyy * lgxz &
|
||||
- lgxy * lgxy * lgzz - lgxx * lgyz * lgyz
|
||||
gxx = gxx * gupzz
|
||||
gxy = gxy * gupzz
|
||||
gxz = gxz * gupzz
|
||||
gyy = gyy * gupzz
|
||||
gyz = gyz * gupzz
|
||||
gzz = gzz * gupzz
|
||||
|
||||
lscale = ONE / ( lscale ** F1o3 )
|
||||
dxx = gxx - ONE
|
||||
dyy = gyy - ONE
|
||||
dzz = gzz - ONE
|
||||
! for A
|
||||
|
||||
lgxx = lgxx * lscale
|
||||
lgxy = lgxy * lscale
|
||||
lgxz = lgxz * lscale
|
||||
lgyy = lgyy * lscale
|
||||
lgyz = lgyz * lscale
|
||||
lgzz = lgzz * lscale
|
||||
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 )
|
||||
|
||||
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
|
||||
trA = gupxx * Axx + gupyy * Ayy + gupzz * Azz &
|
||||
+ TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
|
||||
|
||||
! 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
|
||||
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
|
||||
|
||||
|
||||
@@ -324,6 +324,7 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
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)
|
||||
@@ -349,6 +350,7 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
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)
|
||||
@@ -377,6 +379,7 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
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)
|
||||
@@ -881,18 +884,10 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
|
||||
real*8, dimension(-ord+1:extc(1),-ord+1:extc(2),-ord+1:extc(3)),intent(out):: funcc
|
||||
real*8, dimension(1:3), intent(in) :: SoA
|
||||
|
||||
integer::i,j,k
|
||||
|
||||
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
|
||||
do k=1,extc(3)
|
||||
do j=1,extc(2)
|
||||
do i=1,extc(1)
|
||||
funcc(i,j,k) = func(i,j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
integer::i
|
||||
|
||||
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
|
||||
@@ -917,6 +912,7 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
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)
|
||||
@@ -945,6 +941,7 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
|
||||
|
||||
integer::i
|
||||
|
||||
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)
|
||||
@@ -1122,16 +1119,25 @@ end subroutine d2dump
|
||||
!------------------------------------------------------------------------------
|
||||
|
||||
subroutine polint(xa,ya,x,y,dy,ordn)
|
||||
|
||||
implicit none
|
||||
|
||||
!~~~~~~> Input Parameter:
|
||||
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
|
||||
!~~~~~~> Other parameter:
|
||||
|
||||
integer :: m,n,ns
|
||||
real*8, dimension(ordn) :: c,d,den,ho
|
||||
real*8 :: dif,dift
|
||||
|
||||
!~~~~~~>
|
||||
|
||||
n=ordn
|
||||
m=ordn
|
||||
|
||||
c=ya
|
||||
d=ya
|
||||
@@ -1139,38 +1145,27 @@ end subroutine d2dump
|
||||
|
||||
ns=1
|
||||
dif=abs(x-xa(1))
|
||||
|
||||
do i = 2, ordn
|
||||
dift = abs(x - xa(i))
|
||||
do m=1,n
|
||||
dift=abs(x-xa(m))
|
||||
if(dift < dif) then
|
||||
ns = i
|
||||
ns=m
|
||||
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
|
||||
do m=1,n-1
|
||||
den(1:n-m)=ho(1:n-m)-ho(1+m:n)
|
||||
if (any(den(1:n-m) == 0.0))then
|
||||
write(*,*) 'failure in polint for point',x
|
||||
write(*,*) 'with input points: ',xa
|
||||
stop
|
||||
endif
|
||||
|
||||
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
|
||||
den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
|
||||
d(1:n-m)=ho(1+m:n)*den(1:n-m)
|
||||
c(1:n-m)=ho(1:n-m)*den(1:n-m)
|
||||
if (2*ns < n-m) then
|
||||
dy=c(ns+1)
|
||||
else
|
||||
dy=d(ns)
|
||||
@@ -1180,6 +1175,7 @@ end subroutine d2dump
|
||||
end do
|
||||
|
||||
return
|
||||
|
||||
end subroutine polint
|
||||
!------------------------------------------------------------------------------
|
||||
!
|
||||
@@ -1187,37 +1183,35 @@ end subroutine d2dump
|
||||
!
|
||||
!------------------------------------------------------------------------------
|
||||
subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
|
||||
|
||||
implicit none
|
||||
|
||||
!~~~~~~> Input parameters:
|
||||
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
|
||||
!~~~~~~> Other parameters:
|
||||
|
||||
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
|
||||
|
||||
call polint(x1a,ymtmp,x1,y,dy,ordn)
|
||||
|
||||
return
|
||||
|
||||
end subroutine polin2
|
||||
!------------------------------------------------------------------------------
|
||||
!
|
||||
@@ -1225,15 +1219,18 @@ end subroutine d2dump
|
||||
!
|
||||
!------------------------------------------------------------------------------
|
||||
subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
|
||||
|
||||
implicit none
|
||||
|
||||
!~~~~~~> Input parameters:
|
||||
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
|
||||
!~~~~~~> Other parameters:
|
||||
|
||||
integer :: i,j,m,n
|
||||
real*8, dimension(ordn,ordn) :: yatmp
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
@@ -1242,33 +1239,24 @@ end subroutine d2dump
|
||||
|
||||
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 :: j, k
|
||||
real*8, dimension(ordn,ordn) :: yatmp
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
real*8 :: dy_temp
|
||||
|
||||
do k=1,ordn
|
||||
do j=1,ordn
|
||||
call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
|
||||
end do
|
||||
end do
|
||||
do k=1,ordn
|
||||
call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
|
||||
end do
|
||||
call polint(x3a, ymtmp, x3, y, dy, ordn)
|
||||
#endif
|
||||
|
||||
call polint(x1a,ymtmp,x1,y,dy,ordn)
|
||||
|
||||
return
|
||||
|
||||
end subroutine polin3
|
||||
!--------------------------------------------------------------------------------------
|
||||
! calculate L2norm
|
||||
@@ -1288,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)
|
||||
@@ -1314,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
|
||||
|
||||
@@ -1344,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
|
||||
|
||||
@@ -1409,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
|
||||
|
||||
@@ -1442,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
|
||||
|
||||
@@ -1506,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
|
||||
|
||||
@@ -1709,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
|
||||
@@ -1745,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
|
||||
|
||||
@@ -1789,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
|
||||
@@ -1819,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
|
||||
|
||||
@@ -1857,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
|
||||
@@ -1918,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
|
||||
|
||||
@@ -2151,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
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
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;
|
||||
}
|
||||
|
||||
// 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;
|
||||
swap = b[irow];
|
||||
b[irow] = b[icol];
|
||||
b[icol] = swap;
|
||||
}
|
||||
|
||||
// Then compute inverse from LU factorization
|
||||
info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, a, n, ipiv);
|
||||
indxr[i] = irow;
|
||||
indxc[i] = icol;
|
||||
|
||||
if (info != 0) {
|
||||
cout << "gaussj: Singular Matrix (dgetri info=" << info << ")" << endl;
|
||||
delete[] ipiv;
|
||||
delete[] a_copy;
|
||||
return 1;
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
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
|
||||
|
||||
@@ -159,42 +159,36 @@ integer, parameter :: NO_SYMM=0, OCTANT=2
|
||||
|
||||
call symmetry_bd(3,ex,f,fh,SoA)
|
||||
|
||||
! Interior: all stencil points guaranteed in-bounds
|
||||
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
|
||||
do k=4,ex(3)-3
|
||||
do j=4,ex(2)-3
|
||||
!DIR$ IVDEP
|
||||
do i=4,ex(1)-3
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/cof *( ( &
|
||||
(fh(i-3,j,k)+fh(i+3,j,k)) - &
|
||||
SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
|
||||
FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
|
||||
TWT* fh(i,j,k) )/dX + &
|
||||
( &
|
||||
(fh(i,j-3,k)+fh(i,j+3,k)) - &
|
||||
SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
|
||||
FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
|
||||
TWT* fh(i,j,k) )/dY + &
|
||||
( &
|
||||
(fh(i,j,k-3)+fh(i,j,k+3)) - &
|
||||
SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
|
||||
FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
|
||||
TWT* fh(i,j,k) )/dZ )
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
! Boundary shell: original branching logic for points near edges
|
||||
do k=1,ex(3)
|
||||
do j=1,ex(2)
|
||||
do i=1,ex(1)
|
||||
if(i >= 4 .and. i <= ex(1)-3 .and. &
|
||||
j >= 4 .and. j <= ex(2)-3 .and. &
|
||||
k >= 4 .and. k <= ex(3)-3) cycle
|
||||
|
||||
if(i-3 >= imin .and. i+3 <= imax .and. &
|
||||
j-3 >= jmin .and. j+3 <= jmax .and. &
|
||||
k-3 >= kmin .and. k+3 <= kmax) then
|
||||
#if 0
|
||||
! x direction
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dX/cof * ( &
|
||||
(fh(i-3,j,k)+fh(i+3,j,k)) - &
|
||||
SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
|
||||
FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
|
||||
TWT* fh(i,j,k) )
|
||||
! y direction
|
||||
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dY/cof * ( &
|
||||
(fh(i,j-3,k)+fh(i,j+3,k)) - &
|
||||
SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
|
||||
FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
|
||||
TWT* fh(i,j,k) )
|
||||
! z direction
|
||||
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/dZ/cof * ( &
|
||||
(fh(i,j,k-3)+fh(i,j,k+3)) - &
|
||||
SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
|
||||
FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
|
||||
TWT* fh(i,j,k) )
|
||||
#else
|
||||
! calculation order if important ?
|
||||
f_rhs(i,j,k) = f_rhs(i,j,k) + eps/cof *( ( &
|
||||
(fh(i-3,j,k)+fh(i+3,j,k)) - &
|
||||
SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
|
||||
@@ -210,7 +204,9 @@ integer, parameter :: NO_SYMM=0, OCTANT=2
|
||||
SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
|
||||
FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
|
||||
TWT* fh(i,j,k) )/dZ )
|
||||
#endif
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@@ -233,7 +233,6 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)
|
||||
|
||||
! upper bound set ex-1 only for efficiency,
|
||||
! the loop body will set ex 0 also
|
||||
!$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
|
||||
do k=1,ex(3)-1
|
||||
do j=1,ex(2)-1
|
||||
do i=1,ex(1)-1
|
||||
@@ -483,7 +482,6 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
return
|
||||
|
||||
|
||||
@@ -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_intel_thread -lmkl_core -lifcore -limf -lpthread -lm -ldl -qopenmp
|
||||
## 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 -ipo -qopenmp \
|
||||
-Dfortran3 -Dnewc -I${MKLROOT}/include
|
||||
f90appflags = -O3 -xHost -fp-model fast=2 -fma -ipo -qopenmp \
|
||||
-align array64byte -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