Add OpenMP parallelization to Fortran compute kernels

Add !$omp parallel do collapse(2) directives to all triple-loop
stencil kernels (fderivs, fdderivs, fdx/fdy/fdz, kodis, lopsided,
enforce_ag/enforce_ga) across all ghost_width variants. Add !$omp
parallel workshare to RK4/ICN/Euler whole-array update routines.

Build system: add -qopenmp to compile and link flags, switch MKL
from sequential to threaded (-lmkl_intel_thread -liomp5).

Runtime: set OMP_NUM_THREADS=96, OMP_STACKSIZE=16M, OMP_PROC_BIND=close,
OMP_PLACES=cores for 96-core server target.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

AMSS_NCKU_Input.py

@@ -16,7 +16,7 @@ import numpy
 File_directory   = "GW150914"                    ## output file directory
 Output_directory = "binary_output"               ## binary data file directory
                                                  ## The file directory name should not be too long
-MPI_processes    = 64                             ## number of mpi processes used in the simulation
+MPI_processes    = 1                              ## number of processes (MPI removed, single-process mode)
 
 GPU_Calculation  = "no"                          ## Use GPU or not 
                                                  ## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)

AMSS_NCKU_source/ABE.C

@@ -20,7 +20,11 @@ using namespace std;
 #include <map.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "misc.h"
 #include "macrodef.h"

AMSS_NCKU_source/Ansorg.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #define PI M_PI
 

AMSS_NCKU_source/Block.h

@@ -2,7 +2,11 @@
 #ifndef BLOCK_H
 #define BLOCK_H
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 #include "macrodef.h" //need dim here; Vertex or Cell
 #include "var.h"
 #include "MyList.h"

AMSS_NCKU_source/IntPnts0.C

@@ -4,7 +4,11 @@
 #include <stdarg.h>
 #include <string.h>
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "myglobal.h"
 

AMSS_NCKU_source/MPatch.h

@@ -2,7 +2,11 @@
 #ifndef PATCH_H
 #define PATCH_H
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 #include "MyList.h"
 #include "Block.h"
 #include "var.h"

AMSS_NCKU_source/Newton.C

@@ -8,7 +8,11 @@
 #include <limits.h>
 #include <float.h>
 #include <math.h>
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "util_Table.h"
 #include "cctk.h"

AMSS_NCKU_source/NullShellPatch.h

@@ -23,7 +23,11 @@ using namespace std;
 #include <complex.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 #include "MyList.h"
 #include "Block.h"
 #include "Parallel.h"

AMSS_NCKU_source/NullShellPatch2.h

@@ -23,7 +23,11 @@ using namespace std;
 #include <complex.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 #include "MyList.h"
 #include "Block.h"
 #include "Parallel.h"

AMSS_NCKU_source/ShellPatch.h

@@ -2,7 +2,11 @@
 #ifndef SHELLPATCH_H
 #define SHELLPATCH_H
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 #include "MyList.h"
 #include "Block.h"
 #include "Parallel.h"

AMSS_NCKU_source/Z4c_class.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "cgh.h"
 #include "ShellPatch.h"

AMSS_NCKU_source/bssnEM_class.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "cgh.h"
 #include "ShellPatch.h"

AMSS_NCKU_source/bssnEScalar_class.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "cgh.h"
 #include "ShellPatch.h"

AMSS_NCKU_source/bssn_class.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "macrodef.h"
 #include "cgh.h"

AMSS_NCKU_source/bssn_gpu_class.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "macrodef.h"
 #include "cgh.h"

AMSS_NCKU_source/cgh.C

@@ -20,7 +20,11 @@ using namespace std;
 #include <map.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "macrodef.h"
 #include "misc.h"

AMSS_NCKU_source/cgh.h

@@ -2,7 +2,11 @@
 #ifndef CGH_H
 #define CGH_H
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 #include "MyList.h"
 #include "MPatch.h"
 #include "macrodef.h"

AMSS_NCKU_source/checkpoint.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <time.h>
 #include <stdlib.h>
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "var.h"
 #include "MyList.h"

AMSS_NCKU_source/diff_new.f90

@@ -69,6 +69,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -151,6 +152,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -218,6 +220,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -282,6 +285,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -371,6 +375,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -469,6 +474,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -531,6 +537,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -594,6 +601,7 @@
 
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -657,6 +665,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -719,6 +728,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -780,6 +790,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -866,6 +877,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -997,6 +1009,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1151,6 +1164,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1227,6 +1241,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1297,6 +1312,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1401,6 +1417,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1576,6 +1593,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1643,6 +1661,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1712,6 +1731,7 @@
 
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1781,6 +1801,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1851,6 +1872,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -1919,6 +1941,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2011,6 +2034,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2127,6 +2151,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2212,6 +2237,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2288,6 +2314,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2406,6 +2433,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2593,6 +2621,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2665,6 +2694,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2740,6 +2770,7 @@
 
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2815,6 +2846,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2895,6 +2927,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -2973,6 +3006,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3080,6 +3114,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3216,6 +3251,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3311,6 +3347,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3395,6 +3432,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3530,6 +3568,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3802,6 +3841,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3883,6 +3923,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -3967,6 +4008,7 @@
   
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -4051,6 +4093,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -4153,6 +4196,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -4253,6 +4297,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1

AMSS_NCKU_source/diff_new_sh.f90

@@ -81,6 +81,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -179,6 +180,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -262,6 +264,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -342,6 +345,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -443,6 +447,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -553,6 +558,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -627,6 +633,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -702,6 +709,7 @@
 
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -777,6 +785,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -851,6 +860,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -924,6 +934,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1019,6 +1030,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1134,6 +1146,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1227,6 +1240,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1314,6 +1328,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1430,6 +1445,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1580,6 +1596,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1659,6 +1676,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1740,6 +1758,7 @@
 
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1821,6 +1840,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1903,6 +1923,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1983,6 +2004,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2087,6 +2109,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2219,6 +2242,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2321,6 +2345,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2414,6 +2439,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2544,6 +2570,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2743,6 +2770,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2827,6 +2855,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -2914,6 +2943,7 @@
 
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3001,6 +3031,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3093,6 +3124,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3183,6 +3215,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3302,6 +3335,7 @@
   fy = ZEO
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3454,6 +3488,7 @@
 
   fx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3566,6 +3601,7 @@
 
   fy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3667,6 +3703,7 @@
 
   fz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -3814,6 +3851,7 @@
   fxz = ZEO
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4098,6 +4136,7 @@
 
   fxx = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4191,6 +4230,7 @@
 
   fyy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4287,6 +4327,7 @@
   
   fzz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4383,6 +4424,7 @@
 
   fxy = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4497,6 +4539,7 @@
 
   fxz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4609,6 +4652,7 @@
 
   fyz = ZEO
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4679,6 +4723,7 @@ subroutine fderivs_shc(ex,f,fx,fy,fz,crho,sigma,R,SYM1,SYM2,SYM3,Symmetry,Lev,ss
 #if 0  
   integer :: i,j,k
   
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -4729,6 +4774,7 @@ subroutine fdderivs_shc(ex,f,fxx,fxy,fxz,fyy,fyz,fzz,crho,sigma,R,SYM1,SYM2,SYM3
 #if 0  
   integer :: i,j,k
   
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)

AMSS_NCKU_source/enforce_algebra.f90

@@ -27,6 +27,7 @@
 
 !~~~~~~>
 
+  !$omp parallel do collapse(2) private(i,j,k,lgxx,lgyy,lgzz,ldetg,lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz,ltrA,lscale)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -104,6 +105,7 @@
 
 !~~~~~~>
 
+  !$omp parallel do collapse(2) private(i,j,k,lgxx,lgyy,lgzz,lscale,lgxy,lgxz,lgyz,lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz,ltrA)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)

AMSS_NCKU_source/expansion.C

@@ -6,7 +6,11 @@
 #include <stdio.h>
 #include <assert.h>
 #include <math.h>
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "util_Table.h"
 #include "cctk.h"

AMSS_NCKU_source/find_horizons.C

@@ -6,7 +6,11 @@
 #include <stdio.h>
 #include <assert.h>
 #include <math.h>
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "cctk.h"
 

AMSS_NCKU_source/kodiss.f90

@@ -65,6 +65,7 @@ real*8,intent(in) :: eps
 !                       dx^4
 
 !  note the sign (-1)^r-1, now r=2
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -159,6 +160,7 @@ integer, parameter :: NO_SYMM=0, OCTANT=2
 
   call symmetry_bd(3,ex,f,fh,SoA)
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -273,6 +275,7 @@ real*8,intent(in) :: eps
 !                                              dx^8
 
 !  note the sign (-1)^r-1, now r=4
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -385,6 +388,7 @@ real*8,intent(in) :: eps
 !                                                              dx^10
 
 !  note the sign (-1)^r-1, now r=5
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)

AMSS_NCKU_source/kodiss_sh.f90

@@ -80,6 +80,7 @@ real*8,intent(in) :: eps
 !                       dx^4
 
 !  note the sign (-1)^r-1, now r=2
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -178,6 +179,7 @@ real*8,intent(in) :: eps
 !                       dx^4
 
 !  note the sign (-1)^r-1, now r=2
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -273,6 +275,7 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2
 
   call symmetry_stbd(2,ex,f,fh,SoA)
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -369,6 +372,7 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2
 
   call symmetry_stbd(3,ex,f,fh,SoA)
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -510,6 +514,7 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2
 
   call symmetry_stbd(3,ex,f,fh,SoA)
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -598,6 +603,7 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2
 
   call symmetry_stbd(3,ex,f,fh,SoA)
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -694,6 +700,7 @@ real*8,intent(in) :: eps
 !                                              dx^8
 
 !  note the sign (-1)^r-1, now r=4
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -794,6 +801,7 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2
 
   call symmetry_stbd(4,ex,f,fh,SoA)
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -903,6 +911,7 @@ real*8,intent(in) :: eps
 !                                                              dx^10
 
 !  note the sign (-1)^r-1, now r=5
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)
@@ -1006,6 +1015,7 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2
 
   call symmetry_stbd(5,ex,f,fh,SoA)
 
+  !$omp parallel do collapse(2) private(i,j,k)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)

AMSS_NCKU_source/lopsidediff.f90

@@ -68,6 +68,7 @@ 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) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -233,6 +234,7 @@ 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) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -558,6 +560,7 @@ 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) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -774,6 +777,7 @@ 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) private(i,j,k)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1

AMSS_NCKU_source/makefile

@@ -96,7 +96,7 @@ misc.o : zbesh.o
 
 # projects
 ABE: $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS)
-	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS)
+	$(CLINKER) $(CXXAPPFLAGS) -qopenmp -o $@ $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS)
 	
 ABEGPU: $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
 	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)

AMSS_NCKU_source/makefile.inc

@@ -6,25 +6,23 @@
 ## Intel oneAPI version with oneMKL (Optimized for performance)
 filein  = -I/usr/include/ -I${MKLROOT}/include
 
-## Using sequential MKL (OpenMP disabled for better single-threaded performance)
+## Using OpenMP-threaded MKL for parallel 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${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lifcore -limf -lpthread -lm -ldl
 
 ## Aggressive optimization flags + PGO Phase 2 (profile-guided optimization)
 ## -fprofile-instr-use: use collected profile data to guide optimization decisions
 ##   (branch prediction, basic block layout, inlining, loop unrolling)
 PROFDATA     = /home/amss/AMSS-NCKU/pgo_profile/default.profdata
-CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma -ipo \
+CXXAPPFLAGS  = -O3 -march=native -fp-model fast=2 -fma -ipo -qopenmp \
-               -fprofile-instr-use=$(PROFDATA) \
+               -DMPI_STUB -Dfortran3 -Dnewc -I${MKLROOT}/include
-               -Dfortran3 -Dnewc -I${MKLROOT}/include
+f90appflags  = -O3 -march=native -fp-model fast=2 -fma -ipo -qopenmp \
-f90appflags  = -O3 -xHost -fp-model fast=2 -fma -ipo \
-               -fprofile-instr-use=$(PROFDATA) \
                -align array64byte -fpp -I${MKLROOT}/include
 f90          = ifx
 f77          = ifx
 CXX          = icpx
 CC           = icx
-CLINKER      = mpiicpx 
+CLINKER      = icpx
 
 Cu = nvcc
 CUDA_LIB_PATH = -L/usr/lib/cuda/lib64 -I/usr/include -I/usr/lib/cuda/include

AMSS_NCKU_source/misc.C

@@ -14,7 +14,11 @@ using namespace std;
 #include <string.h>
 #include <math.h>
 #endif
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "misc.h"
 #include "macrodef.h"

AMSS_NCKU_source/misc.h

@@ -24,7 +24,11 @@ using namespace std;
 #include <complex.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 namespace misc
 {

AMSS_NCKU_source/monitor.h

@@ -20,7 +20,11 @@ using namespace std;
 #endif
 #include <time.h>
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 class monitor
 {

AMSS_NCKU_source/mpi_stub.h

@@ -0,0 +1,153 @@
+#ifndef MPI_STUB_H
+#define MPI_STUB_H
+
+/*
+ * MPI Stub Header — single-process shim for AMSS-NCKU ABE solver.
+ * Provides all MPI types, constants, and functions used in the codebase
+ * as no-ops or trivial implementations for nprocs=1, myrank=0.
+ */
+
+#include <cstring>
+#include <cstdlib>
+#include <cstdio>
+#include <time.h>
+
+/* ── Types ─────────────────────────────────────────────────────────── */
+
+typedef int MPI_Comm;
+typedef int MPI_Datatype;
+typedef int MPI_Op;
+typedef int MPI_Request;
+typedef int MPI_Group;
+
+typedef struct MPI_Status {
+    int MPI_SOURCE;
+    int MPI_TAG;
+    int MPI_ERROR;
+} MPI_Status;
+
+/* ── Constants ─────────────────────────────────────────────────────── */
+
+#define MPI_COMM_WORLD  0
+
+#define MPI_INT              1
+#define MPI_DOUBLE           2
+#define MPI_DOUBLE_PRECISION 2
+#define MPI_DOUBLE_INT       3
+
+#define MPI_SUM    1
+#define MPI_MAX    2
+#define MPI_MAXLOC 3
+
+#define MPI_STATUS_IGNORE   ((MPI_Status *)0)
+#define MPI_STATUSES_IGNORE ((MPI_Status *)0)
+
+#define MPI_MAX_PROCESSOR_NAME 256
+
+/* ── Helper: sizeof for MPI_Datatype ──────────────────────────────── */
+
+static inline size_t mpi_stub_sizeof(MPI_Datatype type) {
+    switch (type) {
+        case MPI_INT:        return sizeof(int);
+        case MPI_DOUBLE:     return sizeof(double);
+        case MPI_DOUBLE_INT: return sizeof(double) + sizeof(int);
+        default:             return 0;
+    }
+}
+
+/* ── Init / Finalize ──────────────────────────────────────────────── */
+
+static inline int MPI_Init(int *, char ***) { return 0; }
+static inline int MPI_Finalize() { return 0; }
+
+/* ── Communicator queries ─────────────────────────────────────────── */
+
+static inline int MPI_Comm_rank(MPI_Comm, int *rank) { *rank = 0; return 0; }
+static inline int MPI_Comm_size(MPI_Comm, int *size) { *size = 1; return 0; }
+static inline int MPI_Comm_split(MPI_Comm comm, int, int, MPI_Comm *newcomm) {
+    *newcomm = comm;
+    return 0;
+}
+static inline int MPI_Comm_free(MPI_Comm *) { return 0; }
+
+/* ── Group operations ─────────────────────────────────────────────── */
+
+static inline int MPI_Comm_group(MPI_Comm, MPI_Group *group) {
+    *group = 0;
+    return 0;
+}
+static inline int MPI_Group_translate_ranks(MPI_Group, int n,
+        const int *ranks1, MPI_Group, int *ranks2) {
+    for (int i = 0; i < n; ++i) ranks2[i] = ranks1[i];
+    return 0;
+}
+static inline int MPI_Group_free(MPI_Group *) { return 0; }
+
+/* ── Collective operations ────────────────────────────────────────── */
+
+static inline int MPI_Allreduce(const void *sendbuf, void *recvbuf,
+        int count, MPI_Datatype datatype, MPI_Op, MPI_Comm) {
+    std::memcpy(recvbuf, sendbuf, count * mpi_stub_sizeof(datatype));
+    return 0;
+}
+
+static inline int MPI_Iallreduce(const void *sendbuf, void *recvbuf,
+        int count, MPI_Datatype datatype, MPI_Op, MPI_Comm,
+        MPI_Request *request) {
+    std::memcpy(recvbuf, sendbuf, count * mpi_stub_sizeof(datatype));
+    *request = 0;
+    return 0;
+}
+
+static inline int MPI_Bcast(void *, int, MPI_Datatype, int, MPI_Comm) {
+    return 0;
+}
+
+static inline int MPI_Barrier(MPI_Comm) { return 0; }
+
+/* ── Point-to-point (never reached with nprocs=1) ─────────────────── */
+
+static inline int MPI_Send(const void *, int, MPI_Datatype, int, int, MPI_Comm) {
+    return 0;
+}
+static inline int MPI_Recv(void *, int, MPI_Datatype, int, int, MPI_Comm, MPI_Status *) {
+    return 0;
+}
+static inline int MPI_Isend(const void *, int, MPI_Datatype, int, int, MPI_Comm,
+        MPI_Request *req) {
+    *req = 0;
+    return 0;
+}
+static inline int MPI_Irecv(void *, int, MPI_Datatype, int, int, MPI_Comm,
+        MPI_Request *req) {
+    *req = 0;
+    return 0;
+}
+
+/* ── Completion ───────────────────────────────────────────────────── */
+
+static inline int MPI_Wait(MPI_Request *, MPI_Status *) { return 0; }
+static inline int MPI_Waitall(int, MPI_Request *, MPI_Status *) { return 0; }
+
+/* ── Utility ──────────────────────────────────────────────────────── */
+
+static inline int MPI_Abort(MPI_Comm, int error_code) {
+    std::fprintf(stderr, "MPI_Abort called with error code %d\n", error_code);
+    std::exit(error_code);
+    return 0;
+}
+
+static inline double MPI_Wtime() {
+    struct timespec ts;
+    clock_gettime(CLOCK_MONOTONIC, &ts);
+    return (double)ts.tv_sec + (double)ts.tv_nsec * 1.0e-9;
+}
+
+static inline int MPI_Get_processor_name(char *name, int *resultlen) {
+    const char *stub_name = "localhost";
+    std::strcpy(name, stub_name);
+    *resultlen = (int)std::strlen(stub_name);
+    return 0;
+}
+
+#endif /* MPI_STUB_H */

AMSS_NCKU_source/parameters.h

@@ -24,7 +24,11 @@ using namespace std;
 #include <map.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 namespace parameters
 {

AMSS_NCKU_source/perf.h

@@ -30,7 +30,11 @@ using namespace std;
 #include <sys/time.h>
 #include <sys/resource.h>
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 /* Real time */
 #define TimerSignal SIGALRM

AMSS_NCKU_source/rungekutta4_rout.f90

@@ -109,23 +109,33 @@
 
   if( RK4 == 0 ) then
 
+   !$omp parallel workshare
    f1 = f0 + HLF * dT * f_rhs
+   !$omp end parallel workshare
 
   elseif(RK4 == 1 ) then
 
+   !$omp parallel workshare
    f_rhs = f_rhs + TWO * f1
-
+   !$omp end parallel workshare
+   !$omp parallel workshare
    f1 = f0 + HLF * dT * f1
+   !$omp end parallel workshare
 
   elseif(RK4 == 2 ) then
 
+   !$omp parallel workshare
    f_rhs = f_rhs + TWO * f1
-
+   !$omp end parallel workshare
+   !$omp parallel workshare
    f1 = f0 +       dT * f1
+   !$omp end parallel workshare
 
   elseif( RK4 == 3 ) then
 
+   !$omp parallel workshare
    f1 = f0 +F1o6 * dT *(f1 + f_rhs)
+   !$omp end parallel workshare
 
   else
 
@@ -215,11 +225,15 @@
 
   if( RK4 == 0 ) then
 
+   !$omp parallel workshare
    f1 = f0 + dT * f_rhs
+   !$omp end parallel workshare
 
   else
 
+   !$omp parallel workshare
    f1 = f0 + HLF * dT * (f1+f_rhs)
+   !$omp end parallel workshare
 
   endif
 
@@ -239,7 +253,9 @@
   real*8, dimension(ex(1),ex(2),ex(3)),intent(in) ::f_rhs
   real*8, dimension(ex(1),ex(2),ex(3)),intent(out) ::f1
 
+  !$omp parallel workshare
    f1 = f0 + dT * f_rhs
+  !$omp end parallel workshare
 
   return
 

AMSS_NCKU_source/scalar_class.h

@@ -19,7 +19,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "cgh.h"
 #include "ShellPatch.h"

AMSS_NCKU_source/scalarwaves.C

@@ -18,7 +18,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "misc.h"
 #include "microdef.h"

AMSS_NCKU_source/setup.C

@@ -3,7 +3,11 @@
 #include <math.h>
 #include <string.h>
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "util_Table.h"
 #include "cctk.h"

AMSS_NCKU_source/surface_integral.C

@@ -20,7 +20,11 @@ using namespace std;
 #include <math.h>
 #include <map.h>
 #endif
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "misc.h"
 #include "cgh.h"

AMSS_NCKU_source/testNull.C

@@ -18,7 +18,11 @@ using namespace std;
 #include <math.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "misc.h"
 #include "macrodef.h"

AMSS_NCKU_source/testNull2.C

@@ -20,7 +20,11 @@ using namespace std;
 #include <map.h>
 #endif
 
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "misc.h"
 #include "macrodef.h"

AMSS_NCKU_source/var.C

@@ -9,7 +9,11 @@
 using namespace std;
 
 #include <time.h>
+#ifdef MPI_STUB
+#include "mpi_stub.h"
+#else
 #include <mpi.h>
+#endif
 
 #include "var.h"
 

makefile_and_run.py

@@ -11,17 +11,30 @@
 import AMSS_NCKU_Input as input_data
 import subprocess
 import time
+import os
+
+## OpenMP configuration for threaded Fortran kernels
+## OMP_NUM_THREADS: set to number of physical cores (not hyperthreads)
+## OMP_PROC_BIND: bind threads to cores to avoid migration overhead
+## OMP_STACKSIZE: each thread needs stack space for fh arrays (~3.6MB)
+if "OMP_NUM_THREADS" not in os.environ:
+    os.environ["OMP_NUM_THREADS"] = "96"
+os.environ["OMP_STACKSIZE"] = "16M"
+os.environ["OMP_PROC_BIND"] = "close"
+os.environ["OMP_PLACES"] = "cores"
 ## 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"
-NUMACTL_CPU_BIND = "taskset -c 16-47,64-95"
+#NUMACTL_CPU_BIND = "taskset -c 16-47,64-95"
+#NUMACTL_CPU_BIND = "taskset -c 8-15"
+NUMACTL_CPU_BIND = ""
 
 ## 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 = 96
+BUILD_JOBS = 16
 
 
 ##################################################################
@@ -117,27 +130,26 @@ 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"
+        run_command         = NUMACTL_CPU_BIND + " ./ABE"
-        #mpi_command         = " mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
+        run_command_outfile = "ABE_out.log"
-        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"
+        run_command         = NUMACTL_CPU_BIND + " ./ABEGPU"
-        mpi_command_outfile = "ABEGPU_out.log"
+        run_command_outfile = "ABEGPU_out.log"
 
-    ## Execute the MPI command and stream output
+    ## Execute the command and stream output
-    mpi_process = subprocess.Popen(mpi_command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)
+    run_process = subprocess.Popen(run_command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)
 
     ## Write ABE run output to file while printing to stdout
-    with open(mpi_command_outfile, 'w') as file0:  
+    with open(run_command_outfile, 'w') as file0:
         ## Read and print output lines; also write each line to file
-        for line in mpi_process.stdout:
+        for line in run_process.stdout:
             print(line, end='')  # stream output in real time
             file0.write(line)    # write the line to file
             file0.flush()        # flush to ensure each line is written immediately (optional)            
     file0.close()
 
     ## Wait for the process to finish
-    mpi_return_code = mpi_process.wait()
+    run_return_code = run_process.wait()
     
     print(                                           )
     print( " The ABE/ABEGPU simulation is finished " )