Refresh synced GPU regions incrementally

Reduce staged GPU host-device copies
Cache repeated interpolation plans
2026-04-09 17:07:31 +08:00 · 2026-04-09 16:44:08 +08:00 · 2026-04-09 15:21:01 +08:00 · 2026-04-09 15:06:11 +08:00 · 2026-04-09 14:56:01 +08:00 · 2026-04-09 14:28:36 +08:00
18 changed files with 4664 additions and 1397 deletions
--- a/AMSS_NCKU_source/ABE.C
+++ b/AMSS_NCKU_source/ABE.C
@@ -29,16 +29,11 @@ using namespace std;
 #error "not define ABEtype"
 #endif

-#if (ABEtype == 0)
-
-#ifdef USE_GPU
-#include "bssn_gpu_class.h"
-#else
-#include "bssn_class.h"
-#endif
-
-#elif (ABEtype == 1)
-#include "bssnEScalar_class.h"
+#if (ABEtype == 0)
+#include "bssn_class.h"
+
+#elif (ABEtype == 1)
+#include "bssnEScalar_class.h"

 #elif (ABEtype == 2)
 #include "Z4c_class.h"
--- a/AMSS_NCKU_source/MPatch.C
+++ b/AMSS_NCKU_source/MPatch.C
@@ -2,29 +2,88 @@
 #include <iostream>
 #include <iomanip>
 #include <fstream>
-#include <cstdlib>
-#include <cstdio>
-#include <string>
-#include <cmath>
-#include <new>
-#include <vector>
-using namespace std;
+#include <cstdlib>
+#include <cstdio>
+#include <string>
+#include <cmath>
+#include <new>
+#include <map>
+#include <vector>
+using namespace std;

-#include "misc.h"
-#include "MPatch.h"
-#include "Parallel.h"
-#include "fmisc.h"
-#ifdef INTERP_LB_PROFILE
-#include "interp_lb_profile.h"
-#endif
-
-namespace
-{
-struct InterpBlockView
-{
-  Block *bp;
-  double llb[dim];
-  double uub[dim];
+#include "misc.h"
+#include "MPatch.h"
+#include "Parallel.h"
+#include "fmisc.h"
+#include "bssn_cuda_ops.h"
+#ifdef INTERP_LB_PROFILE
+#include "interp_lb_profile.h"
+#endif
+
+#if defined(__GNUC__) || defined(__clang__)
+extern int bssn_cuda_interp_points_batch(const int *ex,
+                                         const double *X, const double *Y, const double *Z,
+                                         const double *const *fields,
+                                         const double *soa_flat,
+                                         int num_var,
+                                         const double *px, const double *py, const double *pz,
+                                         int num_points,
+                                         int ordn,
+                                         int symmetry,
+                                         double *out) __attribute__((weak));
+#endif
+
+namespace
+{
+struct InterpVarDesc
+{
+  int sgfn;
+  double soa[dim];
+};
+
+struct InterpPlanKey
+{
+  const Patch *patch;
+  const double *x;
+  const double *y;
+  const double *z;
+  int NN;
+  int Symmetry;
+  int myrank;
+};
+
+struct InterpPlanKeyLess
+{
+  bool operator()(const InterpPlanKey &lhs, const InterpPlanKey &rhs) const
+  {
+    if (lhs.patch != rhs.patch) return lhs.patch < rhs.patch;
+    if (lhs.x != rhs.x) return lhs.x < rhs.x;
+    if (lhs.y != rhs.y) return lhs.y < rhs.y;
+    if (lhs.z != rhs.z) return lhs.z < rhs.z;
+    if (lhs.NN != rhs.NN) return lhs.NN < rhs.NN;
+    if (lhs.Symmetry != rhs.Symmetry) return lhs.Symmetry < rhs.Symmetry;
+    return lhs.myrank < rhs.myrank;
+  }
+};
+
+struct CachedInterpPlan
+{
+  int nblocks;
+  vector<int> owner_rank;
+  vector<int> owner_block;
+  vector<vector<int> > block_points;
+  vector<vector<double> > block_px;
+  vector<vector<double> > block_py;
+  vector<vector<double> > block_pz;
+
+  CachedInterpPlan() : nblocks(0) {}
+};
+
+struct InterpBlockView
+{
+  Block *bp;
+  double llb[dim];
+  double uub[dim];
 };

 struct BlockBinIndex
@@ -154,10 +213,10 @@ void build_block_bin_index(Patch *patch, const double *DH, BlockBinIndex &index)
  index.valid = true;
 }

-int find_block_index_for_point(const BlockBinIndex &index, const double *pox, const double *DH)
-{
-  if (!index.valid)
-    return -1;
+int find_block_index_for_point(const BlockBinIndex &index, const double *pox, const double *DH)
+{
+  if (!index.valid)
+    return -1;

  const int bx = coord_to_bin(pox[0], index.lo[0], index.inv[0], index.bins[0]);
  const int by = coord_to_bin(pox[1], index.lo[1], index.inv[1], index.bins[1]);
@@ -175,10 +234,259 @@ int find_block_index_for_point(const BlockBinIndex &index, const double *pox, co
  for (size_t bi = 0; bi < index.views.size(); bi++)
    if (point_in_block_view(index.views[bi], pox, DH))
      return int(bi);
-
-  return -1;
-}
-} // namespace
+
+  return -1;
+}
+
+void collect_interp_vars(MyList<var> *VarList, vector<InterpVarDesc> &vars)
+{
+  vars.clear();
+  MyList<var> *varl = VarList;
+  while (varl)
+  {
+    InterpVarDesc desc;
+    desc.sgfn = varl->data->sgfn;
+    for (int d = 0; d < dim; ++d)
+      desc.soa[d] = varl->data->SoA[d];
+    vars.push_back(desc);
+    varl = varl->next;
+  }
+}
+
+bool should_try_cuda_interp(int ordn, int num_points, int num_var)
+{
+#if defined(__GNUC__) || defined(__clang__)
+  if (!bssn_cuda_interp_points_batch)
+    return false;
+#else
+  return false;
+#endif
+  if (ordn != 6)
+    return false;
+  if (num_points < 32)
+    return false;
+  return num_points * num_var >= 256;
+}
+
+CachedInterpPlan &get_cached_interp_plan(Patch *patch,
+                                         int NN, double **XX,
+                                         int Symmetry, int myrank,
+                                         const double *DH,
+                                         const BlockBinIndex &block_index,
+                                         bool report_bounds_here,
+                                         bool allow_missing_points)
+{
+  static map<InterpPlanKey, CachedInterpPlan, InterpPlanKeyLess> cache;
+
+  InterpPlanKey key;
+  key.patch = patch;
+  key.x = XX[0];
+  key.y = XX[1];
+  key.z = XX[2];
+  key.NN = NN;
+  key.Symmetry = Symmetry;
+  key.myrank = myrank;
+
+  map<InterpPlanKey, CachedInterpPlan, InterpPlanKeyLess>::iterator it = cache.find(key);
+  if (it != cache.end() && it->second.nblocks == static_cast<int>(block_index.views.size()))
+    return it->second;
+
+  CachedInterpPlan &plan = cache[key];
+  plan = CachedInterpPlan();
+  plan.nblocks = static_cast<int>(block_index.views.size());
+  plan.owner_rank.assign(NN, -1);
+  plan.owner_block.assign(NN, -1);
+  plan.block_points.resize(plan.nblocks);
+  plan.block_px.resize(plan.nblocks);
+  plan.block_py.resize(plan.nblocks);
+  plan.block_pz.resize(plan.nblocks);
+
+  for (int j = 0; j < NN; ++j)
+  {
+    double pox[dim];
+    for (int i = 0; i < dim; ++i)
+    {
+      pox[i] = XX[i][j];
+      if (report_bounds_here &&
+          (XX[i][j] < patch->bbox[i] + patch->lli[i] * DH[i] ||
+           XX[i][j] > patch->bbox[dim + i] - patch->uui[i] * DH[i]))
+      {
+        cout << "Patch::Interp_Points: point (";
+        for (int k = 0; k < dim; ++k)
+        {
+          cout << XX[k][j];
+          if (k < dim - 1)
+            cout << ",";
+          else
+            cout << ") is out of current Patch." << endl;
+        }
+        MPI_Abort(MPI_COMM_WORLD, 1);
+      }
+    }
+
+    const int block_i = find_block_index_for_point(block_index, pox, DH);
+    if (block_i >= 0)
+    {
+      Block *BP = block_index.views[block_i].bp;
+      plan.owner_rank[j] = BP->rank;
+      plan.owner_block[j] = block_i;
+      if (BP->rank == myrank)
+      {
+        plan.block_points[block_i].push_back(j);
+        plan.block_px[block_i].push_back(XX[0][j]);
+        plan.block_py[block_i].push_back(XX[1][j]);
+        plan.block_pz[block_i].push_back(XX[2][j]);
+      }
+    }
+  }
+
+  if (!allow_missing_points && report_bounds_here)
+  {
+    for (int j = 0; j < NN; ++j)
+    {
+      if (plan.owner_rank[j] >= 0)
+        continue;
+      cout << "ERROR: Patch::Interp_Points fails to find point (";
+      for (int d = 0; d < dim; ++d)
+      {
+        cout << XX[d][j];
+        if (d < dim - 1)
+          cout << ",";
+        else
+          cout << ")";
+      }
+      cout << " on Patch (";
+      for (int d = 0; d < dim; ++d)
+      {
+        cout << patch->bbox[d] << "+" << patch->lli[d] * DH[d];
+        if (d < dim - 1)
+          cout << ",";
+        else
+          cout << ")--";
+      }
+      cout << "(";
+      for (int d = 0; d < dim; ++d)
+      {
+        cout << patch->bbox[dim + d] << "-" << patch->uui[d] * DH[d];
+        if (d < dim - 1)
+          cout << ",";
+        else
+          cout << ")" << endl;
+      }
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+  }
+
+  return plan;
+}
+
+bool run_cuda_interp_for_block(Block *BP,
+                               const vector<InterpVarDesc> &vars,
+                               const vector<int> &point_ids,
+                               const vector<double> &px,
+                               const vector<double> &py,
+                               const vector<double> &pz,
+                               double *Shellf,
+                               int num_var,
+                               int ordn,
+                               int Symmetry)
+{
+  if (!should_try_cuda_interp(ordn, static_cast<int>(point_ids.size()), num_var))
+    return false;
+
+  vector<const double *> field_ptrs(num_var);
+  vector<double> soa_flat(3 * num_var);
+  for (int v = 0; v < num_var; ++v)
+  {
+    field_ptrs[v] = BP->fgfs[vars[v].sgfn];
+    for (int d = 0; d < dim; ++d)
+      soa_flat[3 * v + d] = vars[v].soa[d];
+  }
+
+  const int npts = static_cast<int>(point_ids.size());
+  vector<double> out(static_cast<size_t>(npts) * static_cast<size_t>(num_var));
+  if (bssn_cuda_interp_points_batch(BP->shape,
+                                    BP->X[0], BP->X[1], BP->X[2],
+                                    field_ptrs.data(),
+                                    soa_flat.data(),
+                                    num_var,
+                                    px.data(), py.data(), pz.data(),
+                                    npts,
+                                    ordn,
+                                    Symmetry,
+                                    out.data()) != 0)
+  {
+    return false;
+  }
+
+  for (int p = 0; p < npts; ++p)
+  {
+    const int j = point_ids[p];
+    memcpy(Shellf + j * num_var, out.data() + p * num_var, sizeof(double) * num_var);
+  }
+  return true;
+}
+
+void run_cpu_interp_for_block(Block *BP,
+                              const vector<InterpVarDesc> &vars,
+                              const vector<int> &point_ids,
+                              const vector<double> &px,
+                              const vector<double> &py,
+                              const vector<double> &pz,
+                              double *Shellf,
+                              int num_var,
+                              int ordn,
+                              int Symmetry)
+{
+  for (size_t p = 0; p < point_ids.size(); ++p)
+  {
+    const int j = point_ids[p];
+    double x = px[p];
+    double y = py[p];
+    double z = pz[p];
+    int ordn_local = ordn;
+    int symmetry_local = Symmetry;
+    for (int v = 0; v < num_var; ++v)
+    {
+      f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2],
+                      BP->fgfs[vars[v].sgfn], Shellf[j * num_var + v],
+                      x, y, z, ordn_local, const_cast<double *>(vars[v].soa), symmetry_local);
+    }
+  }
+}
+
+void interpolate_owned_points(MyList<var> *VarList,
+                              double *Shellf, int Symmetry,
+                              int ordn,
+                              const BlockBinIndex &block_index,
+                              const CachedInterpPlan &plan)
+{
+  vector<InterpVarDesc> vars;
+  collect_interp_vars(VarList, vars);
+  const int num_var = static_cast<int>(vars.size());
+
+  for (size_t bi = 0; bi < plan.block_points.size(); ++bi)
+  {
+    if (plan.block_points[bi].empty())
+      continue;
+
+    Block *BP = block_index.views[bi].bp;
+    bool done = run_cuda_interp_for_block(BP, vars,
+                                          plan.block_points[bi],
+                                          plan.block_px[bi],
+                                          plan.block_py[bi],
+                                          plan.block_pz[bi],
+                                          Shellf, num_var, ordn, Symmetry);
+    if (!done)
+      run_cpu_interp_for_block(BP, vars,
+                               plan.block_points[bi],
+                               plan.block_px[bi],
+                               plan.block_py[bi],
+                               plan.block_pz[bi],
+                               Shellf, num_var, ordn, Symmetry);
+  }
+}
+} // namespace

 Patch::Patch(int DIM, int *shapei, double *bboxi, int levi, bool buflog, int Symmetry) : lev(levi)
 {
@@ -523,60 +831,15 @@ void Patch::Interp_Points(MyList<var> *VarList,

  memset(Shellf, 0, sizeof(double) * NN * num_var);

-  // owner_rank[j] records which MPI rank owns point j
-  // All ranks traverse the same block list so they all agree on ownership
-  int *owner_rank;
-  owner_rank = new int[NN];
-  for (int j = 0; j < NN; j++)
-    owner_rank[j] = -1;
-
-  double DH[dim];
-  for (int i = 0; i < dim; i++)
-    DH[i] = getdX(i);
-  BlockBinIndex block_index;
-  build_block_bin_index(this, DH, block_index);
-
-  for (int j = 0; j < NN; j++) // run along points
-  {
-    double pox[dim];
-    for (int i = 0; i < dim; i++)
-    {
-      pox[i] = XX[i][j];
-      if (myrank == 0 && (XX[i][j] < bbox[i] + lli[i] * DH[i] || XX[i][j] > bbox[dim + i] - uui[i] * DH[i]))
-      {
-        cout << "Patch::Interp_Points: point (";
-        for (int k = 0; k < dim; k++)
-        {
-          cout << XX[k][j];
-          if (k < dim - 1)
-            cout << ",";
-          else
-            cout << ") is out of current Patch." << endl;
-        }
-        MPI_Abort(MPI_COMM_WORLD, 1);
-      }
-    }
-
-    const int block_i = find_block_index_for_point(block_index, pox, DH);
-    if (block_i >= 0)
-    {
-      Block *BP = block_index.views[block_i].bp;
-      owner_rank[j] = BP->rank;
-      if (myrank == BP->rank)
-      {
-        //---> interpolation
-        varl = VarList;
-        int k = 0;
-        while (varl) // run along variables
-        {
-          f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
-                          pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
-          varl = varl->next;
-          k++;
-        }
-      }
-    }
-  }
+  double DH[dim];
+  for (int i = 0; i < dim; i++)
+    DH[i] = getdX(i);
+  BlockBinIndex block_index;
+  build_block_bin_index(this, DH, block_index);
+  CachedInterpPlan &plan = get_cached_interp_plan(this, NN, XX, Symmetry, myrank, DH, block_index, myrank == 0, false);
+  const int *owner_rank = plan.owner_rank.data();
+
+  interpolate_owned_points(VarList, Shellf, Symmetry, ordn, block_index, plan);

  // Replace MPI_Allreduce with per-owner MPI_Bcast:
  // Group consecutive points by owner rank and broadcast each group.
@@ -631,9 +894,8 @@ void Patch::Interp_Points(MyList<var> *VarList,
      MPI_Bcast(Shellf + jstart * num_var, count, MPI_DOUBLE, cur_owner, MPI_COMM_WORLD);
    }
  }
-
-  delete[] owner_rank;
-}
+
+}
 void Patch::Interp_Points(MyList<var> *VarList,
                          int NN, double **XX,
                          double *Shellf, int Symmetry,
@@ -661,102 +923,22 @@ void Patch::Interp_Points(MyList<var> *VarList,

  memset(Shellf, 0, sizeof(double) * NN * num_var);

-  // owner_rank[j] records which MPI rank owns point j
-  int *owner_rank;
-  owner_rank = new int[NN];
-  for (int j = 0; j < NN; j++)
-    owner_rank[j] = -1;
-
-  double DH[dim];
-  for (int i = 0; i < dim; i++)
-    DH[i] = getdX(i);
-  BlockBinIndex block_index;
-  build_block_bin_index(this, DH, block_index);
-
-  // --- Interpolation phase (identical to original) ---
-  for (int j = 0; j < NN; j++)
-  {
-    double pox[dim];
-    for (int i = 0; i < dim; i++)
-    {
-      pox[i] = XX[i][j];
-      if (myrank == 0 && (XX[i][j] < bbox[i] + lli[i] * DH[i] || XX[i][j] > bbox[dim + i] - uui[i] * DH[i]))
-      {
-        cout << "Patch::Interp_Points: point (";
-        for (int k = 0; k < dim; k++)
-        {
-          cout << XX[k][j];
-          if (k < dim - 1)
-            cout << ",";
-          else
-            cout << ") is out of current Patch." << endl;
-        }
-        MPI_Abort(MPI_COMM_WORLD, 1);
-      }
-    }
-
-    const int block_i = find_block_index_for_point(block_index, pox, DH);
-    if (block_i >= 0)
-    {
-      Block *BP = block_index.views[block_i].bp;
-      owner_rank[j] = BP->rank;
-      if (myrank == BP->rank)
-      {
-        varl = VarList;
-        int k = 0;
-        while (varl)
-        {
-          f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
-                          pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
-          varl = varl->next;
-          k++;
-        }
-      }
-    }
-  }
+  double DH[dim];
+  for (int i = 0; i < dim; i++)
+    DH[i] = getdX(i);
+  BlockBinIndex block_index;
+  build_block_bin_index(this, DH, block_index);
+  CachedInterpPlan &plan = get_cached_interp_plan(this, NN, XX, Symmetry, myrank, DH, block_index, myrank == 0, false);
+  const int *owner_rank = plan.owner_rank.data();
+
+  interpolate_owned_points(VarList, Shellf, Symmetry, ordn, block_index, plan);

 #ifdef INTERP_LB_PROFILE
  double t_interp_end = MPI_Wtime();
  double t_interp_local = t_interp_end - t_interp_start;
 #endif

-  // --- Error check for unfound points ---
-  for (int j = 0; j < NN; j++)
-  {
-    if (owner_rank[j] < 0 && myrank == 0)
-    {
-      cout << "ERROR: Patch::Interp_Points fails to find point (";
-      for (int d = 0; d < dim; d++)
-      {
-        cout << XX[d][j];
-        if (d < dim - 1)
-          cout << ",";
-        else
-          cout << ")";
-      }
-      cout << " on Patch (";
-      for (int d = 0; d < dim; d++)
-      {
-        cout << bbox[d] << "+" << lli[d] * DH[d];
-        if (d < dim - 1)
-          cout << ",";
-        else
-          cout << ")--";
-      }
-      cout << "(";
-      for (int d = 0; d < dim; d++)
-      {
-        cout << bbox[dim + d] << "-" << uui[d] * DH[d];
-        if (d < dim - 1)
-          cout << ",";
-        else
-          cout << ")" << endl;
-      }
-      MPI_Abort(MPI_COMM_WORLD, 1);
-    }
-  }
-
-  // --- Targeted point-to-point communication phase ---
+  // --- Targeted point-to-point communication phase ---
  // Compute consumer_rank[j] using the same deterministic formula as surface_integral
  int *consumer_rank = new int[NN];
  {
@@ -873,9 +1055,8 @@ void Patch::Interp_Points(MyList<var> *VarList,
  delete[] send_offset;
  delete[] recv_offset;
  delete[] send_count;
-  delete[] recv_count;
-  delete[] consumer_rank;
-  delete[] owner_rank;
+  delete[] recv_count;
+  delete[] consumer_rank;

 #ifdef INTERP_LB_PROFILE
  {
@@ -923,64 +1104,20 @@ void Patch::Interp_Points(MyList<var> *VarList,

  memset(Shellf, 0, sizeof(double) * NN * num_var);

-  // owner_rank[j] stores the global rank that owns point j
-  int *owner_rank;
-  owner_rank = new int[NN];
-  for (int j = 0; j < NN; j++)
-    owner_rank[j] = -1;
+  // Build global-to-local rank translation for Comm_here
+  MPI_Group world_group, local_group;
+  MPI_Comm_group(MPI_COMM_WORLD, &world_group);
+  MPI_Comm_group(Comm_here, &local_group);

-  // Build global-to-local rank translation for Comm_here
-  MPI_Group world_group, local_group;
-  MPI_Comm_group(MPI_COMM_WORLD, &world_group);
-  MPI_Comm_group(Comm_here, &local_group);
-
-  double DH[dim];
-  for (int i = 0; i < dim; i++)
-    DH[i] = getdX(i);
-  BlockBinIndex block_index;
-  build_block_bin_index(this, DH, block_index);
-
-  for (int j = 0; j < NN; j++) // run along points
-  {
-    double pox[dim];
-    for (int i = 0; i < dim; i++)
-    {
-      pox[i] = XX[i][j];
-      if (lmyrank == 0 && (XX[i][j] < bbox[i] + lli[i] * DH[i] || XX[i][j] > bbox[dim + i] - uui[i] * DH[i]))
-      {
-        cout << "Patch::Interp_Points: point (";
-        for (int k = 0; k < dim; k++)
-        {
-          cout << XX[k][j];
-          if (k < dim - 1)
-            cout << ",";
-          else
-            cout << ") is out of current Patch." << endl;
-        }
-        MPI_Abort(MPI_COMM_WORLD, 1);
-      }
-    }
-
-    const int block_i = find_block_index_for_point(block_index, pox, DH);
-    if (block_i >= 0)
-    {
-      Block *BP = block_index.views[block_i].bp;
-      owner_rank[j] = BP->rank;
-      if (myrank == BP->rank)
-      {
-        //---> interpolation
-        varl = VarList;
-        int k = 0;
-        while (varl) // run along variables
-        {
-          f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
-                          pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
-          varl = varl->next;
-          k++;
-        }
-      }
-    }
-  }
+  double DH[dim];
+  for (int i = 0; i < dim; i++)
+    DH[i] = getdX(i);
+  BlockBinIndex block_index;
+  build_block_bin_index(this, DH, block_index);
+  CachedInterpPlan &plan = get_cached_interp_plan(this, NN, XX, Symmetry, myrank, DH, block_index, lmyrank == 0, true);
+  const int *owner_rank = plan.owner_rank.data();
+
+  interpolate_owned_points(VarList, Shellf, Symmetry, ordn, block_index, plan);

  // Collect unique global owner ranks and translate to local ranks in Comm_here
  // Then broadcast each owner's points via MPI_Bcast on Comm_here
@@ -1008,10 +1145,9 @@ void Patch::Interp_Points(MyList<var> *VarList,
    }
  }

-  MPI_Group_free(&world_group);
-  MPI_Group_free(&local_group);
-  delete[] owner_rank;
-}
+  MPI_Group_free(&world_group);
+  MPI_Group_free(&local_group);
+}
 void Patch::checkBlock()
 {
  int myrank;
--- a/AMSS_NCKU_source/Parallel.C
+++ b/AMSS_NCKU_source/Parallel.C
@@ -1,9 +1,124 @@

-#include "Parallel.h"
-#include "fmisc.h"
-#include "prolongrestrict.h"
-#include "misc.h"
-#include "parameters.h"
+#include "Parallel.h"
+#include "fmisc.h"
+#include "prolongrestrict.h"
+#include "bssn_cuda_ops.h"
+#include "misc.h"
+#include "parameters.h"
+#include <cstring>
+
+#if defined(__GNUC__) || defined(__clang__)
+extern int bssn_cuda_prolong3_pack(int wei,
+                                   const double *llbc, const double *uubc, const int *extc, const double *func,
+                                   const double *llbf, const double *uubf, const int *extf, double *funf,
+                                   const double *llbp, const double *uubp,
+                                   const double *SoA, int symmetry) __attribute__((weak));
+#endif
+
+namespace {
+
+const char *g_parallel_transfer_context = "Parallel::transfer";
+int parallel_next_transfer_tag()
+{
+  const char *ctx = g_parallel_transfer_context ? g_parallel_transfer_context : "Parallel::transfer";
+  unsigned int hash = 5381u;
+  while (*ctx)
+  {
+    hash = ((hash << 5) + hash) + static_cast<unsigned char>(*ctx);
+    ctx++;
+  }
+  return 32 + static_cast<int>(hash % 32000u);
+}
+
+struct ParallelTransferContextGuard
+{
+  const char *prev;
+  explicit ParallelTransferContextGuard(const char *context) : prev(g_parallel_transfer_context)
+  {
+    g_parallel_transfer_context = context;
+  }
+  ~ParallelTransferContextGuard()
+  {
+    g_parallel_transfer_context = prev;
+  }
+};
+
+void parallel_report_mpi_error(const char *context, int errcode, int req_no)
+{
+  char errstr[MPI_MAX_ERROR_STRING];
+  int len = 0;
+  MPI_Error_string(errcode, errstr, &len);
+  int myrank = -1;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  fprintf(stderr, "MPI failure on rank %d in %s (transfer_ctx=%s, req_no=%d): %.*s\n",
+          myrank, context, g_parallel_transfer_context, req_no, len, errstr);
+  fflush(stderr);
+}
+
+void parallel_report_mpi_status_errors(const char *context, int req_no,
+                                       int status_count, MPI_Status *stats)
+{
+  if (!stats)
+    return;
+
+  int myrank = -1;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+  const int count = (status_count >= 0 && status_count <= req_no) ? status_count : req_no;
+  for (int i = 0; i < count; i++)
+  {
+    if (stats[i].MPI_ERROR != MPI_SUCCESS)
+    {
+      char errstr[MPI_MAX_ERROR_STRING];
+      int len = 0;
+      MPI_Error_string(stats[i].MPI_ERROR, errstr, &len);
+      fprintf(stderr,
+              "MPI request failure on rank %d in %s (transfer_ctx=%s, status_idx=%d, source=%d, tag=%d): %.*s\n",
+              myrank, context, g_parallel_transfer_context, i,
+              stats[i].MPI_SOURCE, stats[i].MPI_TAG, len, errstr);
+    }
+  }
+  fflush(stderr);
+}
+
+int parallel_waitsome_checked(int req_no, MPI_Request *reqs, int *outcount,
+                              int *completed, MPI_Status *stats,
+                              const char *context)
+{
+  MPI_Errhandler old_handler;
+  MPI_Comm_get_errhandler(MPI_COMM_WORLD, &old_handler);
+  MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
+  int rc = MPI_Waitsome(req_no, reqs, outcount, completed, stats);
+  MPI_Comm_set_errhandler(MPI_COMM_WORLD, old_handler);
+  MPI_Errhandler_free(&old_handler);
+  if (rc != MPI_SUCCESS)
+  {
+    parallel_report_mpi_status_errors(context, req_no, outcount ? *outcount : req_no, stats);
+    parallel_report_mpi_error(context, rc, req_no);
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+  return rc;
+}
+
+int parallel_waitall_checked(int req_no, MPI_Request *reqs, MPI_Status *stats,
+                             const char *context)
+{
+  MPI_Errhandler old_handler;
+  MPI_Comm_get_errhandler(MPI_COMM_WORLD, &old_handler);
+  MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
+  int rc = MPI_Waitall(req_no, reqs, stats);
+  MPI_Comm_set_errhandler(MPI_COMM_WORLD, old_handler);
+  MPI_Errhandler_free(&old_handler);
+  if (rc != MPI_SUCCESS)
+  {
+    parallel_report_mpi_status_errors(context, req_no, req_no, stats);
+    parallel_report_mpi_error(context, rc, req_no);
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+  return rc;
+}
+
+} // namespace

 int Parallel::partition1(int &nx, int split_size, int min_width, int cpusize, int shape) // special for 1 diemnsion
 {
@@ -3779,11 +3894,18 @@ int Parallel::data_packer(double *data, MyList<Parallel::gridseg> *src, MyList<P
                          src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
                          dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry);
              break;
-            case 3:
-              f_prolong3(DIM, src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
-                         dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
-                         dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry);
-            }
+            case 3:
+              if (!bssn_cuda_prolong3_pack ||
+                  bssn_cuda_prolong3_pack(DIM,
+                                          src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
+                                          dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
+                                          dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry))
+              {
+                f_prolong3(DIM, src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
+                           dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
+                           dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry);
+              }
+            }
          if (dir == UNPACK) // from target data to corresponding grid
            f_copy(DIM, dst->data->Bg->bbox, dst->data->Bg->bbox + dim, dst->data->Bg->shape, dst->data->Bg->fgfs[varld->data->sgfn],
                   dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
@@ -3896,10 +4018,11 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
  MPI_Request *reqs = new MPI_Request[2 * cpusize];
  MPI_Status *stats = new MPI_Status[2 * cpusize];
  int *req_node = new int[2 * cpusize];
-  int *req_is_recv = new int[2 * cpusize];
-  int *completed = new int[2 * cpusize];
-  int req_no = 0;
-  int pending_recv = 0;
+  int *req_is_recv = new int[2 * cpusize];
+  int *completed = new int[2 * cpusize];
+  int req_no = 0;
+  int pending_recv = 0;
+  const int mpi_tag = parallel_next_transfer_tag();

  double **send_data = new double *[cpusize];
  double **rec_data = new double *[cpusize];
@@ -3913,31 +4036,31 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
  }

  // Post receives first so peers can progress rendezvous early.
-  for (node = 0; node < cpusize; node++)
-  {
-    if (node == myrank) continue;
-
-    recv_lengths[node] = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
-    if (recv_lengths[node] > 0)
-    {
+  for (node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    recv_lengths[node] = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+    if (recv_lengths[node] > 0)
+    {
      rec_data[node] = new double[recv_lengths[node]];
      if (!rec_data[node])
      {
        cout << "out of memory when new in short transfer, place 1" << endl;
        MPI_Abort(MPI_COMM_WORLD, 1);
      }
-      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 1;
      req_no++;
      pending_recv++;
    }
-  }
-
-  // Local transfer on this rank.
-  recv_lengths[myrank] = data_packer(0, src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
-  if (recv_lengths[myrank] > 0)
-  {
+  }
+
+  // Local transfer on this rank.
+  recv_lengths[myrank] = data_packer(0, src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  if (recv_lengths[myrank] > 0)
+  {
    rec_data[myrank] = new double[recv_lengths[myrank]];
    if (!rec_data[myrank])
    {
@@ -3947,33 +4070,33 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
    data_packer(rec_data[myrank], src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
  }

-  // Pack and post sends.
-  for (node = 0; node < cpusize; node++)
-  {
-    if (node == myrank) continue;
-
-    send_lengths[node] = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-    if (send_lengths[node] > 0)
-    {
-      send_data[node] = new double[send_lengths[node]];
+  // Pack and post sends.
+  for (node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    send_lengths[node] = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+    if (send_lengths[node] > 0)
+    {
+      send_data[node] = new double[send_lengths[node]];
      if (!send_data[node])
      {
        cout << "out of memory when new in short transfer, place 3" << endl;
        MPI_Abort(MPI_COMM_WORLD, 1);
      }
      data_packer(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 0;
-      req_no++;
-    }
-  }
-
-  // Unpack as soon as receive completes to reduce pure wait time.
-  while (pending_recv > 0)
+      req_no++;
+    }
+  }
+  // Unpack as soon as receive completes to reduce pure wait time.
+  while (pending_recv > 0)
  {
    int outcount = 0;
-    MPI_Waitsome(req_no, reqs, &outcount, completed, stats);
+    parallel_waitsome_checked(req_no, reqs, &outcount, completed, stats,
+                              "Parallel::transfer");
    if (outcount == MPI_UNDEFINED) break;

    for (int i = 0; i < outcount; i++)
@@ -3988,7 +4111,7 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
    }
  }

-  if (req_no > 0) MPI_Waitall(req_no, reqs, stats);
+  if (req_no > 0) parallel_waitall_checked(req_no, reqs, stats, "Parallel::transfer");

  if (rec_data[myrank])
    data_packer(rec_data[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
@@ -4005,12 +4128,12 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
  delete[] stats;
  delete[] req_node;
  delete[] req_is_recv;
-  delete[] completed;
-  delete[] send_data;
-  delete[] rec_data;
-  delete[] send_lengths;
-  delete[] recv_lengths;
-}
+  delete[] completed;
+  delete[] send_data;
+  delete[] rec_data;
+  delete[] send_lengths;
+  delete[] recv_lengths;
+}
 //
 void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gridseg> **dst,
                           MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /*target */,
@@ -4027,8 +4150,9 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
  int *req_node = new int[2 * cpusize];
  int *req_is_recv = new int[2 * cpusize];
  int *completed = new int[2 * cpusize];
-  int req_no = 0;
-  int pending_recv = 0;
+  int req_no = 0;
+  int pending_recv = 0;
+  const int mpi_tag = parallel_next_transfer_tag();

  double **send_data = new double *[cpusize];
  double **rec_data = new double *[cpusize];
@@ -4055,7 +4179,7 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
        cout << "out of memory when new in short transfer, place 1" << endl;
        MPI_Abort(MPI_COMM_WORLD, 1);
      }
-      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 1;
      req_no++;
@@ -4091,7 +4215,7 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
        MPI_Abort(MPI_COMM_WORLD, 1);
      }
      data_packermix(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 0;
      req_no++;
@@ -4102,7 +4226,8 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
  while (pending_recv > 0)
  {
    int outcount = 0;
-    MPI_Waitsome(req_no, reqs, &outcount, completed, stats);
+    parallel_waitsome_checked(req_no, reqs, &outcount, completed, stats,
+                              "Parallel::transfermix");
    if (outcount == MPI_UNDEFINED) break;

    for (int i = 0; i < outcount; i++)
@@ -4117,7 +4242,7 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
    }
  }

-  if (req_no > 0) MPI_Waitall(req_no, reqs, stats);
+  if (req_no > 0) parallel_waitall_checked(req_no, reqs, stats, "Parallel::transfermix");

  if (rec_data[myrank])
    data_packermix(rec_data[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
@@ -4134,16 +4259,21 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
  delete[] stats;
  delete[] req_node;
  delete[] req_is_recv;
-  delete[] completed;
-  delete[] send_data;
-  delete[] rec_data;
-  delete[] send_lengths;
-  delete[] recv_lengths;
-}
-void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry)
-{
-  int cpusize;
-  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
+  delete[] completed;
+  delete[] send_data;
+  delete[] rec_data;
+  delete[] send_lengths;
+  delete[] recv_lengths;
+}
+void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry)
+{
+  Sync(Pat, VarList, Symmetry, "Parallel::Sync(Patch)");
+}
+
+void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry, const char *context)
+{
+  int cpusize;
+  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);

  MyList<Parallel::gridseg> *dst;
  MyList<Parallel::gridseg> **src, **transfer_src, **transfer_dst;
@@ -4159,7 +4289,10 @@ void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry)
                                                                          // but for transfer_dst[node] the data may locate on any node
  }

-  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx(context ? context : "Parallel::Sync(Patch)");
+    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -4177,15 +4310,26 @@ void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry)
  delete[] transfer_src;
  delete[] transfer_dst;
 }
-void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
-{
-  // Patch inner Synch
-  MyList<Patch> *Pp = PatL;
-  while (Pp)
-  {
-    Sync(Pp->data, VarList, Symmetry);
-    Pp = Pp->next;
-  }
+void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
+{
+  Sync(PatL, VarList, Symmetry, "Parallel::Sync(PatchList)");
+}
+
+void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, const char *context)
+{
+  std::string patchlist_context = context ? std::string(context) + " [patchlist]" : "Parallel::Sync(PatchList)";
+
+  // Patch inner Synch
+  MyList<Patch> *Pp = PatL;
+  int patch_index = 0;
+  while (Pp)
+  {
+    std::string patch_context = context ? std::string(context) + " [patch " + std::to_string(patch_index) + "]"
+                                        : "Parallel::Sync(Patch)";
+    Sync(Pp->data, VarList, Symmetry, patch_context.c_str());
+    Pp = Pp->next;
+    patch_index++;
+  }

  // Patch inter Synch
  int cpusize;
@@ -4204,7 +4348,10 @@ void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
    build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx(patchlist_context.c_str());
+    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -4224,10 +4371,15 @@ void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
 }
 // Merged Sync: collect all intra-patch and inter-patch grid segment lists,
 // then issue a single transfer() call instead of N+1 separate ones.
-void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
-{
-  int cpusize;
-  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
+void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
+{
+  Sync_merged(PatL, VarList, Symmetry, "Parallel::Sync_merged");
+}
+
+void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, const char *context)
+{
+  int cpusize;
+  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);

  MyList<Parallel::gridseg> **combined_src = new MyList<Parallel::gridseg> *[cpusize];
  MyList<Parallel::gridseg> **combined_dst = new MyList<Parallel::gridseg> *[cpusize];
@@ -4301,8 +4453,11 @@ void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmet
  if (dst_buffer)
    dst_buffer->destroyList();

-  // Phase C: Single transfer
-  transfer(combined_src, combined_dst, VarList, VarList, Symmetry);
+  // Phase C: Single transfer
+  {
+    ParallelTransferContextGuard transfer_ctx(context ? context : "Parallel::Sync_merged");
+    transfer(combined_src, combined_dst, VarList, VarList, Symmetry);
+  }

  // Phase D: Cleanup
  for (int node = 0; node < cpusize; node++)
@@ -4380,8 +4535,9 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  int cpusize = cache.cpusize;

-  int req_no = 0;
-  int pending_recv = 0;
+  int req_no = 0;
+  int pending_recv = 0;
+  const int mpi_tag = parallel_next_transfer_tag();
  int node;
  int *req_node = cache.tc_req_node;
  int *req_is_recv = cache.tc_req_is_recv;
@@ -4402,7 +4558,7 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
        cache.recv_bufs[node] = new double[rlength];
        cache.recv_buf_caps[node] = rlength;
      }
-      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 1;
      req_no++;
@@ -4440,7 +4596,7 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
        cache.send_buf_caps[node] = slength;
      }
      data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 0;
      req_no++;
@@ -4451,7 +4607,8 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
  while (pending_recv > 0)
  {
    int outcount = 0;
-    MPI_Waitsome(req_no, cache.reqs, &outcount, completed, cache.stats);
+    parallel_waitsome_checked(req_no, cache.reqs, &outcount, completed, cache.stats,
+                              "Parallel::transfer_cached");
    if (outcount == MPI_UNDEFINED) break;

    for (int i = 0; i < outcount; i++)
@@ -4466,11 +4623,11 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
    }
  }

-  if (req_no > 0) MPI_Waitall(req_no, cache.reqs, cache.stats);
+  if (req_no > 0) parallel_waitall_checked(req_no, cache.reqs, cache.stats, "Parallel::transfer_cached");

-  if (self_len > 0)
-    data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
-}
+  if (self_len > 0)
+    data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
+}
 void Parallel::Sync_cached(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, SyncCache &cache)
 {
  if (!cache.valid)
@@ -4669,12 +4826,13 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
  }

  // Now pack and post async MPI operations
-  int myrank;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  int cpusize = cache.cpusize;
-  state.req_no = 0;
-  state.active = true;
-  state.pending_recv = 0;
+  int myrank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  int cpusize = cache.cpusize;
+  state.req_no = 0;
+  state.active = true;
+  state.mpi_tag = parallel_next_transfer_tag();
+  state.pending_recv = 0;
  // Allocate tracking arrays
  delete[] state.req_node; delete[] state.req_is_recv;
  state.req_node = new int[cache.max_reqs];
@@ -4725,7 +4883,7 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
        data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
        state.req_node[state.req_no] = node;
        state.req_is_recv[state.req_no] = 0;
-        MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
+        MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, state.mpi_tag, MPI_COMM_WORLD, cache.reqs + state.req_no++);
      }
      int rlength;
      if (!cache.lengths_valid) {
@@ -4745,7 +4903,7 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
        state.req_node[state.req_no] = node;
        state.req_is_recv[state.req_no] = 1;
        state.pending_recv++;
-        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
+        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, state.mpi_tag, MPI_COMM_WORLD, cache.reqs + state.req_no++);
      }
    }
  }
@@ -4775,7 +4933,8 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
    while (pending > 0)
    {
      int outcount = 0;
-      MPI_Waitsome(state.req_no, cache.reqs, &outcount, completed, cache.stats);
+      parallel_waitsome_checked(state.req_no, cache.reqs, &outcount, completed, cache.stats,
+                                "Parallel::Sync_finish");
      if (outcount == MPI_UNDEFINED) break;
      for (int i = 0; i < outcount; i++)
      {
@@ -4791,8 +4950,8 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
    delete[] completed;
  }

-  // Wait for remaining sends
-  if (state.req_no > 0) MPI_Waitall(state.req_no, cache.reqs, cache.stats);
+  // Wait for remaining sends
+  if (state.req_no > 0) parallel_waitall_checked(state.req_no, cache.reqs, cache.stats, "Parallel::Sync_finish");

  delete[] state.req_node; state.req_node = 0;
  delete[] state.req_is_recv; state.req_is_recv = 0;
@@ -5235,7 +5394,10 @@ void Parallel::PeriodicBD(Patch *Pat, MyList<var> *VarList, int Symmetry)
    build_PhysBD_gstl(Pat, src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx("Parallel::PeriodicBD");
+    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -5506,7 +5668,10 @@ void Parallel::Prolong(Patch *Patc, Patch *Patf,
    build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx("Parallel::Prolong");
+    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -5565,7 +5730,10 @@ void Parallel::Restrict(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
    build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx("Parallel::Restrict");
+    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -5609,7 +5777,10 @@ void Parallel::Restrict_after(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
    build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx("Parallel::Restrict_after");
+    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -5654,7 +5825,10 @@ void Parallel::OutBdLow2Hi(Patch *Patc, Patch *Patf,
    build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx("Parallel::OutBdLow2Hi(Patch)");
+    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -5710,7 +5884,10 @@ void Parallel::OutBdLow2Hi(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
    build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx("Parallel::OutBdLow2Hi(PatchList)");
+    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
+  }

  if (dst)
    dst->destroyList();
@@ -5981,13 +6158,14 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
  int myrank;
  MPI_Comm_size(MPI_COMM_WORLD, &cache.cpusize);
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  int cpusize = cache.cpusize;
-
-  int req_no = 0;
-  int pending_recv = 0;
-  int *req_node = new int[cache.max_reqs];
-  int *req_is_recv = new int[cache.max_reqs];
-  int *completed = new int[cache.max_reqs];
+  int cpusize = cache.cpusize;
+
+  int req_no = 0;
+  int pending_recv = 0;
+  const int mpi_tag = parallel_next_transfer_tag();
+  int *req_node = new int[cache.max_reqs];
+  int *req_is_recv = new int[cache.max_reqs];
+  int *completed = new int[cache.max_reqs];

  // Post receives first so peers can progress rendezvous early.
  for (int node = 0; node < cpusize; node++)
@@ -6004,7 +6182,7 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
        cache.recv_bufs[node] = new double[rlength];
        cache.recv_buf_caps[node] = rlength;
      }
-      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 1;
      req_no++;
@@ -6042,7 +6220,7 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
        cache.send_buf_caps[node] = slength;
      }
      data_packermix(cache.send_bufs[node], cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
      req_node[req_no] = node;
      req_is_recv[req_no] = 0;
      req_no++;
@@ -6053,7 +6231,8 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
  while (pending_recv > 0)
  {
    int outcount = 0;
-    MPI_Waitsome(req_no, cache.reqs, &outcount, completed, cache.stats);
+    parallel_waitsome_checked(req_no, cache.reqs, &outcount, completed, cache.stats,
+                              "Parallel::transfermix_cached");
    if (outcount == MPI_UNDEFINED) break;

    for (int i = 0; i < outcount; i++)
@@ -6068,10 +6247,10 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
    }
  }

-  if (req_no > 0) MPI_Waitall(req_no, cache.reqs, cache.stats);
+  if (req_no > 0) parallel_waitall_checked(req_no, cache.reqs, cache.stats, "Parallel::transfermix_cached");

-  if (self_len > 0)
-    data_packermix(cache.recv_bufs[myrank], cache.combined_src[myrank], cache.combined_dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
+  if (self_len > 0)
+    data_packermix(cache.recv_bufs[myrank], cache.combined_src[myrank], cache.combined_dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);

  delete[] req_node;
  delete[] req_is_recv;
@@ -6787,7 +6966,10 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
    build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
  }

-  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  {
+    ParallelTransferContextGuard transfer_ctx("Parallel::Interp copy");
+    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+  }

  for (int node = 0; node < cpusize; node++)
  {
@@ -6826,7 +7008,10 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
        Sync(PatcL, FutureList, Symmetry);
      }
      //<~~~prolong then
-      transfer(transfer_src, transfer_dst, FutureList, FutureList, Symmetry);
+      {
+        ParallelTransferContextGuard transfer_ctx("Parallel::Interp prolong FutureList");
+        transfer(transfer_src, transfer_dst, FutureList, FutureList, Symmetry);
+      }

      // for StateList
      // time interpolation part
@@ -6842,7 +7027,10 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
        Sync(PatcL, tmList, Symmetry);
      }
      //<~~~prolong then
-      transfer(transfer_src, transfer_dst, tmList, StateList, Symmetry);
+      {
+        ParallelTransferContextGuard transfer_ctx("Parallel::Interp prolong StateList");
+        transfer(transfer_src, transfer_dst, tmList, StateList, Symmetry);
+      }
    }
    else
    {
@@ -6853,7 +7041,10 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
        Sync(PatcL, VarList, Symmetry);
      }
      //<~~~prolong then
-      transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+      {
+        ParallelTransferContextGuard transfer_ctx("Parallel::Interp prolong VarList");
+        transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
+      }
    }

    for (int node = 0; node < cpusize; node++)
--- a/AMSS_NCKU_source/Parallel.h
+++ b/AMSS_NCKU_source/Parallel.h
@@ -89,9 +89,12 @@ namespace Parallel
  void transfermix(MyList<gridseg> **src, MyList<gridseg> **dst,
                   MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /*target */,
                   int Symmetry);
-  void Sync(Patch *Pat, MyList<var> *VarList, int Symmetry);
-  void Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry);
-  void Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry);
+  void Sync(Patch *Pat, MyList<var> *VarList, int Symmetry);
+  void Sync(Patch *Pat, MyList<var> *VarList, int Symmetry, const char *context);
+  void Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry);
+  void Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, const char *context);
+  void Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry);
+  void Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, const char *context);

  struct SyncCache {
    bool valid;
@@ -121,14 +124,15 @@ namespace Parallel
                       MyList<var> *VarList1, MyList<var> *VarList2,
                       int Symmetry, SyncCache &cache);

-  struct AsyncSyncState {
-    int req_no;
-    bool active;
-    int *req_node;
-    int *req_is_recv;
-    int pending_recv;
-    AsyncSyncState() : req_no(0), active(false), req_node(0), req_is_recv(0), pending_recv(0) {}
-  };
+  struct AsyncSyncState {
+    int req_no;
+    bool active;
+    int mpi_tag;
+    int *req_node;
+    int *req_is_recv;
+    int pending_recv;
+    AsyncSyncState() : req_no(0), active(false), mpi_tag(0), req_node(0), req_is_recv(0), pending_recv(0) {}
+  };

  void Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry,
                  SyncCache &cache, AsyncSyncState &state);
--- a/AMSS_NCKU_source/bssn_class.C
+++ b/AMSS_NCKU_source/bssn_class.C
@@ -14,7 +14,8 @@ using namespace std;
 #include <string.h>
 #endif

-#include <time.h>
+#include <time.h>
+#include <unistd.h>

 #include "macrodef.h"
 #include "misc.h"
@@ -36,12 +37,18 @@ using namespace std;
 #include "myglobal.h"
 #endif

-#include "perf.h"
-
-#include "derivatives.h"
-#include "ricci_gamma.h"
-
-//================================================================================================
+#include "perf.h"
+
+#include "derivatives.h"
+#include "ricci_gamma.h"
+
+// Compile-time switch for per-timestep memory usage collection/printing.
+// Default is OFF to reduce overhead in production runs.
+#ifndef BSSN_ENABLE_MEM_USAGE_LOG
+#define BSSN_ENABLE_MEM_USAGE_LOG 0
+#endif
+
+//================================================================================================

 // define bssn_class

@@ -2030,9 +2037,10 @@ void bssn_class::Read_Ansorg()

 void bssn_class::Evolve(int Steps)
 {
-  clock_t prev_clock, curr_clock;
-  double LastDump = 0.0, LastCheck = 0.0, Last2dDump = 0.0;
-  LastAnas = 0;
+  clock_t prev_clock, curr_clock;
+  double LastDump = 0.0, LastCheck = 0.0, Last2dDump = 0.0;
+  LastAnas = 0;
+  LastConsOut = 0;
 #if 0
 //initial checkpoint for special uasge
     {
@@ -2129,8 +2137,10 @@ void bssn_class::Evolve(int Steps)
  #endif
  */
  
-  perf bssn_perf;
-  size_t current_min, current_avg, current_max, peak_min, peak_avg, peak_max;
+#if BSSN_ENABLE_MEM_USAGE_LOG
+  perf bssn_perf;
+  size_t current_min, current_avg, current_max, peak_min, peak_avg, peak_max;
+#endif

  for (int lev = 0; lev < GH->levels; lev++)
    GH->Lt[lev] = PhysTime;
@@ -2224,21 +2234,23 @@ void bssn_class::Evolve(int Steps)
 //		fgt(PhysTime-dT_mon,StartTime,dT_mon/2),ErrorMonitor);
 #endif

-    // Retrieve memory usage information used during computation; master process prints it
-    bssn_perf.MemoryUsage(&current_min, &current_avg, &current_max,
-                          &peak_min, &peak_avg, &peak_max, nprocs);
-    if (myrank == 0)
-    {
-      printf(" Memory usage: current %0.4lg/%0.4lg/%0.4lgMB, "
-             "peak %0.4lg/%0.4lg/%0.4lgMB\n",
-             (double)current_min / (1024.0 * 1024.0),
-             (double)current_avg / (1024.0 * 1024.0),
-             (double)current_max / (1024.0 * 1024.0),
-             (double)peak_min / (1024.0 * 1024.0),
-             (double)peak_avg / (1024.0 * 1024.0),
-             (double)peak_max / (1024.0 * 1024.0));
-      cout << endl;
-    }
+#if BSSN_ENABLE_MEM_USAGE_LOG
+    // Retrieve memory usage information used during computation; master process prints it
+    bssn_perf.MemoryUsage(&current_min, &current_avg, &current_max,
+                          &peak_min, &peak_avg, &peak_max, nprocs);
+    if (myrank == 0)
+    {
+      printf(" Memory usage: current %0.4lg/%0.4lg/%0.4lgMB, "
+             "peak %0.4lg/%0.4lg/%0.4lgMB\n",
+             (double)current_min / (1024.0 * 1024.0),
+             (double)current_avg / (1024.0 * 1024.0),
+             (double)current_max / (1024.0 * 1024.0),
+             (double)peak_min / (1024.0 * 1024.0),
+             (double)peak_avg / (1024.0 * 1024.0),
+             (double)peak_max / (1024.0 * 1024.0));
+      cout << endl;
+    }
+#endif
    
    // Output puncture positions at each step
    if (myrank == 0)
@@ -2286,18 +2298,21 @@ void bssn_class::Evolve(int Steps)
    ////////////////////////////////////////////////////////////

    // When LastCheck >= CheckTime, perform runtime checks and output status data
-    if (LastCheck >= CheckTime)
-    {
-      LastCheck = 0;
-
-      CheckPoint->write_Black_Hole_position(BH_num_input, BH_num, Porg0, Porgbr, Mass);
+    if (LastCheck >= CheckTime)
+    {
+      LastCheck = 0;
+
+      CheckPoint->write_Black_Hole_position(BH_num_input, BH_num, Porg0, Porgbr, Mass);
      CheckPoint->writecheck_cgh(PhysTime, GH);
 #ifdef WithShell
      CheckPoint->writecheck_sh(PhysTime, SH);
-#endif
-      CheckPoint->write_bssn(LastDump, Last2dDump, LastAnas);
-    }
-  }
+#endif
+      CheckPoint->write_bssn(LastDump, Last2dDump, LastAnas);
+    }
+
+    // Keep output/analysis phases aligned across ranks before the next coarse step.
+    MPI_Barrier(MPI_COMM_WORLD);
+  }
  /*
  #ifdef With_AHF
  // final apparent horizon finding
@@ -3016,9 +3031,14 @@ void bssn_class::RecursiveStep(int lev, int num) // in all 2^(lev+1)-1 steps

 #if (PSTR == 0)
 #if 1
-void bssn_class::Step(int lev, int YN)
-{
-  setpbh(BH_num, Porg0, Mass, BH_num_input);
+void bssn_class::Step(int lev, int YN)
+{
+#ifdef USE_GPU
+  Step_MainPath_GPU(lev, YN);
+  return;
+#endif
+
+  setpbh(BH_num, Porg0, Mass, BH_num_input);

  double dT_lev = dT * pow(0.5, Mymax(lev, trfls));

@@ -6238,7 +6258,7 @@ for(int ilev = GH->levels-1;ilev>=lev;ilev--)
 for(int ilev=GH->levels-1;ilev>lev;ilev--)
    RestrictProlong(ilev,1,false,DG_List,DG_List,DG_List);
 #else
-  Parallel::Sync(GH->PatL[lev], DG_List, Symmetry);
+  Parallel::Sync(GH->PatL[lev], DG_List, Symmetry, "bssn_class::Compute_Psi4");
 #endif

 #ifdef WithShell
@@ -6893,10 +6913,10 @@ void bssn_class::AnalysisStuff(int lev, double dT_lev)
 {
  LastAnas += dT_lev;

-  if (LastAnas >= AnasTime)
-  {
-#ifdef Point_Psi4
-#error "not support parallel levels yet"
+  if (LastAnas >= AnasTime)
+  {
+#ifdef Point_Psi4
+#error "not support parallel levels yet"
    // Gam_ijk and R_ij have been calculated in Interp_Constraint()
    double SYM = 1, ANT = -1;
    for (int levh = lev; levh < GH->levels; levh++)
@@ -7240,9 +7260,9 @@ void bssn_class::AnalysisStuff(int lev, double dT_lev)

 //================================================================================================

-void bssn_class::Constraint_Out()
-{
-  LastConsOut += dT * pow(0.5, Mymax(0, trfls));
+void bssn_class::Constraint_Out()
+{
+  LastConsOut += dT * pow(0.5, Mymax(0, trfls));

  if (LastConsOut >= AnasTime)
  // Constraint violation
@@ -7307,7 +7327,7 @@ void bssn_class::Constraint_Out()
          Pp = Pp->next;
        }
      }
-      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry);
+      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry, "bssn_class::Constraint_Out[level]");
    }
 #ifdef WithShell
    if (0) // if the constrait quantities can be reused from the step rhs calculation
@@ -7529,7 +7549,7 @@ void bssn_class::AH_Prepare_derivatives()
      }
      Pp = Pp->next;
    }
-    Parallel::Sync(GH->PatL[lev], AHDList, Symmetry);
+    Parallel::Sync(GH->PatL[lev], AHDList, Symmetry, "bssn_class::AH_Prepare_derivatives");
  }
 }

@@ -7810,7 +7830,7 @@ void bssn_class::Interp_Constraint(bool infg)
          Pp = Pp->next;
        }
      }
-      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry);
+      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry, "bssn_class::Interp_Constraint[level]");
    }
 #ifdef WithShell
    if (0) // if the constrait quantities can be reused from the step rhs calculation
@@ -8068,7 +8088,7 @@ void bssn_class::Compute_Constraint()
        Pp = Pp->next;
      }
    }
-    Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry);
+    Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry, "bssn_class::Compute_Constraint[level]");
  }
  // prolong restrict constraint quantities
  for (lev = GH->levels - 1; lev > 0; lev--)
@@ -8381,12 +8401,18 @@ void bssn_class::Enforce_algcon(int lev, int fg)

 //================================================================================================

-bool bssn_class::check_Stdin_Abort() 
-{
-
-    fd_set readfds;
-
-    struct timeval timeout;
+bool bssn_class::check_Stdin_Abort() 
+{
+    // Non-interactive launches (mpirun via Python/subprocess, batch jobs, redirected stdin)
+    // should not probe stdin. Some MPI runtimes treat stdin as a managed channel and can
+    // fail when rank 0 polls/consumes it.
+    if (!isatty(STDIN_FILENO)) {
+        return false;
+    }
+
+    fd_set readfds;
+
+    struct timeval timeout;

    FD_ZERO(&readfds);
    FD_SET(STDIN_FILENO, &readfds);
@@ -8395,14 +8421,17 @@ bool bssn_class::check_Stdin_Abort()
    timeout.tv_sec = 0;
    timeout.tv_usec = 0;

-    int activity = select(STDIN_FILENO + 1, &readfds, nullptr, nullptr, &timeout);
-
-    if (activity > 0 && FD_ISSET(STDIN_FILENO, &readfds)) {
-        string input_abort;
-        if (cin >> input_abort) {
-            if (input_abort == "stop") {
-                return true;
-            }
+    int activity = select(STDIN_FILENO + 1, &readfds, nullptr, nullptr, &timeout);
+    if (activity <= 0) {
+        return false;
+    }
+
+    if (FD_ISSET(STDIN_FILENO, &readfds)) {
+        string input_abort;
+        if (cin >> input_abort) {
+            if (input_abort == "stop") {
+                return true;
+            }
        }
    }

--- a/AMSS_NCKU_source/bssn_class.h
+++ b/AMSS_NCKU_source/bssn_class.h
@@ -171,16 +171,19 @@ public:
       
       bool check_Stdin_Abort(); 

-       virtual void Setup_Initial_Data_Cao();
-       virtual void Setup_Initial_Data_Lousto();
-       virtual void Initialize();
-       virtual void Read_Ansorg();
-       virtual void Read_Pablo() {};
-       virtual void Compute_Psi4(int lev);
-       virtual void Step(int lev, int YN);
-       virtual void Interp_Constraint(bool infg);
-       virtual void Constraint_Out();
-       virtual void Compute_Constraint();
+       virtual void Setup_Initial_Data_Cao();
+       virtual void Setup_Initial_Data_Lousto();
+       virtual void Initialize();
+       virtual void Read_Ansorg();
+       virtual void Read_Pablo() {};
+       virtual void Compute_Psi4(int lev);
+       virtual void Step(int lev, int YN);
+#ifdef USE_GPU
+       void Step_MainPath_GPU(int lev, int YN);
+#endif
+       virtual void Interp_Constraint(bool infg);
+       virtual void Constraint_Out();
+       virtual void Compute_Constraint();

 #ifdef With_AHF
 protected:
--- a/AMSS_NCKU_source/bssn_cuda_ops.cu
+++ b/AMSS_NCKU_source/bssn_cuda_ops.cu
--- a/AMSS_NCKU_source/bssn_cuda_ops.h
+++ b/AMSS_NCKU_source/bssn_cuda_ops.h
@@ -0,0 +1,45 @@
+#ifndef BSSN_CUDA_OPS_H
+#define BSSN_CUDA_OPS_H
+
+int bssn_cuda_enforce_ga(int *ex,
+                         double *dxx, double *gxy, double *gxz,
+                         double *dyy, double *gyz, double *dzz,
+                         double *Axx, double *Axy, double *Axz,
+                         double *Ayy, double *Ayz, double *Azz);
+
+int bssn_cuda_rk4_boundary_var(int *ex, double dT,
+                               const double *X, const double *Y, const double *Z,
+                               double xmin, double ymin, double zmin,
+                               double xmax, double ymax, double zmax,
+                               const double *state0,
+                               const double *boundary_src,
+                               double *stage_data,
+                               double *rhs_accum,
+                               double propspeed,
+                               const double SoA[3],
+                               int symmetry,
+                               int lev,
+                               int rk_stage,
+                               bool download_to_host = true);
+
+int bssn_cuda_lowerbound(int *ex, double *chi, double tinny, bool download_to_host = true);
+int bssn_cuda_download_buffer(int *ex, double *host_ptr);
+
+int bssn_cuda_prolong3_pack(int wei,
+                            const double *llbc, const double *uubc, const int *extc, const double *func,
+                            const double *llbf, const double *uubf, const int *extf, double *funf,
+                            const double *llbp, const double *uubp,
+                            const double *SoA, int symmetry);
+
+int bssn_cuda_interp_points_batch(const int *ex,
+                                  const double *X, const double *Y, const double *Z,
+                                  const double *const *fields,
+                                  const double *soa_flat,
+                                  int num_var,
+                                  const double *px, const double *py, const double *pz,
+                                  int num_points,
+                                  int ordn,
+                                  int symmetry,
+                                  double *out);
+
+#endif
--- a/AMSS_NCKU_source/bssn_cuda_step.C
+++ b/AMSS_NCKU_source/bssn_cuda_step.C
@@ -0,0 +1,479 @@
+#include "macrodef.h"
+
+#ifdef USE_GPU
+
+#include <cmath>
+#include <vector>
+
+#include "bssn_class.h"
+#include "bssn_cuda_ops.h"
+#include "bssn_gpu.h"
+#include "bssn_macro.h"
+#include "rungekutta4_rout.h"
+
+void bssn_class::Step_MainPath_GPU(int lev, int YN)
+{
+#ifdef WithShell
+#error "Step_MainPath_GPU currently supports Patch grids only."
+#endif
+
+  if (bssn_gpu_bind_process_device(myrank))
+  {
+    cerr << "GPU device bind failure on MPI rank " << myrank << endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+  bssn_gpu_clear_cached_device_buffers();
+
+  setpbh(BH_num, Porg0, Mass, BH_num_input);
+
+  const double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
+
+#if (MAPBH == 1)
+  if (BH_num > 0 && lev == GH->levels - 1)
+  {
+    compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
+    for (int ithBH = 0; ithBH < BH_num; ithBH++)
+    {
+      for (int ith = 0; ith < 3; ith++)
+        Porg1[ithBH][ith] = Porg0[ithBH][ith] + Porg_rhs[ithBH][ith] * dT_lev;
+      if (Symmetry > 0)
+        Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
+      if (Symmetry == 2)
+      {
+        Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
+        Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
+      }
+    }
+  }
+
+  if (lev == a_lev)
+    AnalysisStuff(lev, dT_lev);
+#endif
+
+#ifdef With_AHF
+  AH_Step_Find(lev, dT_lev);
+#endif
+
+  const bool BB = fgt(PhysTime, StartTime, dT_lev / 2);
+  (void)BB;
+  double ndeps = (lev < GH->movls) ? numepsb : numepss;
+  double TRK4 = PhysTime;
+  int iter_count = 0;
+  int pre = 0, cor = 1;
+  int ERROR = 0;
+
+  auto run_stage_on_block =
+      [&](Block *cg, Patch *patch, MyList<var> *state0_list,
+          MyList<var> *boundary_src_list, MyList<var> *stage_data_list,
+          MyList<var> *rhs_list, int rk_stage) {
+        MyList<var> *varl0 = state0_list;
+        MyList<var> *varlb = boundary_src_list;
+        MyList<var> *varls = stage_data_list;
+        MyList<var> *varlr = rhs_list;
+
+        while (varl0)
+        {
+          if (bssn_cuda_rk4_boundary_var(cg->shape, dT_lev,
+                                         cg->X[0], cg->X[1], cg->X[2],
+                                         patch->bbox[0], patch->bbox[1], patch->bbox[2],
+                                         patch->bbox[3], patch->bbox[4], patch->bbox[5],
+                                         cg->fgfs[varl0->data->sgfn],
+                                         cg->fgfs[varlb->data->sgfn],
+                                         cg->fgfs[varls->data->sgfn],
+                                         cg->fgfs[varlr->data->sgfn],
+                                         varl0->data->propspeed,
+                                         varl0->data->SoA,
+                                         Symmetry, lev, rk_stage, false))
+          {
+            cerr << "GPU rk4/boundary failure: lev=" << lev
+                 << " rk_stage=" << rk_stage
+                 << " var=" << varl0->data->name
+                 << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                 << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                 << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+            ERROR = 1;
+            break;
+          }
+          varl0 = varl0->next;
+          varlb = varlb->next;
+          varls = varls->next;
+          varlr = varlr->next;
+        }
+      };
+
+  auto stage_download_var_list =
+      [&](Block *cg, MyList<var> *var_list) {
+        while (var_list)
+        {
+          if (bssn_cuda_download_buffer(cg->shape, cg->fgfs[var_list->data->sgfn]))
+          {
+            cerr << "GPU stage download failure: lev=" << lev
+                 << " var=" << var_list->data->name
+                 << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                 << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                 << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+            ERROR = 1;
+            break;
+          }
+          var_list = var_list->next;
+        }
+      };
+
+  auto stage_upload_var_list =
+      [&](Block *cg, MyList<var> *var_list) {
+        const int n = cg->shape[0] * cg->shape[1] * cg->shape[2];
+        while (var_list)
+        {
+          if (bssn_gpu_stage_upload_buffer(cg->fgfs[var_list->data->sgfn], n))
+          {
+            cerr << "GPU state upload failure: lev=" << lev
+                 << " var=" << var_list->data->name
+                 << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                 << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                 << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+            ERROR = 1;
+            break;
+          }
+          var_list = var_list->next;
+        }
+      };
+
+  auto ensure_stage_device_var_list =
+      [&](Block *cg, MyList<var> *var_list) {
+        const int n = cg->shape[0] * cg->shape[1] * cg->shape[2];
+        while (var_list)
+        {
+          double *host_ptr = cg->fgfs[var_list->data->sgfn];
+          if (!bssn_gpu_find_device_buffer(host_ptr) &&
+              bssn_gpu_stage_upload_buffer(host_ptr, n))
+          {
+            cerr << "GPU state ensure failure: lev=" << lev
+                 << " var=" << var_list->data->name
+                 << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                 << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                 << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+            ERROR = 1;
+            break;
+          }
+          var_list = var_list->next;
+        }
+      };
+
+  auto refresh_synced_device_regions =
+      [&](Block *cg, MyList<var> *var_list, Parallel::SyncCache &cache) {
+        std::vector<Parallel::gridseg *> local_segments;
+        for (int node = 0; node < cache.cpusize; ++node)
+        {
+          MyList<Parallel::gridseg> *seg = cache.combined_dst[node];
+          while (seg)
+          {
+            if (seg->data && seg->data->Bg == cg)
+              local_segments.push_back(seg->data);
+            seg = seg->next;
+          }
+        }
+
+        if (local_segments.empty())
+          return;
+
+        const int n = cg->shape[0] * cg->shape[1] * cg->shape[2];
+        while (var_list)
+        {
+          double *host_ptr = cg->fgfs[var_list->data->sgfn];
+          if (!bssn_gpu_find_device_buffer(host_ptr))
+          {
+            if (bssn_gpu_stage_upload_buffer(host_ptr, n))
+            {
+              cerr << "GPU sync refresh upload failure: lev=" << lev
+                   << " var=" << var_list->data->name
+                   << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                   << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                   << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+              ERROR = 1;
+              break;
+            }
+          }
+          else
+          {
+            for (size_t i = 0; i < local_segments.size(); ++i)
+            {
+              Parallel::gridseg *seg = local_segments[i];
+              if (bssn_gpu_stage_upload_region(host_ptr,
+                                               cg->shape,
+                                               cg->bbox,
+                                               cg->bbox + dim,
+                                               seg->shape,
+                                               seg->llb))
+              {
+                cerr << "GPU sync region refresh failure: lev=" << lev
+                     << " var=" << var_list->data->name
+                     << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                     << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                     << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+                ERROR = 1;
+                break;
+              }
+            }
+            if (ERROR)
+              break;
+          }
+          var_list = var_list->next;
+        }
+      };
+
+  auto refresh_stage_device_after_sync =
+      [&](MyList<var> *var_list, Parallel::SyncCache &cache) {
+        MyList<Patch> *patch_it = GH->PatL[lev];
+        while (patch_it)
+        {
+          MyList<Block> *block_it = patch_it->data->blb;
+          while (block_it)
+          {
+            Block *cg = block_it->data;
+            if (myrank == cg->rank)
+              refresh_synced_device_regions(cg, var_list, cache);
+
+            if (block_it == patch_it->data->ble)
+              break;
+            block_it = block_it->next;
+          }
+          if (ERROR)
+            break;
+          patch_it = patch_it->next;
+        }
+      };
+
+  MyList<Patch> *Pp = GH->PatL[lev];
+  while (Pp)
+  {
+    MyList<Block> *BP = Pp->data->blb;
+    while (BP)
+    {
+      Block *cg = BP->data;
+      if (myrank == cg->rank)
+      {
+        stage_upload_var_list(cg, StateList);
+        if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_FIRST_TIME))
+          ERROR = 1;
+
+        run_stage_on_block(cg, Pp->data, StateList, StateList, SynchList_pre, RHSList, iter_count);
+
+        if (bssn_cuda_lowerbound(cg->shape, cg->fgfs[phi->sgfn], chitiny, false))
+        {
+          cerr << "GPU lowerbound failure: lev=" << lev
+               << " rk_stage=" << iter_count
+               << " var=" << phi->name
+               << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+               << cg->bbox[1] << ":" << cg->bbox[4] << ","
+               << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+          ERROR = 1;
+        }
+        if (!ERROR)
+          stage_download_var_list(cg, SynchList_pre);
+      }
+      if (BP == Pp->data->ble)
+        break;
+      BP = BP->next;
+    }
+    Pp = Pp->next;
+  }
+
+  MPI_Request err_req_pre;
+  {
+    int erh = ERROR;
+    MPI_Iallreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD, &err_req_pre);
+  }
+
+  Parallel::AsyncSyncState async_pre;
+  Parallel::Sync_start(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev], async_pre);
+  Parallel::Sync_finish(sync_cache_pre[lev], async_pre, SynchList_pre, Symmetry);
+  if (!ERROR)
+    refresh_stage_device_after_sync(SynchList_pre, sync_cache_pre[lev]);
+
+  MPI_Wait(&err_req_pre, MPI_STATUS_IGNORE);
+  if (ERROR)
+  {
+    Parallel::Dump_Data(GH->PatL[lev], StateList, 0, PhysTime, dT_lev);
+    if (myrank == 0)
+    {
+      if (ErrorMonitor->outfile)
+        ErrorMonitor->outfile << "find NaN in state variables at t = " << PhysTime
+                              << ", lev = " << lev << endl;
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+  }
+
+#if (MAPBH == 0)
+  if (BH_num > 0 && lev == GH->levels - 1)
+  {
+    compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
+    for (int ithBH = 0; ithBH < BH_num; ithBH++)
+    {
+      f_rungekutta4_scalar(dT_lev, Porg0[ithBH][0], Porg[ithBH][0], Porg_rhs[ithBH][0], iter_count);
+      f_rungekutta4_scalar(dT_lev, Porg0[ithBH][1], Porg[ithBH][1], Porg_rhs[ithBH][1], iter_count);
+      f_rungekutta4_scalar(dT_lev, Porg0[ithBH][2], Porg[ithBH][2], Porg_rhs[ithBH][2], iter_count);
+      if (Symmetry > 0)
+        Porg[ithBH][2] = fabs(Porg[ithBH][2]);
+      if (Symmetry == 2)
+      {
+        Porg[ithBH][0] = fabs(Porg[ithBH][0]);
+        Porg[ithBH][1] = fabs(Porg[ithBH][1]);
+      }
+    }
+  }
+
+  if (lev == a_lev)
+    AnalysisStuff(lev, dT_lev);
+#endif
+
+  for (iter_count = 1; iter_count < 4; iter_count++)
+  {
+    if (iter_count == 1 || iter_count == 3)
+      TRK4 += dT_lev / 2;
+
+    Pp = GH->PatL[lev];
+    while (Pp)
+    {
+      MyList<Block> *BP = Pp->data->blb;
+      while (BP)
+      {
+        Block *cg = BP->data;
+        if (myrank == cg->rank)
+        {
+          ensure_stage_device_var_list(cg, SynchList_pre);
+          if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_THEN))
+            ERROR = 1;
+
+          run_stage_on_block(cg, Pp->data, StateList, SynchList_pre, SynchList_cor, RHSList, iter_count);
+
+          if (bssn_cuda_lowerbound(cg->shape, cg->fgfs[phi1->sgfn], chitiny, false))
+          {
+            cerr << "GPU lowerbound failure: lev=" << lev
+                 << " rk_stage=" << iter_count
+                 << " var=" << phi1->name
+                 << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                 << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                 << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+            ERROR = 1;
+          }
+          if (!ERROR)
+            stage_download_var_list(cg, SynchList_cor);
+        }
+
+        if (BP == Pp->data->ble)
+          break;
+        BP = BP->next;
+      }
+      Pp = Pp->next;
+    }
+
+    MPI_Request err_req_cor;
+    {
+      int erh = ERROR;
+      MPI_Iallreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD, &err_req_cor);
+    }
+
+    Parallel::AsyncSyncState async_cor;
+    Parallel::Sync_start(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev], async_cor);
+    Parallel::Sync_finish(sync_cache_cor[lev], async_cor, SynchList_cor, Symmetry);
+    if (!ERROR && iter_count < 3)
+      refresh_stage_device_after_sync(SynchList_cor, sync_cache_cor[lev]);
+
+    MPI_Wait(&err_req_cor, MPI_STATUS_IGNORE);
+    if (ERROR)
+    {
+      Parallel::Dump_Data(GH->PatL[lev], SynchList_pre, 0, PhysTime, dT_lev);
+      if (myrank == 0)
+      {
+        if (ErrorMonitor->outfile)
+          ErrorMonitor->outfile << "find NaN in RK4 substep#" << iter_count
+                                << " variables at t = " << PhysTime
+                                << ", lev = " << lev << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
+      }
+    }
+
+#if (MAPBH == 0)
+    if (BH_num > 0 && lev == GH->levels - 1)
+    {
+      compute_Porg_rhs(Porg, Porg1, Sfx, Sfy, Sfz, lev);
+      for (int ithBH = 0; ithBH < BH_num; ithBH++)
+      {
+        f_rungekutta4_scalar(dT_lev, Porg0[ithBH][0], Porg1[ithBH][0], Porg_rhs[ithBH][0], iter_count);
+        f_rungekutta4_scalar(dT_lev, Porg0[ithBH][1], Porg1[ithBH][1], Porg_rhs[ithBH][1], iter_count);
+        f_rungekutta4_scalar(dT_lev, Porg0[ithBH][2], Porg1[ithBH][2], Porg_rhs[ithBH][2], iter_count);
+        if (Symmetry > 0)
+          Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
+        if (Symmetry == 2)
+        {
+          Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
+          Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
+        }
+      }
+    }
+#endif
+
+    if (iter_count < 3)
+    {
+      Pp = GH->PatL[lev];
+      while (Pp)
+      {
+        MyList<Block> *BP = Pp->data->blb;
+        while (BP)
+        {
+          BP->data->swapList(SynchList_pre, SynchList_cor, myrank);
+          if (BP == Pp->data->ble)
+            break;
+          BP = BP->next;
+        }
+        Pp = Pp->next;
+      }
+
+#if (MAPBH == 0)
+      if (BH_num > 0 && lev == GH->levels - 1)
+      {
+        for (int ithBH = 0; ithBH < BH_num; ithBH++)
+        {
+          Porg[ithBH][0] = Porg1[ithBH][0];
+          Porg[ithBH][1] = Porg1[ithBH][1];
+          Porg[ithBH][2] = Porg1[ithBH][2];
+        }
+      }
+#endif
+    }
+  }
+
+#if (RPS == 0)
+  RestrictProlong(lev, YN, BB);
+#endif
+
+  bssn_gpu_clear_cached_device_buffers();
+
+  Pp = GH->PatL[lev];
+  while (Pp)
+  {
+    MyList<Block> *BP = Pp->data->blb;
+    while (BP)
+    {
+      Block *cg = BP->data;
+      cg->swapList(StateList, SynchList_cor, myrank);
+      cg->swapList(OldStateList, SynchList_cor, myrank);
+      if (BP == Pp->data->ble)
+        break;
+      BP = BP->next;
+    }
+    Pp = Pp->next;
+  }
+
+  if (BH_num > 0 && lev == GH->levels - 1)
+  {
+    for (int ithBH = 0; ithBH < BH_num; ithBH++)
+    {
+      Porg0[ithBH][0] = Porg1[ithBH][0];
+      Porg0[ithBH][1] = Porg1[ithBH][1];
+      Porg0[ithBH][2] = Porg1[ithBH][2];
+    }
+  }
+}
+
+#endif
--- a/AMSS_NCKU_source/bssn_gpu.cu
+++ b/AMSS_NCKU_source/bssn_gpu.cu
--- a/AMSS_NCKU_source/bssn_gpu.h
+++ b/AMSS_NCKU_source/bssn_gpu.h
@@ -4,10 +4,8 @@
 #include "bssn_macro.h"
 #include "macrodef.fh"

-#define DEVICE_ID 0
-// #define DEVICE_ID_BY_MPI_RANK
-#define GRID_DIM 256
-#define BLOCK_DIM 128
+#define GRID_DIM 256
+#define BLOCK_DIM 128

 #define _FH2_(i, j, k) fh[(i) + (j) * _1D_SIZE[2] + (k) * _2D_SIZE[2]]
 #define _FH3_(i, j, k) fh[(i) + (j) * _1D_SIZE[3] + (k) * _2D_SIZE[3]]
@@ -65,9 +63,22 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,
            double *Gmx_Res, double *Gmy_Res, double *Gmz_Res,
            int &Symmetry, int &Lev, double &eps, int &co);

-int gpu_rhs_ss(RHS_SS_PARA);
-
-/** Init GPU side data in GPUMeta. */
-// void init_fluid_meta_gpu(GPUMeta *gpu_meta);
+int gpu_rhs_ss(RHS_SS_PARA);
+
+int bssn_gpu_bind_process_device(int mpi_rank);
+void bssn_gpu_clear_cached_device_buffers();
+const double *bssn_gpu_find_device_buffer(const double *host_ptr);
+void bssn_gpu_register_device_buffer(const double *host_ptr, const double *device_ptr);
+int bssn_gpu_stage_upload_buffer(const double *host_ptr, int count);
+int bssn_gpu_stage_zero_buffer(const double *host_ptr, int count);
+int bssn_gpu_stage_upload_region(const double *host_ptr,
+                                 const int *full_shape,
+                                 const double *full_llb,
+                                 const double *full_uub,
+                                 const int *region_shape,
+                                 const double *region_llb);
+
+/** Init GPU side data in GPUMeta. */
+// void init_fluid_meta_gpu(GPUMeta *gpu_meta);

 #endif
--- a/AMSS_NCKU_source/bssn_rhs.f90
+++ b/AMSS_NCKU_source/bssn_rhs.f90
@@ -59,9 +59,10 @@
  real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Rxx,Rxy,Rxz,Ryy,Ryz,Rzz
  real*8,intent(in) :: eps
  real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: ham_Res, movx_Res, movy_Res, movz_Res
-  real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gmx_Res, Gmy_Res, Gmz_Res
-!  gont = 0: success; gont = 1: something wrong
-  integer::gont
+  real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gmx_Res, Gmy_Res, Gmz_Res
+!  gont = 0: success; gont = 1: something wrong
+  integer::gont
+  integer :: i,j,k

 !~~~~~~> Other variables:

@@ -83,11 +84,18 @@
  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz
  real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz

-  real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA
-  real*8            :: dX, dY, dZ, PI
-  real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0
-  real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0
-  real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0
+  real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA
+  real*8            :: dX, dY, dZ, PI
+  real*8            :: divb_loc,det_loc
+  real*8            :: gupxx_loc,gupxy_loc,gupxz_loc,gupyy_loc,gupyz_loc,gupzz_loc
+  real*8            :: Rxx_loc,Rxy_loc,Rxz_loc,Ryy_loc,Ryz_loc,Rzz_loc
+  real*8            :: fxx_loc,fxy_loc,fxz_loc
+  real*8            :: Gamxa_loc,Gamya_loc,Gamza_loc
+  real*8            :: f_loc,chin_loc
+  real*8            :: l_fxx,l_fxy,l_fxz,l_fyy,l_fyz,l_fzz,S_loc
+  real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0
+  real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0
+  real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0
  double precision,parameter::FF = 0.75d0,eta=2.d0
  real*8, parameter :: F1o3 = 1.D0/3.D0, F2o3 = 2.D0/3.D0,F3o2=1.5d0, F1o6 = 1.D0/6.D0
  real*8, parameter :: F16=1.6d1,F8=8.d0
@@ -96,11 +104,11 @@
  real*8, dimension(ex(1),ex(2),ex(3)) :: reta
 #endif

-#if (GAUGE == 6 || GAUGE == 7)
-  integer :: BHN,i,j,k
-  real*8, dimension(9) :: Porg
-  real*8, dimension(3) :: Mass
-  real*8 :: r1,r2,M,A,w1,w2,C1,C2
+#if (GAUGE == 6 || GAUGE == 7)
+  integer :: BHN
+  real*8, dimension(9) :: Porg
+  real*8, dimension(3) :: Mass
+  real*8 :: r1,r2,M,A,w1,w2,C1,C2
  real*8, dimension(ex(1),ex(2),ex(3)) :: reta

  call getpbh(BHN,Porg,Mass)
@@ -145,174 +153,204 @@
  dY = Y(2) - Y(1)
  dZ = Z(2) - Z(1)

-  alpn1 = Lap + ONE
-  chin1 = chi + ONE
-  gxx = dxx + ONE
-  gyy = dyy + ONE
-  gzz = dzz + ONE
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+    alpn1(i,j,k) = Lap(i,j,k) + ONE
+    chin1(i,j,k) = chi(i,j,k) + ONE
+    gxx(i,j,k) = dxx(i,j,k) + ONE
+    gyy(i,j,k) = dyy(i,j,k) + ONE
+    gzz(i,j,k) = dzz(i,j,k) + ONE
+  enddo
+  enddo
+  enddo

  call fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev)
  call fderivs(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev)
  call fderivs(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev)
  
-  div_beta = betaxx + betayy + betazz
- 
-  call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
+  call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)

-  chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi
-
-  call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
-  call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev)
-  call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev)
-  call fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
-  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)
-
-  gxx_rhs = - TWO * alpn1 * Axx    -  F2o3 * gxx * div_beta          + &
-              TWO *(  gxx * betaxx +   gxy * betayx +   gxz * betazx)
-
-  gyy_rhs = - TWO * alpn1 * Ayy    -  F2o3 * gyy * div_beta          + &
-              TWO *(  gxy * betaxy +   gyy * betayy +   gyz * betazy)
-
-  gzz_rhs = - TWO * alpn1 * Azz    -  F2o3 * gzz * div_beta          + &
-              TWO *(  gxz * betaxz +   gyz * betayz +   gzz * betazz)
-
-  gxy_rhs = - TWO * alpn1 * Axy    +  F1o3 * gxy    * div_beta       + &
-                      gxx * betaxy                  +   gxz * betazy + &
-                                       gyy * betayx +   gyz * betazx   &
-                                                    -   gxy * betazz
-
-  gyz_rhs = - TWO * alpn1 * Ayz    +  F1o3 * gyz    * div_beta       + &
-                      gxy * betaxz +   gyy * betayz                  + &
-                      gxz * betaxy                  +   gzz * betazy   &
-                                                    -   gyz * betaxx
- 
-  gxz_rhs = - TWO * alpn1 * Axz    +  F1o3 * gxz    * div_beta       + &
-                      gxx * betaxz +   gxy * betayz                  + &
-                                       gyz * betayx +   gzz * betazx   &
-                                                    -   gxz * betayy     !rhs for gij
-
-! invert tilted metric
-  gupzz =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-           gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
-  gupxx =   ( gyy * gzz - gyz * gyz ) / gupzz
-  gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz
-  gupxz =   ( gxy * gyz - gyy * gxz ) / gupzz
-  gupyy =   ( gxx * gzz - gxz * gxz ) / gupzz
-  gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
-  gupzz =   ( gxx * gyy - gxy * gxy ) / gupzz
-
-  if(co == 0)then
-! Gam^i_Res = Gam^i + gup^ij_,j
-  Gmx_Res = Gamx - (gupxx*(gupxx*gxxx+gupxy*gxyx+gupxz*gxzx)&
-                   +gupxy*(gupxx*gxyx+gupxy*gyyx+gupxz*gyzx)&
-                   +gupxz*(gupxx*gxzx+gupxy*gyzx+gupxz*gzzx)&
-                   +gupxx*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)&
-                   +gupxy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)&
-                   +gupxz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)&
-                   +gupxx*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
-                   +gupxy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
-                   +gupxz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
-  Gmy_Res = Gamy - (gupxx*(gupxy*gxxx+gupyy*gxyx+gupyz*gxzx)&
-                   +gupxy*(gupxy*gxyx+gupyy*gyyx+gupyz*gyzx)&
-                   +gupxz*(gupxy*gxzx+gupyy*gyzx+gupyz*gzzx)&
-                   +gupxy*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)&
-                   +gupyy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)&
-                   +gupyz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)&
-                   +gupxy*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
-                   +gupyy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
-                   +gupyz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
-  Gmz_Res = Gamz - (gupxx*(gupxz*gxxx+gupyz*gxyx+gupzz*gxzx)&
-                   +gupxy*(gupxz*gxyx+gupyz*gyyx+gupzz*gyzx)&
-                   +gupxz*(gupxz*gxzx+gupyz*gyzx+gupzz*gzzx)&
-                   +gupxy*(gupxz*gxxy+gupyz*gxyy+gupzz*gxzy)&
-                   +gupyy*(gupxz*gxyy+gupyz*gyyy+gupzz*gyzy)&
-                   +gupyz*(gupxz*gxzy+gupyz*gyzy+gupzz*gzzy)&
-                   +gupxz*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
-                   +gupyz*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
-                   +gupzz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
-  endif
-
-! second kind of connection
-  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 ))
- 
-  Gamxyy =HALF*( gupxx*(TWO*gxyy - gyyx ) + gupxy*gyyy + gupxz*(TWO*gyzy - gyyz ))
-  Gamyyy =HALF*( gupxy*(TWO*gxyy - gyyx ) + gupyy*gyyy + gupyz*(TWO*gyzy - gyyz ))
-  Gamzyy =HALF*( gupxz*(TWO*gxyy - gyyx ) + gupyz*gyyy + gupzz*(TWO*gyzy - gyyz ))
-
-  Gamxzz =HALF*( gupxx*(TWO*gxzz - gzzx ) + gupxy*(TWO*gyzz - gzzy ) + gupxz*gzzz)
-  Gamyzz =HALF*( gupxy*(TWO*gxzz - gzzx ) + gupyy*(TWO*gyzz - gzzy ) + gupyz*gzzz)
-  Gamzzz =HALF*( gupxz*(TWO*gxzz - gzzx ) + gupyz*(TWO*gyzz - gzzy ) + gupzz*gzzz)
-
-  Gamxxy =HALF*( gupxx*gxxy + gupxy*gyyx + gupxz*( gxzy + gyzx - gxyz ) )
-  Gamyxy =HALF*( gupxy*gxxy + gupyy*gyyx + gupyz*( gxzy + gyzx - gxyz ) )
-  Gamzxy =HALF*( gupxz*gxxy + gupyz*gyyx + gupzz*( gxzy + gyzx - gxyz ) )
-
-  Gamxxz =HALF*( gupxx*gxxz + gupxy*( gxyz + gyzx - gxzy ) + gupxz*gzzx )
-  Gamyxz =HALF*( gupxy*gxxz + gupyy*( gxyz + gyzx - gxzy ) + gupyz*gzzx )
-  Gamzxz =HALF*( gupxz*gxxz + gupyz*( gxyz + gyzx - gxzy ) + gupzz*gzzx )
-
-  Gamxyz =HALF*( gupxx*( gxyz + gxzy - gyzx ) + gupxy*gyyz + gupxz*gzzy )
-  Gamyyz =HALF*( gupxy*( gxyz + gxzy - gyzx ) + gupyy*gyyz + gupyz*gzzy )
-  Gamzyz =HALF*( gupxz*( gxyz + gxzy - gyzx ) + gupyz*gyyz + gupzz*gzzy )
-! Raise indices of \tilde A_{ij} and store in R_ij
-
-  Rxx =    gupxx * gupxx * Axx + gupxy * gupxy * Ayy + gupxz * gupxz * Azz + &
-      TWO*(gupxx * gupxy * Axy + gupxx * gupxz * Axz + gupxy * gupxz * Ayz)
-
-  Ryy =    gupxy * gupxy * Axx + gupyy * gupyy * Ayy + gupyz * gupyz * Azz + &
-      TWO*(gupxy * gupyy * Axy + gupxy * gupyz * Axz + gupyy * gupyz * Ayz)
-
-  Rzz =    gupxz * gupxz * Axx + gupyz * gupyz * Ayy + gupzz * gupzz * Azz + &
-      TWO*(gupxz * gupyz * Axy + gupxz * gupzz * Axz + gupyz * gupzz * Ayz)
-
-  Rxy =    gupxx * gupxy * Axx + gupxy * gupyy * Ayy + gupxz * gupyz * Azz + &
-          (gupxx * gupyy       + gupxy * gupxy)* Axy                       + &
-          (gupxx * gupyz       + gupxz * gupxy)* Axz                       + &
-          (gupxy * gupyz       + gupxz * gupyy)* Ayz
-
-  Rxz =    gupxx * gupxz * Axx + gupxy * gupyz * Ayy + gupxz * gupzz * Azz + &
-          (gupxx * gupyz       + gupxy * gupxz)* Axy                       + &
-          (gupxx * gupzz       + gupxz * gupxz)* Axz                       + &
-          (gupxy * gupzz       + gupxz * gupyz)* Ayz
-
-  Ryz =    gupxy * gupxz * Axx + gupyy * gupyz * Ayy + gupyz * gupzz * Azz + &
-          (gupxy * gupyz       + gupyy * gupxz)* Axy                       + &
-          (gupxy * gupzz       + gupyz * gupxz)* Axz                       + &
-          (gupyy * gupzz       + gupyz * gupyz)* Ayz
-
-! Right hand side for Gam^i without shift terms...
-  call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
-  call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
-
-   Gamx_rhs = - TWO * (   Lapx * Rxx +   Lapy * Rxy +   Lapz * Rxz ) + &
-        TWO * alpn1 * (                                                &
-        -F3o2/chin1 * (   chix * Rxx +   chiy * Rxy +   chiz * Rxz ) - &
-              gupxx * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-              gupxy * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-              gupxz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-                        Gamxxx * Rxx + Gamxyy * Ryy + Gamxzz * Rzz   + &
-                TWO * ( Gamxxy * Rxy + Gamxxz * Rxz + Gamxyz * Ryz ) )
-
-   Gamy_rhs = - TWO * (   Lapx * Rxy +   Lapy * Ryy +   Lapz * Ryz ) + &
-        TWO * alpn1 * (                                                &
-        -F3o2/chin1 * (   chix * Rxy +  chiy * Ryy +    chiz * Ryz ) - &
-              gupxy * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-              gupyy * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-              gupyz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-                        Gamyxx * Rxx + Gamyyy * Ryy + Gamyzz * Rzz   + &
-                TWO * ( Gamyxy * Rxy + Gamyxz * Rxz + Gamyyz * Ryz ) )
-
-   Gamz_rhs = - TWO * (   Lapx * Rxz +   Lapy * Ryz +   Lapz * Rzz ) + &
-        TWO * alpn1 * (                                                &
-        -F3o2/chin1 * (   chix * Rxz +  chiy * Ryz +    chiz * Rzz ) - &
-              gupxz * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-              gupyz * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-              gupzz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-                        Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz   + &
-                TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) )
+  call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
+  call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev)
+  call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev)
+  call fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
+  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)
+
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+    divb_loc = betaxx(i,j,k) + betayy(i,j,k) + betazz(i,j,k)
+    div_beta(i,j,k) = divb_loc
+
+    chi_rhs(i,j,k) = F2o3 * chin1(i,j,k) * (alpn1(i,j,k) * trK(i,j,k) - divb_loc)
+
+    gxx_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axx(i,j,k) - F2o3 * gxx(i,j,k) * divb_loc + &
+         TWO * ( gxx(i,j,k) * betaxx(i,j,k) + gxy(i,j,k) * betayx(i,j,k) + gxz(i,j,k) * betazx(i,j,k) )
+
+    gyy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayy(i,j,k) - F2o3 * gyy(i,j,k) * divb_loc + &
+         TWO * ( gxy(i,j,k) * betaxy(i,j,k) + gyy(i,j,k) * betayy(i,j,k) + gyz(i,j,k) * betazy(i,j,k) )
+
+    gzz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Azz(i,j,k) - F2o3 * gzz(i,j,k) * divb_loc + &
+         TWO * ( gxz(i,j,k) * betaxz(i,j,k) + gyz(i,j,k) * betayz(i,j,k) + gzz(i,j,k) * betazz(i,j,k) )
+
+    gxy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axy(i,j,k) + F1o3 * gxy(i,j,k) * divb_loc + &
+         gxx(i,j,k) * betaxy(i,j,k) + gxz(i,j,k) * betazy(i,j,k) + gyy(i,j,k) * betayx(i,j,k) + &
+         gyz(i,j,k) * betazx(i,j,k) - gxy(i,j,k) * betazz(i,j,k)
+
+    gyz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayz(i,j,k) + F1o3 * gyz(i,j,k) * divb_loc + &
+         gxy(i,j,k) * betaxz(i,j,k) + gyy(i,j,k) * betayz(i,j,k) + gxz(i,j,k) * betaxy(i,j,k) + &
+         gzz(i,j,k) * betazy(i,j,k) - gyz(i,j,k) * betaxx(i,j,k)
+
+    gxz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axz(i,j,k) + F1o3 * gxz(i,j,k) * divb_loc + &
+         gxx(i,j,k) * betaxz(i,j,k) + gxy(i,j,k) * betayz(i,j,k) + gyz(i,j,k) * betayx(i,j,k) + &
+         gzz(i,j,k) * betazx(i,j,k) - gxz(i,j,k) * betayy(i,j,k)
+
+    det_loc = gxx(i,j,k) * gyy(i,j,k) * gzz(i,j,k) + gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) + &
+         gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k) * gxz(i,j,k) - &
+         gxy(i,j,k) * gxy(i,j,k) * gzz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) * gyz(i,j,k)
+    gupxx_loc = ( gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gyz(i,j,k) ) / det_loc
+    gupxy_loc = - ( gxy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gxz(i,j,k) ) / det_loc
+    gupxz_loc = ( gxy(i,j,k) * gyz(i,j,k) - gyy(i,j,k) * gxz(i,j,k) ) / det_loc
+    gupyy_loc = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k) * gxz(i,j,k) ) / det_loc
+    gupyz_loc = - ( gxx(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / det_loc
+    gupzz_loc = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k) * gxy(i,j,k) ) / det_loc
+    gupxx(i,j,k) = gupxx_loc
+    gupxy(i,j,k) = gupxy_loc
+    gupxz(i,j,k) = gupxz_loc
+    gupyy(i,j,k) = gupyy_loc
+    gupyz(i,j,k) = gupyz_loc
+    gupzz(i,j,k) = gupzz_loc
+
+    if(co == 0)then
+      Gmx_Res(i,j,k) = Gamx(i,j,k) - ( &
+           gupxx_loc*(gupxx_loc*gxxx(i,j,k)+gupxy_loc*gxyx(i,j,k)+gupxz_loc*gxzx(i,j,k)) + &
+           gupxy_loc*(gupxx_loc*gxyx(i,j,k)+gupxy_loc*gyyx(i,j,k)+gupxz_loc*gyzx(i,j,k)) + &
+           gupxz_loc*(gupxx_loc*gxzx(i,j,k)+gupxy_loc*gyzx(i,j,k)+gupxz_loc*gzzx(i,j,k)) + &
+           gupxx_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + &
+           gupxy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + &
+           gupxz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + &
+           gupxx_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
+           gupxy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
+           gupxz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
+      Gmy_Res(i,j,k) = Gamy(i,j,k) - ( &
+           gupxx_loc*(gupxy_loc*gxxx(i,j,k)+gupyy_loc*gxyx(i,j,k)+gupyz_loc*gxzx(i,j,k)) + &
+           gupxy_loc*(gupxy_loc*gxyx(i,j,k)+gupyy_loc*gyyx(i,j,k)+gupyz_loc*gyzx(i,j,k)) + &
+           gupxz_loc*(gupxy_loc*gxzx(i,j,k)+gupyy_loc*gyzx(i,j,k)+gupyz_loc*gzzx(i,j,k)) + &
+           gupxy_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + &
+           gupyy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + &
+           gupyz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + &
+           gupxy_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
+           gupyy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
+           gupyz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
+      Gmz_Res(i,j,k) = Gamz(i,j,k) - ( &
+           gupxx_loc*(gupxz_loc*gxxx(i,j,k)+gupyz_loc*gxyx(i,j,k)+gupzz_loc*gxzx(i,j,k)) + &
+           gupxy_loc*(gupxz_loc*gxyx(i,j,k)+gupyz_loc*gyyx(i,j,k)+gupzz_loc*gyzx(i,j,k)) + &
+           gupxz_loc*(gupxz_loc*gxzx(i,j,k)+gupyz_loc*gyzx(i,j,k)+gupzz_loc*gzzx(i,j,k)) + &
+           gupxy_loc*(gupxz_loc*gxxy(i,j,k)+gupyz_loc*gxyy(i,j,k)+gupzz_loc*gxzy(i,j,k)) + &
+           gupyy_loc*(gupxz_loc*gxyy(i,j,k)+gupyz_loc*gyyy(i,j,k)+gupzz_loc*gyzy(i,j,k)) + &
+           gupyz_loc*(gupxz_loc*gxzy(i,j,k)+gupyz_loc*gyzy(i,j,k)+gupzz_loc*gzzy(i,j,k)) + &
+           gupxz_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
+           gupyz_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
+           gupzz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
+    endif
+
+    Gamxxx(i,j,k)=HALF*( gupxx_loc*gxxx(i,j,k) + gupxy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupxz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+    Gamyxx(i,j,k)=HALF*( gupxy_loc*gxxx(i,j,k) + gupyy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupyz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+    Gamzxx(i,j,k)=HALF*( gupxz_loc*gxxx(i,j,k) + gupyz_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupzz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+
+    Gamxyy(i,j,k)=HALF*( gupxx_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupxy_loc*gyyy(i,j,k) + gupxz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+    Gamyyy(i,j,k)=HALF*( gupxy_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyy_loc*gyyy(i,j,k) + gupyz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+    Gamzyy(i,j,k)=HALF*( gupxz_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyz_loc*gyyy(i,j,k) + gupzz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+
+    Gamxzz(i,j,k)=HALF*( gupxx_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupxy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupxz_loc*gzzz(i,j,k))
+    Gamyzz(i,j,k)=HALF*( gupxy_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupyz_loc*gzzz(i,j,k))
+    Gamzzz(i,j,k)=HALF*( gupxz_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyz_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupzz_loc*gzzz(i,j,k))
+
+    Gamxxy(i,j,k)=HALF*( gupxx_loc*gxxy(i,j,k) + gupxy_loc*gyyx(i,j,k) + gupxz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+    Gamyxy(i,j,k)=HALF*( gupxy_loc*gxxy(i,j,k) + gupyy_loc*gyyx(i,j,k) + gupyz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+    Gamzxy(i,j,k)=HALF*( gupxz_loc*gxxy(i,j,k) + gupyz_loc*gyyx(i,j,k) + gupzz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+
+    Gamxxz(i,j,k)=HALF*( gupxx_loc*gxxz(i,j,k) + gupxy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupxz_loc*gzzx(i,j,k) )
+    Gamyxz(i,j,k)=HALF*( gupxy_loc*gxxz(i,j,k) + gupyy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupyz_loc*gzzx(i,j,k) )
+    Gamzxz(i,j,k)=HALF*( gupxz_loc*gxxz(i,j,k) + gupyz_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupzz_loc*gzzx(i,j,k) )
+
+    Gamxyz(i,j,k)=HALF*( gupxx_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupxy_loc*gyyz(i,j,k) + gupxz_loc*gzzy(i,j,k) )
+    Gamyyz(i,j,k)=HALF*( gupxy_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyy_loc*gyyz(i,j,k) + gupyz_loc*gzzy(i,j,k) )
+    Gamzyz(i,j,k)=HALF*( gupxz_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyz_loc*gyyz(i,j,k) + gupzz_loc*gzzy(i,j,k) )
+  enddo
+  enddo
+  enddo
+! Raise indices of \tilde A_{ij} and store in R_ij
+
+! Right hand side for Gam^i without shift terms...
+  call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
+  call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+    gupxx_loc = gupxx(i,j,k)
+    gupxy_loc = gupxy(i,j,k)
+    gupxz_loc = gupxz(i,j,k)
+    gupyy_loc = gupyy(i,j,k)
+    gupyz_loc = gupyz(i,j,k)
+    gupzz_loc = gupzz(i,j,k)
+
+    Rxx_loc = gupxx_loc * gupxx_loc * Axx(i,j,k) + gupxy_loc * gupxy_loc * Ayy(i,j,k) + gupxz_loc * gupxz_loc * Azz(i,j,k) + &
+         TWO * (gupxx_loc * gupxy_loc * Axy(i,j,k) + gupxx_loc * gupxz_loc * Axz(i,j,k) + gupxy_loc * gupxz_loc * Ayz(i,j,k))
+    Ryy_loc = gupxy_loc * gupxy_loc * Axx(i,j,k) + gupyy_loc * gupyy_loc * Ayy(i,j,k) + gupyz_loc * gupyz_loc * Azz(i,j,k) + &
+         TWO * (gupxy_loc * gupyy_loc * Axy(i,j,k) + gupxy_loc * gupyz_loc * Axz(i,j,k) + gupyy_loc * gupyz_loc * Ayz(i,j,k))
+    Rzz_loc = gupxz_loc * gupxz_loc * Axx(i,j,k) + gupyz_loc * gupyz_loc * Ayy(i,j,k) + gupzz_loc * gupzz_loc * Azz(i,j,k) + &
+         TWO * (gupxz_loc * gupyz_loc * Axy(i,j,k) + gupxz_loc * gupzz_loc * Axz(i,j,k) + gupyz_loc * gupzz_loc * Ayz(i,j,k))
+    Rxy_loc = gupxx_loc * gupxy_loc * Axx(i,j,k) + gupxy_loc * gupyy_loc * Ayy(i,j,k) + gupxz_loc * gupyz_loc * Azz(i,j,k) + &
+         (gupxx_loc * gupyy_loc + gupxy_loc * gupxy_loc) * Axy(i,j,k) + &
+         (gupxx_loc * gupyz_loc + gupxz_loc * gupxy_loc) * Axz(i,j,k) + &
+         (gupxy_loc * gupyz_loc + gupxz_loc * gupyy_loc) * Ayz(i,j,k)
+    Rxz_loc = gupxx_loc * gupxz_loc * Axx(i,j,k) + gupxy_loc * gupyz_loc * Ayy(i,j,k) + gupxz_loc * gupzz_loc * Azz(i,j,k) + &
+         (gupxx_loc * gupyz_loc + gupxy_loc * gupxz_loc) * Axy(i,j,k) + &
+         (gupxx_loc * gupzz_loc + gupxz_loc * gupxz_loc) * Axz(i,j,k) + &
+         (gupxy_loc * gupzz_loc + gupxz_loc * gupyz_loc) * Ayz(i,j,k)
+    Ryz_loc = gupxy_loc * gupxz_loc * Axx(i,j,k) + gupyy_loc * gupyz_loc * Ayy(i,j,k) + gupyz_loc * gupzz_loc * Azz(i,j,k) + &
+         (gupxy_loc * gupyz_loc + gupyy_loc * gupxz_loc) * Axy(i,j,k) + &
+         (gupxy_loc * gupzz_loc + gupyz_loc * gupxz_loc) * Axz(i,j,k) + &
+         (gupyy_loc * gupzz_loc + gupyz_loc * gupyz_loc) * Ayz(i,j,k)
+    Rxx(i,j,k) = Rxx_loc
+    Ryy(i,j,k) = Ryy_loc
+    Rzz(i,j,k) = Rzz_loc
+    Rxy(i,j,k) = Rxy_loc
+    Rxz(i,j,k) = Rxz_loc
+    Ryz(i,j,k) = Ryz_loc
+
+    Gamx_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxx_loc + Lapy(i,j,k) * Rxy_loc + Lapz(i,j,k) * Rxz_loc) + &
+         TWO * alpn1(i,j,k) * ( &
+         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxx_loc + chiy(i,j,k) * Rxy_loc + chiz(i,j,k) * Rxz_loc) - &
+         gupxx_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
+         gupxy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
+         gupxz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
+         Gamxxx(i,j,k) * Rxx_loc + Gamxyy(i,j,k) * Ryy_loc + Gamxzz(i,j,k) * Rzz_loc + &
+         TWO * (Gamxxy(i,j,k) * Rxy_loc + Gamxxz(i,j,k) * Rxz_loc + Gamxyz(i,j,k) * Ryz_loc))
+
+    Gamy_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxy_loc + Lapy(i,j,k) * Ryy_loc + Lapz(i,j,k) * Ryz_loc) + &
+         TWO * alpn1(i,j,k) * ( &
+         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxy_loc + chiy(i,j,k) * Ryy_loc + chiz(i,j,k) * Ryz_loc) - &
+         gupxy_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
+         gupyy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
+         gupyz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
+         Gamyxx(i,j,k) * Rxx_loc + Gamyyy(i,j,k) * Ryy_loc + Gamyzz(i,j,k) * Rzz_loc + &
+         TWO * (Gamyxy(i,j,k) * Rxy_loc + Gamyxz(i,j,k) * Rxz_loc + Gamyyz(i,j,k) * Ryz_loc))
+
+    Gamz_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxz_loc + Lapy(i,j,k) * Ryz_loc + Lapz(i,j,k) * Rzz_loc) + &
+         TWO * alpn1(i,j,k) * ( &
+         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxz_loc + chiy(i,j,k) * Ryz_loc + chiz(i,j,k) * Rzz_loc) - &
+         gupxz_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
+         gupyz_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
+         gupzz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
+         Gamzxx(i,j,k) * Rxx_loc + Gamzyy(i,j,k) * Ryy_loc + Gamzzz(i,j,k) * Rzz_loc + &
+         TWO * (Gamzxy(i,j,k) * Rxy_loc + Gamzxz(i,j,k) * Rxz_loc + Gamzyz(i,j,k) * Ryz_loc))
+  enddo
+  enddo
+  enddo

  call fdderivs(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,&
                X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev)
@@ -321,38 +359,54 @@
  call fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,&
                X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev)

-  fxx = gxxx + gxyy + gxzz
-  fxy = gxyx + gyyy + gyzz
-  fxz = gxzx + gyzy + gzzz
-
-  Gamxa =       gupxx * Gamxxx + gupyy * Gamxyy + gupzz * Gamxzz + &
-          TWO*( gupxy * Gamxxy + gupxz * Gamxxz + gupyz * Gamxyz )
-  Gamya =       gupxx * Gamyxx + gupyy * Gamyyy + gupzz * Gamyzz + &
-          TWO*( gupxy * Gamyxy + gupxz * Gamyxz + gupyz * Gamyyz )
-  Gamza =       gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz + &
-          TWO*( gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz )
-
-  call fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev)
-  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)
-
-  Gamx_rhs =               Gamx_rhs +  F2o3 *  Gamxa * div_beta        - &
-                     Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz  + &
-             F1o3 * (gupxx * fxx    + gupxy * fxy    + gupxz * fxz    ) + &
-                     gupxx * gxxx   + gupyy * gyyx   + gupzz * gzzx    + &
-              TWO * (gupxy * gxyx   + gupxz * gxzx   + gupyz * gyzx  )
-
-  Gamy_rhs =               Gamy_rhs +  F2o3 *  Gamya * div_beta        - &
-                     Gamxa * betayx - Gamya * betayy - Gamza * betayz  + &
-             F1o3 * (gupxy * fxx    + gupyy * fxy    + gupyz * fxz    ) + &
-                     gupxx * gxxy   + gupyy * gyyy   + gupzz * gzzy    + &
-              TWO * (gupxy * gxyy   + gupxz * gxzy   + gupyz * gyzy  )
-
-  Gamz_rhs =               Gamz_rhs +  F2o3 *  Gamza * div_beta        - &
-                     Gamxa * betazx - Gamya * betazy - Gamza * betazz  + &
-             F1o3 * (gupxz * fxx    + gupyz * fxy    + gupzz * fxz    ) + &
-                     gupxx * gxxz   + gupyy * gyyz   + gupzz * gzzz    + &
-              TWO * (gupxy * gxyz   + gupxz * gxzz   + gupyz * gyzz  )    !rhs for Gam^i
+  call fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev)
+  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)
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+    divb_loc = div_beta(i,j,k)
+    fxx_loc = gxxx(i,j,k) + gxyy(i,j,k) + gxzz(i,j,k)
+    fxy_loc = gxyx(i,j,k) + gyyy(i,j,k) + gyzz(i,j,k)
+    fxz_loc = gxzx(i,j,k) + gyzy(i,j,k) + gzzz(i,j,k)
+
+    gupxx_loc = gupxx(i,j,k)
+    gupxy_loc = gupxy(i,j,k)
+    gupxz_loc = gupxz(i,j,k)
+    gupyy_loc = gupyy(i,j,k)
+    gupyz_loc = gupyz(i,j,k)
+    gupzz_loc = gupzz(i,j,k)
+
+    Gamxa_loc = gupxx_loc * Gamxxx(i,j,k) + gupyy_loc * Gamxyy(i,j,k) + gupzz_loc * Gamxzz(i,j,k) + &
+         TWO * (gupxy_loc * Gamxxy(i,j,k) + gupxz_loc * Gamxxz(i,j,k) + gupyz_loc * Gamxyz(i,j,k))
+    Gamya_loc = gupxx_loc * Gamyxx(i,j,k) + gupyy_loc * Gamyyy(i,j,k) + gupzz_loc * Gamyzz(i,j,k) + &
+         TWO * (gupxy_loc * Gamyxy(i,j,k) + gupxz_loc * Gamyxz(i,j,k) + gupyz_loc * Gamyyz(i,j,k))
+    Gamza_loc = gupxx_loc * Gamzxx(i,j,k) + gupyy_loc * Gamzyy(i,j,k) + gupzz_loc * Gamzzz(i,j,k) + &
+         TWO * (gupxy_loc * Gamzxy(i,j,k) + gupxz_loc * Gamzxz(i,j,k) + gupyz_loc * Gamzyz(i,j,k))
+    Gamxa(i,j,k) = Gamxa_loc
+    Gamya(i,j,k) = Gamya_loc
+    Gamza(i,j,k) = Gamza_loc
+
+    Gamx_rhs(i,j,k) = Gamx_rhs(i,j,k) + F2o3 * Gamxa_loc * divb_loc - &
+         Gamxa_loc * betaxx(i,j,k) - Gamya_loc * betaxy(i,j,k) - Gamza_loc * betaxz(i,j,k) + &
+         F1o3 * (gupxx_loc * fxx_loc + gupxy_loc * fxy_loc + gupxz_loc * fxz_loc) + &
+         gupxx_loc * gxxx(i,j,k) + gupyy_loc * gyyx(i,j,k) + gupzz_loc * gzzx(i,j,k) + &
+         TWO * (gupxy_loc * gxyx(i,j,k) + gupxz_loc * gxzx(i,j,k) + gupyz_loc * gyzx(i,j,k))
+
+    Gamy_rhs(i,j,k) = Gamy_rhs(i,j,k) + F2o3 * Gamya_loc * divb_loc - &
+         Gamxa_loc * betayx(i,j,k) - Gamya_loc * betayy(i,j,k) - Gamza_loc * betayz(i,j,k) + &
+         F1o3 * (gupxy_loc * fxx_loc + gupyy_loc * fxy_loc + gupyz_loc * fxz_loc) + &
+         gupxx_loc * gxxy(i,j,k) + gupyy_loc * gyyy(i,j,k) + gupzz_loc * gzzy(i,j,k) + &
+         TWO * (gupxy_loc * gxyy(i,j,k) + gupxz_loc * gxzy(i,j,k) + gupyz_loc * gyzy(i,j,k))
+
+    Gamz_rhs(i,j,k) = Gamz_rhs(i,j,k) + F2o3 * Gamza_loc * divb_loc - &
+         Gamxa_loc * betazx(i,j,k) - Gamya_loc * betazy(i,j,k) - Gamza_loc * betazz(i,j,k) + &
+         F1o3 * (gupxz_loc * fxx_loc + gupyz_loc * fxy_loc + gupzz_loc * fxz_loc) + &
+         gupxx_loc * gxxz(i,j,k) + gupyy_loc * gyyz(i,j,k) + gupzz_loc * gzzz(i,j,k) + &
+         TWO * (gupxy_loc * gxyz(i,j,k) + gupxz_loc * gxzz(i,j,k) + gupyz_loc * gyzz(i,j,k))
+  enddo
+  enddo
+  enddo

 !first kind of connection stored in gij,k
  gxxx = gxx * Gamxxx + gxy * Gamyxx + gxz * Gamzxx
@@ -601,192 +655,190 @@
            Gamxyz * gxzz + Gamyyz * gyzz + Gamzyz * gzzz  + &
            Gamxzz * gxzy + Gamyzz * gyzy + Gamzzz * gzzy  + &
            Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz )
-!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)
-
-  fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz
-  fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz
-  fxz = fxz - Gamxxz * chix - Gamyxz * chiy - Gamzxz * chiz
-  fyy = fyy - Gamxyy * chix - Gamyyy * chiy - Gamzyy * chiz
-  fyz = fyz - Gamxyz * chix - Gamyyz * chiy - Gamzyz * chiz
-  fzz = fzz - Gamxzz * chix - Gamyzz * chiy - Gamzzz * chiz
-! Store D^l D_l chi - 3/(2*chi) D^l chi D_l chi in f
-
-  f =        gupxx * ( fxx - F3o2/chin1 * chix * chix ) + &
-             gupyy * ( fyy - F3o2/chin1 * chiy * chiy ) + &
-             gupzz * ( fzz - F3o2/chin1 * chiz * chiz ) + &
-       TWO * gupxy * ( fxy - F3o2/chin1 * chix * chiy ) + &
-       TWO * gupxz * ( fxz - F3o2/chin1 * chix * chiz ) + &
-       TWO * gupyz * ( fyz - F3o2/chin1 * chiy * chiz ) 
-! Add chi part to Ricci tensor:
-
-  Rxx = Rxx + (fxx - chix*chix/chin1/TWO + gxx * f)/chin1/TWO
-  Ryy = Ryy + (fyy - chiy*chiy/chin1/TWO + gyy * f)/chin1/TWO
-  Rzz = Rzz + (fzz - chiz*chiz/chin1/TWO + gzz * f)/chin1/TWO
-  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
-
-! 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)
-
-  gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1
-  gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1
-  gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1
-! now get physical second kind of connection
-  Gamxxx = Gamxxx - ( (chix + chix)/chin1 - gxx * gxxx )*HALF
-  Gamyxx = Gamyxx - (                     - gxx * gxxy )*HALF
-  Gamzxx = Gamzxx - (                     - gxx * gxxz )*HALF
-  Gamxyy = Gamxyy - (                     - gyy * gxxx )*HALF
-  Gamyyy = Gamyyy - ( (chiy + chiy)/chin1 - gyy * gxxy )*HALF
-  Gamzyy = Gamzyy - (                     - gyy * gxxz )*HALF
-  Gamxzz = Gamxzz - (                     - gzz * gxxx )*HALF
-  Gamyzz = Gamyzz - (                     - gzz * gxxy )*HALF
-  Gamzzz = Gamzzz - ( (chiz + chiz)/chin1 - gzz * gxxz )*HALF
-  Gamxxy = Gamxxy - (  chiy        /chin1 - gxy * gxxx )*HALF
-  Gamyxy = Gamyxy - (         chix /chin1 - gxy * gxxy )*HALF
-  Gamzxy = Gamzxy - (                     - gxy * gxxz )*HALF
-  Gamxxz = Gamxxz - (  chiz        /chin1 - gxz * gxxx )*HALF
-  Gamyxz = Gamyxz - (                     - gxz * gxxy )*HALF
-  Gamzxz = Gamzxz - (         chix /chin1 - gxz * gxxz )*HALF
-  Gamxyz = Gamxyz - (                     - gyz * gxxx )*HALF
-  Gamyyz = Gamyyz - (  chiz        /chin1 - gyz * gxxy )*HALF
-  Gamzyz = Gamzyz - (         chiy /chin1 - gyz * gxxz )*HALF
-
-  fxx = fxx - Gamxxx*Lapx - Gamyxx*Lapy - Gamzxx*Lapz
-  fyy = fyy - Gamxyy*Lapx - Gamyyy*Lapy - Gamzyy*Lapz
-  fzz = fzz - Gamxzz*Lapx - Gamyzz*Lapy - Gamzzz*Lapz
-  fxy = fxy - Gamxxy*Lapx - Gamyxy*Lapy - Gamzxy*Lapz
-  fxz = fxz - Gamxxz*Lapx - Gamyxz*Lapy - Gamzxz*Lapz
-  fyz = fyz - Gamxyz*Lapx - Gamyyz*Lapy - Gamzyz*Lapz
-
-! store D^i D_i Lap in trK_rhs upto chi
-  trK_rhs =    gupxx * fxx + gupyy * fyy + gupzz * fzz + &
-        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz )
-#if 1        
-!! follow bam code
-  S =  chin1 * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + &
-     TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) )
-  f = F2o3 * trK * trK -(&
-       gupxx * ( &
-       gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + &
-       TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz) ) + &
-       gupyy * ( &
-       gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + &
-       TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz) ) + &
-       gupzz * ( &
-       gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + &
-       TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz) ) + &
-       TWO * ( &
-       gupxy * ( &
-       gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + &
-       gupxy * (Axx * Ayy + Axy * Axy) + &
-       gupxz * (Axx * Ayz + Axz * Axy) + &
-       gupyz * (Axy * Ayz + Axz * Ayy) ) + &
-       gupxz * ( &
-       gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + &
-       gupxy * (Axx * Ayz + Axy * Axz) + &
-       gupxz * (Axx * Azz + Axz * Axz) + &
-       gupyz * (Axy * Azz + Axz * Ayz) ) + &
-       gupyz * ( &
-       gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + &
-       gupxy * (Axy * Ayz + Ayy * Axz) + &
-       gupxz * (Axy * Azz + Ayz * Axz) + &
-       gupyz * (Ayy * Azz + Ayz * Ayz) ) )) -1.6d1*PI*rho + EIGHT * PI * S
-  f = - F1o3 *(  gupxx * fxx + gupyy * fyy + gupzz * fzz + &
-        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + alpn1/chin1*f)
-  
-  fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx
-  fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy
-  fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz
-  fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy
-  fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz
-  fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz
-#else        
-! Add lapse and S_ij parts to Ricci tensor:
-
-  fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx
-  fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy
-  fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz
-  fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy
-  fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz
-  fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz
-
-! Compute trace-free part (note: chi^-1 and chi cancel!):
-
-  f = F1o3 *(  gupxx * fxx + gupyy * fyy + gupzz * fzz + &
-        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) )
-#endif
-
-  Axx_rhs = fxx - gxx * f
-  Ayy_rhs = fyy - gyy * f
-  Azz_rhs = fzz - gzz * f
-  Axy_rhs = fxy - gxy * f
-  Axz_rhs = fxz - gxz * f
-  Ayz_rhs = fyz - gyz * f
-
-! Now: store A_il A^l_j into fij:
-
-  fxx =       gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + &
-       TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz)
-  fyy =       gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + &
-       TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz)
-  fzz =       gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + &
-       TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz)
-  fxy =       gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + &
-              gupxy *(Axx * Ayy + Axy * Axy)                            + &
-              gupxz *(Axx * Ayz + Axz * Axy)                            + &
-              gupyz *(Axy * Ayz + Axz * Ayy)
-  fxz =       gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + &
-              gupxy *(Axx * Ayz + Axy * Axz)                            + &
-              gupxz *(Axx * Azz + Axz * Axz)                            + &
-              gupyz *(Axy * Azz + Axz * Ayz)
-  fyz =       gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + &
-              gupxy *(Axy * Ayz + Ayy * Axz)                            + &
-              gupxz *(Axy * Azz + Ayz * Axz)                            + &
-              gupyz *(Ayy * Azz + Ayz * Ayz)
-
-  f = chin1
-! store D^i D_i Lap in trK_rhs
-  trK_rhs = f*trK_rhs
-          
-  Axx_rhs =           f * Axx_rhs+ alpn1 * (trK * Axx - TWO * fxx)  + &
-           TWO * (  Axx * betaxx +   Axy * betayx +   Axz * betazx )- &
-             F2o3 * Axx * div_beta
-
-  Ayy_rhs =           f * Ayy_rhs+ alpn1 * (trK * Ayy - TWO * fyy)  + &
-           TWO * (  Axy * betaxy +   Ayy * betayy +   Ayz * betazy )- &
-             F2o3 * Ayy * div_beta
-
-  Azz_rhs =           f * Azz_rhs+ alpn1 * (trK * Azz - TWO * fzz)  + &
-           TWO * (  Axz * betaxz +   Ayz * betayz +   Azz * betazz )- &
-             F2o3 * Azz * div_beta
-
-  Axy_rhs =           f * Axy_rhs+ alpn1 *( trK * Axy  - TWO * fxy )+ &
-                    Axx * betaxy                  +   Axz * betazy  + &
-                                     Ayy * betayx +   Ayz * betazx  + &
-             F1o3 * Axy * div_beta                -   Axy * betazz
-
-  Ayz_rhs =           f * Ayz_rhs+ alpn1 *( trK * Ayz  - TWO * fyz )+ &
-                    Axy * betaxz +   Ayy * betayz                   + &
-                    Axz * betaxy                  +   Azz * betazy  + &
-             F1o3 * Ayz * div_beta                -   Ayz * betaxx
- 
-  Axz_rhs =           f * Axz_rhs+ alpn1 *( trK * Axz  - TWO * fxz )+ &
-                    Axx * betaxz +   Axy * betayz                   + &
-                                     Ayz * betayx +   Azz * betazx  + &
-             F1o3 * Axz * div_beta                -   Axz * betayy      !rhs for Aij
-
-! Compute trace of S_ij
-
-  S =  f * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + &
-     TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) )
-
-  trK_rhs = - trK_rhs + alpn1 *( F1o3 * trK * trK         + &
-                gupxx * fxx + gupyy * fyy + gupzz * fzz   + &
-        TWO * ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + &
-       FOUR * PI * ( rho + S ))                                !rhs for trK
+!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)
+
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+    fxx(i,j,k) = fxx(i,j,k) - Gamxxx(i,j,k) * chix(i,j,k) - Gamyxx(i,j,k) * chiy(i,j,k) - Gamzxx(i,j,k) * chiz(i,j,k)
+    fxy(i,j,k) = fxy(i,j,k) - Gamxxy(i,j,k) * chix(i,j,k) - Gamyxy(i,j,k) * chiy(i,j,k) - Gamzxy(i,j,k) * chiz(i,j,k)
+    fxz(i,j,k) = fxz(i,j,k) - Gamxxz(i,j,k) * chix(i,j,k) - Gamyxz(i,j,k) * chiy(i,j,k) - Gamzxz(i,j,k) * chiz(i,j,k)
+    fyy(i,j,k) = fyy(i,j,k) - Gamxyy(i,j,k) * chix(i,j,k) - Gamyyy(i,j,k) * chiy(i,j,k) - Gamzyy(i,j,k) * chiz(i,j,k)
+    fyz(i,j,k) = fyz(i,j,k) - Gamxyz(i,j,k) * chix(i,j,k) - Gamyyz(i,j,k) * chiy(i,j,k) - Gamzyz(i,j,k) * chiz(i,j,k)
+    fzz(i,j,k) = fzz(i,j,k) - Gamxzz(i,j,k) * chix(i,j,k) - Gamyzz(i,j,k) * chiy(i,j,k) - Gamzzz(i,j,k) * chiz(i,j,k)
+
+    chin_loc = chin1(i,j,k)
+    f_loc = gupxx(i,j,k) * (fxx(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chix(i,j,k)) + &
+            gupyy(i,j,k) * (fyy(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiy(i,j,k)) + &
+            gupzz(i,j,k) * (fzz(i,j,k) - F3o2/chin_loc * chiz(i,j,k) * chiz(i,j,k)) + &
+            TWO * gupxy(i,j,k) * (fxy(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiy(i,j,k)) + &
+            TWO * gupxz(i,j,k) * (fxz(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiz(i,j,k)) + &
+            TWO * gupyz(i,j,k) * (fyz(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiz(i,j,k))
+    f(i,j,k) = f_loc
+
+    Rxx(i,j,k) = Rxx(i,j,k) + (fxx(i,j,k) - chix(i,j,k)*chix(i,j,k)/chin_loc/TWO + gxx(i,j,k) * f_loc)/chin_loc/TWO
+    Ryy(i,j,k) = Ryy(i,j,k) + (fyy(i,j,k) - chiy(i,j,k)*chiy(i,j,k)/chin_loc/TWO + gyy(i,j,k) * f_loc)/chin_loc/TWO
+    Rzz(i,j,k) = Rzz(i,j,k) + (fzz(i,j,k) - chiz(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gzz(i,j,k) * f_loc)/chin_loc/TWO
+    Rxy(i,j,k) = Rxy(i,j,k) + (fxy(i,j,k) - chix(i,j,k)*chiy(i,j,k)/chin_loc/TWO + gxy(i,j,k) * f_loc)/chin_loc/TWO
+    Rxz(i,j,k) = Rxz(i,j,k) + (fxz(i,j,k) - chix(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gxz(i,j,k) * f_loc)/chin_loc/TWO
+    Ryz(i,j,k) = Ryz(i,j,k) + (fyz(i,j,k) - chiy(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gyz(i,j,k) * f_loc)/chin_loc/TWO
+  enddo
+  enddo
+  enddo
+
+! 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)
+
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+    chin_loc = chin1(i,j,k)
+    gxxx(i,j,k) = (gupxx(i,j,k) * chix(i,j,k) + gupxy(i,j,k) * chiy(i,j,k) + gupxz(i,j,k) * chiz(i,j,k)) / chin_loc
+    gxxy(i,j,k) = (gupxy(i,j,k) * chix(i,j,k) + gupyy(i,j,k) * chiy(i,j,k) + gupyz(i,j,k) * chiz(i,j,k)) / chin_loc
+    gxxz(i,j,k) = (gupxz(i,j,k) * chix(i,j,k) + gupyz(i,j,k) * chiy(i,j,k) + gupzz(i,j,k) * chiz(i,j,k)) / chin_loc
+
+    Gamxxx(i,j,k) = Gamxxx(i,j,k) - ( (chix(i,j,k) + chix(i,j,k))/chin_loc - gxx(i,j,k) * gxxx(i,j,k) )*HALF
+    Gamyxx(i,j,k) = Gamyxx(i,j,k) - (                                   - gxx(i,j,k) * gxxy(i,j,k) )*HALF
+    Gamzxx(i,j,k) = Gamzxx(i,j,k) - (                                   - gxx(i,j,k) * gxxz(i,j,k) )*HALF
+    Gamxyy(i,j,k) = Gamxyy(i,j,k) - (                                   - gyy(i,j,k) * gxxx(i,j,k) )*HALF
+    Gamyyy(i,j,k) = Gamyyy(i,j,k) - ( (chiy(i,j,k) + chiy(i,j,k))/chin_loc - gyy(i,j,k) * gxxy(i,j,k) )*HALF
+    Gamzyy(i,j,k) = Gamzyy(i,j,k) - (                                   - gyy(i,j,k) * gxxz(i,j,k) )*HALF
+    Gamxzz(i,j,k) = Gamxzz(i,j,k) - (                                   - gzz(i,j,k) * gxxx(i,j,k) )*HALF
+    Gamyzz(i,j,k) = Gamyzz(i,j,k) - (                                   - gzz(i,j,k) * gxxy(i,j,k) )*HALF
+    Gamzzz(i,j,k) = Gamzzz(i,j,k) - ( (chiz(i,j,k) + chiz(i,j,k))/chin_loc - gzz(i,j,k) * gxxz(i,j,k) )*HALF
+    Gamxxy(i,j,k) = Gamxxy(i,j,k) - ( chiy(i,j,k) /chin_loc - gxy(i,j,k) * gxxx(i,j,k) )*HALF
+    Gamyxy(i,j,k) = Gamyxy(i,j,k) - ( chix(i,j,k) /chin_loc - gxy(i,j,k) * gxxy(i,j,k) )*HALF
+    Gamzxy(i,j,k) = Gamzxy(i,j,k) - (                     - gxy(i,j,k) * gxxz(i,j,k) )*HALF
+    Gamxxz(i,j,k) = Gamxxz(i,j,k) - ( chiz(i,j,k) /chin_loc - gxz(i,j,k) * gxxx(i,j,k) )*HALF
+    Gamyxz(i,j,k) = Gamyxz(i,j,k) - (                     - gxz(i,j,k) * gxxy(i,j,k) )*HALF
+    Gamzxz(i,j,k) = Gamzxz(i,j,k) - ( chix(i,j,k) /chin_loc - gxz(i,j,k) * gxxz(i,j,k) )*HALF
+    Gamxyz(i,j,k) = Gamxyz(i,j,k) - (                     - gyz(i,j,k) * gxxx(i,j,k) )*HALF
+    Gamyyz(i,j,k) = Gamyyz(i,j,k) - ( chiz(i,j,k) /chin_loc - gyz(i,j,k) * gxxy(i,j,k) )*HALF
+    Gamzyz(i,j,k) = Gamzyz(i,j,k) - ( chiy(i,j,k) /chin_loc - gyz(i,j,k) * gxxz(i,j,k) )*HALF
+
+    fxx(i,j,k) = fxx(i,j,k) - Gamxxx(i,j,k)*Lapx(i,j,k) - Gamyxx(i,j,k)*Lapy(i,j,k) - Gamzxx(i,j,k)*Lapz(i,j,k)
+    fyy(i,j,k) = fyy(i,j,k) - Gamxyy(i,j,k)*Lapx(i,j,k) - Gamyyy(i,j,k)*Lapy(i,j,k) - Gamzyy(i,j,k)*Lapz(i,j,k)
+    fzz(i,j,k) = fzz(i,j,k) - Gamxzz(i,j,k)*Lapx(i,j,k) - Gamyzz(i,j,k)*Lapy(i,j,k) - Gamzzz(i,j,k)*Lapz(i,j,k)
+    fxy(i,j,k) = fxy(i,j,k) - Gamxxy(i,j,k)*Lapx(i,j,k) - Gamyxy(i,j,k)*Lapy(i,j,k) - Gamzxy(i,j,k)*Lapz(i,j,k)
+    fxz(i,j,k) = fxz(i,j,k) - Gamxxz(i,j,k)*Lapx(i,j,k) - Gamyxz(i,j,k)*Lapy(i,j,k) - Gamzxz(i,j,k)*Lapz(i,j,k)
+    fyz(i,j,k) = fyz(i,j,k) - Gamxyz(i,j,k)*Lapx(i,j,k) - Gamyyz(i,j,k)*Lapy(i,j,k) - Gamzyz(i,j,k)*Lapz(i,j,k)
+
+    trK_rhs(i,j,k) = gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + &
+                     TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k))
+  enddo
+  enddo
+  enddo
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+    divb_loc = div_beta(i,j,k)
+    chin_loc = chin1(i,j,k)
+
+    S_loc = chin_loc * ( gupxx(i,j,k) * Sxx(i,j,k) + gupyy(i,j,k) * Syy(i,j,k) + gupzz(i,j,k) * Szz(i,j,k) + &
+           TWO * (gupxy(i,j,k) * Sxy(i,j,k) + gupxz(i,j,k) * Sxz(i,j,k) + gupyz(i,j,k) * Syz(i,j,k)) )
+    S(i,j,k) = S_loc
+
+    f_loc = F2o3 * trK(i,j,k) * trK(i,j,k) - ( &
+            gupxx(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axx(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + &
+                             gupzz(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + &
+                             TWO * (gupxy(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupxz(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + &
+                                    gupyz(i,j,k) * Axy(i,j,k) * Axz(i,j,k)) ) + &
+            gupyy(i,j,k) * ( gupxx(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayy(i,j,k) + &
+                             gupzz(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + &
+                             TWO * (gupxy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + gupxz(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + &
+                                    gupyz(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k)) ) + &
+            gupzz(i,j,k) * ( gupxx(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + &
+                             gupzz(i,j,k) * Azz(i,j,k) * Azz(i,j,k) + &
+                             TWO * (gupxy(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + gupxz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + &
+                                    gupyz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k)) ) + &
+            TWO * ( gupxy(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + &
+                                     gupzz(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + &
+                                     gupxy(i,j,k) * (Axx(i,j,k) * Ayy(i,j,k) + Axy(i,j,k) * Axy(i,j,k)) + &
+                                     gupxz(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Axy(i,j,k)) + &
+                                     gupyz(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Ayy(i,j,k)) ) + &
+                    gupxz(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + &
+                                     gupzz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + &
+                                     gupxy(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axy(i,j,k) * Axz(i,j,k)) + &
+                                     gupxz(i,j,k) * (Axx(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Axz(i,j,k)) + &
+                                     gupyz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Ayz(i,j,k)) ) + &
+                    gupyz(i,j,k) * ( gupxx(i,j,k) * Axy(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k) + &
+                                     gupzz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k) + &
+                                     gupxy(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Ayy(i,j,k) * Axz(i,j,k)) + &
+                                     gupxz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Axz(i,j,k)) + &
+                                     gupyz(i,j,k) * (Ayy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Ayz(i,j,k)) ) ) ) - &
+            F16 * PI * rho(i,j,k) + EIGHT * PI * S_loc
+
+    f_loc = -F1o3 * ( gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + &
+            TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k)) + &
+            alpn1(i,j,k)/chin_loc * f_loc )
+    f(i,j,k) = f_loc
+
+    l_fxx = alpn1(i,j,k) * (Rxx(i,j,k) - EIGHT * PI * Sxx(i,j,k)) - fxx(i,j,k)
+    l_fxy = alpn1(i,j,k) * (Rxy(i,j,k) - EIGHT * PI * Sxy(i,j,k)) - fxy(i,j,k)
+    l_fxz = alpn1(i,j,k) * (Rxz(i,j,k) - EIGHT * PI * Sxz(i,j,k)) - fxz(i,j,k)
+    l_fyy = alpn1(i,j,k) * (Ryy(i,j,k) - EIGHT * PI * Syy(i,j,k)) - fyy(i,j,k)
+    l_fyz = alpn1(i,j,k) * (Ryz(i,j,k) - EIGHT * PI * Syz(i,j,k)) - fyz(i,j,k)
+    l_fzz = alpn1(i,j,k) * (Rzz(i,j,k) - EIGHT * PI * Szz(i,j,k)) - fzz(i,j,k)
+
+    Axx_rhs(i,j,k) = l_fxx - gxx(i,j,k) * f_loc
+    Ayy_rhs(i,j,k) = l_fyy - gyy(i,j,k) * f_loc
+    Azz_rhs(i,j,k) = l_fzz - gzz(i,j,k) * f_loc
+    Axy_rhs(i,j,k) = l_fxy - gxy(i,j,k) * f_loc
+    Axz_rhs(i,j,k) = l_fxz - gxz(i,j,k) * f_loc
+    Ayz_rhs(i,j,k) = l_fyz - gyz(i,j,k) * f_loc
+
+    fxx(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axx(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + &
+                 gupzz(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + TWO * (gupxy(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + &
+                 gupxz(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyz(i,j,k) * Axy(i,j,k) * Axz(i,j,k))
+    fyy(i,j,k) = gupxx(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayy(i,j,k) + &
+                 gupzz(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + TWO * (gupxy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + &
+                 gupxz(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + gupyz(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k))
+    fzz(i,j,k) = gupxx(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + &
+                 gupzz(i,j,k) * Azz(i,j,k) * Azz(i,j,k) + TWO * (gupxy(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + &
+                 gupxz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + gupyz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k))
+    fxy(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + &
+                 gupzz(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + gupxy(i,j,k) * (Axx(i,j,k) * Ayy(i,j,k) + Axy(i,j,k) * Axy(i,j,k)) + &
+                 gupxz(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Axy(i,j,k)) + &
+                 gupyz(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Ayy(i,j,k))
+    fxz(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + &
+                 gupzz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + gupxy(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axy(i,j,k) * Axz(i,j,k)) + &
+                 gupxz(i,j,k) * (Axx(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Axz(i,j,k)) + &
+                 gupyz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Ayz(i,j,k))
+    fyz(i,j,k) = gupxx(i,j,k) * Axy(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k) + &
+                 gupzz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k) + gupxy(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Ayy(i,j,k) * Axz(i,j,k)) + &
+                 gupxz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Axz(i,j,k)) + &
+                 gupyz(i,j,k) * (Ayy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Ayz(i,j,k))
+
+    trK_rhs(i,j,k) = chin_loc * trK_rhs(i,j,k)
+
+    Axx_rhs(i,j,k) = chin_loc * Axx_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axx(i,j,k) - TWO * fxx(i,j,k)) + &
+                     TWO * (Axx(i,j,k) * betaxx(i,j,k) + Axy(i,j,k) * betayx(i,j,k) + Axz(i,j,k) * betazx(i,j,k)) - &
+                     F2o3 * Axx(i,j,k) * divb_loc
+    Ayy_rhs(i,j,k) = chin_loc * Ayy_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Ayy(i,j,k) - TWO * fyy(i,j,k)) + &
+                     TWO * (Axy(i,j,k) * betaxy(i,j,k) + Ayy(i,j,k) * betayy(i,j,k) + Ayz(i,j,k) * betazy(i,j,k)) - &
+                     F2o3 * Ayy(i,j,k) * divb_loc
+    Azz_rhs(i,j,k) = chin_loc * Azz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Azz(i,j,k) - TWO * fzz(i,j,k)) + &
+                     TWO * (Axz(i,j,k) * betaxz(i,j,k) + Ayz(i,j,k) * betayz(i,j,k) + Azz(i,j,k) * betazz(i,j,k)) - &
+                     F2o3 * Azz(i,j,k) * divb_loc
+    Axy_rhs(i,j,k) = chin_loc * Axy_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axy(i,j,k) - TWO * fxy(i,j,k)) + &
+                     Axx(i,j,k) * betaxy(i,j,k) + Axz(i,j,k) * betazy(i,j,k) + Ayy(i,j,k) * betayx(i,j,k) + &
+                     Ayz(i,j,k) * betazx(i,j,k) + F1o3 * Axy(i,j,k) * divb_loc - Axy(i,j,k) * betazz(i,j,k)
+    Ayz_rhs(i,j,k) = chin_loc * Ayz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Ayz(i,j,k) - TWO * fyz(i,j,k)) + &
+                     Axy(i,j,k) * betaxz(i,j,k) + Ayy(i,j,k) * betayz(i,j,k) + Axz(i,j,k) * betaxy(i,j,k) + &
+                     Azz(i,j,k) * betazy(i,j,k) + F1o3 * Ayz(i,j,k) * divb_loc - Ayz(i,j,k) * betaxx(i,j,k)
+    Axz_rhs(i,j,k) = chin_loc * Axz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axz(i,j,k) - TWO * fxz(i,j,k)) + &
+                     Axx(i,j,k) * betaxz(i,j,k) + Axy(i,j,k) * betayz(i,j,k) + Ayz(i,j,k) * betayx(i,j,k) + &
+                     Azz(i,j,k) * betazx(i,j,k) + F1o3 * Axz(i,j,k) * divb_loc - Axz(i,j,k) * betayy(i,j,k)
+
+    trK_rhs(i,j,k) = - trK_rhs(i,j,k) + alpn1(i,j,k) * ( F1o3 * trK(i,j,k) * trK(i,j,k) + &
+                    gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + &
+                    TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k)) + &
+                    FOUR * PI * (rho(i,j,k) + S_loc) )
+  enddo
+  enddo
+  enddo
  
 !!!! gauge variable part

@@ -948,15 +1000,15 @@
 !!!!!!!!!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.
-
-  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)
+! gxx/gyy/gzz (=dxx/dyy/dzz+1): stencil coefficients sum to zero,
+! so the constant offset has no effect on dissipation.
+
+  call lopsided_kodis(ex,X,Y,Z,dxx,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,dyy,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,dzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps)

  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)
--- a/AMSS_NCKU_source/bssn_rhs_c.C
+++ b/AMSS_NCKU_source/bssn_rhs_c.C
@@ -1139,59 +1139,59 @@ int f_compute_rhs_bssn(int *ex, double &T,
            fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0);
            fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0);
            fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0);
-        }
-        // 7ms //
-        for (int i=0;i<all;i+=1) {
-            gxxx[i] = gxxx[i] - (  Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]
-                                  + Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]) - chix[i]*Axx[i]/chin1[i];
-            gxyx[i] = gxyx[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
-                                  + Gamxxx[i] * Axy[i] + Gamyxx[i] * Ayy[i] + Gamzxx[i] * Ayz[i]) - chix[i]*Axy[i]/chin1[i];
-            gxzx[i] = gxzx[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
-                                  + Gamxxx[i] * Axz[i] + Gamyxx[i] * Ayz[i] + Gamzxx[i] * Azz[i]) - chix[i]*Axz[i]/chin1[i];
-            gyyx[i] = gyyx[i] - (  Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]
-                                  + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chix[i]*Ayy[i]/chin1[i];
-            gyzx[i] = gyzx[i] - (  Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]
-                                  + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chix[i]*Ayz[i]/chin1[i];
-            gzzx[i] = gzzx[i] - (  Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]
-                                  + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chix[i]*Azz[i]/chin1[i];
-            gxxy[i] = gxxy[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
-                                  + Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]) - chiy[i]*Axx[i]/chin1[i];
-            gxyy[i] = gxyy[i] - (  Gamxyy[i] * Axx[i] + Gamyyy[i] * Axy[i] + Gamzyy[i] * Axz[i]
-                                  + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chiy[i]*Axy[i]/chin1[i];
-            gxzy[i] = gxzy[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
-                                  + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chiy[i]*Axz[i]/chin1[i];
-            gyyy[i] = gyyy[i] - (  Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]
-                                  + Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]) - chiy[i]*Ayy[i]/chin1[i];
-            gyzy[i] = gyzy[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
-                                  + Gamxyy[i] * Axz[i] + Gamyyy[i] * Ayz[i] + Gamzyy[i] * Azz[i]) - chiy[i]*Ayz[i]/chin1[i];
-            gzzy[i] = gzzy[i] - (  Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]
-                                  + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiy[i]*Azz[i]/chin1[i];
-            gxxz[i] = gxxz[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
-                                  + Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]) - chiz[i]*Axx[i]/chin1[i];
-            gxyz[i] = gxyz[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
-                                  + Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]) - chiz[i]*Axy[i]/chin1[i];
-            gxzz[i] = gxzz[i] - (  Gamxzz[i] * Axx[i] + Gamyzz[i] * Axy[i] + Gamzzz[i] * Axz[i]
-                                  + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chiz[i]*Axz[i]/chin1[i];
-            gyyz[i] = gyyz[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
-                                  + Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]) - chiz[i]*Ayy[i]/chin1[i];
-            gyzz[i] = gyzz[i] - (  Gamxzz[i] * Axy[i] + Gamyzz[i] * Ayy[i] + Gamzzz[i] * Ayz[i]
-                                  + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiz[i]*Ayz[i]/chin1[i];
-            gzzz[i] = gzzz[i] - (  Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]
-                                  + Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]) - chiz[i]*Azz[i]/chin1[i];
+            // 7ms //
+            for (int i=0;i<all;i+=1) {
+                gxxx[i] = gxxx[i] - (  Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]
+                                      + Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]) - chix[i]*Axx[i]/chin1[i];
+                gxyx[i] = gxyx[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
+                                      + Gamxxx[i] * Axy[i] + Gamyxx[i] * Ayy[i] + Gamzxx[i] * Ayz[i]) - chix[i]*Axy[i]/chin1[i];
+                gxzx[i] = gxzx[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
+                                      + Gamxxx[i] * Axz[i] + Gamyxx[i] * Ayz[i] + Gamzxx[i] * Azz[i]) - chix[i]*Axz[i]/chin1[i];
+                gyyx[i] = gyyx[i] - (  Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]
+                                      + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chix[i]*Ayy[i]/chin1[i];
+                gyzx[i] = gyzx[i] - (  Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]
+                                      + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chix[i]*Ayz[i]/chin1[i];
+                gzzx[i] = gzzx[i] - (  Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]
+                                      + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chix[i]*Azz[i]/chin1[i];
+                gxxy[i] = gxxy[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
+                                      + Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]) - chiy[i]*Axx[i]/chin1[i];
+                gxyy[i] = gxyy[i] - (  Gamxyy[i] * Axx[i] + Gamyyy[i] * Axy[i] + Gamzyy[i] * Axz[i]
+                                      + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chiy[i]*Axy[i]/chin1[i];
+                gxzy[i] = gxzy[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
+                                      + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chiy[i]*Axz[i]/chin1[i];
+                gyyy[i] = gyyy[i] - (  Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]
+                                      + Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]) - chiy[i]*Ayy[i]/chin1[i];
+                gyzy[i] = gyzy[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
+                                      + Gamxyy[i] * Axz[i] + Gamyyy[i] * Ayz[i] + Gamzyy[i] * Azz[i]) - chiy[i]*Ayz[i]/chin1[i];
+                gzzy[i] = gzzy[i] - (  Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]
+                                      + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiy[i]*Azz[i]/chin1[i];
+                gxxz[i] = gxxz[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
+                                      + Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]) - chiz[i]*Axx[i]/chin1[i];
+                gxyz[i] = gxyz[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
+                                      + Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]) - chiz[i]*Axy[i]/chin1[i];
+                gxzz[i] = gxzz[i] - (  Gamxzz[i] * Axx[i] + Gamyzz[i] * Axy[i] + Gamzzz[i] * Axz[i]
+                                      + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chiz[i]*Axz[i]/chin1[i];
+                gyyz[i] = gyyz[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
+                                      + Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]) - chiz[i]*Ayy[i]/chin1[i];
+                gyzz[i] = gyzz[i] - (  Gamxzz[i] * Axy[i] + Gamyzz[i] * Ayy[i] + Gamzzz[i] * Ayz[i]
+                                      + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiz[i]*Ayz[i]/chin1[i];
+                gzzz[i] = gzzz[i] - (  Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]
+                                      + Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]) - chiz[i]*Azz[i]/chin1[i];

-            movx_Res[i] = gupxx[i]*gxxx[i] + gupyy[i]*gxyy[i] + gupzz[i]*gxzz[i]
-                        + gupxy[i]*gxyx[i] + gupxz[i]*gxzx[i] + gupyz[i]*gxzy[i]
-                        + gupxy[i]*gxxy[i] + gupxz[i]*gxxz[i] + gupyz[i]*gxyz[i];
-            movy_Res[i] = gupxx[i]*gxyx[i] + gupyy[i]*gyyy[i] + gupzz[i]*gyzz[i]
-                        + gupxy[i]*gyyx[i] + gupxz[i]*gyzx[i] + gupyz[i]*gyzy[i]
-                        + gupxy[i]*gxyy[i] + gupxz[i]*gxyz[i] + gupyz[i]*gyyz[i];
-            movz_Res[i] = gupxx[i]*gxzx[i] + gupyy[i]*gyzy[i] + gupzz[i]*gzzz[i]
-                        + gupxy[i]*gyzx[i] + gupxz[i]*gzzx[i] + gupyz[i]*gzzy[i]
-                        + gupxy[i]*gxzy[i] + gupxz[i]*gxzz[i] + gupyz[i]*gyzz[i];
+                movx_Res[i] = gupxx[i]*gxxx[i] + gupyy[i]*gxyy[i] + gupzz[i]*gxzz[i]
+                            + gupxy[i]*gxyx[i] + gupxz[i]*gxzx[i] + gupyz[i]*gxzy[i]
+                            + gupxy[i]*gxxy[i] + gupxz[i]*gxxz[i] + gupyz[i]*gxyz[i];
+                movy_Res[i] = gupxx[i]*gxyx[i] + gupyy[i]*gyyy[i] + gupzz[i]*gyzz[i]
+                            + gupxy[i]*gyyx[i] + gupxz[i]*gyzx[i] + gupyz[i]*gyzy[i]
+                            + gupxy[i]*gxyy[i] + gupxz[i]*gxyz[i] + gupyz[i]*gyyz[i];
+                movz_Res[i] = gupxx[i]*gxzx[i] + gupyy[i]*gyzy[i] + gupzz[i]*gzzz[i]
+                            + gupxy[i]*gyzx[i] + gupxz[i]*gzzx[i] + gupyz[i]*gzzy[i]
+                            + gupxy[i]*gxzy[i] + gupxz[i]*gxzz[i] + gupyz[i]*gyzz[i];

-            movx_Res[i] = movx_Res[i] - F2o3*Kx[i] - F8*PI*Sx[i];
-            movy_Res[i] = movy_Res[i] - F2o3*Ky[i] - F8*PI*Sy[i];
-            movz_Res[i] = movz_Res[i] - F2o3*Kz[i] - F8*PI*Sz[i];
+                movx_Res[i] = movx_Res[i] - F2o3*Kx[i] - F8*PI*Sx[i];
+                movy_Res[i] = movy_Res[i] - F2o3*Ky[i] - F8*PI*Sy[i];
+                movz_Res[i] = movz_Res[i] - F2o3*Kz[i] - F8*PI*Sz[i];
+            }
        }


--- a/AMSS_NCKU_source/makefile
+++ b/AMSS_NCKU_source/makefile
@@ -30,8 +30,8 @@ CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma -ipo \
 f90appflags = -O3 -xHost -fp-model fast=2 -fma -ipo \
              -align array64byte -fpp -I${MKLROOT}/include $(POLINT6_FLAG)
 endif
-
-.SUFFIXES: .o .f90 .C .for .cu
+
+.SUFFIXES: .o .f90 .C .for .cu

 .f90.o:
 	$(f90) $(f90appflags) -c $< -o $@
@@ -64,8 +64,8 @@ lopsided_c.o: lopsided_c.C
 lopsided_kodis_c.o: lopsided_kodis_c.C
 	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@

-interp_lb_profile.o: interp_lb_profile.C interp_lb_profile.h
-	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+#interp_lb_profile.o: interp_lb_profile.C interp_lb_profile.h
+#	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@

 ## TwoPunctureABE uses fixed optimal flags with its own PGO profile, independent of CXXAPPFLAGS
 TP_PROFDATA = /home/$(shell whoami)/AMSS-NCKU/pgo_profile/TwoPunctureABE.profdata
@@ -105,13 +105,12 @@ C++FILES = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o\
           Parallel_bam.o scalar_class.o transpbh.o NullShellPatch2.o\
 	   NullShellPatch2_Evo.o writefile_f.o interp_lb_profile.o
 	   
-C++FILES_GPU = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o\
-           cgh.o surface_integral.o ShellPatch.o\
-	   bssnEScalar_class.o perf.o Z4c_class.o NullShellPatch.o\
-	   bssnEM_class.o cpbc_util.o z4c_rhs_point.o checkpoint.o\
-           Parallel_bam.o scalar_class.o transpbh.o NullShellPatch2.o\
-	   NullShellPatch2_Evo.o \
-	   bssn_gpu_class.o bssn_step_gpu.o bssn_macro.o writefile_f.o
+C++FILES_GPU = 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\
+	   bssnEScalar_class.o perf.o Z4c_class.o NullShellPatch.o\
+	   bssnEM_class.o cpbc_util.o z4c_rhs_point.o checkpoint.o\
+           Parallel_bam.o scalar_class.o transpbh.o NullShellPatch2.o\
+	   NullShellPatch2_Evo.o bssn_cuda_step.o writefile_f.o

 F90FILES_BASE = enforce_algebra.o fmisc.o initial_puncture.o prolongrestrict.o\
 	   prolongrestrict_cell.o prolongrestrict_vertex.o\
@@ -143,7 +142,7 @@ initial_guess.o Newton.o Jacobian.o ilucg.o IntPnts0.o IntPnts.o

 TwoPunctureFILES = TwoPunctureABE.o TwoPunctures.o 

-CUDAFILES = bssn_gpu.o bssn_gpu_rhs_ss.o
+CUDAFILES = bssn_gpu.o bssn_cuda_ops.o

 # file dependences
 $(C++FILES) $(C++FILES_GPU) $(F90FILES) $(CFILES) $(AHFDOBJS) $(CUDAFILES): macrodef.fh
--- a/AMSS_NCKU_source/makefile.inc
+++ b/AMSS_NCKU_source/makefile.inc
@@ -9,6 +9,7 @@ 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 -liomp5
+CUDA_LDLIBS = -L/usr/local/cuda-12.9/targets/x86_64-linux/lib -lcudart

 ## Memory allocator switch
 ##   1 (default) : link Intel oneTBB allocator (libtbbmalloc)
@@ -24,6 +25,8 @@ ifeq ($(USE_TBBMALLOC),1)
 LDLIBS := $(TBBMALLOC_LIBS) $(LDLIBS)
 endif

+LDLIBS := $(CUDA_LDLIBS) $(LDLIBS)
+
 ## PGO build mode switch (ABE only; TwoPunctureABE always uses opt flags)
 ##   opt        : (default) maximum performance with PGO profile-guided optimization
 ##   instrument : PGO Phase 1 instrumentation to collect fresh profile data
--- a/AMSS_NCKU_source/surface_integral.C
+++ b/AMSS_NCKU_source/surface_integral.C
@@ -180,19 +180,64 @@ surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry)
 //|============================================================================
 //| Destructor
 //|============================================================================
-surface_integral::~surface_integral()
-{
-  delete[] nx_g;
-  delete[] ny_g;
-  delete[] nz_g;
-  delete[] arcostheta;
-#ifdef GaussInt
-  delete[] wtcostheta;
-#endif
-}
-//|----------------------------------------------------------------
-//  spin weighted spinw component of psi4, general routine
-//  l takes from spinw to maxl; m takes from -l to l
+surface_integral::~surface_integral()
+{
+  release_cached_buffers();
+  delete[] nx_g;
+  delete[] ny_g;
+  delete[] nz_g;
+  delete[] arcostheta;
+#ifdef GaussInt
+  delete[] wtcostheta;
+#endif
+}
+
+void surface_integral::get_surface_points(double rex, double **pox)
+{
+  SpherePointCache &cache = sphere_point_cache[rex];
+  if (!cache.pox[0])
+  {
+    for (int i = 0; i < 3; ++i)
+      cache.pox[i] = new double[n_tot];
+    for (int n = 0; n < n_tot; ++n)
+    {
+      cache.pox[0][n] = rex * nx_g[n];
+      cache.pox[1][n] = rex * ny_g[n];
+      cache.pox[2][n] = rex * nz_g[n];
+    }
+  }
+
+  pox[0] = cache.pox[0];
+  pox[1] = cache.pox[1];
+  pox[2] = cache.pox[2];
+}
+
+double *surface_integral::get_shellf_buffer(int num_var)
+{
+  double *&buffer = shellf_cache[num_var];
+  if (!buffer)
+    buffer = new double[n_tot * num_var];
+  return buffer;
+}
+
+void surface_integral::release_cached_buffers()
+{
+  for (map<double, SpherePointCache>::iterator it = sphere_point_cache.begin(); it != sphere_point_cache.end(); ++it)
+  {
+    delete[] it->second.pox[0];
+    delete[] it->second.pox[1];
+    delete[] it->second.pox[2];
+    it->second.pox[0] = it->second.pox[1] = it->second.pox[2] = 0;
+  }
+  sphere_point_cache.clear();
+
+  for (map<int, double *>::iterator it = shellf_cache.begin(); it != shellf_cache.end(); ++it)
+    delete[] it->second;
+  shellf_cache.clear();
+}
+//|----------------------------------------------------------------
+//  spin weighted spinw component of psi4, general routine
+//  l takes from spinw to maxl; m takes from -l to l
 //|----------------------------------------------------------------
 void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4,
                                 int spinw, int maxl, int NN, double *RP, double *IP,
@@ -209,16 +254,9 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
  MyList<var> *DG_List = new MyList<var>(Rpsi4);
  DG_List->insert(Ipsi4);

-  int n;
-  double *pox[3];
-  for (int i = 0; i < 3; i++)
-    pox[i] = new double[n_tot];
-  for (n = 0; n < n_tot; n++)
-  {
-    pox[0][n] = rex * nx_g[n];
-    pox[1][n] = rex * ny_g[n];
-    pox[2][n] = rex * nz_g[n];
-  }
+  int n;
+  double *pox[3];
+  get_surface_points(rex, pox);

  int mp, Lp, Nmin, Nmax;
  mp = n_tot / cpusize;
@@ -234,8 +272,7 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
    Nmax = Nmin + mp - 1;
  }

-  double *shellf;
-  shellf = new double[n_tot * InList];
+  double *shellf = get_shellf_buffer(InList);

  GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax);

@@ -375,14 +412,10 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *

  //|------= Free memory.

-  delete[] pox[0];
-  delete[] pox[1];
-  delete[] pox[2];
-  delete[] shellf;
-  delete[] RP_out;
-  delete[] IP_out;
-  DG_List->clearList();
-}
+  delete[] RP_out;
+  delete[] IP_out;
+  DG_List->clearList();
+}
 void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4,
                                 int spinw, int maxl, int NN, double *RP, double *IP,
                                 monitor *Monitor, MPI_Comm Comm_here) // NN is the length of RP and IP
@@ -402,19 +435,11 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
  MyList<var> *DG_List = new MyList<var>(Rpsi4);
  DG_List->insert(Ipsi4);

-  int n;
-  double *pox[3];
-  for (int i = 0; i < 3; i++)
-    pox[i] = new double[n_tot];
-  for (n = 0; n < n_tot; n++)
-  {
-    pox[0][n] = rex * nx_g[n];
-    pox[1][n] = rex * ny_g[n];
-    pox[2][n] = rex * nz_g[n];
-  }
-
-  double *shellf;
-  shellf = new double[n_tot * InList];
+  int n;
+  double *pox[3];
+  get_surface_points(rex, pox);
+
+  double *shellf = get_shellf_buffer(InList);

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

@@ -577,14 +602,10 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *

  //|------= Free memory.

-  delete[] pox[0];
-  delete[] pox[1];
-  delete[] pox[2];
-  delete[] shellf;
-  delete[] RP_out;
-  delete[] IP_out;
-  DG_List->clearList();
-}
+  delete[] RP_out;
+  delete[] IP_out;
+  DG_List->clearList();
+}
 //|----------------------------------------------------------------
 //  for shell patch
 //|----------------------------------------------------------------
@@ -597,19 +618,11 @@ void surface_integral::surf_Wave(double rex, int lev, ShellPatch *GH, var *Rpsi4
  MyList<var> *DG_List = new MyList<var>(Rpsi4);
  DG_List->insert(Ipsi4);

-  int n;
-  double *pox[3];
-  for (int i = 0; i < 3; i++)
-    pox[i] = new double[n_tot];
-  for (n = 0; n < n_tot; n++)
-  {
-    pox[0][n] = rex * nx_g[n];
-    pox[1][n] = rex * ny_g[n];
-    pox[2][n] = rex * nz_g[n];
-  }
+  int n;
+  double *pox[3];
+  get_surface_points(rex, pox);

-  double *shellf;
-  shellf = new double[n_tot * InList];
+  double *shellf = get_shellf_buffer(InList);

  GH->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry);

@@ -2570,12 +2583,8 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var
  Rout[5] = sy;
  Rout[6] = sz;

-  delete[] pox[0];
-  delete[] pox[1];
-  delete[] pox[2];
-  delete[] shellf;
-  DG_List->clearList();
-}
+  DG_List->clearList();
+}
 void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var *trK,
                                     var *gxx, var *gxy, var *gxz, var *gyy, var *gyz, var *gzz,
                                     var *Axx, var *Axy, var *Axz, var *Ayy, var *Ayz, var *Azz,
@@ -2637,19 +2646,11 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var
  DG_List->insert(Ayz);
  DG_List->insert(Azz);

-  int n;
-  double *pox[3];
-  for (int i = 0; i < 3; i++)
-    pox[i] = new double[n_tot];
-  for (n = 0; n < n_tot; n++)
-  {
-    pox[0][n] = rex * nx_g[n];
-    pox[1][n] = rex * ny_g[n];
-    pox[2][n] = rex * nz_g[n];
-  }
-
-  double *shellf;
-  shellf = new double[n_tot * InList];
+  int n;
+  double *pox[3];
+  get_surface_points(rex, pox);
+
+  double *shellf = get_shellf_buffer(InList);

  // we have assumed there is only one box on this level,
  // so we do not need loop boxes
@@ -2839,12 +2840,8 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var
  Rout[5] = sy;
  Rout[6] = sz;

-  delete[] pox[0];
-  delete[] pox[1];
-  delete[] pox[2];
-  delete[] shellf;
-  DG_List->clearList();
-}
+  DG_List->clearList();
+}
 //|----------------------------------------------------------------
 //  for shell patch
 //|----------------------------------------------------------------
--- a/AMSS_NCKU_source/surface_integral.h
+++ b/AMSS_NCKU_source/surface_integral.h
@@ -20,25 +20,41 @@ using namespace std;
 #include "cgh.h"
 #include "ShellPatch.h"
 #include "NullShellPatch.h"
-#include "NullShellPatch2.h"
-#include "var.h"
-#include "monitor.h"
+#include "NullShellPatch2.h"
+#include "var.h"
+#include "monitor.h"
+#include <map>

 class surface_integral
 {

-private:
-	int Symmetry, factor;
-	int N_theta, N_phi; // Number of points in Theta & Phi directions
-	double dphi, dcostheta;
-	double *arcostheta, *wtcostheta;
-	int n_tot; // size of arrays
-
-	double *nx_g, *ny_g, *nz_g; // global list of unit normals
-	int myrank, cpusize;
-
-public:
-	surface_integral(int iSymmetry);
+private:
+	struct SpherePointCache
+	{
+		double *pox[3];
+		SpherePointCache()
+		{
+			pox[0] = pox[1] = pox[2] = 0;
+		}
+	};
+
+	int Symmetry, factor;
+	int N_theta, N_phi; // Number of points in Theta & Phi directions
+	double dphi, dcostheta;
+	double *arcostheta, *wtcostheta;
+	int n_tot; // size of arrays
+
+	double *nx_g, *ny_g, *nz_g; // global list of unit normals
+	int myrank, cpusize;
+	map<double, SpherePointCache> sphere_point_cache;
+	map<int, double *> shellf_cache;
+
+	void get_surface_points(double rex, double **pox);
+	double *get_shellf_buffer(int num_var);
+	void release_cached_buffers();
+
+public:
+	surface_integral(int iSymmetry);
 	~surface_integral();

 	void surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *Ipsi4,
--- a/makefile_and_run.py
+++ b/makefile_and_run.py
@@ -9,6 +9,7 @@


 import AMSS_NCKU_Input as input_data
+import os
 import subprocess
 import time

@@ -57,6 +58,48 @@ BUILD_JOBS = 64
 ##################################################################


+##################################################################
+
+def prepare_gpu_runtime_env():
+    """
+    Create a user-private CUDA MPS environment for GPU runs.
+
+    On shared machines another user's daemon may already occupy the default
+    /tmp/nvidia-mps pipe directory, which makes plain cudaSetDevice/cudaMalloc
+    fail with cudaErrorMpsConnectionFailed.  Binding AMSS-NCKU to a private
+    pipe directory avoids cross-user interference.
+    """
+    env = os.environ.copy()
+
+    pipe_dir = env.get("CUDA_MPS_PIPE_DIRECTORY", f"/tmp/amss-ncku-mps-{os.getuid()}")
+    log_dir  = env.get("CUDA_MPS_LOG_DIRECTORY",  f"/tmp/amss-ncku-mps-log-{os.getuid()}")
+
+    os.makedirs(pipe_dir, exist_ok=True)
+    os.makedirs(log_dir, exist_ok=True)
+
+    env["CUDA_MPS_PIPE_DIRECTORY"] = pipe_dir
+    env["CUDA_MPS_LOG_DIRECTORY"] = log_dir
+
+    control_socket = os.path.join(pipe_dir, "control")
+    if not os.path.exists(control_socket):
+        start = subprocess.run(
+            ["nvidia-cuda-mps-control", "-d"],
+            env=env,
+            stdout=subprocess.DEVNULL,
+            stderr=subprocess.DEVNULL,
+        )
+        if start.returncode != 0:
+            print(f" Warning: failed to start private CUDA MPS daemon in {pipe_dir}")
+        else:
+            print(f" Using private CUDA MPS pipe directory: {pipe_dir}")
+    else:
+        print(f" Using existing private CUDA MPS pipe directory: {pipe_dir}")
+
+    return env
+
+##################################################################
+
+

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

@@ -146,16 +189,29 @@ def run_ABE():

    ## Define the command to run; cast other values to strings as needed
    
+    run_env = None
+
    if (input_data.GPU_Calculation == "no"):
        mpi_command         = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
        #mpi_command         = " mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
        mpi_command_outfile = "ABE_out.log"
    elif (input_data.GPU_Calculation == "yes"):
-        mpi_command         = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + " ./ABEGPU"
+        run_env = prepare_gpu_runtime_env()
+        if int(input_data.MPI_processes) == 1:
+            mpi_command = "./ABEGPU"
+        else:
+            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)
+    mpi_process = subprocess.Popen(
+        mpi_command,
+        shell=True,
+        stdout=subprocess.PIPE,
+        stderr=subprocess.STDOUT,
+        text=True,
+        env=run_env,
+    )

    ## Write ABE run output to file while printing to stdout
    with open(mpi_command_outfile, 'w') as file0:
Author	SHA1	Message	Date
CGH0S7	3b16795e78	Refresh synced GPU regions incrementally	2026-04-09 17:07:31 +08:00
CGH0S7	5b00d49070	Reduce staged GPU host-device copies	2026-04-09 16:44:08 +08:00
CGH0S7	42e851d19a	Cache repeated interpolation plans	2026-04-09 15:21:01 +08:00
CGH0S7	06fa643365	Refine batched CUDA interpolation kernel	2026-04-09 15:06:11 +08:00
CGH0S7	c47349b7a9	Add batched CUDA patch interpolation path	2026-04-09 14:56:01 +08:00
CGH0S7	ad999e4c5a	Add guarded GPU prolong3 path scaffold	2026-04-09 14:28:36 +08:00
CGH0S7	e1e3b4a448	Reduce GPU RK4 transfer overhead	2026-04-09 12:11:40 +08:00
CGH0S7	49409645c0	Stabilize GPU output path and MPI sync	2026-04-09 10:57:49 +08:00
CGH0S7	4e3946a4f0	Persist GPU RK4 stage caches	2026-04-08 20:59:15 +08:00
CGH0S7	a0af9b8804	Trim GPU main-path transfer overhead	2026-04-08 20:16:25 +08:00
CGH0S7	01ac1f9250	Cache GPU main-path device buffers	2026-04-08 19:43:17 +08:00
CGH0S7	ea470737db	Add runnable GPU main-path prototype	2026-04-08 19:14:37 +08:00
CGH0S7	8c1f4d8108	迁移C算子的循环融合和临时量消除	2026-03-03 16:20:15 +08:00
CGH0S7	d310ef918b	bssn_rhs(fortran): migrate C kernel loop-fusion optimizations	2026-03-03 16:20:15 +08:00
CGH0S7	b35e1b289f	设置开关关闭内存打印统计	2026-03-03 16:17:47 +08:00
CGH0S7	05851b2c59	关闭静态负载	2026-03-03 16:17:47 +08:00
ianchb	3b39583d67	fix(bssn_rhs)	2026-03-03 16:06:33 +08:00