From 8da2fbbd59bf3db10fdc552c302c00c20ba5a229 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sat, 7 Apr 2018 09:40:28 +0200 Subject: tmp --- src/Cartoon2DBC.c | 172 ++++++++++++++++++++++++++---------------------------- src/RegisterSym.c | 2 +- src/SetGrid.c | 4 +- 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]; -- cgit v1.2.3