fix(bssn_rhs)

OutBdLow2Hi_cached 读的是 coarse owned 区域（非 coarse ghost/buffer）
回退旧行为：编译时定义 RP_SYNC_COARSE_AFTER_RESTRICT=1

AMSS_NCKU_source/ABE.C

@@ -24,16 +24,18 @@ using namespace std;
 
 #include "misc.h"
 #include "macrodef.h"
-#ifdef USE_GPU
-extern void bssn_cuda_dump_stage_profile();
-#endif
 
 #ifndef ABEtype
 #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"
@@ -472,13 +474,10 @@ int main(int argc, char *argv[])
             cout << endl;
       }
 
-	      ADM->Evolve(Steps);
+      ADM->Evolve(Steps);
-#ifdef USE_GPU
-	      bssn_cuda_dump_stage_profile();
-#endif
 
-	      if (myrank == 0)
+      if (myrank == 0)
-	      {
+      {
             cout << endl;
             cout << " Total Evolve Time: "  << MPI_Wtime() - End_clock   << " seconds!" << endl;
             cout << " Total Running Time: " << MPI_Wtime() - Begin_clock << " seconds!" << endl;

AMSS_NCKU_source/Block.C

@@ -11,10 +11,6 @@ using namespace std;
 
 #include "Block.h"
 #include "misc.h"
-#ifdef USE_GPU
-#include "bssn_gpu.h"
-#include "bssn_cuda_ops.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)
 {
@@ -105,11 +101,6 @@ Block::~Block()
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
   if (myrank == rank)
   {
-#ifdef USE_GPU
-    bssn_gpu_clear_cached_device_buffers();
-    bssn_cuda_release_rk4_caches();
-    bssn_cuda_release_interp_caches();
-#endif
     for (int i = 0; i < dim; i++)
       delete[] X[i];
     for (int i = 0; i < ingfs; i++)

AMSS_NCKU_source/MPatch.C

@@ -7,7 +7,6 @@
 #include <string>
 #include <cmath>
 #include <new>
-#include <map>
 #include <vector>
 using namespace std;
 
@@ -15,82 +14,12 @@ using namespace std;
 #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 CachedInterpPlanEntry
-{
-  bool valid;
-  InterpPlanKey key;
-  vector<double> xvals;
-  vector<double> yvals;
-  vector<double> zvals;
-  CachedInterpPlan plan;
-
-  CachedInterpPlanEntry() : valid(false) {}
-};
-
 struct InterpBlockView
 {
   Block *bp;
@@ -249,309 +178,8 @@ int find_block_index_for_point(const BlockBinIndex &index, const double *pox, co
 
   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;
-}
-
-CachedInterpPlanEntry &interp_plan_cache_entry()
-{
-  static CachedInterpPlanEntry cache;
-  return cache;
-}
-
-bool same_interp_plan_key(const InterpPlanKey &lhs, const InterpPlanKey &rhs)
-{
-  return lhs.patch == rhs.patch &&
-         lhs.NN == rhs.NN &&
-         lhs.Symmetry == rhs.Symmetry &&
-         lhs.myrank == rhs.myrank;
-}
-
-bool same_interp_plan_points(const CachedInterpPlanEntry &cache, int NN, double **XX)
-{
-  if (static_cast<int>(cache.xvals.size()) != NN ||
-      static_cast<int>(cache.yvals.size()) != NN ||
-      static_cast<int>(cache.zvals.size()) != NN)
-    return false;
-
-  for (int j = 0; j < NN; ++j)
-  {
-    if (cache.xvals[j] != XX[0][j] ||
-        cache.yvals[j] != XX[1][j] ||
-        cache.zvals[j] != XX[2][j])
-      return false;
-  }
-  return true;
-}
-
-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)
-{
-  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;
-
-  CachedInterpPlanEntry &cache = interp_plan_cache_entry();
-  if (cache.valid &&
-      same_interp_plan_key(cache.key, key) &&
-      same_interp_plan_points(cache, NN, XX) &&
-      cache.plan.nblocks == static_cast<int>(block_index.views.size()))
-    return cache.plan;
-
-  cache.valid = true;
-  cache.key = key;
-  cache.xvals.assign(XX[0], XX[0] + NN);
-  cache.yvals.assign(XX[1], XX[1] + NN);
-  cache.zvals.assign(XX[2], XX[2] + NN);
-  cache.plan = CachedInterpPlan();
-  CachedInterpPlan &plan = cache.plan;
-  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;
-}
-
-void release_interp_plan_cache_internal()
-{
-  CachedInterpPlanEntry &cache = interp_plan_cache_entry();
-  cache.valid = false;
-  cache.xvals.clear();
-  cache.yvals.clear();
-  cache.zvals.clear();
-  cache.plan = CachedInterpPlan();
-}
-
-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
 
-void patch_release_interp_plan_cache()
-{
-  release_interp_plan_cache_internal();
-}
-
 Patch::Patch(int DIM, int *shapei, double *bboxi, int levi, bool buflog, int Symmetry) : lev(levi)
 {
 
@@ -895,15 +523,60 @@ 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);
-  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);
+  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++;
+        }
+      }
+    }
+  }
 
   // Replace MPI_Allreduce with per-owner MPI_Bcast:
   // Group consecutive points by owner rank and broadcast each group.
@@ -959,6 +632,7 @@ void Patch::Interp_Points(MyList<var> *VarList,
     }
   }
 
+  delete[] owner_rank;
 }
 void Patch::Interp_Points(MyList<var> *VarList,
                           int NN, double **XX,
@@ -987,21 +661,101 @@ 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);
-  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);
+  // --- 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++;
+        }
+      }
+    }
+  }
 
 #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 ---
   // Compute consumer_rank[j] using the same deterministic formula as surface_integral
   int *consumer_rank = new int[NN];
@@ -1121,6 +875,7 @@ void Patch::Interp_Points(MyList<var> *VarList,
   delete[] send_count;
   delete[] recv_count;
   delete[] consumer_rank;
+  delete[] owner_rank;
 
 #ifdef INTERP_LB_PROFILE
   {
@@ -1168,6 +923,12 @@ 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);
@@ -1178,10 +939,48 @@ void Patch::Interp_Points(MyList<var> *VarList,
     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);
+  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++;
+        }
+      }
+    }
+  }
 
   // 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
@@ -1211,6 +1010,7 @@ void Patch::Interp_Points(MyList<var> *VarList,
 
   MPI_Group_free(&world_group);
   MPI_Group_free(&local_group);
+  delete[] owner_rank;
 }
 void Patch::checkBlock()
 {

AMSS_NCKU_source/MPatch.h

@@ -52,6 +52,4 @@ public:
                      double *Shellf, MPI_Comm Comm_here);
 };
 
-void patch_release_interp_plan_cache();
-
 #endif /* PATCH_H */

AMSS_NCKU_source/Parallel.C

@@ -2,298 +2,8 @@
 #include "Parallel.h"
 #include "fmisc.h"
 #include "prolongrestrict.h"
-#include "bssn_cuda_ops.h"
-#include "bssn_gpu.h"
 #include "misc.h"
 #include "parameters.h"
-#include <cstring>
-
-#if defined(__GNUC__) || defined(__clang__)
-extern int bssn_cuda_prolong3_pack(int wei,
-                                   const double *llbc, const double *uubc, const int *extc, const double *func,
-                                   const double *llbf, const double *uubf, const int *extf, double *funf,
-                                   const double *llbp, const double *uubp,
-                                   const double *SoA, int symmetry) __attribute__((weak));
-extern int bssn_cuda_restrict3_pack(int wei,
-                                    const double *llbc, const double *uubc, const int *extc, double *func,
-                                    const double *llbf, const double *uubf, const int *extf, const double *funf,
-                                    const double *llbr, const double *uubr,
-                                    const double *SoA, int symmetry) __attribute__((weak));
-#endif
-
-namespace {
-
-const char *g_parallel_transfer_context = "Parallel::transfer";
-int parallel_next_transfer_tag()
-{
-  const char *ctx = g_parallel_transfer_context ? g_parallel_transfer_context : "Parallel::transfer";
-  unsigned int hash = 5381u;
-  while (*ctx)
-  {
-    hash = ((hash << 5) + hash) + static_cast<unsigned char>(*ctx);
-    ctx++;
-  }
-  return 32 + static_cast<int>(hash % 32000u);
-}
-
-struct ParallelTransferContextGuard
-{
-  const char *prev;
-  explicit ParallelTransferContextGuard(const char *context) : prev(g_parallel_transfer_context)
-  {
-    g_parallel_transfer_context = context;
-  }
-  ~ParallelTransferContextGuard()
-  {
-    g_parallel_transfer_context = prev;
-  }
-};
-
-bool parallel_can_gpu_pack_segments(MyList<Parallel::gridseg> *src, MyList<Parallel::gridseg> *dst,
-                                    int rank_in, MyList<var> *VarLists, MyList<var> *VarListd)
-{
-  int myrank;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  if (!src || !dst)
-    return false;
-
-  if (src->data->Bg->lev != dst->data->Bg->lev)
-    return false;
-
-  while (src && dst)
-  {
-    if ((dst->data->Bg->rank == rank_in) && (src->data->Bg->rank == myrank))
-    {
-      MyList<var> *varls = VarLists;
-      MyList<var> *varld = VarListd;
-      while (varls && varld)
-      {
-        (void)varld;
-        if (!bssn_gpu_find_device_buffer(src->data->Bg->fgfs[varls->data->sgfn]))
-          return false;
-        varls = varls->next;
-        varld = varld->next;
-      }
-      if (varls || varld)
-        return false;
-    }
-    src = src->next;
-    dst = dst->next;
-  }
-  return true;
-}
-
-bool parallel_gpu_pack_segments(double *data,
-                                MyList<Parallel::gridseg> *src, MyList<Parallel::gridseg> *dst,
-                                int rank_in, MyList<var> *VarLists, MyList<var> *VarListd)
-{
-  if (!data)
-    return false;
-
-  int myrank;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  if (!src || !dst)
-    return false;
-  if (src->data->Bg->lev != dst->data->Bg->lev)
-    return false;
-
-  int size_out = 0;
-  while (src && dst)
-  {
-    if ((dst->data->Bg->rank == rank_in) && (src->data->Bg->rank == myrank))
-    {
-      MyList<var> *varls = VarLists;
-      MyList<var> *varld = VarListd;
-      while (varls && varld)
-      {
-        (void)varld;
-        if (bssn_gpu_stage_download_region_to_buffer(src->data->Bg->fgfs[varls->data->sgfn],
-                                                     src->data->Bg->shape,
-                                                     src->data->Bg->bbox,
-                                                     src->data->Bg->bbox + dim,
-                                                     dst->data->shape,
-                                                     dst->data->llb,
-                                                     data + size_out))
-          return false;
-        size_out += dst->data->shape[0] * dst->data->shape[1] * dst->data->shape[2];
-        varls = varls->next;
-        varld = varld->next;
-      }
-      if (varls || varld)
-        return false;
-    }
-    src = src->next;
-    dst = dst->next;
-  }
-  return true;
-}
-
-bool parallel_can_gpu_unpack_segments(MyList<Parallel::gridseg> *src, MyList<Parallel::gridseg> *dst,
-                                      int rank_in, MyList<var> *VarLists, MyList<var> *VarListd)
-{
-  int myrank;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  if (!src || !dst)
-    return false;
-
-  if (src->data->Bg->lev != dst->data->Bg->lev)
-    return false;
-
-  while (src && dst)
-  {
-    if ((src->data->Bg->rank == rank_in) && (dst->data->Bg->rank == myrank))
-    {
-      MyList<var> *varls = VarLists;
-      MyList<var> *varld = VarListd;
-      while (varls && varld)
-      {
-        (void)varls;
-        if (!bssn_gpu_find_device_buffer(dst->data->Bg->fgfs[varld->data->sgfn]))
-          return false;
-        varls = varls->next;
-        varld = varld->next;
-      }
-      if (varls || varld)
-        return false;
-    }
-    src = src->next;
-    dst = dst->next;
-  }
-  return true;
-}
-
-bool parallel_gpu_unpack_segments(const double *data,
-                                  MyList<Parallel::gridseg> *src, MyList<Parallel::gridseg> *dst,
-                                  int rank_in, MyList<var> *VarLists, MyList<var> *VarListd)
-{
-  if (!data)
-    return false;
-
-  int myrank;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  if (!src || !dst)
-    return false;
-  if (src->data->Bg->lev != dst->data->Bg->lev)
-    return false;
-
-  int size_out = 0;
-  while (src && dst)
-  {
-    if ((src->data->Bg->rank == rank_in) && (dst->data->Bg->rank == myrank))
-    {
-      MyList<var> *varls = VarLists;
-      MyList<var> *varld = VarListd;
-      while (varls && varld)
-      {
-        (void)varls;
-        if (bssn_gpu_stage_upload_buffer_to_region(data + size_out,
-                                                   dst->data->Bg->fgfs[varld->data->sgfn],
-                                                   dst->data->Bg->shape,
-                                                   dst->data->Bg->bbox,
-                                                   dst->data->Bg->bbox + dim,
-                                                   dst->data->shape,
-                                                   dst->data->llb))
-          return false;
-        size_out += dst->data->shape[0] * dst->data->shape[1] * dst->data->shape[2];
-        varls = varls->next;
-        varld = varld->next;
-      }
-      if (varls || varld)
-        return false;
-    }
-    src = src->next;
-    dst = dst->next;
-  }
-  return true;
-}
-
-int parallel_var_list_count(MyList<var> *var_list)
-{
-  int count = 0;
-  while (var_list)
-  {
-    count++;
-    var_list = var_list->next;
-  }
-  return count;
-}
-
-void parallel_report_mpi_error(const char *context, int errcode, int req_no)
-{
-  char errstr[MPI_MAX_ERROR_STRING];
-  int len = 0;
-  MPI_Error_string(errcode, errstr, &len);
-  int myrank = -1;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-  fprintf(stderr, "MPI failure on rank %d in %s (transfer_ctx=%s, req_no=%d): %.*s\n",
-          myrank, context, g_parallel_transfer_context, req_no, len, errstr);
-  fflush(stderr);
-}
-
-void parallel_report_mpi_status_errors(const char *context, int req_no,
-                                       int status_count, MPI_Status *stats)
-{
-  if (!stats)
-    return;
-
-  int myrank = -1;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-
-  const int count = (status_count >= 0 && status_count <= req_no) ? status_count : req_no;
-  for (int i = 0; i < count; i++)
-  {
-    if (stats[i].MPI_ERROR != MPI_SUCCESS)
-    {
-      char errstr[MPI_MAX_ERROR_STRING];
-      int len = 0;
-      MPI_Error_string(stats[i].MPI_ERROR, errstr, &len);
-      fprintf(stderr,
-              "MPI request failure on rank %d in %s (transfer_ctx=%s, status_idx=%d, source=%d, tag=%d): %.*s\n",
-              myrank, context, g_parallel_transfer_context, i,
-              stats[i].MPI_SOURCE, stats[i].MPI_TAG, len, errstr);
-    }
-  }
-  fflush(stderr);
-}
-
-int parallel_waitsome_checked(int req_no, MPI_Request *reqs, int *outcount,
-                              int *completed, MPI_Status *stats,
-                              const char *context)
-{
-  MPI_Errhandler old_handler;
-  MPI_Comm_get_errhandler(MPI_COMM_WORLD, &old_handler);
-  MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
-  int rc = MPI_Waitsome(req_no, reqs, outcount, completed, stats);
-  MPI_Comm_set_errhandler(MPI_COMM_WORLD, old_handler);
-  MPI_Errhandler_free(&old_handler);
-  if (rc != MPI_SUCCESS)
-  {
-    parallel_report_mpi_status_errors(context, req_no, outcount ? *outcount : req_no, stats);
-    parallel_report_mpi_error(context, rc, req_no);
-    MPI_Abort(MPI_COMM_WORLD, 1);
-  }
-  return rc;
-}
-
-int parallel_waitall_checked(int req_no, MPI_Request *reqs, MPI_Status *stats,
-                             const char *context)
-{
-  MPI_Errhandler old_handler;
-  MPI_Comm_get_errhandler(MPI_COMM_WORLD, &old_handler);
-  MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
-  int rc = MPI_Waitall(req_no, reqs, stats);
-  MPI_Comm_set_errhandler(MPI_COMM_WORLD, old_handler);
-  MPI_Errhandler_free(&old_handler);
-  if (rc != MPI_SUCCESS)
-  {
-    parallel_report_mpi_status_errors(context, req_no, req_no, stats);
-    parallel_report_mpi_error(context, rc, req_no);
-    MPI_Abort(MPI_COMM_WORLD, 1);
-  }
-  return rc;
-}
-
-} // namespace
 
 int Parallel::partition1(int &nx, int split_size, int min_width, int cpusize, int shape) // special for 1 diemnsion
 {
@@ -4065,28 +3775,14 @@ int Parallel::data_packer(double *data, MyList<Parallel::gridseg> *src, MyList<P
                      dst->data->llb, dst->data->uub);
               break;
             case 2:
-              if (!bssn_cuda_restrict3_pack ||
+              f_restrict3(DIM, dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
-                  bssn_cuda_restrict3_pack(DIM,
+                          src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
-                                           dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
+                          dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry);
-                                           src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
-                                           dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry))
-              {
-                f_restrict3(DIM, dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
-                            src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
-                            dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry);
-              }
               break;
             case 3:
-              if (!bssn_cuda_prolong3_pack ||
+              f_prolong3(DIM, src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
-                  bssn_cuda_prolong3_pack(DIM,
+                         dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
-                                          src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
+                         dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry);
-                                          dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
-                                          dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry))
-              {
-                f_prolong3(DIM, src->data->Bg->bbox, src->data->Bg->bbox + dim, src->data->Bg->shape, src->data->Bg->fgfs[varls->data->sgfn],
-                           dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
-                           dst->data->llb, dst->data->uub, varls->data->SoA, Symmetry);
-              }
             }
           if (dir == UNPACK) // from target data to corresponding grid
             f_copy(DIM, dst->data->Bg->bbox, dst->data->Bg->bbox + dim, dst->data->Bg->shape, dst->data->Bg->fgfs[varld->data->sgfn],
@@ -4204,7 +3900,6 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
   int *completed = new int[2 * cpusize];
   int req_no = 0;
   int pending_recv = 0;
-  const int mpi_tag = parallel_next_transfer_tag();
 
   double **send_data = new double *[cpusize];
   double **rec_data = new double *[cpusize];
@@ -4231,7 +3926,7 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
         cout << "out of memory when new in short transfer, place 1" << endl;
         MPI_Abort(MPI_COMM_WORLD, 1);
       }
-      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 1;
       req_no++;
@@ -4267,18 +3962,18 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
         MPI_Abort(MPI_COMM_WORLD, 1);
       }
       data_packer(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 0;
       req_no++;
     }
   }
+
   // Unpack as soon as receive completes to reduce pure wait time.
   while (pending_recv > 0)
   {
     int outcount = 0;
-    parallel_waitsome_checked(req_no, reqs, &outcount, completed, stats,
+    MPI_Waitsome(req_no, reqs, &outcount, completed, stats);
-                              "Parallel::transfer");
     if (outcount == MPI_UNDEFINED) break;
 
     for (int i = 0; i < outcount; i++)
@@ -4293,7 +3988,7 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
     }
   }
 
-  if (req_no > 0) parallel_waitall_checked(req_no, reqs, stats, "Parallel::transfer");
+  if (req_no > 0) MPI_Waitall(req_no, reqs, stats);
 
   if (rec_data[myrank])
     data_packer(rec_data[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
@@ -4334,7 +4029,6 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
   int *completed = new int[2 * cpusize];
   int req_no = 0;
   int pending_recv = 0;
-  const int mpi_tag = parallel_next_transfer_tag();
 
   double **send_data = new double *[cpusize];
   double **rec_data = new double *[cpusize];
@@ -4361,7 +4055,7 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
         cout << "out of memory when new in short transfer, place 1" << endl;
         MPI_Abort(MPI_COMM_WORLD, 1);
       }
-      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 1;
       req_no++;
@@ -4397,7 +4091,7 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
         MPI_Abort(MPI_COMM_WORLD, 1);
       }
       data_packermix(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, reqs + req_no);
+      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 0;
       req_no++;
@@ -4408,8 +4102,7 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
   while (pending_recv > 0)
   {
     int outcount = 0;
-    parallel_waitsome_checked(req_no, reqs, &outcount, completed, stats,
+    MPI_Waitsome(req_no, reqs, &outcount, completed, stats);
-                              "Parallel::transfermix");
     if (outcount == MPI_UNDEFINED) break;
 
     for (int i = 0; i < outcount; i++)
@@ -4424,7 +4117,7 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
     }
   }
 
-  if (req_no > 0) parallel_waitall_checked(req_no, reqs, stats, "Parallel::transfermix");
+  if (req_no > 0) MPI_Waitall(req_no, reqs, stats);
 
   if (rec_data[myrank])
     data_packermix(rec_data[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
@@ -4448,11 +4141,6 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
   delete[] recv_lengths;
 }
 void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry)
-{
-  Sync(Pat, VarList, Symmetry, "Parallel::Sync(Patch)");
-}
-
-void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry, const char *context)
 {
   int cpusize;
   MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
@@ -4471,10 +4159,7 @@ void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry, const char *
                                                                           // but for transfer_dst[node] the data may locate on any node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-    ParallelTransferContextGuard transfer_ctx(context ? context : "Parallel::Sync(Patch)");
-    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -4494,23 +4179,12 @@ void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry, const char *
 }
 void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
 {
-  Sync(PatL, VarList, Symmetry, "Parallel::Sync(PatchList)");
-}
-
-void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, const char *context)
-{
-  std::string patchlist_context = context ? std::string(context) + " [patchlist]" : "Parallel::Sync(PatchList)";
-
   // Patch inner Synch
   MyList<Patch> *Pp = PatL;
-  int patch_index = 0;
   while (Pp)
   {
-    std::string patch_context = context ? std::string(context) + " [patch " + std::to_string(patch_index) + "]"
+    Sync(Pp->data, VarList, Symmetry);
-                                        : "Parallel::Sync(Patch)";
-    Sync(Pp->data, VarList, Symmetry, patch_context.c_str());
     Pp = Pp->next;
-    patch_index++;
   }
 
   // Patch inter Synch
@@ -4530,10 +4204,7 @@ void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, con
     build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-    ParallelTransferContextGuard transfer_ctx(patchlist_context.c_str());
-    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -4554,11 +4225,6 @@ void Parallel::Sync(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, con
 // Merged Sync: collect all intra-patch and inter-patch grid segment lists,
 // then issue a single transfer() call instead of N+1 separate ones.
 void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry)
-{
-  Sync_merged(PatL, VarList, Symmetry, "Parallel::Sync_merged");
-}
-
-void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, const char *context)
 {
   int cpusize;
   MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
@@ -4636,10 +4302,7 @@ void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmet
     dst_buffer->destroyList();
 
   // Phase C: Single transfer
-  {
+  transfer(combined_src, combined_dst, VarList, VarList, Symmetry);
-    ParallelTransferContextGuard transfer_ctx(context ? context : "Parallel::Sync_merged");
-    transfer(combined_src, combined_dst, VarList, VarList, Symmetry);
-  }
 
   // Phase D: Cleanup
   for (int node = 0; node < cpusize; node++)
@@ -4657,8 +4320,7 @@ Parallel::SyncCache::SyncCache()
     : valid(false), cpusize(0), combined_src(0), combined_dst(0),
       send_lengths(0), recv_lengths(0), send_bufs(0), recv_bufs(0),
       send_buf_caps(0), recv_buf_caps(0), reqs(0), stats(0), max_reqs(0),
-      lengths_valid(false), lengths_var_count(-1),
+      lengths_valid(false), tc_req_node(0), tc_req_is_recv(0), tc_completed(0)
-      tc_req_node(0), tc_req_is_recv(0), tc_completed(0)
 {
 }
 // SyncCache invalidate: free grid segment lists but keep buffers
@@ -4677,7 +4339,6 @@ void Parallel::SyncCache::invalidate()
   }
   valid = false;
   lengths_valid = false;
-  lengths_var_count = -1;
 }
 // SyncCache destroy: free everything
 void Parallel::SyncCache::destroy()
@@ -4708,8 +4369,6 @@ void Parallel::SyncCache::destroy()
   reqs = 0; stats = 0;
   tc_req_node = 0; tc_req_is_recv = 0; tc_completed = 0;
   cpusize = 0; max_reqs = 0;
-  lengths_valid = false;
-  lengths_var_count = -1;
 }
 // transfer_cached: reuse pre-allocated buffers from SyncCache
 void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel::gridseg> **dst,
@@ -4723,9 +4382,6 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
 
   int req_no = 0;
   int pending_recv = 0;
-  const int mpi_tag = parallel_next_transfer_tag();
-  const int current_var_count = parallel_var_list_count(VarList1);
-  const bool lengths_match = cache.lengths_valid && cache.lengths_var_count == current_var_count;
   int node;
   int *req_node = cache.tc_req_node;
   int *req_is_recv = cache.tc_req_is_recv;
@@ -4736,14 +4392,8 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
   {
     if (node == myrank) continue;
 
-    int rlength;
+    int rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
-    if (!lengths_match)
+    cache.recv_lengths[node] = rlength;
-    {
-      rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
-      cache.recv_lengths[node] = rlength;
-    }
-    else
-      rlength = cache.recv_lengths[node];
     if (rlength > 0)
     {
       if (rlength > cache.recv_buf_caps[node])
@@ -4752,7 +4402,7 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
         cache.recv_bufs[node] = new double[rlength];
         cache.recv_buf_caps[node] = rlength;
       }
-      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 1;
       req_no++;
@@ -4761,14 +4411,8 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
   }
 
   // Local transfer on this rank.
-  int self_len;
+  int self_len = data_packer(0, src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
-  if (!lengths_match)
+  cache.recv_lengths[myrank] = self_len;
-  {
-    self_len = data_packer(0, src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
-    cache.recv_lengths[myrank] = self_len;
-  }
-  else
-    self_len = cache.recv_lengths[myrank];
   if (self_len > 0)
   {
     if (self_len > cache.recv_buf_caps[myrank])
@@ -4785,14 +4429,8 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
   {
     if (node == myrank) continue;
 
-    int slength;
+    int slength = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-    if (!lengths_match)
+    cache.send_lengths[node] = slength;
-    {
-      slength = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      cache.send_lengths[node] = slength;
-    }
-    else
-      slength = cache.send_lengths[node];
     if (slength > 0)
     {
       if (slength > cache.send_buf_caps[node])
@@ -4802,22 +4440,18 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
         cache.send_buf_caps[node] = slength;
       }
       data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 0;
       req_no++;
     }
   }
 
-  cache.lengths_valid = true;
-  cache.lengths_var_count = current_var_count;
-
   // Unpack as soon as receive completes to reduce pure wait time.
   while (pending_recv > 0)
   {
     int outcount = 0;
-    parallel_waitsome_checked(req_no, cache.reqs, &outcount, completed, cache.stats,
+    MPI_Waitsome(req_no, cache.reqs, &outcount, completed, cache.stats);
-                              "Parallel::transfer_cached");
     if (outcount == MPI_UNDEFINED) break;
 
     for (int i = 0; i < outcount; i++)
@@ -4832,7 +4466,7 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
     }
   }
 
-  if (req_no > 0) parallel_waitall_checked(req_no, cache.reqs, cache.stats, "Parallel::transfer_cached");
+  if (req_no > 0) MPI_Waitall(req_no, cache.reqs, cache.stats);
 
   if (self_len > 0)
     data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
@@ -5038,11 +4672,8 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
   int myrank;
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
   int cpusize = cache.cpusize;
-  const int current_var_count = parallel_var_list_count(VarList);
-  const bool lengths_match = cache.lengths_valid && cache.lengths_var_count == current_var_count;
   state.req_no = 0;
   state.active = true;
-  state.mpi_tag = parallel_next_transfer_tag();
   state.pending_recv = 0;
   // Allocate tracking arrays
   delete[] state.req_node; delete[] state.req_is_recv;
@@ -5057,7 +4688,7 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
     if (node == myrank)
     {
       int length;
-      if (!lengths_match) {
+      if (!cache.lengths_valid) {
         length = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
         cache.recv_lengths[node] = length;
       } else {
@@ -5071,16 +4702,13 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
           cache.recv_bufs[node] = new double[length];
           cache.recv_buf_caps[node] = length;
         }
-        bssn_gpu_prepare_host_buffer(cache.recv_bufs[node], length);
+        data_packer(cache.recv_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
-        if (!parallel_can_gpu_pack_segments(src[myrank], dst[myrank], node, VarList, VarList) ||
-            !parallel_gpu_pack_segments(cache.recv_bufs[node], src[myrank], dst[myrank], node, VarList, VarList))
-          data_packer(cache.recv_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
       }
     }
     else
     {
       int slength;
-      if (!lengths_match) {
+      if (!cache.lengths_valid) {
         slength = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
         cache.send_lengths[node] = slength;
       } else {
@@ -5094,16 +4722,13 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
           cache.send_bufs[node] = new double[slength];
           cache.send_buf_caps[node] = slength;
         }
-        bssn_gpu_prepare_host_buffer(cache.send_bufs[node], slength);
+        data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
-        if (!parallel_can_gpu_pack_segments(src[myrank], dst[myrank], node, VarList, VarList) ||
-            !parallel_gpu_pack_segments(cache.send_bufs[node], src[myrank], dst[myrank], node, VarList, VarList))
-          data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
         state.req_node[state.req_no] = node;
         state.req_is_recv[state.req_no] = 0;
-        MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, state.mpi_tag, MPI_COMM_WORLD, cache.reqs + state.req_no++);
+        MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
       }
       int rlength;
-      if (!lengths_match) {
+      if (!cache.lengths_valid) {
         rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList, VarList, Symmetry);
         cache.recv_lengths[node] = rlength;
       } else {
@@ -5120,16 +4745,15 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
         state.req_node[state.req_no] = node;
         state.req_is_recv[state.req_no] = 1;
         state.pending_recv++;
-        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, state.mpi_tag, MPI_COMM_WORLD, cache.reqs + state.req_no++);
+        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
       }
     }
   }
   cache.lengths_valid = true;
-  cache.lengths_var_count = current_var_count;
 }
 // Sync_finish: progressive unpack as receives complete, then wait for sends
 void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
-                           MyList<var> *VarList, int Symmetry, bool unpack_to_host)
+                           MyList<var> *VarList, int Symmetry)
 {
   if (!state.active)
     return;
@@ -5141,12 +4765,7 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
 
   // Unpack local data first (no MPI needed)
   if (cache.recv_bufs[myrank] && cache.recv_lengths[myrank] > 0)
-  {
+    data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList, VarList, Symmetry);
-    if (unpack_to_host ||
-        !parallel_can_gpu_unpack_segments(src[myrank], dst[myrank], myrank, VarList, VarList) ||
-        !parallel_gpu_unpack_segments(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, VarList, VarList))
-      data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList, VarList, Symmetry);
-  }
 
   // Progressive unpack of remote receives
   if (state.pending_recv > 0 && state.req_no > 0)
@@ -5156,8 +4775,7 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
     while (pending > 0)
     {
       int outcount = 0;
-      parallel_waitsome_checked(state.req_no, cache.reqs, &outcount, completed, cache.stats,
+      MPI_Waitsome(state.req_no, cache.reqs, &outcount, completed, cache.stats);
-                                "Parallel::Sync_finish");
       if (outcount == MPI_UNDEFINED) break;
       for (int i = 0; i < outcount; i++)
       {
@@ -5165,10 +4783,7 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
         if (idx >= 0 && state.req_is_recv[idx])
         {
           int recv_node = state.req_node[idx];
-          if (unpack_to_host ||
+          data_packer(cache.recv_bufs[recv_node], src[recv_node], dst[recv_node], recv_node, UNPACK, VarList, VarList, Symmetry);
-              !parallel_can_gpu_unpack_segments(src[recv_node], dst[recv_node], recv_node, VarList, VarList) ||
-              !parallel_gpu_unpack_segments(cache.recv_bufs[recv_node], src[recv_node], dst[recv_node], recv_node, VarList, VarList))
-            data_packer(cache.recv_bufs[recv_node], src[recv_node], dst[recv_node], recv_node, UNPACK, VarList, VarList, Symmetry);
           pending--;
         }
       }
@@ -5177,7 +4792,7 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
   }
 
   // Wait for remaining sends
-  if (state.req_no > 0) parallel_waitall_checked(state.req_no, cache.reqs, cache.stats, "Parallel::Sync_finish");
+  if (state.req_no > 0) MPI_Waitall(state.req_no, cache.reqs, cache.stats);
 
   delete[] state.req_node; state.req_node = 0;
   delete[] state.req_is_recv; state.req_is_recv = 0;
@@ -5620,10 +5235,7 @@ void Parallel::PeriodicBD(Patch *Pat, MyList<var> *VarList, int Symmetry)
     build_PhysBD_gstl(Pat, src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-    ParallelTransferContextGuard transfer_ctx("Parallel::PeriodicBD");
-    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -5894,10 +5506,7 @@ void Parallel::Prolong(Patch *Patc, Patch *Patf,
     build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-    ParallelTransferContextGuard transfer_ctx("Parallel::Prolong");
-    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -5956,10 +5565,7 @@ void Parallel::Restrict(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
     build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-    ParallelTransferContextGuard transfer_ctx("Parallel::Restrict");
-    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -6003,10 +5609,7 @@ void Parallel::Restrict_after(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
     build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-    ParallelTransferContextGuard transfer_ctx("Parallel::Restrict_after");
-    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -6051,10 +5654,7 @@ void Parallel::OutBdLow2Hi(Patch *Patc, Patch *Patf,
     build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-    ParallelTransferContextGuard transfer_ctx("Parallel::OutBdLow2Hi(Patch)");
-    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -6110,10 +5710,7 @@ void Parallel::OutBdLow2Hi(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
     build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-    ParallelTransferContextGuard transfer_ctx("Parallel::OutBdLow2Hi(PatchList)");
-    transfer(transfer_src, transfer_dst, VarList1, VarList2, Symmetry);
-  }
 
   if (dst)
     dst->destroyList();
@@ -6388,26 +5985,17 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
 
   int req_no = 0;
   int pending_recv = 0;
-  const int mpi_tag = parallel_next_transfer_tag();
+  int *req_node = new int[cache.max_reqs];
-  const int current_var_count = parallel_var_list_count(VarList1);
+  int *req_is_recv = new int[cache.max_reqs];
-  const bool lengths_match = cache.lengths_valid && cache.lengths_var_count == current_var_count;
+  int *completed = new int[cache.max_reqs];
-  int *req_node = cache.tc_req_node;
-  int *req_is_recv = cache.tc_req_is_recv;
-  int *completed = cache.tc_completed;
 
   // Post receives first so peers can progress rendezvous early.
   for (int node = 0; node < cpusize; node++)
   {
     if (node == myrank) continue;
 
-    int rlength;
+    int rlength = data_packermix(0, cache.combined_src[node], cache.combined_dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
-    if (!lengths_match)
+    cache.recv_lengths[node] = rlength;
-    {
-      rlength = data_packermix(0, cache.combined_src[node], cache.combined_dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
-      cache.recv_lengths[node] = rlength;
-    }
-    else
-      rlength = cache.recv_lengths[node];
     if (rlength > 0)
     {
       if (rlength > cache.recv_buf_caps[node])
@@ -6416,7 +6004,7 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
         cache.recv_bufs[node] = new double[rlength];
         cache.recv_buf_caps[node] = rlength;
       }
-      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 1;
       req_no++;
@@ -6425,14 +6013,8 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
   }
 
   // Local transfer on this rank.
-  int self_len;
+  int self_len = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
-  if (!lengths_match)
+  cache.recv_lengths[myrank] = self_len;
-  {
-    self_len = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
-    cache.recv_lengths[myrank] = self_len;
-  }
-  else
-    self_len = cache.recv_lengths[myrank];
   if (self_len > 0)
   {
     if (self_len > cache.recv_buf_caps[myrank])
@@ -6449,14 +6031,8 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
   {
     if (node == myrank) continue;
 
-    int slength;
+    int slength = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-    if (!lengths_match)
+    cache.send_lengths[node] = slength;
-    {
-      slength = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      cache.send_lengths[node] = slength;
-    }
-    else
-      slength = cache.send_lengths[node];
     if (slength > 0)
     {
       if (slength > cache.send_buf_caps[node])
@@ -6466,22 +6042,18 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
         cache.send_buf_caps[node] = slength;
       }
       data_packermix(cache.send_bufs[node], cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, mpi_tag, MPI_COMM_WORLD, cache.reqs + req_no);
+      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 0;
       req_no++;
     }
   }
 
-  cache.lengths_valid = true;
-  cache.lengths_var_count = current_var_count;
-
   // Unpack as soon as receive completes to reduce pure wait time.
   while (pending_recv > 0)
   {
     int outcount = 0;
-    parallel_waitsome_checked(req_no, cache.reqs, &outcount, completed, cache.stats,
+    MPI_Waitsome(req_no, cache.reqs, &outcount, completed, cache.stats);
-                              "Parallel::transfermix_cached");
     if (outcount == MPI_UNDEFINED) break;
 
     for (int i = 0; i < outcount; i++)
@@ -6496,10 +6068,14 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
     }
   }
 
-  if (req_no > 0) parallel_waitall_checked(req_no, cache.reqs, cache.stats, "Parallel::transfermix_cached");
+  if (req_no > 0) MPI_Waitall(req_no, cache.reqs, cache.stats);
 
   if (self_len > 0)
     data_packermix(cache.recv_bufs[myrank], cache.combined_src[myrank], cache.combined_dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
+
+  delete[] req_node;
+  delete[] req_is_recv;
+  delete[] completed;
 }
 
 // collect all buffer grid segments or blocks for given patch
@@ -7211,10 +6787,7 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
     build_gstl(src[node], dst, &transfer_src[node], &transfer_dst[node]); // for transfer[node], data locate on cpu#node
   }
 
-  {
+  transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-    ParallelTransferContextGuard transfer_ctx("Parallel::Interp copy");
-    transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-  }
 
   for (int node = 0; node < cpusize; node++)
   {
@@ -7253,10 +6826,7 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
         Sync(PatcL, FutureList, Symmetry);
       }
       //<~~~prolong then
-      {
+      transfer(transfer_src, transfer_dst, FutureList, FutureList, Symmetry);
-        ParallelTransferContextGuard transfer_ctx("Parallel::Interp prolong FutureList");
-        transfer(transfer_src, transfer_dst, FutureList, FutureList, Symmetry);
-      }
 
       // for StateList
       // time interpolation part
@@ -7272,10 +6842,7 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
         Sync(PatcL, tmList, Symmetry);
       }
       //<~~~prolong then
-      {
+      transfer(transfer_src, transfer_dst, tmList, StateList, Symmetry);
-        ParallelTransferContextGuard transfer_ctx("Parallel::Interp prolong StateList");
-        transfer(transfer_src, transfer_dst, tmList, StateList, Symmetry);
-      }
     }
     else
     {
@@ -7286,10 +6853,7 @@ void Parallel::fill_level_data(MyList<Patch> *PatLd, MyList<Patch> *PatLs, MyLis
         Sync(PatcL, VarList, Symmetry);
       }
       //<~~~prolong then
-      {
+      transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-        ParallelTransferContextGuard transfer_ctx("Parallel::Interp prolong VarList");
-        transfer(transfer_src, transfer_dst, VarList, VarList, Symmetry);
-      }
     }
 
     for (int node = 0; node < cpusize; node++)

AMSS_NCKU_source/Parallel.h

@@ -90,11 +90,8 @@ namespace Parallel
                    MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /*target */,
                    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;
@@ -111,7 +108,6 @@ namespace Parallel
     MPI_Status *stats;
     int max_reqs;
     bool lengths_valid;
-    int lengths_var_count;
     int *tc_req_node;
     int *tc_req_is_recv;
     int *tc_completed;
@@ -128,17 +124,16 @@ namespace Parallel
   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) {}
+    AsyncSyncState() : req_no(0), active(false), 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, bool unpack_to_host = true);
+                   MyList<var> *VarList, int Symmetry);
   void OutBdLow2Hi(Patch *Patc, Patch *Patf,
                    MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
                    int Symmetry);

AMSS_NCKU_source/bssn_class.C

@@ -15,7 +15,6 @@ using namespace std;
 #endif
 
 #include <time.h>
-#include <unistd.h>
 
 #include "macrodef.h"
 #include "misc.h"
@@ -31,10 +30,6 @@ using namespace std;
 #include "getnp4.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"
@@ -46,12 +41,6 @@ using namespace std;
 #include "derivatives.h"
 #include "ricci_gamma.h"
 
-// Compile-time switch for per-timestep memory usage collection/printing.
-// Default is OFF to reduce overhead in production runs.
-#ifndef BSSN_ENABLE_MEM_USAGE_LOG
-#define BSSN_ENABLE_MEM_USAGE_LOG 0
-#endif
-
 //================================================================================================
 
 // define bssn_class
@@ -749,7 +738,6 @@ void bssn_class::Initialize()
   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];
 }
 
 //================================================================================================
@@ -1026,24 +1014,6 @@ bssn_class::~bssn_class()
       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
@@ -1079,23 +1049,6 @@ bssn_class::~bssn_class()
   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();
-  }
-}
-
 //================================================================================================
 
 
@@ -2080,7 +2033,6 @@ void bssn_class::Evolve(int Steps)
   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
      {
@@ -2177,10 +2129,8 @@ void bssn_class::Evolve(int Steps)
   #endif
   */
   
-#if BSSN_ENABLE_MEM_USAGE_LOG
   perf bssn_perf;
   size_t current_min, current_avg, current_max, peak_min, peak_avg, peak_max;
-#endif
 
   for (int lev = 0; lev < GH->levels; lev++)
     GH->Lt[lev] = PhysTime;
@@ -2265,7 +2215,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);
-    InvalidateSyncCaches();
+    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(); }
 #endif
 
 #if (REGLEV == 0 && (PSTR == 1 || PSTR == 2))
@@ -2274,7 +2224,6 @@ void bssn_class::Evolve(int Steps)
 //		fgt(PhysTime-dT_mon,StartTime,dT_mon/2),ErrorMonitor);
 #endif
 
-#if BSSN_ENABLE_MEM_USAGE_LOG
     // Retrieve memory usage information used during computation; master process prints it
     bssn_perf.MemoryUsage(&current_min, &current_avg, &current_max,
                           &peak_min, &peak_avg, &peak_max, nprocs);
@@ -2290,7 +2239,6 @@ void bssn_class::Evolve(int Steps)
              (double)peak_max / (1024.0 * 1024.0));
       cout << endl;
     }
-#endif
     
     // Output puncture positions at each step
     if (myrank == 0)
@@ -2349,9 +2297,6 @@ void bssn_class::Evolve(int Steps)
 #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
@@ -2486,7 +2431,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))
-  InvalidateSyncCaches();
+  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(); }
 #endif
 }
 
@@ -2665,7 +2610,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))
-  InvalidateSyncCaches();
+  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(); }
 #endif
 }
 
@@ -2832,7 +2777,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))
-        InvalidateSyncCaches();
+        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(); }
 
         //               a_stream.clear();
         //               a_stream.str("");
@@ -2847,7 +2792,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))
-      InvalidateSyncCaches();
+      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(); }
 
       //               a_stream.clear();
       //               a_stream.str("");
@@ -2866,7 +2811,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))
-          InvalidateSyncCaches();
+          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(); }
 
           //               a_stream.clear();
           //               a_stream.str("");
@@ -2882,7 +2827,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))
-          InvalidateSyncCaches();
+          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(); }
 
           //               a_stream.clear();
           //               a_stream.str("");
@@ -3073,11 +3018,6 @@ void bssn_class::RecursiveStep(int lev, int num) // in all 2^(lev+1)-1 steps
 #if 1
 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));
@@ -5808,6 +5748,12 @@ void bssn_class::SHStep()
 // 0: do not use mixing two levels data for OutBD; 1: do use
 
 #define MIXOUTB 0
+// In the cached Restrict->OutBdLow2Hi path, coarse Sync is usually redundant:
+// OutBdLow2Hi_cached reads coarse owned cells (build_owned_gsl type-4), not coarse ghost/buffer cells.
+// Keep a switch to restore the old behavior if needed for debugging.
+#ifndef RP_SYNC_COARSE_AFTER_RESTRICT
+#define RP_SYNC_COARSE_AFTER_RESTRICT 0
+#endif
 void bssn_class::RestrictProlong(int lev, int YN, bool BB,
                                  MyList<var> *SL, MyList<var> *OL, MyList<var> *corL)
 // we assume
@@ -5871,7 +5817,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 //       misc::tillherecheck(GH->Commlev[GH->mylev],GH->start_rank[GH->mylev],a_stream.str());
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (PSTR == 1 || PSTR == 2)
 //       a_stream.clear();
@@ -5922,7 +5870,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 //       misc::tillherecheck(GH->Commlev[GH->mylev],GH->start_rank[GH->mylev],a_stream.str());
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (PSTR == 1 || PSTR == 2)
 //       a_stream.clear();
@@ -6008,7 +5958,9 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -6030,7 +5982,9 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -6095,7 +6049,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, SynchList_pre, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -6119,7 +6075,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -6196,13 +6154,16 @@ void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
 #else
       Parallel::Restrict_after(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, Symmetry);
 #endif
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
     }
 
     Parallel::Sync_cached(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_rp_fine[lev]);
   }
 }
 #undef MIXOUTB
+#undef RP_SYNC_COARSE_AFTER_RESTRICT
 
 //================================================================================================
 
@@ -6298,7 +6259,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_cached(GH->PatL[lev], DG_List, Symmetry, sync_cache_psi4[lev]);
+  Parallel::Sync(GH->PatL[lev], DG_List, Symmetry);
 #endif
 
 #ifdef WithShell
@@ -7325,9 +7286,6 @@ void bssn_class::Constraint_Out()
             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], 
@@ -7362,7 +7320,6 @@ void bssn_class::Constraint_Out()
                                  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;
@@ -7371,7 +7328,7 @@ void bssn_class::Constraint_Out()
           Pp = Pp->next;
         }
       }
-      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry, "bssn_class::Constraint_Out[level]");
+      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry);
     }
 #ifdef WithShell
     if (0) // if the constrait quantities can be reused from the step rhs calculation
@@ -7593,7 +7550,7 @@ void bssn_class::AH_Prepare_derivatives()
       }
       Pp = Pp->next;
     }
-    Parallel::Sync(GH->PatL[lev], AHDList, Symmetry, "bssn_class::AH_Prepare_derivatives");
+    Parallel::Sync(GH->PatL[lev], AHDList, Symmetry);
   }
 }
 
@@ -7832,9 +7789,6 @@ void bssn_class::Interp_Constraint(bool infg)
             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], 
@@ -7869,7 +7823,6 @@ void bssn_class::Interp_Constraint(bool infg)
                                  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;
@@ -7878,7 +7831,7 @@ void bssn_class::Interp_Constraint(bool infg)
           Pp = Pp->next;
         }
       }
-      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry, "bssn_class::Interp_Constraint[level]");
+      Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry);
     }
 #ifdef WithShell
     if (0) // if the constrait quantities can be reused from the step rhs calculation
@@ -8136,7 +8089,7 @@ void bssn_class::Compute_Constraint()
         Pp = Pp->next;
       }
     }
-    Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry, "bssn_class::Compute_Constraint[level]");
+    Parallel::Sync(GH->PatL[lev], ConstraintList, Symmetry);
   }
   // prolong restrict constraint quantities
   for (lev = GH->levels - 1; lev > 0; lev--)
@@ -8451,12 +8404,6 @@ void bssn_class::Enforce_algcon(int lev, int fg)
 
 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;
 
@@ -8470,11 +8417,8 @@ bool bssn_class::check_Stdin_Abort()
     timeout.tv_usec = 0;
 
     int activity = select(STDIN_FILENO + 1, &readfds, nullptr, nullptr, &timeout);
-    if (activity <= 0) {
-        return false;
-    }
 
-    if (FD_ISSET(STDIN_FILENO, &readfds)) {
+    if (activity > 0 && FD_ISSET(STDIN_FILENO, &readfds)) {
         string input_abort;
         if (cin >> input_abort) {
             if (input_abort == "stop") {

AMSS_NCKU_source/bssn_class.h

@@ -132,7 +132,6 @@ public:
        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;
@@ -177,12 +176,8 @@ public:
        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();

AMSS_NCKU_source/bssn_cuda_ops.h

@@ -1,68 +0,0 @@
-#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 *phi_field,
-                               const double *lap_field,
-                               const double *boundary_src,
-                               double *stage_data,
-                               double *rhs_accum,
-                               double propspeed,
-                               const double SoA[3],
-                               int symmetry,
-                               int lev,
-                               int rk_stage,
-                               bool force_host_boundary_fix,
-                               bool download_to_host = true);
-
-int bssn_cuda_rk4_boundary_batch(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,
-                                 int symmetry,
-                                 const double *const *state0_list,
-                                 double *const *stage_data_list,
-                                 double *const *rhs_accum_list,
-                                 int num_var,
-                                 int rk_stage,
-                                 bool download_to_host = false);
-
-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);
-void bssn_cuda_release_rk4_caches();
-void bssn_cuda_release_interp_caches();
-
-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_restrict3_pack(int wei,
-                             const double *llbc, const double *uubc, const int *extc, double *func,
-                             const double *llbf, const double *uubf, const int *extf, const double *funf,
-                             const double *llbr, const double *uubr,
-                             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

@@ -1,936 +0,0 @@
-#include "macrodef.h"
-
-#ifdef USE_GPU
-
-#include <algorithm>
-#include <cmath>
-#include <cstring>
-#include <cstdlib>
-#include <iomanip>
-#include <vector>
-
-#include "bssn_class.h"
-#include "bssn_cuda_ops.h"
-#include "bssn_gpu.h"
-#include "bssn_macro.h"
-
-namespace
-{
-enum StageProfileMetric
-{
-  STAGE_PROFILE_TOTAL = 0,
-  STAGE_PROFILE_RHS,
-  STAGE_PROFILE_RUN_STAGE,
-  STAGE_PROFILE_RUN_STAGE_DEVICE,
-  STAGE_PROFILE_RUN_STAGE_HOST_FIX,
-  STAGE_PROFILE_LOWERBOUND,
-  STAGE_PROFILE_ENSURE,
-  STAGE_PROFILE_DOWNLOAD,
-  STAGE_PROFILE_CLEAR_CACHE,
-  STAGE_PROFILE_SYNC_START,
-  STAGE_PROFILE_SYNC_FINISH,
-  STAGE_PROFILE_REFRESH,
-  STAGE_PROFILE_COUNT
-};
-
-static const int kStageProfileMaxLevels = 32;
-
-struct StageProfileStore
-{
-  bool env_checked;
-  bool enabled;
-  int calls[kStageProfileMaxLevels];
-  double metric[kStageProfileMaxLevels][STAGE_PROFILE_COUNT];
-};
-
-StageProfileStore &stage_profile_store()
-{
-  static StageProfileStore store = {};
-  return store;
-}
-
-bool stage_profile_enabled()
-{
-  StageProfileStore &store = stage_profile_store();
-  if (!store.env_checked)
-  {
-    const char *env = getenv("AMSS_GPU_STAGE_TIMING");
-    store.enabled = (env && env[0] && strcmp(env, "0") != 0);
-    store.env_checked = true;
-  }
-  return store.enabled;
-}
-
-void stage_profile_note_call(int lev)
-{
-  if (lev >= 0 && lev < kStageProfileMaxLevels)
-    stage_profile_store().calls[lev]++;
-}
-
-void stage_profile_add(int lev, StageProfileMetric metric, double seconds)
-{
-  if (lev >= 0 && lev < kStageProfileMaxLevels)
-    stage_profile_store().metric[lev][metric] += seconds;
-}
-
-const char *stage_profile_metric_name(StageProfileMetric metric)
-{
-  switch (metric)
-  {
-  case STAGE_PROFILE_TOTAL:
-    return "total";
-  case STAGE_PROFILE_RHS:
-    return "rhs";
-  case STAGE_PROFILE_RUN_STAGE:
-    return "run_stage";
-  case STAGE_PROFILE_RUN_STAGE_DEVICE:
-    return "run_stage_dev";
-  case STAGE_PROFILE_RUN_STAGE_HOST_FIX:
-    return "run_stage_host";
-  case STAGE_PROFILE_LOWERBOUND:
-    return "lower";
-  case STAGE_PROFILE_ENSURE:
-    return "ensure";
-  case STAGE_PROFILE_DOWNLOAD:
-    return "download";
-  case STAGE_PROFILE_CLEAR_CACHE:
-    return "clear_cache";
-  case STAGE_PROFILE_SYNC_START:
-    return "sync_start";
-  case STAGE_PROFILE_SYNC_FINISH:
-    return "sync_finish";
-  case STAGE_PROFILE_REFRESH:
-    return "refresh";
-  default:
-    return "unknown";
-  }
-}
-} // namespace
-
-void bssn_cuda_dump_stage_profile()
-{
-  if (!stage_profile_enabled())
-    return;
-
-  int myrank = 0;
-  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-
-  StageProfileStore &store = stage_profile_store();
-  int global_calls_sum[kStageProfileMaxLevels] = {};
-  double global_metric_sum[kStageProfileMaxLevels][STAGE_PROFILE_COUNT] = {};
-  double global_metric_max[kStageProfileMaxLevels][STAGE_PROFILE_COUNT] = {};
-
-  MPI_Reduce(store.calls, global_calls_sum, kStageProfileMaxLevels, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
-  MPI_Reduce(store.metric[0], global_metric_sum[0],
-             kStageProfileMaxLevels * STAGE_PROFILE_COUNT,
-             MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
-  MPI_Reduce(store.metric[0], global_metric_max[0],
-             kStageProfileMaxLevels * STAGE_PROFILE_COUNT,
-             MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
-
-  if (myrank != 0)
-    return;
-
-  cout << endl;
-  cout << " GPU stage timing summary (sum/max over MPI ranks) " << endl;
-  cout << " lev  calls";
-  for (int metric = 0; metric < STAGE_PROFILE_COUNT; ++metric)
-    cout << "  " << setw(22) << stage_profile_metric_name(static_cast<StageProfileMetric>(metric));
-  cout << endl;
-
-  for (int lev = 0; lev < kStageProfileMaxLevels; ++lev)
-  {
-    if (global_calls_sum[lev] == 0)
-      continue;
-
-    cout << setw(4) << lev << "  " << setw(5) << global_calls_sum[lev];
-    for (int metric = 0; metric < STAGE_PROFILE_COUNT; ++metric)
-    {
-      cout << "  "
-           << setw(10) << setprecision(6) << fixed << global_metric_sum[lev][metric]
-           << "/"
-           << setw(10) << setprecision(6) << fixed << global_metric_max[lev][metric];
-    }
-    cout << endl;
-  }
-  cout << endl;
-}
-
-void bssn_class::Step_MainPath_GPU(int lev, int YN)
-{
-#ifdef WithShell
-#error "Step_MainPath_GPU currently supports Patch grids only."
-#endif
-
-  const bool profile_enabled = stage_profile_enabled();
-  const double step_total_begin = profile_enabled ? MPI_Wtime() : 0.0;
-  if (profile_enabled)
-    stage_profile_note_call(lev);
-
-  if (bssn_gpu_bind_process_device(myrank))
-  {
-    cerr << "GPU device bind failure on MPI rank " << myrank << endl;
-    MPI_Abort(MPI_COMM_WORLD, 1);
-  }
-  if (profile_enabled)
-  {
-    const double t0 = MPI_Wtime();
-    bssn_gpu_clear_cached_device_buffers();
-    stage_profile_add(lev, STAGE_PROFILE_CLEAR_CACHE, MPI_Wtime() - t0);
-  }
-  else
-    bssn_gpu_clear_cached_device_buffers();
-
-  setpbh(BH_num, Porg0, Mass, BH_num_input);
-
-  const double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
-
-#if (MAPBH == 1)
-  if (BH_num > 0 && lev == GH->levels - 1)
-  {
-    compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
-    for (int ithBH = 0; ithBH < BH_num; ithBH++)
-    {
-      for (int ith = 0; ith < 3; ith++)
-        Porg1[ithBH][ith] = Porg0[ithBH][ith] + Porg_rhs[ithBH][ith] * dT_lev;
-      if (Symmetry > 0)
-        Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
-      if (Symmetry == 2)
-      {
-        Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
-        Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
-      }
-    }
-  }
-
-  if (lev == a_lev)
-    AnalysisStuff(lev, dT_lev);
-#endif
-
-#ifdef With_AHF
-  AH_Step_Find(lev, dT_lev);
-#endif
-
-  const bool BB = fgt(PhysTime, StartTime, dT_lev / 2);
-  (void)BB;
-  double ndeps = (lev < GH->movls) ? numepsb : numepss;
-  double TRK4 = PhysTime;
-  int iter_count = 0;
-  int pre = 0, cor = 1;
-  int ERROR = 0;
-  const bool keep_stage_sync_on_device = (RPS == 1) && (MAPBH == 1) && (REGLEV == 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;
-        std::vector<const double *> batch_state0;
-        std::vector<double *> batch_stage;
-        std::vector<double *> batch_rhs;
-
-        while (varl0)
-        {
-          const bool force_host_boundary_fix = false;
-          const bool can_batch_device_path = (lev > 0) && !force_host_boundary_fix;
-          if (can_batch_device_path)
-          {
-            batch_state0.push_back(cg->fgfs[varl0->data->sgfn]);
-            batch_stage.push_back(cg->fgfs[varls->data->sgfn]);
-            batch_rhs.push_back(cg->fgfs[varlr->data->sgfn]);
-            varl0 = varl0->next;
-            varlb = varlb->next;
-            varls = varls->next;
-            varlr = varlr->next;
-            continue;
-          }
-
-          const double var_begin = profile_enabled ? MPI_Wtime() : 0.0;
-          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[phi0->sgfn],
-                                         cg->fgfs[Lap0->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,
-                                         force_host_boundary_fix, 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;
-          }
-          if (profile_enabled)
-          {
-            stage_profile_add(lev,
-                              force_host_boundary_fix ? STAGE_PROFILE_RUN_STAGE_HOST_FIX
-                                                      : STAGE_PROFILE_RUN_STAGE_DEVICE,
-                              MPI_Wtime() - var_begin);
-          }
-          varl0 = varl0->next;
-          varlb = varlb->next;
-          varls = varls->next;
-          varlr = varlr->next;
-        }
-
-        if (!ERROR && !batch_state0.empty())
-        {
-          const double batch_begin = profile_enabled ? MPI_Wtime() : 0.0;
-          if (bssn_cuda_rk4_boundary_batch(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],
-                                           Symmetry,
-                                           &batch_state0[0],
-                                           &batch_stage[0],
-                                           &batch_rhs[0],
-                                           static_cast<int>(batch_state0.size()),
-                                           rk_stage, false))
-          {
-            cerr << "GPU rk4/boundary batch failure: lev=" << lev
-                 << " rk_stage=" << rk_stage
-                 << " vars=" << batch_state0.size()
-                 << " bbox=(" << cg->bbox[0] << ":" << cg->bbox[3] << ","
-                 << cg->bbox[1] << ":" << cg->bbox[4] << ","
-                 << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
-            ERROR = 1;
-          }
-          else if (profile_enabled)
-          {
-            stage_profile_add(lev, STAGE_PROFILE_RUN_STAGE_DEVICE, MPI_Wtime() - batch_begin);
-          }
-        }
-      };
-
-  auto stage_download_var_list =
-      [&](Block *cg, MyList<var> *var_list, bool skip_unmapped) {
-        std::vector<double *> batch_host_ptrs;
-        std::vector<MyList<var> *> batch_vars;
-        while (var_list)
-        {
-          double *host_ptr = cg->fgfs[var_list->data->sgfn];
-          if (skip_unmapped && !bssn_gpu_find_device_buffer(host_ptr))
-          {
-            var_list = var_list->next;
-            continue;
-          }
-          batch_host_ptrs.push_back(host_ptr);
-          batch_vars.push_back(var_list);
-          var_list = var_list->next;
-        }
-        if (!batch_host_ptrs.empty() &&
-            bssn_gpu_download_buffer_batch(cg->shape, &batch_host_ptrs[0],
-                                           static_cast<int>(batch_host_ptrs.size())))
-        {
-          for (size_t i = 0; i < batch_host_ptrs.size(); ++i)
-          {
-            if (bssn_cuda_download_buffer(cg->shape, batch_host_ptrs[i]))
-            {
-              cerr << "GPU stage download failure: lev=" << lev
-                   << " var=" << batch_vars[i]->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;
-            }
-          }
-        }
-      };
-
-  auto stage_download_patch_list =
-      [&](MyList<var> *var_list, bool skip_unmapped) {
-        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)
-              stage_download_var_list(cg, var_list, skip_unmapped);
-
-            if (block_it == patch_it->data->ble)
-              break;
-            block_it = block_it->next;
-          }
-          if (ERROR)
-            break;
-          patch_it = patch_it->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)
-      {
-        double t0 = 0.0;
-        if (profile_enabled)
-          t0 = MPI_Wtime();
-        if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_FIRST_TIME))
-          ERROR = 1;
-        if (profile_enabled)
-          stage_profile_add(lev, STAGE_PROFILE_RHS, MPI_Wtime() - t0);
-
-        if (profile_enabled)
-          t0 = MPI_Wtime();
-        run_stage_on_block(cg, Pp->data, StateList, StateList, SynchList_pre, RHSList, iter_count);
-        if (profile_enabled)
-          stage_profile_add(lev, STAGE_PROFILE_RUN_STAGE, MPI_Wtime() - t0);
-
-        if (profile_enabled)
-          t0 = MPI_Wtime();
-        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 (profile_enabled)
-          stage_profile_add(lev, STAGE_PROFILE_LOWERBOUND, MPI_Wtime() - t0);
-      }
-      if (BP == Pp->data->ble)
-        break;
-      BP = BP->next;
-    }
-    Pp = Pp->next;
-  }
-
-  if (!ERROR)
-  {
-    if (!keep_stage_sync_on_device)
-    {
-      double t0 = 0.0;
-      if (profile_enabled)
-        t0 = MPI_Wtime();
-      stage_download_patch_list(SynchList_pre, false);
-      if (profile_enabled)
-        stage_profile_add(lev, STAGE_PROFILE_DOWNLOAD, MPI_Wtime() - t0);
-      if (!ERROR)
-      {
-        if (profile_enabled)
-          t0 = MPI_Wtime();
-        bssn_gpu_clear_cached_device_buffers();
-        if (profile_enabled)
-          stage_profile_add(lev, STAGE_PROFILE_CLEAR_CACHE, MPI_Wtime() - t0);
-      }
-    }
-  }
-
-  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;
-  if (profile_enabled)
-  {
-    const double t0 = MPI_Wtime();
-    Parallel::Sync_start(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev], async_pre);
-    stage_profile_add(lev, STAGE_PROFILE_SYNC_START, MPI_Wtime() - t0);
-  }
-  else
-    Parallel::Sync_start(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev], async_pre);
-  if (profile_enabled)
-  {
-    const double t0 = MPI_Wtime();
-    Parallel::Sync_finish(sync_cache_pre[lev], async_pre, SynchList_pre, Symmetry,
-                          !keep_stage_sync_on_device);
-    stage_profile_add(lev, STAGE_PROFILE_SYNC_FINISH, MPI_Wtime() - t0);
-  }
-  else
-    Parallel::Sync_finish(sync_cache_pre[lev], async_pre, SynchList_pre, Symmetry,
-                          !keep_stage_sync_on_device);
-  if (!ERROR && !keep_stage_sync_on_device)
-  {
-    if (profile_enabled)
-    {
-      const double t0 = MPI_Wtime();
-      refresh_stage_device_after_sync(SynchList_pre, sync_cache_pre[lev]);
-      stage_profile_add(lev, STAGE_PROFILE_REFRESH, MPI_Wtime() - t0);
-    }
-    else
-      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)
-        {
-          double t0 = 0.0;
-          if (profile_enabled)
-            t0 = MPI_Wtime();
-          ensure_stage_device_var_list(cg, SynchList_pre);
-          if (profile_enabled)
-            stage_profile_add(lev, STAGE_PROFILE_ENSURE, MPI_Wtime() - t0);
-
-          if (profile_enabled)
-            t0 = MPI_Wtime();
-          if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_THEN))
-            ERROR = 1;
-          if (profile_enabled)
-            stage_profile_add(lev, STAGE_PROFILE_RHS, MPI_Wtime() - t0);
-
-          if (profile_enabled)
-            t0 = MPI_Wtime();
-          run_stage_on_block(cg, Pp->data, StateList, SynchList_pre, SynchList_cor, RHSList, iter_count);
-          if (profile_enabled)
-            stage_profile_add(lev, STAGE_PROFILE_RUN_STAGE, MPI_Wtime() - t0);
-
-          if (profile_enabled)
-            t0 = MPI_Wtime();
-          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 (profile_enabled)
-            stage_profile_add(lev, STAGE_PROFILE_LOWERBOUND, MPI_Wtime() - t0);
-        }
-
-        if (BP == Pp->data->ble)
-          break;
-        BP = BP->next;
-      }
-      Pp = Pp->next;
-    }
-
-    if (!ERROR)
-    {
-      if (!keep_stage_sync_on_device)
-      {
-        double t0 = 0.0;
-        if (profile_enabled)
-          t0 = MPI_Wtime();
-        stage_download_patch_list(SynchList_cor, false);
-        if (profile_enabled)
-          stage_profile_add(lev, STAGE_PROFILE_DOWNLOAD, MPI_Wtime() - t0);
-        if (!ERROR)
-        {
-          if (profile_enabled)
-            t0 = MPI_Wtime();
-          bssn_gpu_clear_cached_device_buffers();
-          if (profile_enabled)
-            stage_profile_add(lev, STAGE_PROFILE_CLEAR_CACHE, MPI_Wtime() - t0);
-        }
-      }
-    }
-
-    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;
-    if (profile_enabled)
-    {
-      const double t0 = MPI_Wtime();
-      Parallel::Sync_start(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev], async_cor);
-      stage_profile_add(lev, STAGE_PROFILE_SYNC_START, MPI_Wtime() - t0);
-    }
-    else
-      Parallel::Sync_start(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev], async_cor);
-    if (profile_enabled)
-    {
-      const double t0 = MPI_Wtime();
-      Parallel::Sync_finish(sync_cache_cor[lev], async_cor, SynchList_cor, Symmetry,
-                            !keep_stage_sync_on_device);
-      stage_profile_add(lev, STAGE_PROFILE_SYNC_FINISH, MPI_Wtime() - t0);
-    }
-    else
-      Parallel::Sync_finish(sync_cache_cor[lev], async_cor, SynchList_cor, Symmetry,
-                            !keep_stage_sync_on_device);
-    if (!ERROR && !keep_stage_sync_on_device && iter_count < 3)
-    {
-      if (profile_enabled)
-      {
-        const double t0 = MPI_Wtime();
-        refresh_stage_device_after_sync(SynchList_cor, sync_cache_cor[lev]);
-        stage_profile_add(lev, STAGE_PROFILE_REFRESH, MPI_Wtime() - t0);
-      }
-      else
-        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
-
-  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 (!ERROR && keep_stage_sync_on_device)
-  {
-    // After the swaps above, only StateList points at arrays updated during this step.
-    // OldStateList/SynchList_cor remain valid on host because their backing arrays were
-    // read-only during the RK step, and SynchList_pre is reused only as scratch later.
-    const double t0 = profile_enabled ? MPI_Wtime() : 0.0;
-    stage_download_patch_list(StateList, true);
-    if (profile_enabled)
-      stage_profile_add(lev, STAGE_PROFILE_DOWNLOAD, MPI_Wtime() - t0);
-  }
-
-  if (profile_enabled)
-  {
-    const double t0 = MPI_Wtime();
-    bssn_gpu_clear_cached_device_buffers();
-    stage_profile_add(lev, STAGE_PROFILE_CLEAR_CACHE, MPI_Wtime() - t0);
-  }
-  else
-    bssn_gpu_clear_cached_device_buffers();
-
-  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];
-    }
-  }
-
-  if (profile_enabled)
-    stage_profile_add(lev, STAGE_PROFILE_TOTAL, MPI_Wtime() - step_total_begin);
-}
-
-#endif

AMSS_NCKU_source/bssn_gpu.h

@@ -4,6 +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
 
@@ -65,42 +67,6 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,
 
 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);
-int bssn_gpu_download_buffer_batch(const int *ex, double **host_ptrs, int num_buffers);
-
 /** Init GPU side data in GPUMeta. */
 // void init_fluid_meta_gpu(GPUMeta *gpu_meta);
 

AMSS_NCKU_source/bssn_macro.h

@@ -68,7 +68,6 @@ if(TIME_COUNT_EACH_RANK == 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.f90

@@ -62,7 +62,6 @@
   real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gmx_Res, Gmy_Res, Gmz_Res
 !  gont = 0: success; gont = 1: something wrong
   integer::gont
-  integer :: i,j,k
 
 !~~~~~~> Other variables:
 
@@ -86,13 +85,6 @@
 
   real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA
   real*8            :: dX, dY, dZ, PI
-  real*8            :: divb_loc,det_loc
-  real*8            :: gupxx_loc,gupxy_loc,gupxz_loc,gupyy_loc,gupyz_loc,gupzz_loc
-  real*8            :: Rxx_loc,Rxy_loc,Rxz_loc,Ryy_loc,Ryz_loc,Rzz_loc
-  real*8            :: fxx_loc,fxy_loc,fxz_loc
-  real*8            :: Gamxa_loc,Gamya_loc,Gamza_loc
-  real*8            :: f_loc,chin_loc
-  real*8            :: l_fxx,l_fxy,l_fxz,l_fyy,l_fyz,l_fzz,S_loc
   real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0
   real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0
   real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0
@@ -105,7 +97,7 @@
 #endif
 
 #if (GAUGE == 6 || GAUGE == 7)
-  integer :: BHN
+  integer :: BHN,i,j,k
   real*8, dimension(9) :: Porg
   real*8, dimension(3) :: Mass
   real*8 :: r1,r2,M,A,w1,w2,C1,C2
@@ -153,24 +145,22 @@
   dY = Y(2) - Y(1)
   dZ = Z(2) - Z(1)
 
-  do k=1,ex(3)
+  alpn1 = Lap + ONE
-  do j=1,ex(2)
+  chin1 = chi + ONE
-  do i=1,ex(1)
+  gxx = dxx + ONE
-    alpn1(i,j,k) = Lap(i,j,k) + ONE
+  gyy = dyy + ONE
-    chin1(i,j,k) = chi(i,j,k) + ONE
+  gzz = dzz + ONE
-    gxx(i,j,k) = dxx(i,j,k) + ONE
-    gyy(i,j,k) = dyy(i,j,k) + ONE
-    gzz(i,j,k) = dzz(i,j,k) + ONE
-  enddo
-  enddo
-  enddo
 
   call fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev)
   call fderivs(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev)
   call fderivs(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev)
   
+  div_beta = betaxx + betayy + betazz
+ 
   call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
 
+  chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi
+
   call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
   call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev)
   call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev)
@@ -178,179 +168,151 @@
   call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev)
   call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
 
-  do k=1,ex(3)
+  gxx_rhs = - TWO * alpn1 * Axx    -  F2o3 * gxx * div_beta          + &
-  do j=1,ex(2)
+              TWO *(  gxx * betaxx +   gxy * betayx +   gxz * betazx)
-  do i=1,ex(1)
-    divb_loc = betaxx(i,j,k) + betayy(i,j,k) + betazz(i,j,k)
-    div_beta(i,j,k) = divb_loc
 
-    chi_rhs(i,j,k) = F2o3 * chin1(i,j,k) * (alpn1(i,j,k) * trK(i,j,k) - divb_loc)
+  gyy_rhs = - TWO * alpn1 * Ayy    -  F2o3 * gyy * div_beta          + &
+              TWO *(  gxy * betaxy +   gyy * betayy +   gyz * betazy)
 
-    gxx_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axx(i,j,k) - F2o3 * gxx(i,j,k) * divb_loc + &
+  gzz_rhs = - TWO * alpn1 * Azz    -  F2o3 * gzz * div_beta          + &
-         TWO * ( gxx(i,j,k) * betaxx(i,j,k) + gxy(i,j,k) * betayx(i,j,k) + gxz(i,j,k) * betazx(i,j,k) )
+              TWO *(  gxz * betaxz +   gyz * betayz +   gzz * betazz)
 
-    gyy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayy(i,j,k) - F2o3 * gyy(i,j,k) * divb_loc + &
+  gxy_rhs = - TWO * alpn1 * Axy    +  F1o3 * gxy    * div_beta       + &
-         TWO * ( gxy(i,j,k) * betaxy(i,j,k) + gyy(i,j,k) * betayy(i,j,k) + gyz(i,j,k) * betazy(i,j,k) )
+                      gxx * betaxy                  +   gxz * betazy + &
+                                       gyy * betayx +   gyz * betazx   &
+                                                    -   gxy * betazz
 
-    gzz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Azz(i,j,k) - F2o3 * gzz(i,j,k) * divb_loc + &
+  gyz_rhs = - TWO * alpn1 * Ayz    +  F1o3 * gyz    * div_beta       + &
-         TWO * ( gxz(i,j,k) * betaxz(i,j,k) + gyz(i,j,k) * betayz(i,j,k) + gzz(i,j,k) * betazz(i,j,k) )
+                      gxy * betaxz +   gyy * betayz                  + &
+                      gxz * betaxy                  +   gzz * betazy   &
+                                                    -   gyz * betaxx
  
-    gxy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axy(i,j,k) + F1o3 * gxy(i,j,k) * divb_loc + &
+  gxz_rhs = - TWO * alpn1 * Axz    +  F1o3 * gxz    * div_beta       + &
-         gxx(i,j,k) * betaxy(i,j,k) + gxz(i,j,k) * betazy(i,j,k) + gyy(i,j,k) * betayx(i,j,k) + &
+                      gxx * betaxz +   gxy * betayz                  + &
-         gyz(i,j,k) * betazx(i,j,k) - gxy(i,j,k) * betazz(i,j,k)
+                                       gyz * betayx +   gzz * betazx   &
+                                                    -   gxz * betayy     !rhs for gij
 
-    gyz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayz(i,j,k) + F1o3 * gyz(i,j,k) * divb_loc + &
+! invert tilted metric
-         gxy(i,j,k) * betaxz(i,j,k) + gyy(i,j,k) * betayz(i,j,k) + gxz(i,j,k) * betaxy(i,j,k) + &
+  gupzz =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-         gzz(i,j,k) * betazy(i,j,k) - gyz(i,j,k) * betaxx(i,j,k)
+           gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
+  gupxx =   ( gyy * gzz - gyz * gyz ) / gupzz
+  gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz
+  gupxz =   ( gxy * gyz - gyy * gxz ) / gupzz
+  gupyy =   ( gxx * gzz - gxz * gxz ) / gupzz
+  gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
+  gupzz =   ( gxx * gyy - gxy * gxy ) / gupzz
 
-    gxz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axz(i,j,k) + F1o3 * gxz(i,j,k) * divb_loc + &
+  if(co == 0)then
-         gxx(i,j,k) * betaxz(i,j,k) + gxy(i,j,k) * betayz(i,j,k) + gyz(i,j,k) * betayx(i,j,k) + &
+! Gam^i_Res = Gam^i + gup^ij_,j
-         gzz(i,j,k) * betazx(i,j,k) - gxz(i,j,k) * betayy(i,j,k)
+  Gmx_Res = Gamx - (gupxx*(gupxx*gxxx+gupxy*gxyx+gupxz*gxzx)&
+                   +gupxy*(gupxx*gxyx+gupxy*gyyx+gupxz*gyzx)&
+                   +gupxz*(gupxx*gxzx+gupxy*gyzx+gupxz*gzzx)&
+                   +gupxx*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)&
+                   +gupxy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)&
+                   +gupxz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)&
+                   +gupxx*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
+                   +gupxy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
+                   +gupxz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
+  Gmy_Res = Gamy - (gupxx*(gupxy*gxxx+gupyy*gxyx+gupyz*gxzx)&
+                   +gupxy*(gupxy*gxyx+gupyy*gyyx+gupyz*gyzx)&
+                   +gupxz*(gupxy*gxzx+gupyy*gyzx+gupyz*gzzx)&
+                   +gupxy*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)&
+                   +gupyy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)&
+                   +gupyz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)&
+                   +gupxy*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
+                   +gupyy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
+                   +gupyz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
+  Gmz_Res = Gamz - (gupxx*(gupxz*gxxx+gupyz*gxyx+gupzz*gxzx)&
+                   +gupxy*(gupxz*gxyx+gupyz*gyyx+gupzz*gyzx)&
+                   +gupxz*(gupxz*gxzx+gupyz*gyzx+gupzz*gzzx)&
+                   +gupxy*(gupxz*gxxy+gupyz*gxyy+gupzz*gxzy)&
+                   +gupyy*(gupxz*gxyy+gupyz*gyyy+gupzz*gyzy)&
+                   +gupyz*(gupxz*gxzy+gupyz*gyzy+gupzz*gzzy)&
+                   +gupxz*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
+                   +gupyz*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
+                   +gupzz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
+  endif
 
-    det_loc = gxx(i,j,k) * gyy(i,j,k) * gzz(i,j,k) + gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) + &
+! second kind of connection
-         gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k) * gxz(i,j,k) - &
+  Gamxxx =HALF*( gupxx*gxxx + gupxy*(TWO*gxyx - gxxy ) + gupxz*(TWO*gxzx - gxxz ))
-         gxy(i,j,k) * gxy(i,j,k) * gzz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) * gyz(i,j,k)
+  Gamyxx =HALF*( gupxy*gxxx + gupyy*(TWO*gxyx - gxxy ) + gupyz*(TWO*gxzx - gxxz ))
-    gupxx_loc = ( gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gyz(i,j,k) ) / det_loc
+  Gamzxx =HALF*( gupxz*gxxx + gupyz*(TWO*gxyx - gxxy ) + gupzz*(TWO*gxzx - gxxz ))
-    gupxy_loc = - ( gxy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupxz_loc = ( gxy(i,j,k) * gyz(i,j,k) - gyy(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupyy_loc = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupyz_loc = - ( gxx(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupzz_loc = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k) * gxy(i,j,k) ) / det_loc
-    gupxx(i,j,k) = gupxx_loc
-    gupxy(i,j,k) = gupxy_loc
-    gupxz(i,j,k) = gupxz_loc
-    gupyy(i,j,k) = gupyy_loc
-    gupyz(i,j,k) = gupyz_loc
-    gupzz(i,j,k) = gupzz_loc
  
-    if(co == 0)then
+  Gamxyy =HALF*( gupxx*(TWO*gxyy - gyyx ) + gupxy*gyyy + gupxz*(TWO*gyzy - gyyz ))
-      Gmx_Res(i,j,k) = Gamx(i,j,k) - ( &
+  Gamyyy =HALF*( gupxy*(TWO*gxyy - gyyx ) + gupyy*gyyy + gupyz*(TWO*gyzy - gyyz ))
-           gupxx_loc*(gupxx_loc*gxxx(i,j,k)+gupxy_loc*gxyx(i,j,k)+gupxz_loc*gxzx(i,j,k)) + &
+  Gamzyy =HALF*( gupxz*(TWO*gxyy - gyyx ) + gupyz*gyyy + gupzz*(TWO*gyzy - gyyz ))
-           gupxy_loc*(gupxx_loc*gxyx(i,j,k)+gupxy_loc*gyyx(i,j,k)+gupxz_loc*gyzx(i,j,k)) + &
-           gupxz_loc*(gupxx_loc*gxzx(i,j,k)+gupxy_loc*gyzx(i,j,k)+gupxz_loc*gzzx(i,j,k)) + &
-           gupxx_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + &
-           gupxy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + &
-           gupxz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + &
-           gupxx_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
-           gupxy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
-           gupxz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
-      Gmy_Res(i,j,k) = Gamy(i,j,k) - ( &
-           gupxx_loc*(gupxy_loc*gxxx(i,j,k)+gupyy_loc*gxyx(i,j,k)+gupyz_loc*gxzx(i,j,k)) + &
-           gupxy_loc*(gupxy_loc*gxyx(i,j,k)+gupyy_loc*gyyx(i,j,k)+gupyz_loc*gyzx(i,j,k)) + &
-           gupxz_loc*(gupxy_loc*gxzx(i,j,k)+gupyy_loc*gyzx(i,j,k)+gupyz_loc*gzzx(i,j,k)) + &
-           gupxy_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + &
-           gupyy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + &
-           gupyz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + &
-           gupxy_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
-           gupyy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
-           gupyz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
-      Gmz_Res(i,j,k) = Gamz(i,j,k) - ( &
-           gupxx_loc*(gupxz_loc*gxxx(i,j,k)+gupyz_loc*gxyx(i,j,k)+gupzz_loc*gxzx(i,j,k)) + &
-           gupxy_loc*(gupxz_loc*gxyx(i,j,k)+gupyz_loc*gyyx(i,j,k)+gupzz_loc*gyzx(i,j,k)) + &
-           gupxz_loc*(gupxz_loc*gxzx(i,j,k)+gupyz_loc*gyzx(i,j,k)+gupzz_loc*gzzx(i,j,k)) + &
-           gupxy_loc*(gupxz_loc*gxxy(i,j,k)+gupyz_loc*gxyy(i,j,k)+gupzz_loc*gxzy(i,j,k)) + &
-           gupyy_loc*(gupxz_loc*gxyy(i,j,k)+gupyz_loc*gyyy(i,j,k)+gupzz_loc*gyzy(i,j,k)) + &
-           gupyz_loc*(gupxz_loc*gxzy(i,j,k)+gupyz_loc*gyzy(i,j,k)+gupzz_loc*gzzy(i,j,k)) + &
-           gupxz_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
-           gupyz_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
-           gupzz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
-    endif
 
-    Gamxxx(i,j,k)=HALF*( gupxx_loc*gxxx(i,j,k) + gupxy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupxz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+  Gamxzz =HALF*( gupxx*(TWO*gxzz - gzzx ) + gupxy*(TWO*gyzz - gzzy ) + gupxz*gzzz)
-    Gamyxx(i,j,k)=HALF*( gupxy_loc*gxxx(i,j,k) + gupyy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupyz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+  Gamyzz =HALF*( gupxy*(TWO*gxzz - gzzx ) + gupyy*(TWO*gyzz - gzzy ) + gupyz*gzzz)
-    Gamzxx(i,j,k)=HALF*( gupxz_loc*gxxx(i,j,k) + gupyz_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupzz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+  Gamzzz =HALF*( gupxz*(TWO*gxzz - gzzx ) + gupyz*(TWO*gyzz - gzzy ) + gupzz*gzzz)
 
-    Gamxyy(i,j,k)=HALF*( gupxx_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupxy_loc*gyyy(i,j,k) + gupxz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+  Gamxxy =HALF*( gupxx*gxxy + gupxy*gyyx + gupxz*( gxzy + gyzx - gxyz ) )
-    Gamyyy(i,j,k)=HALF*( gupxy_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyy_loc*gyyy(i,j,k) + gupyz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+  Gamyxy =HALF*( gupxy*gxxy + gupyy*gyyx + gupyz*( gxzy + gyzx - gxyz ) )
-    Gamzyy(i,j,k)=HALF*( gupxz_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyz_loc*gyyy(i,j,k) + gupzz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+  Gamzxy =HALF*( gupxz*gxxy + gupyz*gyyx + gupzz*( gxzy + gyzx - gxyz ) )
 
-    Gamxzz(i,j,k)=HALF*( gupxx_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupxy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupxz_loc*gzzz(i,j,k))
+  Gamxxz =HALF*( gupxx*gxxz + gupxy*( gxyz + gyzx - gxzy ) + gupxz*gzzx )
-    Gamyzz(i,j,k)=HALF*( gupxy_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupyz_loc*gzzz(i,j,k))
+  Gamyxz =HALF*( gupxy*gxxz + gupyy*( gxyz + gyzx - gxzy ) + gupyz*gzzx )
-    Gamzzz(i,j,k)=HALF*( gupxz_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyz_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupzz_loc*gzzz(i,j,k))
+  Gamzxz =HALF*( gupxz*gxxz + gupyz*( gxyz + gyzx - gxzy ) + gupzz*gzzx )
 
-    Gamxxy(i,j,k)=HALF*( gupxx_loc*gxxy(i,j,k) + gupxy_loc*gyyx(i,j,k) + gupxz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+  Gamxyz =HALF*( gupxx*( gxyz + gxzy - gyzx ) + gupxy*gyyz + gupxz*gzzy )
-    Gamyxy(i,j,k)=HALF*( gupxy_loc*gxxy(i,j,k) + gupyy_loc*gyyx(i,j,k) + gupyz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+  Gamyyz =HALF*( gupxy*( gxyz + gxzy - gyzx ) + gupyy*gyyz + gupyz*gzzy )
-    Gamzxy(i,j,k)=HALF*( gupxz_loc*gxxy(i,j,k) + gupyz_loc*gyyx(i,j,k) + gupzz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+  Gamzyz =HALF*( gupxz*( gxyz + gxzy - gyzx ) + gupyz*gyyz + gupzz*gzzy )
-
-    Gamxxz(i,j,k)=HALF*( gupxx_loc*gxxz(i,j,k) + gupxy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupxz_loc*gzzx(i,j,k) )
-    Gamyxz(i,j,k)=HALF*( gupxy_loc*gxxz(i,j,k) + gupyy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupyz_loc*gzzx(i,j,k) )
-    Gamzxz(i,j,k)=HALF*( gupxz_loc*gxxz(i,j,k) + gupyz_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupzz_loc*gzzx(i,j,k) )
-
-    Gamxyz(i,j,k)=HALF*( gupxx_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupxy_loc*gyyz(i,j,k) + gupxz_loc*gzzy(i,j,k) )
-    Gamyyz(i,j,k)=HALF*( gupxy_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyy_loc*gyyz(i,j,k) + gupyz_loc*gzzy(i,j,k) )
-    Gamzyz(i,j,k)=HALF*( gupxz_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyz_loc*gyyz(i,j,k) + gupzz_loc*gzzy(i,j,k) )
-  enddo
-  enddo
-  enddo
 ! Raise indices of \tilde A_{ij} and store in R_ij
 
+  Rxx =    gupxx * gupxx * Axx + gupxy * gupxy * Ayy + gupxz * gupxz * Azz + &
+      TWO*(gupxx * gupxy * Axy + gupxx * gupxz * Axz + gupxy * gupxz * Ayz)
+
+  Ryy =    gupxy * gupxy * Axx + gupyy * gupyy * Ayy + gupyz * gupyz * Azz + &
+      TWO*(gupxy * gupyy * Axy + gupxy * gupyz * Axz + gupyy * gupyz * Ayz)
+
+  Rzz =    gupxz * gupxz * Axx + gupyz * gupyz * Ayy + gupzz * gupzz * Azz + &
+      TWO*(gupxz * gupyz * Axy + gupxz * gupzz * Axz + gupyz * gupzz * Ayz)
+
+  Rxy =    gupxx * gupxy * Axx + gupxy * gupyy * Ayy + gupxz * gupyz * Azz + &
+          (gupxx * gupyy       + gupxy * gupxy)* Axy                       + &
+          (gupxx * gupyz       + gupxz * gupxy)* Axz                       + &
+          (gupxy * gupyz       + gupxz * gupyy)* Ayz
+
+  Rxz =    gupxx * gupxz * Axx + gupxy * gupyz * Ayy + gupxz * gupzz * Azz + &
+          (gupxx * gupyz       + gupxy * gupxz)* Axy                       + &
+          (gupxx * gupzz       + gupxz * gupxz)* Axz                       + &
+          (gupxy * gupzz       + gupxz * gupyz)* Ayz
+
+  Ryz =    gupxy * gupxz * Axx + gupyy * gupyz * Ayy + gupyz * gupzz * Azz + &
+          (gupxy * gupyz       + gupyy * gupxz)* Axy                       + &
+          (gupxy * gupzz       + gupyz * gupxz)* Axz                       + &
+          (gupyy * gupzz       + gupyz * gupyz)* Ayz
+
 ! Right hand side for Gam^i without shift terms...
   call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
   call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-    gupxx_loc = gupxx(i,j,k)
-    gupxy_loc = gupxy(i,j,k)
-    gupxz_loc = gupxz(i,j,k)
-    gupyy_loc = gupyy(i,j,k)
-    gupyz_loc = gupyz(i,j,k)
-    gupzz_loc = gupzz(i,j,k)
 
-    Rxx_loc = gupxx_loc * gupxx_loc * Axx(i,j,k) + gupxy_loc * gupxy_loc * Ayy(i,j,k) + gupxz_loc * gupxz_loc * Azz(i,j,k) + &
+   Gamx_rhs = - TWO * (   Lapx * Rxx +   Lapy * Rxy +   Lapz * Rxz ) + &
-         TWO * (gupxx_loc * gupxy_loc * Axy(i,j,k) + gupxx_loc * gupxz_loc * Axz(i,j,k) + gupxy_loc * gupxz_loc * Ayz(i,j,k))
+        TWO * alpn1 * (                                                &
-    Ryy_loc = gupxy_loc * gupxy_loc * Axx(i,j,k) + gupyy_loc * gupyy_loc * Ayy(i,j,k) + gupyz_loc * gupyz_loc * Azz(i,j,k) + &
+        -F3o2/chin1 * (   chix * Rxx +   chiy * Rxy +   chiz * Rxz ) - &
-         TWO * (gupxy_loc * gupyy_loc * Axy(i,j,k) + gupxy_loc * gupyz_loc * Axz(i,j,k) + gupyy_loc * gupyz_loc * Ayz(i,j,k))
+              gupxx * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-    Rzz_loc = gupxz_loc * gupxz_loc * Axx(i,j,k) + gupyz_loc * gupyz_loc * Ayy(i,j,k) + gupzz_loc * gupzz_loc * Azz(i,j,k) + &
+              gupxy * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-         TWO * (gupxz_loc * gupyz_loc * Axy(i,j,k) + gupxz_loc * gupzz_loc * Axz(i,j,k) + gupyz_loc * gupzz_loc * Ayz(i,j,k))
+              gupxz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-    Rxy_loc = gupxx_loc * gupxy_loc * Axx(i,j,k) + gupxy_loc * gupyy_loc * Ayy(i,j,k) + gupxz_loc * gupyz_loc * Azz(i,j,k) + &
+                        Gamxxx * Rxx + Gamxyy * Ryy + Gamxzz * Rzz   + &
-         (gupxx_loc * gupyy_loc + gupxy_loc * gupxy_loc) * Axy(i,j,k) + &
+                TWO * ( Gamxxy * Rxy + Gamxxz * Rxz + Gamxyz * Ryz ) )
-         (gupxx_loc * gupyz_loc + gupxz_loc * gupxy_loc) * Axz(i,j,k) + &
-         (gupxy_loc * gupyz_loc + gupxz_loc * gupyy_loc) * Ayz(i,j,k)
-    Rxz_loc = gupxx_loc * gupxz_loc * Axx(i,j,k) + gupxy_loc * gupyz_loc * Ayy(i,j,k) + gupxz_loc * gupzz_loc * Azz(i,j,k) + &
-         (gupxx_loc * gupyz_loc + gupxy_loc * gupxz_loc) * Axy(i,j,k) + &
-         (gupxx_loc * gupzz_loc + gupxz_loc * gupxz_loc) * Axz(i,j,k) + &
-         (gupxy_loc * gupzz_loc + gupxz_loc * gupyz_loc) * Ayz(i,j,k)
-    Ryz_loc = gupxy_loc * gupxz_loc * Axx(i,j,k) + gupyy_loc * gupyz_loc * Ayy(i,j,k) + gupyz_loc * gupzz_loc * Azz(i,j,k) + &
-         (gupxy_loc * gupyz_loc + gupyy_loc * gupxz_loc) * Axy(i,j,k) + &
-         (gupxy_loc * gupzz_loc + gupyz_loc * gupxz_loc) * Axz(i,j,k) + &
-         (gupyy_loc * gupzz_loc + gupyz_loc * gupyz_loc) * Ayz(i,j,k)
-    Rxx(i,j,k) = Rxx_loc
-    Ryy(i,j,k) = Ryy_loc
-    Rzz(i,j,k) = Rzz_loc
-    Rxy(i,j,k) = Rxy_loc
-    Rxz(i,j,k) = Rxz_loc
-    Ryz(i,j,k) = Ryz_loc
 
-    Gamx_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxx_loc + Lapy(i,j,k) * Rxy_loc + Lapz(i,j,k) * Rxz_loc) + &
+   Gamy_rhs = - TWO * (   Lapx * Rxy +   Lapy * Ryy +   Lapz * Ryz ) + &
-         TWO * alpn1(i,j,k) * ( &
+        TWO * alpn1 * (                                                &
-         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxx_loc + chiy(i,j,k) * Rxy_loc + chiz(i,j,k) * Rxz_loc) - &
+        -F3o2/chin1 * (   chix * Rxy +  chiy * Ryy +    chiz * Ryz ) - &
-         gupxx_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
+              gupxy * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-         gupxy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
+              gupyy * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-         gupxz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
+              gupyz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-         Gamxxx(i,j,k) * Rxx_loc + Gamxyy(i,j,k) * Ryy_loc + Gamxzz(i,j,k) * Rzz_loc + &
+                        Gamyxx * Rxx + Gamyyy * Ryy + Gamyzz * Rzz   + &
-         TWO * (Gamxxy(i,j,k) * Rxy_loc + Gamxxz(i,j,k) * Rxz_loc + Gamxyz(i,j,k) * Ryz_loc))
+                TWO * ( Gamyxy * Rxy + Gamyxz * Rxz + Gamyyz * Ryz ) )
 
-    Gamy_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxy_loc + Lapy(i,j,k) * Ryy_loc + Lapz(i,j,k) * Ryz_loc) + &
+   Gamz_rhs = - TWO * (   Lapx * Rxz +   Lapy * Ryz +   Lapz * Rzz ) + &
-         TWO * alpn1(i,j,k) * ( &
+        TWO * alpn1 * (                                                &
-         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxy_loc + chiy(i,j,k) * Ryy_loc + chiz(i,j,k) * Ryz_loc) - &
+        -F3o2/chin1 * (   chix * Rxz +  chiy * Ryz +    chiz * Rzz ) - &
-         gupxy_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
+              gupxz * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-         gupyy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
+              gupyz * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-         gupyz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
+              gupzz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-         Gamyxx(i,j,k) * Rxx_loc + Gamyyy(i,j,k) * Ryy_loc + Gamyzz(i,j,k) * Rzz_loc + &
+                        Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz   + &
-         TWO * (Gamyxy(i,j,k) * Rxy_loc + Gamyxz(i,j,k) * Rxz_loc + Gamyyz(i,j,k) * Ryz_loc))
+                TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) )
-
-    Gamz_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxz_loc + Lapy(i,j,k) * Ryz_loc + Lapz(i,j,k) * Rzz_loc) + &
-         TWO * alpn1(i,j,k) * ( &
-         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxz_loc + chiy(i,j,k) * Ryz_loc + chiz(i,j,k) * Rzz_loc) - &
-         gupxz_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
-         gupyz_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
-         gupzz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
-         Gamzxx(i,j,k) * Rxx_loc + Gamzyy(i,j,k) * Ryy_loc + Gamzzz(i,j,k) * Rzz_loc + &
-         TWO * (Gamzxy(i,j,k) * Rxy_loc + Gamzxz(i,j,k) * Rxz_loc + Gamzyz(i,j,k) * Ryz_loc))
-  enddo
-  enddo
-  enddo
 
   call fdderivs(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,&
                 X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev)
@@ -359,54 +321,38 @@
   call fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,&
                 X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev)
 
+  fxx = gxxx + gxyy + gxzz
+  fxy = gxyx + gyyy + gyzz
+  fxz = gxzx + gyzy + gzzz
+
+  Gamxa =       gupxx * Gamxxx + gupyy * Gamxyy + gupzz * Gamxzz + &
+          TWO*( gupxy * Gamxxy + gupxz * Gamxxz + gupyz * Gamxyz )
+  Gamya =       gupxx * Gamyxx + gupyy * Gamyyy + gupzz * Gamyzz + &
+          TWO*( gupxy * Gamyxy + gupxz * Gamyxz + gupyz * Gamyyz )
+  Gamza =       gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz + &
+          TWO*( gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz )
+
   call fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev)
   call fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev)
   call fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev)
-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-    divb_loc = div_beta(i,j,k)
-    fxx_loc = gxxx(i,j,k) + gxyy(i,j,k) + gxzz(i,j,k)
-    fxy_loc = gxyx(i,j,k) + gyyy(i,j,k) + gyzz(i,j,k)
-    fxz_loc = gxzx(i,j,k) + gyzy(i,j,k) + gzzz(i,j,k)
 
-    gupxx_loc = gupxx(i,j,k)
+  Gamx_rhs =               Gamx_rhs +  F2o3 *  Gamxa * div_beta        - &
-    gupxy_loc = gupxy(i,j,k)
+                     Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz  + &
-    gupxz_loc = gupxz(i,j,k)
+             F1o3 * (gupxx * fxx    + gupxy * fxy    + gupxz * fxz    ) + &
-    gupyy_loc = gupyy(i,j,k)
+                     gupxx * gxxx   + gupyy * gyyx   + gupzz * gzzx    + &
-    gupyz_loc = gupyz(i,j,k)
+              TWO * (gupxy * gxyx   + gupxz * gxzx   + gupyz * gyzx  )
-    gupzz_loc = gupzz(i,j,k)
 
-    Gamxa_loc = gupxx_loc * Gamxxx(i,j,k) + gupyy_loc * Gamxyy(i,j,k) + gupzz_loc * Gamxzz(i,j,k) + &
+  Gamy_rhs =               Gamy_rhs +  F2o3 *  Gamya * div_beta        - &
-         TWO * (gupxy_loc * Gamxxy(i,j,k) + gupxz_loc * Gamxxz(i,j,k) + gupyz_loc * Gamxyz(i,j,k))
+                     Gamxa * betayx - Gamya * betayy - Gamza * betayz  + &
-    Gamya_loc = gupxx_loc * Gamyxx(i,j,k) + gupyy_loc * Gamyyy(i,j,k) + gupzz_loc * Gamyzz(i,j,k) + &
+             F1o3 * (gupxy * fxx    + gupyy * fxy    + gupyz * fxz    ) + &
-         TWO * (gupxy_loc * Gamyxy(i,j,k) + gupxz_loc * Gamyxz(i,j,k) + gupyz_loc * Gamyyz(i,j,k))
+                     gupxx * gxxy   + gupyy * gyyy   + gupzz * gzzy    + &
-    Gamza_loc = gupxx_loc * Gamzxx(i,j,k) + gupyy_loc * Gamzyy(i,j,k) + gupzz_loc * Gamzzz(i,j,k) + &
+              TWO * (gupxy * gxyy   + gupxz * gxzy   + gupyz * gyzy  )
-         TWO * (gupxy_loc * Gamzxy(i,j,k) + gupxz_loc * Gamzxz(i,j,k) + gupyz_loc * Gamzyz(i,j,k))
-    Gamxa(i,j,k) = Gamxa_loc
-    Gamya(i,j,k) = Gamya_loc
-    Gamza(i,j,k) = Gamza_loc
 
-    Gamx_rhs(i,j,k) = Gamx_rhs(i,j,k) + F2o3 * Gamxa_loc * divb_loc - &
+  Gamz_rhs =               Gamz_rhs +  F2o3 *  Gamza * div_beta        - &
-         Gamxa_loc * betaxx(i,j,k) - Gamya_loc * betaxy(i,j,k) - Gamza_loc * betaxz(i,j,k) + &
+                     Gamxa * betazx - Gamya * betazy - Gamza * betazz  + &
-         F1o3 * (gupxx_loc * fxx_loc + gupxy_loc * fxy_loc + gupxz_loc * fxz_loc) + &
+             F1o3 * (gupxz * fxx    + gupyz * fxy    + gupzz * fxz    ) + &
-         gupxx_loc * gxxx(i,j,k) + gupyy_loc * gyyx(i,j,k) + gupzz_loc * gzzx(i,j,k) + &
+                     gupxx * gxxz   + gupyy * gyyz   + gupzz * gzzz    + &
-         TWO * (gupxy_loc * gxyx(i,j,k) + gupxz_loc * gxzx(i,j,k) + gupyz_loc * gyzx(i,j,k))
+              TWO * (gupxy * gxyz   + gupxz * gxzz   + gupyz * gyzz  )    !rhs for Gam^i
-
-    Gamy_rhs(i,j,k) = Gamy_rhs(i,j,k) + F2o3 * Gamya_loc * divb_loc - &
-         Gamxa_loc * betayx(i,j,k) - Gamya_loc * betayy(i,j,k) - Gamza_loc * betayz(i,j,k) + &
-         F1o3 * (gupxy_loc * fxx_loc + gupyy_loc * fxy_loc + gupyz_loc * fxz_loc) + &
-         gupxx_loc * gxxy(i,j,k) + gupyy_loc * gyyy(i,j,k) + gupzz_loc * gzzy(i,j,k) + &
-         TWO * (gupxy_loc * gxyy(i,j,k) + gupxz_loc * gxzy(i,j,k) + gupyz_loc * gyzy(i,j,k))
-
-    Gamz_rhs(i,j,k) = Gamz_rhs(i,j,k) + F2o3 * Gamza_loc * divb_loc - &
-         Gamxa_loc * betazx(i,j,k) - Gamya_loc * betazy(i,j,k) - Gamza_loc * betazz(i,j,k) + &
-         F1o3 * (gupxz_loc * fxx_loc + gupyz_loc * fxy_loc + gupzz_loc * fxz_loc) + &
-         gupxx_loc * gxxz(i,j,k) + gupyy_loc * gyyz(i,j,k) + gupzz_loc * gzzz(i,j,k) + &
-         TWO * (gupxy_loc * gxyz(i,j,k) + gupxz_loc * gxzz(i,j,k) + gupyz_loc * gyzz(i,j,k))
-  enddo
-  enddo
-  enddo
 
 !first kind of connection stored in gij,k
   gxxx = gxx * Gamxxx + gxy * Gamyxx + gxz * Gamzxx
@@ -658,187 +604,189 @@
 !covariant second derivative of chi respect to tilted metric
   call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
 
-  do k=1,ex(3)
+  fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz
-  do j=1,ex(2)
+  fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz
-  do i=1,ex(1)
+  fxz = fxz - Gamxxz * chix - Gamyxz * chiy - Gamzxz * chiz
-    fxx(i,j,k) = fxx(i,j,k) - Gamxxx(i,j,k) * chix(i,j,k) - Gamyxx(i,j,k) * chiy(i,j,k) - Gamzxx(i,j,k) * chiz(i,j,k)
+  fyy = fyy - Gamxyy * chix - Gamyyy * chiy - Gamzyy * chiz
-    fxy(i,j,k) = fxy(i,j,k) - Gamxxy(i,j,k) * chix(i,j,k) - Gamyxy(i,j,k) * chiy(i,j,k) - Gamzxy(i,j,k) * chiz(i,j,k)
+  fyz = fyz - Gamxyz * chix - Gamyyz * chiy - Gamzyz * chiz
-    fxz(i,j,k) = fxz(i,j,k) - Gamxxz(i,j,k) * chix(i,j,k) - Gamyxz(i,j,k) * chiy(i,j,k) - Gamzxz(i,j,k) * chiz(i,j,k)
+  fzz = fzz - Gamxzz * chix - Gamyzz * chiy - Gamzzz * chiz
-    fyy(i,j,k) = fyy(i,j,k) - Gamxyy(i,j,k) * chix(i,j,k) - Gamyyy(i,j,k) * chiy(i,j,k) - Gamzyy(i,j,k) * chiz(i,j,k)
+! Store D^l D_l chi - 3/(2*chi) D^l chi D_l chi in f
-    fyz(i,j,k) = fyz(i,j,k) - Gamxyz(i,j,k) * chix(i,j,k) - Gamyyz(i,j,k) * chiy(i,j,k) - Gamzyz(i,j,k) * chiz(i,j,k)
-    fzz(i,j,k) = fzz(i,j,k) - Gamxzz(i,j,k) * chix(i,j,k) - Gamyzz(i,j,k) * chiy(i,j,k) - Gamzzz(i,j,k) * chiz(i,j,k)
 
-    chin_loc = chin1(i,j,k)
+  f =        gupxx * ( fxx - F3o2/chin1 * chix * chix ) + &
-    f_loc = gupxx(i,j,k) * (fxx(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chix(i,j,k)) + &
+             gupyy * ( fyy - F3o2/chin1 * chiy * chiy ) + &
-            gupyy(i,j,k) * (fyy(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiy(i,j,k)) + &
+             gupzz * ( fzz - F3o2/chin1 * chiz * chiz ) + &
-            gupzz(i,j,k) * (fzz(i,j,k) - F3o2/chin_loc * chiz(i,j,k) * chiz(i,j,k)) + &
+       TWO * gupxy * ( fxy - F3o2/chin1 * chix * chiy ) + &
-            TWO * gupxy(i,j,k) * (fxy(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiy(i,j,k)) + &
+       TWO * gupxz * ( fxz - F3o2/chin1 * chix * chiz ) + &
-            TWO * gupxz(i,j,k) * (fxz(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiz(i,j,k)) + &
+       TWO * gupyz * ( fyz - F3o2/chin1 * chiy * chiz ) 
-            TWO * gupyz(i,j,k) * (fyz(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiz(i,j,k))
+! Add chi part to Ricci tensor:
-    f(i,j,k) = f_loc
 
-    Rxx(i,j,k) = Rxx(i,j,k) + (fxx(i,j,k) - chix(i,j,k)*chix(i,j,k)/chin_loc/TWO + gxx(i,j,k) * f_loc)/chin_loc/TWO
+  Rxx = Rxx + (fxx - chix*chix/chin1/TWO + gxx * f)/chin1/TWO
-    Ryy(i,j,k) = Ryy(i,j,k) + (fyy(i,j,k) - chiy(i,j,k)*chiy(i,j,k)/chin_loc/TWO + gyy(i,j,k) * f_loc)/chin_loc/TWO
+  Ryy = Ryy + (fyy - chiy*chiy/chin1/TWO + gyy * f)/chin1/TWO
-    Rzz(i,j,k) = Rzz(i,j,k) + (fzz(i,j,k) - chiz(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gzz(i,j,k) * f_loc)/chin_loc/TWO
+  Rzz = Rzz + (fzz - chiz*chiz/chin1/TWO + gzz * f)/chin1/TWO
-    Rxy(i,j,k) = Rxy(i,j,k) + (fxy(i,j,k) - chix(i,j,k)*chiy(i,j,k)/chin_loc/TWO + gxy(i,j,k) * f_loc)/chin_loc/TWO
+  Rxy = Rxy + (fxy - chix*chiy/chin1/TWO + gxy * f)/chin1/TWO
-    Rxz(i,j,k) = Rxz(i,j,k) + (fxz(i,j,k) - chix(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gxz(i,j,k) * f_loc)/chin_loc/TWO
+  Rxz = Rxz + (fxz - chix*chiz/chin1/TWO + gxz * f)/chin1/TWO
-    Ryz(i,j,k) = Ryz(i,j,k) + (fyz(i,j,k) - chiy(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gyz(i,j,k) * f_loc)/chin_loc/TWO
+  Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO
-  enddo
-  enddo
-  enddo
 
 ! covariant second derivatives of the lapse respect to physical metric
   call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
                 SYM,SYM,SYM,symmetry,Lev)
 
-  do k=1,ex(3)
+  gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1
-  do j=1,ex(2)
+  gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1
-  do i=1,ex(1)
+  gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1
-    chin_loc = chin1(i,j,k)
+! now get physical second kind of connection
-    gxxx(i,j,k) = (gupxx(i,j,k) * chix(i,j,k) + gupxy(i,j,k) * chiy(i,j,k) + gupxz(i,j,k) * chiz(i,j,k)) / chin_loc
+  Gamxxx = Gamxxx - ( (chix + chix)/chin1 - gxx * gxxx )*HALF
-    gxxy(i,j,k) = (gupxy(i,j,k) * chix(i,j,k) + gupyy(i,j,k) * chiy(i,j,k) + gupyz(i,j,k) * chiz(i,j,k)) / chin_loc
+  Gamyxx = Gamyxx - (                     - gxx * gxxy )*HALF
-    gxxz(i,j,k) = (gupxz(i,j,k) * chix(i,j,k) + gupyz(i,j,k) * chiy(i,j,k) + gupzz(i,j,k) * chiz(i,j,k)) / chin_loc
+  Gamzxx = Gamzxx - (                     - gxx * gxxz )*HALF
+  Gamxyy = Gamxyy - (                     - gyy * gxxx )*HALF
+  Gamyyy = Gamyyy - ( (chiy + chiy)/chin1 - gyy * gxxy )*HALF
+  Gamzyy = Gamzyy - (                     - gyy * gxxz )*HALF
+  Gamxzz = Gamxzz - (                     - gzz * gxxx )*HALF
+  Gamyzz = Gamyzz - (                     - gzz * gxxy )*HALF
+  Gamzzz = Gamzzz - ( (chiz + chiz)/chin1 - gzz * gxxz )*HALF
+  Gamxxy = Gamxxy - (  chiy        /chin1 - gxy * gxxx )*HALF
+  Gamyxy = Gamyxy - (         chix /chin1 - gxy * gxxy )*HALF
+  Gamzxy = Gamzxy - (                     - gxy * gxxz )*HALF
+  Gamxxz = Gamxxz - (  chiz        /chin1 - gxz * gxxx )*HALF
+  Gamyxz = Gamyxz - (                     - gxz * gxxy )*HALF
+  Gamzxz = Gamzxz - (         chix /chin1 - gxz * gxxz )*HALF
+  Gamxyz = Gamxyz - (                     - gyz * gxxx )*HALF
+  Gamyyz = Gamyyz - (  chiz        /chin1 - gyz * gxxy )*HALF
+  Gamzyz = Gamzyz - (         chiy /chin1 - gyz * gxxz )*HALF
 
-    Gamxxx(i,j,k) = Gamxxx(i,j,k) - ( (chix(i,j,k) + chix(i,j,k))/chin_loc - gxx(i,j,k) * gxxx(i,j,k) )*HALF
+  fxx = fxx - Gamxxx*Lapx - Gamyxx*Lapy - Gamzxx*Lapz
-    Gamyxx(i,j,k) = Gamyxx(i,j,k) - (                                   - gxx(i,j,k) * gxxy(i,j,k) )*HALF
+  fyy = fyy - Gamxyy*Lapx - Gamyyy*Lapy - Gamzyy*Lapz
-    Gamzxx(i,j,k) = Gamzxx(i,j,k) - (                                   - gxx(i,j,k) * gxxz(i,j,k) )*HALF
+  fzz = fzz - Gamxzz*Lapx - Gamyzz*Lapy - Gamzzz*Lapz
-    Gamxyy(i,j,k) = Gamxyy(i,j,k) - (                                   - gyy(i,j,k) * gxxx(i,j,k) )*HALF
+  fxy = fxy - Gamxxy*Lapx - Gamyxy*Lapy - Gamzxy*Lapz
-    Gamyyy(i,j,k) = Gamyyy(i,j,k) - ( (chiy(i,j,k) + chiy(i,j,k))/chin_loc - gyy(i,j,k) * gxxy(i,j,k) )*HALF
+  fxz = fxz - Gamxxz*Lapx - Gamyxz*Lapy - Gamzxz*Lapz
-    Gamzyy(i,j,k) = Gamzyy(i,j,k) - (                                   - gyy(i,j,k) * gxxz(i,j,k) )*HALF
+  fyz = fyz - Gamxyz*Lapx - Gamyyz*Lapy - Gamzyz*Lapz
-    Gamxzz(i,j,k) = Gamxzz(i,j,k) - (                                   - gzz(i,j,k) * gxxx(i,j,k) )*HALF
-    Gamyzz(i,j,k) = Gamyzz(i,j,k) - (                                   - gzz(i,j,k) * gxxy(i,j,k) )*HALF
-    Gamzzz(i,j,k) = Gamzzz(i,j,k) - ( (chiz(i,j,k) + chiz(i,j,k))/chin_loc - gzz(i,j,k) * gxxz(i,j,k) )*HALF
-    Gamxxy(i,j,k) = Gamxxy(i,j,k) - ( chiy(i,j,k) /chin_loc - gxy(i,j,k) * gxxx(i,j,k) )*HALF
-    Gamyxy(i,j,k) = Gamyxy(i,j,k) - ( chix(i,j,k) /chin_loc - gxy(i,j,k) * gxxy(i,j,k) )*HALF
-    Gamzxy(i,j,k) = Gamzxy(i,j,k) - (                     - gxy(i,j,k) * gxxz(i,j,k) )*HALF
-    Gamxxz(i,j,k) = Gamxxz(i,j,k) - ( chiz(i,j,k) /chin_loc - gxz(i,j,k) * gxxx(i,j,k) )*HALF
-    Gamyxz(i,j,k) = Gamyxz(i,j,k) - (                     - gxz(i,j,k) * gxxy(i,j,k) )*HALF
-    Gamzxz(i,j,k) = Gamzxz(i,j,k) - ( chix(i,j,k) /chin_loc - gxz(i,j,k) * gxxz(i,j,k) )*HALF
-    Gamxyz(i,j,k) = Gamxyz(i,j,k) - (                     - gyz(i,j,k) * gxxx(i,j,k) )*HALF
-    Gamyyz(i,j,k) = Gamyyz(i,j,k) - ( chiz(i,j,k) /chin_loc - gyz(i,j,k) * gxxy(i,j,k) )*HALF
-    Gamzyz(i,j,k) = Gamzyz(i,j,k) - ( chiy(i,j,k) /chin_loc - gyz(i,j,k) * gxxz(i,j,k) )*HALF
 
-    fxx(i,j,k) = fxx(i,j,k) - Gamxxx(i,j,k)*Lapx(i,j,k) - Gamyxx(i,j,k)*Lapy(i,j,k) - Gamzxx(i,j,k)*Lapz(i,j,k)
+! store D^i D_i Lap in trK_rhs upto chi
-    fyy(i,j,k) = fyy(i,j,k) - Gamxyy(i,j,k)*Lapx(i,j,k) - Gamyyy(i,j,k)*Lapy(i,j,k) - Gamzyy(i,j,k)*Lapz(i,j,k)
+  trK_rhs =    gupxx * fxx + gupyy * fyy + gupzz * fzz + &
-    fzz(i,j,k) = fzz(i,j,k) - Gamxzz(i,j,k)*Lapx(i,j,k) - Gamyzz(i,j,k)*Lapy(i,j,k) - Gamzzz(i,j,k)*Lapz(i,j,k)
+        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz )
-    fxy(i,j,k) = fxy(i,j,k) - Gamxxy(i,j,k)*Lapx(i,j,k) - Gamyxy(i,j,k)*Lapy(i,j,k) - Gamzxy(i,j,k)*Lapz(i,j,k)
+#if 1        
-    fxz(i,j,k) = fxz(i,j,k) - Gamxxz(i,j,k)*Lapx(i,j,k) - Gamyxz(i,j,k)*Lapy(i,j,k) - Gamzxz(i,j,k)*Lapz(i,j,k)
+!! follow bam code
-    fyz(i,j,k) = fyz(i,j,k) - Gamxyz(i,j,k)*Lapx(i,j,k) - Gamyyz(i,j,k)*Lapy(i,j,k) - Gamzyz(i,j,k)*Lapz(i,j,k)
+  S =  chin1 * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + &
+     TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) )
+  f = F2o3 * trK * trK -(&
+       gupxx * ( &
+       gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + &
+       TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz) ) + &
+       gupyy * ( &
+       gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + &
+       TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz) ) + &
+       gupzz * ( &
+       gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + &
+       TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz) ) + &
+       TWO * ( &
+       gupxy * ( &
+       gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + &
+       gupxy * (Axx * Ayy + Axy * Axy) + &
+       gupxz * (Axx * Ayz + Axz * Axy) + &
+       gupyz * (Axy * Ayz + Axz * Ayy) ) + &
+       gupxz * ( &
+       gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + &
+       gupxy * (Axx * Ayz + Axy * Axz) + &
+       gupxz * (Axx * Azz + Axz * Axz) + &
+       gupyz * (Axy * Azz + Axz * Ayz) ) + &
+       gupyz * ( &
+       gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + &
+       gupxy * (Axy * Ayz + Ayy * Axz) + &
+       gupxz * (Axy * Azz + Ayz * Axz) + &
+       gupyz * (Ayy * Azz + Ayz * Ayz) ) )) -1.6d1*PI*rho + EIGHT * PI * S
+  f = - F1o3 *(  gupxx * fxx + gupyy * fyy + gupzz * fzz + &
+        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + alpn1/chin1*f)
   
-    trK_rhs(i,j,k) = gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + &
+  fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx
-                     TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k))
+  fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy
-  enddo
+  fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz
-  enddo
+  fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy
-  enddo
+  fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz
-  do k=1,ex(3)
+  fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz
-  do j=1,ex(2)
+#else        
-  do i=1,ex(1)
+! Add lapse and S_ij parts to Ricci tensor:
-    divb_loc = div_beta(i,j,k)
-    chin_loc = chin1(i,j,k)
 
-    S_loc = chin_loc * ( gupxx(i,j,k) * Sxx(i,j,k) + gupyy(i,j,k) * Syy(i,j,k) + gupzz(i,j,k) * Szz(i,j,k) + &
+  fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx
-           TWO * (gupxy(i,j,k) * Sxy(i,j,k) + gupxz(i,j,k) * Sxz(i,j,k) + gupyz(i,j,k) * Syz(i,j,k)) )
+  fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy
-    S(i,j,k) = S_loc
+  fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz
+  fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy
+  fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz
+  fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz
 
-    f_loc = F2o3 * trK(i,j,k) * trK(i,j,k) - ( &
+! Compute trace-free part (note: chi^-1 and chi cancel!):
-            gupxx(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axx(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + &
-                             gupzz(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + &
-                             TWO * (gupxy(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupxz(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + &
-                                    gupyz(i,j,k) * Axy(i,j,k) * Axz(i,j,k)) ) + &
-            gupyy(i,j,k) * ( gupxx(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayy(i,j,k) + &
-                             gupzz(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + &
-                             TWO * (gupxy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + gupxz(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + &
-                                    gupyz(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k)) ) + &
-            gupzz(i,j,k) * ( gupxx(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + &
-                             gupzz(i,j,k) * Azz(i,j,k) * Azz(i,j,k) + &
-                             TWO * (gupxy(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + gupxz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + &
-                                    gupyz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k)) ) + &
-            TWO * ( gupxy(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + &
-                                     gupzz(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + &
-                                     gupxy(i,j,k) * (Axx(i,j,k) * Ayy(i,j,k) + Axy(i,j,k) * Axy(i,j,k)) + &
-                                     gupxz(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Axy(i,j,k)) + &
-                                     gupyz(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Ayy(i,j,k)) ) + &
-                    gupxz(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + &
-                                     gupzz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + &
-                                     gupxy(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axy(i,j,k) * Axz(i,j,k)) + &
-                                     gupxz(i,j,k) * (Axx(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Axz(i,j,k)) + &
-                                     gupyz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Ayz(i,j,k)) ) + &
-                    gupyz(i,j,k) * ( gupxx(i,j,k) * Axy(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k) + &
-                                     gupzz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k) + &
-                                     gupxy(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Ayy(i,j,k) * Axz(i,j,k)) + &
-                                     gupxz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Axz(i,j,k)) + &
-                                     gupyz(i,j,k) * (Ayy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Ayz(i,j,k)) ) ) ) - &
-            F16 * PI * rho(i,j,k) + EIGHT * PI * S_loc
 
-    f_loc = -F1o3 * ( gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + &
+  f = F1o3 *(  gupxx * fxx + gupyy * fyy + gupzz * fzz + &
-            TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k)) + &
+        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) )
-            alpn1(i,j,k)/chin_loc * f_loc )
+#endif
-    f(i,j,k) = f_loc
 
-    l_fxx = alpn1(i,j,k) * (Rxx(i,j,k) - EIGHT * PI * Sxx(i,j,k)) - fxx(i,j,k)
+  Axx_rhs = fxx - gxx * f
-    l_fxy = alpn1(i,j,k) * (Rxy(i,j,k) - EIGHT * PI * Sxy(i,j,k)) - fxy(i,j,k)
+  Ayy_rhs = fyy - gyy * f
-    l_fxz = alpn1(i,j,k) * (Rxz(i,j,k) - EIGHT * PI * Sxz(i,j,k)) - fxz(i,j,k)
+  Azz_rhs = fzz - gzz * f
-    l_fyy = alpn1(i,j,k) * (Ryy(i,j,k) - EIGHT * PI * Syy(i,j,k)) - fyy(i,j,k)
+  Axy_rhs = fxy - gxy * f
-    l_fyz = alpn1(i,j,k) * (Ryz(i,j,k) - EIGHT * PI * Syz(i,j,k)) - fyz(i,j,k)
+  Axz_rhs = fxz - gxz * f
-    l_fzz = alpn1(i,j,k) * (Rzz(i,j,k) - EIGHT * PI * Szz(i,j,k)) - fzz(i,j,k)
+  Ayz_rhs = fyz - gyz * f
 
-    Axx_rhs(i,j,k) = l_fxx - gxx(i,j,k) * f_loc
+! Now: store A_il A^l_j into fij:
-    Ayy_rhs(i,j,k) = l_fyy - gyy(i,j,k) * f_loc
-    Azz_rhs(i,j,k) = l_fzz - gzz(i,j,k) * f_loc
-    Axy_rhs(i,j,k) = l_fxy - gxy(i,j,k) * f_loc
-    Axz_rhs(i,j,k) = l_fxz - gxz(i,j,k) * f_loc
-    Ayz_rhs(i,j,k) = l_fyz - gyz(i,j,k) * f_loc
 
-    fxx(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axx(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + &
+  fxx =       gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + &
-                 gupzz(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + TWO * (gupxy(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + &
+       TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz)
-                 gupxz(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyz(i,j,k) * Axy(i,j,k) * Axz(i,j,k))
+  fyy =       gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + &
-    fyy(i,j,k) = gupxx(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayy(i,j,k) + &
+       TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz)
-                 gupzz(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + TWO * (gupxy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + &
+  fzz =       gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + &
-                 gupxz(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + gupyz(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k))
+       TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz)
-    fzz(i,j,k) = gupxx(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + &
+  fxy =       gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + &
-                 gupzz(i,j,k) * Azz(i,j,k) * Azz(i,j,k) + TWO * (gupxy(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + &
+              gupxy *(Axx * Ayy + Axy * Axy)                            + &
-                 gupxz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + gupyz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k))
+              gupxz *(Axx * Ayz + Axz * Axy)                            + &
-    fxy(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + &
+              gupyz *(Axy * Ayz + Axz * Ayy)
-                 gupzz(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + gupxy(i,j,k) * (Axx(i,j,k) * Ayy(i,j,k) + Axy(i,j,k) * Axy(i,j,k)) + &
+  fxz =       gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + &
-                 gupxz(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Axy(i,j,k)) + &
+              gupxy *(Axx * Ayz + Axy * Axz)                            + &
-                 gupyz(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Ayy(i,j,k))
+              gupxz *(Axx * Azz + Axz * Axz)                            + &
-    fxz(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + &
+              gupyz *(Axy * Azz + Axz * Ayz)
-                 gupzz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + gupxy(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axy(i,j,k) * Axz(i,j,k)) + &
+  fyz =       gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + &
-                 gupxz(i,j,k) * (Axx(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Axz(i,j,k)) + &
+              gupxy *(Axy * Ayz + Ayy * Axz)                            + &
-                 gupyz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Ayz(i,j,k))
+              gupxz *(Axy * Azz + Ayz * Axz)                            + &
-    fyz(i,j,k) = gupxx(i,j,k) * Axy(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k) + &
+              gupyz *(Ayy * Azz + Ayz * Ayz)
-                 gupzz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k) + gupxy(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Ayy(i,j,k) * Axz(i,j,k)) + &
-                 gupxz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Axz(i,j,k)) + &
-                 gupyz(i,j,k) * (Ayy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Ayz(i,j,k))
 
-    trK_rhs(i,j,k) = chin_loc * trK_rhs(i,j,k)
+  f = chin1
+! store D^i D_i Lap in trK_rhs
+  trK_rhs = f*trK_rhs
           
-    Axx_rhs(i,j,k) = chin_loc * Axx_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axx(i,j,k) - TWO * fxx(i,j,k)) + &
+  Axx_rhs =           f * Axx_rhs+ alpn1 * (trK * Axx - TWO * fxx)  + &
-                     TWO * (Axx(i,j,k) * betaxx(i,j,k) + Axy(i,j,k) * betayx(i,j,k) + Axz(i,j,k) * betazx(i,j,k)) - &
+           TWO * (  Axx * betaxx +   Axy * betayx +   Axz * betazx )- &
-                     F2o3 * Axx(i,j,k) * divb_loc
+             F2o3 * Axx * div_beta
-    Ayy_rhs(i,j,k) = chin_loc * Ayy_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Ayy(i,j,k) - TWO * fyy(i,j,k)) + &
-                     TWO * (Axy(i,j,k) * betaxy(i,j,k) + Ayy(i,j,k) * betayy(i,j,k) + Ayz(i,j,k) * betazy(i,j,k)) - &
-                     F2o3 * Ayy(i,j,k) * divb_loc
-    Azz_rhs(i,j,k) = chin_loc * Azz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Azz(i,j,k) - TWO * fzz(i,j,k)) + &
-                     TWO * (Axz(i,j,k) * betaxz(i,j,k) + Ayz(i,j,k) * betayz(i,j,k) + Azz(i,j,k) * betazz(i,j,k)) - &
-                     F2o3 * Azz(i,j,k) * divb_loc
-    Axy_rhs(i,j,k) = chin_loc * Axy_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axy(i,j,k) - TWO * fxy(i,j,k)) + &
-                     Axx(i,j,k) * betaxy(i,j,k) + Axz(i,j,k) * betazy(i,j,k) + Ayy(i,j,k) * betayx(i,j,k) + &
-                     Ayz(i,j,k) * betazx(i,j,k) + F1o3 * Axy(i,j,k) * divb_loc - Axy(i,j,k) * betazz(i,j,k)
-    Ayz_rhs(i,j,k) = chin_loc * Ayz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Ayz(i,j,k) - TWO * fyz(i,j,k)) + &
-                     Axy(i,j,k) * betaxz(i,j,k) + Ayy(i,j,k) * betayz(i,j,k) + Axz(i,j,k) * betaxy(i,j,k) + &
-                     Azz(i,j,k) * betazy(i,j,k) + F1o3 * Ayz(i,j,k) * divb_loc - Ayz(i,j,k) * betaxx(i,j,k)
-    Axz_rhs(i,j,k) = chin_loc * Axz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axz(i,j,k) - TWO * fxz(i,j,k)) + &
-                     Axx(i,j,k) * betaxz(i,j,k) + Axy(i,j,k) * betayz(i,j,k) + Ayz(i,j,k) * betayx(i,j,k) + &
-                     Azz(i,j,k) * betazx(i,j,k) + F1o3 * Axz(i,j,k) * divb_loc - Axz(i,j,k) * betayy(i,j,k)
 
-    trK_rhs(i,j,k) = - trK_rhs(i,j,k) + alpn1(i,j,k) * ( F1o3 * trK(i,j,k) * trK(i,j,k) + &
+  Ayy_rhs =           f * Ayy_rhs+ alpn1 * (trK * Ayy - TWO * fyy)  + &
-                    gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + &
+           TWO * (  Axy * betaxy +   Ayy * betayy +   Ayz * betazy )- &
-                    TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k)) + &
+             F2o3 * Ayy * div_beta
-                    FOUR * PI * (rho(i,j,k) + S_loc) )
+
-  enddo
+  Azz_rhs =           f * Azz_rhs+ alpn1 * (trK * Azz - TWO * fzz)  + &
-  enddo
+           TWO * (  Axz * betaxz +   Ayz * betayz +   Azz * betazz )- &
-  enddo
+             F2o3 * Azz * div_beta
+
+  Axy_rhs =           f * Axy_rhs+ alpn1 *( trK * Axy  - TWO * fxy )+ &
+                    Axx * betaxy                  +   Axz * betazy  + &
+                                     Ayy * betayx +   Ayz * betazx  + &
+             F1o3 * Axy * div_beta                -   Axy * betazz
+
+  Ayz_rhs =           f * Ayz_rhs+ alpn1 *( trK * Ayz  - TWO * fyz )+ &
+                    Axy * betaxz +   Ayy * betayz                   + &
+                    Axz * betaxy                  +   Azz * betazy  + &
+             F1o3 * Ayz * div_beta                -   Ayz * betaxx
+ 
+  Axz_rhs =           f * Axz_rhs+ alpn1 *( trK * Axz  - TWO * fxz )+ &
+                    Axx * betaxz +   Axy * betayz                   + &
+                                     Ayz * betayx +   Azz * betazx  + &
+             F1o3 * Axz * div_beta                -   Axz * betayy      !rhs for Aij
+
+! Compute trace of S_ij
+
+  S =  f * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + &
+     TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) )
+
+  trK_rhs = - trK_rhs + alpn1 *( F1o3 * trK * trK         + &
+                gupxx * fxx + gupyy * fyy + gupzz * fzz   + &
+        TWO * ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + &
+       FOUR * PI * ( rho + S ))                                !rhs for trK
   
 !!!! gauge variable part
 
@@ -1000,15 +948,15 @@
 !!!!!!!!!advection term + Kreiss-Oliger dissipation (merged for cache efficiency)
 ! lopsided_kodis shares the symmetry_bd buffer between advection and
 ! dissipation, eliminating redundant full-grid copies. For metric variables
-! gxx/gyy/gzz (=dxx/dyy/dzz+1): stencil coefficients sum to zero,
+! gxx/gyy/gzz (=dxx/dyy/dzz+1): kodis stencil coefficients sum to zero,
 ! so the constant offset has no effect on dissipation.
 
-  call lopsided_kodis(ex,X,Y,Z,dxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
   call lopsided_kodis(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,eps)
   call lopsided_kodis(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,eps)
-  call lopsided_kodis(ex,X,Y,Z,dyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
   call lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps)
-  call lopsided_kodis(ex,X,Y,Z,dzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
 
   call lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
   call lopsided_kodis(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,eps)

AMSS_NCKU_source/bssn_rhs_c.C

@@ -716,7 +716,6 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
         // 24ms //
         fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);
-        fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);
 
         // 6ms //
         for (int i=0;i<all;i+=1) {
@@ -1014,24 +1013,17 @@ int f_compute_rhs_bssn(int *ex, double &T,
             betaz_rhs[i] = FF * dtSfz[i];
 
             reta[i] =
-                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
+                gupxx[i] * chix[i] * chix[i]
-                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
+                + gupyy[i] * chiy[i] * chiy[i]
-                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
+                + gupzz[i] * chiz[i] * chiz[i]
-                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
+                + TWO * ( gupxy[i] * chix[i] * chiy[i]
-                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
+                        + gupxz[i] * chix[i] * chiz[i]
-                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
+                        + gupyz[i] * chiy[i] * chiz[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];
@@ -1039,24 +1031,17 @@ int f_compute_rhs_bssn(int *ex, double &T,
             dtSfz_rhs[i] = Gamz_rhs[i] - reta[i] * dtSfz[i];
             #elif (GAUGE == 4 || GAUGE == 5)
             reta[i] =
-                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
+                gupxx[i] * chix[i] * chix[i]
-                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
+                + gupyy[i] * chiy[i] * chiy[i]
-                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
+                + gupzz[i] * chiz[i] * chiz[i]
-                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
+                + TWO * ( gupxy[i] * chix[i] * chiy[i]
-                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
+                        + gupxz[i] * chix[i] * chiz[i]
-                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
+                        + gupyz[i] * chiy[i] * chiz[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

@@ -27,10 +27,6 @@ using namespace std;
 #include "cgh.h"
 #include "Parallel.h"
 #include "parameters.h"
-#ifdef USE_GPU
-#include "bssn_gpu.h"
-#include "bssn_cuda_ops.h"
-#endif
 
 //================================================================================================
 
@@ -889,13 +885,6 @@ void cgh::recompose_cgh(int nprocs, bool *lev_flag,
       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();
-      bssn_cuda_release_rk4_caches();
-      bssn_cuda_release_interp_caches();
-      patch_release_interp_plan_cache();
-#endif
       Parallel::KillBlocks(PatL[lev]);
       PatL[lev]->destroyList();
       PatL[lev] = tmPat;
@@ -925,13 +914,6 @@ void cgh::recompose_cgh(int nprocs, bool *lev_flag,
       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();
-      bssn_cuda_release_rk4_caches();
-      bssn_cuda_release_interp_caches();
-      patch_release_interp_plan_cache();
-#endif
       Parallel::KillBlocks(PatL[lev]);
       PatL[lev]->destroyList();
       PatL[lev] = tmPat;
@@ -1540,13 +1522,6 @@ void cgh::recompose_cgh_Onelevel(int nprocs, int lev,
   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();
-  bssn_cuda_release_rk4_caches();
-  bssn_cuda_release_interp_caches();
-  patch_release_interp_plan_cache();
-#endif
   Parallel::KillBlocks(PatL[lev]);
   PatL[lev]->destroyList();
   PatL[lev] = tmPat;
@@ -1570,13 +1545,6 @@ void cgh::recompose_cgh_Onelevel(int nprocs, int lev,
   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();
-  bssn_cuda_release_rk4_caches();
-  bssn_cuda_release_interp_caches();
-  patch_release_interp_plan_cache();
-#endif
   Parallel::KillBlocks(PatL[lev]);
   PatL[lev]->destroyList();
   PatL[lev] = tmPat;

AMSS_NCKU_source/makefile

@@ -64,8 +64,8 @@ lopsided_c.o: lopsided_c.C
 lopsided_kodis_c.o: lopsided_kodis_c.C
 	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
 
-#interp_lb_profile.o: interp_lb_profile.C interp_lb_profile.h
+interp_lb_profile.o: interp_lb_profile.C interp_lb_profile.h
-#	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
 
 ## TwoPunctureABE uses fixed optimal flags with its own PGO profile, independent of CXXAPPFLAGS
 TP_PROFDATA = /home/$(shell whoami)/AMSS-NCKU/pgo_profile/TwoPunctureABE.profdata
@@ -106,11 +106,12 @@ C++FILES = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.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 bssn_class.o surface_integral.o ShellPatch.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_cuda_step.o writefile_f.o
+	   NullShellPatch2_Evo.o \
+	   bssn_gpu_class.o bssn_step_gpu.o bssn_macro.o writefile_f.o
 
 F90FILES_BASE = enforce_algebra.o fmisc.o initial_puncture.o prolongrestrict.o\
 	   prolongrestrict_cell.o prolongrestrict_vertex.o\
@@ -142,7 +143,7 @@ initial_guess.o Newton.o Jacobian.o ilucg.o IntPnts0.o IntPnts.o
 
 TwoPunctureFILES = TwoPunctureABE.o TwoPunctures.o 
 
-CUDAFILES = bssn_gpu.o bssn_cuda_ops.o
+CUDAFILES = bssn_gpu.o bssn_gpu_rhs_ss.o
 
 # file dependences
 $(C++FILES) $(C++FILES_GPU) $(F90FILES) $(CFILES) $(AHFDOBJS) $(CUDAFILES): macrodef.fh

AMSS_NCKU_source/makefile.inc

@@ -9,7 +9,6 @@ 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)
@@ -25,8 +24,6 @@ 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

@@ -182,7 +182,6 @@ surface_integral::surface_integral(int iSymmetry) : Symmetry(iSymmetry)
 //|============================================================================
 surface_integral::~surface_integral()
 {
-  release_cached_buffers();
   delete[] nx_g;
   delete[] ny_g;
   delete[] nz_g;
@@ -191,50 +190,6 @@ surface_integral::~surface_integral()
   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
@@ -256,7 +211,14 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
 
   int n;
   double *pox[3];
-  get_surface_points(rex, pox);
+  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 mp, Lp, Nmin, Nmax;
   mp = n_tot / cpusize;
@@ -272,7 +234,8 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
     Nmax = Nmin + mp - 1;
   }
 
-  double *shellf = get_shellf_buffer(InList);
+  double *shellf;
+  shellf = new double[n_tot * InList];
 
   GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax);
 
@@ -412,6 +375,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();
@@ -437,9 +404,17 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
 
   int n;
   double *pox[3];
-  get_surface_points(rex, pox);
+  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 = get_shellf_buffer(InList);
+  double *shellf;
+  shellf = new double[n_tot * InList];
 
   //    misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Interp_Points");
 
@@ -602,6 +577,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();
@@ -620,9 +599,17 @@ void surface_integral::surf_Wave(double rex, int lev, ShellPatch *GH, var *Rpsi4
 
   int n;
   double *pox[3];
-  get_surface_points(rex, pox);
+  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 = get_shellf_buffer(InList);
+  double *shellf;
+  shellf = new double[n_tot * InList];
 
   GH->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry);
 
@@ -2583,6 +2570,10 @@ 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();
 }
 void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var *trK,
@@ -2648,9 +2639,17 @@ void surface_integral::surf_MassPAng(double rex, int lev, cgh *GH, var *chi, var
 
   int n;
   double *pox[3];
-  get_surface_points(rex, pox);
+  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 = get_shellf_buffer(InList);
+  double *shellf;
+  shellf = new double[n_tot * InList];
 
   // we have assumed there is only one box on this level,
   // so we do not need loop boxes
@@ -2840,6 +2839,10 @@ 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();
 }
 //|----------------------------------------------------------------

AMSS_NCKU_source/surface_integral.h

@@ -23,21 +23,11 @@ using namespace std;
 #include "NullShellPatch2.h"
 #include "var.h"
 #include "monitor.h"
-#include <map>
 
 class surface_integral
 {
 
 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;
@@ -46,12 +36,6 @@ private:
 
 	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);

makefile_and_run.py

@@ -9,7 +9,6 @@
 
 
 import AMSS_NCKU_Input as input_data
-import os
 import subprocess
 import time
 
@@ -58,48 +57,6 @@ 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
-
-##################################################################
-
-
 
 ##################################################################
 
@@ -189,29 +146,16 @@ 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"):
-        run_env = prepare_gpu_runtime_env()
+        mpi_command         = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + " ./ABEGPU"
-        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_process = subprocess.Popen(mpi_command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)
-        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: