Add thread-safe ShellPatch::setupintintstuff with OpenMP

Split prolongpointstru into search-only (prolongpointstru_search) and append-only (prolongpointstru_append) functions. The search is read-only and thread-safe; each thread builds private linked lists via prolongpointstru_append, merged after the parallel loop. This eliminates critical-section contention and delivers ~2.2x speedup: setupintintstuff: 511s -> 252s, total init: 592s -> 267s. Also add -qopenmp to ShellPatch.o compilation via makefile override rule and <omp.h> include with _OPENMP guards + fallback stubs. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
2026-05-10 02:10:20 +08:00
parent 26f5add3fb
commit fae4093863
3 changed files with 3738 additions and 3489 deletions
--- a/AMSS_NCKU_source/ShellPatch.C
+++ b/AMSS_NCKU_source/ShellPatch.C
@@ -8,6 +8,12 @@
 #include <cmath>
 #include <new>
 #include <map>
+#ifdef _OPENMP
+#include <omp.h>
+#else
+static inline int omp_get_max_threads() { return 1; }
+static inline int omp_get_thread_num()  { return 0; }
+#endif
 #include <vector>
 #include <thread>
 #include <algorithm>
@@ -2546,6 +2552,176 @@ void ShellPatch::prolongpointstru(MyList<pointstru> *&psul, MyList<ss_patch> *sP
  }
 }
 // used by _src construction, so these x,y,z do not coinside with grid point
+ShellPatch::PointSearchResult ShellPatch::prolongpointstru_search(bool ssyn, int tsst, MyList<ss_patch> *sPp,
+    double DH[dim], MyList<Patch> *Pp, double CDH[dim], double x, double y, double z, int Symmetry, int rank_in)
+{
+  PointSearchResult sr;
+  sr.found = false;
+  sr.Bg = nullptr;
+  sr.gx = x; sr.gy = y; sr.gz = z;
+  sr.lx = sr.ly = sr.lz = 0.0;
+  sr.ssst = -1;
+
+  MyList<Block> *Bgl;
+  Block *Bg;
+  double llb[dim], uub[dim];
+
+  if (ssyn)
+  {
+    int sst;
+    double lx, ly, lz;
+    getlocalpox(x, y, z, sst, lx, ly, lz);
+    while (sPp)
+    {
+      if (sPp->data->sst == sst)
+      {
+        Bgl = sPp->data->blb;
+        while (Bgl)
+        {
+          Bg = Bgl->data;
+          if (Bg->rank == rank_in)
+          {
+            for (int j = 0; j < 2; j++)
+            {
+              if (feq(Bg->bbox[j], -PI / 4 - overghost * DH[j], DH[j] / 2))
+                llb[j] = -PI / 4;
+              else if (feq(Bg->bbox[j], sPp->data->bbox[j], DH[j] / 2))
+                llb[j] = Bg->bbox[j];
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+              else
+                llb[j] = Bg->bbox[j] + (ghost_width - 1) * DH[j];
+#else
+#ifdef Cell
+              else
+                llb[j] = Bg->bbox[j] + ghost_width * DH[j];
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+              if (feq(Bg->bbox[dim + j], PI / 4 + overghost * DH[j], DH[j] / 2))
+                uub[j] = PI / 4;
+              else if (feq(Bg->bbox[dim + j], sPp->data->bbox[dim + j], DH[j] / 2))
+                uub[j] = Bg->bbox[dim + j];
+              else
+                uub[j] = Bg->bbox[dim + j] - ghost_width * DH[j];
+            }
+            if (feq(Bg->bbox[2], sPp->data->bbox[2], DH[2] / 2))
+              llb[2] = Bg->bbox[2];
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+            else
+              llb[2] = Bg->bbox[2] + (ghost_width - 1) * DH[2];
+#else
+#ifdef Cell
+            else
+              llb[2] = Bg->bbox[2] + ghost_width * DH[2];
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+            if (feq(Bg->bbox[dim + 2], sPp->data->bbox[dim + 2], DH[2] / 2))
+              uub[2] = Bg->bbox[dim + 2];
+            else
+              uub[2] = Bg->bbox[dim + 2] - ghost_width * DH[2];
+            if (lx > llb[0] - 0.0001 * DH[0] && lx < uub[0] + 0.0001 * DH[0] &&
+                ly > llb[1] - 0.0001 * DH[1] && ly < uub[1] + 0.0001 * DH[1] &&
+                lz > llb[2] - 0.0001 * DH[2] && lz < uub[2] + 0.0001 * DH[2])
+            {
+              sr.found = true;
+              sr.Bg = Bg;
+              sr.lx = lx; sr.ly = ly; sr.lz = lz;
+              sr.ssst = sPp->data->sst;
+              return sr;
+            }
+          }
+          if (Bgl == sPp->data->ble) break;
+          Bgl = Bgl->next;
+        }
+      }
+      sPp = sPp->next;
+    }
+  }
+  else
+  {
+    while (Pp)
+    {
+      Bgl = Pp->data->blb;
+      while (Bgl)
+      {
+        Bg = Bgl->data;
+        if (Bg->rank == rank_in)
+        {
+          for (int j = 0; j < dim; j++)
+          {
+            if (feq(Bg->bbox[j], Pp->data->bbox[j], CDH[j] / 2))
+              llb[j] = Bg->bbox[j];
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+            else
+              llb[j] = Bg->bbox[j] + (ghost_width - 1) * CDH[j];
+#else
+#ifdef Cell
+            else
+              llb[j] = Bg->bbox[j] + ghost_width * CDH[j];
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+            if (feq(Bg->bbox[dim + j], Pp->data->bbox[dim + j], CDH[j] / 2))
+              uub[j] = Bg->bbox[dim + j];
+            else
+              uub[j] = Bg->bbox[dim + j] - ghost_width * CDH[j];
+          }
+          if (x > llb[0] - 0.0001 * CDH[0] && x < uub[0] + 0.0001 * CDH[0] &&
+              y > llb[1] - 0.0001 * CDH[1] && y < uub[1] + 0.0001 * CDH[1] &&
+              z > llb[2] - 0.0001 * CDH[2] && z < uub[2] + 0.0001 * CDH[2])
+          {
+            sr.found = true;
+            sr.Bg = Bg;
+            sr.lx = x; sr.ly = y; sr.lz = z;
+            sr.ssst = -1;
+            return sr;
+          }
+        }
+        if (Bgl == Pp->data->ble) break;
+        Bgl = Bgl->next;
+      }
+      Pp = Pp->next;
+    }
+  }
+  return sr;
+}
+
+void ShellPatch::prolongpointstru_append(MyList<pointstru> *&psul, const PointSearchResult &sr, int tsst)
+{
+  MyList<pointstru> *ps = new MyList<pointstru>;
+  ps->data = new pointstru;
+  ps->data->Bg = sr.Bg;
+  ps->data->gpox[0] = sr.gx;
+  ps->data->gpox[1] = sr.gy;
+  ps->data->gpox[2] = sr.gz;
+  ps->data->lpox[0] = sr.lx;
+  ps->data->lpox[1] = sr.ly;
+  ps->data->lpox[2] = sr.lz;
+  ps->data->ssst = sr.ssst;
+  ps->data->tsst = tsst;
+  ps->data->dumyd = getdumydimension(ps->data->ssst, ps->data->tsst);
+  ps->data->coef = 0;
+  ps->data->sind = 0;
+  ps->next = 0;
+  if (psul)
+    psul->catList(ps);
+  else
+    psul = ps;
+}
+
 bool ShellPatch::prolongpointstru(MyList<pointstru> *&psul, bool ssyn, int tsst, MyList<ss_patch> *sPp, double DH[dim],
                                  MyList<Patch> *Pp, double CDH[dim], double x, double y, double z, int Symmetry, int rank_in)
 {
@@ -2763,12 +2939,26 @@ void ShellPatch::setupintintstuff(int cpusize, MyList<Patch> *CPatL, int Symmetr
  }

  sPp = PatL;
+  // Thread-local lists to avoid critical-section contention
+  int nt = omp_get_max_threads();
+  MyList<pointstru> ***tl_csats = new MyList<pointstru>**[nt];
+  MyList<pointstru> ***tl_ss    = new MyList<pointstru>**[nt];
+  for (int t = 0; t < nt; t++) {
+    tl_csats[t] = new MyList<pointstru>*[cpusize];
+    tl_ss[t]    = new MyList<pointstru>*[cpusize];
+    for (int i = 0; i < cpusize; i++) {
+      tl_csats[t][i] = nullptr;
+      tl_ss[t][i]    = nullptr;
+    }
+  }
  while (sPp)
  {
+#pragma omp parallel for collapse(3) schedule(guided)
    for (int iz = 0; iz < sPp->data->shape[2]; iz++)
      for (int is = 0; is < sPp->data->shape[1]; is++)
        for (int ir = 0; ir < sPp->data->shape[0]; ir++)
        {
+          double x, y, z;
 #ifdef Vertex
 #ifdef Cell
 #error Both Cell and Vertex are defined
@@ -2785,44 +2975,46 @@ void ShellPatch::setupintintstuff(int cpusize, MyList<Patch> *CPatL, int Symmetr
 #error Not define Vertex nor Cell
 #endif
 #endif
+          int tid = omp_get_thread_num();
          if (z < sPp->data->bbox[2] + (SC_width + 0.0001) * DH[2])
          {
            double gx, gy, gz;
            getglobalpox(gx, gy, gz, sPp->data->sst, x, y, z);
-            bool flag = false;
+            PointSearchResult sr;
+            int found_rank = -1;
            for (int i = 0; i < cpusize; i++)
            {
-              flag = prolongpointstru(csats_src[i], false, sPp->data->sst, PatL, DH, CPatL, CDH, gx, gy, gz, Symmetry, i);
-              if (flag)
-                break;
+              sr = prolongpointstru_search(false, sPp->data->sst, PatL, DH, CPatL, CDH, gx, gy, gz, Symmetry, i);
+              if (sr.found) { found_rank = i; break; }
            }
-            if (!flag)
+            if (found_rank >= 0)
+              prolongpointstru_append(tl_csats[tid][found_rank], sr, sPp->data->sst);
+            else
            {
-              CPatL->data->checkBlock();
              if (myrank == 0)
              {
                cout << "ShellPatch::prolongpointstru fail to find cardisian source point for" << endl;
                cout << "sst = " << sPp->data->sst << " lx,ly,lz = " << x << "," << y << "," << z << endl;
                cout << "x,y,z = " << gx << "," << gy << "," << gz << endl;
-                MPI_Abort(MPI_COMM_WORLD, 1);
              }
+              MPI_Abort(MPI_COMM_WORLD, 1);
            }
          }
-          //	       else if(x<-PI/4-(overghost-ghost_width-0.0001)*DH[0] || x>PI/4+(overghost-ghost_width-0.0001)*DH[0] ||
-          //	               y<-PI/4-(overghost-ghost_width-0.0001)*DH[1] || y>PI/4+(overghost-ghost_width-0.0001)*DH[1]   ) //0.0001 is for vertex center
          if (x < -PI / 4 - (overghost - ghost_width - 0.0001) * DH[0] || x > PI / 4 + (overghost - ghost_width - 0.0001) * DH[0] ||
              y < -PI / 4 - (overghost - ghost_width - 0.0001) * DH[1] || y > PI / 4 + (overghost - ghost_width - 0.0001) * DH[1])
          {
            double gx, gy, gz;
            getglobalpox(gx, gy, gz, sPp->data->sst, x, y, z);
-            bool flag = false;
+            PointSearchResult sr;
+            int found_rank = -1;
            for (int i = 0; i < cpusize; i++)
            {
-              flag = prolongpointstru(ss_src[i], true, sPp->data->sst, PatL, DH, CPatL, CDH, gx, gy, gz, Symmetry, i);
-              if (flag)
-                break;
+              sr = prolongpointstru_search(true, sPp->data->sst, PatL, DH, CPatL, CDH, gx, gy, gz, Symmetry, i);
+              if (sr.found) { found_rank = i; break; }
            }
-            if (!flag)
+            if (found_rank >= 0)
+              prolongpointstru_append(tl_ss[tid][found_rank], sr, sPp->data->sst);
+            else
            {
              if (myrank == 0)
              {
@@ -2830,17 +3022,41 @@ void ShellPatch::setupintintstuff(int cpusize, MyList<Patch> *CPatL, int Symmetr
                cout << "sst = " << sPp->data->sst << " lx,ly,lz = " << x << "," << y << "," << z << endl;
                if (sPp->data->sst == -1)
                  cout << "your angular resolution for shell is too coarse?" << endl;
-                MPI_Abort(MPI_COMM_WORLD, 1);
              }
+              MPI_Abort(MPI_COMM_WORLD, 1);
            }
          }
        }
    sPp = sPp->next;
  }
+  // Merge thread-local lists into main lists
+  for (int t = 0; t < nt; t++) {
+    for (int i = 0; i < cpusize; i++) {
+      if (tl_csats[t][i]) {
+        if (csats_src[i]) csats_src[i]->catList(tl_csats[t][i]);
+        else csats_src[i] = tl_csats[t][i];
+      }
+      if (tl_ss[t][i]) {
+        if (ss_src[i]) ss_src[i]->catList(tl_ss[t][i]);
+        else ss_src[i] = tl_ss[t][i];
+      }
+    }
+    delete[] tl_csats[t];
+    delete[] tl_ss[t];
+  }
+  delete[] tl_csats;
+  delete[] tl_ss;
  if (myrank == 0)
    cout << " ShellPatch::setup interpatch interpolation stuffs ss_src completes" << endl;

  Pp = CPatL;
+  // Thread-local lists for csatc_src
+  nt = omp_get_max_threads();
+  MyList<pointstru> ***tl_csatc = new MyList<pointstru>**[nt];
+  for (int t = 0; t < nt; t++) {
+    tl_csatc[t] = new MyList<pointstru>*[cpusize];
+    for (int i = 0; i < cpusize; i++) tl_csatc[t][i] = nullptr;
+  }
  while (Pp)
  {
    double llb[dim], uub[dim];
@@ -2857,10 +3073,12 @@ void ShellPatch::setupintintstuff(int cpusize, MyList<Patch> *CPatL, int Symmetr
        llb[j] = Pp->data->bbox[j] + (CS_width + 0.0001) * CDH[j];
      uub[j] = Pp->data->bbox[dim + j] - (CS_width + 0.0001) * CDH[j];
    }
+#pragma omp parallel for collapse(3) schedule(guided)
    for (int iz = 0; iz < Pp->data->shape[2]; iz++)
      for (int iy = 0; iy < Pp->data->shape[1]; iy++)
        for (int ix = 0; ix < Pp->data->shape[0]; ix++)
        {
+          double x, y, z;
 #ifdef Vertex
 #ifdef Cell
 #error Both Cell and Vertex are defined
@@ -2881,29 +3099,40 @@ void ShellPatch::setupintintstuff(int cpusize, MyList<Patch> *CPatL, int Symmetr
              y < llb[1] || y > uub[1] ||
              z < llb[2] || z > uub[2])
          {
-            int sst;
-            double lx, ly, lz;
-            bool flag = false;
-            getlocalpox(x, y, z, sst, lx, ly, lz);
+            int tid = omp_get_thread_num();
+            PointSearchResult sr;
+            int found_rank = -1;
            for (int i = 0; i < cpusize; i++)
            {
-              flag = prolongpointstru(csatc_src[i], true, -1, PatL, DH, CPatL, CDH, x, y, z, Symmetry, i);
-              if (flag)
-                break;
+              sr = prolongpointstru_search(true, -1, PatL, DH, CPatL, CDH, x, y, z, Symmetry, i);
+              if (sr.found) { found_rank = i; break; }
            }
-            if (!flag)
+            if (found_rank >= 0)
+              prolongpointstru_append(tl_csatc[tid][found_rank], sr, -1);
+            else
            {
              if (myrank == 0)
              {
                cout << "ShellPatch::prolongpointstru fail to find shell source point for" << endl;
                cout << "sst = -1, x,y,z = " << x << "," << y << "," << z << endl;
-                MPI_Abort(MPI_COMM_WORLD, 1);
              }
+              MPI_Abort(MPI_COMM_WORLD, 1);
            }
          }
        }
    Pp = Pp->next;
  }
+  // Merge thread-local csatc_src lists
+  for (int t = 0; t < nt; t++) {
+    for (int i = 0; i < cpusize; i++) {
+      if (tl_csatc[t][i]) {
+        if (csatc_src[i]) csatc_src[i]->catList(tl_csatc[t][i]);
+        else csatc_src[i] = tl_csatc[t][i];
+      }
+    }
+    delete[] tl_csatc[t];
+  }
+  delete[] tl_csatc;
  if (myrank == 0)
    cout << " ShellPatch::setup interpatch interpolation stuffs csatc_src and csats_src completes" << endl;

--- a/AMSS_NCKU_source/ShellPatch.h
+++ b/AMSS_NCKU_source/ShellPatch.h
@@ -102,6 +102,16 @@ public:
                 //-1: means no dumy dimension at all; 0: means rho; 1: means sigma
   };

+   // Thread-safe search result (no pointers to shared mutable state)
+   struct PointSearchResult
+   {
+      bool found;
+      Block *Bg;
+      double gx, gy, gz; // global Cartesian coordinates
+      double lx, ly, lz; // local coordinates within the found block
+      int ssst;          // source shell-patch type (-1 = Cartesian)
+   };
+
   int myrank;
   int shape[dim];   // for (rho, sigma, R), for rho and sigma means number of points for every pi/2
   double Rrange[2]; // for Rmin and Rmax
@@ -175,6 +185,12 @@ public:
                         MyList<Patch> *Pp, double CDH[dim], MyList<pointstru> *pss);
   bool prolongpointstru(MyList<pointstru> *&psul, bool ssyn, int tsst, MyList<ss_patch> *sPp, double DH[dim],
                         MyList<Patch> *Pp, double CDH[dim], double x, double y, double z, int Symmetry, int rank_in);
+   // Read-only point search — thread-safe (no shared mutable state modified)
+   PointSearchResult prolongpointstru_search(bool ssyn, int tsst, MyList<ss_patch> *sPp, double DH[dim],
+                                             MyList<Patch> *Pp, double CDH[dim], double x, double y, double z,
+                                             int Symmetry, int rank_in);
+   // Append a search result to a linked list — use inside omp critical section
+   void prolongpointstru_append(MyList<pointstru> *&psul, const PointSearchResult &sr, int tsst);
   void setupintintstuff(int cpusize, MyList<Patch> *CPatL, int Symmetry);
   void intertransfer(MyList<pointstru> **src, MyList<pointstru> **dst,
                      MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /*target */,
--- a/AMSS_NCKU_source/makefile
+++ b/AMSS_NCKU_source/makefile
@@ -73,6 +73,10 @@ endif
 .C.o:
 	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@

+# ShellPatch.C uses OpenMP for setupintintstuff search loops
+ShellPatch.o: ShellPatch.C
+	${CXX} $(CXXAPPFLAGS) $(OMP_FLAG) -c $< $(filein) -o $@
+
 .for.o:
 	$(f77) -c $< -o $@