From e06a5b337476d4f772e4e90892118734146f281a Mon Sep 17 00:00:00 2001
From: Liwei Ji <jiliwei.phys@gmail.com>
Date: Fri, 28 Jun 2024 04:27:22 +0000
Subject: [PATCH] Multipole: only integrate on root rank

---
 Multipole/src/multipole.cxx |  3 +-
 Multipole/src/surface.cxx   | 73 ++++++++++++++++++++-----------------
 Multipole/src/surface.hxx   |  2 +-
 3 files changed, 42 insertions(+), 36 deletions(-)

diff --git a/Multipole/src/multipole.cxx b/Multipole/src/multipole.cxx
index e0cf602c..f4e146dd 100644
--- a/Multipole/src/multipole.cxx
+++ b/Multipole/src/multipole.cxx
@@ -155,7 +155,8 @@ extern "C" void Multipole_Calc(CCTK_ARGUMENTS) {
       // Intergate of conj(sYlm)*F*sin(theta) over the sphere at radiusr[i]
       for (int l = 0; l <= l_max; l++) {
         for (int m = -l; m <= l; m++) {
-          g_sphere->integrate(g_sphere->getRealY()[si][l][m + l],
+          g_sphere->integrate(CCTK_PASS_CTOC,
+                              g_sphere->getRealY()[si][l][m + l],
                               g_sphere->getImagY()[si][l][m + l],
                               g_sphere->getRealF(), g_sphere->getImagF(),
                               &modes(v, i, l, m, 0), &modes(v, i, l, m, 1));
diff --git a/Multipole/src/surface.cxx b/Multipole/src/surface.cxx
index d35a0220..413da778 100644
--- a/Multipole/src/surface.cxx
+++ b/Multipole/src/surface.cxx
@@ -66,50 +66,55 @@ void Surface::interpolate(CCTK_ARGUMENTS, int realFieldIndex,
 }
 
 // Take the integral of conj(array1)*array2*sin(th)
-void Surface::integrate(const std::vector<CCTK_REAL> &array1r,
+void Surface::integrate(CCTK_ARGUMENTS, const std::vector<CCTK_REAL> &array1r,
                         const std::vector<CCTK_REAL> &array1i,
                         const std::vector<CCTK_REAL> &array2r,
                         const std::vector<CCTK_REAL> &array2i, CCTK_REAL *outRe,
                         CCTK_REAL *outIm) {
   DECLARE_CCTK_PARAMETERS;
 
-  std::vector<CCTK_REAL> fReal(theta_.size());
-  std::vector<CCTK_REAL> fImag(theta_.size());
+  if (CCTK_MyProc(cctkGH) == 0) {
 
-  // integrand: conj(array1)*array2*sin(th)
-  for (size_t i = 0; i < theta_.size(); ++i) {
-    fReal[i] = (array1r[i] * array2r[i] + array1i[i] * array2i[i]) *
-               std::sin(theta_[i]);
-    fImag[i] = (array1r[i] * array2i[i] - array1i[i] * array2r[i]) *
-               std::sin(theta_[i]);
-  }
+    std::vector<CCTK_REAL> fReal(theta_.size());
+    std::vector<CCTK_REAL> fImag(theta_.size());
 
-  if (CCTK_Equals(integration_method, "midpoint")) {
-    *outRe = Midpoint2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-    *outIm = Midpoint2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-  } else if (CCTK_Equals(integration_method, "trapezoidal")) {
-    *outRe =
-        Trapezoidal2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-    *outIm =
-        Trapezoidal2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-  } else if (CCTK_Equals(integration_method, "Simpson")) {
-    if (nPhi_ % 2 != 0 || nTheta_ % 2 != 0) {
-      CCTK_WARN(CCTK_WARN_ABORT, "The Simpson integration method requires even "
-                                 "nTheta_ and even nPhi_");
+    // integrand: conj(array1)*array2*sin(th)
+    for (size_t i = 0; i < theta_.size(); ++i) {
+      fReal[i] = (array1r[i] * array2r[i] + array1i[i] * array2i[i]) *
+                 std::sin(theta_[i]);
+      fImag[i] = (array1r[i] * array2i[i] - array1i[i] * array2r[i]) *
+                 std::sin(theta_[i]);
     }
-    *outRe = Simpson2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-    *outIm = Simpson2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-  } else if (CCTK_Equals(integration_method, "DriscollHealy")) {
-    if (nTheta_ % 2 != 0) {
-      CCTK_WARN(CCTK_WARN_ABORT,
-                "The Driscoll&Healy integration method requires even nTheta_");
+
+    if (CCTK_Equals(integration_method, "midpoint")) {
+      *outRe = Midpoint2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+      *outIm = Midpoint2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+    } else if (CCTK_Equals(integration_method, "trapezoidal")) {
+      *outRe =
+          Trapezoidal2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+      *outIm =
+          Trapezoidal2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+    } else if (CCTK_Equals(integration_method, "Simpson")) {
+      if (nPhi_ % 2 != 0 || nTheta_ % 2 != 0) {
+        CCTK_WARN(CCTK_WARN_ABORT,
+                  "The Simpson integration method requires even "
+                  "nTheta_ and even nPhi_");
+      }
+      *outRe = Simpson2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+      *outIm = Simpson2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+    } else if (CCTK_Equals(integration_method, "DriscollHealy")) {
+      if (nTheta_ % 2 != 0) {
+        CCTK_WARN(
+            CCTK_WARN_ABORT,
+            "The Driscoll&Healy integration method requires even nTheta_");
+      }
+      *outRe =
+          DriscollHealy2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+      *outIm =
+          DriscollHealy2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
+    } else {
+      CCTK_WARN(CCTK_WARN_ABORT, "internal error");
     }
-    *outRe =
-        DriscollHealy2DIntegral(fReal.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-    *outIm =
-        DriscollHealy2DIntegral(fImag.data(), nTheta_, nPhi_, dTheta_, dPhi_);
-  } else {
-    CCTK_WARN(CCTK_WARN_ABORT, "internal error");
   }
 }
 
diff --git a/Multipole/src/surface.hxx b/Multipole/src/surface.hxx
index 45bd0139..c95cda84 100644
--- a/Multipole/src/surface.hxx
+++ b/Multipole/src/surface.hxx
@@ -60,7 +60,7 @@ public:
   void interpolate(CCTK_ARGUMENTS, int realFieldIndex, int imagFieldIndex);
 
   // Take the integral of conj(array1)*array2*sin(th)
-  void integrate(const std::vector<CCTK_REAL> &array1r,
+  void integrate(CCTK_ARGUMENTS, const std::vector<CCTK_REAL> &array1r,
                  const std::vector<CCTK_REAL> &array1i,
                  const std::vector<CCTK_REAL> &array2r,
                  const std::vector<CCTK_REAL> &array2i, CCTK_REAL *outre,