tmpwip

3 files changed, 87 insertions, 91 deletions

diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index 40b3503..7955be7 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -116,7 +116,6 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi)
     int y_idx;
     int gi, var0;
     int n_vars;
-    int timelevel;
 
     gi     = CCTK_GroupIndexFromVarI (vi);
     n_vars = CCTK_NumVarsInGroupI (gi);
@@ -160,8 +159,6 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi)
 
     s = cctkGH->cctk_nghostzones[0];
 
-    timelevel = 0;
-
     /* find the index for which y == 0 */
     for (i = 0; i < cctkGH->cctk_lsh[1]; i++) {
         if (fabs(y[CCTK_GFINDEX3D(cctkGH, 0, i, 0)]) < 1e-12) {
@@ -173,100 +170,99 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi)
         CCTK_WARN(0, "y is not 0 anywhere");
 
     /* make sure that the input data is synced */
-    CCTK_SyncGroupI(cctkGH, gi);
+    //CCTK_SyncGroupI(cctkGH, gi);
 
 #pragma omp parallel for
   for (int k = 0; k < cctkGH->cctk_lsh[2]; k++)
-      for (int jj = 0; jj < cctkGH->cctk_lsh[1]; jj++)
-          for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) {
-              CCTK_REAL f[100];
-              CCTK_REAL rho_val;
-              int idx_dst, ijk;
-
-              if (jj != y_idx || i >= s)
-                  continue;
-              //if (jj == y_idx && i >= s)
-              //    continue;
-
-              /* index into linear memory */
-              ijk     = CCTK_GFINDEX3D(cctkGH, i, y_idx, k);
-              idx_dst = CCTK_GFINDEX3D(cctkGH, i, jj, k);
-
-              rho_val = rho[idx_dst];
+      for (int i = 0; i < s; i++) {
+          CCTK_REAL f[100];
+          CCTK_REAL rho_val;
+          int idx_dst, idx_src;
+
+          //if (jj == y_idx && i >= s)
+          //    continue;
+
+          /* index into linear memory */
+          idx_src = CCTK_GFINDEX3D(cctkGH, 2 * s -i, y_idx, k);
+          idx_dst = CCTK_GFINDEX3D(cctkGH, i, y_idx, k);
+
+          rho_val = rho[idx_dst];
+
+          /* interpolate grid functions */
+          for (int n = 0; n < n_vars; n++) {
+              double *src = cctkGH->data[var0 + n][0];
+              //f[n] = cartoon2d_interp(cctkGH, order, cctkGH->data[var0+n][0],
+              //                        i, ijk, rho_val);
+              f[n] = src[idx_src];
+          }
 
-              /* interpolate grid functions */
+          /* rotate grid tensor by matrix multiplication */
+          if (tensortype == TENSORTYPE_SCALAR) {
+              /* Scalar groups can have arbitrary many components; these
+                 "components" are then all scalars.  */
               for (int n = 0; n < n_vars; n++) {
-                  f[n] = cartoon2d_interp(cctkGH, order, cctkGH->data[var0+n][timelevel],
-                                          i, ijk, rho_val);
-              }
-
-              /* rotate grid tensor by matrix multiplication */
-              if (tensortype == TENSORTYPE_SCALAR) {
-                  /* Scalar groups can have arbitrary many components; these
-                     "components" are then all scalars.  */
-                  for (int n = 0; n < n_vars; n++) {
-                      CCTK_REAL *t = cctkGH->data[var0 + n][timelevel];
-                      t[idx_dst] = f[n];
-                  }
-              } else if (tensortype == TENSORTYPE_U) {
-                  CCTK_REAL *tx = cctkGH->data[vi][timelevel];
-                  CCTK_REAL *ty = cctkGH->data[vi+1][timelevel];
-                  CCTK_REAL *tz = cctkGH->data[vi+2][timelevel];
-
-                  CCTK_REAL rxx_val = rxx[idx_dst];
-                  CCTK_REAL ryy_val = ryy[idx_dst];
-                  CCTK_REAL rxy_val = rxy[idx_dst];
-                  CCTK_REAL ryx_val = ryx[idx_dst];
-
-                  CCTK_REAL fx = f[0];
-                  CCTK_REAL fy = f[1];
-                  CCTK_REAL fz = f[2];
-
-                  tx[idx_dst] = rxx_val * fx + rxy_val * fy;
-                  ty[idx_dst] = ryx_val * fx + ryy_val * fy;
-                  tz[idx_dst] = fz;
-
-                  //tx[ijk-dj] = sxx * fx + sxy * fy;
-                  //ty[ijk-dj] = syx * fx + syy * fy;
-                  //tz[ijk-dj] = fz;
-              } else if (tensortype == TENSORTYPE_DDSYM) {
-                  CCTK_REAL *txx = cctkGH->data[vi][timelevel];
-                  CCTK_REAL *txy = cctkGH->data[vi+1][timelevel];
-                  CCTK_REAL *txz = cctkGH->data[vi+2][timelevel];
-                  CCTK_REAL *tyy = cctkGH->data[vi+3][timelevel];
-                  CCTK_REAL *tyz = cctkGH->data[vi+4][timelevel];
-                  CCTK_REAL *tzz = cctkGH->data[vi+5][timelevel];
-
-                  CCTK_REAL sxx_val = sxx[idx_dst];
-                  CCTK_REAL syy_val = syy[idx_dst];
-                  CCTK_REAL sxy_val = sxy[idx_dst];
-                  CCTK_REAL syx_val = syx[idx_dst];
-
-                  CCTK_REAL fxx = f[0];
-                  CCTK_REAL fxy = f[1];
-                  CCTK_REAL fxz = f[2];
-                  CCTK_REAL fyy = f[3];
-                  CCTK_REAL fyz = f[4];
-                  CCTK_REAL fzz = f[5];
-
-                  txx[idx_dst] = fxx * sxx_val * sxx_val + 2 * fxy * sxx_val * syx_val + fyy * syx_val * syx_val;
-                  tyy[idx_dst] = fxx * sxy_val * sxy_val + 2 * fxy * sxy_val * syy_val + fyy * syy_val * syy_val;
-                  txy[idx_dst] = fxx * sxx_val * sxy_val +     fxy * (sxy_val * syx_val + sxx_val * syy_val) + fyy * syx_val * syy_val;
-                  txz[idx_dst] = fxz * sxx_val + fyz * syx_val;
-                  tyz[idx_dst] = fxz * sxy_val + fyz * syy_val;
-                  tzz[idx_dst] = fzz;
-
-                  //txx[ijk-dj] = fxx*rxx*rxx + 2*fxy*rxx*ryx + fyy*ryx*ryx;
-                  //tyy[ijk-dj] = fxx*rxy*rxy + 2*fxy*rxy*ryy + fyy*ryy*ryy;
-                  //txy[ijk-dj] = fxx*rxx*rxy + fxy*(rxy*ryx+rxx*ryy) + fyy*ryx*ryy;
-                  //txz[ijk-dj] = fxz*rxx + fyz*ryx;
-                  //tyz[ijk-dj] = fxz*rxy + fyz*ryy;
-                  //tzz[ijk-dj] = fzz;
+                  CCTK_REAL *t = cctkGH->data[var0 + n][0];
+                  t[idx_dst] = f[n];
               }
+          } else if (tensortype == TENSORTYPE_U) {
+              CCTK_REAL *tx = cctkGH->data[vi][0];
+              CCTK_REAL *ty = cctkGH->data[vi+1][0];
+              CCTK_REAL *tz = cctkGH->data[vi+2][0];
+
+              CCTK_REAL rxx_val = rxx[idx_dst];
+              CCTK_REAL ryy_val = ryy[idx_dst];
+              CCTK_REAL rxy_val = rxy[idx_dst];
+              CCTK_REAL ryx_val = ryx[idx_dst];
+
+              CCTK_REAL fx = f[0];
+              CCTK_REAL fy = f[1];
+              CCTK_REAL fz = f[2];
+
+              tx[idx_dst] = rxx_val * fx + rxy_val * fy;
+              ty[idx_dst] = ryx_val * fx + ryy_val * fy;
+              tz[idx_dst] = fz;
+
+              //tx[ijk-dj] = sxx * fx + sxy * fy;
+              //ty[ijk-dj] = syx * fx + syy * fy;
+              //tz[ijk-dj] = fz;
+          } else if (tensortype == TENSORTYPE_DDSYM) {
+              CCTK_REAL *txx = cctkGH->data[vi][0];
+              CCTK_REAL *txy = cctkGH->data[vi+1][0];
+              CCTK_REAL *txz = cctkGH->data[vi+2][0];
+              CCTK_REAL *tyy = cctkGH->data[vi+3][0];
+              CCTK_REAL *tyz = cctkGH->data[vi+4][0];
+              CCTK_REAL *tzz = cctkGH->data[vi+5][0];
+
+              CCTK_REAL sxx_val = sxx[idx_dst];
+              CCTK_REAL syy_val = syy[idx_dst];
+              CCTK_REAL sxy_val = sxy[idx_dst];
+              CCTK_REAL syx_val = syx[idx_dst];
+
+              CCTK_REAL fxx = f[0];
+              CCTK_REAL fxy = f[1];
+              CCTK_REAL fxz = f[2];
+              CCTK_REAL fyy = f[3];
+              CCTK_REAL fyz = f[4];
+              CCTK_REAL fzz = f[5];
+
+              txx[idx_dst] = fxx * sxx_val * sxx_val + 2 * fxy * sxx_val * syx_val + fyy * syx_val * syx_val;
+              tyy[idx_dst] = fxx * sxy_val * sxy_val + 2 * fxy * sxy_val * syy_val + fyy * syy_val * syy_val;
+              txy[idx_dst] = fxx * sxx_val * sxy_val +     fxy * (sxy_val * syx_val + sxx_val * syy_val) + fyy * syx_val * syy_val;
+              txz[idx_dst] = fxz * sxx_val + fyz * syx_val;
+              tyz[idx_dst] = fxz * sxy_val + fyz * syy_val;
+              tzz[idx_dst] = fzz;
+
+              //txx[ijk-dj] = fxx*rxx*rxx + 2*fxy*rxx*ryx + fyy*ryx*ryx;
+              //tyy[ijk-dj] = fxx*rxy*rxy + 2*fxy*rxy*ryy + fyy*ryy*ryy;
+              //txy[ijk-dj] = fxx*rxx*rxy + fxy*(rxy*ryx+rxx*ryy) + fyy*ryx*ryy;
+              //txz[ijk-dj] = fxz*rxx + fyz*ryx;
+              //tyz[ijk-dj] = fxz*rxy + fyz*ryy;
+              //tzz[ijk-dj] = fzz;
           }
+      }
 
     /* syncs needed after interpolation (actually only for x direction) */
-    CCTK_SyncGroupI(cctkGH, gi);
+    //CCTK_SyncGroupI(cctkGH, gi);
 
     return 0;
 }
diff --git a/src/RegisterSym.c b/src/RegisterSym.c
index dbed7be..13783ad 100644
--- a/src/RegisterSym.c
+++ b/src/RegisterSym.c
@@ -47,7 +47,7 @@ void Cartoon2D_RegisterSymmetries (CCTK_ARGUMENTS)
   faces[5] = 0;
   
   for (f=0; f<6; ++f) {
-    width[f] = nboundaryzones[f/2];
+    width[f] = f < 2 ? nboundaryzones[f/2] : 0;
   }
   
   handle = SymmetryRegister ("cartoon");
diff --git a/src/SetGrid.c b/src/SetGrid.c
index 28ff2f7..74f3e29 100644
--- a/src/SetGrid.c
+++ b/src/SetGrid.c
@@ -284,8 +284,8 @@ void Cartoon2D_setup_coord(CCTK_ARGUMENTS)
                 rho[idx] = sqrt(x_val * x_val + y_val * y_val);
                 zero[idx] = 0.0;
 
-                rxx[idx] = x_val / rho[idx];
-                ryx[idx] = y_val / rho[idx];
+                rxx[idx] = fabs(rho[idx]) > 1e-8 ? x_val / rho[idx] : 1.0;
+                ryx[idx] = fabs(rho[idx]) > 1e-8 ? y_val / rho[idx] : 0.0;
                 rxy[idx] = -ryx[idx];
                 ryy[idx] = rxx[idx];