aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/ApplyCartoon.c56
-rw-r--r--src/Cartoon2D.h2
-rw-r--r--src/Cartoon2DBC.c167
-rw-r--r--src/Cartoon2D_tensors.h4
-rw-r--r--src/interpolate.c55
5 files changed, 240 insertions, 44 deletions
diff --git a/src/ApplyCartoon.c b/src/ApplyCartoon.c
index a203105..61424e8 100644
--- a/src/ApplyCartoon.c
+++ b/src/ApplyCartoon.c
@@ -25,6 +25,7 @@ CCTK_FILEVERSION(BetaThorns_Cartoon2D_ApplyCartoon_c);
********************************************************************/
/* #define DEBUG 1 */
#define TENSORTYPE_BUFF_SIZE 7
+#define PROLONG_BUFF_SIZE 1000
/********************************************************************
********************* Local Routine Prototypes *********************
@@ -77,11 +78,15 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS);
@endreturndesc
@@*/
-void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) {
+void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS)
+{
+
DECLARE_CCTK_ARGUMENTS;
int num_vars, err, i, gi, group_tags_table;
CCTK_INT * vars;
char tensortype[TENSORTYPE_BUFF_SIZE];
+ char prolongtype[PROLONG_BUFF_SIZE];
+ int prolongmethod;
/* Allocate memory to hold selected bcs */
num_vars = Boundary_SelectedGVs(cctkGH, 0, NULL, NULL, NULL, NULL, NULL);
@@ -145,6 +150,49 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) {
#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_LAGRANGE; /* 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.
@@ -154,13 +202,13 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) {
macro */
if (CCTK_Equals(tensortype, "scalar"))
{
- BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, vars[i]);
+ BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, prolongmethod, vars[i]);
} else if (CCTK_Equals(tensortype, "u"))
{
- BndCartoon2DVI(cctkGH, TENSORTYPE_U, vars[i]);
+ BndCartoon2DVI(cctkGH, TENSORTYPE_U, prolongmethod, vars[i]);
} else if (CCTK_Equals(tensortype, "dd_sym"))
{
- BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, vars[i]);
+ BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, prolongmethod, vars[i]);
} else
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
diff --git a/src/Cartoon2D.h b/src/Cartoon2D.h
index 930c647..476e985 100644
--- a/src/Cartoon2D.h
+++ b/src/Cartoon2D.h
@@ -3,7 +3,7 @@
extern "C" {
#endif
-int BndCartoon2DVI(const cGH *GH, int tensortype, int vi);
+int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi);
#ifdef __cplusplus
}
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;
+}
diff --git a/src/Cartoon2D_tensors.h b/src/Cartoon2D_tensors.h
index 4d0092b..29c72ec 100644
--- a/src/Cartoon2D_tensors.h
+++ b/src/Cartoon2D_tensors.h
@@ -10,3 +10,7 @@
#define N_TENSORTYPE_SCALAR 1
#define N_TENSORTYPE_U 3
#define N_TENSORTYPE_DDSYM 6
+
+/* Prolongation methods */
+#define PROLONG_LAGRANGE 1
+#define PROLONG_ENO 2
diff --git a/src/interpolate.c b/src/interpolate.c
index 8a39c44..d48da6e 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -16,6 +16,7 @@
*/
#include "cctk.h"
+#include "math.h"
#define PRIVATE
#define PUBLIC
@@ -40,6 +41,9 @@ PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
CCTK_REAL y[],
CCTK_REAL x);
+PRIVATE CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+ CCTK_REAL y[],
+ CCTK_REAL x);
/*****************************************************************************/
@@ -225,3 +229,54 @@ CCTK_REAL xr = (x - xc) / dx;
return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))) );
}
+
+/*****************************************************************************/
+
+CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+ CCTK_REAL y[],
+ CCTK_REAL x)
+{
+
+ CCTK_REAL diff[4];
+ CCTK_INT seed = 0;
+ CCTK_INT j,r;
+
+ /* Find seed index */
+ while (x > x0 + dx * ((float)seed+0.5)) seed++;
+
+ if (seed!=2) {
+ /* Not enough stencil, only perform linear interpolation */
+ seed = 0;
+ while (x > x0 + dx * ((float)(seed+1))) seed++;
+ if (seed==4) seed=3;
+ return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(float)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 */
+ double c0 = y[r+1];
+ double c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2];
+ double c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2];
+
+ double xc = x0 + dx * (double)(r+1);
+ double xr = (x - xc) / dx;
+
+ return( c0 + xr*(c1 + xr*c2) );
+
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 06:38:37 +0000