From 140eaa74a9e30828636c5bedafbda6723db758de Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 26 Feb 2016 17:15:35 +0100 Subject: Largely rewrite to work with mesh refinement. As the actual work is now mostly done in the evolution code, this really only handles the symmetry condition on the x-axis. --- src/ApplyCartoon.c | 247 ++++++-------------- src/Cartoon2D.h | 2 +- src/Cartoon2DBC.c | 672 ++++++++++++++++------------------------------------- src/RegisterSym.c | 10 +- src/SetGrid.c | 28 +++ src/interpolate.c | 454 +++++++++--------------------------- 6 files changed, 411 insertions(+), 1002 deletions(-) (limited to 'src') diff --git a/src/ApplyCartoon.c b/src/ApplyCartoon.c index b97c289..e5f3bc0 100644 --- a/src/ApplyCartoon.c +++ b/src/ApplyCartoon.c @@ -84,189 +84,86 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS); void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; + int num_vars, err, i; + CCTK_INT * vars; - int num_vars, err, i, gi, group_tags_table; - CCTK_INT * vars; - char tensortype[TENSORTYPE_BUFF_SIZE]; - char prolongtype[PROLONG_BUFF_SIZE]; - int prolongmethod = -1; + //return; - int whiskycartoon, len; + /* Allocate memory to hold selected bcs */ + num_vars = Boundary_SelectedGVs(cctkGH, 0, NULL, NULL, NULL, NULL, NULL); + vars = malloc(num_vars*sizeof *vars); - - - - /* Check grid size */ - if(cctk_gsh[1] != 2*cctk_nghostzones[1]+1) - { - CCTK_WARN(0, "Grid size in y direction inappropriate."); - } - - /* Allocate memory to hold selected bcs */ - num_vars = Boundary_SelectedGVs(cctkGH, 0, NULL, NULL, NULL, NULL, NULL); -#ifdef DEBUG - printf("Cartoon_ApplyBoundaries: num_vars is %d\n",num_vars); -#endif - vars = malloc(num_vars*sizeof *vars); - - /* get all selected vars */ - err = Boundary_SelectedGVs(cctkGH, num_vars, vars, NULL, NULL, NULL, NULL); - if (err != num_vars) - { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Boundary_SelectedGVs returned %d selected variables, but %d " - "expected", err, num_vars); + /* get all selected vars */ + err = Boundary_SelectedGVs(cctkGH, num_vars, vars, NULL, NULL, NULL, NULL); + if (err != num_vars) { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Boundary_SelectedGVs returned %d selected variables, but %d " + "expected", err, num_vars); } - /* Apply CartoonBC to each of them. */ - /* One should probably check to see that the entire group has been - * selected, since Cartoon operates on an entire group at a time. - * For now I'll just skip a variable if it is not a 'group leader'. - */ - for (i=0; i= 0) - { - char* value = malloc (len + 1); - Util_TableGetString (group_tags_table, len + 1, value, "whiskycartoon"); - if (CCTK_Equals(value, "yes")) - { - whiskycartoon = 1; - } - else if (CCTK_Equals(value, "no")) - { - whiskycartoon = 0; - } - else - { + int gi = CCTK_GroupIndexFromVarI(vars[i]); + if (gi < 0) { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Tags table entry \"whiskycartoon\" for variable group %s should be either \"yes\" or \"no\", but is \"%s\"", - CCTK_GroupName(gi), value); - } - free (value); - } - - if (whiskycartoon) - { - - /*###################################################################*/ - - - /* Get tensor type from group tags table */ - err = Util_TableGetString(group_tags_table, TENSORTYPE_BUFF_SIZE, - tensortype, "tensortypealias"); - if (err<0) - { - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Error in TAGS table for variable group %s", - CCTK_GroupName(gi)); - } -#ifdef DEBUG - printf("Cartoon_ApplyBoundaries: tensor type is %s\n",tensortype); -#endif - - /* Get prolongation type from group tags table */ - - err = Util_TableGetString(group_tags_table, PROLONG_BUFF_SIZE, - prolongtype, "Prolongation"); - - if (err == UTIL_ERROR_TABLE_NO_SUCH_KEY) - { - /* Use the default */ - - prolongmethod = PROLONG_LAGRANGE; - } - else if (err < 0) - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Error (%d) in TAGS table for variable group %s", err, - CCTK_GroupName(gi)); - } - else - { - if (CCTK_Equals(prolongtype, "None")) - { - prolongmethod = PROLONG_NONE; /* But why? */ - } - else if (CCTK_Equals(prolongtype, "Lagrange")) - { - prolongmethod = PROLONG_LAGRANGE; - } - else if (CCTK_Equals(prolongtype, "TVD")) - { - prolongmethod = PROLONG_ENO; - } - else if (CCTK_Equals(prolongtype, "ENO")) - { - prolongmethod = PROLONG_ENO; - } - else - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Error in TAGS table for variable group %s", - CCTK_GroupName(gi)); - } - } - - /* Here one should check that the group's size is correct for the - * specified tensortype. - */ + "Invalid variable index %d selected for a boundary " + "condition", gi); + } + if (vars[i] != CCTK_FirstVarIndexI(gi)) { + /* not a group leader -- skip to next var index */ + continue; + } + + /* Get table handle for group tags table */ + group_tags_table = CCTK_GroupTagsTableI(gi); + if (group_tags_table < 0) { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Tags table for variable group %s is %d", CCTK_GroupName(gi), + group_tags_table); + } + + /* Get tensor type from group tags table */ + ret = Util_TableGetString(group_tags_table, TENSORTYPE_BUFF_SIZE, + tensortype, "tensortypealias"); + if (ret < 0) { + CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in TAGS table for variable group %s", + CCTK_GroupName(gi)); + } + + /* Get prolongation type from group tags table */ + ret = Util_TableGetString(group_tags_table, PROLONG_BUFF_SIZE, + prolongtype, "Prolongation"); + + if (!(ret == UTIL_ERROR_TABLE_NO_SUCH_KEY || + (ret >= 0 && CCTK_Equals(prolongtype, "Lagrange")))) + continue; + + /* Call BndCartoon2DVI, passing the appropriate tensor type integer + macro */ + if (CCTK_Equals(tensortype, "scalar")) { + BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, vars[i]); + } else if ((CCTK_Equals(tensortype, "u")) || (CCTK_Equals(tensortype, "d"))) { + BndCartoon2DVI(cctkGH, TENSORTYPE_U, vars[i]); + } else if (CCTK_Equals(tensortype, "dd_sym")) { + BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, vars[i]); + } else { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "invalid tensor type for group %s", CCTK_GroupName(gi)); + } + } - /* Call BndCartoon2DVI, passing the appropriate tensor type integer - macro */ - if (CCTK_Equals(tensortype, "scalar")) - { - BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, prolongmethod, vars[i]); - } else if ((CCTK_Equals(tensortype, "u")) || - (CCTK_Equals(tensortype, "d"))) - { - BndCartoon2DVI(cctkGH, TENSORTYPE_U, prolongmethod, vars[i]); - } else if (CCTK_Equals(tensortype, "dd_sym")) - { - BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, prolongmethod, vars[i]); - } else - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "invalid tensor type for group %s", CCTK_GroupName(gi)); - } - } - } - /* Free data */ - free(vars); + /* Free data */ + free(vars); } diff --git a/src/Cartoon2D.h b/src/Cartoon2D.h index b4fadc3..d371729 100644 --- a/src/Cartoon2D.h +++ b/src/Cartoon2D.h @@ -7,7 +7,7 @@ extern "C" { #endif -int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi); +int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi); CCTK_INT Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_, 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; diff --git a/src/RegisterSym.c b/src/RegisterSym.c index 2145df0..dbed7be 100644 --- a/src/RegisterSym.c +++ b/src/RegisterSym.c @@ -60,10 +60,10 @@ void Cartoon2D_RegisterSymmetries (CCTK_ARGUMENTS) CCTK_WARN (0, "Could not register the Cartoon boundaries -- probably some other thorn has already registered the same boundary faces for a different symmetry"); } - ierr = SymmetryRegisterGridInterpolator - (cctkGH, handle, Cartoon2D_SymmetryInterpolate); - if (ierr < 0) { - CCTK_WARN (0, "Could not register the symmetry interpolator"); - } + //ierr = SymmetryRegisterGridInterpolator + // (cctkGH, handle, Cartoon2D_SymmetryInterpolate); + //if (ierr < 0) { + // CCTK_WARN (0, "Could not register the symmetry interpolator"); + //} } diff --git a/src/SetGrid.c b/src/SetGrid.c index 954d301..28ff2f7 100644 --- a/src/SetGrid.c +++ b/src/SetGrid.c @@ -9,10 +9,12 @@ @version $Id$ @@*/ +#include #include #include #include "cctk.h" +#include "cctk_Arguments.h" #include "cctk_Parameter.h" #include "cctk_Parameters.h" @@ -267,3 +269,29 @@ int Cartoon2D_SetGrid(void) return 0; } + +void Cartoon2D_setup_coord(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + +#pragma omp parallel for + for (int k = 0; k < cctkGH->cctk_lsh[2]; k++) + for (int j = 0; j < cctkGH->cctk_lsh[1]; j++) + for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) { + int idx = CCTK_GFINDEX3D(cctkGH, i, j, k); + CCTK_REAL x_val = x[idx], y_val = y[idx]; + 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]; + rxy[idx] = -ryx[idx]; + ryy[idx] = rxx[idx]; + + sxx[idx] = ryy[idx]; + syx[idx] = -ryx[idx]; + sxy[idx] = -rxy[idx]; + syy[idx] = rxx[idx]; + } +} diff --git a/src/interpolate.c b/src/interpolate.c index 3ae6b3f..151f08c 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -23,128 +23,7 @@ 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[]); - -/* prototypes for private functions defined in this file */ -static CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], - CCTK_REAL x); -static CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], - CCTK_REAL x); -static CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], - CCTK_REAL x); -static CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], - CCTK_REAL x); -static CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], - CCTK_REAL x); - -/*****************************************************************************/ - -/* - * This function does local Lagrange polynomial interpolation on a - * uniform grid. That is, given order+1 uniformly spaced data points - * (x0 + i*dx, y[i]) , i = 0...order, it interpolates a degree-(order) - * (Lagrange) polynomial through them and evaluates this polynomial at - * a specified x coordinate. In general the error in this computed - * function value is O(h^(order+1)) . - * - * Except for possible end effects (see end_action (below)), the local - * interpolation is internally centered to improve its conditioning. Even - * so, the interpolation is highly ill-conditioned for small dx and/or - * large order and/or x outside the domain of the data points. - * - * The interpolation formulas were (are) all derived via Maple, see - * "./interpolate.in" and "./interpolate.out". - * - * Arguments: - * order = (in) The order of the local interpolation, i.e. the degree - * of the interpolated polynomial. - * x0 = (in) The x coordinate corresponding to y[0]. - * dx = (in) The x spacing between the data points. - * y[0...order] = (in) The y data array. - * x = (in) The x coordinate for the interpolation. - */ -CCTK_REAL interpolate_local(int order, - CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], - CCTK_REAL x, - CCTK_INT excised_points[]) -{ - - /* - We assume that the excised region has one of the following 4 forms: - - o o o o o (No excision) - x o o o o (left excision) - o o o o x (right excision) - x o o o x (left and right excision) - - Anything else (e.g., x o x o x) will go horribly wrong. - */ - - int start, end, available_order; - - start = 0; - end = order; - - while ((start < order + 1)&&(excised_points[start])) - { - ++start; - } - while ((end > -1)&&(excised_points[end])) - { - --end; - } - - if (end < start) - { - /* The whole block is excised */ - return y[0]; - } - - /* The possible interpolation order that can be used */ - - available_order = end - start; - - if (available_order == 0) - { - /* There is only one non-excised point */ - return y[start]; - } - - if (dx * start + x0 > x) - { - /* We would have to extrapolate to the left */ - return y[start]; - } - - if (dx * end + x0 < x) - { - /* We would have to extrapolate to the right */ - return y[end]; - } - - switch (available_order) - { - case 1: return interpolate_local_order1(x0 + start * dx, dx, &y[start], x); - case 2: return interpolate_local_order2(x0 + start * dx, dx, &y[start], x); - case 3: return interpolate_local_order3(x0 + start * dx, dx, &y[start], x); - case 4: return interpolate_local_order4(x0 + start * dx, dx, &y[start], x); - case 5: return interpolate_local_order5(x0 + start * dx, dx, &y[start], x); - default: CCTK_WARN(0, "Interpolation order not supported"); - return 0; /* Avoid warning */ - } -} + CCTK_REAL x); /*****************************************************************************/ @@ -156,13 +35,13 @@ static CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { -CCTK_REAL c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1]; -CCTK_REAL c1 = - y[0] + + y[1]; + CCTK_REAL c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1]; + CCTK_REAL c1 = - y[0] + + y[1]; -CCTK_REAL xc = x0 + 0.5*dx; -CCTK_REAL xr = (x - xc) / dx; + CCTK_REAL xc = x0 + 0.5*dx; + CCTK_REAL xr = (x - xc) / dx; -return( c0 + xr*c1 ); + return c0 + xr * c1; } /*****************************************************************************/ @@ -175,14 +54,14 @@ static CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { -CCTK_REAL c0 = y[1]; -CCTK_REAL c1 = (-1.0/2.0)*y[0] + (+1.0/2.0)*y[2]; -CCTK_REAL c2 = (+1.0/2.0)*y[0] - y[1] + (+1.0/2.0)*y[2]; + CCTK_REAL c0 = y[1]; + CCTK_REAL c1 = (-1.0/2.0)*y[0] + (+1.0/2.0)*y[2]; + CCTK_REAL c2 = (+1.0/2.0)*y[0] - y[1] + (+1.0/2.0)*y[2]; -CCTK_REAL xc = x0 + 1.0*dx; -CCTK_REAL xr = (x - xc) / dx; + CCTK_REAL xc = x0 + 1.0*dx; + CCTK_REAL xr = (x - xc) / dx; -return( c0 + xr*(c1 + xr*c2) ); + return c0 + xr*(c1 + xr*c2); } /*****************************************************************************/ @@ -195,19 +74,19 @@ static CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { -CCTK_REAL c0 = (-1.0/16.0)*y[0] + (+9.0/16.0)*y[1] - + (+9.0/16.0)*y[2] + (-1.0/16.0)*y[3]; -CCTK_REAL c1 = (+1.0/24.0)*y[0] + (-9.0/ 8.0)*y[1] - + (+9.0/ 8.0)*y[2] + (-1.0/24.0)*y[3]; -CCTK_REAL c2 = (+1.0/ 4.0)*y[0] + (-1.0/ 4.0)*y[1] - + (-1.0/ 4.0)*y[2] + (+1.0/ 4.0)*y[3]; -CCTK_REAL c3 = (-1.0/ 6.0)*y[0] + (+1.0/ 2.0)*y[1] - + (-1.0/ 2.0)*y[2] + ( 1.0/ 6.0)*y[3]; - -CCTK_REAL xc = x0 + 1.5*dx; -CCTK_REAL xr = (x - xc) / dx; - -return( c0 + xr*(c1 + xr*(c2 + xr*c3)) ); + CCTK_REAL c0 = (-1.0/16.0)*y[0] + (+9.0/16.0)*y[1] + + (+9.0/16.0)*y[2] + (-1.0/16.0)*y[3]; + CCTK_REAL c1 = (+1.0/24.0)*y[0] + (-9.0/ 8.0)*y[1] + + (+9.0/ 8.0)*y[2] + (-1.0/24.0)*y[3]; + CCTK_REAL c2 = (+1.0/ 4.0)*y[0] + (-1.0/ 4.0)*y[1] + + (-1.0/ 4.0)*y[2] + (+1.0/ 4.0)*y[3]; + CCTK_REAL c3 = (-1.0/ 6.0)*y[0] + (+1.0/ 2.0)*y[1] + + (-1.0/ 2.0)*y[2] + ( 1.0/ 6.0)*y[3]; + + CCTK_REAL xc = x0 + 1.5*dx; + CCTK_REAL xr = (x - xc) / dx; + + return c0 + xr*(c1 + xr*(c2 + xr*c3)); } /*****************************************************************************/ @@ -220,20 +99,20 @@ static CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { -CCTK_REAL c0 = y[2]; -CCTK_REAL c1 = (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1] - + (+2.0/ 3.0)*y[3] + (-1.0/12.0)*y[4]; -CCTK_REAL c2 = + (-1.0/24.0)*y[0] + (+2.0/ 3.0)*y[1] + (-5.0/ 4.0)*y[2] - + (+2.0/ 3.0)*y[3] + (-1.0/24.0)*y[4]; -CCTK_REAL c3 = + (-1.0/12.0)*y[0] + (+1.0/ 6.0)*y[1] - + (-1.0/ 6.0)*y[3] + (+1.0/12.0)*y[4]; -CCTK_REAL c4 = + (+1.0/24.0)*y[0] + (-1.0/ 6.0)*y[1] + (+1.0/ 4.0)*y[2] - + (-1.0/ 6.0)*y[3] + (+1.0/24.0)*y[4]; - -CCTK_REAL xc = x0 + 2.0*dx; -CCTK_REAL xr = (x - xc) / dx; - -return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4))) ); + CCTK_REAL c0 = y[2]; + CCTK_REAL c1 = (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1] + + (+2.0/ 3.0)*y[3] + (-1.0/12.0)*y[4]; + CCTK_REAL c2 = + (-1.0/24.0)*y[0] + (+2.0/ 3.0)*y[1] + (-5.0/ 4.0)*y[2] + + (+2.0/ 3.0)*y[3] + (-1.0/24.0)*y[4]; + CCTK_REAL c3 = + (-1.0/12.0)*y[0] + (+1.0/ 6.0)*y[1] + + (-1.0/ 6.0)*y[3] + (+1.0/12.0)*y[4]; + CCTK_REAL c4 = + (+1.0/24.0)*y[0] + (-1.0/ 6.0)*y[1] + (+1.0/ 4.0)*y[2] + + (-1.0/ 6.0)*y[3] + (+1.0/24.0)*y[4]; + + CCTK_REAL xc = x0 + 2.0*dx; + CCTK_REAL xr = (x - xc) / dx; + + return c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4))); } /*****************************************************************************/ @@ -246,200 +125,81 @@ static CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { -CCTK_REAL c0 = (+ 3.0/256.0)*y[0] + (-25.0/256.0)*y[1] - + (+75.0/128.0)*y[2] + (+75.0/128.0)*y[3] - + (-25.0/256.0)*y[4] + (+ 3.0/256.0)*y[5]; -CCTK_REAL c1 = (- 3.0/640.0)*y[0] + (+25.0/384.0)*y[1] - + (-75.0/ 64.0)*y[2] + (+75.0/ 64.0)*y[3] - + (-25.0/384.0)*y[4] + (+ 3.0/640.0)*y[5]; -CCTK_REAL c2 = (- 5.0/ 96.0)*y[0] + (+13.0/ 32.0)*y[1] - + (-17.0/ 48.0)*y[2] + (-17.0/ 48.0)*y[3] - + (+13.0/ 32.0)*y[4] + (- 5.0/ 96.0)*y[5]; -CCTK_REAL c3 = (+ 1.0/ 48.0)*y[0] + (-13.0/ 48.0)*y[1] - + (+17.0/ 24.0)*y[2] + (-17.0/ 24.0)*y[3] - + (+13.0/ 48.0)*y[4] + (- 1.0/ 48.0)*y[5]; -CCTK_REAL c4 = (+ 1.0/ 48.0)*y[0] + (- 1.0/ 16.0)*y[1] - + (+ 1.0/ 24.0)*y[2] + (+ 1.0/ 24.0)*y[3] - + (- 1.0/ 16.0)*y[4] + (+ 1.0/ 48.0)*y[5]; -CCTK_REAL c5 = (- 1.0/120.0)*y[0] + (+ 1.0/ 24.0)*y[1] - + (- 1.0/ 12.0)*y[2] + (+ 1.0/ 12.0)*y[3] - + (- 1.0/ 24.0)*y[4] + (+ 1.0/120.0)*y[5]; - -CCTK_REAL xc = x0 + 2.5*dx; -CCTK_REAL xr = (x - xc) / dx; - -return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))) ); + CCTK_REAL c0 = (+ 3.0/256.0)*y[0] + (-25.0/256.0)*y[1] + + (+75.0/128.0)*y[2] + (+75.0/128.0)*y[3] + + (-25.0/256.0)*y[4] + (+ 3.0/256.0)*y[5]; + CCTK_REAL c1 = (- 3.0/640.0)*y[0] + (+25.0/384.0)*y[1] + + (-75.0/ 64.0)*y[2] + (+75.0/ 64.0)*y[3] + + (-25.0/384.0)*y[4] + (+ 3.0/640.0)*y[5]; + CCTK_REAL c2 = (- 5.0/ 96.0)*y[0] + (+13.0/ 32.0)*y[1] + + (-17.0/ 48.0)*y[2] + (-17.0/ 48.0)*y[3] + + (+13.0/ 32.0)*y[4] + (- 5.0/ 96.0)*y[5]; + CCTK_REAL c3 = (+ 1.0/ 48.0)*y[0] + (-13.0/ 48.0)*y[1] + + (+17.0/ 24.0)*y[2] + (-17.0/ 24.0)*y[3] + + (+13.0/ 48.0)*y[4] + (- 1.0/ 48.0)*y[5]; + CCTK_REAL c4 = (+ 1.0/ 48.0)*y[0] + (- 1.0/ 16.0)*y[1] + + (+ 1.0/ 24.0)*y[2] + (+ 1.0/ 24.0)*y[3] + + (- 1.0/ 16.0)*y[4] + (+ 1.0/ 48.0)*y[5]; + CCTK_REAL c5 = (- 1.0/120.0)*y[0] + (+ 1.0/ 24.0)*y[1] + + (- 1.0/ 12.0)*y[2] + (+ 1.0/ 12.0)*y[3] + + (- 1.0/ 24.0)*y[4] + (+ 1.0/120.0)*y[5]; + + CCTK_REAL xc = x0 + 2.5*dx; + CCTK_REAL xr = (x - xc) / dx; + + return c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))); } - -/*****************************************************************************/ - -CCTK_REAL interpolate_eno(CCTK_INT order, - CCTK_REAL x0, CCTK_REAL dx, - CCTK_REAL y[], - CCTK_REAL x, - CCTK_INT excised_points[]) +/* + * This function does local Lagrange polynomial interpolation on a + * uniform grid. That is, given order+1 uniformly spaced data points + * (x0 + i*dx, y[i]) , i = 0...order, it interpolates a degree-(order) + * (Lagrange) polynomial through them and evaluates this polynomial at + * a specified x coordinate. In general the error in this computed + * function value is O(h^(order+1)) . + * + * Except for possible end effects (see end_action (below)), the local + * interpolation is internally centered to improve its conditioning. Even + * so, the interpolation is highly ill-conditioned for small dx and/or + * large order and/or x outside the domain of the data points. + * + * The interpolation formulas were (are) all derived via Maple, see + * "./interpolate.in" and "./interpolate.out". + * + * Arguments: + * order = (in) The order of the local interpolation, i.e. the degree + * of the interpolated polynomial. + * x0 = (in) The x coordinate corresponding to y[0]. + * dx = (in) The x spacing between the data points. + * y[0...order] = (in) The y data array. + * x = (in) The x coordinate for the interpolation. + */ +CCTK_REAL interpolate_local(int order, + CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x) { -/* CCTK_REAL diff[4]; */ -/* CCTK_INT seed = 0; */ -/* CCTK_INT j,r; */ - -/* CCTK_REAL c0; */ -/* CCTK_REAL c1; */ -/* CCTK_REAL c2; */ - -/* CCTK_REAL xc; */ -/* CCTK_REAL xr; */ - - CCTK_REAL result; - CCTK_REAL undiv_diff [13][6]; - CCTK_INT diff, i, ii, r; - CCTK_REAL fp_i, fp_ii, x_rel; - - /* Find seed index */ -/* while (x > x0 + dx * ((CCTK_REAL)seed+0.5)) seed++; */ - -/* if (seed!=2) { */ - /* Not enough stencil, only perform linear interpolation */ -/* seed = 0; */ -/* while (x > x0 + dx * ((CCTK_REAL)(seed+1))) seed++; */ -/* if (seed==4) seed=3; */ -/* return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(CCTK_REAL)seed*dx); */ -/* } */ - - /* Calculate the undivided differences */ -/* for (j=0; j<=3; j++) */ -/* diff[j] = y[j+1] - y[j]; */ - - /* Find the stencil */ -/* if ( fabs(diff[1]) < fabs(diff[2]) ) { */ -/* if ( fabs(diff[1]-diff[0]) < fabs(diff[2]-diff[1]) ) */ -/* r = 0; */ -/* else */ -/* r = 1; */ -/* } else { */ -/* if ( fabs(diff[2]-diff[1]) < fabs(diff[3]-diff[2]) ) */ -/* r = 1; */ -/* else */ -/* r = 2; */ -/* } */ - - /* Interpolate second order */ -/* c0 = y[r+1]; */ -/* c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2]; */ -/* c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2]; */ - -/* xc = x0 + dx * (CCTK_REAL)(r+1); */ -/* xr = (x - xc) / dx; */ -/* result = ( c0 + xr*(c1 + xr*c2) ); */ - - for (i = 1; i < 2 * order + 2; ++i) - { - undiv_diff[i][0] = y[i-1]; - } - for (i = 0; i < 6; ++i) - { - undiv_diff[0][i] = 1e10; - undiv_diff[2 * order + 2][i] = 1e10; - } - - for (diff = 1; diff < order + 1; ++diff) - { - for (i = 1; i < 2 * order + 2; ++i) - { - undiv_diff[i][diff] = undiv_diff[i][diff-1] - undiv_diff[i-1][diff-1]; - } - undiv_diff[0][diff] = -10 * undiv_diff[1][diff]; - undiv_diff[2 * order + 2][diff]= -10 * undiv_diff[2 * order + 1][diff]; - } + int start, end; - fp_i = (x - x0) / dx; - fp_ii = floor(fp_i); - x_rel = fp_i - fp_ii; - - ii = (CCTK_INT)(fp_ii) + 1; - - if (ii < 1) - { - ii = 1; - x_rel = fp_i - (CCTK_REAL)(ii) + 1.0; - } - else if (ii > 2 * order + 1) - { - ii = 2 * order + 1; - x_rel = fp_i - (CCTK_REAL)(ii) + 1.0; - } + start = 0; + end = order; - r = 0; - for (diff = 1; diff < order + 1; ++diff) - { - if (fabs(undiv_diff[ii-r][diff]) < fabs(undiv_diff[ii-r+1][diff]) ) - { - ++r; - } - - if (ii - r < diff + 1) - { - --r; - } - else if (ii - r + diff > 2 * order + 1) - { - ++r; - } - } - - result = undiv_diff[ii-r][0]; + /* The possible interpolation order that can be used */ - switch (order) - { - case 1: - result += (x_rel + r - 0.0) / 1.0 * undiv_diff[ii-r+1][1]; - break; - case 2: - result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] + - (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2])); - break; - case 3: - result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] + - (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] + - (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3]))); - break; - case 4: - result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] + - (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] + - (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3] + - (x_rel + r - 3.0) / 4.0 * (undiv_diff[ii-r+4][4])))); - break; - case 5: - result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] + - (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] + - (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3] + - (x_rel + r - 3.0) / 4.0 * (undiv_diff[ii-r+4][4] + - (x_rel + r - 4.0) / 5.0 * (undiv_diff[ii-r+5][5]))))); - break; - } + if (x0 > x) + return y[start]; -/* #define CARTOON_ENO_DBG 1 */ -#ifdef CARTOON_ENO_DBG - - printf("Cartoon ENO interpolation.\n" - "Input: order %d x0 %g dx %g x %g.\n", - order, x0, dx, x); - printf("Data is\n"); - for (i = 0; i < 2 * order + 1; ++i) - { - printf(" %g", y[i]); + if (dx * end + x0 < x) + return y[end]; + + switch (order) { + case 1: return interpolate_local_order1(x0, dx, y, x); + case 2: return interpolate_local_order2(x0, dx, y, x); + case 3: return interpolate_local_order3(x0, dx, y, x); + case 4: return interpolate_local_order4(x0, dx, y, x); + case 5: return interpolate_local_order5(x0, dx, y, x); + default: CCTK_WARN(0, "Interpolation order not supported"); + return 0; /* Avoid warning */ } - printf("\nEntries used were\n"); - printf("0: %g 1: %g 2: %g\n", - undiv_diff[ii-r][0], undiv_diff[ii-r+1][1], undiv_diff[ii-r+2][2]); - printf("Parameters are fp_i %g fp_ii %g x_rel %g ii %d r %d\n", - fp_i, fp_ii, x_rel, ii, r); - printf("Result is %g\n", result); - -#endif - - return result; - } -- cgit v1.2.3