aboutsummaryrefslogtreecommitdiff
path: root/src/Cartoon2DBC.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/Cartoon2DBC.c')
-rw-r--r--src/Cartoon2DBC.c167
1 files changed, 128 insertions, 39 deletions
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index 840c042..2added6 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -19,6 +19,7 @@
#include <math.h>
#include "cctk.h"
+#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_FortranString.h"
@@ -30,11 +31,16 @@ 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 vi);
+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);
CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx,
CCTK_REAL y[], CCTK_REAL x);
+CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+ 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,
@@ -63,10 +69,11 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
group, using the tensortype argument to determine how many there
are and how to treat them. */
-int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
+int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi)
{
- DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
/*** FIXME: can CCTK_SyncGroup() have a 'const cGH *' parameter ?? ***/
union
{
@@ -74,14 +81,14 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
cGH *non_const_ptr;
} GH_fake_const;
- static int xindx=-1;
+/* static int xindx=-1; */
int i, j, k, di, dj, ijk;
int jj, n, s;
int lnx, lnz, ny;
int n_tensortype;
int timelevel;
- CCTK_REAL *xp;
+/* CCTK_REAL *xp; */
CCTK_REAL x, y, rho;
CCTK_REAL lx0, dx0, dy0;
CCTK_REAL rxx, rxy, ryx, ryy;
@@ -105,18 +112,26 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
n_tensortype = 0;
}
- s = GH->cctk_nghostzones[0];
- lnx = GH->cctk_lsh[0];
- lnz = GH->cctk_lsh[2];
- ny = GH->cctk_gsh[1];
+ s = cctkGH->cctk_nghostzones[0];
+ lnx = cctkGH->cctk_lsh[0];
+ lnz = cctkGH->cctk_lsh[2];
+ ny = cctkGH->cctk_gsh[1];
+ /*
xindx = (xindx<0)?CCTK_CoordIndex(-1,"x","cart3d"):xindx;
- xp = GH->data[xindx][0];
+ xp = cctkGH->data[xindx][0];
+ */
/* xp = GH->data[CCTK_CoordIndex(-1,"x","cart3d")][0];*/
+ /*
lx0 = xp[0];
- dx0 = GH->cctk_delta_space[0]/GH->cctk_levfac[0];
- dy0 = GH->cctk_delta_space[1]/GH->cctk_levfac[1];
+ dx0 = cctkGH->cctk_delta_space[0]/cctkGH->cctk_levfac[0];
+ dy0 = cctkGH->cctk_delta_space[1]/cctkGH->cctk_levfac[1];
+ */
+
+ lx0 = CCTK_ORIGIN_SPACE(0);
+ dx0 = CCTK_DELTA_SPACE(0);
+ dy0 = CCTK_DELTA_SPACE(1);
timelevel = 0;
@@ -131,7 +146,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
/* this union helps us to avoid compiler warning about discarding
the const qualifier from a pointer target type */
- GH_fake_const.const_ptr = GH;
+ GH_fake_const.const_ptr = cctkGH;
/* make sure that the input data is synced */
CCTK_SyncGroup(GH_fake_const.non_const_ptr, CCTK_GroupNameFromVarI(vi));
@@ -145,7 +160,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
for (i = -s; i < lnx; i++) {
/* index into linear memory */
- ijk = CCTK_GFINDEX3D(GH,i,j,k);
+ 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;}
@@ -171,19 +186,32 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
syy = rxx;
/* interpolate grid functions */
- for (n = 0; n < n_tensortype; n++)
- f[n] = Cartoon2DInterp(GH, GH->data[vi+n][timelevel], i, di, ijk, rho);
+ switch (prolongtype)
+ {
+ case PROLONG_LAGRANGE:
+ for (n = 0; n < n_tensortype; n++)
+ f[n] = Cartoon2DInterp(cctkGH, cctkGH->data[vi+n][timelevel],
+ i, di, ijk, rho);
+ break;
+ case PROLONG_ENO:
+ for (n = 0; n < n_tensortype; n++)
+ f[n] = Cartoon2DInterp_ENO(cctkGH, cctkGH->data[vi+n][timelevel], i,
+ di, ijk, rho);
+ break;
+ default:
+ CCTK_WARN(0, "Prolongation type not recognized by Cartoon2D.");
+ }
/* rotate grid tensor by matrix multiplication */
if (tensortype == TENSORTYPE_SCALAR) {
- t = GH->data[vi][timelevel];
+ t = cctkGH->data[vi][timelevel];
t[ijk+dj] = f[0];
t[ijk-dj] = f[0];
}
else if (tensortype == TENSORTYPE_U) {
- tx = GH->data[vi][timelevel];
- ty = GH->data[vi+1][timelevel];
- tz = GH->data[vi+2][timelevel];
+ 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];
@@ -197,12 +225,12 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
tz[ijk-dj] = fz;
}
else if (tensortype == TENSORTYPE_DDSYM) {
- txx = GH->data[vi][timelevel];
- txy = GH->data[vi+1][timelevel];
- txz = GH->data[vi+2][timelevel];
- tyy = GH->data[vi+3][timelevel];
- tyz = GH->data[vi+4][timelevel];
- tzz = GH->data[vi+5][timelevel];
+ 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];
@@ -234,7 +262,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
/* this union helps us to avoid compiler warning about discarding
the const qualifier from a pointer target type */
- GH_fake_const.const_ptr = GH;
+ GH_fake_const.const_ptr = cctkGH;
/* syncs needed after interpolation (actually only for x direction) */
CCTK_SyncGroup(GH_fake_const.non_const_ptr, CCTK_GroupNameFromVarI(vi));
@@ -244,7 +272,7 @@ int BndCartoon2DVI(const cGH *GH, int tensortype, int vi)
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)
(int *retval, const cGH **GH, int *tensortype, int *vi)
{
- *retval = BndCartoon2DVI(*GH, *tensortype, *vi);
+ *retval = BndCartoon2DVI(*GH, *tensortype, PROLONG_LAGRANGE, *vi);
}
int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
@@ -256,7 +284,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, vi));
+ return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE, vi));
}
void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)
@@ -299,7 +327,8 @@ int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group)
"%s should have only %d component(s).", group, n_tensortype);
}
- return(BndCartoon2DVI(GH, tensortype, CCTK_FirstVarIndex(group)));
+ return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE,
+ CCTK_FirstVarIndex(group)));
}
void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
@@ -311,31 +340,38 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
}
/* interpolation on x-axis */
-static CCTK_REAL Cartoon2DInterp(const cGH *GH, CCTK_REAL *f, int i, int di,
- int ijk, CCTK_REAL x)
+static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f,
+ int i, int di, int ijk, CCTK_REAL x)
{
- DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
- static int xindx=-1;
+/* static int xindx=-1; */
CCTK_REAL x0;
CCTK_REAL lx0, dx0;
CCTK_REAL y[6], r;
- CCTK_REAL *xp;
+/* CCTK_REAL *xp; */
int lnx;
int n, offset;
(void) (di + 0);
- lnx = GH->cctk_lsh[0];
+ lnx = cctkGH->cctk_lsh[0];
+ /*
xindx = (xindx<0)?CCTK_CoordIndex(-1,"x","cart3d"):xindx;
- xp = GH->data[xindx][0];
-
+ xp = cctkGH->data[xindx][0];
+ */
/* xp = GH->data[CCTK_CoordIndex(-1,"x","cart3d")][0];*/
+/*
lx0 = xp[0];
- dx0 = GH->cctk_delta_space[0]/GH->cctk_levfac[0];
+ dx0 = cctkGH->cctk_delta_space[0]/cctkGH->cctk_levfac[0];
+*/
+
+ lx0 = CCTK_ORIGIN_SPACE(0);
+ dx0 = CCTK_DELTA_SPACE(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) */
@@ -366,3 +402,56 @@ static CCTK_REAL Cartoon2DInterp(const cGH *GH, CCTK_REAL *f, int i, int di,
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;
+
+/* static int xindx=-1; */
+
+ CCTK_REAL x0;
+ CCTK_REAL lx0, dx0;
+ CCTK_REAL y[5], r;
+/* CCTK_REAL *xp; */
+ int lnx;
+ int n, offset;
+
+ lnx = cctkGH->cctk_lsh[0];
+
+ lx0 = CCTK_ORIGIN_SPACE(0);
+ dx0 = CCTK_DELTA_SPACE(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 = 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+4 >= lnx) offset = i+5-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 <= 4; n++)
+ {
+ y[n] = f[ijk-offset+n];
+ }
+
+ /* call interpolator */
+ r = interpolate_eno(x0, dx0, y, x);
+
+ return r;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 04:56:40 +0000