aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--interface.ccl10
-rw-r--r--schedule.ccl14
-rw-r--r--src/ApplyCartoon.c247
-rw-r--r--src/Cartoon2D.h2
-rw-r--r--src/Cartoon2DBC.c672
-rw-r--r--src/RegisterSym.c10
-rw-r--r--src/SetGrid.c28
-rw-r--r--src/interpolate.c454
8 files changed, 435 insertions, 1002 deletions
diff --git a/interface.ccl b/interface.ccl
index f1df4ec..da70202 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,6 +2,7 @@
# $Header$
implements: cartoon2d
+inherits: grid
INCLUDES HEADER: Cartoon2D_tensors.h in Cartoon2D_tensors.h
INCLUDES HEADER: Cartoon2D.h in Cartoon2D.h
@@ -81,6 +82,15 @@ REQUIRES FUNCTION Boundary_SelectedGVs
private:
+CCTK_REAL coord type=GF
+{
+ rho
+ zero
+ rxx, ryy, rxy, ryx
+ sxx, syy, sxy, syx
+}
+
+
int excision_variables type=SCALAR
{
excision_active
diff --git a/schedule.ccl b/schedule.ccl
index 0ca570f..47a06eb 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -5,6 +5,20 @@ if (cartoon_active)
{
STORAGE: excision_variables
+ STORAGE: coord
+
+ schedule Cartoon2D_setup_coord at CCTK_POSTREGRID {
+ LANG: C
+ } ""
+
+ schedule Cartoon2D_setup_coord at CCTK_POSTREGRIDINITIAL {
+ LANG: C
+ } ""
+
+ schedule Cartoon2D_setup_coord at CCTK_BASEGRID {
+ LANG: C
+ } ""
+
schedule Cartoon2D_CheckTensorTypes at CCTK_PARAMCHECK
{
LANG: C
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<num_vars; ++i) {
-#ifdef DEBUG
- printf("Cartoon_ApplyBoundaries: i=%d applying cartoon to vi %d\n",i,
- vars[i]);
-#endif
- gi = CCTK_GroupIndexFromVarI(vars[i]);
- if (gi<0) {
- CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "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;
- }
-
- /* Here one should check that the entire group is registered,
- * using CCTK_NumVarsInGroupI.
+ /* 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 < num_vars; i++) {
+ char tensortype[TENSORTYPE_BUFF_SIZE];
+ char prolongtype[PROLONG_BUFF_SIZE];
+ int group_tags_table, ret;
- /* 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);
- }
-
-
- /* This lines are writen for Whisky2D to avoid to treat the hydrovariables with Cartoon */
-
- len = Util_TableGetString(group_tags_table,0,NULL,"whiskycartoon");
-
- whiskycartoon = 1; /* default is yes */
- if (len >= 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 <math.h>
#include <stdio.h>
#include <string.h>
#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;
-
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-02 00:12:00 +0000