Switch to NVIDIA HPC Toolchain

Rewritten done by Codex.
This still has errors, do not pick this one now.

.gitignore

@@ -1,6 +1,6 @@
 __pycache__
 GW150914
-GW150914-origin
+GW150914*
 docs
 *.tmp
-
+.codex

.idea/vcs.xml

@@ -1,6 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<project version="4">
-  <component name="VcsDirectoryMappings">
-    <mapping directory="" vcs="Git" />
-  </component>
-</project>

AMSS_NCKU_Input.py

@@ -16,9 +16,9 @@ import numpy
 File_directory   = "GW150914"                    ## output file directory
 Output_directory = "binary_output"               ## binary data file directory
                                                  ## The file directory name should not be too long
-MPI_processes    = 64                             ## number of mpi processes used in the simulation
+MPI_processes    = 8                             ## number of mpi processes used in the simulation
 
-GPU_Calculation  = "no"                          ## Use GPU or not 
+GPU_Calculation  = "yes"                          ## Use GPU or not 
                                                  ## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)
 CPU_Part         = 1.0
 GPU_Part         = 0.0

AMSS_NCKU_Program.py

@@ -126,6 +126,11 @@ setup.generate_AMSSNCKU_input()
 #inputvalue = input()  ## Wait for user input (press Enter) to proceed
 #print()
 
+setup.print_puncture_information()
+
+
+##################################################################
+
 ## Generate AMSS-NCKU program input files based on the configured parameters
 
 print(                                                                                   )
@@ -253,7 +258,7 @@ print()
 if (input_data.GPU_Calculation == "no"):
     ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE")
 elif (input_data.GPU_Calculation == "yes"):
-    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABEGPU")
+    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE_CUDA")
 
 if not os.path.exists( ABE_file ):
     print(                                                                                                  )
@@ -307,7 +312,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
     
     import generate_TwoPuncture_input
     
-    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input(numerical_grid.puncture_data)
+    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
     
     print(                                                                                              )
     print( " The input parfile for the TwoPunctureABE executable has been generated. " )
@@ -349,7 +354,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
 
 import renew_puncture_parameter
     
-renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory, numerical_grid.puncture_data)
+renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory)
 
 
 ## Generated AMSS-NCKU input filename

AMSS_NCKU_source/BSSN_GPU/bssn_gpu.h

@@ -1,73 +0,0 @@
-
-#ifndef BSSN_GPU_H_
-#define BSSN_GPU_H_
-#include "bssn_macro.h"
-#include "macrodef.fh"
-
-#define DEVICE_ID 0
-// #define DEVICE_ID_BY_MPI_RANK
-#define GRID_DIM 256
-#define BLOCK_DIM 128
-
-#define _FH2_(i, j, k) fh[(i) + (j) * _1D_SIZE[2] + (k) * _2D_SIZE[2]]
-#define _FH3_(i, j, k) fh[(i) + (j) * _1D_SIZE[3] + (k) * _2D_SIZE[3]]
-#define pow2(x) ((x) * (x))
-#define TimeBetween(a, b) ((b.tv_sec - a.tv_sec) + (b.tv_usec - a.tv_usec) / 1000000.0f)
-#define M_ metac.
-#define Mh_ meta->
-#define Ms_ metassc.
-#define Msh_ metass->
-
-// #define TIMING
-
-#define RHS_SS_PARA int calledby, int mpi_rank, int *ex, double &T, double *crho, double *sigma, double *R, double *X, double *Y, double *Z, double *drhodx, double *drhody, double *drhodz, double *dsigmadx, double *dsigmady, double *dsigmadz, double *dRdx, double *dRdy, double *dRdz, double *drhodxx, double *drhodxy, double *drhodxz, double *drhodyy, double *drhodyz, double *drhodzz, double *dsigmadxx, double *dsigmadxy, double *dsigmadxz, double *dsigmadyy, double *dsigmadyz, double *dsigmadzz, double *dRdxx, double *dRdxy, double *dRdxz, double *dRdyy, double *dRdyz, double *dRdzz, double *chi, double *trK, double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz, double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz, double *Gamx, double *Gamy, double *Gamz, double *Lap, double *betax, double *betay, double *betaz, double *dtSfx, double *dtSfy, double *dtSfz, double *chi_rhs, double *trK_rhs, double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs, double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs, double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs, double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs, double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs, double *rho, double *Sx, double *Sy, double *Sz, double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz, double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz, double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz, double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz, double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz, double *ham_Res, double *movx_Res, double *movy_Res, double *movz_Res, double *Gmx_Res, double *Gmy_Res, double *Gmz_Res, int &Symmetry, int &Lev, double &eps, int &sst, int &co
-
-/**  main function */
-int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,
-            double *X, double *Y, double *Z,
-
-            double *chi, double *trK,
-
-            double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
-
-            double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
-
-            double *Gamx, double *Gamy, double *Gamz,
-
-            double *Lap, double *betax, double *betay, double *betaz,
-
-            double *dtSfx, double *dtSfy, double *dtSfz,
-
-            double *chi_rhs, double *trK_rhs,
-
-            double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
-
-            double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
-
-            double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
-
-            double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
-
-            double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
-
-            double *rho, double *Sx, double *Sy, double *Sz, double *Sxx,
-            double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
-
-            double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
-
-            double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
-
-            double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
-
-            double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
-
-            double *ham_Res, double *movx_Res, double *movy_Res, double *movz_Res,
-            double *Gmx_Res, double *Gmy_Res, double *Gmz_Res,
-            int &Symmetry, int &Lev, double &eps, int &co);
-
-int gpu_rhs_ss(RHS_SS_PARA);
-
-/** Init GPU side data in GPUMeta. */
-// void init_fluid_meta_gpu(GPUMeta *gpu_meta);
-
-#endif

AMSS_NCKU_source/BSSN_GPU/bssn_gpu_class.h

@@ -1,210 +0,0 @@
-
-#ifndef BSSN_GPU_CLASS_H
-#define BSSN_GPU_CLASS_H
-
-#ifdef newc
-#include <iostream>
-#include <iomanip>
-#include <fstream>
-#include <cstdlib>
-#include <string>
-#include <cmath>
-using namespace std;
-#else
-#include <iostream.h>
-#include <iomanip.h>
-#include <fstream.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#endif
-
-#include <mpi.h>
-
-#include "macrodef.h"
-#include "cgh.h"
-#include "ShellPatch.h"
-#include "misc.h"
-#include "var.h"
-#include "MyList.h"
-#include "monitor.h"
-#include "surface_integral.h"
-#include "checkpoint.h"
-
-// added by yangquan
-#include "bssn_macro.h"
-
-extern void setpbh(int iBHN, double **iPBH, double *iMass, int rBHN);
-
-class bssn_class
-{
-public:
-       // added by yangquan
-       //----------------------
-       int gpu_num_mynode;
-       int cpu_core_num_mynode;
-       int mpi_process_num_mynode;
-       int my_sequence_mynode;
-       int mynode_id;
-       int use_gpu;
-
-       virtual void Step_GPU(int lev, int YN);
-       virtual void Get_runtime_envirment();
-       // virtual void Step_OPENMP(int lev,int YN);
-       //----------------------
-
-       int ngfs;
-       int nprocs, myrank;
-       cgh *GH;
-       ShellPatch *SH;
-       double PhysTime;
-
-       int checkrun;
-       char checkfilename[50];
-       int Steps;
-       double StartTime, TotalTime;
-       double AnasTime, DumpTime, d2DumpTime, CheckTime;
-       double LastAnas, LastConsOut;
-       double Courant;
-       double numepss, numepsb, numepsh;
-       int Symmetry;
-       int maxl, decn;
-       double maxrex, drex;
-       int trfls, a_lev;
-
-       double dT;
-       double chitiny;
-
-       double **Porg0, **Porgbr, **Porg, **Porg1, **Porg_rhs;
-       int BH_num, BH_num_input;
-       double *Mass, *Pmom, *Spin;
-       double ADMMass;
-
-       var *phio, *trKo;
-       var *gxxo, *gxyo, *gxzo, *gyyo, *gyzo, *gzzo;
-       var *Axxo, *Axyo, *Axzo, *Ayyo, *Ayzo, *Azzo;
-       var *Gmxo, *Gmyo, *Gmzo;
-       var *Lapo, *Sfxo, *Sfyo, *Sfzo;
-       var *dtSfxo, *dtSfyo, *dtSfzo;
-
-       var *phi0, *trK0;
-       var *gxx0, *gxy0, *gxz0, *gyy0, *gyz0, *gzz0;
-       var *Axx0, *Axy0, *Axz0, *Ayy0, *Ayz0, *Azz0;
-       var *Gmx0, *Gmy0, *Gmz0;
-       var *Lap0, *Sfx0, *Sfy0, *Sfz0;
-       var *dtSfx0, *dtSfy0, *dtSfz0;
-
-       var *phi, *trK;
-       var *gxx, *gxy, *gxz, *gyy, *gyz, *gzz;
-       var *Axx, *Axy, *Axz, *Ayy, *Ayz, *Azz;
-       var *Gmx, *Gmy, *Gmz;
-       var *Lap, *Sfx, *Sfy, *Sfz;
-       var *dtSfx, *dtSfy, *dtSfz;
-
-       var *phi1, *trK1;
-       var *gxx1, *gxy1, *gxz1, *gyy1, *gyz1, *gzz1;
-       var *Axx1, *Axy1, *Axz1, *Ayy1, *Ayz1, *Azz1;
-       var *Gmx1, *Gmy1, *Gmz1;
-       var *Lap1, *Sfx1, *Sfy1, *Sfz1;
-       var *dtSfx1, *dtSfy1, *dtSfz1;
-
-       var *phi_rhs, *trK_rhs;
-       var *gxx_rhs, *gxy_rhs, *gxz_rhs, *gyy_rhs, *gyz_rhs, *gzz_rhs;
-       var *Axx_rhs, *Axy_rhs, *Axz_rhs, *Ayy_rhs, *Ayz_rhs, *Azz_rhs;
-       var *Gmx_rhs, *Gmy_rhs, *Gmz_rhs;
-       var *Lap_rhs, *Sfx_rhs, *Sfy_rhs, *Sfz_rhs;
-       var *dtSfx_rhs, *dtSfy_rhs, *dtSfz_rhs;
-
-       var *rho, *Sx, *Sy, *Sz, *Sxx, *Sxy, *Sxz, *Syy, *Syz, *Szz;
-
-       var *Gamxxx, *Gamxxy, *Gamxxz, *Gamxyy, *Gamxyz, *Gamxzz;
-       var *Gamyxx, *Gamyxy, *Gamyxz, *Gamyyy, *Gamyyz, *Gamyzz;
-       var *Gamzxx, *Gamzxy, *Gamzxz, *Gamzyy, *Gamzyz, *Gamzzz;
-
-       var *Rxx, *Rxy, *Rxz, *Ryy, *Ryz, *Rzz;
-
-       var *Rpsi4, *Ipsi4;
-       var *t1Rpsi4, *t1Ipsi4, *t2Rpsi4, *t2Ipsi4;
-
-       var *Cons_Ham, *Cons_Px, *Cons_Py, *Cons_Pz, *Cons_Gx, *Cons_Gy, *Cons_Gz;
-
-#ifdef Point_Psi4
-       var *phix, *phiy, *phiz;
-       var *trKx, *trKy, *trKz;
-       var *Axxx, *Axxy, *Axxz;
-       var *Axyx, *Axyy, *Axyz;
-       var *Axzx, *Axzy, *Axzz;
-       var *Ayyx, *Ayyy, *Ayyz;
-       var *Ayzx, *Ayzy, *Ayzz;
-       var *Azzx, *Azzy, *Azzz;
-#endif
-       // FIXME: uc = StateList, up = OldStateList, upp = SynchList_cor; so never touch these three data
-       MyList<var> *StateList, *SynchList_pre, *SynchList_cor, *RHSList;
-       MyList<var> *OldStateList, *DumpList;
-       MyList<var> *ConstraintList;
-
-       monitor *ErrorMonitor, *Psi4Monitor, *BHMonitor, *MAPMonitor;
-       monitor *ConVMonitor;
-       surface_integral *Waveshell;
-       checkpoint *CheckPoint;
-
-public:
-       bssn_class(double Couranti, double StartTimei, double TotalTimei, double DumpTimei, double d2DumpTimei, double CheckTimei, double AnasTimei,
-                  int Symmetryi, int checkruni, char *checkfilenamei, double numepssi, double numepsbi, double numepshi,
-                  int a_levi, int maxli, int decni, double maxrexi, double drexi);
-       ~bssn_class();
-
-       void Evolve(int Steps);
-       void RecursiveStep(int lev);
-#if (PSTR == 1)
-       void ParallelStep();
-       void SHStep();
-#endif
-       void RestrictProlong(int lev, int YN, bool BB, MyList<var> *SL, MyList<var> *OL, MyList<var> *corL);
-       void RestrictProlong_aux(int lev, int YN, bool BB, MyList<var> *SL, MyList<var> *OL, MyList<var> *corL);
-       void RestrictProlong(int lev, int YN, bool BB);
-       void ProlongRestrict(int lev, int YN, bool BB);
-       void Setup_Black_Hole_position();
-       void compute_Porg_rhs(double **BH_PS, double **BH_RHS, var *forx, var *fory, var *forz, int lev);
-       bool read_Pablo_file(int *ext, double *datain, char *filename);
-       void write_Pablo_file(int *ext, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax,
-                             char *filename);
-       void AnalysisStuff(int lev, double dT_lev);
-       void Setup_KerrSchild();
-       void Enforce_algcon(int lev, int fg);
-
-       void testRestrict();
-       void testOutBd();
-
-       virtual void Setup_Initial_Data_Lousto();
-       virtual void Setup_Initial_Data_Cao();
-       virtual void Initialize();
-       virtual void Read_Ansorg();
-       virtual void Read_Pablo() {};
-       virtual void Compute_Psi4(int lev);
-       virtual void Step(int lev, int YN);
-       virtual void Interp_Constraint(bool infg);
-       virtual void Constraint_Out();
-       virtual void Compute_Constraint();
-
-#ifdef With_AHF
-protected:
-       MyList<var> *AHList, *AHDList, *GaugeList;
-       int AHfindevery;
-       double AHdumptime;
-       int *lastahdumpid, HN_num; // number of possible horizons
-       int *findeveryl;
-       double *xc, *yc, *zc, *xr, *yr, *zr;
-       bool *trigger;
-       double *dTT;
-       int *dumpid;
-
-public:
-       void AH_Prepare_derivatives();
-       bool AH_Interp_Points(MyList<var> *VarList,
-                             int NN, double **XX,
-                             double *Shellf, int Symmetryi);
-       void AH_Step_Find(int lev, double dT_lev);
-#endif
-};
-#endif /* BSSN_GPU_CLASS_H */

AMSS_NCKU_source/Parallel.C

@@ -4,6 +4,273 @@
 #include "prolongrestrict.h"
 #include "misc.h"
 #include "parameters.h"
+#include <cstdlib>
+#include <cstdio>
+
+#ifndef USE_CUDA_Z4C
+#define USE_CUDA_Z4C 0
+#endif
+
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+#include <cuda_runtime_api.h>
+#endif
+#if USE_CUDA_BSSN
+#include "bssn_rhs_cuda.h"
+#endif
+#if USE_CUDA_Z4C
+#include "z4c_rhs_cuda.h"
+#endif
+
+namespace {
+
+struct SyncProfileStats
+{
+  long long start_calls;
+  long long finish_calls;
+  double start_sec;
+  double finish_sec;
+  double direct_pack_sec;
+  double direct_unpack_sec;
+  double wait_sec;
+};
+
+SyncProfileStats &sync_profile_stats()
+{
+  static SyncProfileStats stats = {0, 0, 0.0, 0.0, 0.0, 0.0, 0.0};
+  return stats;
+}
+
+bool sync_profile_enabled()
+{
+  static int enabled = -1;
+  if (enabled < 0)
+  {
+    const char *env = getenv("AMSS_PROFILE_SYNC");
+    enabled = (env && atoi(env) != 0) ? 1 : 0;
+  }
+  return enabled != 0;
+}
+
+int sync_profile_every()
+{
+  static int every = -1;
+  if (every < 0)
+  {
+    const char *env = getenv("AMSS_PROFILE_SYNC_EVERY");
+    every = (env && atoi(env) > 0) ? atoi(env) : 100;
+  }
+  return every;
+}
+
+void sync_profile_maybe_log()
+{
+  if (!sync_profile_enabled())
+    return;
+  SyncProfileStats &stats = sync_profile_stats();
+  if (stats.finish_calls <= 0 || stats.finish_calls % sync_profile_every() != 0)
+    return;
+  int rank = 0;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  fprintf(stderr,
+          "[AMSS-SYNC][rank %d] start=%lld finish=%lld avg_start=%.6f s avg_finish=%.6f s avg_wait=%.6f s avg_cuda_pack=%.6f s avg_cuda_unpack=%.6f s\n",
+          rank,
+          stats.start_calls,
+          stats.finish_calls,
+          stats.start_calls ? stats.start_sec / (double)stats.start_calls : 0.0,
+          stats.finish_calls ? stats.finish_sec / (double)stats.finish_calls : 0.0,
+          stats.finish_calls ? stats.wait_sec / (double)stats.finish_calls : 0.0,
+          stats.finish_calls ? stats.direct_pack_sec / (double)stats.finish_calls : 0.0,
+          stats.finish_calls ? stats.direct_unpack_sec / (double)stats.finish_calls : 0.0);
+  fflush(stderr);
+}
+
+bool cuda_sync_pinned_enabled()
+{
+  static int enabled = -1;
+  if (enabled < 0)
+  {
+    const char *env = getenv("AMSS_CUDA_PINNED_SYNC");
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+    enabled = (!env || atoi(env) != 0) ? 1 : 0;
+#else
+    enabled = 0;
+#endif
+  }
+  return enabled != 0;
+}
+
+void free_comm_buffer(double *&ptr, unsigned char &is_pinned)
+{
+  if (!ptr)
+    return;
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+  if (is_pinned)
+    cudaFreeHost(ptr);
+  else
+    delete[] ptr;
+#else
+  delete[] ptr;
+#endif
+  ptr = 0;
+  is_pinned = 0;
+}
+
+double *alloc_comm_buffer(int length, unsigned char &is_pinned)
+{
+  is_pinned = 0;
+  if (length <= 0)
+    return 0;
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+  if (cuda_sync_pinned_enabled())
+  {
+    double *ptr = 0;
+    cudaError_t err = cudaMallocHost((void **)&ptr, (size_t)length * sizeof(double));
+    if (err == cudaSuccess)
+    {
+      is_pinned = 1;
+      return ptr;
+    }
+  }
+#endif
+  return new double[length];
+}
+
+void ensure_comm_buffer(double **buffers, unsigned char *pinned_flags, int *caps, int idx, int length)
+{
+  if (length <= caps[idx])
+    return;
+  free_comm_buffer(buffers[idx], pinned_flags[idx]);
+  buffers[idx] = alloc_comm_buffer(length, pinned_flags[idx]);
+  if (!buffers[idx])
+  {
+    fprintf(stderr, "Parallel: failed to allocate communication buffer (%d doubles)\n", length);
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+  caps[idx] = length;
+}
+
+int cuda_seg_begin(const Parallel::gridseg *seg, Block *bg, int dir)
+{
+  const double dx = bg->getdX(dir);
+  return (int)floor((seg->llb[dir] - bg->bbox[dir]) / dx + 0.5);
+}
+
+int cuda_state_var_count(MyList<var> *src_vars, MyList<var> *dst_vars)
+{
+  int count = 0;
+  while (src_vars && dst_vars)
+  {
+    ++count;
+    src_vars = src_vars->next;
+    dst_vars = dst_vars->next;
+  }
+  return (src_vars || dst_vars) ? -1 : count;
+}
+
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+bool cuda_state_count_direct_supported(int state_count)
+{
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  return state_count == Z4C_CUDA_STATE_COUNT;
+#elif USE_CUDA_BSSN
+  return state_count > 0 && state_count <= BSSN_CUDA_STATE_COUNT;
+#else
+  (void)state_count;
+  return false;
+#endif
+}
+
+bool cuda_can_direct_pack(const Parallel::gridseg *src, const Parallel::gridseg *dst, int type)
+{
+  if (type != 1 || !src || !dst || !src->Bg)
+    return false;
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  return z4c_cuda_has_resident_state(src->Bg) != 0;
+#elif USE_CUDA_BSSN
+  return bssn_cuda_has_resident_state(src->Bg) != 0;
+#else
+  return false;
+#endif
+}
+
+bool cuda_can_direct_unpack(const Parallel::gridseg *dst, int type)
+{
+  if (type != 1 || !dst || !dst->Bg)
+    return false;
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  return z4c_cuda_has_resident_state(dst->Bg) != 0;
+#elif USE_CUDA_BSSN
+  return bssn_cuda_has_resident_state(dst->Bg) != 0;
+#else
+  return false;
+#endif
+}
+
+bool cuda_direct_pack_segment(double *buffer,
+                              const Parallel::gridseg *src,
+                              const Parallel::gridseg *dst,
+                              int state_count)
+{
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  if (state_count != Z4C_CUDA_STATE_COUNT)
+    return false;
+#elif USE_CUDA_BSSN
+  if (state_count <= 0 || state_count > BSSN_CUDA_STATE_COUNT)
+    return false;
+#else
+  return false;
+#endif
+  const double t0 = sync_profile_enabled() ? MPI_Wtime() : 0.0;
+  const int i0 = cuda_seg_begin(dst, src->Bg, 0);
+  const int j0 = cuda_seg_begin(dst, src->Bg, 1);
+  const int k0 = cuda_seg_begin(dst, src->Bg, 2);
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  const bool ok = z4c_cuda_pack_state_batch_to_host_buffer(src->Bg, state_count, buffer, src->Bg->shape,
+                                                           i0, j0, k0,
+                                                           dst->shape[0], dst->shape[1], dst->shape[2]) == 0;
+#else
+  const bool ok = bssn_cuda_pack_state_batch_to_host_buffer(src->Bg, state_count, buffer, src->Bg->shape,
+                                                            i0, j0, k0,
+                                                            dst->shape[0], dst->shape[1], dst->shape[2]) == 0;
+#endif
+  if (sync_profile_enabled())
+    sync_profile_stats().direct_pack_sec += MPI_Wtime() - t0;
+  return ok;
+}
+
+bool cuda_direct_unpack_segment(double *buffer,
+                                const Parallel::gridseg *dst,
+                                int state_count)
+{
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  if (state_count != Z4C_CUDA_STATE_COUNT)
+    return false;
+#elif USE_CUDA_BSSN
+  if (state_count <= 0 || state_count > BSSN_CUDA_STATE_COUNT)
+    return false;
+#else
+  return false;
+#endif
+  const double t0 = sync_profile_enabled() ? MPI_Wtime() : 0.0;
+  const int i0 = cuda_seg_begin(dst, dst->Bg, 0);
+  const int j0 = cuda_seg_begin(dst, dst->Bg, 1);
+  const int k0 = cuda_seg_begin(dst, dst->Bg, 2);
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  const bool ok = z4c_cuda_unpack_state_batch_from_host_buffer(dst->Bg, state_count, buffer, dst->Bg->shape,
+                                                               i0, j0, k0,
+                                                               dst->shape[0], dst->shape[1], dst->shape[2]) == 0;
+#else
+  const bool ok = bssn_cuda_unpack_state_batch_from_host_buffer(dst->Bg, state_count, buffer, dst->Bg->shape,
+                                                                i0, j0, k0,
+                                                                dst->shape[0], dst->shape[1], dst->shape[2]) == 0;
+#endif
+  if (sync_profile_enabled())
+    sync_profile_stats().direct_unpack_sec += MPI_Wtime() - t0;
+  return ok;
+}
+#endif
+
+} // namespace
 
 int Parallel::partition1(int &nx, int split_size, int min_width, int cpusize, int shape) // special for 1 diemnsion
 {
@@ -3732,15 +3999,8 @@ int Parallel::data_packer(double *data, MyList<Parallel::gridseg> *src, MyList<P
 
   MyList<var> *varls, *varld;
 
-  varls = VarLists;
+  const int state_count = cuda_state_var_count(VarLists, VarListd);
-  varld = VarListd;
+  if (state_count < 0)
-  while (varls && varld)
-  {
-    varls = varls->next;
-    varld = varld->next;
-  }
-
-  if (varls || varld)
   {
     cout << "error in short data packer, var lists does not match." << endl;
     MPI_Abort(MPI_COMM_WORLD, 1);
@@ -3761,10 +4021,36 @@ int Parallel::data_packer(double *data, MyList<Parallel::gridseg> *src, MyList<P
     {
       varls = VarLists;
       varld = VarListd;
+      int state_idx = 0;
       while (varls && varld)
       {
         if (data)
         {
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+          bool handled_by_cuda = false;
+          if (dir == PACK && cuda_state_count_direct_supported(state_count) &&
+              cuda_can_direct_pack(src->data, dst->data, type))
+          {
+            handled_by_cuda = cuda_direct_pack_segment(data + size_out, src->data, dst->data, state_count);
+            if (!handled_by_cuda)
+            {
+              cout << "Parallel::data_packer: CUDA direct pack failed." << endl;
+              MPI_Abort(MPI_COMM_WORLD, 1);
+            }
+          }
+          else if (dir == UNPACK && cuda_state_count_direct_supported(state_count) &&
+                   cuda_can_direct_unpack(dst->data, type))
+          {
+            handled_by_cuda = cuda_direct_unpack_segment(data + size_out, dst->data, state_count);
+            if (!handled_by_cuda)
+            {
+              cout << "Parallel::data_packer: CUDA direct unpack failed." << endl;
+              MPI_Abort(MPI_COMM_WORLD, 1);
+            }
+          }
+          if (!handled_by_cuda)
+          {
+#endif
           if (dir == PACK)
             switch (type)
             {
@@ -3788,10 +4074,24 @@ int Parallel::data_packer(double *data, MyList<Parallel::gridseg> *src, MyList<P
             f_copy(DIM, dst->data->Bg->bbox, dst->data->Bg->bbox + dim, dst->data->Bg->shape, dst->data->Bg->fgfs[varld->data->sgfn],
                    dst->data->llb, dst->data->uub, dst->data->shape, data + size_out,
                    dst->data->llb, dst->data->uub);
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+          }
+          else
+          {
+            size_out += (state_count - 1) * dst->data->shape[0] * dst->data->shape[1] * dst->data->shape[2];
+            while (varls->next && varld->next)
+            {
+              varls = varls->next;
+              varld = varld->next;
+              ++state_idx;
+            }
+          }
+#endif
         }
         size_out += dst->data->shape[0] * dst->data->shape[1] * dst->data->shape[2];
         varls = varls->next;
         varld = varld->next;
+        ++state_idx;
       }
     }
     dst = dst->next;
@@ -4319,7 +4619,8 @@ void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmet
 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),
+      send_buf_caps(0), recv_buf_caps(0), send_buf_pinned(0), recv_buf_pinned(0),
+      reqs(0), stats(0), max_reqs(0),
       lengths_valid(false), tc_req_node(0), tc_req_is_recv(0), tc_completed(0)
 {
 }
@@ -4352,11 +4653,27 @@ void Parallel::SyncCache::destroy()
   if (recv_buf_caps) delete[] recv_buf_caps;
   for (int i = 0; i < cpusize; i++)
   {
-    if (send_bufs && send_bufs[i]) delete[] send_bufs[i];
+    if (send_bufs && send_bufs[i])
-    if (recv_bufs && recv_bufs[i]) delete[] recv_bufs[i];
+    {
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+      free_comm_buffer(send_bufs[i], send_buf_pinned[i]);
+#else
+      delete[] send_bufs[i];
+#endif
+    }
+    if (recv_bufs && recv_bufs[i])
+    {
+#if USE_CUDA_BSSN || USE_CUDA_Z4C
+      free_comm_buffer(recv_bufs[i], recv_buf_pinned[i]);
+#else
+      delete[] recv_bufs[i];
+#endif
+    }
   }
   if (send_bufs) delete[] send_bufs;
   if (recv_bufs) delete[] recv_bufs;
+  if (send_buf_pinned) delete[] send_buf_pinned;
+  if (recv_buf_pinned) delete[] recv_buf_pinned;
   if (reqs) delete[] reqs;
   if (stats) delete[] stats;
   if (tc_req_node) delete[] tc_req_node;
@@ -4366,6 +4683,7 @@ void Parallel::SyncCache::destroy()
   send_lengths = recv_lengths = 0;
   send_buf_caps = recv_buf_caps = 0;
   send_bufs = recv_bufs = 0;
+  send_buf_pinned = recv_buf_pinned = 0;
   reqs = 0; stats = 0;
   tc_req_node = 0; tc_req_is_recv = 0; tc_completed = 0;
   cpusize = 0; max_reqs = 0;
@@ -4396,12 +4714,7 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
     cache.recv_lengths[node] = rlength;
     if (rlength > 0)
     {
-      if (rlength > cache.recv_buf_caps[node])
+      ensure_comm_buffer(cache.recv_bufs, cache.recv_buf_pinned, cache.recv_buf_caps, node, rlength);
-      {
-        if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-        cache.recv_bufs[node] = new double[rlength];
-        cache.recv_buf_caps[node] = rlength;
-      }
       MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 1;
@@ -4415,12 +4728,7 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
   cache.recv_lengths[myrank] = self_len;
   if (self_len > 0)
   {
-    if (self_len > cache.recv_buf_caps[myrank])
+    ensure_comm_buffer(cache.recv_bufs, cache.recv_buf_pinned, cache.recv_buf_caps, myrank, self_len);
-    {
-      if (cache.recv_bufs[myrank]) delete[] cache.recv_bufs[myrank];
-      cache.recv_bufs[myrank] = new double[self_len];
-      cache.recv_buf_caps[myrank] = self_len;
-    }
     data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
   }
 
@@ -4433,12 +4741,7 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
     cache.send_lengths[node] = slength;
     if (slength > 0)
     {
-      if (slength > cache.send_buf_caps[node])
+      ensure_comm_buffer(cache.send_bufs, cache.send_buf_pinned, cache.send_buf_caps, node, slength);
-      {
-        if (cache.send_bufs[node]) delete[] cache.send_bufs[node];
-        cache.send_bufs[node] = new double[slength];
-        cache.send_buf_caps[node] = slength;
-      }
       data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
       MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;
@@ -4471,82 +4774,57 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
   if (self_len > 0)
     data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
 }
-void Parallel::Sync_cached(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, SyncCache &cache)
+void Parallel::Sync_ensure_cache(MyList<Patch> *PatL, int Symmetry, SyncCache &cache)
 {
-  if (!cache.valid)
+  if (cache.valid)
+    return;
+
+  int cpusize;
+  MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
+  cache.cpusize = cpusize;
+
+  if (!cache.combined_src)
   {
-    int cpusize;
+    cache.combined_src = new MyList<Parallel::gridseg> *[cpusize];
-    MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
+    cache.combined_dst = new MyList<Parallel::gridseg> *[cpusize];
-    cache.cpusize = cpusize;
+    cache.send_lengths = new int[cpusize];
-
+    cache.recv_lengths = new int[cpusize];
-    // Allocate cache arrays if needed
+    cache.send_bufs = new double *[cpusize];
-    if (!cache.combined_src)
+    cache.recv_bufs = new double *[cpusize];
+    cache.send_buf_caps = new int[cpusize];
+    cache.recv_buf_caps = new int[cpusize];
+    cache.send_buf_pinned = new unsigned char[cpusize];
+    cache.recv_buf_pinned = new unsigned char[cpusize];
+    for (int i = 0; i < cpusize; i++)
     {
-      cache.combined_src = new MyList<Parallel::gridseg> *[cpusize];
+      cache.send_bufs[i] = cache.recv_bufs[i] = 0;
-      cache.combined_dst = new MyList<Parallel::gridseg> *[cpusize];
+      cache.send_buf_caps[i] = cache.recv_buf_caps[i] = 0;
-      cache.send_lengths = new int[cpusize];
+      cache.send_buf_pinned[i] = cache.recv_buf_pinned[i] = 0;
-      cache.recv_lengths = new int[cpusize];
-      cache.send_bufs = new double *[cpusize];
-      cache.recv_bufs = new double *[cpusize];
-      cache.send_buf_caps = new int[cpusize];
-      cache.recv_buf_caps = new int[cpusize];
-      for (int i = 0; i < cpusize; i++)
-      {
-        cache.send_bufs[i] = cache.recv_bufs[i] = 0;
-        cache.send_buf_caps[i] = cache.recv_buf_caps[i] = 0;
-      }
-      cache.max_reqs = 2 * cpusize;
-      cache.reqs = new MPI_Request[cache.max_reqs];
-      cache.stats = new MPI_Status[cache.max_reqs];
-      cache.tc_req_node = new int[cache.max_reqs];
-      cache.tc_req_is_recv = new int[cache.max_reqs];
-      cache.tc_completed = new int[cache.max_reqs];
     }
+    cache.max_reqs = 2 * cpusize;
+    cache.reqs = new MPI_Request[cache.max_reqs];
+    cache.stats = new MPI_Status[cache.max_reqs];
+    cache.tc_req_node = new int[cache.max_reqs];
+    cache.tc_req_is_recv = new int[cache.max_reqs];
+    cache.tc_completed = new int[cache.max_reqs];
+  }
 
+  for (int node = 0; node < cpusize; node++)
+  {
+    cache.combined_src[node] = cache.combined_dst[node] = 0;
+    cache.send_lengths[node] = cache.recv_lengths[node] = 0;
+  }
+
+  MyList<Patch> *Pp = PatL;
+  while (Pp)
+  {
+    Patch *Pat = Pp->data;
+    MyList<Parallel::gridseg> *dst_ghost = build_ghost_gsl(Pat);
     for (int node = 0; node < cpusize; node++)
     {
-      cache.combined_src[node] = cache.combined_dst[node] = 0;
+      MyList<Parallel::gridseg> *src_owned = build_owned_gsl0(Pat, node);
-      cache.send_lengths[node] = cache.recv_lengths[node] = 0;
-    }
-
-    // Build intra-patch segments (same as Sync_merged Phase A)
-    MyList<Patch> *Pp = PatL;
-    while (Pp)
-    {
-      Patch *Pat = Pp->data;
-      MyList<Parallel::gridseg> *dst_ghost = build_ghost_gsl(Pat);
-      for (int node = 0; node < cpusize; node++)
-      {
-        MyList<Parallel::gridseg> *src_owned = build_owned_gsl0(Pat, node);
-        MyList<Parallel::gridseg> *tsrc = 0, *tdst = 0;
-        build_gstl(src_owned, dst_ghost, &tsrc, &tdst);
-        if (tsrc)
-        {
-          if (cache.combined_src[node])
-            cache.combined_src[node]->catList(tsrc);
-          else
-            cache.combined_src[node] = tsrc;
-        }
-        if (tdst)
-        {
-          if (cache.combined_dst[node])
-            cache.combined_dst[node]->catList(tdst);
-          else
-            cache.combined_dst[node] = tdst;
-        }
-        if (src_owned) src_owned->destroyList();
-      }
-      if (dst_ghost) dst_ghost->destroyList();
-      Pp = Pp->next;
-    }
-
-    // Build inter-patch segments (same as Sync_merged Phase B)
-    MyList<Parallel::gridseg> *dst_buffer = build_buffer_gsl(PatL);
-    for (int node = 0; node < cpusize; node++)
-    {
-      MyList<Parallel::gridseg> *src_owned = build_owned_gsl(PatL, node, 5, Symmetry);
       MyList<Parallel::gridseg> *tsrc = 0, *tdst = 0;
-      build_gstl(src_owned, dst_buffer, &tsrc, &tdst);
+      build_gstl(src_owned, dst_ghost, &tsrc, &tdst);
       if (tsrc)
       {
         if (cache.combined_src[node])
@@ -4563,11 +4841,40 @@ void Parallel::Sync_cached(MyList<Patch> *PatL, MyList<var> *VarList, int Symmet
       }
       if (src_owned) src_owned->destroyList();
     }
-    if (dst_buffer) dst_buffer->destroyList();
+    if (dst_ghost) dst_ghost->destroyList();
-
+    Pp = Pp->next;
-    cache.valid = true;
   }
 
+  MyList<Parallel::gridseg> *dst_buffer = build_buffer_gsl(PatL);
+  for (int node = 0; node < cpusize; node++)
+  {
+    MyList<Parallel::gridseg> *src_owned = build_owned_gsl(PatL, node, 5, Symmetry);
+    MyList<Parallel::gridseg> *tsrc = 0, *tdst = 0;
+    build_gstl(src_owned, dst_buffer, &tsrc, &tdst);
+    if (tsrc)
+    {
+      if (cache.combined_src[node])
+        cache.combined_src[node]->catList(tsrc);
+      else
+        cache.combined_src[node] = tsrc;
+    }
+    if (tdst)
+    {
+      if (cache.combined_dst[node])
+        cache.combined_dst[node]->catList(tdst);
+      else
+        cache.combined_dst[node] = tdst;
+    }
+    if (src_owned) src_owned->destroyList();
+  }
+  if (dst_buffer) dst_buffer->destroyList();
+
+  cache.valid = true;
+}
+void Parallel::Sync_cached(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, SyncCache &cache)
+{
+  Sync_ensure_cache(PatL, Symmetry, cache);
+
   // Use cached lists with buffer-reusing transfer
   transfer_cached(cache.combined_src, cache.combined_dst, VarList, VarList, Symmetry, cache);
 }
@@ -4575,98 +4882,8 @@ void Parallel::Sync_cached(MyList<Patch> *PatL, MyList<var> *VarList, int Symmet
 void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry,
                           SyncCache &cache, AsyncSyncState &state)
 {
-  // Ensure cache is built
+  const double t_start = sync_profile_enabled() ? MPI_Wtime() : 0.0;
-  if (!cache.valid)
+  Sync_ensure_cache(PatL, Symmetry, cache);
-  {
-    // Build cache (same logic as Sync_cached)
-    int cpusize;
-    MPI_Comm_size(MPI_COMM_WORLD, &cpusize);
-    cache.cpusize = cpusize;
-
-    if (!cache.combined_src)
-    {
-      cache.combined_src = new MyList<Parallel::gridseg> *[cpusize];
-      cache.combined_dst = new MyList<Parallel::gridseg> *[cpusize];
-      cache.send_lengths = new int[cpusize];
-      cache.recv_lengths = new int[cpusize];
-      cache.send_bufs = new double *[cpusize];
-      cache.recv_bufs = new double *[cpusize];
-      cache.send_buf_caps = new int[cpusize];
-      cache.recv_buf_caps = new int[cpusize];
-      for (int i = 0; i < cpusize; i++)
-      {
-        cache.send_bufs[i] = cache.recv_bufs[i] = 0;
-        cache.send_buf_caps[i] = cache.recv_buf_caps[i] = 0;
-      }
-      cache.max_reqs = 2 * cpusize;
-      cache.reqs = new MPI_Request[cache.max_reqs];
-      cache.stats = new MPI_Status[cache.max_reqs];
-      cache.tc_req_node = new int[cache.max_reqs];
-      cache.tc_req_is_recv = new int[cache.max_reqs];
-      cache.tc_completed = new int[cache.max_reqs];
-    }
-
-    for (int node = 0; node < cpusize; node++)
-    {
-      cache.combined_src[node] = cache.combined_dst[node] = 0;
-      cache.send_lengths[node] = cache.recv_lengths[node] = 0;
-    }
-
-    MyList<Patch> *Pp = PatL;
-    while (Pp)
-    {
-      Patch *Pat = Pp->data;
-      MyList<Parallel::gridseg> *dst_ghost = build_ghost_gsl(Pat);
-      for (int node = 0; node < cpusize; node++)
-      {
-        MyList<Parallel::gridseg> *src_owned = build_owned_gsl0(Pat, node);
-        MyList<Parallel::gridseg> *tsrc = 0, *tdst = 0;
-        build_gstl(src_owned, dst_ghost, &tsrc, &tdst);
-        if (tsrc)
-        {
-          if (cache.combined_src[node])
-            cache.combined_src[node]->catList(tsrc);
-          else
-            cache.combined_src[node] = tsrc;
-        }
-        if (tdst)
-        {
-          if (cache.combined_dst[node])
-            cache.combined_dst[node]->catList(tdst);
-          else
-            cache.combined_dst[node] = tdst;
-        }
-        if (src_owned) src_owned->destroyList();
-      }
-      if (dst_ghost) dst_ghost->destroyList();
-      Pp = Pp->next;
-    }
-
-    MyList<Parallel::gridseg> *dst_buffer = build_buffer_gsl(PatL);
-    for (int node = 0; node < cpusize; node++)
-    {
-      MyList<Parallel::gridseg> *src_owned = build_owned_gsl(PatL, node, 5, Symmetry);
-      MyList<Parallel::gridseg> *tsrc = 0, *tdst = 0;
-      build_gstl(src_owned, dst_buffer, &tsrc, &tdst);
-      if (tsrc)
-      {
-        if (cache.combined_src[node])
-          cache.combined_src[node]->catList(tsrc);
-        else
-          cache.combined_src[node] = tsrc;
-      }
-      if (tdst)
-      {
-        if (cache.combined_dst[node])
-          cache.combined_dst[node]->catList(tdst);
-        else
-          cache.combined_dst[node] = tdst;
-      }
-      if (src_owned) src_owned->destroyList();
-    }
-    if (dst_buffer) dst_buffer->destroyList();
-    cache.valid = true;
-  }
 
   // Now pack and post async MPI operations
   int myrank;
@@ -4683,6 +4900,27 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
   MyList<Parallel::gridseg> **src = cache.combined_src;
   MyList<Parallel::gridseg> **dst = cache.combined_dst;
 
+  for (int node = 0; node < cpusize; node++)
+  {
+    if (node == myrank)
+      continue;
+    int rlength;
+    if (!cache.lengths_valid) {
+      rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList, VarList, Symmetry);
+      cache.recv_lengths[node] = rlength;
+    } else {
+      rlength = cache.recv_lengths[node];
+    }
+    if (rlength > 0)
+    {
+      ensure_comm_buffer(cache.recv_bufs, cache.recv_buf_pinned, cache.recv_buf_caps, node, rlength);
+      state.req_node[state.req_no] = node;
+      state.req_is_recv[state.req_no] = 1;
+      state.pending_recv++;
+      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
+    }
+  }
+
   for (int node = 0; node < cpusize; node++)
   {
     if (node == myrank)
@@ -4696,12 +4934,7 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
       }
       if (length > 0)
       {
-        if (length > cache.recv_buf_caps[node])
+        ensure_comm_buffer(cache.recv_bufs, cache.recv_buf_pinned, cache.recv_buf_caps, node, length);
-        {
-          if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-          cache.recv_bufs[node] = new double[length];
-          cache.recv_buf_caps[node] = length;
-        }
         data_packer(cache.recv_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
       }
     }
@@ -4716,40 +4949,21 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
       }
       if (slength > 0)
       {
-        if (slength > cache.send_buf_caps[node])
+        ensure_comm_buffer(cache.send_bufs, cache.send_buf_pinned, cache.send_buf_caps, node, slength);
-        {
-          if (cache.send_bufs[node]) delete[] cache.send_bufs[node];
-          cache.send_bufs[node] = new double[slength];
-          cache.send_buf_caps[node] = slength;
-        }
         data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
         state.req_node[state.req_no] = node;
         state.req_is_recv[state.req_no] = 0;
         MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
       }
-      int rlength;
-      if (!cache.lengths_valid) {
-        rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList, VarList, Symmetry);
-        cache.recv_lengths[node] = rlength;
-      } else {
-        rlength = cache.recv_lengths[node];
-      }
-      if (rlength > 0)
-      {
-        if (rlength > cache.recv_buf_caps[node])
-        {
-          if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-          cache.recv_bufs[node] = new double[rlength];
-          cache.recv_buf_caps[node] = rlength;
-        }
-        state.req_node[state.req_no] = node;
-        state.req_is_recv[state.req_no] = 1;
-        state.pending_recv++;
-        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
-      }
     }
   }
   cache.lengths_valid = true;
+  if (sync_profile_enabled())
+  {
+    SyncProfileStats &stats = sync_profile_stats();
+    stats.start_calls++;
+    stats.start_sec += MPI_Wtime() - t_start;
+  }
 }
 // Sync_finish: progressive unpack as receives complete, then wait for sends
 void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
@@ -4762,6 +4976,8 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
   MyList<Parallel::gridseg> **src = cache.combined_src;
   MyList<Parallel::gridseg> **dst = cache.combined_dst;
+  const double t_finish = sync_profile_enabled() ? MPI_Wtime() : 0.0;
+  double wait_sec = 0.0;
 
   // Unpack local data first (no MPI needed)
   if (cache.recv_bufs[myrank] && cache.recv_lengths[myrank] > 0)
@@ -4771,15 +4987,17 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
   if (state.pending_recv > 0 && state.req_no > 0)
   {
     int pending = state.pending_recv;
-    int *completed = new int[cache.max_reqs];
     while (pending > 0)
     {
       int outcount = 0;
-      MPI_Waitsome(state.req_no, cache.reqs, &outcount, completed, cache.stats);
+      const double t_wait = sync_profile_enabled() ? MPI_Wtime() : 0.0;
+      MPI_Waitsome(state.req_no, cache.reqs, &outcount, cache.tc_completed, cache.stats);
+      if (sync_profile_enabled())
+        wait_sec += MPI_Wtime() - t_wait;
       if (outcount == MPI_UNDEFINED) break;
       for (int i = 0; i < outcount; i++)
       {
-        int idx = completed[i];
+        int idx = cache.tc_completed[i];
         if (idx >= 0 && state.req_is_recv[idx])
         {
           int recv_node = state.req_node[idx];
@@ -4788,15 +5006,28 @@ void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
         }
       }
     }
-    delete[] completed;
   }
 
   // Wait for remaining sends
-  if (state.req_no > 0) MPI_Waitall(state.req_no, cache.reqs, cache.stats);
+  if (state.req_no > 0)
+  {
+    const double t_wait = sync_profile_enabled() ? MPI_Wtime() : 0.0;
+    MPI_Waitall(state.req_no, cache.reqs, cache.stats);
+    if (sync_profile_enabled())
+      wait_sec += MPI_Wtime() - t_wait;
+  }
 
   delete[] state.req_node; state.req_node = 0;
   delete[] state.req_is_recv; state.req_is_recv = 0;
   state.active = false;
+  if (sync_profile_enabled())
+  {
+    SyncProfileStats &stats = sync_profile_stats();
+    stats.finish_calls++;
+    stats.finish_sec += MPI_Wtime() - t_finish;
+    stats.wait_sec += wait_sec;
+    sync_profile_maybe_log();
+  }
 }
 // collect buffer grid segments or blocks for the periodic boundary condition of given patch
 // ---------------------------------------------------
@@ -5924,10 +6155,13 @@ void Parallel::Restrict_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
       cache.recv_bufs = new double *[cpusize];
       cache.send_buf_caps = new int[cpusize];
       cache.recv_buf_caps = new int[cpusize];
+      cache.send_buf_pinned = new unsigned char[cpusize];
+      cache.recv_buf_pinned = new unsigned char[cpusize];
       for (int i = 0; i < cpusize; i++)
       {
         cache.send_bufs[i] = cache.recv_bufs[i] = 0;
         cache.send_buf_caps[i] = cache.recv_buf_caps[i] = 0;
+        cache.send_buf_pinned[i] = cache.recv_buf_pinned[i] = 0;
       }
       cache.max_reqs = 2 * cpusize;
       cache.reqs = new MPI_Request[cache.max_reqs];
@@ -5973,10 +6207,13 @@ void Parallel::OutBdLow2Hi_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
       cache.recv_bufs = new double *[cpusize];
       cache.send_buf_caps = new int[cpusize];
       cache.recv_buf_caps = new int[cpusize];
+      cache.send_buf_pinned = new unsigned char[cpusize];
+      cache.recv_buf_pinned = new unsigned char[cpusize];
       for (int i = 0; i < cpusize; i++)
       {
         cache.send_bufs[i] = cache.recv_bufs[i] = 0;
         cache.send_buf_caps[i] = cache.recv_buf_caps[i] = 0;
+        cache.send_buf_pinned[i] = cache.recv_buf_pinned[i] = 0;
       }
       cache.max_reqs = 2 * cpusize;
       cache.reqs = new MPI_Request[cache.max_reqs];
@@ -6022,10 +6259,13 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
       cache.recv_bufs = new double *[cpusize];
       cache.send_buf_caps = new int[cpusize];
       cache.recv_buf_caps = new int[cpusize];
+      cache.send_buf_pinned = new unsigned char[cpusize];
+      cache.recv_buf_pinned = new unsigned char[cpusize];
       for (int i = 0; i < cpusize; i++)
       {
         cache.send_bufs[i] = cache.recv_bufs[i] = 0;
         cache.send_buf_caps[i] = cache.recv_buf_caps[i] = 0;
+        cache.send_buf_pinned[i] = cache.recv_buf_pinned[i] = 0;
       }
       cache.max_reqs = 2 * cpusize;
       cache.reqs = new MPI_Request[cache.max_reqs];
@@ -6068,12 +6308,7 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
     cache.recv_lengths[node] = rlength;
     if (rlength > 0)
     {
-      if (rlength > cache.recv_buf_caps[node])
+      ensure_comm_buffer(cache.recv_bufs, cache.recv_buf_pinned, cache.recv_buf_caps, node, rlength);
-      {
-        if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-        cache.recv_bufs[node] = new double[rlength];
-        cache.recv_buf_caps[node] = rlength;
-      }
       MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;
       req_is_recv[req_no] = 1;
@@ -6087,12 +6322,7 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
   cache.recv_lengths[myrank] = self_len;
   if (self_len > 0)
   {
-    if (self_len > cache.recv_buf_caps[myrank])
+    ensure_comm_buffer(cache.recv_bufs, cache.recv_buf_pinned, cache.recv_buf_caps, myrank, self_len);
-    {
-      if (cache.recv_bufs[myrank]) delete[] cache.recv_bufs[myrank];
-      cache.recv_bufs[myrank] = new double[self_len];
-      cache.recv_buf_caps[myrank] = self_len;
-    }
     data_packermix(cache.recv_bufs[myrank], cache.combined_src[myrank], cache.combined_dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
   }
 
@@ -6105,12 +6335,7 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
     cache.send_lengths[node] = slength;
     if (slength > 0)
     {
-      if (slength > cache.send_buf_caps[node])
+      ensure_comm_buffer(cache.send_bufs, cache.send_buf_pinned, cache.send_buf_caps, node, slength);
-      {
-        if (cache.send_bufs[node]) delete[] cache.send_bufs[node];
-        cache.send_bufs[node] = new double[slength];
-        cache.send_buf_caps[node] = slength;
-      }
       data_packermix(cache.send_bufs[node], cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
       MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
       req_node[req_no] = node;

AMSS_NCKU_source/Parallel.h

@@ -104,6 +104,8 @@ namespace Parallel
     double **recv_bufs;
     int *send_buf_caps;
     int *recv_buf_caps;
+    unsigned char *send_buf_pinned;
+    unsigned char *recv_buf_pinned;
     MPI_Request *reqs;
     MPI_Status *stats;
     int max_reqs;
@@ -117,6 +119,7 @@ namespace Parallel
   };
 
   void Sync_cached(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, SyncCache &cache);
+  void Sync_ensure_cache(MyList<Patch> *PatL, int Symmetry, SyncCache &cache);
   void transfer_cached(MyList<gridseg> **src, MyList<gridseg> **dst,
                        MyList<var> *VarList1, MyList<var> *VarList2,
                        int Symmetry, SyncCache &cache);

AMSS_NCKU_source/Z4c_class.C

@@ -28,6 +28,14 @@ using namespace std;
 #include "kodiss.h"
 #include "parameters.h"
 
+#ifndef USE_CUDA_Z4C
+#define USE_CUDA_Z4C 0
+#endif
+
+#if USE_CUDA_Z4C && (ABEtype == 2)
+#include "z4c_rhs_cuda.h"
+#endif
+
 #ifdef With_AHF
 #include "derivatives.h"
 #include "myglobal.h"
@@ -170,8 +178,550 @@ Z4c_class::~Z4c_class()
 #ifndef CPBC
 // for sommerfeld boundary
 
+#if USE_CUDA_Z4C && (ABEtype == 2)
+#ifdef WithShell
+#error "USE_CUDA_Z4C resident path currently supports Cartesian non-shell Z4C only"
+#endif
+#if (MRBD == 2)
+#error "USE_CUDA_Z4C resident path does not support MRBD == 2"
+#endif
+
+namespace {
+
+static const int k_z4c_cuda_bh_state_indices[3] = {18, 19, 20};
+
+bool fill_z4c_cuda_views(Block *cg, MyList<var> *vars,
+                         double **host_views,
+                         double *propspeeds = 0,
+                         double *soa_flat = 0)
+{
+  int idx = 0;
+  while (vars && idx < Z4C_CUDA_STATE_COUNT)
+  {
+    host_views[idx] = cg->fgfs[vars->data->sgfn];
+    if (propspeeds)
+      propspeeds[idx] = vars->data->propspeed;
+    if (soa_flat)
+    {
+      soa_flat[3 * idx + 0] = vars->data->SoA[0];
+      soa_flat[3 * idx + 1] = vars->data->SoA[1];
+      soa_flat[3 * idx + 2] = vars->data->SoA[2];
+    }
+    vars = vars->next;
+    ++idx;
+  }
+  return idx == Z4C_CUDA_STATE_COUNT && vars == 0;
+}
+
+void z4c_cuda_download_level_state(MyList<Patch> *PatL, MyList<var> *vars, int myrank, bool release_ctx)
+{
+  MyList<Patch> *Pp = PatL;
+  while (Pp)
+  {
+    MyList<Block> *BP = Pp->data->blb;
+    while (BP)
+    {
+      Block *cg = BP->data;
+      if (myrank == cg->rank && z4c_cuda_has_resident_state(cg))
+      {
+        double *state_out[Z4C_CUDA_STATE_COUNT];
+        if (!fill_z4c_cuda_views(cg, vars, state_out))
+        {
+          cout << "CUDA Z4C state list mismatch on resident state download" << endl;
+          MPI_Abort(MPI_COMM_WORLD, 1);
+        }
+        if (z4c_cuda_download_resident_state(cg, cg->shape, state_out))
+        {
+          cout << "CUDA Z4C resident state download failed" << endl;
+          MPI_Abort(MPI_COMM_WORLD, 1);
+        }
+        if (release_ctx)
+          z4c_cuda_release_step_ctx(cg);
+      }
+      if (BP == Pp->data->ble)
+        break;
+      BP = BP->next;
+    }
+    Pp = Pp->next;
+  }
+}
+
+bool z4c_cuda_patch_contains_point(Patch *patch, const double *point)
+{
+  if (!patch)
+    return false;
+  for (int d = 0; d < dim; d++)
+  {
+    const double h = patch->getdX(d);
+    const double lo = patch->bbox[d] + patch->lli[d] * h;
+    const double hi = patch->bbox[dim + d] - patch->uui[d] * h;
+    if (point[d] < lo || point[d] > hi)
+      return false;
+  }
+  return true;
+}
+
+bool z4c_cuda_point_in_block(Patch *patch, Block *block,
+                             const double *point, const double *DH)
+{
+  if (!patch || !block)
+    return false;
+  for (int d = 0; d < dim; d++)
+  {
+    double llb;
+    double uub;
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+    llb = (feq(block->bbox[d], patch->bbox[d], DH[d] / 2))
+              ? block->bbox[d] + patch->lli[d] * DH[d]
+              : block->bbox[d] + (ghost_width - 0.5) * DH[d];
+    uub = (feq(block->bbox[dim + d], patch->bbox[dim + d], DH[d] / 2))
+              ? block->bbox[dim + d] - patch->uui[d] * DH[d]
+              : block->bbox[dim + d] - (ghost_width - 0.5) * DH[d];
+#else
+#ifdef Cell
+    llb = (feq(block->bbox[d], patch->bbox[d], DH[d] / 2))
+              ? block->bbox[d] + patch->lli[d] * DH[d]
+              : block->bbox[d] + ghost_width * DH[d];
+    uub = (feq(block->bbox[dim + d], patch->bbox[dim + d], DH[d] / 2))
+              ? block->bbox[dim + d] - patch->uui[d] * DH[d]
+              : block->bbox[dim + d] - ghost_width * DH[d];
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+    if (point[d] - llb < -DH[d] / 2 || point[d] - uub > DH[d] / 2)
+      return false;
+  }
+  return true;
+}
+
+int z4c_cuda_interp_tile_start(const double *coords, int n, double x, double dx, int ordn)
+{
+  if (!coords || n <= ordn)
+    return 0;
+  int cxi = int((x - coords[0]) / dx + 0.4) + 1;
+  int start = cxi - ordn / 2;
+  if (start < 0)
+    start = 0;
+  const int max_start = n - ordn;
+  if (start > max_start)
+    start = max_start;
+  return start;
+}
+
+bool z4c_cuda_interp_bh_point_resident(MyList<Patch> *PatL,
+                                       int myrank,
+                                       const double *point,
+                                       var *forx, var *fory, var *forz,
+                                       int Symmetry,
+                                       double *shellf)
+{
+  const int ordn = 2 * ghost_width;
+  int owner_rank = -1;
+
+  shellf[0] = shellf[1] = shellf[2] = 0.0;
+
+  MyList<Patch> *PL = PatL;
+  while (PL)
+  {
+    Patch *patch = PL->data;
+    if (!z4c_cuda_patch_contains_point(patch, point))
+    {
+      PL = PL->next;
+      continue;
+    }
+
+    double DH[dim];
+    for (int d = 0; d < dim; d++)
+      DH[d] = patch->getdX(d);
+
+    MyList<Block> *BP = patch->blb;
+    while (BP)
+    {
+      Block *block = BP->data;
+      if (z4c_cuda_point_in_block(patch, block, point, DH))
+      {
+        owner_rank = block->rank;
+        if (myrank == owner_rank)
+        {
+          int interp_ordn = ordn;
+          int interp_sym = Symmetry;
+          double x = point[0];
+          double y = point[1];
+          double z = point[2];
+
+          if (z4c_cuda_has_resident_state(block) &&
+              block->shape[0] >= ordn && block->shape[1] >= ordn && block->shape[2] >= ordn)
+          {
+            const int sx = ordn;
+            const int sy = ordn;
+            const int sz = ordn;
+            const int region_all = sx * sy * sz;
+            const int i0 = z4c_cuda_interp_tile_start(block->X[0], block->shape[0], x, DH[0], ordn);
+            const int j0 = z4c_cuda_interp_tile_start(block->X[1], block->shape[1], y, DH[1], ordn);
+            const int k0 = z4c_cuda_interp_tile_start(block->X[2], block->shape[2], z, DH[2], ordn);
+            double *packed_fields = new double[3 * region_all];
+            var *vars[3] = {forx, fory, forz};
+            for (int f = 0; f < 3; f++)
+            {
+              if (z4c_cuda_pack_state_region_to_host_buffer(block,
+                                                            k_z4c_cuda_bh_state_indices[f],
+                                                            packed_fields + f * region_all,
+                                                            block->shape,
+                                                            i0, j0, k0,
+                                                            sx, sy, sz) != 0)
+              {
+                delete[] packed_fields;
+                cout << "CUDA Z4C BH tile download failed" << endl;
+                MPI_Abort(MPI_COMM_WORLD, 1);
+              }
+              int tile_shape[3] = {sx, sy, sz};
+              f_global_interp(tile_shape,
+                              block->X[0] + i0,
+                              block->X[1] + j0,
+                              block->X[2] + k0,
+                              packed_fields + f * region_all,
+                              shellf[f],
+                              x, y, z,
+                              interp_ordn,
+                              vars[f]->SoA,
+                              interp_sym);
+            }
+            delete[] packed_fields;
+          }
+          else
+          {
+            f_global_interp(block->shape, block->X[0], block->X[1], block->X[2],
+                            block->fgfs[forx->sgfn], shellf[0],
+                            x, y, z, interp_ordn, forx->SoA, interp_sym);
+            f_global_interp(block->shape, block->X[0], block->X[1], block->X[2],
+                            block->fgfs[fory->sgfn], shellf[1],
+                            x, y, z, interp_ordn, fory->SoA, interp_sym);
+            f_global_interp(block->shape, block->X[0], block->X[1], block->X[2],
+                            block->fgfs[forz->sgfn], shellf[2],
+                            x, y, z, interp_ordn, forz->SoA, interp_sym);
+          }
+        }
+        break;
+      }
+      if (BP == patch->ble)
+        break;
+      BP = BP->next;
+    }
+
+    if (owner_rank >= 0)
+      break;
+    PL = PL->next;
+  }
+
+  if (owner_rank < 0)
+    return false;
+
+  MPI_Bcast(shellf, 3, MPI_DOUBLE, owner_rank, MPI_COMM_WORLD);
+  return true;
+}
+
+bool z4c_cuda_compute_porg_rhs_resident(cgh *GH,
+                                        int ilev,
+                                        int myrank,
+                                        int BH_num,
+                                        double **BH_PS,
+                                        double **BH_RHS,
+                                        var *forx, var *fory, var *forz,
+                                        int Symmetry)
+{
+  for (int n = 0; n < BH_num; n++)
+  {
+    double shellf[3] = {0.0, 0.0, 0.0};
+    int lev = ilev;
+    while (lev >= 0 &&
+           !z4c_cuda_interp_bh_point_resident(GH->PatL[lev], myrank, BH_PS[n],
+                                              forx, fory, forz, Symmetry, shellf))
+    {
+      --lev;
+    }
+    if (lev < 0)
+      return false;
+    BH_RHS[n][0] = -shellf[0];
+    BH_RHS[n][1] = -shellf[1];
+    BH_RHS[n][2] = -shellf[2];
+  }
+  return true;
+}
+
+} // namespace
+#endif
+
 void Z4c_class::Step(int lev, int YN)
 {
+#if USE_CUDA_Z4C && (ABEtype == 2)
+  double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
+#ifdef With_AHF
+  AH_Step_Find(lev, dT_lev);
+#endif
+  bool BB = fgt(PhysTime, StartTime, dT_lev / 2);
+  double ndeps = numepss;
+  if (lev < GH->movls)
+    ndeps = numepsb;
+  double TRK4 = PhysTime;
+  int iter_count = 0;
+  int pre = 0, cor = 1;
+  int ERROR = 0;
+
+  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 *state_in[Z4C_CUDA_STATE_COUNT];
+        double *state_out[Z4C_CUDA_STATE_COUNT];
+        double propspeed[Z4C_CUDA_STATE_COUNT];
+        double soa_flat[3 * Z4C_CUDA_STATE_COUNT];
+        if (!fill_z4c_cuda_views(cg, StateList, state_in, propspeed, soa_flat) ||
+            !fill_z4c_cuda_views(cg, SynchList_pre, state_out))
+        {
+          cout << "CUDA Z4C state list mismatch on predictor step" << endl;
+          MPI_Abort(MPI_COMM_WORLD, 1);
+        }
+        int apply_bam_bc = 0;
+#if (MRBD == 0)
+#if (SommerType == 0)
+        apply_bam_bc = (lev == 0) ? 1 : 0;
+#endif
+#elif (MRBD == 1)
+        apply_bam_bc = 1;
+#endif
+        int keep_resident_state = 1;
+        int apply_enforce_ga = 0;
+#if (AGM == 0)
+        apply_enforce_ga = 1;
+#endif
+        if (z4c_cuda_rk4_substep(cg,
+                                 cg->shape, cg->X[0], cg->X[1], cg->X[2],
+                                 state_in, state_out,
+                                 propspeed, soa_flat, Pp->data->bbox,
+                                 dT_lev, TRK4, iter_count, apply_bam_bc,
+                                 Symmetry, lev, ndeps, pre,
+                                 keep_resident_state, apply_enforce_ga, chitiny))
+        {
+          cout << "CUDA Z4C predictor substep failed in domain: ("
+               << cg->bbox[0] << ":" << cg->bbox[3] << ","
+               << cg->bbox[1] << ":" << cg->bbox[4] << ","
+               << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+          ERROR = 1;
+        }
+      }
+      if (BP == Pp->data->ble)
+        break;
+      BP = BP->next;
+    }
+    Pp = Pp->next;
+  }
+
+  {
+    int erh = ERROR;
+    MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+  }
+  if (ERROR)
+  {
+    if (myrank == 0 && ErrorMonitor->outfile)
+      ErrorMonitor->outfile << "CUDA Z4C failed in predictor at t = " << PhysTime
+                            << ", lev = " << lev << endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+
+  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
+
+  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) && (LastAnas + dT_lev >= AnasTime))
+    z4c_cuda_download_level_state(GH->PatL[lev], SynchList_pre, myrank, false);
+  if (lev == a_lev)
+    AnalysisStuff(lev, dT_lev);
+
+  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 *state_in[Z4C_CUDA_STATE_COUNT];
+          double *state_out[Z4C_CUDA_STATE_COUNT];
+          double propspeed[Z4C_CUDA_STATE_COUNT];
+          double soa_flat[3 * Z4C_CUDA_STATE_COUNT];
+          if (!fill_z4c_cuda_views(cg, SynchList_pre, state_in, propspeed, soa_flat) ||
+              !fill_z4c_cuda_views(cg, SynchList_cor, state_out))
+          {
+            cout << "CUDA Z4C state list mismatch on corrector step" << endl;
+            MPI_Abort(MPI_COMM_WORLD, 1);
+          }
+          int apply_bam_bc = 0;
+#if (MRBD == 0)
+#if (SommerType == 0)
+          apply_bam_bc = (lev == 0) ? 1 : 0;
+#endif
+#elif (MRBD == 1)
+          apply_bam_bc = 1;
+#endif
+          int keep_resident_state = 1;
+          int apply_enforce_ga = 0;
+#if (AGM == 0)
+          apply_enforce_ga = 1;
+#elif (AGM == 1)
+          apply_enforce_ga = (iter_count == 3) ? 1 : 0;
+#endif
+          if (z4c_cuda_rk4_substep(cg,
+                                   cg->shape, cg->X[0], cg->X[1], cg->X[2],
+                                   state_in, state_out,
+                                   propspeed, soa_flat, Pp->data->bbox,
+                                   dT_lev, TRK4, iter_count, apply_bam_bc,
+                                   Symmetry, lev, ndeps, cor,
+                                   keep_resident_state, apply_enforce_ga, chitiny))
+          {
+            cout << "CUDA Z4C corrector substep failed in domain: ("
+                 << cg->bbox[0] << ":" << cg->bbox[3] << ","
+                 << cg->bbox[1] << ":" << cg->bbox[4] << ","
+                 << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
+            ERROR = 1;
+          }
+        }
+        if (BP == Pp->data->ble)
+          break;
+        BP = BP->next;
+      }
+      Pp = Pp->next;
+    }
+
+    {
+      int erh = ERROR;
+      MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+    }
+    if (ERROR)
+    {
+      if (myrank == 0 && ErrorMonitor->outfile)
+        ErrorMonitor->outfile << "CUDA Z4C failed in RK4 substep#" << iter_count
+                              << " at t = " << PhysTime
+                              << ", lev = " << lev << endl;
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+
+    Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
+
+    if (BH_num > 0 && lev == GH->levels - 1)
+    {
+      if (!z4c_cuda_compute_porg_rhs_resident(GH, lev, myrank, BH_num,
+                                              Porg, Porg1,
+                                              Sfx, Sfy, Sfz, Symmetry))
+      {
+        if (myrank == 0 && ErrorMonitor->outfile)
+          ErrorMonitor->outfile << "CUDA Z4C failed to interpolate black-hole shift at t = "
+                                << PhysTime << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
+      }
+      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]);
+        }
+      }
+    }
+
+    if (iter_count < 3)
+    {
+      Pp = GH->PatL[lev];
+      while (Pp)
+      {
+        MyList<Block> *BP = Pp->data->blb;
+        while (BP)
+        {
+          Block *cg = BP->data;
+          cg->swapList(SynchList_pre, SynchList_cor, myrank);
+          if (BP == Pp->data->ble)
+            break;
+          BP = BP->next;
+        }
+        Pp = Pp->next;
+      }
+      if (BH_num > 0 && lev == GH->levels - 1)
+      {
+        for (int ithBH = 0; ithBH < BH_num; ithBH++)
+        {
+          Porg[ithBH][0] = Porg1[ithBH][0];
+          Porg[ithBH][1] = Porg1[ithBH][1];
+          Porg[ithBH][2] = Porg1[ithBH][2];
+        }
+      }
+    }
+  }
+
+  z4c_cuda_download_level_state(GH->PatL[lev], SynchList_cor, myrank, true);
+
+#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 (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];
+    }
+  }
+#else
   double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
 #ifdef With_AHF
   AH_Step_Find(lev, dT_lev);
@@ -1042,9 +1592,13 @@ void Z4c_class::Step(int lev, int YN)
       Porg0[ithBH][2] = Porg1[ithBH][2];
     }
   }
+#endif
 }
 #else
 // for constraint preserving boundary (CPBC)
+#if USE_CUDA_Z4C && (ABEtype == 2)
+#error "USE_CUDA_Z4C resident path does not support CPBC"
+#endif
 #ifndef WithShell
 #error "CPBC only supports Shell"
 #endif

AMSS_NCKU_source/Z4c_rhs.f90

@@ -94,29 +94,31 @@
                Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon,                    &            
                Symmetry,Lev,eps,co)
 
+  if (co == 0) then
 #if (ABV == 0)  
-  call ricci_gamma(ex, X, Y, Z,                                      &
+    call ricci_gamma(ex, X, Y, Z,                                      &
-               chi,                                                  &
+                 chi,                                                  &
-               dxx    ,   gxy    ,   gxz    ,   dyy    ,   gyz    ,   dzz,&
+                 dxx    ,   gxy    ,   gxz    ,   dyy    ,   gyz    ,   dzz,&
-               Gamx   ,  Gamy    ,  Gamz    , &
+                 Gamx   ,  Gamy    ,  Gamz    , &
-               Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
+                 Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
-               Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
+                 Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
-               Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
+                 Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
-               Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
+                 Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
-               Symmetry)
+                 Symmetry)
 #endif
-  call constraint_bssn(ex, X, Y, Z,&
+    call constraint_bssn(ex, X, Y, Z,&
-               chi,trK, &
+                 chi,trK, &
-               dxx,gxy,gxz,dyy,gyz,dzz, &
+                 dxx,gxy,gxz,dyy,gyz,dzz, &
-               Axx,Axy,Axz,Ayy,Ayz,Azz, &
+                 Axx,Axy,Axz,Ayy,Ayz,Azz, &
-               Gamx,Gamy,Gamz,&
+                 Gamx,Gamy,Gamz,&
-               Lap,betax,betay,betaz,rho,Sx,Sy,Sz,&
+                 Lap,betax,betay,betaz,rho,Sx,Sy,Sz,&
-               Gamxxx, Gamxxy, Gamxxz,Gamxyy, Gamxyz, Gamxzz, &
+                 Gamxxx, Gamxxy, Gamxxz,Gamxyy, Gamxyz, Gamxzz, &
-               Gamyxx, Gamyxy, Gamyxz,Gamyyy, Gamyyz, Gamyzz, &
+                 Gamyxx, Gamyxy, Gamyxz,Gamyyy, Gamyyz, Gamyzz, &
-               Gamzxx, Gamzxy, Gamzxz,Gamzyy, Gamzyz, Gamzzz, &
+                 Gamzxx, Gamzxy, Gamzxz,Gamzyy, Gamzyz, Gamzzz, &
-               Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, &
+                 Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, &
-               Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon, &
+                 Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon, &
-               Symmetry)
+                 Symmetry)
+  endif
 
   return
 
@@ -227,6 +229,7 @@
   call get_Z4cparameters(kappa1,kappa2,kappa3,FF,eta)
 
 !!! sanity check
+#ifdef DEBUG
   dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) &
       +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz)                   &
       +sum(Gamx)+sum(Gamy)+sum(Gamz)                                           &
@@ -261,6 +264,7 @@
      gont = 1
      return
   endif
+#endif
 
   PI = dacos(-ONE)
 
@@ -1263,30 +1267,32 @@
 
   endif
 
+  if (co == 0) then
 #if (ABV == 0)  
-  call ricci_gamma(ex, X, Y, Z,                                      &
+    call ricci_gamma(ex, X, Y, Z,                                      &
-               chi,                                                  &
+                 chi,                                                  &
-               dxx    ,   gxy    ,   gxz    ,   dyy    ,   gyz    ,   dzz,&
+                 dxx    ,   gxy    ,   gxz    ,   dyy    ,   gyz    ,   dzz,&
-               Gamx   ,  Gamy    ,  Gamz    , &
+                 Gamx   ,  Gamy    ,  Gamz    , &
-               Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
+                 Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
-               Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
+                 Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
-               Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
+                 Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
-               Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
+                 Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
-               Symmetry)
+                 Symmetry)
 #endif
 
-  call constraint_bssn(ex, X, Y, Z,&
+    call constraint_bssn(ex, X, Y, Z,&
-               chi,trK, &
+                 chi,trK, &
-               dxx,gxy,gxz,dyy,gyz,dzz, &
+                 dxx,gxy,gxz,dyy,gyz,dzz, &
-               Axx,Axy,Axz,Ayy,Ayz,Azz, &
+                 Axx,Axy,Axz,Ayy,Ayz,Azz, &
-               Gamx,Gamy,Gamz,&
+                 Gamx,Gamy,Gamz,&
-               Lap,betax,betay,betaz,rho,Sx,Sy,Sz,&
+                 Lap,betax,betay,betaz,rho,Sx,Sy,Sz,&
-               Gamxxx, Gamxxy, Gamxxz,Gamxyy, Gamxyz, Gamxzz, &
+                 Gamxxx, Gamxxy, Gamxxz,Gamxyy, Gamxyz, Gamxzz, &
-               Gamyxx, Gamyxy, Gamyxz,Gamyyy, Gamyyz, Gamyzz, &
+                 Gamyxx, Gamyxy, Gamyxz,Gamyyy, Gamyyz, Gamyzz, &
-               Gamzxx, Gamzxy, Gamzxz,Gamzyy, Gamzyz, Gamzzz, &
+                 Gamzxx, Gamzxy, Gamzxz,Gamzyy, Gamzyz, Gamzzz, &
-               Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, &
+                 Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, &
-               Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon, &
+                 Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon, &
-               Symmetry)
+                 Symmetry)
+  endif
 
   gont = 0
 

AMSS_NCKU_source/Z4c_rhs_ss.f90

@@ -122,6 +122,7 @@
   call get_Z4cparameters(kappa1,kappa2,kappa3,FF,eta)
 
 !!! sanity check
+#ifdef DEBUG
   dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) &
       +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz)                   &
       +sum(Gamx)+sum(Gamy)+sum(Gamz)                                           &
@@ -156,6 +157,7 @@
      gont = 1
      return
   endif
+#endif
 
   PI = dacos(-ONE)
 
@@ -1388,41 +1390,43 @@
   call kodis_sh(ex,crho,sigma,R,TZ,TZ_rhs,SSS,Symmetry,eps,sst)
   endif
 
+  if (co == 0) then
 #if (ABV == 1)  
-  call ricci_gamma_ss(ex,crho,sigma,R,X, Y, Z,                                 &
+    call ricci_gamma_ss(ex,crho,sigma,R,X, Y, Z,                               &
-               drhodx, drhody, drhodz,                                         &
+                 drhodx, drhody, drhodz,                                       &
-               dsigmadx,dsigmady,dsigmadz,                                     &
+                 dsigmadx,dsigmady,dsigmadz,                                   &
-               dRdx,dRdy,dRdz,                                                 &
+                 dRdx,dRdy,dRdz,                                               &
-               drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz,                &
+                 drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz,              &
-               dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz,    &
+                 dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz,  &
-               dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz,                            &
+                 dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz,                          &
-               chi,                                                  &
+                 chi,                                                          &
-               dxx    ,   gxy    ,   gxz    ,   dyy    ,   gyz    ,   dzz,&
+                 dxx    ,   gxy    ,   gxz    ,   dyy    ,   gyz    ,   dzz,&
-               Gamx   ,  Gamy    ,  Gamz    , &
+                 Gamx   ,  Gamy    ,  Gamz    , &
-               Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
+                 Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
-               Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
+                 Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
-               Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
+                 Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
-               Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
+                 Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
-               Symmetry,Lev,sst)
+                 Symmetry,Lev,sst)
-  call constraint_bssn_ss(ex,crho,sigma,R,X, Y, Z,  &
-               drhodx, drhody, drhodz,                                         &
-               dsigmadx,dsigmady,dsigmadz,                                     &
-               dRdx,dRdy,dRdz,                                                 &
-               drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz,                &
-               dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz,    &
-               dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz,                            &
-               chi,trK, &
-               dxx,gxy,gxz,dyy,gyz,dzz, &
-               Axx,Axy,Axz,Ayy,Ayz,Azz, &
-               Gamx,Gamy,Gamz,&
-               Lap,betax,betay,betaz,rho,Sx,Sy,Sz,&
-               Gamxxx, Gamxxy, Gamxxz,Gamxyy, Gamxyz, Gamxzz, &
-               Gamyxx, Gamyxy, Gamyxz,Gamyyy, Gamyyz, Gamyzz, &
-               Gamzxx, Gamzxy, Gamzxz,Gamzyy, Gamzyz, Gamzzz, &
-               Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, &
-               Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon, &
-               Symmetry,Lev,sst)
 #endif
+    call constraint_bssn_ss(ex,crho,sigma,R,X, Y, Z,  &
+                 drhodx, drhody, drhodz,                                       &
+                 dsigmadx,dsigmady,dsigmadz,                                   &
+                 dRdx,dRdy,dRdz,                                               &
+                 drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz,              &
+                 dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz,  &
+                 dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz,                          &
+                 chi,trK, &
+                 dxx,gxy,gxz,dyy,gyz,dzz, &
+                 Axx,Axy,Axz,Ayy,Ayz,Azz, &
+                 Gamx,Gamy,Gamz,&
+                 Lap,betax,betay,betaz,rho,Sx,Sy,Sz,&
+                 Gamxxx, Gamxxy, Gamxxz,Gamxyy, Gamxyz, Gamxzz, &
+                 Gamyxx, Gamyxy, Gamyxz,Gamyyy, Gamyyz, Gamyzz, &
+                 Gamzxx, Gamzxy, Gamzxz,Gamzyy, Gamzyz, Gamzzz, &
+                 Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, &
+                 Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon, &
+                 Symmetry,Lev,sst)
+  endif
 
   gont = 0
 

AMSS_NCKU_source/bssn_class.h

@@ -48,6 +48,7 @@ public:
        double StartTime, TotalTime;
        double AnasTime, DumpTime, d2DumpTime, CheckTime;
        double LastAnas, LastConsOut;
+       bool cuda_level0_constraint_cache_valid;
        int *ConstraintRefreshLevels;
        double Courant;
        double numepss, numepsb, numepsh;

AMSS_NCKU_source/bssn_rhs_c.C

@@ -1098,12 +1098,12 @@ 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 );
@@ -1116,12 +1116,12 @@ 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 );

AMSS_NCKU_source/bssn_rhs_cuda.h

@@ -0,0 +1,127 @@
+#ifndef BSSN_RHS_CUDA_H
+#define BSSN_RHS_CUDA_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+enum {
+    BSSN_CUDA_STATE_COUNT = 24,
+    BSSN_CUDA_MATTER_COUNT = 10
+};
+
+int f_compute_rhs_bssn(int *ex, double &T,
+                       double *X, double *Y, double *Z,
+                       double *chi, double *trK,
+                       double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
+                       double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
+                       double *Gamx, double *Gamy, double *Gamz,
+                       double *Lap, double *betax, double *betay, double *betaz,
+                       double *dtSfx, double *dtSfy, double *dtSfz,
+                       double *chi_rhs, double *trK_rhs,
+                       double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
+                       double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
+                       double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
+                       double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
+                       double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
+                       double *rho, double *Sx, double *Sy, double *Sz,
+                       double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
+                       double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
+                       double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
+                       double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
+                       double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
+                       double *ham_Res, double *movx_Res, double *movy_Res, double *movz_Res,
+                       double *Gmx_Res, double *Gmy_Res, double *Gmz_Res,
+                       int &Symmetry, int &Lev, double &eps, int &co);
+
+int bssn_cuda_rk4_substep(void *block_tag,
+                          int *ex, double *X, double *Y, double *Z,
+                          double **state_host_in,
+                          double **state_host_out,
+                          double **matter_host,
+                          const double *propspeed,
+                          const double *soa_flat,
+                          const double *bbox,
+                          double &dT,
+                          double &T,
+                          int &RK4,
+                          int &apply_bam_bc,
+                          int &Symmetry,
+                          int &Lev,
+                          double &eps,
+                          int &co,
+                          int &use_zero_matter,
+                          int &keep_resident_state,
+                          int &apply_enforce_ga,
+                          double &chitiny);
+
+int bssn_cuda_copy_state_region_to_host(void *block_tag,
+                                        int state_index,
+                                        double *host_state,
+                                        int *ex,
+                                        int i0, int j0, int k0,
+                                        int sx, int sy, int sz);
+
+int bssn_cuda_copy_state_region_from_host(void *block_tag,
+                                          int state_index,
+                                          double *host_state,
+                                          int *ex,
+                                          int i0, int j0, int k0,
+                                          int sx, int sy, int sz);
+
+int bssn_cuda_download_resident_state(void *block_tag,
+                                      int *ex,
+                                      double **state_host_out);
+
+int bssn_cuda_download_constraint_outputs(int *ex,
+                                          double **constraint_host_out);
+
+int bssn_cuda_pack_state_region_to_host_buffer(void *block_tag,
+                                               int state_index,
+                                               double *host_buffer,
+                                               int *ex,
+                                               int i0, int j0, int k0,
+                                               int sx, int sy, int sz);
+
+int bssn_cuda_unpack_state_region_from_host_buffer(void *block_tag,
+                                                   int state_index,
+                                                   double *host_buffer,
+                                                   int *ex,
+                                                   int i0, int j0, int k0,
+                                                   int sx, int sy, int sz);
+
+int bssn_cuda_pack_state_batch_to_host_buffer(void *block_tag,
+                                              int state_count,
+                                              double *host_buffer,
+                                              int *ex,
+                                              int i0, int j0, int k0,
+                                              int sx, int sy, int sz);
+
+int bssn_cuda_unpack_state_batch_from_host_buffer(void *block_tag,
+                                                  int state_count,
+                                                  double *host_buffer,
+                                                  int *ex,
+                                                  int i0, int j0, int k0,
+                                                  int sx, int sy, int sz);
+
+int bssn_cuda_download_state_subset(void *block_tag,
+                                    int *ex,
+                                    int subset_count,
+                                    const int *state_indices,
+                                    double **state_host_out);
+
+int bssn_cuda_upload_state_subset(void *block_tag,
+                                  int *ex,
+                                  int subset_count,
+                                  const int *state_indices,
+                                  double **state_host_in);
+
+int bssn_cuda_has_resident_state(void *block_tag);
+
+void bssn_cuda_release_step_ctx(void *block_tag);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif