修改block划分，对负载高的rank所在block进行划分，添加到空rank，空rank是平移得到的

Rank [  20]: Calc 0.000003 s

Rank [  35]: Calc 0.000003 s

Rank [  10]: Calc 0.000010 s

Rank [  17]: Calc 0.000005 s

Rank [  32]: Calc 0.000003 s，而且rank不固定，一般就是rank 10 和 rank 52；
但尽管有很多，比前者时间还是少很多
对于第二个Surf_Wave 调用 Interp_points，我发现前四个rank时间最长，比较固定，就是下面四个rank

Rank [  27]: Calc 0.331978 s

Rank [  35]: Calc 0.242219 s

Rank [  28]: Calc 0.242132 s

Rank [  36]: Calc 0.197024 s
因此下面surf_wave是核心

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)
+
 
 ##################################################################
 
@@ -58,7 +66,8 @@ if os.path.exists(File_directory):
     ## Prompt whether to overwrite the existing directory
     while True:
         try:
-            inputvalue = input()
+            ## inputvalue = input()
+            inputvalue = "continue"
             ## If the user agrees to overwrite, proceed and remove the existing directory
             if ( inputvalue == "continue" ):
                 print( " Continue the calculation !!! " )
@@ -424,26 +433,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/MPatch.C

@@ -442,7 +442,6 @@ void Patch::Interp_Points(MyList<var> *VarList,
       Bp = Bp->next;
     }
   }
-
   // Replace MPI_Allreduce with per-owner MPI_Bcast:
   // Group consecutive points by owner rank and broadcast each group.
   // Since each point's data is non-zero only on the owner rank,
@@ -507,6 +506,9 @@ void Patch::Interp_Points(MyList<var> *VarList,
   // Targeted point-to-point overload: each owner sends each point only to
   // the one rank that needs it for integration (consumer), reducing
   // communication volume by ~nprocs times compared to the Bcast version.
+  /*
+  double t_calc_end, t_calc_total = 0;
+  double t_calc_start = MPI_Wtime();*/
   int myrank, nprocs;
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
@@ -607,7 +609,9 @@ void Patch::Interp_Points(MyList<var> *VarList,
       Bp = Bp->next;
     }
   }
-
+  /*
+      t_calc_end = MPI_Wtime();
+      t_calc_total = t_calc_end - t_calc_start;*/
   // --- Error check for unfound points ---
   for (int j = 0; j < NN; j++)
   {
@@ -764,6 +768,63 @@ void Patch::Interp_Points(MyList<var> *VarList,
   delete[] recv_count;
   delete[] consumer_rank;
   delete[] owner_rank;
+  /*
+  // 4. 汇总并输出真正干活最慢的 Top 4
+  struct RankStats {
+    int rank;
+    double calc_time; // 净计算时间
+  };
+
+  // 创建当前进程的统计数据
+  RankStats local_stat;
+  local_stat.rank = myrank;
+  local_stat.calc_time = t_calc_total;
+
+  // 为所有进程的统计数据分配内存
+  RankStats *all_stats = nullptr;
+  if (myrank == 0) {
+    all_stats = new RankStats[nprocs];
+  }
+
+  // 使用MPI_Gather收集所有进程的数据到rank 0
+  MPI_Gather(&local_stat, sizeof(RankStats), MPI_BYTE,
+             all_stats, sizeof(RankStats), MPI_BYTE,
+             0, MPI_COMM_WORLD);
+
+  // 准备输出前4个rank的信息（所有rank都参与，确保广播后一致）
+  int top10_ranks[10] = { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 };
+  double top10_times[10] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+  int num_top10 = 0;
+
+  if (myrank == 0) {
+    // 按 calc_time（净计算时间）排序
+    std::sort(all_stats, all_stats + nprocs, [](const RankStats& a, const RankStats& b) {
+        return a.calc_time > b.calc_time;
+    });
+
+    // 取前4个
+    num_top10 = (nprocs < 10) ? nprocs : 10;
+    for (int i = 0; i < num_top10; i++) {
+      top10_ranks[i] = all_stats[i].rank;
+      top10_times[i] = all_stats[i].calc_time;
+    }
+
+    printf("\n--- Top %d Ranks by ACTIVE COMPUTATION (CPU Time) ---\n", num_top10);
+    for (int i = 0; i < num_top10; i++) {
+      printf("Rank [%4d]: Calc %.6f s\n", top10_ranks[i], top10_times[i]);
+    }
+
+    // 清理分配的内存
+    delete[] all_stats;
+  }
+
+  // 广播前4个rank的信息给所有进程
+  MPI_Bcast(&num_top10, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  if (num_top10 > 0) {
+    MPI_Bcast(top10_ranks, 10, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Bcast(top10_times, 10, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  }
+*/
 }
 void Patch::Interp_Points(MyList<var> *VarList,
                           int NN, double **XX,

AMSS_NCKU_source/NullShellPatch.h

@@ -24,6 +24,7 @@ using namespace std;
 #endif
 
 #include <mpi.h>
+#include <memory.h>
 #include "MyList.h"
 #include "Block.h"
 #include "Parallel.h"

AMSS_NCKU_source/Parallel.C

@@ -4,6 +4,7 @@
 #include "prolongrestrict.h"
 #include "misc.h"
 #include "parameters.h"
+#include <set>
 
 int Parallel::partition1(int &nx, int split_size, int min_width, int cpusize, int shape) // special for 1 diemnsion
 {
@@ -115,7 +116,7 @@ int Parallel::partition3(int *nxyz, int split_size, int *min_width, int cpusize,
   return nx * ny * nz;
 #undef SEARCH_SIZE
 }
-#elif 1 // Zhihui's idea one on 2013-09-25
+#elif 0 // Zhihui's idea one on 2013-09-25
 {
   int nx, ny, nz;
   int hmin_width;
@@ -150,7 +151,7 @@ int Parallel::partition3(int *nxyz, int split_size, int *min_width, int cpusize,
 
   return nx * ny * nz;
 }
-#elif 1 // Zhihui's idea two on 2013-09-25
+#elif 0 // Zhihui's idea two on 2013-09-25
 {
   int nx, ny, nz;
   const int hmin_width = 8; // for example we use 8
@@ -500,6 +501,428 @@ MyList<Block> *Parallel::distribute(MyList<Patch> *PatchLIST, int cpusize, int i
 
   return BlL;
 }
+MyList<Block> *Parallel::distribute_hard(MyList<Patch> *PatchLIST, int cpusize, int ingfsi, int fngfsi,
+                                    bool periodic, int nodes)
+{
+#ifdef USE_GPU_DIVIDE
+  double cpu_part, gpu_part;
+  map<string, double>::iterator iter;
+  iter = parameters::dou_par.find("cpu part");
+  if (iter != parameters::dou_par.end())
+  {
+    cpu_part = iter->second;
+  }
+  else
+  {
+    int myrank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    // read parameter from file
+    const int LEN = 256;
+    char pline[LEN];
+    string str, sgrp, skey, sval;
+    int sind;
+    char pname[50];
+    {
+      map<string, string>::iterator iter = parameters::str_par.find("inputpar");
+      if (iter != parameters::str_par.end())
+      {
+        strcpy(pname, (iter->second).c_str());
+      }
+      else
+      {
+        cout << "Error inputpar" << endl;
+        exit(0);
+      }
+    }
+    ifstream inf(pname, ifstream::in);
+    if (!inf.good() && myrank == 0)
+    {
+      cout << "Can not open parameter file " << pname << endl;
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+
+    for (int i = 1; inf.good(); i++)
+    {
+      inf.getline(pline, LEN);
+      str = pline;
+
+      int status = misc::parse_parts(str, sgrp, skey, sval, sind);
+      if (status == -1)
+      {
+        cout << "error reading parameter file " << pname << " in line " << i << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
+      }
+      else if (status == 0)
+        continue;
+
+      if (sgrp == "ABE")
+      {
+        if (skey == "cpu part")
+          cpu_part = atof(sval.c_str());
+      }
+    }
+    inf.close();
+
+    parameters::dou_par.insert(map<string, double>::value_type("cpu part", cpu_part));
+  }
+  iter = parameters::dou_par.find("gpu part");
+  if (iter != parameters::dou_par.end())
+  {
+    gpu_part = iter->second;
+  }
+  else
+  {
+    int myrank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    // read parameter from file
+    const int LEN = 256;
+    char pline[LEN];
+    string str, sgrp, skey, sval;
+    int sind;
+    char pname[50];
+    {
+      map<string, string>::iterator iter = parameters::str_par.find("inputpar");
+      if (iter != parameters::str_par.end())
+      {
+        strcpy(pname, (iter->second).c_str());
+      }
+      else
+      {
+        cout << "Error inputpar" << endl;
+        exit(0);
+      }
+    }
+    ifstream inf(pname, ifstream::in);
+    if (!inf.good() && myrank == 0)
+    {
+      cout << "Can not open parameter file " << pname << endl;
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+
+    for (int i = 1; inf.good(); i++)
+    {
+      inf.getline(pline, LEN);
+      str = pline;
+
+      int status = misc::parse_parts(str, sgrp, skey, sval, sind);
+      if (status == -1)
+      {
+        cout << "error reading parameter file " << pname << " in line " << i << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
+      }
+      else if (status == 0)
+        continue;
+
+      if (sgrp == "ABE")
+      {
+        if (skey == "gpu part")
+          gpu_part = atof(sval.c_str());
+      }
+    }
+    inf.close();
+
+    parameters::dou_par.insert(map<string, double>::value_type("gpu part", gpu_part));
+  }
+
+  if (nodes == 0)
+    nodes = cpusize / 2;
+#else
+  if (nodes == 0)
+    nodes = cpusize;
+#endif
+
+  if (dim != 3)
+  {
+    cout << "distrivute: now we only support 3-dimension" << endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+
+  MyList<Block> *BlL = 0;
+  int split_size, min_size, block_size = 0;
+
+  int min_width = 2 * Mymax(ghost_width, buffer_width);
+  int nxyz[dim], mmin_width[dim], min_shape[dim];
+
+  MyList<Patch> *PLi = PatchLIST;
+  for (int i = 0; i < dim; i++)
+    min_shape[i] = PLi->data->shape[i];
+  int lev = PLi->data->lev;
+  PLi = PLi->next;
+  while (PLi)
+  {
+    Patch *PP = PLi->data;
+    for (int i = 0; i < dim; i++)
+      min_shape[i] = Mymin(min_shape[i], PP->shape[i]);
+    if (lev != PLi->data->lev)
+      cout << "Parallel::distribute CAUSTION: meet Patches for different level: " << lev << " and " << PLi->data->lev << endl;
+    PLi = PLi->next;
+  }
+
+  for (int i = 0; i < dim; i++)
+    mmin_width[i] = Mymin(min_width, min_shape[i]);
+
+  min_size = mmin_width[0];
+  for (int i = 1; i < dim; i++)
+    min_size = min_size * mmin_width[i];
+
+  PLi = PatchLIST;
+  while (PLi)
+  {
+    Patch *PP = PLi->data;
+    //    PP->checkPatch(true);
+    int bs = PP->shape[0];
+    for (int i = 1; i < dim; i++)
+      bs = bs * PP->shape[i];
+    block_size = block_size + bs;
+    PLi = PLi->next;
+  }
+  split_size = Mymax(min_size, block_size / nodes);
+  split_size = Mymax(1, split_size);
+
+  int n_rank = 0;
+  PLi = PatchLIST;
+  int reacpu = 0;
+  int current_block_id = 0;
+  while (PLi) {
+    Block *ng0, *ng;
+    bool first_block_in_patch = true; 
+    Patch *PP = PLi->data;
+    reacpu += partition3(nxyz, split_size, mmin_width, nodes, PP->shape);
+
+    for (int i = 0; i < nxyz[0]; i++)
+    for (int j = 0; j < nxyz[1]; j++)
+    for (int k = 0; k < nxyz[2]; k++)
+    {
+        // --- 1. 定义局部变量 ---
+        int ibbox_here[6], shape_here[3];
+        double bbox_here[6], dd;
+        Block *current_ng_start = nullptr; // 本次循环产生的第一个(或唯一一个)块
+
+        // --- 2. 核心逻辑分支 ---
+        if (current_block_id == 27 || current_block_id == 28 ||
+            current_block_id == 35 || current_block_id == 36)
+        {
+            // A. 计算原始索引 (不带 Ghost)
+            int ib0 = (PP->shape[0] * i) / nxyz[0];
+            int ib3 = (PP->shape[0] * (i + 1)) / nxyz[0] - 1;
+            int jb1 = (PP->shape[1] * j) / nxyz[1];
+            int jb4 = (PP->shape[1] * (j + 1)) / nxyz[1] - 1;
+            int kb2 = (PP->shape[2] * k) / nxyz[2];
+            int kb5 = (PP->shape[2] * (k + 1)) / nxyz[2] - 1;
+
+            int r_l, r_r;
+            if(current_block_id == 27)      { r_l = 26; r_r = 27; }
+            else if(current_block_id == 28) { r_l = 28; r_r = 29; }
+            else if(current_block_id == 35) { r_l = 34; r_r = 35; }
+            else                            { r_l = 36; r_r = 37; }
+            Block * split_first_block = nullptr;
+            Block * split_last_block = nullptr; 
+            // 拆分逻辑：该函数应更新类成员变量 split_first_block 和 split_last_block
+            splitHotspotBlock(BlL, dim, ib0, ib3, jb1, jb4, kb2, kb5, 
+                              PP, r_l, r_r, ingfsi, fngfsi, periodic,split_first_block,split_last_block);
+            
+            current_ng_start = split_first_block;
+            ng = split_last_block; 
+        }
+        else 
+        {
+            // B. 普通块逻辑 (含 Ghost 扩张)
+            ibbox_here[0] = (PP->shape[0] * i) / nxyz[0];
+            ibbox_here[3] = (PP->shape[0] * (i + 1)) / nxyz[0] - 1;
+            ibbox_here[1] = (PP->shape[1] * j) / nxyz[1];
+            ibbox_here[4] = (PP->shape[1] * (j + 1)) / nxyz[1] - 1;
+            ibbox_here[2] = (PP->shape[2] * k) / nxyz[2];
+            ibbox_here[5] = (PP->shape[2] * (k + 1)) / nxyz[2] - 1;
+
+            if (periodic) {
+                for(int d=0; d<3; d++) {
+                    ibbox_here[d] -= ghost_width;
+                    ibbox_here[d+3] += ghost_width;
+                }
+            } else {
+                ibbox_here[0] = Mymax(0, ibbox_here[0] - ghost_width);
+                ibbox_here[3] = Mymin(PP->shape[0] - 1, ibbox_here[3] + ghost_width);
+                ibbox_here[1] = Mymax(0, ibbox_here[1] - ghost_width);
+                ibbox_here[4] = Mymin(PP->shape[1] - 1, ibbox_here[4] + ghost_width);
+                ibbox_here[2] = Mymax(0, ibbox_here[2] - ghost_width);
+                ibbox_here[5] = Mymin(PP->shape[2] - 1, ibbox_here[5] + ghost_width);
+            }
+
+            for(int d=0; d<3; d++) shape_here[d] = ibbox_here[d+3] - ibbox_here[d] + 1;
+
+            // 物理坐标计算 (根据你的宏定义 Cell/Vertex)
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+          // 0--4, 5--10
+          dd = (PP->bbox[3] - PP->bbox[0]) / (PP->shape[0] - 1);
+          bbox_here[0] = PP->bbox[0] + ibbox_here[0] * dd;
+          bbox_here[3] = PP->bbox[0] + ibbox_here[3] * dd;
+
+          dd = (PP->bbox[4] - PP->bbox[1]) / (PP->shape[1] - 1);
+          bbox_here[1] = PP->bbox[1] + ibbox_here[1] * dd;
+          bbox_here[4] = PP->bbox[1] + ibbox_here[4] * dd;
+
+          dd = (PP->bbox[5] - PP->bbox[2]) / (PP->shape[2] - 1);
+          bbox_here[2] = PP->bbox[2] + ibbox_here[2] * dd;
+          bbox_here[5] = PP->bbox[2] + ibbox_here[5] * dd;
+#else
+#ifdef Cell
+          // 0--5, 5--10
+          dd = (PP->bbox[3] - PP->bbox[0]) / PP->shape[0];
+          bbox_here[0] = PP->bbox[0] + (ibbox_here[0]) * dd;
+          bbox_here[3] = PP->bbox[0] + (ibbox_here[3] + 1) * dd;
+
+          dd = (PP->bbox[4] - PP->bbox[1]) / PP->shape[1];
+          bbox_here[1] = PP->bbox[1] + (ibbox_here[1]) * dd;
+          bbox_here[4] = PP->bbox[1] + (ibbox_here[4] + 1) * dd;
+
+          dd = (PP->bbox[5] - PP->bbox[2]) / PP->shape[2];
+          bbox_here[2] = PP->bbox[2] + (ibbox_here[2]) * dd;
+          bbox_here[5] = PP->bbox[2] + (ibbox_here[5] + 1) * dd;
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+            ng = createMappedBlock(BlL, dim, shape_here, bbox_here, current_block_id, ingfsi, fngfsi, PP->lev);
+            current_ng_start = ng;
+        }
+
+        // --- 3. 统一处理 Patch 起始 Block 指针 ---
+        if (first_block_in_patch) {
+            ng0 = current_ng_start;
+            
+            // 立即设置 PP->blb，避免后续循环覆盖 ng0
+            MyList<Block> *Bp_start = BlL;
+            while (Bp_start && Bp_start->data != ng0) Bp_start = Bp_start->next;
+            PP->blb = Bp_start;
+            
+            first_block_in_patch = false;
+        }
+        
+        current_block_id++;
+    }
+
+    // --- 4. 设置 Patch 结束 Block 指针 ---
+    MyList<Block> *Bp_end = BlL;
+    while (Bp_end && Bp_end->data != ng) Bp_end = Bp_end->next;
+    PP->ble = Bp_end;
+
+    PLi = PLi->next;
+    first_block_in_patch = true; 
+  }
+  if (reacpu < nodes * 2 / 3)
+  {
+    int myrank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    if (myrank == 0)
+      cout << "Parallel::distribute CAUSTION: level#" << lev << " uses essencially " << reacpu << " processors vs " << nodes << " nodes run, your scientific computation scale is not as large as you estimate." << endl;
+  }
+
+  return BlL;
+}
+
+ /**
+ * @brief 将当前 Block 几何二等分并存入列表
+ * @param axis 拆分轴：0-x, 1-y, 2-z (建议选最长轴)
+ */
+Block* Parallel::splitHotspotBlock(MyList<Block>* &BlL, int _dim, 
+                                 int ib0_orig, int ib3_orig, 
+                                 int jb1_orig, int jb4_orig, 
+                                 int kb2_orig, int kb5_orig, 
+                                 Patch* PP, int r_left, int r_right, 
+                                 int ingfsi, int fngfsi, bool periodic,
+                                 Block* &split_first_block, Block* &split_last_block)
+{
+    // 1. 索引二分 (基于无 ghost 的原始索引)
+    int mid = (ib0_orig + ib3_orig) / 2;
+
+    // 左块原始索引: [ib0, mid], 右块原始索引: [mid+1, ib3]
+    int indices_L[6] = {ib0_orig, jb1_orig, kb2_orig, mid, jb4_orig, kb5_orig};
+    int indices_R[6] = {mid + 1, jb1_orig, kb2_orig, ib3_orig, jb4_orig, kb5_orig};
+
+    // 2. 内部处理逻辑 (复刻原 distribute 逻辑)
+    auto createSubBlock = [&](int* ib_raw, int target_rank) {
+        int ib_final[6];
+        int sh_here[3];
+        double bb_here[6], dd;
+
+        // --- 逻辑 A: Ghost 扩张 ---
+        if (periodic) {
+            ib_final[0] = ib_raw[0] - ghost_width;
+            ib_final[3] = ib_raw[3] + ghost_width;
+            ib_final[1] = ib_raw[1] - ghost_width;
+            ib_final[4] = ib_raw[4] + ghost_width;
+            ib_final[2] = ib_raw[2] - ghost_width;
+            ib_final[5] = ib_raw[5] + ghost_width;
+        } else {
+            ib_final[0] = Mymax(0, ib_raw[0] - ghost_width);
+            ib_final[3] = Mymin(PP->shape[0] - 1, ib_raw[3] + ghost_width);
+            ib_final[1] = Mymax(0, ib_raw[1] - ghost_width);
+            ib_final[4] = Mymin(PP->shape[1] - 1, ib_raw[4] + ghost_width);
+            ib_final[2] = Mymax(0, ib_raw[2] - ghost_width);
+            ib_final[5] = Mymin(PP->shape[2] - 1, ib_raw[5] + ghost_width);
+        }
+
+        sh_here[0] = ib_final[3] - ib_final[0] + 1;
+        sh_here[1] = ib_final[4] - ib_final[1] + 1;
+        sh_here[2] = ib_final[5] - ib_final[2] + 1;
+
+        // --- 逻辑 B: 物理坐标计算 (严格匹配 Cell 模式) ---
+        // X 方向
+        dd = (PP->bbox[3] - PP->bbox[0]) / PP->shape[0];
+        bb_here[0] = PP->bbox[0] + ib_final[0] * dd;
+        bb_here[3] = PP->bbox[0] + (ib_final[3] + 1) * dd;
+
+        // Y 方向
+        dd = (PP->bbox[4] - PP->bbox[1]) / PP->shape[1];
+        bb_here[1] = PP->bbox[1] + ib_final[1] * dd;
+        bb_here[4] = PP->bbox[1] + (ib_final[4] + 1) * dd;
+
+        // Z 方向
+        dd = (PP->bbox[5] - PP->bbox[2]) / PP->shape[2];
+        bb_here[2] = PP->bbox[2] + ib_final[2] * dd;
+        bb_here[5] = PP->bbox[2] + (ib_final[5] + 1) * dd;
+
+        Block* Bg = new Block(dim, sh_here, bb_here, target_rank, ingfsi, fngfsi, PP->lev);
+        if (BlL) BlL->insert(Bg);
+        else     BlL = new MyList<Block>(Bg);
+
+       
+        return Bg;
+    };
+
+    // 执行创建
+    split_first_block = createSubBlock(indices_L, r_left);
+    split_last_block  = createSubBlock(indices_R, r_right);
+}
+
+
+/**
+ * @brief 创建映射后的 Block
+ */
+  Block* Parallel::createMappedBlock(MyList<Block>* &BlL, int _dim, int* shape, double* bbox,
+                        int block_id, int ingfsi, int fngfsi, int lev)
+  {
+      // 映射表逻辑
+      int target_rank = block_id;
+      if (block_id == 26)      target_rank = 25;
+      else if (block_id == 29) target_rank = 30;
+      else if (block_id == 34) target_rank = 33;
+      else if (block_id == 37) target_rank = 38;
+
+      Block* ng = new Block(dim, shape, bbox, target_rank, ingfsi, fngfsi, lev);
+
+      if (BlL) BlL->insert(ng);
+      else     BlL = new MyList<Block>(ng);
+
+      return ng;
+  }
+
+
+
+
 #elif (PSTR == 1 || PSTR == 2 || PSTR == 3)
 MyList<Block> *Parallel::distribute(MyList<Patch> *PatchLIST, int cpusize, int ingfsi, int fngfsi,
                                     bool periodic, int start_rank, int end_rank, int nodes)
@@ -3853,7 +4276,8 @@ void Parallel::Sync_merged(MyList<Patch> *PatL, MyList<var> *VarList, int Symmet
 Parallel::SyncCache::SyncCache()
     : valid(false), cpusize(0), combined_src(0), combined_dst(0),
       send_lengths(0), recv_lengths(0), send_bufs(0), recv_bufs(0),
-      send_buf_caps(0), recv_buf_caps(0), reqs(0), stats(0), max_reqs(0)
+      send_buf_caps(0), recv_buf_caps(0), reqs(0), stats(0), max_reqs(0),
+      lengths_valid(false)
 {
 }
 // SyncCache invalidate: free grid segment lists but keep buffers
@@ -3871,6 +4295,7 @@ void Parallel::SyncCache::invalidate()
     send_lengths[i] = recv_lengths[i] = 0;
   }
   valid = false;
+  lengths_valid = false;
 }
 // SyncCache destroy: free everything
 void Parallel::SyncCache::destroy()
@@ -4172,8 +4597,13 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
   {
     if (node == myrank)
     {
-      int length = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
-      cache.recv_lengths[node] = length;
+      int length;
+      if (!cache.lengths_valid) {
+        length = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
+        cache.recv_lengths[node] = length;
+      } else {
+        length = cache.recv_lengths[node];
+      }
       if (length > 0)
       {
         if (length > cache.recv_buf_caps[node])
@@ -4187,8 +4617,13 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
     }
     else
     {
-      int slength = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
-      cache.send_lengths[node] = slength;
+      int slength;
+      if (!cache.lengths_valid) {
+        slength = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
+        cache.send_lengths[node] = slength;
+      } else {
+        slength = cache.send_lengths[node];
+      }
       if (slength > 0)
       {
         if (slength > cache.send_buf_caps[node])
@@ -4200,8 +4635,13 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
         data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList, VarList, Symmetry);
         MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 2, MPI_COMM_WORLD, cache.reqs + state.req_no++);
       }
-      int rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList, VarList, Symmetry);
-      cache.recv_lengths[node] = rlength;
+      int rlength;
+      if (!cache.lengths_valid) {
+        rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList, VarList, Symmetry);
+        cache.recv_lengths[node] = rlength;
+      } else {
+        rlength = cache.recv_lengths[node];
+      }
       if (rlength > 0)
       {
         if (rlength > cache.recv_buf_caps[node])
@@ -4214,6 +4654,7 @@ void Parallel::Sync_start(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetr
       }
     }
   }
+  cache.lengths_valid = true;
 }
 // Sync_finish: wait for async MPI operations and unpack
 void Parallel::Sync_finish(SyncCache &cache, AsyncSyncState &state,
@@ -5268,6 +5709,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)
 {
@@ -6267,3 +6905,224 @@ void Parallel::checkpatchlist(MyList<Patch> *PatL, bool buflog)
     PL = PL->next;
   }
 }
+// Check if load balancing is needed based on interpolation times
+bool Parallel::check_load_balance_need(double *rank_times, int nprocs, int &num_heavy, int *heavy_ranks)
+{
+  // Calculate average time
+  double avg_time = 0;
+  for (int r = 0; r < nprocs; r++)
+  {
+    avg_time += rank_times[r];
+  }
+  avg_time /= nprocs;
+
+  // Identify heavy ranks (time > 1.5x average)
+  std::vector<std::pair<int, double>> rank_times_vec;
+  for (int r = 0; r < nprocs; r++)
+  {
+    if (rank_times[r] > avg_time * 1.5)
+    {
+      rank_times_vec.push_back(std::make_pair(r, rank_times[r]));
+    }
+  }
+
+  // Sort by time (descending)
+  std::sort(rank_times_vec.begin(), rank_times_vec.end(),
+            [](const std::pair<int, double>& a, const std::pair<int, double>& b) {
+              return a.second > b.second;
+            });
+
+  // Take top 4 heavy ranks
+  num_heavy = std::min(4, (int)rank_times_vec.size());
+  if (num_heavy > 0)
+  {
+    for (int i = 0; i < num_heavy; i++)
+    {
+      heavy_ranks[i] = rank_times_vec[i].first;
+    }
+    return true;  // Load balancing is needed
+  }
+
+  return false;  // No load balancing needed
+}
+
+// Split blocks belonging to heavy ranks to improve load balancing
+// Strategy: Split heavy rank blocks in half, merge 8 light ranks to free 4 ranks
+void Parallel::split_heavy_blocks(MyList<Patch> *PatL, int *heavy_ranks, int num_heavy,
+                                  int split_factor, int cpusize, int ingfsi, int fngfsi)
+{
+  int myrank, nprocs;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+
+  if (myrank != 0) return; // Only rank 0 performs the analysis
+
+  cout << "\n=== Load Balancing Strategy ===" << endl;
+  cout << "Heavy ranks to split (in half): " << num_heavy << endl;
+  for (int i = 0; i < num_heavy; i++)
+    cout << "  Heavy rank " << heavy_ranks[i] << endl;
+
+  // Step 1: Identify all blocks and their ranks
+  std::vector<int> all_ranks;
+  std::map<int, std::vector<Block*>> rank_to_blocks;
+
+  MyList<Patch> *PL = PatL;
+  while (PL)
+  {
+    Patch *PP = PL->data;
+    MyList<Block> *BP = PP->blb;
+    while (BP)
+    {
+      Block *block = BP->data;
+      all_ranks.push_back(block->rank);
+      rank_to_blocks[block->rank].push_back(block);
+      BP = BP->next;
+    }
+    PL = PL->next;
+  }
+
+  // Step 2: Identify light ranks (not in heavy_ranks list)
+  std::set<int> heavy_set(heavy_ranks, heavy_ranks + num_heavy);
+  std::vector<int> light_ranks;
+  for (int r : all_ranks)
+  {
+    if (heavy_set.find(r) == heavy_set.end())
+    {
+      light_ranks.push_back(r);
+    }
+  }
+
+  // Remove duplicates from light_ranks
+  std::sort(light_ranks.begin(), light_ranks.end());
+  light_ranks.erase(std::unique(light_ranks.begin(), light_ranks.end()), light_ranks.end());
+
+  cout << "Found " << light_ranks.size() << " light ranks (candidates for merging)" << endl;
+
+  // Step 3: Select 8 light ranks to merge (those with smallest workload)
+  // For now, we select the first 8 light ranks
+  int num_to_merge = 8;
+  if (light_ranks.size() < num_to_merge)
+  {
+    cout << "WARNING: Not enough light ranks to merge. Found " << light_ranks.size()
+         << ", need " << num_to_merge << endl;
+    num_to_merge = light_ranks.size();
+  }
+
+  std::vector<int> ranks_to_merge(light_ranks.begin(), light_ranks.begin() + num_to_merge);
+
+  cout << "Light ranks to merge (8 -> 4 merged ranks):" << endl;
+  for (int i = 0; i < num_to_merge; i++)
+    cout << "  Rank " << ranks_to_merge[i] << endl;
+
+  // Step 4: Analyze blocks that need to be split
+  cout << "\n=== Analyzing blocks for splitting ===" << endl;
+
+  struct BlockSplitInfo {
+    Block *original_block;
+    int split_dim;
+    int split_point;
+  };
+
+  std::vector<BlockSplitInfo> blocks_to_split;
+
+  PL = PatL;
+  while (PL)
+  {
+    Patch *PP = PL->data;
+    MyList<Block> *BP = PP->blb;
+    while (BP)
+    {
+      Block *block = BP->data;
+
+      // Check if this block belongs to a heavy rank
+      for (int i = 0; i < num_heavy; i++)
+      {
+        if (block->rank == heavy_ranks[i])
+        {
+          // Find the largest dimension for splitting
+          int max_dim = 0;
+          int max_size = block->shape[0];
+          for (int d = 1; d < dim; d++)
+          {
+            if (block->shape[d] > max_size)
+            {
+              max_size = block->shape[d];
+              max_dim = d;
+            }
+          }
+
+          int split_point = max_size / 2;
+
+          BlockSplitInfo info;
+          info.original_block = block;
+          info.split_dim = max_dim;
+          info.split_point = split_point;
+          blocks_to_split.push_back(info);
+
+          cout << "Block at rank " << block->rank << " will be split" << endl;
+          cout << "  Shape: [" << block->shape[0] << ", " << block->shape[1] << ", " << block->shape[2] << "]" << endl;
+          cout << "  Split along dimension " << max_dim << " at index " << split_point << endl;
+          break;
+        }
+      }
+
+      BP = BP->next;
+    }
+    PL = PL->next;
+  }
+
+  cout << "\nTotal blocks to split: " << blocks_to_split.size() << endl;
+
+  // Step 5: Calculate new rank assignments
+  // Strategy:
+  // - For each heavy rank, its blocks are split in half
+  // - First half keeps the original rank
+  // - Second half gets a new rank (from the freed light ranks)
+  // - 8 light ranks are merged into 4 ranks, freeing up 4 ranks
+
+  std::vector<int> freed_ranks;
+  for (size_t i = 0; i < ranks_to_merge.size(); i += 2)
+  {
+    // Merge pairs of light ranks: (ranks_to_merge[i], ranks_to_merge[i+1]) -> ranks_to_merge[i]
+    // This frees up ranks_to_merge[i+1]
+    if (i + 1 < ranks_to_merge.size())
+    {
+      freed_ranks.push_back(ranks_to_merge[i + 1]);
+      cout << "Merging ranks " << ranks_to_merge[i] << " and " << ranks_to_merge[i + 1]
+           << " -> keeping rank " << ranks_to_merge[i] << ", freeing rank " << ranks_to_merge[i + 1] << endl;
+    }
+  }
+
+  cout << "\nFreed ranks available for split blocks: ";
+  for (int r : freed_ranks)
+    cout << r << " ";
+  cout << endl;
+
+  // Step 6: Assign new ranks to split blocks
+  int freed_idx = 0;
+  for (size_t i = 0; i < blocks_to_split.size(); i++)
+  {
+    BlockSplitInfo &info = blocks_to_split[i];
+    Block *original = info.original_block;
+
+    if (freed_idx < freed_ranks.size())
+    {
+      cout << "\nSplitting block at rank " << original->rank << endl;
+      cout << "  First half: keeps rank " << original->rank << endl;
+      cout << "  Second half: gets new rank " << freed_ranks[freed_idx] << endl;
+      freed_idx++;
+    }
+    else
+    {
+      cout << "WARNING: Not enough freed ranks for all split blocks!" << endl;
+      break;
+    }
+  }
+
+  cout << "\n=== Load Balancing Analysis Complete ===" << endl;
+  cout << "Next steps:" << endl;
+  cout << "  1. Recompose the grid with new rank assignments" << endl;
+  cout << "  2. Data migration will be handled by recompose_cgh" << endl;
+  cout << "  3. Ghost zone communication will be updated automatically" << endl;
+}
+

AMSS_NCKU_source/Parallel.h

@@ -11,7 +11,7 @@
 #include <cmath>
 #include <new>
 using namespace std;
-
+#include <memory.h>
 #include "Parallel_bam.h"
 #include "var.h"
 #include "MPatch.h"
@@ -32,6 +32,16 @@ namespace Parallel
   int partition2(int *nxy, int split_size, int *min_width, int cpusize, int *shape); // special for 2 diemnsions
   int partition3(int *nxyz, int split_size, int *min_width, int cpusize, int *shape);
   MyList<Block> *distribute(MyList<Patch> *PatchLIST, int cpusize, int ingfsi, int fngfs, bool periodic, int nodes = 0); // produce corresponding Blocks
+  MyList<Block> *distribute_hard(MyList<Patch> *PatchLIST, int cpusize, int ingfsi, int fngfs, bool periodic, int nodes = 0); // produce corresponding Blocks
+  Block* splitHotspotBlock(MyList<Block>* &BlL, int _dim, 
+                                 int ib0_orig, int ib3_orig, 
+                                 int jb1_orig, int jb4_orig, 
+                                 int kb2_orig, int kb5_orig, 
+                                 Patch* PP, int r_left, int r_right, 
+                                 int ingfsi, int fngfsi, bool periodic,
+                                 Block* &split_first_block, Block* &split_last_block);
+  Block* createMappedBlock(MyList<Block>* &BlL, int _dim, int* shape, double* bbox,
+                        int block_id, int ingfsi, int fngfsi, int lev); 
   void KillBlocks(MyList<Patch> *PatchLIST);
 
   void setfunction(MyList<Block> *BlL, var *vn, double func(double x, double y, double z));
@@ -97,6 +107,7 @@ namespace Parallel
     MPI_Request *reqs;
     MPI_Status *stats;
     int max_reqs;
+    bool lengths_valid;
     SyncCache();
     void invalidate();
     void destroy();
@@ -129,6 +140,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);
@@ -198,6 +218,18 @@ namespace Parallel
 #if (PSTR == 1 || PSTR == 2 || PSTR == 3)
   MyList<Block> *distribute(MyList<Patch> *PatchLIST, int cpusize, int ingfsi, int fngfsi,
                             bool periodic, int start_rank, int end_rank, int nodes = 0);
+  
+  // Redistribute blocks with time statistics for load balancing
+  MyList<Block> *distribute(MyList<Patch> *PatchLIST, MyList<Block> *OldBlockL,
+                            int cpusize, int ingfsi, int fngfsi,
+                            bool periodic, int start_rank, int end_rank, int nodes = 0);
 #endif
-}
-#endif /*PARALLEL_H */
+
+  // Dynamic load balancing: split blocks for heavy ranks
+    void split_heavy_blocks(MyList<Patch> *PatL, int *heavy_ranks, int num_heavy,
+                            int split_factor, int cpusize, int ingfsi, int fngfsi);
+
+    // Check if load balancing is needed based on interpolation times
+    bool check_load_balance_need(double *rank_times, int nprocs, int &num_heavy, int *heavy_ranks);
+  }
+  #endif /*PARALLEL_H */

AMSS_NCKU_source/bssn_class.C

@@ -5819,21 +5819,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 +5870,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 +5949,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 +5971,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 +6036,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 +6060,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 +6101,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 +6114,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

@@ -43,6 +43,14 @@ cgh::cgh(int ingfsi, int fngfsi, int Symmetry, char *filename, int checkrun,
   end_rank = 0;
 #endif
 
+  // Initialize load balancing variables
+  enable_load_balance = false;
+  load_balance_check_interval = 10;  // Check every 10 time steps
+  current_time_step = 0;
+  rank_interp_times = nullptr;
+  heavy_ranks = nullptr;
+  num_heavy_ranks = 0;
+
   if (!checkrun)
   {
     read_bbox(Symmetry, filename);
@@ -113,6 +121,12 @@ cgh::~cgh()
     delete[] Porgls[lev];
   }
   delete[] Porgls;
+
+  // Clean up load balancing memory
+  if (rank_interp_times)
+    delete[] rank_interp_times;
+  if (heavy_ranks)
+    delete[] heavy_ranks;
 }
 
 //================================================================================================
@@ -130,7 +144,7 @@ void cgh::compose_cgh(int nprocs)
   for (int lev = 0; lev < levels; lev++)
   {
     checkPatchList(PatL[lev], false);
-    Parallel::distribute(PatL[lev], nprocs, ingfs, fngfs, false);
+    Parallel::distribute_hard(PatL[lev], nprocs, ingfs, fngfs, false);
 #if (RPB == 1)
     // we need distributed box of PatL[lev] and PatL[lev-1]
     if (lev > 0)
@@ -1705,3 +1719,121 @@ void cgh::settrfls(const int lev)
 {
   trfls = lev;
 }
+
+//================================================================================================
+// Load Balancing Functions
+//================================================================================================
+
+// Initialize load balancing
+void cgh::init_load_balance(int nprocs)
+{
+  if (rank_interp_times)
+    delete[] rank_interp_times;
+  if (heavy_ranks)
+    delete[] heavy_ranks;
+
+  rank_interp_times = new double[nprocs];
+  heavy_ranks = new int[4];  // Maximum 4 heavy ranks
+  num_heavy_ranks = 0;
+
+  for (int i = 0; i < nprocs; i++)
+    rank_interp_times[i] = 0.0;
+}
+
+// Update interpolation time for a rank
+void cgh::update_interp_time(int rank, double time)
+{
+  if (rank_interp_times && rank >= 0)
+  {
+    rank_interp_times[rank] = time;
+  }
+}
+
+// Check and perform load balancing if needed
+bool cgh::check_and_rebalance(int nprocs, int lev,
+                               MyList<var> *OldList, MyList<var> *StateList,
+                               MyList<var> *FutureList, MyList<var> *tmList,
+                               int Symmetry, bool BB)
+{
+  int myrank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+  // Only check at specified intervals
+  current_time_step++;
+  if (current_time_step % load_balance_check_interval != 0)
+    return false;
+
+  if (myrank == 0)
+  {
+    cout << "\n=== Checking load balance at time step " << current_time_step << " ===" << endl;
+  }
+
+  // Collect all rank times on rank 0
+  double *all_times = nullptr;
+  if (myrank == 0)
+  {
+    all_times = new double[nprocs];
+  }
+
+  MPI_Gather(rank_interp_times, 1, MPI_DOUBLE, all_times, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+
+  bool need_rebalance = false;
+
+  if (myrank == 0)
+  {
+    // Check if load balancing is needed
+    need_rebalance = Parallel::check_load_balance_need(all_times, nprocs, num_heavy_ranks, heavy_ranks);
+
+    if (need_rebalance)
+    {
+      cout << "=== Load imbalance detected! Need to rebalance ===" << endl;
+      cout << "Top " << num_heavy_ranks << " heavy ranks: ";
+      for (int i = 0; i < num_heavy_ranks; i++)
+      {
+        cout << heavy_ranks[i] << " (" << all_times[heavy_ranks[i]] << " s) ";
+      }
+      cout << endl;
+
+      // Analyze blocks that need to be split
+      Parallel::split_heavy_blocks(PatL[lev], heavy_ranks, num_heavy_ranks, 2, nprocs, ingfs, fngfs);
+
+      // Set lev_flag to trigger recompose_cgh
+      cout << "=== Triggering recompose_cgh for level " << lev << " ===" << endl;
+    }
+    else
+    {
+      cout << "=== Load is balanced, no rebalancing needed ===" << endl;
+    }
+
+    delete[] all_times;
+  }
+
+  // Broadcast the decision to all ranks
+  MPI_Bcast(&need_rebalance, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
+
+  if (need_rebalance)
+  {
+    // Broadcast heavy ranks information
+    MPI_Bcast(&num_heavy_ranks, 1, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Bcast(heavy_ranks, num_heavy_ranks, MPI_INT, 0, MPI_COMM_WORLD);
+
+    // Perform recompose_cgh on the specified level
+    if (myrank == 0)
+    {
+      cout << "=== Performing recompose_cgh ===" << endl;
+    }
+
+    // Call recompose_cgh_Onelevel for the specified level
+    bool *lev_flag = new bool[1];
+    lev_flag[0] = true;
+    recompose_cgh_Onelevel(nprocs, lev, OldList, StateList, FutureList, tmList, Symmetry, BB);
+    delete[] lev_flag;
+
+    // Reset time counter after rebalancing
+    current_time_step = 0;
+
+    return true;
+  }
+
+  return false;
+}

AMSS_NCKU_source/cgh.h

@@ -87,6 +87,21 @@ public:
 #if (PSTR == 1 || PSTR == 2 || PSTR == 3)
    void construct_mylev(int nprocs);
 #endif
+
+   // Load balancing support
+   bool enable_load_balance;         // Enable load balancing
+   int load_balance_check_interval;  // Check interval (in time steps)
+   int current_time_step;            // Current time step counter
+   double *rank_interp_times;        // Store interpolation times for each rank
+   int *heavy_ranks;                 // Store heavy rank numbers
+   int num_heavy_ranks;              // Number of heavy ranks
+
+   void init_load_balance(int nprocs);
+   void update_interp_time(int rank, double time);
+   bool check_and_rebalance(int nprocs, int lev,
+                           MyList<var> *OldList, MyList<var> *StateList,
+                           MyList<var> *FutureList, MyList<var> *tmList,
+                           int Symmetry, bool BB);
 };
 
 #endif /* CGH_H */

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/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/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/makefile.inc

@@ -13,7 +13,7 @@ LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore
 ## 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
+PROFDATA     = ../../pgo_profile/default.profdata
 CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma -ipo \
                -fprofile-instr-use=$(PROFDATA) \
                -Dfortran3 -Dnewc -I${MKLROOT}/include

AMSS_NCKU_source/surface_integral.C

@@ -11,6 +11,8 @@
 #include <strstream>
 #include <cmath>
 #include <map>
+#include <vector>
+#include <algorithm>
 using namespace std;
 #else
 #include <iostream.h>
@@ -238,9 +240,6 @@ void surface_integral::surf_Wave(double rex, int lev, cgh *GH, var *Rpsi4, var *
   shellf = new double[n_tot * InList];
 
   GH->PatL[lev]->data->Interp_Points(DG_List, n_tot, pox, shellf, Symmetry, Nmin, Nmax);
-
-  //|~~~~~> Integrate the dot product of Dphi with the surface normal.
-
   double *RP_out, *IP_out;
   RP_out = new double[NN];
   IP_out = new double[NN];

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 " )