diff options
author | hawke <hawke@eec4d7dc-71c2-46d6-addf-10296150bf52> | 2005-09-14 13:01:59 +0000 |
---|---|---|
committer | hawke <hawke@eec4d7dc-71c2-46d6-addf-10296150bf52> | 2005-09-14 13:01:59 +0000 |
commit | 2b42ff6020d0ab36b11b65466ef6850f7ec68242 (patch) | |
tree | 3a389a7d05130250c70d1f5bdc8590bf845a4781 | |
parent | 49ebff9fcba20d579d938f5142903178f22d1387 (diff) |
Modify the ENO interpolators so that it takes an "order" parameter.
Add knowledge of the PROLONG_NONE type to lower level bits of the
code, in case that case propagates.
git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@94 eec4d7dc-71c2-46d6-addf-10296150bf52
-rw-r--r-- | param.ccl | 7 | ||||
-rw-r--r-- | src/Cartoon2DBC.c | 22 | ||||
-rw-r--r-- | src/Cartoon2D_tensors.h | 1 | ||||
-rw-r--r-- | src/interpolate.c | 186 |
4 files changed, 173 insertions, 43 deletions
@@ -20,6 +20,11 @@ INT order "Cartoon's interpolation order" 1:5 :: "From linear to fifth order." } 4 +INT eno_order "The interpolation order applied to the ENO interpolator" +{ + 1:5 :: "From linear to fifth order." +} 4 + BOOLEAN allow_grid_resize "Allow grid to be resized in a cartoon-compatible way" { -} "no"
\ No newline at end of file +} "no" diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c index ba995d0..f74c30f 100644 --- a/src/Cartoon2DBC.c +++ b/src/Cartoon2DBC.c @@ -39,7 +39,7 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, 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 interpolate_eno(CCTK_INT order, 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); @@ -198,6 +198,10 @@ int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi) /* interpolate grid functions */ switch (prolongtype) { + case PROLONG_NONE: + /* no-op */ + return (0); + break; case PROLONG_LAGRANGE: for (n = 0; n < n_vars; n++) f[n] = Cartoon2DInterp(cctkGH, cctkGH->data[var0+n][timelevel], @@ -383,9 +387,10 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, CCTK_REAL x0; CCTK_REAL lx0, dx0; - CCTK_REAL y[5], r; + CCTK_REAL y[11], r; int lnx; int n, offset; + int tmp_order = 2; lnx = cctkGH->cctk_lsh[0]; @@ -400,25 +405,30 @@ static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, /* offset of stencil, note that rotation leads to x close to x0 for large x */ - offset = 2; + 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+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 <= 4; n++) +/* for (n = 0; n < 2 * tmp_order + 1; ++n) */ + for (n = 0; n < 2 * eno_order + 1; ++n) { y[n] = f[ijk-offset+n]; } /* call interpolator */ - r = interpolate_eno(x0, dx0, y, x); +/* r = interpolate_eno(tmp_order, x0, dx0, y, x); */ + r = interpolate_eno(eno_order, x0, dx0, y, x); return r; } diff --git a/src/Cartoon2D_tensors.h b/src/Cartoon2D_tensors.h index bbb4fa8..7670eab 100644 --- a/src/Cartoon2D_tensors.h +++ b/src/Cartoon2D_tensors.h @@ -9,3 +9,4 @@ /* Prolongation methods */ #define PROLONG_LAGRANGE 1 #define PROLONG_ENO 2 +#define PROLONG_NONE 3 diff --git a/src/interpolate.c b/src/interpolate.c index d990cf5..1d774a2 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -17,6 +17,7 @@ */ #include <math.h> +#include <stdio.h> #include "cctk.h" @@ -24,7 +25,7 @@ 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 interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); @@ -214,58 +215,171 @@ return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))) ); /*****************************************************************************/ -CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx, +CCTK_REAL interpolate_eno(CCTK_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 diff[4]; */ +/* CCTK_INT seed = 0; */ +/* CCTK_INT j,r; */ - CCTK_REAL c0; - CCTK_REAL c1; - CCTK_REAL c2; +/* CCTK_REAL c0; */ +/* CCTK_REAL c1; */ +/* CCTK_REAL c2; */ - CCTK_REAL xc; - CCTK_REAL xr; +/* 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++; +/* while (x > x0 + dx * ((CCTK_REAL)seed+0.5)) seed++; */ - if (seed!=2) { +/* 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); - } +/* 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]; +/* 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; - } +/* 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]; +/* 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; +/* xc = x0 + dx * (CCTK_REAL)(r+1); */ +/* xr = (x - xc) / dx; */ +/* result = ( c0 + xr*(c1 + xr*c2) ); */ - return( 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]; + } + + 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; + } + + 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]; + + 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; + } + +/* #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]); + } + 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; } |