Add OpenMP parallelization to BSSN RHS hot-path stencil routines

Enable OpenMP threading for the dominant computational kernels: - makefile.inc: add -qopenmp to f90appflags - diff_new.f90: split fderivs/fdderivs into OpenMP interior + serial boundary - kodiss.f90: split kodis into OpenMP interior + serial boundary - lopsidediff.f90: add OMP PARALLEL DO COLLAPSE(2) to lopsided - fmisc.f90: parallelize symmetry_bd bulk array copy - bssn_rhs.f90: add OMP WORKSHARE to array-syntax operations Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-02-07 13:58:55 +08:00
19 changed files with 598 additions and 1906 deletions
--- a/AMSS_NCKU_Program.py
+++ b/AMSS_NCKU_Program.py
@@ -8,14 +8,6 @@
 ##
 ##################################################################

-## Guard against re-execution by multiprocessing child processes.
-## Without this, using 'spawn' or 'forkserver' context would cause every
-## worker to re-run the entire script, spawning exponentially more
-## workers (fork bomb).
-if __name__ != '__main__':
-    import sys as _sys
-    _sys.exit(0)
-

 ##################################################################

@@ -432,31 +424,26 @@ print(

 import plot_xiaoqu
 import plot_GW_strain_amplitude_xiaoqu
-from parallel_plot_helper import run_plot_tasks_parallel
-
-plot_tasks = []

 ## Plot black hole trajectory
-plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot,   (binary_results_directory, figure_directory) ) )
-plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot3D, (binary_results_directory, figure_directory) ) )
+plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directory, figure_directory )
+plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directory, figure_directory )

 ## Plot black hole separation vs. time
-plot_tasks.append( ( plot_xiaoqu.generate_puncture_distence_plot, (binary_results_directory, figure_directory) ) )
+plot_xiaoqu.generate_puncture_distence_plot( binary_results_directory, figure_directory )

 ## Plot gravitational waveforms (psi4 and strain amplitude)
 for i in range(input_data.Detector_Number):
-    plot_tasks.append( ( plot_xiaoqu.generate_gravitational_wave_psi4_plot, (binary_results_directory, figure_directory, i) ) )
-    plot_tasks.append( ( plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot, (binary_results_directory, figure_directory, i) ) )
+    plot_xiaoqu.generate_gravitational_wave_psi4_plot( binary_results_directory, figure_directory, i )
+    plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot( binary_results_directory, figure_directory, i )

 ## Plot ADM mass evolution
 for i in range(input_data.Detector_Number):
-    plot_tasks.append( ( plot_xiaoqu.generate_ADMmass_plot, (binary_results_directory, figure_directory, i) ) )
+    plot_xiaoqu.generate_ADMmass_plot( binary_results_directory, figure_directory, i )

 ## Plot Hamiltonian constraint violation over time
 for i in range(input_data.grid_level):
-    plot_tasks.append( ( plot_xiaoqu.generate_constraint_check_plot, (binary_results_directory, figure_directory, i) ) )
-
-run_plot_tasks_parallel(plot_tasks)
+    plot_xiaoqu.generate_constraint_check_plot( binary_results_directory, figure_directory, i )

 ## Plot stored binary data
 plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )
--- a/AMSS_NCKU_source/Parallel.C
+++ b/AMSS_NCKU_source/Parallel.C
@@ -3756,358 +3756,6 @@ void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
  delete[] transfer_src;
  delete[] transfer_dst;
 }
-//
-// Async Sync: split into SyncBegin (initiate MPI) and SyncEnd (wait + unpack)
-// This allows overlapping MPI communication with computation.
-//
-static void transfer_begin(Parallel::TransferState *ts)
-{
-  int myrank;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  int cpusize = ts->cpusize;
-
-  ts->reqs = new MPI_Request[2 * cpusize];
-  ts->stats = new MPI_Status[2 * cpusize];
-  ts->req_no = 0;
-  ts->send_data = new double *[cpusize];
-  ts->rec_data = new double *[cpusize];
-  int length;
-
-  for (int node = 0; node < cpusize; node++)
-  {
-    ts->send_data[node] = ts->rec_data[node] = 0;
-    if (node == myrank)
-    {
-      // Local copy: pack then immediately unpack (no MPI needed)
-      if ((length = Parallel::data_packer(0, ts->transfer_src[myrank], ts->transfer_dst[myrank],
-                                          node, PACK, ts->VarList1, ts->VarList2, ts->Symmetry)))
-      {
-        double *local_data = new double[length];
-        if (!local_data)
-        {
-          cout << "out of memory in transfer_begin, local copy" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        Parallel::data_packer(local_data, ts->transfer_src[myrank], ts->transfer_dst[myrank],
-                              node, PACK, ts->VarList1, ts->VarList2, ts->Symmetry);
-        Parallel::data_packer(local_data, ts->transfer_src[node], ts->transfer_dst[node],
-                              node, UNPACK, ts->VarList1, ts->VarList2, ts->Symmetry);
-        delete[] local_data;
-      }
-    }
-    else
-    {
-      // send from this cpu to cpu#node
-      if ((length = Parallel::data_packer(0, ts->transfer_src[myrank], ts->transfer_dst[myrank],
-                                          node, PACK, ts->VarList1, ts->VarList2, ts->Symmetry)))
-      {
-        ts->send_data[node] = new double[length];
-        if (!ts->send_data[node])
-        {
-          cout << "out of memory in transfer_begin, send" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        Parallel::data_packer(ts->send_data[node], ts->transfer_src[myrank], ts->transfer_dst[myrank],
-                              node, PACK, ts->VarList1, ts->VarList2, ts->Symmetry);
-        MPI_Isend((void *)ts->send_data[node], length, MPI_DOUBLE, node, 1, MPI_COMM_WORLD,
-                  ts->reqs + ts->req_no++);
-      }
-      // receive from cpu#node to this cpu
-      if ((length = Parallel::data_packer(0, ts->transfer_src[node], ts->transfer_dst[node],
-                                          node, UNPACK, ts->VarList1, ts->VarList2, ts->Symmetry)))
-      {
-        ts->rec_data[node] = new double[length];
-        if (!ts->rec_data[node])
-        {
-          cout << "out of memory in transfer_begin, recv" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        MPI_Irecv((void *)ts->rec_data[node], length, MPI_DOUBLE, node, 1, MPI_COMM_WORLD,
-                  ts->reqs + ts->req_no++);
-      }
-    }
-  }
-  // NOTE: MPI_Waitall is NOT called here - that happens in transfer_end
-}
-//
-static void transfer_end(Parallel::TransferState *ts)
-{
-  // Wait for all pending MPI operations
-  MPI_Waitall(ts->req_no, ts->reqs, ts->stats);
-
-  // Unpack received data from remote ranks
-  for (int node = 0; node < ts->cpusize; node++)
-    if (ts->rec_data[node])
-      Parallel::data_packer(ts->rec_data[node], ts->transfer_src[node], ts->transfer_dst[node],
-                            node, UNPACK, ts->VarList1, ts->VarList2, ts->Symmetry);
-
-  // Cleanup MPI buffers
-  for (int node = 0; node < ts->cpusize; node++)
-  {
-    if (ts->send_data[node])
-      delete[] ts->send_data[node];
-    if (ts->rec_data[node])
-      delete[] ts->rec_data[node];
-  }
-  delete[] ts->reqs;
-  delete[] ts->stats;
-  delete[] ts->send_data;
-  delete[] ts->rec_data;
-}
-//
-Parallel::SyncHandle *Parallel::SyncBegin(Patch *Pat, MyList<var> *VarList, int Symmetry)
-{
-  int cpusize;
-  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
-
-  SyncHandle *handle = new SyncHandle;
-  handle->num_states = 1;
-  handle->states = new TransferState[1];
-
-  TransferState *ts = &handle->states[0];
-  ts->cpusize = cpusize;
-  ts->VarList1 = VarList;
-  ts->VarList2 = VarList;
-  ts->Symmetry = Symmetry;
-  ts->owns_gsl = true;
-
-  ts->dst = build_ghost_gsl(Pat);
-  ts->src = new MyList<Parallel::gridseg> *[cpusize];
-  ts->transfer_src = new MyList<Parallel::gridseg> *[cpusize];
-  ts->transfer_dst = new MyList<Parallel::gridseg> *[cpusize];
-  for (int node = 0; node < cpusize; node++)
-  {
-    ts->src[node] = build_owned_gsl0(Pat, node);
-    build_gstl(ts->src[node], ts->dst, &ts->transfer_src[node], &ts->transfer_dst[node]);
-  }
-
-  transfer_begin(ts);
-
-  return handle;
-}
-//
-Parallel::SyncHandle *Parallel::SyncBegin(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
-{
-  int cpusize;
-  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
-
-  // Count patches
-  int num_patches = 0;
-  MyList<Patch> *Pp = PatL;
-  while (Pp) { num_patches++; Pp = Pp->next; }
-
-  SyncHandle *handle = new SyncHandle;
-  handle->num_states = num_patches + 1; // intra-patch transfers + 1 inter-patch transfer
-  handle->states = new TransferState[handle->num_states];
-
-  // Intra-patch sync: for each patch, build ghost lists and initiate transfer
-  int idx = 0;
-  Pp = PatL;
-  while (Pp)
-  {
-    TransferState *ts = &handle->states[idx];
-    ts->cpusize = cpusize;
-    ts->VarList1 = VarList;
-    ts->VarList2 = VarList;
-    ts->Symmetry = Symmetry;
-    ts->owns_gsl = true;
-
-    ts->dst = build_ghost_gsl(Pp->data);
-    ts->src = new MyList<Parallel::gridseg> *[cpusize];
-    ts->transfer_src = new MyList<Parallel::gridseg> *[cpusize];
-    ts->transfer_dst = new MyList<Parallel::gridseg> *[cpusize];
-    for (int node = 0; node < cpusize; node++)
-    {
-      ts->src[node] = build_owned_gsl0(Pp->data, node);
-      build_gstl(ts->src[node], ts->dst, &ts->transfer_src[node], &ts->transfer_dst[node]);
-    }
-
-    transfer_begin(ts);
-
-    idx++;
-    Pp = Pp->next;
-  }
-
-  // Inter-patch sync: buffer zone exchange between patches
-  {
-    TransferState *ts = &handle->states[idx];
-    ts->cpusize = cpusize;
-    ts->VarList1 = VarList;
-    ts->VarList2 = VarList;
-    ts->Symmetry = Symmetry;
-    ts->owns_gsl = true;
-
-    ts->dst = build_buffer_gsl(PatL);
-    ts->src = new MyList<Parallel::gridseg> *[cpusize];
-    ts->transfer_src = new MyList<Parallel::gridseg> *[cpusize];
-    ts->transfer_dst = new MyList<Parallel::gridseg> *[cpusize];
-    for (int node = 0; node < cpusize; node++)
-    {
-      ts->src[node] = build_owned_gsl(PatL, node, 5, Symmetry);
-      build_gstl(ts->src[node], ts->dst, &ts->transfer_src[node], &ts->transfer_dst[node]);
-    }
-
-    transfer_begin(ts);
-  }
-
-  return handle;
-}
-//
-void Parallel::SyncEnd(SyncHandle *handle)
-{
-  if (!handle)
-    return;
-
-  // Wait for all pending transfers and unpack
-  for (int i = 0; i < handle->num_states; i++)
-  {
-    TransferState *ts = &handle->states[i];
-    transfer_end(ts);
-
-    // Cleanup grid segment lists only if this state owns them
-    if (ts->owns_gsl)
-    {
-      if (ts->dst)
-        ts->dst->destroyList();
-      for (int node = 0; node < ts->cpusize; node++)
-      {
-        if (ts->src[node])
-          ts->src[node]->destroyList();
-        if (ts->transfer_src[node])
-          ts->transfer_src[node]->destroyList();
-        if (ts->transfer_dst[node])
-          ts->transfer_dst[node]->destroyList();
-      }
-      delete[] ts->src;
-      delete[] ts->transfer_src;
-      delete[] ts->transfer_dst;
-    }
-  }
-
-  delete[] handle->states;
-  delete handle;
-}
-//
-// SyncPreparePlan: Pre-build grid segment lists for a patch list.
-// The plan can be reused across multiple SyncBeginWithPlan calls
-// as long as the mesh topology does not change (no regridding).
-//
-Parallel::SyncPlan *Parallel::SyncPreparePlan(MyList<Patch> *PatL, int Symmetry)
-{
-  int cpusize;
-  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
-
-  // Count patches
-  int num_patches = 0;
-  MyList<Patch> *Pp = PatL;
-  while (Pp) { num_patches++; Pp = Pp->next; }
-
-  SyncPlan *plan = new SyncPlan;
-  plan->num_entries = num_patches + 1; // intra-patch + 1 inter-patch
-  plan->Symmetry = Symmetry;
-  plan->entries = new SyncPlanEntry[plan->num_entries];
-
-  // Intra-patch entries: ghost zone exchange within each patch
-  int idx = 0;
-  Pp = PatL;
-  while (Pp)
-  {
-    SyncPlanEntry *pe = &plan->entries[idx];
-    pe->cpusize = cpusize;
-    pe->dst = build_ghost_gsl(Pp->data);
-    pe->src = new MyList<Parallel::gridseg> *[cpusize];
-    pe->transfer_src = new MyList<Parallel::gridseg> *[cpusize];
-    pe->transfer_dst = new MyList<Parallel::gridseg> *[cpusize];
-    for (int node = 0; node < cpusize; node++)
-    {
-      pe->src[node] = build_owned_gsl0(Pp->data, node);
-      build_gstl(pe->src[node], pe->dst, &pe->transfer_src[node], &pe->transfer_dst[node]);
-    }
-    idx++;
-    Pp = Pp->next;
-  }
-
-  // Inter-patch entry: buffer zone exchange between patches
-  {
-    SyncPlanEntry *pe = &plan->entries[idx];
-    pe->cpusize = cpusize;
-    pe->dst = build_buffer_gsl(PatL);
-    pe->src = new MyList<Parallel::gridseg> *[cpusize];
-    pe->transfer_src = new MyList<Parallel::gridseg> *[cpusize];
-    pe->transfer_dst = new MyList<Parallel::gridseg> *[cpusize];
-    for (int node = 0; node < cpusize; node++)
-    {
-      pe->src[node] = build_owned_gsl(PatL, node, 5, Symmetry);
-      build_gstl(pe->src[node], pe->dst, &pe->transfer_src[node], &pe->transfer_dst[node]);
-    }
-  }
-
-  return plan;
-}
-//
-void Parallel::SyncFreePlan(SyncPlan *plan)
-{
-  if (!plan)
-    return;
-
-  for (int i = 0; i < plan->num_entries; i++)
-  {
-    SyncPlanEntry *pe = &plan->entries[i];
-    if (pe->dst)
-      pe->dst->destroyList();
-    for (int node = 0; node < pe->cpusize; node++)
-    {
-      if (pe->src[node])
-        pe->src[node]->destroyList();
-      if (pe->transfer_src[node])
-        pe->transfer_src[node]->destroyList();
-      if (pe->transfer_dst[node])
-        pe->transfer_dst[node]->destroyList();
-    }
-    delete[] pe->src;
-    delete[] pe->transfer_src;
-    delete[] pe->transfer_dst;
-  }
-  delete[] plan->entries;
-  delete plan;
-}
-//
-// SyncBeginWithPlan: Use pre-built GSLs from a SyncPlan to initiate async transfer.
-// This avoids the O(cpusize * blocks^2) cost of rebuilding GSLs on every call.
-//
-Parallel::SyncHandle *Parallel::SyncBeginWithPlan(SyncPlan *plan, MyList<var> *VarList)
-{
-  return SyncBeginWithPlan(plan, VarList, VarList);
-}
-//
-Parallel::SyncHandle *Parallel::SyncBeginWithPlan(SyncPlan *plan, MyList<var> *VarList1, MyList<var> *VarList2)
-{
-  SyncHandle *handle = new SyncHandle;
-  handle->num_states = plan->num_entries;
-  handle->states = new TransferState[handle->num_states];
-
-  for (int i = 0; i < plan->num_entries; i++)
-  {
-    SyncPlanEntry *pe = &plan->entries[i];
-    TransferState *ts = &handle->states[i];
-
-    ts->cpusize = pe->cpusize;
-    ts->VarList1 = VarList1;
-    ts->VarList2 = VarList2;
-    ts->Symmetry = plan->Symmetry;
-    ts->owns_gsl = false; // GSLs are owned by the plan, not this handle
-
-    // Borrow GSL pointers from the plan (do NOT free them in SyncEnd)
-    ts->transfer_src = pe->transfer_src;
-    ts->transfer_dst = pe->transfer_dst;
-    ts->src = pe->src;
-    ts->dst = pe->dst;
-
-    transfer_begin(ts);
-  }
-
-  return handle;
-}
 // collect buffer grid segments or blocks for the periodic boundary condition of given patch
 // ---------------------------------------------------
 // |con |                                       |con |
--- a/AMSS_NCKU_source/Parallel.h
+++ b/AMSS_NCKU_source/Parallel.h
@@ -81,53 +81,6 @@ namespace Parallel
                   int Symmetry);
  void Sync(Patch *Pat, MyList<var> *VarList, int Symmetry);
  void Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry);
-
-  // Async Sync: overlap MPI communication with computation
-  struct TransferState
-  {
-    MPI_Request *reqs;
-    MPI_Status *stats;
-    int req_no;
-    double **send_data;
-    double **rec_data;
-    int cpusize;
-    MyList<gridseg> **transfer_src;
-    MyList<gridseg> **transfer_dst;
-    MyList<gridseg> **src;
-    MyList<gridseg> *dst;
-    MyList<var> *VarList1;
-    MyList<var> *VarList2;
-    int Symmetry;
-    bool owns_gsl; // true if this state owns and should free the GSLs
-  };
-  struct SyncHandle
-  {
-    TransferState *states;
-    int num_states;
-  };
-  SyncHandle *SyncBegin(Patch *Pat, MyList<var> *VarList, int Symmetry);
-  SyncHandle *SyncBegin(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry);
-  void SyncEnd(SyncHandle *handle);
-
-  // Cached GSL plan: pre-build grid segment lists once, reuse across multiple Sync calls
-  struct SyncPlanEntry
-  {
-    int cpusize;
-    MyList<gridseg> **transfer_src;
-    MyList<gridseg> **transfer_dst;
-    MyList<gridseg> **src;
-    MyList<gridseg> *dst;
-  };
-  struct SyncPlan
-  {
-    SyncPlanEntry *entries;
-    int num_entries;
-    int Symmetry;
-  };
-  SyncPlan *SyncPreparePlan(MyList<Patch> *PatL, int Symmetry);
-  void SyncFreePlan(SyncPlan *plan);
-  SyncHandle *SyncBeginWithPlan(SyncPlan *plan, MyList<var> *VarList);
-  SyncHandle *SyncBeginWithPlan(SyncPlan *plan, MyList<var> *VarList1, MyList<var> *VarList2);
  void OutBdLow2Hi(Patch *Patc, Patch *Patf,
                   MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
                   int Symmetry);
--- a/AMSS_NCKU_source/TwoPunctures.C
+++ b/AMSS_NCKU_source/TwoPunctures.C
--- a/AMSS_NCKU_source/TwoPunctures.h
+++ b/AMSS_NCKU_source/TwoPunctures.h
@@ -1,8 +1,7 @@
+
 #ifndef TWO_PUNCTURES_H
 #define TWO_PUNCTURES_H

-#include <omp.h>
-
 #define StencilSize 19
 #define N_PlaneRelax 1
 #define NRELAX 200
@@ -43,17 +42,32 @@ private:

       int ntotal;

-       // ===== Precomputed spectral derivative matrices =====
-       double *D1_A, *D2_A;
-       double *D1_B, *D2_B;
-       double *DF1_phi, *DF2_phi;
-
-       // ===== Pre-allocated workspace for LineRelax (per-thread) =====
-       int max_threads;
-       double **ws_diag_be, **ws_e_be, **ws_f_be, **ws_b_be, **ws_x_be;
-       double **ws_l_be, **ws_u_be, **ws_d_be, **ws_y_be;
-       double **ws_diag_al, **ws_e_al, **ws_f_al, **ws_b_al, **ws_x_al;
-       double **ws_l_al, **ws_u_al, **ws_d_al, **ws_y_al;
+       // 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
       {
@@ -71,28 +85,6 @@ public:
                    int Newtonmaxit);
       ~TwoPunctures();

-       // 02/07: New/modified methods
-       void allocate_workspace();
-       void free_workspace();
-       void precompute_derivative_matrices();
-       void build_cheb_deriv_matrices(int n, double *D1, double *D2);
-       void build_fourier_deriv_matrices(int N, double *DF1, double *DF2);
-       void Derivatives_AB3_MatMul(int nvar, int n1, int n2, int n3, derivs v);
-       void ThomasAlgorithm_ws(int N, double *b, double *a, double *c, double *x, double *q,
-                                double *l, double *u_ws, double *d, double *y);
-       void LineRelax_be_omp(double *dv,
-                         int const i, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols, int **cols,
-                         double **JFD, int tid);
-       void LineRelax_al_omp(double *dv,
-                         int const j, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols,
-                         int **cols, double **JFD, int tid);
-       void relax_omp(double *dv, int const nvar, int const n1, int const n2, int const n3,
-                  double const *rhs, int const *ncols, int **cols, double **JFD);
-
       void Solve();
       void set_initial_guess(derivs v);
       int index(int i, int j, int k, int l, int a, int b, int c, int d);
@@ -151,11 +143,23 @@ public:
       double BY_KKofxyz(double x, double y, double z);
       void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix);
       void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u);
+       void relax(double *dv, int const nvar, int const n1, int const n2, int const n3,
+                  double const *rhs, int const *ncols, int **cols, double **JFD);
+       void LineRelax_be(double *dv,
+                         int const i, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols, int **cols,
+                         double **JFD);
       void JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
                         int n3, derivs dv, derivs u, double *values);
       void LinEquations(double A, double B, double X, double R,
                         double x, double r, double phi,
                         double y, double z, derivs dU, derivs U, double *values);
+       void LineRelax_al(double *dv,
+                         int const j, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols,
+                         int **cols, double **JFD);
       void ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q);
       void Save(char *fname);
       // provided by Vasileios Paschalidis (vpaschal@illinois.edu)
--- a/AMSS_NCKU_source/Z4c_class.C
+++ b/AMSS_NCKU_source/Z4c_class.C
@@ -186,12 +186,6 @@ void Z4c_class::Step(int lev, int YN)
  int ERROR = 0;

  MyList<ss_patch> *sPp;
-
-  // Pre-build grid segment lists once for this level's patches.
-  // These are reused across predictor + 3 corrector SyncBegin calls,
-  // avoiding O(cpusize * blocks^2) rebuild each time.
-  Parallel::SyncPlan *sync_plan = Parallel::SyncPreparePlan(GH->PatL[lev], Symmetry);
-
  // Predictor
  MyList<Patch> *Pp = GH->PatL[lev];
  while (Pp)
@@ -327,17 +321,13 @@ void Z4c_class::Step(int lev, int YN)
    }
    Pp = Pp->next;
  }
-  // Start async ghost zone exchange - overlaps with error check and Shell computation
-  Parallel::SyncHandle *sync_pre = Parallel::SyncBeginWithPlan(sync_plan, SynchList_pre);
-
-  // check error information (overlaps with MPI transfer)
+  // check error information
  {
    int erh = ERROR;
    MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
  }
  if (ERROR)
  {
-    Parallel::SyncEnd(sync_pre); sync_pre = 0;
    Parallel::Dump_Data(GH->PatL[lev], StateList, 0, PhysTime, dT_lev);
    if (myrank == 0)
    {
@@ -485,7 +475,6 @@ void Z4c_class::Step(int lev, int YN)
  }
  if (ERROR)
  {
-    Parallel::SyncEnd(sync_pre); sync_pre = 0;
    SH->Dump_Data(StateList, 0, PhysTime, dT_lev);
    if (myrank == 0)
    {
@@ -496,8 +485,7 @@ void Z4c_class::Step(int lev, int YN)
  }
 #endif

-  // Complete async ghost zone exchange
-  if (sync_pre) Parallel::SyncEnd(sync_pre);
+  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);

 #ifdef WithShell
  if (lev == 0)
@@ -705,17 +693,13 @@ void Z4c_class::Step(int lev, int YN)
      Pp = Pp->next;
    }

-    // Start async ghost zone exchange - overlaps with error check and Shell computation
-    Parallel::SyncHandle *sync_cor = Parallel::SyncBeginWithPlan(sync_plan, SynchList_cor);
-
-    // check error information (overlaps with MPI transfer)
+    // check error information
    {
      int erh = ERROR;
      MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    }
    if (ERROR)
    {
-      Parallel::SyncEnd(sync_cor); sync_cor = 0;
      Parallel::Dump_Data(GH->PatL[lev], SynchList_pre, 0, PhysTime, dT_lev);
      if (myrank == 0)
      {
@@ -873,7 +857,6 @@ void Z4c_class::Step(int lev, int YN)
    }
    if (ERROR)
    {
-      Parallel::SyncEnd(sync_cor); sync_cor = 0;
      SH->Dump_Data(SynchList_pre, 0, PhysTime, dT_lev);
      if (myrank == 0)
      {
@@ -885,8 +868,7 @@ void Z4c_class::Step(int lev, int YN)
    }
 #endif

-    // Complete async ghost zone exchange
-    if (sync_cor) Parallel::SyncEnd(sync_cor);
+    Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);

 #ifdef WithShell
    if (lev == 0)
@@ -1060,8 +1042,6 @@ void Z4c_class::Step(int lev, int YN)
      Porg0[ithBH][2] = Porg1[ithBH][2];
    }
  }
-
-  Parallel::SyncFreePlan(sync_plan);
 }
 #else
 // for constraint preserving boundary (CPBC)
@@ -1095,10 +1075,6 @@ void Z4c_class::Step(int lev, int YN)
  int ERROR = 0;

  MyList<ss_patch> *sPp;
-
-  // Pre-build grid segment lists once for this level's patches.
-  Parallel::SyncPlan *sync_plan = Parallel::SyncPreparePlan(GH->PatL[lev], Symmetry);
-
  // Predictor
  MyList<Patch> *Pp = GH->PatL[lev];
  while (Pp)
@@ -1566,17 +1542,13 @@ void Z4c_class::Step(int lev, int YN)
  }
 #endif
  }
-  // Start async ghost zone exchange - overlaps with error check
-  Parallel::SyncHandle *sync_pre = Parallel::SyncBeginWithPlan(sync_plan, SynchList_pre);
-
-  // check error information (overlaps with MPI transfer)
+  // check error information
  {
    int erh = ERROR;
    MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
  }
  if (ERROR)
  {
-    Parallel::SyncEnd(sync_pre); sync_pre = 0;
    SH->Dump_Data(StateList, 0, PhysTime, dT_lev);
    if (myrank == 0)
    {
@@ -1586,8 +1558,7 @@ void Z4c_class::Step(int lev, int YN)
    }
  }

-  // Complete async ghost zone exchange
-  if (sync_pre) Parallel::SyncEnd(sync_pre);
+  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);

  if (lev == 0)
  {
@@ -2132,17 +2103,13 @@ void Z4c_class::Step(int lev, int YN)
        sPp = sPp->next;
      }
    }
-    // Start async ghost zone exchange - overlaps with error check
-    Parallel::SyncHandle *sync_cor = Parallel::SyncBeginWithPlan(sync_plan, SynchList_cor);
-
-    // check error information (overlaps with MPI transfer)
+    // check error information
    {
      int erh = ERROR;
      MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    }
    if (ERROR)
    {
-      Parallel::SyncEnd(sync_cor); sync_cor = 0;
      SH->Dump_Data(SynchList_pre, 0, PhysTime, dT_lev);
      if (myrank == 0)
      {
@@ -2153,8 +2120,7 @@ void Z4c_class::Step(int lev, int YN)
      }
    }

-    // Complete async ghost zone exchange
-    if (sync_cor) Parallel::SyncEnd(sync_cor);
+    Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);

    if (lev == 0)
    {
@@ -2380,8 +2346,6 @@ void Z4c_class::Step(int lev, int YN)
 	  DG_List->clearList();
 	}
 #endif
-
-  Parallel::SyncFreePlan(sync_plan);
 }
 #endif
 #undef MRBD
--- a/AMSS_NCKU_source/bssn_class.C
+++ b/AMSS_NCKU_source/bssn_class.C
@@ -3035,12 +3035,6 @@ void bssn_class::Step(int lev, int YN)
  int ERROR = 0;

  MyList<ss_patch> *sPp;
-
-  // Pre-build grid segment lists once for this level's patches.
-  // These are reused across predictor + 3 corrector SyncBegin calls,
-  // avoiding O(cpusize * blocks^2) rebuild each time.
-  Parallel::SyncPlan *sync_plan = Parallel::SyncPreparePlan(GH->PatL[lev], Symmetry);
-
  // Predictor
  MyList<Patch> *Pp = GH->PatL[lev];
  while (Pp)
@@ -3164,18 +3158,13 @@ void bssn_class::Step(int lev, int YN)
    }
    Pp = Pp->next;
  }
-
-  // Start async ghost zone exchange - overlaps with error check and Shell computation
-  Parallel::SyncHandle *sync_pre = Parallel::SyncBeginWithPlan(sync_plan, SynchList_pre);
-
-  // check error information (overlaps with MPI transfer)
+  // check error information
  {
    int erh = ERROR;
    MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
  }
  if (ERROR)
  {
-    Parallel::SyncEnd(sync_pre); sync_pre = 0;
    Parallel::Dump_Data(GH->PatL[lev], StateList, 0, PhysTime, dT_lev);
    if (myrank == 0)
    {
@@ -3335,7 +3324,6 @@ void bssn_class::Step(int lev, int YN)

  if (ERROR)
  {
-    Parallel::SyncEnd(sync_pre); sync_pre = 0;
    SH->Dump_Data(StateList, 0, PhysTime, dT_lev);
    if (myrank == 0)
    {
@@ -3346,8 +3334,7 @@ void bssn_class::Step(int lev, int YN)
  }
 #endif

-  // Complete async ghost zone exchange
-  if (sync_pre) Parallel::SyncEnd(sync_pre);
+  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);

 #ifdef WithShell
  if (lev == 0)
@@ -3541,10 +3528,7 @@ void bssn_class::Step(int lev, int YN)
      Pp = Pp->next;
    }

-    // Start async ghost zone exchange - overlaps with error check and Shell computation
-    Parallel::SyncHandle *sync_cor = Parallel::SyncBeginWithPlan(sync_plan, SynchList_cor);
-
-    // check error information (overlaps with MPI transfer)
+    // check error information
    {
      int erh = ERROR;
      MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
@@ -3552,7 +3536,6 @@ void bssn_class::Step(int lev, int YN)

    if (ERROR)
    {
-      Parallel::SyncEnd(sync_cor); sync_cor = 0;
      Parallel::Dump_Data(GH->PatL[lev], SynchList_pre, 0, PhysTime, dT_lev);
      if (myrank == 0)
      {
@@ -3709,7 +3692,6 @@ void bssn_class::Step(int lev, int YN)
    }
    if (ERROR)
    {
-      Parallel::SyncEnd(sync_cor); sync_cor = 0;
      SH->Dump_Data(SynchList_pre, 0, PhysTime, dT_lev);
      if (myrank == 0)
      {
@@ -3722,8 +3704,7 @@ void bssn_class::Step(int lev, int YN)
    }
 #endif

-    // Complete async ghost zone exchange
-    if (sync_cor) Parallel::SyncEnd(sync_cor);
+    Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);

 #ifdef WithShell
    if (lev == 0)
@@ -3914,8 +3895,6 @@ void bssn_class::Step(int lev, int YN)
      Porg0[ithBH][2] = Porg1[ithBH][2];
    }
  }
-
-  Parallel::SyncFreePlan(sync_plan);
 }

 //================================================================================================
@@ -4838,12 +4817,6 @@ void bssn_class::Step(int lev, int YN)
  int ERROR = 0;

  MyList<ss_patch> *sPp;
-
-  // Pre-build grid segment lists once for this level's patches.
-  // These are reused across predictor + 3 corrector SyncBegin calls,
-  // avoiding O(cpusize * blocks^2) rebuild each time.
-  Parallel::SyncPlan *sync_plan = Parallel::SyncPreparePlan(GH->PatL[lev], Symmetry);
-
  // Predictor
  MyList<Patch> *Pp = GH->PatL[lev];
  while (Pp)
@@ -4970,17 +4943,13 @@ void bssn_class::Step(int lev, int YN)

  //   misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after Predictor rhs calculation");

-  // Start async ghost zone exchange - overlaps with error check and BH position
-  Parallel::SyncHandle *sync_pre = Parallel::SyncBeginWithPlan(sync_plan, SynchList_pre);
-
-  // check error information (overlaps with MPI transfer)
+  // check error information
  {
    int erh = ERROR;
    MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, GH->Commlev[lev]);
  }
  if (ERROR)
  {
-    Parallel::SyncEnd(sync_pre); sync_pre = 0;
    Parallel::Dump_Data(GH->PatL[lev], StateList, 0, PhysTime, dT_lev);
    if (myrank == 0)
    {
@@ -4992,8 +4961,7 @@ void bssn_class::Step(int lev, int YN)

  //   misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Predictor sync");

-  // Complete async ghost zone exchange
-  if (sync_pre) Parallel::SyncEnd(sync_pre);
+  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);

 #if (MAPBH == 0)
  // for black hole position
@@ -5172,17 +5140,13 @@ void bssn_class::Step(int lev, int YN)

    //   misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Corrector error check");

-    // Start async ghost zone exchange - overlaps with error check and BH position
-    Parallel::SyncHandle *sync_cor = Parallel::SyncBeginWithPlan(sync_plan, SynchList_cor);
-
-    // check error information (overlaps with MPI transfer)
+    // check error information
    {
      int erh = ERROR;
      MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, GH->Commlev[lev]);
    }
    if (ERROR)
    {
-      Parallel::SyncEnd(sync_cor); sync_cor = 0;
      Parallel::Dump_Data(GH->PatL[lev], SynchList_pre, 0, PhysTime, dT_lev);
      if (myrank == 0)
      {
@@ -5196,8 +5160,7 @@ void bssn_class::Step(int lev, int YN)

    //    misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Corrector sync");

-    // Complete async ghost zone exchange
-    if (sync_cor) Parallel::SyncEnd(sync_cor);
+    Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);

    //    misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after Corrector sync");

@@ -5313,8 +5276,6 @@ void bssn_class::Step(int lev, int YN)

  //     if(myrank==GH->start_rank[lev]) cout<<GH->mylev<<endl;
  //     misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"complet GH Step");
-
-  Parallel::SyncFreePlan(sync_plan);
 }

 //================================================================================================
--- a/AMSS_NCKU_source/bssn_rhs.f90
+++ b/AMSS_NCKU_source/bssn_rhs.f90
@@ -168,6 +168,8 @@
  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)

@@ -201,6 +203,8 @@
  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
@@ -234,6 +238,8 @@
  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 ))
@@ -282,6 +288,8 @@
          (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)
@@ -336,6 +344,8 @@
  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    ) + &
@@ -375,6 +385,8 @@
  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)
@@ -401,6 +413,8 @@
   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  + &
@@ -601,9 +615,13 @@
            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
@@ -626,11 +644,15 @@
  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
@@ -791,6 +813,8 @@
 !!!! 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
@@ -945,60 +969,103 @@
  SSA(2)=SYM
  SSA(3)=ANTI

-!!!!!!!!!advection term + Kreiss-Oliger dissipation (merged for cache efficiency)
-! lopsided_kodis shares the symmetry_bd buffer between advection and
-! dissipation, eliminating redundant full-grid copies. For metric variables
-! gxx/gyy/gzz (=dxx/dyy/dzz+1): kodis stencil coefficients sum to zero,
-! so the constant offset has no effect on dissipation.
+!!!!!!!!!advection term part

-  call lopsided_kodis(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
-  call lopsided_kodis(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,eps)
-  call lopsided_kodis(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,eps)
-  call lopsided_kodis(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
-  call lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps)
-  call lopsided_kodis(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS)
+  call lopsided(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA)
+  call lopsided(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA)
+  call lopsided(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS)

-  call lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
-  call lopsided_kodis(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,eps)
-  call lopsided_kodis(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA,eps)
-  call lopsided_kodis(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
-  call lopsided_kodis(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA,eps)
-  call lopsided_kodis(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS)
+  call lopsided(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA)
+  call lopsided(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA)
+  call lopsided(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS)

-  call lopsided_kodis(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS,eps)
-  call lopsided_kodis(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS)

-  call lopsided_kodis(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS,eps)
-  call lopsided_kodis(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS,eps)
-  call lopsided_kodis(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA,eps)
-
-#if 1 
-!! bam does not apply dissipation on gauge variables
-  call lopsided_kodis(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS,eps)
-#if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7)
-  call lopsided_kodis(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS,eps)
-  call lopsided_kodis(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS,eps)
-  call lopsided_kodis(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,eps)
-#endif
-#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
-  call lopsided_kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,eps)
-  call lopsided_kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,eps)
-  call lopsided_kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA,eps)
-#endif
-#else
-! No dissipation on gauge variables (advection only)
+  call lopsided(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS)
+  call lopsided(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS)
+  call lopsided(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA)
+!!
  call lopsided(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS)
+
 #if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7)
  call lopsided(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS)
  call lopsided(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS)
  call lopsided(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA)
 #endif
+
 #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
  call lopsided(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS)
  call lopsided(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS)
  call lopsided(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA)
 #endif
+
+  if(eps>0)then 
+! usual Kreiss-Oliger dissipation      
+  call kodis(ex,X,Y,Z,chi,chi_rhs,SSS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,trK,trK_rhs,SSS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,dxx,gxx_rhs,SSS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,gxy,gxy_rhs,AAS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,gxz,gxz_rhs,ASA,Symmetry,eps)
+  call kodis(ex,X,Y,Z,dyy,gyy_rhs,SSS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,gyz,gyz_rhs,SAA,Symmetry,eps)
+  call kodis(ex,X,Y,Z,dzz,gzz_rhs,SSS,Symmetry,eps)
+#if 0
+#define i 42
+#define j 40
+#define k 40
+if(Lev == 1)then
+write(*,*) X(i),Y(j),Z(k)
+write(*,*) "before",Axx_rhs(i,j,k)
+endif
+#undef i
+#undef j
+#undef k
+!!stop
 #endif
+  call kodis(ex,X,Y,Z,Axx,Axx_rhs,SSS,Symmetry,eps)
+#if 0
+#define i 42
+#define j 40
+#define k 40
+if(Lev == 1)then
+write(*,*) X(i),Y(j),Z(k)
+write(*,*) "after",Axx_rhs(i,j,k)
+endif
+#undef i
+#undef j
+#undef k
+!!stop
+#endif
+  call kodis(ex,X,Y,Z,Axy,Axy_rhs,AAS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,Axz,Axz_rhs,ASA,Symmetry,eps)
+  call kodis(ex,X,Y,Z,Ayy,Ayy_rhs,SSS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,Ayz,Ayz_rhs,SAA,Symmetry,eps)
+  call kodis(ex,X,Y,Z,Azz,Azz_rhs,SSS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,Gamx,Gamx_rhs,ASS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,Gamy,Gamy_rhs,SAS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,Gamz,Gamz_rhs,SSA,Symmetry,eps)
+
+#if 1 
+!! bam does not apply dissipation on gauge variables
+  call kodis(ex,X,Y,Z,Lap,Lap_rhs,SSS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,betax,betax_rhs,ASS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,betay,betay_rhs,SAS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,betaz,betaz_rhs,SSA,Symmetry,eps)
+#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
+  call kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,ASS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,SAS,Symmetry,eps)
+  call kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,SSA,Symmetry,eps)
+#endif
+#endif
+
+  endif

  if(co == 0)then
 ! ham_Res = trR + 2/3 * K^2 - A_ij * A^ij - 16 * PI * rho
--- a/AMSS_NCKU_source/diff_new.f90
+++ b/AMSS_NCKU_source/diff_new.f90
@@ -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,9 +1040,13 @@

 ! set kmax and kmin 0
    endif
+  enddo
+  enddo
+  enddo
 #elif 0
-! x direction   
-        if(i+2 <= imax .and. i-2 >= imin)then
+  do k=1,ex(3)-1
+  do j=1,ex(2)-1
+  do i=1,ex(1)-1
 !
 !              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
 !  fx(i) = ---------------------------------------------
@@ -1079,8 +1083,32 @@

 ! set kmax and kmin 0
    endif
+  enddo
+  enddo
+  enddo
 #else
-! for bam comparison
+! 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
   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
@@ -1094,10 +1122,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

@@ -1401,10 +1429,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
 !
@@ -1482,8 +1510,47 @@
   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
+! 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
   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
@@ -1518,10 +1585,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

--- a/AMSS_NCKU_source/fmisc.f90
+++ b/AMSS_NCKU_source/fmisc.f90
@@ -881,9 +881,18 @@ 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
+  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

-  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
--- a/AMSS_NCKU_source/kodiss.f90
+++ b/AMSS_NCKU_source/kodiss.f90
@@ -159,36 +159,12 @@ integer, parameter :: NO_SYMM=0, OCTANT=2

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

-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-
-  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 ?
+! 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)) + &
@@ -204,9 +180,37 @@ 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
+  !$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
+   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 )
+  endif
  enddo
  enddo
  enddo
--- a/AMSS_NCKU_source/lopsidediff.f90
+++ b/AMSS_NCKU_source/lopsidediff.f90
@@ -233,6 +233,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) SCHEDULE(static) PRIVATE(i,j,k)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
@@ -482,206 +483,12 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)
  enddo
  enddo
  enddo
+  !$OMP END PARALLEL DO

  return

  end subroutine lopsided

-!-----------------------------------------------------------------------------
-! Combined advection (lopsided) + Kreiss-Oliger dissipation (kodis)
-! Shares the symmetry_bd buffer fh, eliminating one full-grid copy per call.
-! Mathematically identical to calling lopsided then kodis separately.
-!-----------------------------------------------------------------------------
-subroutine lopsided_kodis(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA,eps)
-  implicit none
-
-!~~~~~~> Input parameters:
-
-  integer, intent(in)  :: ex(1:3),Symmetry
-  real*8,  intent(in)  :: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3))
-  real*8,dimension(ex(1),ex(2),ex(3)),intent(in)   :: f,Sfx,Sfy,Sfz
-
-  real*8,dimension(ex(1),ex(2),ex(3)),intent(inout):: f_rhs
-  real*8,dimension(3),intent(in) ::SoA
-  real*8,intent(in) :: eps
-
-!~~~~~~> local variables:
-! note index -2,-1,0, so we have 3 extra points
-  real*8,dimension(-2:ex(1),-2:ex(2),-2:ex(3))   :: fh
-  integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
-  real*8 :: dX,dY,dZ
-  real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
-  real*8,  parameter :: ZEO=0.d0,ONE=1.d0, F3=3.d0
-  real*8,  parameter :: TWO=2.d0,F6=6.0d0,F18=1.8d1
-  real*8,  parameter :: F12=1.2d1, F10=1.d1,EIT=8.d0
-  integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
-! kodis parameters
-  real*8, parameter :: SIX=6.d0,FIT=1.5d1,TWT=2.d1
-  real*8, parameter :: cof=6.4d1   ! 2^6
-
-  dX = X(2)-X(1)
-  dY = Y(2)-Y(1)
-  dZ = Z(2)-Z(1)
-
-  d12dx = ONE/F12/dX
-  d12dy = ONE/F12/dY
-  d12dz = ONE/F12/dZ
-
-  d2dx = ONE/TWO/dX
-  d2dy = ONE/TWO/dY
-  d2dz = ONE/TWO/dZ
-
-  imax = ex(1)
-  jmax = ex(2)
-  kmax = ex(3)
-
-  imin = 1
-  jmin = 1
-  kmin = 1
-  if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -2
-  if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -2
-  if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -2
-
-! Single symmetry_bd call shared by both advection and dissipation
-  call symmetry_bd(3,ex,f,fh,SoA)
-
-! ---- Advection (lopsided) loop ----
-! upper bound set ex-1 only for efficiency, 
-! the loop body will set ex 0 also
-  do k=1,ex(3)-1
-  do j=1,ex(2)-1
-  do i=1,ex(1)-1
-! x direction   
-    if(Sfx(i,j,k) > ZEO)then
-      if(i+3 <= imax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
-                  Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) &
-                                    -F6*fh(i+2,j,k)+    fh(i+3,j,k))
-     elseif(i+2 <= imax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
-                  Sfx(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))
-
-     elseif(i+1 <= imax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
-                  Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) &
-                                    -F6*fh(i-2,j,k)+    fh(i-3,j,k))
-     endif
-   elseif(Sfx(i,j,k) < ZEO)then
-      if(i-3 >= imin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
-                  Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) &
-                                    -F6*fh(i-2,j,k)+    fh(i-3,j,k))
-     elseif(i-2 >= imin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
-                  Sfx(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))
-
-     elseif(i-1 >= imin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
-                  Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) &
-                                    -F6*fh(i+2,j,k)+    fh(i+3,j,k))
-     endif
-   endif
-
-! y direction   
-    if(Sfy(i,j,k) > ZEO)then
-      if(j+3 <= jmax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
-                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) &
-                                    -F6*fh(i,j+2,k)+    fh(i,j+3,k))
-     elseif(j+2 <= jmax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
-                  Sfy(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))
-
-     elseif(j+1 <= jmax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
-                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) &
-                                    -F6*fh(i,j-2,k)+    fh(i,j-3,k))
-     endif
-   elseif(Sfy(i,j,k) < ZEO)then
-      if(j-3 >= jmin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
-                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) &
-                                    -F6*fh(i,j-2,k)+    fh(i,j-3,k))
-     elseif(j-2 >= jmin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
-                  Sfy(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))
-
-     elseif(j-1 >= jmin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
-                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) &
-                                    -F6*fh(i,j+2,k)+    fh(i,j+3,k))
-     endif
-   endif
-
-! z direction   
-    if(Sfz(i,j,k) > ZEO)then
-      if(k+3 <= kmax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
-                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) &
-                                    -F6*fh(i,j,k+2)+    fh(i,j,k+3))
-     elseif(k+2 <= kmax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
-                  Sfz(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))
-
-     elseif(k+1 <= kmax)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
-                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) &
-                                    -F6*fh(i,j,k-2)+    fh(i,j,k-3))
-     endif
-   elseif(Sfz(i,j,k) < ZEO)then
-      if(k-3 >= kmin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
-                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) &
-                                    -F6*fh(i,j,k-2)+    fh(i,j,k-3))
-     elseif(k-2 >= kmin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
-                  Sfz(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))
-
-     elseif(k-1 >= kmin)then
-     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
-                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) &
-                                    -F6*fh(i,j,k+2)+    fh(i,j,k+3))
-     endif
-   endif
-  enddo
-  enddo
-  enddo
-
-! ---- Dissipation (kodis) loop ----
-  if(eps > ZEO) then
-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-
-  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
-   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 )
-  endif
-
-  enddo
-  enddo
-  enddo
-  endif
-
-  return
-
-  end subroutine lopsided_kodis
-
 #elif (ghost_width == 4)
 ! sixth order code
 ! Compute advection terms in right hand sides of field equations
--- a/AMSS_NCKU_source/makefile
+++ b/AMSS_NCKU_source/makefile
@@ -16,12 +16,6 @@ include makefile.inc
 .cu.o:
 	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)

-TwoPunctures.o: TwoPunctures.C
-	${CXX} $(CXXAPPFLAGS) -qopenmp -c $< -o $@
-
-TwoPunctureABE.o: TwoPunctureABE.C
-	${CXX} $(CXXAPPFLAGS) -qopenmp -c $< -o $@
-
 # Input files
 C++FILES = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o\
           cgh.o bssn_class.o surface_integral.o ShellPatch.o\
@@ -102,7 +96,7 @@ ABEGPU: $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
 	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)

 TwoPunctureABE: $(TwoPunctureFILES)
-	$(CLINKER) $(CXXAPPFLAGS) -qopenmp -o $@ $(TwoPunctureFILES) $(LDLIBS)
+	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(TwoPunctureFILES) $(LDLIBS)

 clean:
 	rm *.o ABE ABEGPU TwoPunctureABE make.log -f
--- a/AMSS_NCKU_source/makefile.inc
+++ b/AMSS_NCKU_source/makefile.inc
@@ -8,16 +8,17 @@ filein  = -I/usr/include/ -I${MKLROOT}/include

 ## Using sequential MKL (OpenMP disabled for better single-threaded performance)
 ## Added -lifcore for Intel Fortran runtime and -limf for Intel math library
-LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl
+LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lifcore -limf -lpthread -lm -ldl -qopenmp

 ## 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
-CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma -ipo \
+## 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 \
+f90appflags  = -O3 -xHost -fp-model fast=2 -fma -ipo -qopenmp \
               -align array64byte -fpp -I${MKLROOT}/include
 f90          = ifx
 f77          = ifx
--- a/makefile_and_run.py
+++ b/makefile_and_run.py
@@ -10,7 +10,7 @@

 import AMSS_NCKU_Input as input_data
 import subprocess
-import time
+
 ## 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)
@@ -152,7 +152,7 @@ def run_ABE():
 ## Run the AMSS-NCKU TwoPuncture program TwoPunctureABE

 def run_TwoPunctureABE():
-    tp_time1=time.time()
+
    print(                                                          )
    print( " Running the AMSS-NCKU executable file TwoPunctureABE " ) 
    print(                                                          )
@@ -179,9 +179,7 @@ def run_TwoPunctureABE():
    print(                                               )
    print( " The TwoPunctureABE simulation is finished " ) 
    print(                                               )
-    tp_time2=time.time()
-    et=tp_time2-tp_time1
-    print(f"Used time: {et}")
+    
    return

 ##################################################################
--- a/parallel_plot_helper.py
+++ b/parallel_plot_helper.py
@@ -1,29 +0,0 @@
-import multiprocessing
-
-def run_plot_task(task):
-    """Execute a single plotting task.
-    
-    Parameters
-    ----------
-    task : tuple
-        A tuple of (function, args_tuple) where function is a callable
-        plotting function and args_tuple contains its arguments.
-    """
-    func, args = task
-    return func(*args)
-
-
-def run_plot_tasks_parallel(plot_tasks):
-    """Execute a list of independent plotting tasks in parallel.
-
-    Uses the 'fork' context to create worker processes so that the main
-    script is NOT re-imported/re-executed in child processes.
-
-    Parameters
-    ----------
-    plot_tasks : list of tuples
-        Each element is (function, args_tuple).
-    """
-    ctx = multiprocessing.get_context('fork')
-    with ctx.Pool() as pool:
-        pool.map(run_plot_task, plot_tasks)
--- a/plot_GW_strain_amplitude_xiaoqu.py
+++ b/plot_GW_strain_amplitude_xiaoqu.py
@@ -11,8 +11,6 @@
 import numpy                               ## numpy for array operations
 import scipy                               ## scipy for interpolation and signal processing
 import math
-import matplotlib
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt     ## matplotlib for plotting
 import os                                  ## os for system/file operations

--- a/plot_binary_data.py
+++ b/plot_binary_data.py
@@ -8,23 +8,16 @@
 ##
 #################################################

-## Restrict OpenMP to one thread per process so that running
-## many workers in parallel does not create an O(workers * BLAS_threads)
-## thread explosion.  The variable MUST be set before numpy/scipy
-## are imported, because the BLAS library reads them only at load time.
-import os
-os.environ.setdefault("OMP_NUM_THREADS",        "1")
-
 import numpy
 import scipy
-import matplotlib
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt
 from   matplotlib.colors    import LogNorm
 from   mpl_toolkits.mplot3d import Axes3D
 ## import torch
 import AMSS_NCKU_Input      as input_data

+import os
+

 #########################################################################################

@@ -199,19 +192,3 @@ def get_data_xy( Rmin, Rmax, n, data0, time, figure_title, figure_outdir ):

 ####################################################################################

-
-####################################################################################
-## Allow this module to be run as a standalone script so that each
-## binary-data plot can be executed in a fresh subprocess whose BLAS
-## environment variables (set above) take effect before numpy loads.
-##
-## Usage:  python3 plot_binary_data.py <filename> <binary_outdir> <figure_outdir>
-####################################################################################
-
-if __name__ == '__main__':
-    import sys
-    if len(sys.argv) != 4:
-        print(f"Usage: {sys.argv[0]} <filename> <binary_outdir> <figure_outdir>")
-        sys.exit(1)
-    plot_binary_data(sys.argv[1], sys.argv[2], sys.argv[3])
-
--- a/plot_xiaoqu.py
+++ b/plot_xiaoqu.py
@@ -8,8 +8,6 @@
 #################################################

 import numpy                               ## numpy for array operations
-import matplotlib
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt     ## matplotlib for plotting
 from   mpl_toolkits.mplot3d import Axes3D  ## needed for 3D plots
 import glob
@@ -17,9 +15,6 @@ import os                                  ## operating system utilities

 import plot_binary_data
 import AMSS_NCKU_Input as input_data
-import subprocess
-import sys
-import multiprocessing

 # plt.rcParams['text.usetex'] = True  ## enable LaTeX fonts in plots

@@ -55,40 +50,10 @@ def generate_binary_data_plot( binary_outdir, figure_outdir ):
        file_list.append(x)
        print(x)

-    ## Plot each file in parallel using subprocesses.
-    ## Each subprocess is a fresh Python process where the BLAS thread-count
-    ## environment variables (set at the top of plot_binary_data.py) take
-    ## effect before numpy is imported.  This avoids the thread explosion
-    ## that occurs when multiprocessing.Pool with 'fork' context inherits
-    ## already-initialized multi-threaded BLAS from the parent.
-    script = os.path.join( os.path.dirname(__file__), "plot_binary_data.py" )
-    max_workers = min( multiprocessing.cpu_count(), len(file_list) ) if file_list else 0
-
-    running = []
-    failed  = []
+    ## Plot each file in the list
    for filename in file_list:
        print(filename)
-        proc = subprocess.Popen(
-            [sys.executable, script, filename, binary_outdir, figure_outdir],
-        )
-        running.append( (proc, filename) )
-        ## Keep at most max_workers subprocesses active at a time
-        if len(running) >= max_workers:
-            p, fn = running.pop(0)
-            p.wait()
-            if p.returncode != 0:
-                failed.append(fn)
-
-    ## Wait for all remaining subprocesses to finish
-    for p, fn in running:
-        p.wait()
-        if p.returncode != 0:
-            failed.append(fn)
-
-    if failed:
-        print( " WARNING: the following binary data plots failed:" )
-        for fn in failed:
-            print( "   ", fn )
+        plot_binary_data.plot_binary_data(filename, binary_outdir, figure_outdir)

    print(                        )
    print( " Binary Data Plot Has been Finished " )