补全C算子和Fortran算子的数学差异

AMSS_NCKU_source/bssn_rhs_c.C

@@ -70,10 +70,34 @@ int f_compute_rhs_bssn(int *ex, double &T,
     const double FF = 0.75, eta = 2.0;
     const double F1o3 = 1.0/3.0, F2o3 = 2.0/3.0, F3o2 = 1.5, F1o6 = 1.0/6.0;
     const double F16 = 16.0, F8 = 8.0;
-    #if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5)
+    #if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7)
     double reta[all];
     /* 使用时：reta[idx]，其中 idx = i + nx*(j + ny*k)  (Fortran列主序) */
     #endif
+
+    #if (GAUGE == 6 || GAUGE == 7)
+    int BHN = 0;
+    double Porg[9] = {0.0};
+    double Mass[3] = {0.0};
+
+    extern "C" {
+        #ifdef fortran1
+        void getpbh(int &, double *, double *);
+        #elif defined(fortran2)
+        void GETPBH(int &, double *, double *);
+        #else
+        void getpbh_(int &, double *, double *);
+        #endif
+    }
+
+    #ifdef fortran1
+    getpbh(BHN, Porg, Mass);
+    #elif defined(fortran2)
+    GETPBH(BHN, Porg, Mass);
+    #else
+    getpbh_(BHN, Porg, Mass);
+    #endif
+    #endif
     PI = acos(-1.0); 
     dX = X[1] - X[0];
     dY = Y[1] - Y[0];
@@ -966,11 +990,67 @@ int f_compute_rhs_bssn(int *ex, double &T,
         }
 
         // 1ms //
+        #if (GAUGE == 6 || GAUGE == 7)
+        if (BHN == 2) {
+            const double M = Mass[0] + Mass[1];
+            const double A = TWO / M;
+            const double w1 = 1.2e1;
+            const double w2 = w1;
+            const double C1 = ONE / Mass[0] - A;
+            const double C2 = ONE / Mass[1] - A;
+            const double denom =
+                (Porg[0] - Porg[3]) * (Porg[0] - Porg[3]) +
+                (Porg[1] - Porg[4]) * (Porg[1] - Porg[4]) +
+                (Porg[2] - Porg[5]) * (Porg[2] - Porg[5]);
+
+            for (int k0 = 0; k0 < nz; ++k0) {
+                for (int j0 = 0; j0 < ny; ++j0) {
+                    for (int i0 = 0; i0 < nx; ++i0) {
+                        const size_t p = idx_ex(i0, j0, k0, ex);
+                        const double dx1 = Porg[0] - X[i0];
+                        const double dy1 = Porg[1] - Y[j0];
+                        const double dz1 = Porg[2] - Z[k0];
+                        const double dx2 = Porg[3] - X[i0];
+                        const double dy2 = Porg[4] - Y[j0];
+                        const double dz2 = Porg[5] - Z[k0];
+                        const double r1 = (dx1 * dx1 + dy1 * dy1 + dz1 * dz1) / denom;
+                        const double r2 = (dx2 * dx2 + dy2 * dy2 + dz2 * dz2) / denom;
+                        #if (GAUGE == 6)
+                        reta[p] = A + C1 / (ONE + w1 * r1) + C2 / (ONE + w2 * r2);
+                        #else
+                        reta[p] = A + C1 * exp(-w1 * r1) + C2 * exp(-w2 * r2);
+                        #endif
+                    }
+                }
+            }
+        } else {
+            printf("not support BH_num in Jason's form %d %d\n", (GAUGE == 6) ? 1 : 2, BHN);
+            for (int i = 0; i < all; ++i) reta[i] = ZEO;
+        }
+        #endif
+
         for (int i = 0; i < all; i += 1) {
+            #if (GAUGE == 0)
             betax_rhs[i] = FF * dtSfx[i];
             betay_rhs[i] = FF * dtSfy[i];
             betaz_rhs[i] = FF * dtSfz[i];
-#if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5)
+
+            dtSfx_rhs[i] = Gamx_rhs[i] - eta * dtSfx[i];
+            dtSfy_rhs[i] = Gamy_rhs[i] - eta * dtSfy[i];
+            dtSfz_rhs[i] = Gamz_rhs[i] - eta * dtSfz[i];
+            #elif (GAUGE == 1)
+            betax_rhs[i] = Gamx[i] - eta * betax[i];
+            betay_rhs[i] = Gamy[i] - eta * betay[i];
+            betaz_rhs[i] = Gamz[i] - eta * betaz[i];
+
+            dtSfx_rhs[i] = ZEO;
+            dtSfy_rhs[i] = ZEO;
+            dtSfz_rhs[i] = ZEO;
+            #elif (GAUGE == 2 || GAUGE == 3)
+            betax_rhs[i] = FF * dtSfx[i];
+            betay_rhs[i] = FF * dtSfy[i];
+            betaz_rhs[i] = FF * dtSfz[i];
+
             reta[i] =
                 gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
                 + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
@@ -979,15 +1059,45 @@ int f_compute_rhs_bssn(int *ex, double &T,
                         + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
                         + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
 
-            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (1.0 - sqrt(chin1[i])), 2.0 );
+            #if (GAUGE == 2)
+            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
+            #else
+            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - chin1[i]), 2.0 );
+            #endif
 
             dtSfx_rhs[i] = Gamx_rhs[i] - reta[i] * dtSfx[i];
             dtSfy_rhs[i] = Gamy_rhs[i] - reta[i] * dtSfy[i];
             dtSfz_rhs[i] = Gamz_rhs[i] - reta[i] * dtSfz[i];
+            #elif (GAUGE == 4 || GAUGE == 5)
+            reta[i] =
+                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
+                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
+                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
+                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
+                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
+                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
+
+            #if (GAUGE == 4)
+            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
             #else
-            dtSfx_rhs[i] = Gamx_rhs[i] - eta * dtSfx[i];
+            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - chin1[i]), 2.0 );
-            dtSfy_rhs[i] = Gamy_rhs[i] - eta * dtSfy[i];
+            #endif
-            dtSfz_rhs[i] = Gamz_rhs[i] - eta * dtSfz[i];
+
+            betax_rhs[i] = FF * Gamx[i] - reta[i] * betax[i];
+            betay_rhs[i] = FF * Gamy[i] - reta[i] * betay[i];
+            betaz_rhs[i] = FF * Gamz[i] - reta[i] * betaz[i];
+
+            dtSfx_rhs[i] = ZEO;
+            dtSfy_rhs[i] = ZEO;
+            dtSfz_rhs[i] = ZEO;
+            #elif (GAUGE == 6 || GAUGE == 7)
+            betax_rhs[i] = FF * dtSfx[i];
+            betay_rhs[i] = FF * dtSfy[i];
+            betaz_rhs[i] = FF * dtSfz[i];
+
+            dtSfx_rhs[i] = Gamx_rhs[i] - reta[i] * dtSfx[i];
+            dtSfy_rhs[i] = Gamy_rhs[i] - reta[i] * dtSfy[i];
+            dtSfz_rhs[i] = Gamz_rhs[i] - reta[i] * dtSfz[i];
             #endif
         }
         // 26ms //