aboutsummaryrefslogtreecommitdiff
path: root/src/Cartoon2DBC.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/Cartoon2DBC.c')
-rw-r--r--src/Cartoon2DBC.c672
1 files changed, 198 insertions, 474 deletions
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index 2dbca7f..40b3503 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -32,18 +32,9 @@ CCTK_FILEVERSION(BetaThorns_Cartoon2D_Cartoon2DBC_c);
int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group);
int BndCartoon2DVN(const cGH *GH, int tensortype, const char *var);
-int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi);
-static CCTK_REAL Cartoon2DInterp(const cGH *GH, CCTK_REAL *f, int i, int di,
- int ijk, CCTK_REAL x);
-static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
- CCTK_REAL *f, int i, int di,
- int ijk, CCTK_REAL x);
+int BndCartoon2DVI(const cGH *GH, int tensortype, int vi);
CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx,
- CCTK_REAL y[], CCTK_REAL x,
- CCTK_INT excised_points[]);
-CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx,
- CCTK_REAL y[], CCTK_REAL x,
- CCTK_INT excised_points[]);
+ CCTK_REAL y[], CCTK_REAL x);
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)(int *retval, const cGH **GH,
int *tensortype, int *vi);
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)(int *retval, const cGH **GH,
@@ -55,6 +46,46 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS);
void SetupExcisionVars(CCTK_ARGUMENTS);
+/* interpolation on x-axis */
+static CCTK_REAL cartoon2d_interp(const cGH *cctkGH, int order, CCTK_REAL *f,
+ int i, int ijk, CCTK_REAL x_val)
+{
+ CCTK_REAL x0;
+ CCTK_REAL lx0, dx0;
+ int lnx;
+ int n, offset;
+
+ lnx = cctkGH->cctk_lsh[0];
+
+ dx0 = cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0];
+ lx0 = cctkGH->cctk_origin_space[0] + cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0] * cctkGH->cctk_levoff[0] / cctkGH->cctk_levoffdenom[0] + dx0 * cctkGH->cctk_lbnd[0];
+
+
+ /* find i such that x(i) < x <= x(i+1)
+ for rotation on entry always x > x(i), but sometimes also x > x(i+1) */
+ while (x_val > lx0 + dx0 * (i+1)) {i++; ijk++;}
+
+ /* offset of stencil, note that rotation leads to x close to x0
+ for large x */
+ offset = order/2;
+
+ /* shift stencil at boundaries */
+ /* note that for simplicity we don't distinguish between true boundaries
+ and decomposition boundaries: the sync fixes that! */
+ if (i-offset < 0) offset = i;
+ while (lx0 + dx0 * (i - offset) < 0) offset--;
+
+ if (i-offset+order >= lnx) offset = i+order-lnx+1;
+
+ x0 = lx0 + dx0 * (i-offset);
+ if(x0 < 0)
+ fprintf(stderr, "x0 %g\n", x0);
+
+ /* call interpolator */
+ return interpolate_local(order, x0, dx0, f + ijk - offset, x_val);
+}
+
+
/* set boundaries of a grid tensor assuming axisymmetry
- handles lower boundary in x
- does not touch other boundaries
@@ -75,241 +106,175 @@ void SetupExcisionVars(CCTK_ARGUMENTS);
group, using the tensortype argument to determine how many there
are and how to treat them. */
-int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi)
+int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi)
{
- DECLARE_CCTK_ARGUMENTS;
- DECLARE_CCTK_PARAMETERS;
-
- char * groupname;
-
- int i, j, k, di, dj, ijk;
- int jj, n, s;
- int lnx, lnz, ny;
- int gi, var0;
- int n_vars;
- int timelevel;
- CCTK_REAL x, y, rho;
- CCTK_REAL lx0, dx0, dy0;
- CCTK_REAL rxx, rxy, ryx, ryy;
- CCTK_REAL sxx, sxy, syx, syy;
- CCTK_REAL f[100], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz;
- CCTK_REAL *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz;
-
- gi = CCTK_GroupIndexFromVarI (vi);
- assert (gi >= 0);
- n_vars = CCTK_NumVarsInGroupI (gi);
- assert (n_vars > 0);
- var0 = CCTK_FirstVarIndexI (gi);
- assert (var0 >= 0);
-
- if (n_vars > 100)
- {
- groupname = CCTK_GroupName (gi);
- CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Group \"%s\" has more than 100 variables -- cannot apply Cartoon boundary condition",
- groupname);
- free (groupname);
- return -1;
- }
-
- switch (tensortype) {
- case TENSORTYPE_SCALAR:
- break;
- case TENSORTYPE_U:
- if (n_vars != 3)
- {
- groupname = CCTK_GroupName (gi);
- CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Group \"%s\" is declared with the tensor type U, but does not contain 3 variables -- cannot apply Cartoon boundary condition",
- groupname);
- free (groupname);
- return -1;
- }
- break;
- case TENSORTYPE_DDSYM:
- if (n_vars != 6)
- {
- groupname = CCTK_GroupName (gi);
- CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Group \"%s\" is declared with the tensor type DDSYM, but does not contain 6 variables -- cannot apply Cartoon boundary condition",
- groupname);
- free (groupname);
- return -1;
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
+
+ int i, s;
+ int y_idx;
+ int gi, var0;
+ int n_vars;
+ int timelevel;
+
+ gi = CCTK_GroupIndexFromVarI (vi);
+ n_vars = CCTK_NumVarsInGroupI (gi);
+ var0 = CCTK_FirstVarIndexI (gi);
+
+ if (n_vars > 100) {
+ char *groupname = CCTK_GroupName (gi);
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Group \"%s\" has more than 100 variables -- cannot apply Cartoon boundary condition",
+ groupname);
+ free (groupname);
+ return -1;
}
- break;
- default:
- CCTK_WARN(0,"Tensor type not recognized by Cartoon2D.");
- }
-
- s = cctkGH->cctk_nghostzones[0];
- lnx = cctkGH->cctk_lsh[0];
- lnz = cctkGH->cctk_lsh[2];
- ny = cctkGH->cctk_gsh[1];
-
- dx0 = CCTK_DELTA_SPACE(0);
- dy0 = CCTK_DELTA_SPACE(1);
- lx0 = CCTK_ORIGIN_SPACE(0) + dx0 * cctk_lbnd[0];
-
- timelevel = 0;
-
- /* Cactus loop in C: grid points with y = 0 */
- /* compare thorn_BAM_Elliptic */
- /* strides used in stencils, should be provided by cactus */
- di = 1;
- dj = lnx;
-
- /* y = 0 */
- assert (ny % 2 == 1);
- j = ny/2;
-
- /* make sure that the input data is synced */
- CCTK_SyncGroupI(cctkGH, gi);
-
- /* z-direction: include lower and upper boundary */
- for (k = 0; k < lnz; k++)
-
- /* y-direction: as many zombies as the evolution stencil needs */
-#if 0
- for (jj = 1, dj = jj*lnx; jj <= ny/2; jj++, dj = jj*lnx)
-#else
- for (jj = 0, dj = jj*lnx; jj <= ny/2; jj++, dj = jj*lnx)
-#endif
-
- /* x-direction: zombie for x < 0, including upper boundary for extrapol */
-#if 0
- for (i = -s; i < lnx; i++) {
-#else
- for (i = 0; i < lnx; i++)
- if (! (i >= s && jj == 0)) {
-#endif
-
- /* index into linear memory */
- ijk = CCTK_GFINDEX3D(cctkGH,i,j,k);
-
- /* what a neat way to hack in the zombie for x < 0, y = 0 */
- if (i < 0) {i += s; ijk += s; dj = 0;}
-
- /* compute coordinates (could also use Cactus grid function) */
- x = lx0 + dx0 * i;
-#if 0
- y = (dj) ? dy0 * jj : 0;
-#else
- y = dy0 * jj;
-#endif
- rho = sqrt(x*x + y*y);
-
- /* compute rotation matrix
- note that this also works for x <= 0 (at lower boundary in x)
- note that rho is nonzero by definition if cactus checks dy ...
- */
- rxx = x/rho;
- ryx = y/rho;
- rxy = -ryx;
- ryy = rxx;
-
- /* inverse rotation matrix, assuming detr = 1 */
- sxx = ryy;
- syx = -ryx;
- sxy = -rxy;
- syy = rxx;
-
- /* interpolate grid functions */
- switch (prolongtype)
- {
- case PROLONG_NONE:
- /* no-op */
- return (0);
+
+ switch (tensortype) {
+ case TENSORTYPE_SCALAR:
break;
- case PROLONG_LAGRANGE:
- for (n = 0; n < n_vars; n++)
- f[n] = Cartoon2DInterp(cctkGH, cctkGH->data[var0+n][timelevel],
- i, di, ijk, rho);
+ case TENSORTYPE_U:
+ if (n_vars != 3) {
+ char *groupname = CCTK_GroupName (gi);
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Group \"%s\" is declared with the tensor type U, but does not contain 3 variables -- cannot apply Cartoon boundary condition",
+ groupname);
+ free (groupname);
+ return -1;
+ }
break;
- case PROLONG_ENO:
- for (n = 0; n < n_vars; n++)
- f[n] = Cartoon2DInterp_ENO(cctkGH, cctkGH->data[var0+n][timelevel], i,
- di, ijk, rho);
+ case TENSORTYPE_DDSYM:
+ if (n_vars != 6) {
+ char *groupname = CCTK_GroupName (gi);
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Group \"%s\" is declared with the tensor type DDSYM, but does not contain 6 variables -- cannot apply Cartoon boundary condition",
+ groupname);
+ free (groupname);
+ return -1;
+ }
break;
- default:
- CCTK_WARN(0, "Prolongation type not recognized by Cartoon2D.");
- }
-
- /* rotate grid tensor by matrix multiplication */
- if (tensortype == TENSORTYPE_SCALAR) {
- /* Scalar groups can have arbitrary many components; these
- "components" are then all scalars. */
- for (n = 0; n < n_vars; n++)
- {
- t = cctkGH->data[var0 + n][timelevel];
- t[ijk+dj] = f[n];
- t[ijk-dj] = f[n];
- }
- }
- else if (tensortype == TENSORTYPE_U) {
- tx = cctkGH->data[vi][timelevel];
- ty = cctkGH->data[vi+1][timelevel];
- tz = cctkGH->data[vi+2][timelevel];
- fx = f[0];
- fy = f[1];
- fz = f[2];
-
- tx[ijk+dj] = rxx * fx + rxy * fy;
- ty[ijk+dj] = ryx * fx + ryy * fy;
- tz[ijk+dj] = fz;
-
- tx[ijk-dj] = sxx * fx + sxy * fy;
- ty[ijk-dj] = syx * fx + syy * fy;
- tz[ijk-dj] = fz;
- }
- else if (tensortype == TENSORTYPE_DDSYM) {
- txx = cctkGH->data[vi][timelevel];
- txy = cctkGH->data[vi+1][timelevel];
- txz = cctkGH->data[vi+2][timelevel];
- tyy = cctkGH->data[vi+3][timelevel];
- tyz = cctkGH->data[vi+4][timelevel];
- tzz = cctkGH->data[vi+5][timelevel];
- fxx = f[0];
- fxy = f[1];
- fxz = f[2];
- fyy = f[3];
- fyz = f[4];
- fzz = f[5];
-
- txx[ijk+dj] = fxx*sxx*sxx + 2*fxy*sxx*syx + fyy*syx*syx;
- tyy[ijk+dj] = fxx*sxy*sxy + 2*fxy*sxy*syy + fyy*syy*syy;
- txy[ijk+dj] = fxx*sxx*sxy + fxy*(sxy*syx+sxx*syy) + fyy*syx*syy;
- txz[ijk+dj] = fxz*sxx + fyz*syx;
- tyz[ijk+dj] = fxz*sxy + fyz*syy;
- tzz[ijk+dj] = 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;
- }
- else {
- return(-1);
+ default:
+ CCTK_WARN(0,"Tensor type not recognized by Cartoon2D.");
}
-#if 0
- /* what a neat way to hack out the zombies for x < 0, y = 0 */
- if (dj == 0) {i -= s; ijk -= s; dj = jj*lnx;}
-#endif
- }
+ s = cctkGH->cctk_nghostzones[0];
- /* syncs needed after interpolation (actually only for x direction) */
- CCTK_SyncGroupI(cctkGH, gi);
+ timelevel = 0;
- return(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) {
+ y_idx = i;
+ break;
+ }
+ }
+ if (i == cctkGH->cctk_lsh[1])
+ CCTK_WARN(0, "y is not 0 anywhere");
+
+ /* make sure that the input data is synced */
+ 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];
+
+ /* interpolate grid functions */
+ 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;
+ }
+ }
+
+ /* syncs needed after interpolation (actually only for x direction) */
+ CCTK_SyncGroupI(cctkGH, gi);
+
+ return 0;
}
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)
(int *retval, const cGH **GH, int *tensortype, int *vi)
{
- *retval = BndCartoon2DVI(*GH, *tensortype, PROLONG_LAGRANGE, *vi);
+ *retval = BndCartoon2DVI(*GH, *tensortype, *vi);
}
int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
@@ -321,7 +286,7 @@ int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
if (verbose) printf("Cartoon2D called for %s\n", inpvarname);
vi = CCTK_VarIndex(inpvarname);
- return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE, vi));
+ return(BndCartoon2DVI(GH, tensortype, vi));
}
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)
@@ -338,8 +303,7 @@ int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group)
if (verbose && GH->cctk_iteration==1) CCTK_VInfo(CCTK_THORNSTRING,"Cartoon2D called for %s (message appears only once)", group);
- return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE,
- CCTK_FirstVarIndex(group)));
+ return(BndCartoon2DVI(GH, tensortype, CCTK_FirstVarIndex(group)));
}
void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
@@ -350,246 +314,6 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
free(group);
}
-/* interpolation on x-axis */
-static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f,
- int i, int di, int ijk, CCTK_REAL x)
-{
- DECLARE_CCTK_ARGUMENTS;
- DECLARE_CCTK_PARAMETERS;
-
- CCTK_REAL x0;
- CCTK_REAL lx0, dx0;
- CCTK_REAL y[6], r;
- int lnx;
- int n, offset;
-
- CCTK_REAL *old_mask;
- CCTK_INT *new_mask;
-
- CCTK_INT excised_points[6];
-
- if (*excision_active < 0)
- {
- SetupExcisionVars((cGH*)cctkGH);
- }
-
- if (old_excision)
- {
- old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
- }
- else
- {
- old_mask = NULL;
- }
-
- if (new_excision)
- {
- new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
- }
- else
- {
- new_mask = NULL;
- }
-
- /*
- If this boundary point is excised then we return
- without doing anything
- */
-
- if (*excision_active)
- {
- if ((*excision_active == 1)||(*excision_active == 3))
- {
- if (old_mask[ijk] < 0.75)
- {
- return f[ijk];
- }
- }
- if ((*excision_active == 2)||(*excision_active == 3))
- {
- if (SpaceMask_CheckStateBits(new_mask, (ijk),
- *new_excision_field,
- *new_excision_descriptor))
- {
- return f[ijk];
- }
- }
- }
-
- lnx = cctkGH->cctk_lsh[0];
-
- dx0 = CCTK_DELTA_SPACE(0);
- lx0 = CCTK_ORIGIN_SPACE(0) + dx0 * cctk_lbnd[0];
-
- /* find i such that x(i) < x <= x(i+1)
- for rotation on entry always x > x(i), but sometimes also x > x(i+1) */
- while (x > lx0 + dx0 * (i+1)) {i++; ijk++;}
-
- /* first attempt to interface to JT's interpolator */
-
- /* offset of stencil, note that rotation leads to x close to x0
- for large x */
- offset = order/2;
-
- /* shift stencil at boundaries */
- /* note that for simplicity we don't distinguish between true boundaries
- and decomposition boundaries: the sync fixes that! */
- if (i-offset < 0) offset = i;
- if (i-offset+order >= lnx) offset = i+order-lnx+1;
-
- /* fill in info */
- /* fills in old data near axis, but as long as stencil at axis is
- centered, no interpolation happens anyway */
- x0 = lx0 + dx0 * (i-offset);
- for (n = 0; n <= order; n++) {
- y[n] = f[ijk-offset+n];
- excised_points[n] = 0;
- if (*excision_active)
- {
- if ((*excision_active == 1)||(*excision_active == 3))
- {
- excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 0 : 1;
- }
- if ((*excision_active == 2)||(*excision_active == 3))
- {
- excised_points[n] |=
- (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
- *new_excision_field,
- *new_excision_descriptor));
- }
- }
- }
-
- /* call interpolator */
- r = interpolate_local(order, x0, dx0, y, x, excised_points);
-
- return r;
-}
-
-static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH,
- CCTK_REAL *f, int i, int di,
- int ijk, CCTK_REAL x)
-{
-
- DECLARE_CCTK_PARAMETERS;
- DECLARE_CCTK_ARGUMENTS;
-
- CCTK_REAL x0;
- CCTK_REAL lx0, dx0;
- CCTK_REAL y[11], r;
- int lnx;
- int n, offset;
- /* int tmp_order = 2; */
-
- CCTK_REAL *old_mask;
- CCTK_INT *new_mask;
-
- CCTK_INT excised_points[11];
-
- if (*excision_active < 0)
- {
- SetupExcisionVars((cGH*)cctkGH);
- }
-
- if (old_excision)
- {
- old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
- }
- else
- {
- old_mask = NULL;
- }
-
- if (new_excision)
- {
- new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
- }
- else
- {
- new_mask = NULL;
- }
-
- /*
- If this boundary point is excised then we return
- without doing anything
- */
-
- if (*excision_active)
- {
- if ((*excision_active == 1)||(*excision_active == 3))
- {
- if (old_mask[ijk] < 0.75)
- {
- return f[ijk];
- }
- }
- if ((*excision_active == 2)||(*excision_active == 3))
- {
- if (SpaceMask_CheckStateBits(new_mask, (ijk),
- *new_excision_field,
- *new_excision_descriptor))
- {
- return f[ijk];
- }
- }
- }
-
- lnx = cctkGH->cctk_lsh[0];
-
- dx0 = CCTK_DELTA_SPACE(0);
- lx0 = CCTK_ORIGIN_SPACE(0) + dx0 * cctk_lbnd[0];
-
- /* find i such that x(i) < x <= x(i+1)
- for rotation on entry always x > x(i), but sometimes also x > x(i+1) */
- while (x > lx0 + dx0 * (i+1)) {i++; ijk++;}
-
- /* first attempt to interface to JT's interpolator */
-
- /* offset of stencil, note that rotation leads to x close to x0
- for large x */
- offset = eno_order;
-/* offset = tmp_order; */
-
- /* shift stencil at boundaries */
- /* note that for simplicity we don't distinguish between true boundaries
- and decomposition boundaries: the sync fixes that! */
- if (i-offset < 0) offset = i;
-/* if (i-offset+4 >= lnx) offset = i+5-lnx; */
-/* if (i-offset+2*tmp_order >= lnx) offset = i+2*tmp_order+1-lnx; */
- if (i-offset+2*eno_order >= lnx) offset = i+2*eno_order+1-lnx;
-
- /* fill in info */
- /* fills in old data near axis, but as long as stencil at axis is
- centered, no interpolation happens anyway */
- x0 = lx0 + dx0 * (i-offset);
-/* for (n = 0; n < 2 * tmp_order + 1; ++n) */
- for (n = 0; n < 2 * eno_order + 1; ++n)
- {
- y[n] = f[ijk-offset+n];
- excised_points[n] = 0;
- if (*excision_active)
- {
- if ((*excision_active == 1)||(*excision_active == 3))
- {
- excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 1 : 0;
- }
- if ((*excision_active == 2)||(*excision_active == 3))
- {
- excised_points[n] |=
- (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
- *new_excision_field,
- *new_excision_descriptor));
- }
- }
- }
-
- /* call interpolator */
-/* r = interpolate_eno(tmp_order, x0, dx0, y, x); */
- r = interpolate_eno(eno_order, x0, dx0, y, x, excised_points);
-
- return r;
-}
-
void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-15 19:13:30 +0000