Reduce redundant GPU host buffer preparation

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"

AMSS_NCKU_source/Block.C

@@ -9,8 +9,11 @@
 #include <new>
 using namespace std;
 
-#include "Block.h"
-#include "misc.h"
+#include "Block.h"
+#include "misc.h"
+#ifdef USE_GPU
+#include "bssn_gpu.h"
+#endif
 
 Block::Block(int DIM, int *shapei, double *bboxi, int ranki, int ingfsi, int fngfsi, int levi, const int cgpui) : rank(ranki), ingfs(ingfsi), fngfs(fngfsi), lev(levi), cgpu(cgpui)
 {
@@ -95,14 +98,17 @@ Block::Block(int DIM, int *shapei, double *bboxi, int ranki, int ingfsi, int fng
   }
 #endif
 }
-Block::~Block()
-{
-  int myrank;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  if (myrank == rank)
-  {
-    for (int i = 0; i < dim; i++)
-      delete[] X[i];
+Block::~Block()
+{
+  int myrank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  if (myrank == rank)
+  {
+#ifdef USE_GPU
+    bssn_gpu_clear_cached_device_buffers();
+#endif
+    for (int i = 0; i < dim; i++)
+      delete[] X[i];
     for (int i = 0; i < ingfs; i++)
       free(igfs[i]);
     delete[] igfs;

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;

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,19 +124,20 @@ 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);
-  void Sync_finish(SyncCache &cache, AsyncSyncState &state,
-                   MyList<var> *VarList, int Symmetry);
+  void Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry,
+                  SyncCache &cache, AsyncSyncState &state);
+  void Sync_finish(SyncCache &cache, AsyncSyncState &state,
+                   MyList<var> *VarList, int Symmetry, bool unpack_to_host = true);
   void OutBdLow2Hi(Patch *Patc, Patch *Patf,
                    MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
                    int Symmetry);

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"
@@ -28,8 +29,12 @@ using namespace std;
 #include "rungekutta4_rout.h"
 #include "sommerfeld_rout.h"
 #include "getnp4.h"
-#include "shellfunctions.h"
-#include "parameters.h"
+#include "shellfunctions.h"
+#include "parameters.h"
+#ifdef USE_GPU
+#include "bssn_macro.h"
+#include "bssn_gpu.h"
+#endif
 
 #ifdef With_AHF
 #include "derivatives.h"
@@ -740,11 +745,12 @@ void bssn_class::Initialize()
   // Initialize sync caches (per-level, for predictor and corrector)
   sync_cache_pre = new Parallel::SyncCache[GH->levels];
   sync_cache_cor = new Parallel::SyncCache[GH->levels];
-  sync_cache_rp_coarse = new Parallel::SyncCache[GH->levels];
-  sync_cache_rp_fine = new Parallel::SyncCache[GH->levels];
-  sync_cache_restrict = new Parallel::SyncCache[GH->levels];
-  sync_cache_outbd = new Parallel::SyncCache[GH->levels];
-}
+  sync_cache_rp_coarse = new Parallel::SyncCache[GH->levels];
+  sync_cache_rp_fine = new Parallel::SyncCache[GH->levels];
+  sync_cache_restrict = new Parallel::SyncCache[GH->levels];
+  sync_cache_outbd = new Parallel::SyncCache[GH->levels];
+  sync_cache_psi4 = new Parallel::SyncCache[GH->levels];
+}
 
 //================================================================================================
 
@@ -756,8 +762,8 @@ void bssn_class::Initialize()
 
 //================================================================================================
 
-bssn_class::~bssn_class()
-{
+bssn_class::~bssn_class()
+{
 #ifdef With_AHF
   AHList->clearList();
   AHDList->clearList();
@@ -1014,12 +1020,30 @@ bssn_class::~bssn_class()
       sync_cache_rp_coarse[i].destroy();
     delete[] sync_cache_rp_coarse;
   }
-  if (sync_cache_rp_fine)
-  {
-    for (int i = 0; i < GH->levels; i++)
-      sync_cache_rp_fine[i].destroy();
-    delete[] sync_cache_rp_fine;
-  }
+  if (sync_cache_rp_fine)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_rp_fine[i].destroy();
+    delete[] sync_cache_rp_fine;
+  }
+  if (sync_cache_restrict)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_restrict[i].destroy();
+    delete[] sync_cache_restrict;
+  }
+  if (sync_cache_outbd)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_outbd[i].destroy();
+    delete[] sync_cache_outbd;
+  }
+  if (sync_cache_psi4)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_psi4[i].destroy();
+    delete[] sync_cache_psi4;
+  }
 
   delete GH;
 #ifdef WithShell
@@ -1052,8 +1076,25 @@ bssn_class::~bssn_class()
   delete ConVMonitor;
   delete Waveshell;
 
-  delete CheckPoint;
-}
+  delete CheckPoint;
+}
+
+void bssn_class::InvalidateSyncCaches()
+{
+  if (!GH)
+    return;
+
+  for (int il = 0; il < GH->levels; il++)
+  {
+    sync_cache_pre[il].invalidate();
+    sync_cache_cor[il].invalidate();
+    sync_cache_rp_coarse[il].invalidate();
+    sync_cache_rp_fine[il].invalidate();
+    sync_cache_restrict[il].invalidate();
+    sync_cache_outbd[il].invalidate();
+    sync_cache_psi4[il].invalidate();
+  }
+}
 
 //================================================================================================
 
@@ -2036,9 +2077,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
      {
@@ -2223,7 +2265,7 @@ void bssn_class::Evolve(int Steps)
     GH->Regrid(Symmetry, BH_num, Porgbr, Porg0,
                SynchList_cor, OldStateList, StateList, SynchList_pre,
                fgt(PhysTime - dT_mon, StartTime, dT_mon / 2), ErrorMonitor);
-    for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+    InvalidateSyncCaches();
 #endif
 
 #if (REGLEV == 0 && (PSTR == 1 || PSTR == 2))
@@ -2296,18 +2338,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
@@ -2441,7 +2486,7 @@ void bssn_class::RecursiveStep(int lev)
   if (GH->Regrid_Onelevel(lev, Symmetry, BH_num, Porgbr, Porg0,
                       SynchList_cor, OldStateList, StateList, SynchList_pre,
                       fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
-  for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+  InvalidateSyncCaches();
 #endif
 }
 
@@ -2620,7 +2665,7 @@ void bssn_class::ParallelStep()
   if (GH->Regrid_Onelevel(GH->mylev, Symmetry, BH_num, Porgbr, Porg0,
                       SynchList_cor, OldStateList, StateList, SynchList_pre,
                       fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
-  for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+  InvalidateSyncCaches();
 #endif
 }
 
@@ -2787,7 +2832,7 @@ void bssn_class::ParallelStep()
         if (GH->Regrid_Onelevel(lev + 1, Symmetry, BH_num, Porgbr, Porg0,
                             SynchList_cor, OldStateList, StateList, SynchList_pre,
                             fgt(PhysTime - dT_levp1, StartTime, dT_levp1 / 2), ErrorMonitor))
-        for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+        InvalidateSyncCaches();
 
         //               a_stream.clear();
         //               a_stream.str("");
@@ -2802,7 +2847,7 @@ void bssn_class::ParallelStep()
       if (GH->Regrid_Onelevel(lev, Symmetry, BH_num, Porgbr, Porg0,
                           SynchList_cor, OldStateList, StateList, SynchList_pre,
                           fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
-      for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+      InvalidateSyncCaches();
 
       //               a_stream.clear();
       //               a_stream.str("");
@@ -2821,7 +2866,7 @@ void bssn_class::ParallelStep()
           if (GH->Regrid_Onelevel(lev - 1, Symmetry, BH_num, Porgbr, Porg0,
                               SynchList_cor, OldStateList, StateList, SynchList_pre,
                               fgt(PhysTime - dT_lev, StartTime, dT_levm1 / 2), ErrorMonitor))
-          for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+          InvalidateSyncCaches();
 
           //               a_stream.clear();
           //               a_stream.str("");
@@ -2837,7 +2882,7 @@ void bssn_class::ParallelStep()
           if (GH->Regrid_Onelevel(lev - 1, Symmetry, BH_num, Porgbr, Porg0,
                               SynchList_cor, OldStateList, StateList, SynchList_pre,
                               fgt(PhysTime - dT_lev, StartTime, dT_levm1 / 2), ErrorMonitor))
-          for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+          InvalidateSyncCaches();
 
           //               a_stream.clear();
           //               a_stream.str("");
@@ -3026,9 +3071,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));
 
@@ -6248,7 +6298,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_cached(GH->PatL[lev], DG_List, Symmetry, sync_cache_psi4[lev]);
 #endif
 
 #ifdef WithShell
@@ -6903,10 +6953,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++)
@@ -7250,9 +7300,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
@@ -7272,12 +7322,15 @@ void bssn_class::Constraint_Out()
           MyList<Block> *BP = Pp->data->blb;
           while (BP)
           {
-            Block *cg = BP->data;
-            if (myrank == cg->rank)
-            {
-              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
-                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
-                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
+            Block *cg = BP->data;
+            if (myrank == cg->rank)
+            {
+#ifdef USE_GPU
+              gpu_rhs(CALLED_BY_CONSTRAINT_CONS_ONLY, myrank, RHS_PARA_CALLED_Constraint_Out);
+#else
+              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
+                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
+                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
                                  cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
                                  cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], 
                                  cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn],
@@ -7305,11 +7358,12 @@ void bssn_class::Constraint_Out()
                                  cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn],
                                  cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], 
                                  cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn],
-                                 cg->fgfs[Cons_Ham->sgfn],
-                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
-                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
-                                 Symmetry, lev, ndeps, pre);
-            }
+                                 cg->fgfs[Cons_Ham->sgfn],
+                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
+                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
+                                 Symmetry, lev, ndeps, pre);
+#endif
+            }
             if (BP == Pp->data->ble)
               break;
             BP = BP->next;
@@ -7317,7 +7371,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
@@ -7539,7 +7593,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");
   }
 }
 
@@ -7775,12 +7829,15 @@ void bssn_class::Interp_Constraint(bool infg)
           MyList<Block> *BP = Pp->data->blb;
           while (BP)
           {
-            Block *cg = BP->data;
-            if (myrank == cg->rank)
-            {
-              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
-                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
-                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
+            Block *cg = BP->data;
+            if (myrank == cg->rank)
+            {
+#ifdef USE_GPU
+              gpu_rhs(CALLED_BY_CONSTRAINT_CONS_ONLY, myrank, RHS_PARA_CALLED_Interp_Constraint);
+#else
+              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
+                                 cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
+                                 cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
                                  cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
                                  cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], 
                                  cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn],
@@ -7808,11 +7865,12 @@ void bssn_class::Interp_Constraint(bool infg)
                                  cg->fgfs[Gamzyy->sgfn], cg->fgfs[Gamzyz->sgfn], cg->fgfs[Gamzzz->sgfn],
                                  cg->fgfs[Rxx->sgfn], cg->fgfs[Rxy->sgfn], cg->fgfs[Rxz->sgfn], 
                                  cg->fgfs[Ryy->sgfn], cg->fgfs[Ryz->sgfn], cg->fgfs[Rzz->sgfn],
-                                 cg->fgfs[Cons_Ham->sgfn],
-                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
-                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
-                                 Symmetry, lev, ndeps, pre);
-            }
+                                 cg->fgfs[Cons_Ham->sgfn],
+                                 cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
+                                 cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
+                                 Symmetry, lev, ndeps, pre);
+#endif
+            }
             if (BP == Pp->data->ble)
               break;
             BP = BP->next;
@@ -7820,7 +7878,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
@@ -8078,7 +8136,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--)
@@ -8391,12 +8449,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);
@@ -8405,14 +8469,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;
+            }
         }
     }
 

AMSS_NCKU_source/bssn_class.h

@@ -128,10 +128,11 @@ public:
 
        Parallel::SyncCache *sync_cache_pre;  // per-level cache for predictor sync
        Parallel::SyncCache *sync_cache_cor;  // per-level cache for corrector sync
-       Parallel::SyncCache *sync_cache_rp_coarse;  // RestrictProlong sync on PatL[lev-1]
-       Parallel::SyncCache *sync_cache_rp_fine;    // RestrictProlong sync on PatL[lev]
-       Parallel::SyncCache *sync_cache_restrict;   // cached Restrict in RestrictProlong
-       Parallel::SyncCache *sync_cache_outbd;      // cached OutBdLow2Hi in RestrictProlong
+       Parallel::SyncCache *sync_cache_rp_coarse;  // RestrictProlong sync on PatL[lev-1]
+       Parallel::SyncCache *sync_cache_rp_fine;    // RestrictProlong sync on PatL[lev]
+       Parallel::SyncCache *sync_cache_restrict;   // cached Restrict in RestrictProlong
+       Parallel::SyncCache *sync_cache_outbd;      // cached OutBdLow2Hi in RestrictProlong
+       Parallel::SyncCache *sync_cache_psi4;       // cached Psi4 sync on PatL[lev]
 
        monitor *ErrorMonitor, *Psi4Monitor, *BHMonitor, *MAPMonitor;
        monitor *ConVMonitor;
@@ -171,16 +172,20 @@ 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() {};
+       void InvalidateSyncCaches();
+       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:

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

AMSS_NCKU_source/bssn_cuda_step.C

@@ -0,0 +1,555 @@
+#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;
+#if (MAPBH == 0)
+  const bool need_host_stage_sync = (BH_num > 0 && lev == GH->levels - 1);
+#else
+  const bool need_host_stage_sync = false;
+#endif
+  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 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;
+        }
+      };
+
+  auto refresh_stage_host_before_sync =
+      [&](MyList<var> *var_list, Parallel::SyncCache &cache) -> bool {
+        if (!cache.valid || !cache.combined_src || myrank < 0 || myrank >= cache.cpusize)
+          return false;
+
+        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)
+            {
+              std::vector<Parallel::gridseg *> local_segments;
+              MyList<Parallel::gridseg> *seg = cache.combined_src[myrank];
+              while (seg)
+              {
+                if (seg->data && seg->data->Bg == cg)
+                  local_segments.push_back(seg->data);
+                seg = seg->next;
+              }
+
+              if (!local_segments.empty())
+              {
+                MyList<var> *var_it = var_list;
+                while (var_it)
+                {
+                  double *host_ptr = cg->fgfs[var_it->data->sgfn];
+                  for (size_t i = 0; i < local_segments.size(); ++i)
+                  {
+                    Parallel::gridseg *src_seg = local_segments[i];
+                    if (bssn_gpu_stage_download_region(host_ptr,
+                                                       cg->shape,
+                                                       cg->bbox,
+                                                       cg->bbox + dim,
+                                                       src_seg->shape,
+                                                       src_seg->llb))
+                    {
+                      cerr << "GPU sync region download failure: lev=" << lev
+                           << " var=" << var_it->data->name
+                           << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                           << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                           << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+                      ERROR = 1;
+                      return true;
+                    }
+                  }
+                  var_it = var_it->next;
+                }
+              }
+            }
+
+            if (block_it == patch_it->data->ble)
+              break;
+            block_it = block_it->next;
+          }
+          patch_it = patch_it->next;
+        }
+
+        return true;
+      };
+
+  auto can_pack_sync_from_device =
+      [&](MyList<var> *var_list, Parallel::SyncCache &cache) -> bool {
+        if (!cache.valid || !cache.combined_src || myrank < 0 || myrank >= cache.cpusize)
+          return false;
+
+        MyList<Parallel::gridseg> *seg = cache.combined_src[myrank];
+        while (seg)
+        {
+          MyList<var> *var_it = var_list;
+          while (var_it)
+          {
+            if (!bssn_gpu_find_device_buffer(seg->data->Bg->fgfs[var_it->data->sgfn]))
+              return false;
+            var_it = var_it->next;
+          }
+          seg = seg->next;
+        }
+        return true;
+      };
+
+  MyList<Patch> *Pp = GH->PatL[lev];
+  while (Pp)
+  {
+    MyList<Block> *BP = Pp->data->blb;
+    while (BP)
+    {
+      Block *cg = BP->data;
+      if (myrank == cg->rank)
+      {
+        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 && !sync_cache_pre[lev].valid)
+          stage_download_var_list(cg, SynchList_pre);
+      }
+      if (BP == Pp->data->ble)
+        break;
+      BP = BP->next;
+    }
+    Pp = Pp->next;
+  }
+
+  if (!ERROR && sync_cache_pre[lev].valid && !can_pack_sync_from_device(SynchList_pre, sync_cache_pre[lev]))
+    refresh_stage_host_before_sync(SynchList_pre, sync_cache_pre[lev]);
+
+  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, need_host_stage_sync);
+  if (!ERROR && need_host_stage_sync)
+    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 && (!sync_cache_cor[lev].valid || iter_count == 3))
+            stage_download_var_list(cg, SynchList_cor);
+        }
+
+        if (BP == Pp->data->ble)
+          break;
+        BP = BP->next;
+      }
+      Pp = Pp->next;
+    }
+
+    if (!ERROR && sync_cache_cor[lev].valid && iter_count < 3 &&
+        !can_pack_sync_from_device(SynchList_cor, sync_cache_cor[lev]))
+      refresh_stage_host_before_sync(SynchList_cor, sync_cache_cor[lev]);
+
+    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);
+    const bool unpack_cor_to_host = (iter_count == 3) || need_host_stage_sync;
+    Parallel::Sync_finish(sync_cache_cor[lev], async_cor, SynchList_cor, Symmetry, unpack_cor_to_host);
+    if (!ERROR && iter_count < 3 && unpack_cor_to_host)
+      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

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,44 @@ 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();
+void bssn_gpu_release_pinned_host_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);
+void bssn_gpu_prepare_host_buffer(const double *host_ptr, int count);
+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);
+int bssn_gpu_stage_download_region(double *host_ptr,
+                                   const int *full_shape,
+                                   const double *full_llb,
+                                   const double *full_uub,
+                                   const int *region_shape,
+                                   const double *region_llb);
+int bssn_gpu_stage_download_region_to_buffer(const double *host_src_ptr,
+                                             const int *full_shape,
+                                             const double *full_llb,
+                                             const double *full_uub,
+                                             const int *region_shape,
+                                             const double *region_llb,
+                                             double *host_dst_ptr);
+int bssn_gpu_stage_upload_buffer_to_region(const double *host_src_ptr,
+                                           double *host_dst_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

AMSS_NCKU_source/bssn_macro.h

@@ -65,9 +65,10 @@ if(TIME_COUNT_EACH_RANK == 1){\
      }\
 }
 
-//3---------------------GPU---------------------
-#define CALLED_BY_STEP 0
-#define CALLED_BY_CONSTRAINT 1
+//3---------------------GPU---------------------
+#define CALLED_BY_STEP 0
+#define CALLED_BY_CONSTRAINT 1
+#define CALLED_BY_CONSTRAINT_CONS_ONLY 2
 
 #define RHS_PARA_CALLED_FIRST_TIME cg->shape,TRK4,cg->X[0],cg->X[1],cg->X[2],cg->fgfs[phi0->sgfn],cg->fgfs[trK0->sgfn],cg->fgfs[gxx0->sgfn],cg->fgfs[gxy0->sgfn],cg->fgfs[gxz0->sgfn],cg->fgfs[gyy0->sgfn],cg->fgfs[gyz0->sgfn],cg->fgfs[gzz0->sgfn],cg->fgfs[Axx0->sgfn],cg->fgfs[Axy0->sgfn],cg->fgfs[Axz0->sgfn],cg->fgfs[Ayy0->sgfn],cg->fgfs[Ayz0->sgfn],cg->fgfs[Azz0->sgfn],cg->fgfs[Gmx0->sgfn],cg->fgfs[Gmy0->sgfn],cg->fgfs[Gmz0->sgfn],cg->fgfs[Lap0->sgfn],cg->fgfs[Sfx0->sgfn],cg->fgfs[Sfy0->sgfn],cg->fgfs[Sfz0->sgfn],cg->fgfs[dtSfx0->sgfn],cg->fgfs[dtSfy0->sgfn],cg->fgfs[dtSfz0->sgfn],cg->fgfs[phi_rhs->sgfn],cg->fgfs[trK_rhs->sgfn],cg->fgfs[gxx_rhs->sgfn],cg->fgfs[gxy_rhs->sgfn],cg->fgfs[gxz_rhs->sgfn],cg->fgfs[gyy_rhs->sgfn],cg->fgfs[gyz_rhs->sgfn],cg->fgfs[gzz_rhs->sgfn],cg->fgfs[Axx_rhs->sgfn],cg->fgfs[Axy_rhs->sgfn],cg->fgfs[Axz_rhs->sgfn],cg->fgfs[Ayy_rhs->sgfn],cg->fgfs[Ayz_rhs->sgfn],cg->fgfs[Azz_rhs->sgfn],cg->fgfs[Gmx_rhs->sgfn],cg->fgfs[Gmy_rhs->sgfn],cg->fgfs[Gmz_rhs->sgfn],cg->fgfs[Lap_rhs->sgfn],cg->fgfs[Sfx_rhs->sgfn],cg->fgfs[Sfy_rhs->sgfn],cg->fgfs[Sfz_rhs->sgfn],cg->fgfs[dtSfx_rhs->sgfn],cg->fgfs[dtSfy_rhs->sgfn],cg->fgfs[dtSfz_rhs->sgfn],cg->fgfs[rho->sgfn],cg->fgfs[Sx->sgfn],cg->fgfs[Sy->sgfn],cg->fgfs[Sz->sgfn],cg->fgfs[Sxx->sgfn],cg->fgfs[Sxy->sgfn],cg->fgfs[Sxz->sgfn],cg->fgfs[Syy->sgfn],cg->fgfs[Syz->sgfn],cg->fgfs[Szz->sgfn],cg->fgfs[Gamxxx->sgfn],cg->fgfs[Gamxxy->sgfn],cg->fgfs[Gamxxz->sgfn],cg->fgfs[Gamxyy->sgfn],cg->fgfs[Gamxyz->sgfn],cg->fgfs[Gamxzz->sgfn],cg->fgfs[Gamyxx->sgfn],cg->fgfs[Gamyxy->sgfn],cg->fgfs[Gamyxz->sgfn],cg->fgfs[Gamyyy->sgfn],cg->fgfs[Gamyyz->sgfn],cg->fgfs[Gamyzz->sgfn],cg->fgfs[Gamzxx->sgfn],cg->fgfs[Gamzxy->sgfn],cg->fgfs[Gamzxz->sgfn],cg->fgfs[Gamzyy->sgfn],cg->fgfs[Gamzyz->sgfn],cg->fgfs[Gamzzz->sgfn],cg->fgfs[Rxx->sgfn],cg->fgfs[Rxy->sgfn],cg->fgfs[Rxz->sgfn],cg->fgfs[Ryy->sgfn],cg->fgfs[Ryz->sgfn],cg->fgfs[Rzz->sgfn],cg->fgfs[Cons_Ham->sgfn],cg->fgfs[Cons_Px->sgfn],cg->fgfs[Cons_Py->sgfn],cg->fgfs[Cons_Pz->sgfn],cg->fgfs[Cons_Gx->sgfn],cg->fgfs[Cons_Gy->sgfn],cg->fgfs[Cons_Gz->sgfn],Symmetry,lev,ndeps,pre
 

AMSS_NCKU_source/bssn_rhs_c.C

@@ -1022,9 +1022,16 @@ int f_compute_rhs_bssn(int *ex, double &T,
                         + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
 
             #if (GAUGE == 2)
-            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
+            {
+                const double chi_sqrt = sqrt(chin1[i]);
+                const double damping = ONE - chi_sqrt;
+                reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / (damping * damping);
+            }
             #else
-            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - chin1[i]), 2.0 );
+            {
+                const double damping = ONE - chin1[i];
+                reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / (damping * damping);
+            }
             #endif
 
             dtSfx_rhs[i] = Gamx_rhs[i] - reta[i] * dtSfx[i];
@@ -1040,9 +1047,16 @@ int f_compute_rhs_bssn(int *ex, double &T,
                         + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
 
             #if (GAUGE == 4)
-            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
+            {
+                const double chi_sqrt = sqrt(chin1[i]);
+                const double damping = ONE - chi_sqrt;
+                reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / (damping * damping);
+            }
             #else
-            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - chin1[i]), 2.0 );
+            {
+                const double damping = ONE - chin1[i];
+                reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / (damping * damping);
+            }
             #endif
 
             betax_rhs[i] = FF * Gamx[i] - reta[i] * betax[i];

AMSS_NCKU_source/cgh.C

@@ -23,10 +23,13 @@ using namespace std;
 #include <mpi.h>
 
 #include "macrodef.h"
-#include "misc.h"
-#include "cgh.h"
-#include "Parallel.h"
-#include "parameters.h"
+#include "misc.h"
+#include "cgh.h"
+#include "Parallel.h"
+#include "parameters.h"
+#ifdef USE_GPU
+#include "bssn_gpu.h"
+#endif
 
 //================================================================================================
 
@@ -881,13 +884,17 @@ void cgh::recompose_cgh(int nprocs, bool *lev_flag,
       tmPat = construct_patchlist(lev, Symmetry);
       // tmPat construction completes
       Parallel::distribute(tmPat, nprocs, ingfs, fngfs, false);
-      //    checkPatchList(tmPat,true);
-      bool CC = (lev > trfls);
-      Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
-
-      Parallel::KillBlocks(PatL[lev]);
-      PatL[lev]->destroyList();
-      PatL[lev] = tmPat;
+      //    checkPatchList(tmPat,true);
+      bool CC = (lev > trfls);
+      Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
+
+#ifdef USE_GPU
+      bssn_gpu_clear_cached_device_buffers();
+      bssn_gpu_release_pinned_host_buffers();
+#endif
+      Parallel::KillBlocks(PatL[lev]);
+      PatL[lev]->destroyList();
+      PatL[lev] = tmPat;
 #if (RPB == 1)
       Parallel::destroypsuList_bam(bdsul[lev]);
       Parallel::destroypsuList_bam(rsul[lev]);
@@ -910,13 +917,17 @@ void cgh::recompose_cgh(int nprocs, bool *lev_flag,
       tmPat = construct_patchlist(lev, Symmetry);
       // tmPat construction completes
       Parallel::distribute(tmPat, end_rank[lev] - start_rank[lev] + 1, ingfs, fngfs, false, start_rank[lev], end_rank[lev]);
-      //    checkPatchList(tmPat,true);
-      bool CC = (lev > trfls);
-      Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
-
-      Parallel::KillBlocks(PatL[lev]);
-      PatL[lev]->destroyList();
-      PatL[lev] = tmPat;
+      //    checkPatchList(tmPat,true);
+      bool CC = (lev > trfls);
+      Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
+
+#ifdef USE_GPU
+      bssn_gpu_clear_cached_device_buffers();
+      bssn_gpu_release_pinned_host_buffers();
+#endif
+      Parallel::KillBlocks(PatL[lev]);
+      PatL[lev]->destroyList();
+      PatL[lev] = tmPat;
 #if (RPB == 1)
 #error "not support yet"
 #endif
@@ -1518,13 +1529,17 @@ void cgh::recompose_cgh_Onelevel(int nprocs, int lev,
   tmPat = construct_patchlist(lev, Symmetry);
   // tmPat construction completes
   Parallel::distribute(tmPat, nprocs, ingfs, fngfs, false);
-  //    checkPatchList(tmPat,true);
-  bool CC = (lev > trfls);
-  Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
-
-  Parallel::KillBlocks(PatL[lev]);
-  PatL[lev]->destroyList();
-  PatL[lev] = tmPat;
+  //    checkPatchList(tmPat,true);
+  bool CC = (lev > trfls);
+  Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
+
+#ifdef USE_GPU
+  bssn_gpu_clear_cached_device_buffers();
+  bssn_gpu_release_pinned_host_buffers();
+#endif
+  Parallel::KillBlocks(PatL[lev]);
+  PatL[lev]->destroyList();
+  PatL[lev] = tmPat;
 }
 #elif (PSTR == 1 || PSTR == 2 || PSTR == 3)
 #warning "recompose_cgh_Onelevel is not implimented yet"
@@ -1540,14 +1555,18 @@ void cgh::recompose_cgh_Onelevel(int nprocs, int lev,
   // tmPat construction completes
   Parallel::distribute(tmPat, end_rank[lev] - start_rank[lev] + 1, ingfs, fngfs, false, start_rank[lev], end_rank[lev]);
   misc::tillherecheck(Commlev[lev], start_rank[lev], "after distribute");
-  //    checkPatchList(tmPat,true);
-  bool CC = (lev > trfls);
-  Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
-  misc::tillherecheck(Commlev[lev], start_rank[lev], "after fill_level_data");
-
-  Parallel::KillBlocks(PatL[lev]);
-  PatL[lev]->destroyList();
-  PatL[lev] = tmPat;
+  //    checkPatchList(tmPat,true);
+  bool CC = (lev > trfls);
+  Parallel::fill_level_data(tmPat, PatL[lev], PatL[lev - 1], OldList, StateList, FutureList, tmList, Symmetry, BB, CC);
+  misc::tillherecheck(Commlev[lev], start_rank[lev], "after fill_level_data");
+
+#ifdef USE_GPU
+  bssn_gpu_clear_cached_device_buffers();
+  bssn_gpu_release_pinned_host_buffers();
+#endif
+  Parallel::KillBlocks(PatL[lev]);
+  PatL[lev]->destroyList();
+  PatL[lev] = tmPat;
 }
 
 

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 $@
@@ -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

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

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
 //|----------------------------------------------------------------

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,

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: