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>
Remove MPI dependency, replace with single-process stub for non-MPI builds
2026-02-10 23:40:17 +08:00 · 2026-02-10 22:51:11 +08:00
43 changed files with 463 additions and 49 deletions
--- a/AMSS_NCKU_Input.py
+++ b/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)
--- a/AMSS_NCKU_source/ABE.C
+++ b/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"
--- a/AMSS_NCKU_source/Ansorg.h
+++ b/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

--- a/AMSS_NCKU_source/Block.h
+++ b/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"
--- a/AMSS_NCKU_source/IntPnts0.C
+++ b/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"

--- a/AMSS_NCKU_source/MPatch.h
+++ b/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"
--- a/AMSS_NCKU_source/Newton.C
+++ b/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"
--- a/AMSS_NCKU_source/NullShellPatch.h
+++ b/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"
--- a/AMSS_NCKU_source/NullShellPatch2.h
+++ b/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"
--- a/AMSS_NCKU_source/ShellPatch.h
+++ b/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"
--- a/AMSS_NCKU_source/Z4c_class.h
+++ b/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"
--- a/AMSS_NCKU_source/bssnEM_class.h
+++ b/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"
--- a/AMSS_NCKU_source/bssnEScalar_class.h
+++ b/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"
--- a/AMSS_NCKU_source/bssn_class.h
+++ b/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"
--- a/AMSS_NCKU_source/bssn_gpu_class.h
+++ b/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"
--- a/AMSS_NCKU_source/cgh.C
+++ b/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"
--- a/AMSS_NCKU_source/cgh.h
+++ b/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"
--- a/AMSS_NCKU_source/checkpoint.h
+++ b/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"
--- a/AMSS_NCKU_source/diff_new.f90
+++ b/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
--- a/AMSS_NCKU_source/diff_new_sh.f90
+++ b/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)
--- a/AMSS_NCKU_source/enforce_algebra.f90
+++ b/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)
--- a/AMSS_NCKU_source/expansion.C
+++ b/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"
--- a/AMSS_NCKU_source/find_horizons.C
+++ b/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"

--- a/AMSS_NCKU_source/kodiss.f90
+++ b/AMSS_NCKU_source/kodiss.f90
@@ -65,7 +65,8 @@ real*8,intent(in) :: eps
 !                       dx^4

 !  note the sign (-1)^r-1, now r=2
-  do k=1,ex(3)
+  !$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,7 +160,8 @@ integer, parameter :: NO_SYMM=0, OCTANT=2

  call symmetry_bd(3,ex,f,fh,SoA)

-  do k=1,ex(3)
+  !$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,7 +275,8 @@ real*8,intent(in) :: eps
 !                                              dx^8

 !  note the sign (-1)^r-1, now r=4
-  do k=1,ex(3)
+  !$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,7 +388,8 @@ real*8,intent(in) :: eps
 !                                                              dx^10

 !  note the sign (-1)^r-1, now r=5
-  do k=1,ex(3)
+  !$omp parallel do collapse(2) private(i,j,k)
+  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)

--- a/AMSS_NCKU_source/kodiss_sh.f90
+++ b/AMSS_NCKU_source/kodiss_sh.f90
@@ -80,7 +80,8 @@ real*8,intent(in) :: eps
 !                       dx^4

 !  note the sign (-1)^r-1, now r=2
-  do k=1,ex(3)
+  !$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,7 +179,8 @@ real*8,intent(in) :: eps
 !                       dx^4

 !  note the sign (-1)^r-1, now r=2
-  do k=1,ex(3)
+  !$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,7 +275,8 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2

  call symmetry_stbd(2,ex,f,fh,SoA)

-  do k=1,ex(3)
+  !$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,7 +372,8 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2

  call symmetry_stbd(3,ex,f,fh,SoA)

-  do k=1,ex(3)
+  !$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,7 +514,8 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2

  call symmetry_stbd(3,ex,f,fh,SoA)

-  do k=1,ex(3)
+  !$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,7 +603,8 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2

  call symmetry_stbd(3,ex,f,fh,SoA)

-  do k=1,ex(3)
+  !$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,7 +700,8 @@ real*8,intent(in) :: eps
 !                                              dx^8

 !  note the sign (-1)^r-1, now r=4
-  do k=1,ex(3)
+  !$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,7 +801,8 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2

  call symmetry_stbd(4,ex,f,fh,SoA)

-  do k=1,ex(3)
+  !$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,7 +911,8 @@ real*8,intent(in) :: eps
 !                                                              dx^10

 !  note the sign (-1)^r-1, now r=5
-  do k=1,ex(3)
+  !$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,7 +1015,8 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2

  call symmetry_stbd(5,ex,f,fh,SoA)

-  do k=1,ex(3)
+  !$omp parallel do collapse(2) private(i,j,k)
+  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)

--- a/AMSS_NCKU_source/lopsidediff.f90
+++ b/AMSS_NCKU_source/lopsidediff.f90
@@ -68,7 +68,8 @@ 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
-  do k=1,ex(3)-1
+  !$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
 ! x direction   
@@ -233,7 +234,8 @@ 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
-  do k=1,ex(3)-1
+  !$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
 #if 0  
@@ -558,7 +560,8 @@ 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
-  do k=1,ex(3)-1
+  !$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
 ! x direction   
@@ -774,7 +777,8 @@ 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
-  do k=1,ex(3)-1
+  !$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
 ! x direction   
--- a/AMSS_NCKU_source/makefile
+++ b/AMSS_NCKU_source/makefile
@@ -95,8 +95,8 @@ $(CUDAFILES): bssn_gpu.h gpu_mem.h gpu_rhsSS_mem.h
 misc.o : zbesh.o

 # projects
-ABE: $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) 
-	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS)
+ABE: $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS)
+	$(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)
--- a/AMSS_NCKU_source/makefile.inc
+++ b/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 \
-               -fprofile-instr-use=$(PROFDATA) \
-               -Dfortran3 -Dnewc -I${MKLROOT}/include
-f90appflags  = -O3 -xHost -fp-model fast=2 -fma -ipo \
-               -fprofile-instr-use=$(PROFDATA) \
+CXXAPPFLAGS  = -O3 -march=native -fp-model fast=2 -fma -ipo -qopenmp \
+               -DMPI_STUB -Dfortran3 -Dnewc -I${MKLROOT}/include
+f90appflags  = -O3 -march=native -fp-model fast=2 -fma -ipo -qopenmp \
               -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
--- a/AMSS_NCKU_source/misc.C
+++ b/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"
--- a/AMSS_NCKU_source/misc.h
+++ b/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
 {
--- a/AMSS_NCKU_source/monitor.h
+++ b/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
 {
--- a/AMSS_NCKU_source/mpi_stub.h
+++ b/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 */
--- a/AMSS_NCKU_source/parameters.h
+++ b/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
 {
--- a/AMSS_NCKU_source/perf.h
+++ b/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
--- a/AMSS_NCKU_source/rungekutta4_rout.f90
+++ b/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

@@ -134,7 +144,7 @@

  endif

-  return   
+  return

  end subroutine rungekutta4_rout
 !-----------------------------------------------------------------------------
@@ -215,15 +225,19 @@

  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

-  return   
+  return

  end subroutine icn_rout
 !~~~~~~~~~~~~~~~~~~  
@@ -239,8 +253,10 @@
  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   
+  return

  end subroutine euler_rout
--- a/AMSS_NCKU_source/scalar_class.h
+++ b/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"
--- a/AMSS_NCKU_source/scalarwaves.C
+++ b/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"
--- a/AMSS_NCKU_source/setup.C
+++ b/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"
--- a/AMSS_NCKU_source/surface_integral.C
+++ b/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"
--- a/AMSS_NCKU_source/testNull.C
+++ b/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"
--- a/AMSS_NCKU_source/testNull2.C
+++ b/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"
--- a/AMSS_NCKU_source/var.C
+++ b/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"

--- a/makefile_and_run.py
+++ b/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


 ##################################################################
@@ -115,29 +128,28 @@ def run_ABE():
    print(                                                      )

    ## 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"
+        run_command         = NUMACTL_CPU_BIND + " ./ABE"
+        run_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_outfile = "ABEGPU_out.log"
- 
-    ## Execute the MPI command and stream output
-    mpi_process = subprocess.Popen(mpi_command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)
+        run_command         = NUMACTL_CPU_BIND + " ./ABEGPU"
+        run_command_outfile = "ABEGPU_out.log"
+
+    ## Execute the command and stream output
+    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 " )