diff options
-rw-r--r-- | src/ApplyCartoon.c | 56 | ||||
-rw-r--r-- | src/Cartoon2D.h | 2 | ||||
-rw-r--r-- | src/Cartoon2DBC.c | 167 | ||||
-rw-r--r-- | src/Cartoon2D_tensors.h | 4 | ||||
-rw-r--r-- | src/interpolate.c | 55 |
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) ); + +} |