skip redundant MPI ghost cell syncs for stages 0, 1 & 2

BSSN 每个 RK4 时间步执行 4 次 MPI ghost zone 同步：
Stage 0（预测）结束后：Parallel::Sync(SynchList_pre)
Stage 1（校正 1）结束后：Parallel::Sync(SynchList_cor)
Stage 2（校正 2）结束后：Parallel::Sync(SynchList_cor)
Stage 3（校正 3）结束后：Parallel::Sync(SynchList_cor) ← 必要（为下一步提供 ghost）

bssnEM_class.C、Z4c_class.C 结构相同，一起修改了

AMSS_NCKU_Program.py

@@ -8,6 +8,14 @@
 ##
 ##################################################################
 
+## Guard against re-execution by multiprocessing child processes.
+## Without this, using 'spawn' or 'forkserver' context would cause every
+## worker to re-run the entire script, spawning exponentially more
+## workers (fork bomb).
+if __name__ != '__main__':
+    import sys as _sys
+    _sys.exit(0)
+
 
 ##################################################################
 
@@ -424,26 +432,31 @@ print(
 
 import plot_xiaoqu
 import plot_GW_strain_amplitude_xiaoqu
+from parallel_plot_helper import run_plot_tasks_parallel
+
+plot_tasks = []
 
 ## Plot black hole trajectory
-plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directory, figure_directory )
-plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot,   (binary_results_directory, figure_directory) ) )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot3D, (binary_results_directory, figure_directory) ) )
 
 ## Plot black hole separation vs. time
-plot_xiaoqu.generate_puncture_distence_plot( binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_distence_plot, (binary_results_directory, figure_directory) ) )
 
 ## Plot gravitational waveforms (psi4 and strain amplitude)
 for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_gravitational_wave_psi4_plot( binary_results_directory, figure_directory, i )
-    plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_gravitational_wave_psi4_plot, (binary_results_directory, figure_directory, i) ) )
+    plot_tasks.append( ( plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot, (binary_results_directory, figure_directory, i) ) )
 
 ## Plot ADM mass evolution
 for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_ADMmass_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_ADMmass_plot, (binary_results_directory, figure_directory, i) ) )
 
 ## Plot Hamiltonian constraint violation over time
 for i in range(input_data.grid_level):
-    plot_xiaoqu.generate_constraint_check_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_constraint_check_plot, (binary_results_directory, figure_directory, i) ) )
+
+run_plot_tasks_parallel(plot_tasks)
 
 ## Plot stored binary data
 plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )

AMSS_NCKU_source/Parallel.C

@@ -5286,6 +5286,203 @@ void Parallel::OutBdLow2Himix(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
   delete[] transfer_src;
   delete[] transfer_dst;
 }
+
+// Restrict_cached: cache grid segment lists, reuse buffers via transfer_cached
+void Parallel::Restrict_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
+                               MyList<var> *VarList1, MyList<var> *VarList2,
+                               int Symmetry, SyncCache &cache)
+{
+  if (!cache.valid)
+  {
+    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];
+    }
+
+    MyList<Parallel::gridseg> *dst = build_complete_gsl(PatcL);
+    for (int node = 0; node < cpusize; node++)
+    {
+      MyList<Parallel::gridseg> *src_owned = build_owned_gsl(PatfL, node, 2, Symmetry);
+      build_gstl(src_owned, dst, &cache.combined_src[node], &cache.combined_dst[node]);
+      if (src_owned) src_owned->destroyList();
+    }
+    if (dst) dst->destroyList();
+
+    cache.valid = true;
+  }
+
+  transfer_cached(cache.combined_src, cache.combined_dst, VarList1, VarList2, Symmetry, cache);
+}
+
+// OutBdLow2Hi_cached: cache grid segment lists, reuse buffers via transfer_cached
+void Parallel::OutBdLow2Hi_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
+                                  MyList<var> *VarList1, MyList<var> *VarList2,
+                                  int Symmetry, SyncCache &cache)
+{
+  if (!cache.valid)
+  {
+    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];
+    }
+
+    MyList<Parallel::gridseg> *dst = build_buffer_gsl(PatfL);
+    for (int node = 0; node < cpusize; node++)
+    {
+      MyList<Parallel::gridseg> *src_owned = build_owned_gsl(PatcL, node, 4, Symmetry);
+      build_gstl(src_owned, dst, &cache.combined_src[node], &cache.combined_dst[node]);
+      if (src_owned) src_owned->destroyList();
+    }
+    if (dst) dst->destroyList();
+
+    cache.valid = true;
+  }
+
+  transfer_cached(cache.combined_src, cache.combined_dst, VarList1, VarList2, Symmetry, cache);
+}
+
+// OutBdLow2Himix_cached: same as OutBdLow2Hi_cached but uses transfermix for unpacking
+void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
+                                     MyList<var> *VarList1, MyList<var> *VarList2,
+                                     int Symmetry, SyncCache &cache)
+{
+  if (!cache.valid)
+  {
+    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];
+    }
+
+    MyList<Parallel::gridseg> *dst = build_buffer_gsl(PatfL);
+    for (int node = 0; node < cpusize; node++)
+    {
+      MyList<Parallel::gridseg> *src_owned = build_owned_gsl(PatcL, node, 4, Symmetry);
+      build_gstl(src_owned, dst, &cache.combined_src[node], &cache.combined_dst[node]);
+      if (src_owned) src_owned->destroyList();
+    }
+    if (dst) dst->destroyList();
+
+    cache.valid = true;
+  }
+
+  // Use transfermix instead of transfer for mix-mode interpolation
+  int myrank;
+  MPI_Comm_size(MPI_COMM_WORLD, &cache.cpusize);
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  int cpusize = cache.cpusize;
+
+  int req_no = 0;
+  for (int node = 0; node < cpusize; node++)
+  {
+    if (node == myrank)
+    {
+      int length = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+      cache.recv_lengths[node] = length;
+      if (length > 0)
+      {
+        if (length > cache.recv_buf_caps[node])
+        {
+          if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
+          cache.recv_bufs[node] = new double[length];
+          cache.recv_buf_caps[node] = length;
+        }
+        data_packermix(cache.recv_bufs[node], cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+      }
+    }
+    else
+    {
+      int slength = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+      cache.send_lengths[node] = slength;
+      if (slength > 0)
+      {
+        if (slength > cache.send_buf_caps[node])
+        {
+          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++);
+      }
+      int rlength = data_packermix(0, cache.combined_src[node], cache.combined_dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+      cache.recv_lengths[node] = rlength;
+      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;
+        }
+        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no++);
+      }
+    }
+  }
+
+  MPI_Waitall(req_no, cache.reqs, cache.stats);
+
+  for (int node = 0; node < cpusize; node++)
+    if (cache.recv_bufs[node] && cache.recv_lengths[node] > 0)
+      data_packermix(cache.recv_bufs[node], cache.combined_src[node], cache.combined_dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+}
+
 // collect all buffer grid segments or blocks for given patch
 MyList<Parallel::gridseg> *Parallel::build_buffer_gsl(Patch *Pat)
 {

AMSS_NCKU_source/Parallel.h

@@ -130,6 +130,15 @@ namespace Parallel
   void OutBdLow2Himix(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
                       MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
                       int Symmetry);
+  void Restrict_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
+                       MyList<var> *VarList1, MyList<var> *VarList2,
+                       int Symmetry, SyncCache &cache);
+  void OutBdLow2Hi_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
+                          MyList<var> *VarList1, MyList<var> *VarList2,
+                          int Symmetry, SyncCache &cache);
+  void OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
+                             MyList<var> *VarList1, MyList<var> *VarList2,
+                             int Symmetry, SyncCache &cache);
   void Prolong(Patch *Patc, Patch *Patf,
                MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
                int Symmetry);

AMSS_NCKU_source/Z4c_class.C

@@ -485,25 +485,7 @@ void Z4c_class::Step(int lev, int YN)
   }
 #endif
 
-  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
-
-#ifdef WithShell
-  if (lev == 0)
-  {
-    clock_t prev_clock, curr_clock;
-    if (myrank == 0)
-      curr_clock = clock();
-    SH->Synch(SynchList_pre, Symmetry);
-    if (myrank == 0)
-    {
-      prev_clock = curr_clock;
-      curr_clock = clock();
-      cout << " Shell stuff synchronization used " 
-           << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) 
-           << " seconds! " << endl;
-    }
-  }
-#endif
+  // CA-RK4: skip post-prediction sync (redundant; ghost cells computable locally)
 
   // for black hole position
   if (BH_num > 0 && lev == GH->levels - 1)
@@ -868,6 +850,8 @@ void Z4c_class::Step(int lev, int YN)
     }
 #endif
 
+    // CA-RK4: only sync after last corrector (iter_count == 3); stages 1 & 2 are redundant
+    if (iter_count == 3) {
       Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
 
 #ifdef WithShell
@@ -887,6 +871,7 @@ void Z4c_class::Step(int lev, int YN)
         }
       }
 #endif
+    } // end CA-RK4 guard
     // for black hole position
     if (BH_num > 0 && lev == GH->levels - 1)
     {
@@ -1558,7 +1543,7 @@ void Z4c_class::Step(int lev, int YN)
     }
   }
 
-  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
+  // CA-RK4: skip post-prediction MPI ghost sync (redundant; ghost cells computable locally)
 
   if (lev == 0)
   {
@@ -2120,6 +2105,8 @@ void Z4c_class::Step(int lev, int YN)
       }
     }
 
+    // CA-RK4: only MPI sync after last corrector (iter_count == 3); stages 1 & 2 are redundant
+    if (iter_count == 3)
       Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
 
     if (lev == 0)

AMSS_NCKU_source/bssnEM_class.C

@@ -1221,25 +1221,7 @@ void bssnEM_class::Step(int lev, int YN)
   }
 #endif
 
-  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
-
-#ifdef WithShell
-  if (lev == 0)
-  {
-    clock_t prev_clock, curr_clock;
-    if (myrank == 0)
-      curr_clock = clock();
-    SH->Synch(SynchList_pre, Symmetry);
-    if (myrank == 0)
-    {
-      prev_clock = curr_clock;
-      curr_clock = clock();
-      cout << " Shell stuff synchronization used " 
-           << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) 
-           << " seconds! " << endl;
-    }
-  }
-#endif
+  // CA-RK4: skip post-prediction sync (redundant; ghost cells computable locally)
 
   // for black hole position
   if (BH_num > 0 && lev == GH->levels - 1)
@@ -1683,6 +1665,8 @@ void bssnEM_class::Step(int lev, int YN)
     }
 #endif
 
+    // CA-RK4: only sync after last corrector (iter_count == 3); stages 1 & 2 are redundant
+    if (iter_count == 3) {
       Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
 
 #ifdef WithShell
@@ -1702,6 +1686,7 @@ void bssnEM_class::Step(int lev, int YN)
         }
       }
 #endif
+    } // end CA-RK4 guard
     // for black hole position
     if (BH_num > 0 && lev == GH->levels - 1)
     {

AMSS_NCKU_source/bssn_class.C

@@ -2426,9 +2426,9 @@ void bssn_class::RecursiveStep(int lev)
 #endif
 
 #if (REGLEV == 0)
-  GH->Regrid_Onelevel(lev, Symmetry, BH_num, Porgbr, Porg0,
+  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);
+                      fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
   for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); }
 #endif
 }
@@ -2605,9 +2605,9 @@ void bssn_class::ParallelStep()
   delete[] tporg;
   delete[] tporgo;
 #if (REGLEV == 0)
-  GH->Regrid_Onelevel(GH->mylev, Symmetry, BH_num, Porgbr, Porg0,
+  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);
+                      fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
   for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); }
 #endif
 }
@@ -2772,9 +2772,9 @@ void bssn_class::ParallelStep()
       if (lev + 1 >= GH->movls)
       {
         //	       GH->Regrid_Onelevel_aux(lev,Symmetry,BH_num,Porgbr,Porg0,
-        GH->Regrid_Onelevel(lev + 1, Symmetry, BH_num, Porgbr, Porg0,
+        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);
+                            fgt(PhysTime - dT_levp1, StartTime, dT_levp1 / 2), ErrorMonitor))
         for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); }
 
         //               a_stream.clear();
@@ -2787,9 +2787,9 @@ void bssn_class::ParallelStep()
     // for this level
     if (YN == 1)
     {
-      GH->Regrid_Onelevel(lev, Symmetry, BH_num, Porgbr, Porg0,
+      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);
+                          fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
       for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); }
 
       //               a_stream.clear();
@@ -2806,9 +2806,9 @@ void bssn_class::ParallelStep()
         if (YN == 1)
         {
           //	   GH->Regrid_Onelevel_aux(lev-2,Symmetry,BH_num,Porgbr,Porg0,
-          GH->Regrid_Onelevel(lev - 1, Symmetry, BH_num, Porgbr, Porg0,
+          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);
+                              fgt(PhysTime - dT_lev, StartTime, dT_levm1 / 2), ErrorMonitor))
           for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); }
 
           //               a_stream.clear();
@@ -2822,9 +2822,9 @@ void bssn_class::ParallelStep()
         if (i % 4 == 3)
         {
           //	   GH->Regrid_Onelevel_aux(lev-2,Symmetry,BH_num,Porgbr,Porg0,
-          GH->Regrid_Onelevel(lev - 1, Symmetry, BH_num, Porgbr, Porg0,
+          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);
+                              fgt(PhysTime - dT_lev, StartTime, dT_levm1 / 2), ErrorMonitor))
           for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); }
 
           //               a_stream.clear();
@@ -3349,27 +3349,7 @@ void bssn_class::Step(int lev, int YN)
   }
 #endif
 
-  Parallel::AsyncSyncState async_pre;
-  Parallel::Sync_start(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev], async_pre);
-
-#ifdef WithShell
-  if (lev == 0)
-  {
-    clock_t prev_clock, curr_clock;
-    if (myrank == 0)
-      curr_clock = clock();
-    SH->Synch(SynchList_pre, Symmetry);
-    if (myrank == 0)
-    {
-      prev_clock = curr_clock;
-      curr_clock = clock();
-      cout << " Shell stuff synchronization used "
-           << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC)
-           << " seconds! " << endl;
-    }
-  }
-#endif
-  Parallel::Sync_finish(sync_cache_pre[lev], async_pre, SynchList_pre, Symmetry);
+  // CA-RK4: skip post-prediction sync (redundant; ghost cells computable locally)
 
 #ifdef WithShell
   // Complete non-blocking error reduction and check
@@ -3709,6 +3689,8 @@ void bssn_class::Step(int lev, int YN)
     }
 #endif
 
+    // CA-RK4: only sync after last corrector (iter_count == 3); stages 1 & 2 are redundant
+    if (iter_count == 3) {
       Parallel::AsyncSyncState async_cor;
       Parallel::Sync_start(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev], async_cor);
 
@@ -3730,6 +3712,7 @@ void bssn_class::Step(int lev, int YN)
       }
 #endif
       Parallel::Sync_finish(sync_cache_cor[lev], async_cor, SynchList_cor, Symmetry);
+    } // end CA-RK4 guard
 
 #ifdef WithShell
     // Complete non-blocking error reduction and check
@@ -5819,21 +5802,11 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 #endif
 
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],SynchList_pre,SL,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, GH->bdsul[lev], Symmetry);
@@ -5880,21 +5853,11 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 #endif
 
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SL, SL, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SL, SL, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],SL,SL,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->bdsul[lev], Symmetry);
@@ -5969,21 +5932,11 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],SynchList_pre,SL,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, GH->bdsul[lev], Symmetry);
@@ -6001,21 +5954,11 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
       Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SL, SL, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SL, SL, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],SL,SL,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->bdsul[lev], Symmetry);
@@ -6076,21 +6019,11 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SynchList_pre, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SynchList_pre, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],SynchList_pre,SynchList_cor,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, GH->bdsul[lev], Symmetry);
@@ -6110,21 +6043,11 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, StateList, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, StateList, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],StateList,SynchList_cor,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, GH->bdsul[lev], Symmetry);
@@ -6161,21 +6084,11 @@ void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
       }
 
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SynchList_pre, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SynchList_pre, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],SynchList_pre,SynchList_cor,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, GH->bdsul[lev], Symmetry);
@@ -6184,21 +6097,11 @@ void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
     else // no time refinement levels and for all same time levels
     {
 #if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
 #if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, StateList, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Hi(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry);
 #elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, StateList, SynchList_cor, Symmetry);
+      Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry);
 #endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
 #elif (RPB == 1)
       //       Parallel::OutBdLow2Hi_bam(GH->PatL[lev-1],GH->PatL[lev],StateList,SynchList_cor,Symmetry);
       Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, GH->bdsul[lev], Symmetry);

AMSS_NCKU_source/bssn_rhs.f90

@@ -945,103 +945,60 @@
   SSA(2)=SYM
   SSA(3)=ANTI
 
-!!!!!!!!!advection term part
+!!!!!!!!!advection term + Kreiss-Oliger dissipation (merged for cache efficiency)
+! lopsided_kodis shares the symmetry_bd buffer between advection and
+! dissipation, eliminating redundant full-grid copies. For metric variables
+! gxx/gyy/gzz (=dxx/dyy/dzz+1): kodis stencil coefficients sum to zero,
+! so the constant offset has no effect on dissipation.
 
-  call lopsided(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS)
-  call lopsided(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS)
-  call lopsided(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA)
-  call lopsided(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS)
-  call lopsided(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA)
-  call lopsided(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_kodis(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,eps)
+  call lopsided_kodis(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps)
+  call lopsided_kodis(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
 
-  call lopsided(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS)
-  call lopsided(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS)
-  call lopsided(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA)
-  call lopsided(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS)
-  call lopsided(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA)
-  call lopsided(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,eps)
+  call lopsided_kodis(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA,eps)
+  call lopsided_kodis(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA,eps)
+  call lopsided_kodis(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
 
-  call lopsided(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS)
-  call lopsided(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_kodis(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS,eps)
 
-  call lopsided(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS)
-  call lopsided(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS)
-  call lopsided(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA)
-!!
+  call lopsided_kodis(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS,eps)
+  call lopsided_kodis(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS,eps)
+  call lopsided_kodis(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA,eps)
+
+#if 1 
+!! bam does not apply dissipation on gauge variables
+  call lopsided_kodis(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+#if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7)
+  call lopsided_kodis(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS,eps)
+  call lopsided_kodis(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS,eps)
+  call lopsided_kodis(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,eps)
+#endif
+#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
+  call lopsided_kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,eps)
+  call lopsided_kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,eps)
+  call lopsided_kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA,eps)
+#endif
+#else
+! No dissipation on gauge variables (advection only)
   call lopsided(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS)
-
 #if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7)
   call lopsided(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS)
   call lopsided(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS)
   call lopsided(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA)
 #endif
-
 #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
   call lopsided(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS)
   call lopsided(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS)
   call lopsided(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA)
 #endif
-
-  if(eps>0)then 
-! usual Kreiss-Oliger dissipation      
-  call kodis(ex,X,Y,Z,chi,chi_rhs,SSS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,trK,trK_rhs,SSS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,dxx,gxx_rhs,SSS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,gxy,gxy_rhs,AAS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,gxz,gxz_rhs,ASA,Symmetry,eps)
-  call kodis(ex,X,Y,Z,dyy,gyy_rhs,SSS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,gyz,gyz_rhs,SAA,Symmetry,eps)
-  call kodis(ex,X,Y,Z,dzz,gzz_rhs,SSS,Symmetry,eps)
-#if 0
-#define i 42
-#define j 40
-#define k 40
-if(Lev == 1)then
-write(*,*) X(i),Y(j),Z(k)
-write(*,*) "before",Axx_rhs(i,j,k)
-endif
-#undef i
-#undef j
-#undef k
-!!stop
 #endif
-  call kodis(ex,X,Y,Z,Axx,Axx_rhs,SSS,Symmetry,eps)
-#if 0
-#define i 42
-#define j 40
-#define k 40
-if(Lev == 1)then
-write(*,*) X(i),Y(j),Z(k)
-write(*,*) "after",Axx_rhs(i,j,k)
-endif
-#undef i
-#undef j
-#undef k
-!!stop
-#endif
-  call kodis(ex,X,Y,Z,Axy,Axy_rhs,AAS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,Axz,Axz_rhs,ASA,Symmetry,eps)
-  call kodis(ex,X,Y,Z,Ayy,Ayy_rhs,SSS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,Ayz,Ayz_rhs,SAA,Symmetry,eps)
-  call kodis(ex,X,Y,Z,Azz,Azz_rhs,SSS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,Gamx,Gamx_rhs,ASS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,Gamy,Gamy_rhs,SAS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,Gamz,Gamz_rhs,SSA,Symmetry,eps)
-
-#if 1 
-!! bam does not apply dissipation on gauge variables
-  call kodis(ex,X,Y,Z,Lap,Lap_rhs,SSS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,betax,betax_rhs,ASS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,betay,betay_rhs,SAS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,betaz,betaz_rhs,SSA,Symmetry,eps)
-#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
-  call kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,ASS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,SAS,Symmetry,eps)
-  call kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,SSA,Symmetry,eps)
-#endif
-#endif
-
-  endif
 
   if(co == 0)then
 ! ham_Res = trR + 2/3 * K^2 - A_ij * A^ij - 16 * PI * rho

AMSS_NCKU_source/cgh.C

@@ -1301,13 +1301,13 @@ bool cgh::Interp_One_Point(MyList<var> *VarList,
 }
 
 
-void cgh::Regrid_Onelevel(int lev, int Symmetry, int BH_num, double **Porgbr, double **Porg0,
+bool cgh::Regrid_Onelevel(int lev, int Symmetry, int BH_num, double **Porgbr, double **Porg0,
                           MyList<var> *OldList, MyList<var> *StateList,
                           MyList<var> *FutureList, MyList<var> *tmList, bool BB,
                           monitor *ErrorMonitor)
 {
   if (lev < movls)
-    return;
+    return false;
 
 #if (0)
   // #if (PSTR == 1 || PSTR == 2)
@@ -1396,7 +1396,7 @@ void cgh::Regrid_Onelevel(int lev, int Symmetry, int BH_num, double **Porgbr, do
       for (bhi = 0; bhi < BH_num; bhi++)
         delete[] tmpPorg[bhi];
       delete[] tmpPorg;
-      return;
+      return false;
     }
     // x direction
     rr = (Porg0[bhi][0] - handle[lev][grd][0]) / dX;
@@ -1500,6 +1500,7 @@ void cgh::Regrid_Onelevel(int lev, int Symmetry, int BH_num, double **Porgbr, do
   for (int bhi = 0; bhi < BH_num; bhi++)
     delete[] tmpPorg[bhi];
   delete[] tmpPorg;
+  return tot_flag;
 }
 
 

AMSS_NCKU_source/cgh.h

@@ -74,7 +74,7 @@ public:
                                MyList<var> *OldList, MyList<var> *StateList,
                                MyList<var> *FutureList, MyList<var> *tmList,
                                int Symmetry, bool BB);
-   void Regrid_Onelevel(int lev, int Symmetry, int BH_num, double **Porgbr, double **Porg0,
+   bool Regrid_Onelevel(int lev, int Symmetry, int BH_num, double **Porgbr, double **Porg0,
                         MyList<var> *OldList, MyList<var> *StateList,
                         MyList<var> *FutureList, MyList<var> *tmList, bool BB,
                         monitor *ErrorMonitor);

AMSS_NCKU_source/diff_new.f90

@@ -69,6 +69,8 @@
   fy = ZEO
   fz = ZEO
 
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
+!DIR$ UNROLL PARTIAL(4)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1
@@ -371,6 +373,8 @@
   fxz = ZEO
   fyz = ZEO
 
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
+!DIR$ UNROLL PARTIAL(4)
   do k=1,ex(3)-1
   do j=1,ex(2)-1
   do i=1,ex(1)-1

AMSS_NCKU_source/fdderivs_c.C

@@ -0,0 +1,268 @@
+#include "tool.h"
+void fdderivs(const int ex[3],
+              const double *f,
+              double *fxx, double *fxy, double *fxz,
+              double *fyy, double *fyz, double *fzz,
+              const double *X, const double *Y, const double *Z,
+              double SYM1, double SYM2, double SYM3,
+              int Symmetry, int onoff)
+{
+    (void)onoff;
+
+    const int NO_SYMM = 0, EQ_SYMM = 1;
+    const double ZEO = 0.0, ONE = 1.0, TWO = 2.0;
+    const double F1o4   = 2.5e-1;          // 1/4
+    const double F8     = 8.0;
+    const double F16    = 16.0;
+    const double F30    = 30.0;
+    const double F1o12  = ONE / 12.0;
+    const double F1o144 = ONE / 144.0;
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    int iminF = 1, jminF = 1, kminF = 1;
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -1;
+    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -1;
+    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -1;
+
+    const double SoA[3] = { SYM1, SYM2, SYM3 };
+
+    /* fh: (ex1+2)*(ex2+2)*(ex3+2) because ord=2 */
+    const size_t nx = (size_t)ex1 + 2;
+    const size_t ny = (size_t)ex2 + 2;
+    const size_t nz = (size_t)ex3 + 2;
+    const size_t fh_size = nx * ny * nz;
+
+    static double *fh = NULL;
+    static size_t cap = 0;
+
+    if (fh_size > cap) {
+        free(fh);
+        fh = (double*)aligned_alloc(64, fh_size * sizeof(double));
+        cap = fh_size;
+    }
+    // double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return;
+
+    symmetry_bd(2, ex, f, fh, SoA);
+
+    /* 系数：按 Fortran 原式 */
+    const double Sdxdx = ONE / (dX * dX);
+    const double Sdydy = ONE / (dY * dY);
+    const double Sdzdz = ONE / (dZ * dZ);
+
+    const double Fdxdx = F1o12 / (dX * dX);
+    const double Fdydy = F1o12 / (dY * dY);
+    const double Fdzdz = F1o12 / (dZ * dZ);
+
+    const double Sdxdy = F1o4 / (dX * dY);
+    const double Sdxdz = F1o4 / (dX * dZ);
+    const double Sdydz = F1o4 / (dY * dZ);
+
+    const double Fdxdy = F1o144 / (dX * dY);
+    const double Fdxdz = F1o144 / (dX * dZ);
+    const double Fdydz = F1o144 / (dY * dZ);
+
+    /* 输出清零：fxx,fyy,fzz,fxy,fxz,fyz = 0 */
+    const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3;
+    for (size_t p = 0; p < all; ++p) {
+        fxx[p] = ZEO; fyy[p] = ZEO; fzz[p] = ZEO;
+        fxy[p] = ZEO; fxz[p] = ZEO; fyz[p] = ZEO;
+    }
+
+    /*
+     * Fortran:
+     * do k=1,ex3-1
+     * do j=1,ex2-1
+     * do i=1,ex1-1
+     */
+    
+    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
+                const int iF = i0 + 1;
+                const size_t p = idx_ex(i0, j0, k0, ex);
+
+                /* 高阶分支：i±2,j±2,k±2 都在范围内 */
+                if ((iF + 2) <= imaxF && (iF - 2) >= iminF &&
+                    (jF + 2) <= jmaxF && (jF - 2) >= jminF &&
+                    (kF + 2) <= kmaxF && (kF - 2) >= kminF)
+                {
+                    fxx[p] = Fdxdx * (
+                        -fh[idx_fh_F_ord2(iF - 2, jF,     kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] -
+                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF + 2, jF,     kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]
+                    );
+
+                    fyy[p] = Fdydy * (
+                        -fh[idx_fh_F_ord2(iF,     jF - 2, kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] -
+                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF,     jF + 2, kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]
+                    );
+
+                    fzz[p] = Fdzdz * (
+                        -fh[idx_fh_F_ord2(iF,     jF,     kF - 2, ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] -
+                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF,     jF,     kF + 2, ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]
+                    );
+
+                    /* fxy 高阶：完全照搬 Fortran 的括号结构 */
+                    {
+                        const double t_jm2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF - 2, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF - 2, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF - 2, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF - 2, kF, ex)] );
+
+                        const double t_jm1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF - 1, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF - 1, kF, ex)] );
+
+                        const double t_jp1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF + 1, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF + 1, kF, ex)] );
+
+                        const double t_jp2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF + 2, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF + 2, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF + 2, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF + 2, kF, ex)] );
+
+                        fxy[p] = Fdxdy * ( t_jm2 - F8 * t_jm1 + F8 * t_jp1 - t_jp2 );
+                    }
+
+                    /* fxz 高阶 */
+                    {
+                        const double t_km2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 2, ex)] );
+
+                        const double t_km1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 1, ex)] );
+
+                        const double t_kp1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 1, ex)] );
+
+                        const double t_kp2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 2, ex)] );
+
+                        fxz[p] = Fdxdz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
+                    }
+
+                    /* fyz 高阶 */
+                    {
+                        const double t_km2 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 2, ex)] );
+
+                        const double t_km1 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 1, ex)] );
+
+                        const double t_kp1 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 1, ex)] );
+
+                        const double t_kp2 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 2, ex)] );
+
+                        fyz[p] = Fdydz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
+                    }
+                }
+                /* 二阶分支：i±1,j±1,k±1 在范围内 */
+                else if ((iF + 1) <= imaxF && (iF - 1) >= iminF &&
+                         (jF + 1) <= jmaxF && (jF - 1) >= jminF &&
+                         (kF + 1) <= kmaxF && (kF - 1) >= kminF)
+                {
+                    fxx[p] = Sdxdx * (
+                        fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] -
+                        TWO * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] +
+                        fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]
+                    );
+
+                    fyy[p] = Sdydy * (
+                        fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] -
+                        TWO * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] +
+                        fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]
+                    );
+
+                    fzz[p] = Sdzdz * (
+                        fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] -
+                        TWO * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] +
+                        fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]
+                    );
+
+                    fxy[p] = Sdxdy * (
+                        fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)] -
+                        fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)] -
+                        fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)] +
+                        fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)]
+                    );
+
+                    fxz[p] = Sdxdz * (
+                        fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)] -
+                        fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)] -
+                        fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)] +
+                        fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)]
+                    );
+
+                    fyz[p] = Sdydz * (
+                        fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)] -
+                        fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)] -
+                        fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)] +
+                        fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)]
+                    );
+                }else{
+                    fxx[p] = 0.0;
+                    fyy[p] = 0.0;
+                    fzz[p] = 0.0;
+                    fxy[p] = 0.0;
+                    fxz[p] = 0.0;
+                    fyz[p] = 0.0;
+                }
+            }
+        }
+    }
+
+    // free(fh);
+}

AMSS_NCKU_source/fderivs_c.C

@@ -0,0 +1,150 @@
+#include "tool.h"
+
+/*
+ * C 版 fderivs
+ *
+ * Fortran:
+ * subroutine fderivs(ex,f,fx,fy,fz,X,Y,Z,SYM1,SYM2,SYM3,symmetry,onoff)
+ *
+ * 约定：
+ *   f, fx, fy, fz: ex1*ex2*ex3，按 idx_ex 布局
+ *   X: ex1, Y: ex2, Z: ex3
+ */
+void fderivs(const int ex[3],
+             const double *f,
+             double *fx, double *fy, double *fz,
+             const double *X, const double *Y, const double *Z,
+             double SYM1, double SYM2, double SYM3,
+             int Symmetry, int onoff)
+{
+    (void)onoff; // Fortran 里没用到
+
+    const double ZEO = 0.0, ONE = 1.0;
+    const double TWO = 2.0, EIT = 8.0;
+    const double F12 = 12.0;
+
+    const int NO_SYMM = 0, EQ_SYMM = 1; // OCTANT=2 在本子程序里不直接用
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    // dX = X(2)-X(1) -> C: X[1]-X[0]
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+
+    // Fortran 1-based bounds
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    int iminF = 1, jminF = 1, kminF = 1;
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -1;
+    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -1;
+    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -1;
+
+    // SoA(1:3) = SYM1,SYM2,SYM3
+    const double SoA[3] = { SYM1, SYM2, SYM3 };
+
+    // fh: (ex1+2)*(ex2+2)*(ex3+2) because ord=2
+    const size_t nx = (size_t)ex1 + 2;
+    const size_t ny = (size_t)ex2 + 2;
+    const size_t nz = (size_t)ex3 + 2;
+    const size_t fh_size = nx * ny * nz;
+    static double *fh = NULL;
+    static size_t cap = 0;
+
+    if (fh_size > cap) {
+        free(fh);
+        fh = (double*)aligned_alloc(64, fh_size * sizeof(double));
+        cap = fh_size;
+    }
+    // double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return;
+
+    // call symmetry_bd(2,ex,f,fh,SoA)
+    symmetry_bd(2, ex, f, fh, SoA);
+
+    const double d12dx = ONE / F12 / dX;
+    const double d12dy = ONE / F12 / dY;
+    const double d12dz = ONE / F12 / dZ;
+
+    const double d2dx  = ONE / TWO / dX;
+    const double d2dy  = ONE / TWO / dY;
+    const double d2dz  = ONE / TWO / dZ;
+
+    // fx = fy = fz = 0
+    const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3;
+    for (size_t p = 0; p < all; ++p) {
+        fx[p] = ZEO;
+        fy[p] = ZEO;
+        fz[p] = ZEO;
+    }
+
+    /*
+     * Fortran loops:
+     * do k=1,ex3-1
+     * do j=1,ex2-1
+     * do i=1,ex1-1
+     *
+     * C: k0=0..ex3-2, j0=0..ex2-2, i0=0..ex1-2
+     */
+    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
+                const int iF = i0 + 1;
+                const size_t p = idx_ex(i0, j0, k0, ex);
+
+                // if(i+2 <= imax .and. i-2 >= imin ... )  (全是 Fortran 索引)
+                if ((iF + 2) <= imaxF && (iF - 2) >= iminF &&
+                    (jF + 2) <= jmaxF && (jF - 2) >= jminF &&
+                    (kF + 2) <= kmaxF && (kF - 2) >= kminF)
+                {
+                    fx[p] = d12dx * (
+                        fh[idx_fh_F_ord2(iF - 2, jF,     kF,     ex)] -
+                        EIT * fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] +
+                        EIT * fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF + 2, jF,     kF,     ex)]
+                    );
+
+                    fy[p] = d12dy * (
+                        fh[idx_fh_F_ord2(iF,     jF - 2, kF,     ex)] -
+                        EIT * fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] +
+                        EIT * fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF,     jF + 2, kF,     ex)]
+                    );
+
+                    fz[p] = d12dz * (
+                        fh[idx_fh_F_ord2(iF,     jF,     kF - 2, ex)] -
+                        EIT * fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] +
+                        EIT * fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)] -
+                        fh[idx_fh_F_ord2(iF,     jF,     kF + 2, ex)]
+                    );
+                }
+                // elseif(i+1 <= imax .and. i-1 >= imin ...)
+                else if ((iF + 1) <= imaxF && (iF - 1) >= iminF &&
+                         (jF + 1) <= jmaxF && (jF - 1) >= jminF &&
+                         (kF + 1) <= kmaxF && (kF - 1) >= kminF)
+                {
+                    fx[p] = d2dx * (
+                        -fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] +
+                         fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]
+                    );
+
+                    fy[p] = d2dy * (
+                        -fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] +
+                         fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]
+                    );
+
+                    fz[p] = d2dz * (
+                        -fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] +
+                         fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]
+                    );
+                }
+            }
+        }
+    }
+
+    // free(fh);
+}

AMSS_NCKU_source/fmisc.f90

@@ -883,13 +883,17 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
 
   integer::i
 
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
   funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
    do i=0,ord-1
       funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
    enddo
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
    do i=0,ord-1
       funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
    enddo
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
    do i=0,ord-1
       funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
    enddo
@@ -1112,6 +1116,7 @@ end subroutine d2dump
 ! Lagrangian polynomial interpolation
 !------------------------------------------------------------------------------
 
+!DIR$ ATTRIBUTES FORCEINLINE :: polint
   subroutine polint(xa, ya, x, y, dy, ordn)
   implicit none
 

AMSS_NCKU_source/kodiss.f90

@@ -65,6 +65,8 @@ real*8,intent(in) :: eps
 !                       dx^4
 
 !  note the sign (-1)^r-1, now r=2
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
+!DIR$ UNROLL PARTIAL(4)
   do k=1,ex(3)
   do j=1,ex(2)
   do i=1,ex(1)

AMSS_NCKU_source/kodiss_c.C

@@ -0,0 +1,109 @@
+#include "tool.h"
+
+/*
+ * C 版 kodis
+ *
+ * Fortran signature:
+ * subroutine kodis(ex,X,Y,Z,f,f_rhs,SoA,Symmetry,eps)
+ *
+ * 约定：
+ *   X: ex1, Y: ex2, Z: ex3
+ *   f, f_rhs: ex1*ex2*ex3 按 idx_ex 布局
+ *   SoA[3]
+ *   eps: double
+ */
+void kodis(const int ex[3],
+           const double *X, const double *Y, const double *Z,
+           const double *f, double *f_rhs,
+           const double SoA[3],
+           int Symmetry, double eps)
+{
+    const double ONE = 1.0, SIX = 6.0, FIT = 15.0, TWT = 20.0;
+    const double cof = 64.0;             // 2^6
+    const int NO_SYMM = 0, OCTANT = 2;
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    // Fortran: dX = X(2)-X(1) -> C: X[1]-X[0]
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+    (void)ONE; // ONE 在原 Fortran 里只是参数，这里不一定用得上
+
+    // Fortran: imax=ex(1) 等是 1-based 上界
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    // Fortran: imin=jmin=kmin=1，某些对称情况变 -2
+    int iminF = 1, jminF = 1, kminF = 1;
+
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -2;
+    if (Symmetry == OCTANT && fabs(X[0]) < dX) iminF = -2;
+    if (Symmetry == OCTANT && fabs(Y[0]) < dY) jminF = -2;
+
+    // 分配 fh：大小 (ex1+3)*(ex2+3)*(ex3+3)，对应 ord=3
+    const size_t nx = (size_t)ex1 + 3;
+    const size_t ny = (size_t)ex2 + 3;
+    const size_t nz = (size_t)ex3 + 3;
+    const size_t fh_size = nx * ny * nz;
+
+    double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return;
+
+    // Fortran: call symmetry_bd(3,ex,f,fh,SoA)
+    symmetry_bd(3, ex, f, fh, SoA);
+
+    /*
+     * Fortran loops:
+     * do k=1,ex3
+     * do j=1,ex2
+     * do i=1,ex1
+     *
+     * C: k0=0..ex3-1, j0=0..ex2-1, i0=0..ex1-1
+     * 并定义 Fortran index: iF=i0+1, ...
+     */
+    for (int k0 = 0; k0 < ex3; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 < ex2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 < ex1; ++i0) {
+                const int iF = i0 + 1;
+
+                // Fortran if 条件：
+                // i-3 >= imin .and. i+3 <= imax  等（都是 Fortran 索引）
+                if ((iF - 3) >= iminF && (iF + 3) <= imaxF &&
+                    (jF - 3) >= jminF && (jF + 3) <= jmaxF &&
+                    (kF - 3) >= kminF && (kF + 3) <= kmaxF)
+                {
+                    const size_t p = idx_ex(i0, j0, k0, ex);
+
+                    // 三个方向各一份同型的 7 点组合（实际上是对称的 6th-order dissipation/filter 核）
+                    const double Dx_term =
+                        ( (fh[idx_fh_F(iF - 3, jF, kF, ex)] + fh[idx_fh_F(iF + 3, jF, kF, ex)]) -
+                          SIX * (fh[idx_fh_F(iF - 2, jF, kF, ex)] + fh[idx_fh_F(iF + 2, jF, kF, ex)]) +
+                          FIT * (fh[idx_fh_F(iF - 1, jF, kF, ex)] + fh[idx_fh_F(iF + 1, jF, kF, ex)]) -
+                          TWT *  fh[idx_fh_F(iF    , jF, kF, ex)] ) / dX;
+
+                    const double Dy_term =
+                        ( (fh[idx_fh_F(iF, jF - 3, kF, ex)] + fh[idx_fh_F(iF, jF + 3, kF, ex)]) -
+                          SIX * (fh[idx_fh_F(iF, jF - 2, kF, ex)] + fh[idx_fh_F(iF, jF + 2, kF, ex)]) +
+                          FIT * (fh[idx_fh_F(iF, jF - 1, kF, ex)] + fh[idx_fh_F(iF, jF + 1, kF, ex)]) -
+                          TWT *  fh[idx_fh_F(iF, jF    , kF, ex)] ) / dY;
+
+                    const double Dz_term =
+                        ( (fh[idx_fh_F(iF, jF, kF - 3, ex)] + fh[idx_fh_F(iF, jF, kF + 3, ex)]) -
+                          SIX * (fh[idx_fh_F(iF, jF, kF - 2, ex)] + fh[idx_fh_F(iF, jF, kF + 2, ex)]) +
+                          FIT * (fh[idx_fh_F(iF, jF, kF - 1, ex)] + fh[idx_fh_F(iF, jF, kF + 1, ex)]) -
+                          TWT *  fh[idx_fh_F(iF, jF, kF    , ex)] ) / dZ;
+
+                    // Fortran:
+                    // f_rhs(i,j,k) = f_rhs(i,j,k) + eps/cof*(Dx_term + Dy_term + Dz_term)
+                    f_rhs[p] += (eps / cof) * (Dx_term + Dy_term + Dz_term);
+                }
+            }
+        }
+    }
+
+    free(fh);
+}

AMSS_NCKU_source/lopsided_c.C

@@ -0,0 +1,255 @@
+#include "tool.h"
+/*
+ * 你需要提供 symmetry_bd 的 C 版本（或 Fortran 绑到 C 的接口）。
+ * Fortran: call symmetry_bd(3,ex,f,fh,SoA)
+ *
+ * 约定：
+ *   nghost = 3
+ *   ex[3]  = {ex1,ex2,ex3}
+ *   f      = 原始网格 (ex1*ex2*ex3)
+ *   fh     = 扩展网格 ((ex1+3)*(ex2+3)*(ex3+3))，对应 Fortran 的 (-2:ex1, ...)
+ *   SoA[3] = 输入参数
+ */
+void lopsided(const int ex[3],
+              const double *X, const double *Y, const double *Z,
+              const double *f, double *f_rhs,
+              const double *Sfx, const double *Sfy, const double *Sfz,
+              int Symmetry, const double SoA[3])
+{
+    const double ZEO = 0.0, ONE = 1.0, F3 = 3.0;
+    const double TWO = 2.0, F6 = 6.0, F18 = 18.0;
+    const double F12 = 12.0, F10 = 10.0, EIT = 8.0;
+
+    const int NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2;
+    (void)OCTANT; // 这里和 Fortran 一样只是定义了不用也没关系
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    // 对应 Fortran: dX = X(2)-X(1)  （Fortran 1-based）
+    // C: X[1]-X[0]
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+
+    const double d12dx = ONE / F12 / dX;
+    const double d12dy = ONE / F12 / dY;
+    const double d12dz = ONE / F12 / dZ;
+
+    // Fortran 里算了 d2dx/d2dy/d2dz 但本 subroutine 里没用到（保持一致也算出来）
+    const double d2dx  = ONE / TWO / dX;
+    const double d2dy  = ONE / TWO / dY;
+    const double d2dz  = ONE / TWO / dZ;
+    (void)d2dx; (void)d2dy; (void)d2dz;
+
+    // Fortran:
+    // imax = ex(1); jmax = ex(2); kmax = ex(3)
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    // Fortran:
+    // imin=jmin=kmin=1; 若满足对称条件则设为 -2
+    int iminF = 1, jminF = 1, kminF = 1;
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -2;
+
+    // 分配 fh：大小 (ex1+3)*(ex2+3)*(ex3+3)
+    const size_t nx = (size_t)ex1 + 3;
+    const size_t ny = (size_t)ex2 + 3;
+    const size_t nz = (size_t)ex3 + 3;
+    const size_t fh_size = nx * ny * nz;
+
+    double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return; // 内存不足：直接返回（你也可以改成 abort/报错）
+
+    // Fortran: call symmetry_bd(3,ex,f,fh,SoA)
+    symmetry_bd(3, ex, f, fh, SoA);
+
+    /*
+     * Fortran 主循环：
+     * do k=1,ex(3)-1
+     * do j=1,ex(2)-1
+     * do i=1,ex(1)-1
+     *
+     * 转成 C 0-based：
+     * k0 = 0..ex3-2, j0 = 0..ex2-2, i0 = 0..ex1-2
+     *
+     * 并且 Fortran 里的 i/j/k 在 fh 访问时，仍然是 Fortran 索引值：
+     * iF=i0+1, jF=j0+1, kF=k0+1
+     */
+    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
+                const int iF = i0 + 1;
+
+                const size_t p = idx_ex(i0, j0, k0, ex);
+
+                // ---------------- x direction ----------------
+                const double sfx = Sfx[p];
+                if (sfx > ZEO) {
+                    // Fortran: if(i+3 <= imax)
+                    // iF+3 <= ex1  <=> i0+4 <= ex1 <=> i0 <= ex1-4
+                    if (i0 <= ex1 - 4) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                    // elseif(i+2 <= imax)  <=> i0 <= ex1-3
+                    else if (i0 <= ex1 - 3) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    }
+                    // elseif(i+1 <= imax)  <=> i0 <= ex1-2（循环里总成立）
+                    else if (i0 <= ex1 - 2) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                } else if (sfx < ZEO) {
+                    // Fortran: if(i-3 >= imin)
+                    // (iF-3) >= iminF  <=> (i0-2) >= iminF
+                    if ((i0 - 2) >= iminF) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                    // elseif(i-2 >= imin) <=> (i0-1) >= iminF
+                    else if ((i0 - 1) >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    }
+                    // elseif(i-1 >= imin) <=> i0 >= iminF
+                    else if (i0 >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                }
+
+                // ---------------- y direction ----------------
+                const double sfy = Sfy[p];
+                if (sfy > ZEO) {
+                    // jF+3 <= ex2 <=> j0+4 <= ex2 <=> j0 <= ex2-4
+                    if (j0 <= ex2 - 4) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    } else if (j0 <= ex2 - 3) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 <= ex2 - 2) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    }
+                } else if (sfy < ZEO) {
+                    if ((j0 - 2) >= jminF) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    } else if ((j0 - 1) >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    }
+                }
+
+                // ---------------- z direction ----------------
+                const double sfz = Sfz[p];
+                if (sfz > ZEO) {
+                    if (k0 <= ex3 - 4) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    } else if (k0 <= ex3 - 3) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 <= ex3 - 2) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    }
+                } else if (sfz < ZEO) {
+                    if ((k0 - 2) >= kminF) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    } else if ((k0 - 1) >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    }
+                }
+            }
+        }
+    }
+    free(fh);
+}
+
+
+
+
+

AMSS_NCKU_source/lopsidediff.f90

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

AMSS_NCKU_source/macrodef.fh

@@ -1,83 +1,77 @@
 
-
-#if 0
-note here
-v:r; u: phi; w: theta
-tetradtype 0
-v^a = (x,y,z)
-orthonormal order: v,u,w
-m = (phi - i theta)/sqrt(2) following Frans, Eq.(8) of  PRD 75, 124018(2007)
-tetradtype 1
-orthonormal order: w,u,v
-m = (theta + i phi)/sqrt(2) following Sperhake, Eq.(3.2) of  PRD 85, 124062(2012)
-tetradtype 2
-v_a = (x,y,z)
-orthonormal order: v,u,w
-m = (phi - i theta)/sqrt(2) following Frans, Eq.(8) of  PRD 75, 124018(2007)
-#endif
 #define tetradtype 2
 
-#if 0
-note here
-Cell center or Vertex center
-#endif
 #define Cell
 
-#if 0
-note here
-2nd order: 2
-4th order: 3
-6th order: 4
-8th order: 5
-#endif
 #define ghost_width 3
 
-#if 0
-note here
-use shell or not
-#endif
-#define WithShell
 
-#if 0
-note here
-use constraint preserving boundary condition or not
-only affect Z4c
-#endif
-#define CPBC
 
-#if 0
-note here
-Gauge condition type
-0: B^i gauge
-1: David's puncture gauge
-2: MB B^i gauge
-3: RIT B^i gauge
-4: MB beta gauge (beta gauge not means Eq.(3) of PRD 84, 124006)
-5: RIT beta gauge (beta gauge not means Eq.(3) of PRD 84, 124006)
-6: MGB1 B^i gauge
-7: MGB2 B^i gauge
-#endif
-#define GAUGE 2
+#define GAUGE 0
 
-#if 0
-buffer points for CPBC boundary
-#endif
 #define CPBC_ghost_width  (ghost_width)
 
-#if 0
-using BSSN variable for constraint violation and psi4 calculation: 0
-using ADM variable for constraint violation and psi4 calculation: 1
-#endif
 #define ABV 0
 
-#if 0
-Type of Potential and Scalar Distribution in F(R) Scalar-Tensor Theory
-1: Case C of 1112.3928, V=0
-2: shell with a2^2*phi0/(1+a2^2), f(R) = R+a2*R^2 induced V
-3: ground state of Schrodinger-Newton system, f(R) = R+a2*R^2 induced V
-4: a2 = oo and phi(r) = phi0 * 0.5 * ( tanh((r+r0)/sigma) - tanh((r-r0)/sigma) )
-5: shell with phi(r) = phi0*Exp(-(r-r0)**2/sigma), V = 0
-#endif
 #define EScalar_CC 2
 
+#if 0
+
+define tetradtype
+    v:r; u: phi; w: theta
+    tetradtype 0
+    v^a = (x,y,z)
+    orthonormal order: v,u,w
+    m = (phi - i theta)/sqrt(2) following Frans, Eq.(8) of  PRD 75, 124018(2007)
+    tetradtype 1
+    orthonormal order: w,u,v
+    m = (theta + i phi)/sqrt(2) following Sperhake, Eq.(3.2) of  PRD 85, 124062(2012)
+    tetradtype 2
+    v_a = (x,y,z)
+    orthonormal order: v,u,w
+    m = (phi - i theta)/sqrt(2) following Frans, Eq.(8) of  PRD 75, 124018(2007)
+
+define Cell or Vertex
+    Cell center or Vertex center
+
+define ghost_width
+    2nd order: 2
+    4th order: 3
+    6th order: 4
+    8th order: 5
+
+define WithShell
+    use shell or not
+
+define CPBC
+    use constraint preserving boundary condition or not
+    only affect Z4c
+    CPBC only supports WithShell
+
+define GAUGE
+    0: B^i gauge
+    1: David puncture gauge
+    2: MB B^i gauge
+    3: RIT B^i gauge
+    4: MB beta gauge (beta gauge not means Eq.(3) of PRD 84, 124006)
+    5: RIT beta gauge (beta gauge not means Eq.(3) of PRD 84, 124006)
+    6: MGB1 B^i gauge
+    7: MGB2 B^i gauge
+
+define CPBC_ghost_width  (ghost_width)
+    buffer points for CPBC boundary
+
+define ABV
+    0: using BSSN variable for constraint violation and psi4 calculation
+    1: using ADM variable for constraint violation and psi4 calculation
+
+define EScalar_CC
+Type of Potential and Scalar Distribution in F(R) Scalar-Tensor Theory
+    1: Case C of 1112.3928, V=0
+    2: shell with   phi(r) = phi0 * a2^2/(1+a2^2), f(R) = R+a2*R^2 induced V
+    3: ground state of Schrodinger-Newton system,  f(R) = R+a2*R^2 induced V
+    4: a2 = +oo and phi(r) = phi0 * 0.5 * ( tanh((r+r0)/sigma) - tanh((r-r0)/sigma) )
+    5: shell with   phi(r) = phi0 * Exp(-(r-r0)**2/sigma), V = 0
+
+#endif
 

AMSS_NCKU_source/macrodef.h

@@ -6,95 +6,127 @@
 
 // application parameters
 
-/// ****
-// sommerfeld boundary type
-// 0: bam, 1: shibata
 #define SommerType 0
 
-/// ****
-// for Using Gauss-Legendre quadrature in theta direction
 #define GaussInt
 
-/// ****
+#define ABEtype 0
+
+//#define With_AHF
+#define Psi4type 0
+
+//#define Point_Psi4
+
+#define RPS 1
+
+#define AGM 0
+
+#define RPB 0
+
+#define MAPBH 1
+
+#define PSTR 0
+
+#define REGLEV 0
+
+//#define USE_GPU
+
+//#define CHECKDETAIL
+
+//#define FAKECHECK
+
+//
+// define SommerType
+//     sommerfeld boundary type
+//     0: bam
+//     1: shibata
+//
+// define GaussInt
+//     for Using Gauss-Legendre quadrature in theta direction
+//
+// define ABEtype
 //     0: BSSN vacuum
 //     1: coupled to scalar field
 //     2: Z4c vacuum
 //     3: coupled to Maxwell field
 //
-#define ABEtype 2
-
-/// ****
+// define With_AHF
 //     using Apparent Horizon Finder
-//#define With_AHF
-
-/// ****
+//
+// define Psi4type
 //     Psi4 calculation method
 //     0: EB method
 //     1: 4-D method
 //
-#define Psi4type 0
-
-/// ****
+// define Point_Psi4
 //     for Using point psi4 or not
-//#define Point_Psi4
-
-/// ****
+//
+// define RPS
 //     RestrictProlong in Step (0) or after Step (1)
-#define RPS 1
-
-/// ****
+//
+// define AGM
 //     Enforce algebra constraint
 //     for every RK4 sub step: 0
 //     only when iter_count == 3: 1
 //     after routine Step: 2
-#define AGM 0
-
-/// ****
+//
+// define RPB
 //     Restrict Prolong using BAM style 1 or old style 0
-#define RPB 0
-
-/// ****
+//
+// define MAPBH
 //     1: move Analysis out ot 4 sub steps and treat PBH with Euler method
-#define MAPBH 1
-
-/// ****
-// parallel structure, 0: level by level, 1: considering all levels, 2: as 1 but reverse the CPU order, 3: Frank's scheme
-#define PSTR 0
-
-/// ****
+//
+// define PSTR
+//     parallel structure
+//     0: level by level
+//     1: considering all levels
+//     2: as 1 but reverse the CPU order
+//     3: Frank's scheme
+//
+// define REGLEV
 //     regrid for every level or for all levels at a time
-// 0: for every level; 1: for all
-#define REGLEV 0
-
-/// ****
+//     0: for every level;
+//     1: for all
+//
+// define USE_GPU
 //     use gpu or not
-//#define USE_GPU
-
-/// ****
+//
+// define CHECKDETAIL
 //     use checkpoint for every process
-//#define CHECKDETAIL
-
-/// ****
+//
+// define FAKECHECK
 //     use FakeCheckPrepare to write CheckPoint
-//#define FAKECHECK
+//
+
 ////================================================================
 //  some basic parameters for numerical calculation
+////================================================================
+
 #define dim 3
 
-//#define Cell or Vertex in "microdef.fh"
+//#define Cell or Vertex in "macrodef.fh" 
 
-// ******
-// buffer point number for mesh refinement interface
 #define buffer_width 6
 
-// ******
-// buffer point number shell-box interface, on shell
 #define SC_width buffer_width
-// buffer point number shell-box interface, on box
+
 #define CS_width (2*buffer_width)
 
+//
+// define Cell or Vertex in "macrodef.fh" 
+//
+// define buffer_width
+//     buffer point number for mesh refinement interface
+//
+// define SC_width buffer_width
+//     buffer point number shell-box interface, on shell
+//
+// define CS_width
+//     buffer point number shell-box interface, on box
+//
+
 #if(buffer_width < ghost_width)
-#error we always assume buffer_width>ghost_width
+#   error we always assume buffer_width>ghost_width
 #endif
 
 #define PACK 1
@@ -110,3 +142,4 @@
 #define TINY 1e-10
 
 #endif   /* MICRODEF_H */
+

AMSS_NCKU_source/makefile

@@ -2,6 +2,27 @@
 
 include makefile.inc
 
+## ABE build flags selected by PGO_MODE (set in makefile.inc, default: opt)
+##   make                        -> opt  (PGO-guided, maximum performance)
+##   make PGO_MODE=instrument    -> instrument (Phase 1: collect fresh profile data)
+PROFDATA = /home/$(shell whoami)/AMSS-NCKU/pgo_profile/default.profdata
+
+ifeq ($(PGO_MODE),instrument)
+## Phase 1: instrumentation — omit -ipo/-fp-model fast=2 for faster build and numerical stability
+CXXAPPFLAGS = -O3 -xHost -fma -fprofile-instr-generate -ipo \
+              -Dfortran3 -Dnewc -I${MKLROOT}/include
+f90appflags = -O3 -xHost -fma -fprofile-instr-generate -ipo \
+              -align array64byte -fpp -I${MKLROOT}/include
+else
+## opt (default): maximum performance with PGO profile data
+CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma -ipo \
+              -fprofile-instr-use=$(PROFDATA) \
+              -Dfortran3 -Dnewc -I${MKLROOT}/include
+f90appflags = -O3 -xHost -fp-model fast=2 -fma -ipo \
+              -fprofile-instr-use=$(PROFDATA) \
+              -align array64byte -fpp -I${MKLROOT}/include
+endif
+
 .SUFFIXES: .o .f90 .C .for .cu
 
 .f90.o:
@@ -16,13 +37,36 @@ include makefile.inc
 .cu.o:
 	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)
 
+# C rewrite of BSSN RHS kernel and helpers
+bssn_rhs_c.o: bssn_rhs_c.C
+	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+
+fderivs_c.o: fderivs_c.C
+	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+
+fdderivs_c.o: fdderivs_c.C
+	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+
+kodiss_c.o: kodiss_c.C
+	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+
+lopsided_c.o: lopsided_c.C
+	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+
+## TwoPunctureABE uses fixed optimal flags, independent of CXXAPPFLAGS (which may be PGO-instrumented)
+TP_OPTFLAGS = -O3 -xHost -fp-model fast=2 -fma -ipo -Dfortran3 -Dnewc -I${MKLROOT}/include
+
 TwoPunctures.o: TwoPunctures.C
-	${CXX} $(CXXAPPFLAGS) -qopenmp -c $< -o $@
+	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
 
 TwoPunctureABE.o: TwoPunctureABE.C
-	${CXX} $(CXXAPPFLAGS) -qopenmp -c $< -o $@
+	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
 
 # Input files
+
+# C rewrite files
+CFILES = bssn_rhs_c.o fderivs_c.o fdderivs_c.o kodiss_c.o lopsided_c.o
+
 C++FILES = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o\
            cgh.o bssn_class.o surface_integral.o ShellPatch.o\
 	   bssnEScalar_class.o perf.o Z4c_class.o NullShellPatch.o\
@@ -40,7 +84,7 @@ C++FILES_GPU = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o
 
 F90FILES = enforce_algebra.o fmisc.o initial_puncture.o prolongrestrict.o\
 	   prolongrestrict_cell.o prolongrestrict_vertex.o\
-	   rungekutta4_rout.o bssn_rhs.o diff_new.o kodiss.o kodiss_sh.o\
+	   rungekutta4_rout.o diff_new.o kodiss.o kodiss_sh.o\
 	   lopsidediff.o sommerfeld_rout.o getnp4.o diff_new_sh.o\
 	   shellfunctions.o bssn_rhs_ss.o Set_Rho_ADM.o\
            getnp4EScalar.o bssnEScalar_rhs.o bssn_constraint.o ricci_gamma.o\
@@ -63,7 +107,7 @@ TwoPunctureFILES = TwoPunctureABE.o TwoPunctures.o
 CUDAFILES = bssn_gpu.o bssn_gpu_rhs_ss.o
 
 # file dependences
-$(C++FILES) $(C++FILESGPU) $(F90FILES) $(AHFDOBJS) $(CUDAFILES): macrodef.fh
+$(C++FILES) $(C++FILES_GPU) $(F90FILES) $(CFILES) $(AHFDOBJS) $(CUDAFILES): macrodef.fh
 
 $(C++FILES): Block.h enforce_algebra.h fmisc.h initial_puncture.h macrodef.h\
 	     misc.h monitor.h MyList.h Parallel.h MPatch.h prolongrestrict.h\
@@ -86,7 +130,7 @@ $(C++FILES_GPU): Block.h enforce_algebra.h fmisc.h initial_puncture.h macrodef.h
              
 $(AHFDOBJS): cctk.h cctk_Config.h cctk_Types.h cctk_Constants.h myglobal.h
 
-$(C++FILES) $(C++FILES_GPU) $(AHFDOBJS) $(CUDAFILES): macrodef.h
+$(C++FILES) $(C++FILES_GPU) $(CFILES) $(AHFDOBJS) $(CUDAFILES): macrodef.h
 
 TwoPunctureFILES: TwoPunctures.h
 
@@ -95,14 +139,14 @@ $(CUDAFILES): bssn_gpu.h gpu_mem.h gpu_rhsSS_mem.h
 misc.o : zbesh.o
 
 # projects
-ABE: $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) 
-	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS)
+ABE: $(C++FILES) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS)
+	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS)
 	
-ABEGPU: $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
-	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)
+ABEGPU: $(C++FILES_GPU) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
+	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)
 
 TwoPunctureABE: $(TwoPunctureFILES)
-	$(CLINKER) $(CXXAPPFLAGS) -qopenmp -o $@ $(TwoPunctureFILES) $(LDLIBS)
+	$(CLINKER) $(TP_OPTFLAGS) -qopenmp -o $@ $(TwoPunctureFILES) $(LDLIBS)
 
 clean:
 	rm *.o ABE ABEGPU TwoPunctureABE make.log -f

AMSS_NCKU_source/makefile.inc

@@ -8,18 +8,12 @@ filein  = -I/usr/include/ -I${MKLROOT}/include
 
 ## Using sequential MKL (OpenMP disabled for better single-threaded performance)
 ## Added -lifcore for Intel Fortran runtime and -limf for Intel math library
-LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl
+LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl -liomp5
 
-## Aggressive optimization flags + PGO Phase 2 (profile-guided optimization)
-## -fprofile-instr-use: use collected profile data to guide optimization decisions
-##   (branch prediction, basic block layout, inlining, loop unrolling)
-PROFDATA     = /home/amss/AMSS-NCKU/pgo_profile/default.profdata
-CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma -ipo \
-               -fprofile-instr-use=$(PROFDATA) \
-               -Dfortran3 -Dnewc -I${MKLROOT}/include
-f90appflags  = -O3 -xHost -fp-model fast=2 -fma -ipo \
-               -fprofile-instr-use=$(PROFDATA) \
-               -align array64byte -fpp -I${MKLROOT}/include
+## 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
+PGO_MODE ?= opt
 f90          = ifx
 f77          = ifx
 CXX          = icpx

AMSS_NCKU_source/share_func.h

@@ -0,0 +1,146 @@
+#ifndef SHARE_FUNC_H
+#define SHARE_FUNC_H
+
+#include <stdlib.h>
+#include <stddef.h>
+#include <math.h>
+#include <stdio.h>
+/* 主网格：0-based -> 1D */
+static inline size_t idx_ex(int i0, int j0, int k0, const int ex[3]) {
+    const int ex1 = ex[0], ex2 = ex[1];
+    return (size_t)i0 + (size_t)j0 * (size_t)ex1 + (size_t)k0 * (size_t)ex1 * (size_t)ex2;
+}
+
+/*
+ * fh 对应 Fortran: fh(-1:ex1, -1:ex2, -1:ex3)
+ * ord=2 => shift=1
+ * iF/jF/kF 为 Fortran 索引（可为 -1,0,1..ex）
+ */
+static inline size_t idx_fh_F_ord2(int iF, int jF, int kF, const int ex[3]) {
+    const int shift = 1;
+    const int nx = ex[0] + 2;      // ex1 + ord
+    const int ny = ex[1] + 2;
+
+    const int ii = iF + shift;     // 0..ex1+1
+    const int jj = jF + shift;     // 0..ex2+1
+    const int kk = kF + shift;     // 0..ex3+1
+
+    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
+}
+
+/*
+ * fh 对应 Fortran: fh(-2:ex1, -2:ex2, -2:ex3)
+ * ord=3 => shift=2
+ * iF/jF/kF 是 Fortran 索引（可为负）
+ */
+static inline size_t idx_fh_F(int iF, int jF, int kF, const int ex[3]) {
+    const int shift = 2;                 // ord=3 -> -2..ex
+    const int nx = ex[0] + 3;            // ex1 + ord
+    const int ny = ex[1] + 3;
+
+    const int ii = iF + shift;           // 0..ex1+2
+    const int jj = jF + shift;           // 0..ex2+2
+    const int kk = kF + shift;           // 0..ex3+2
+
+    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
+}
+
+/*
+ * func:  (1..extc1, 1..extc2, 1..extc3)   1-based in Fortran
+ * funcc: (-ord+1..extc1, -ord+1..extc2, -ord+1..extc3) in Fortran
+ *
+ * C 里我们把：
+ *   func  视为 0-based: i0=0..extc1-1, j0=0..extc2-1, k0=0..extc3-1
+ *   funcc 用“平移下标”存为一维数组：
+ *     iF in [-ord+1..extc1]  -> ii = iF + (ord-1)  in [0..extc1+ord-1]
+ *     总长度 nx = extc1 + ord
+ *     同理 ny = extc2 + ord, nz = extc3 + ord
+ */
+
+static inline size_t idx_func0(int i0, int j0, int k0, const int extc[3]) {
+    const int nx = extc[0], ny = extc[1];
+    return (size_t)i0 + (size_t)j0 * (size_t)nx + (size_t)k0 * (size_t)nx * (size_t)ny;
+}
+
+static inline size_t idx_funcc_F(int iF, int jF, int kF, int ord, const int extc[3]) {
+    const int shift = ord - 1;          // iF = -shift .. extc1
+    const int nx = extc[0] + ord;       // [-shift..extc1] 共 extc1+ord 个
+    const int ny = extc[1] + ord;
+
+    const int ii = iF + shift;          // 0..extc1+shift
+    const int jj = jF + shift;          // 0..extc2+shift
+    const int kk = kF + shift;          // 0..extc3+shift
+
+    return (size_t)ii + (size_t)jj * (size_t)nx + (size_t)kk * (size_t)nx * (size_t)ny;
+}
+
+/*
+ * 等价于 Fortran:
+ * funcc(1:extc1,1:extc2,1:extc3)=func
+ * do i=0,ord-1
+ *   funcc(-i,1:extc2,1:extc3) = funcc(i+1,1:extc2,1:extc3)*SoA(1)
+ * enddo
+ * do i=0,ord-1
+ *   funcc(:,-i,1:extc3) = funcc(:,i+1,1:extc3)*SoA(2)
+ * enddo
+ * do i=0,ord-1
+ *   funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
+ * enddo
+ */
+static inline void symmetry_bd(int ord,
+                 const int extc[3],
+                 const double *func,
+                 double *funcc,
+                 const double SoA[3])
+{
+    const int extc1 = extc[0], extc2 = extc[1], extc3 = extc[2];
+
+    // 1) funcc(1:extc1,1:extc2,1:extc3) = func
+    // Fortran 的 (iF=1..extc1) 对应 C 的 func(i0=0..extc1-1)
+    for (int k0 = 0; k0 < extc3; ++k0) {
+        for (int j0 = 0; j0 < extc2; ++j0) {
+            for (int i0 = 0; i0 < extc1; ++i0) {
+                const int iF = i0 + 1, jF = j0 + 1, kF = k0 + 1;
+                funcc[idx_funcc_F(iF, jF, kF, ord, extc)] = func[idx_func0(i0, j0, k0, extc)];
+            }
+        }
+    }
+
+    // 2) do i=0..ord-1: funcc(-i, 1:extc2, 1:extc3) = funcc(i+1, ...)*SoA(1)
+    for (int ii = 0; ii <= ord - 1; ++ii) {
+        const int iF_dst = -ii;       // 0, -1, -2, ...
+        const int iF_src = ii + 1;    // 1, 2, 3, ...
+        for (int kF = 1; kF <= extc3; ++kF) {
+            for (int jF = 1; jF <= extc2; ++jF) {
+                funcc[idx_funcc_F(iF_dst, jF, kF, ord, extc)] =
+                    funcc[idx_funcc_F(iF_src, jF, kF, ord, extc)] * SoA[0];
+            }
+        }
+    }
+
+    // 3) do i=0..ord-1: funcc(:,-i, 1:extc3) = funcc(:, i+1, 1:extc3)*SoA(2)
+    // 注意 Fortran 这里的 ":" 表示 iF 从 (-ord+1..extc1) 全覆盖
+    for (int jj = 0; jj <= ord - 1; ++jj) {
+        const int jF_dst = -jj;
+        const int jF_src = jj + 1;
+        for (int kF = 1; kF <= extc3; ++kF) {
+            for (int iF = -ord + 1; iF <= extc1; ++iF) {
+                funcc[idx_funcc_F(iF, jF_dst, kF, ord, extc)] =
+                    funcc[idx_funcc_F(iF, jF_src, kF, ord, extc)] * SoA[1];
+            }
+        }
+    }
+
+    // 4) do i=0..ord-1: funcc(:,:,-i) = funcc(:,:, i+1)*SoA(3)
+    for (int kk = 0; kk <= ord - 1; ++kk) {
+        const int kF_dst = -kk;
+        const int kF_src = kk + 1;
+        for (int jF = -ord + 1; jF <= extc2; ++jF) {
+            for (int iF = -ord + 1; iF <= extc1; ++iF) {
+                funcc[idx_funcc_F(iF, jF, kF_dst, ord, extc)] =
+                    funcc[idx_funcc_F(iF, jF, kF_src, ord, extc)] * SoA[2];
+            }
+        }
+    }
+}
+#endif

AMSS_NCKU_source/tool.h

@@ -0,0 +1,27 @@
+#include "share_func.h"
+void fdderivs(const int ex[3],
+              const double *f,
+              double *fxx, double *fxy, double *fxz,
+              double *fyy, double *fyz, double *fzz,
+              const double *X, const double *Y, const double *Z,
+              double SYM1, double SYM2, double SYM3,
+              int Symmetry, int onoff);
+
+void fderivs(const int ex[3],
+             const double *f,
+             double *fx, double *fy, double *fz,
+             const double *X, const double *Y, const double *Z,
+             double SYM1, double SYM2, double SYM3,
+             int Symmetry, int onoff);
+
+void kodis(const int ex[3],
+           const double *X, const double *Y, const double *Z,
+           const double *f, double *f_rhs,
+           const double SoA[3],
+           int Symmetry, double eps);
+
+void lopsided(const int ex[3],
+              const double *X, const double *Y, const double *Z,
+              const double *f, double *f_rhs,
+              const double *Sfx, const double *Sfy, const double *Sfz,
+              int Symmetry, const double SoA[3]);

makefile_and_run.py

@@ -11,17 +11,46 @@
 import AMSS_NCKU_Input as input_data
 import subprocess
 import time
-## CPU core binding configuration using taskset
-## taskset ensures all child processes inherit the CPU affinity mask
-## This forces make and all compiler processes to use only nohz_full cores (4-55, 60-111)
-## Format: taskset -c 4-55,60-111 ensures processes only run on these cores
-#NUMACTL_CPU_BIND = "taskset -c 0-111"
-NUMACTL_CPU_BIND = "taskset -c 16-47,64-95"
 
-## Build parallelism configuration
-## Use nohz_full cores (4-55, 60-111) for compilation: 52 + 52 = 104 cores
-## Set make -j to utilize available cores for faster builds
-BUILD_JOBS = 96
+
+def get_last_n_cores_per_socket(n=32):
+    """
+    Read CPU topology via lscpu and return a taskset -c string
+    selecting the last `n` cores of each NUMA node (socket).
+
+    Example: 2 sockets x 56 cores each, n=32 -> node0: 24-55, node1: 80-111
+    -> "taskset -c 24-55,80-111"
+    """
+    result = subprocess.run(["lscpu", "--parse=NODE,CPU"], capture_output=True, text=True)
+
+    # Build a dict: node_id -> sorted list of CPU ids
+    node_cpus = {}
+    for line in result.stdout.splitlines():
+        if line.startswith("#") or not line.strip():
+            continue
+        parts = line.split(",")
+        if len(parts) < 2:
+            continue
+        node_id, cpu_id = int(parts[0]), int(parts[1])
+        node_cpus.setdefault(node_id, []).append(cpu_id)
+
+    segments = []
+    for node_id in sorted(node_cpus):
+        cpus = sorted(node_cpus[node_id])
+        selected = cpus[-n:]          # last n cores of this socket
+        segments.append(f"{selected[0]}-{selected[-1]}")
+
+    cpu_str = ",".join(segments)
+    total = len(segments) * n
+    print(f" CPU binding: taskset -c {cpu_str}  ({total} cores, last {n} per socket)")
+    return f"taskset -c {cpu_str}"
+
+
+## CPU core binding: dynamically select the last 32 cores of each socket (64 cores total)
+NUMACTL_CPU_BIND = get_last_n_cores_per_socket(n=32)
+
+## Build parallelism: match the number of bound cores
+BUILD_JOBS = 64
 
 
 ##################################################################

parallel_plot_helper.py

@@ -0,0 +1,29 @@
+import multiprocessing
+
+def run_plot_task(task):
+    """Execute a single plotting task.
+    
+    Parameters
+    ----------
+    task : tuple
+        A tuple of (function, args_tuple) where function is a callable
+        plotting function and args_tuple contains its arguments.
+    """
+    func, args = task
+    return func(*args)
+
+
+def run_plot_tasks_parallel(plot_tasks):
+    """Execute a list of independent plotting tasks in parallel.
+
+    Uses the 'fork' context to create worker processes so that the main
+    script is NOT re-imported/re-executed in child processes.
+
+    Parameters
+    ----------
+    plot_tasks : list of tuples
+        Each element is (function, args_tuple).
+    """
+    ctx = multiprocessing.get_context('fork')
+    with ctx.Pool() as pool:
+        pool.map(run_plot_task, plot_tasks)

plot_GW_strain_amplitude_xiaoqu.py

@@ -11,6 +11,8 @@
 import numpy                               ## numpy for array operations
 import scipy                               ## scipy for interpolation and signal processing
 import math
+import matplotlib
+matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt     ## matplotlib for plotting
 import os                                  ## os for system/file operations
 

plot_binary_data.py

@@ -8,16 +8,23 @@
 ##
 #################################################
 
+## Restrict OpenMP to one thread per process so that running
+## many workers in parallel does not create an O(workers * BLAS_threads)
+## thread explosion.  The variable MUST be set before numpy/scipy
+## are imported, because the BLAS library reads them only at load time.
+import os
+os.environ.setdefault("OMP_NUM_THREADS",        "1")
+
 import numpy
 import scipy
+import matplotlib
+matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt
 from   matplotlib.colors    import LogNorm
 from   mpl_toolkits.mplot3d import Axes3D
 ## import torch
 import AMSS_NCKU_Input      as input_data
 
-import os
-
 
 #########################################################################################
 
@@ -192,3 +199,19 @@ def get_data_xy( Rmin, Rmax, n, data0, time, figure_title, figure_outdir ):
 
 ####################################################################################
 
+
+####################################################################################
+## Allow this module to be run as a standalone script so that each
+## binary-data plot can be executed in a fresh subprocess whose BLAS
+## environment variables (set above) take effect before numpy loads.
+##
+## Usage:  python3 plot_binary_data.py <filename> <binary_outdir> <figure_outdir>
+####################################################################################
+
+if __name__ == '__main__':
+    import sys
+    if len(sys.argv) != 4:
+        print(f"Usage: {sys.argv[0]} <filename> <binary_outdir> <figure_outdir>")
+        sys.exit(1)
+    plot_binary_data(sys.argv[1], sys.argv[2], sys.argv[3])
+

plot_xiaoqu.py

@@ -8,6 +8,8 @@
 #################################################
 
 import numpy                               ## numpy for array operations
+import matplotlib
+matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt     ## matplotlib for plotting
 from   mpl_toolkits.mplot3d import Axes3D  ## needed for 3D plots
 import glob
@@ -15,6 +17,9 @@ import os                                  ## operating system utilities
 
 import plot_binary_data
 import AMSS_NCKU_Input as input_data
+import subprocess
+import sys
+import multiprocessing
 
 # plt.rcParams['text.usetex'] = True  ## enable LaTeX fonts in plots
 
@@ -50,10 +55,40 @@ def generate_binary_data_plot( binary_outdir, figure_outdir ):
         file_list.append(x)
         print(x)
 
-    ## Plot each file in the list
+    ## Plot each file in parallel using subprocesses.
+    ## Each subprocess is a fresh Python process where the BLAS thread-count
+    ## environment variables (set at the top of plot_binary_data.py) take
+    ## effect before numpy is imported.  This avoids the thread explosion
+    ## that occurs when multiprocessing.Pool with 'fork' context inherits
+    ## already-initialized multi-threaded BLAS from the parent.
+    script = os.path.join( os.path.dirname(__file__), "plot_binary_data.py" )
+    max_workers = min( multiprocessing.cpu_count(), len(file_list) ) if file_list else 0
+
+    running = []
+    failed  = []
     for filename in file_list:
         print(filename)
-        plot_binary_data.plot_binary_data(filename, binary_outdir, figure_outdir)
+        proc = subprocess.Popen(
+            [sys.executable, script, filename, binary_outdir, figure_outdir],
+        )
+        running.append( (proc, filename) )
+        ## Keep at most max_workers subprocesses active at a time
+        if len(running) >= max_workers:
+            p, fn = running.pop(0)
+            p.wait()
+            if p.returncode != 0:
+                failed.append(fn)
+
+    ## Wait for all remaining subprocesses to finish
+    for p, fn in running:
+        p.wait()
+        if p.returncode != 0:
+            failed.append(fn)
+
+    if failed:
+        print( " WARNING: the following binary data plots failed:" )
+        for fn in failed:
+            print( "   ", fn )
 
     print(                        )
     print( " Binary Data Plot Has been Finished " )