diff options
-rw-r--r-- | interface.ccl | 2 | ||||
-rw-r--r-- | schedule.ccl | 5 | ||||
-rw-r--r-- | src/Cartoon2DBC.c | 243 | ||||
-rw-r--r-- | src/SetBoundAxisym.c | 226 | ||||
-rw-r--r-- | src/error_exit.c | 63 | ||||
-rw-r--r-- | src/interpolate.c | 1072 | ||||
-rw-r--r-- | src/make.code.defn | 2 |
7 files changed, 310 insertions, 1303 deletions
diff --git a/interface.ccl b/interface.ccl index 52065c5..e83a4fc 100644 --- a/interface.ccl +++ b/interface.ccl @@ -1,2 +1,4 @@ # Interface definition for thorn Cartoon2D # $Header$ + +implements: cartoon2d diff --git a/schedule.ccl b/schedule.ccl index 730af80..0012a25 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,2 +1,7 @@ # Schedule definitions for thorn Cartoon2D # $Header$ + +schedule Cartoon2D_CheckParameters at CCTK_PARAMCHECK +{ + LANG: C +} "Check Cartoon2D parameters" diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c index f1e1ebf..f7dfe70 100644 --- a/src/Cartoon2DBC.c +++ b/src/Cartoon2DBC.c @@ -4,46 +4,261 @@ @author Sai Iyer @desc Apply Cartoon2D boundary conditions - Adapted from Bernd Bruegmann's cartoon_2d_tensorbound.c - and set_bound_axisym.c + An implementation of Steve Brandt's idea for doing + axisymmetry with 3d stencils. + Cactus 4 version of Bernd Bruegmann's code. @enddesc @@*/ static char *rcsid = "$Id$" +#include <stdio.h> +#include <assert.h> +#include <stdlib.h> +#include <ctype.h> +#include <stdarg.h> +#include <string.h> +#include <math.h> + #include "cctk.h" -#include "cctk_arguments.h" #include "cctk_parameters.h" +#include "cctk_FortranString.h" #include "cartoon2D.h" -int Cartoon2DBCVarI(cGH *GH, int vi, int tensortype) { +/* Number of components */ +CCTK_INT tensor_type_n[N_TENSOR_TYPES] = {1, 3, 6}; + + +CCTK_REAL Cartoon2DInterp(cGH *GH, CCTK_REAL *f, CCTK_INT i, CCTK_INT di, + CCTK_INT ijk, CCTK_REAL x); +CCTK_INT Cartoon2DBCVarI(cGH *GH, CCTK_INT vi, CCTK_INT tensor_type); +CCTK_INT Cartoon2DBCVar(cGH *GH, const char *impvarname, CCTK_INT tensortype); + + +void FMODIFIER FORTRAN_NAME(Cartoon2DBCVarI)(CCTK_INT *retval, cGH *GH, + CCTK_INT *vi, CCTK_INT *tensortype); +void FMODIFIER FORTRAN_NAME(Cartoon2DBCVar)(CCTK_INT *retval, cGH *GH, + ONE_FORTSTRING_ARG, CCTK_INT *tensortype); + +/* set boundaries of a grid tensor assuming axisymmetry + - handles lower boundary in x + - does not touch other boundaries + - coordinates and rotation coefficients are independent of z and should + be precomputed + - this is also true for the constants in interpolator, but this may + be too messy + - minimizes conceptual warpage, not computationally optimized +*/ +/* uses rotation matrix and its inverse as linear transformation on + arbitrary tensor indices -- I consider this a good compromise + between doing index loops versus using explicit formulas in + cos(phi) etc with messy signs +*/ + +CCTK_INT Cartoon2DBCVarI(cGH *GH, CCTK_INT vi, CCTK_INT tensor_type) + +{ DECLARE_CCTK_PARAMETERS - static int firstcall = 1, pr = 0; + CCTK_INT i, j, k, di, dj, dk, ijk; + CCTK_INT jj, n, s; + CCTK_INT lnx, lny, lnz, ny; + CCTK_REAL lx0, dx0, dy0; + CCTK_REAL rxx, rxy, ryx, ryy; + CCTK_REAL sxx, sxy, syx, syy; + CCTK_REAL f[6], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz; + CCTK_REAL *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz; - if (firstcall) { - firstcall = 0; - pr = CCTK_Equals(verbose, "yes"); + s = GH->cctk_nghostzones[0]; + lnx = GH->cctk_lsh[0]; + lny = GH->cctk_lsh[1]; + lnz = GH->cctk_lsh[2]; + ny = GH->cctk_gsh[1]; + lx0 = GH->data[CCTK_CoordIndex("x")][0][0]; + dx0 = GH->cctk_delta_space[0]/GH->cctk_levfac[0]; + dy0 = GH->cctk_delta_space[1]/GH->cctk_levfac[1]; + + /* 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; + dk = lnx*lny; + + /* y = 0 */ + j = ny/2; + + /* z-direction: include lower and upper boundary */ + for (k = 0; k < lnz; k++) + + /* y-direction: as many zombies as the evolution stencil needs */ + for (jj = 1, dj = jj*lnx; jj <= ny/2; jj++, dj = jj*lnx) + + /* x-direction: zombie for x < 0, including upper boundary for extrapol */ + for (i = -s; i < lnx; i++) { + + /* index into linear memory */ + ijk = CCTK_GFINDEX3D(GH,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; + y = (dj) ? dy0 * jj : 0; + 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 */ + for (n = 0; n < tensor_type_n[tensor_type]; n++) + f[n] = axisym_interpolate(GH, GH->data[vi+n][0], i, di, ijk, rho); + + /* rotate grid tensor by matrix multiplication */ + if (tensor_type == TENSOR_TYPE_SCALAR) { + t = GH->data[vi][0]; + t[ijk+dj] = f[0]; + t[ijk-dj] = f[0]; + } + else if (tensor_type == TENSOR_TYPE_U) { + tx = GH->data[vi][0]; + ty = GH->data[vi+1][0]; + tz = GH->data[vi+2][0]; + 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 (tensor_type == TENSOR_TYPE_DDSYM) { + txx = GH->data[vi][0]; + txy = GH->data[vi+1][0]; + txz = GH->data[vi+2][0]; + tyy = GH->data[vi+3][0]; + tyz = GH->data[vi+4][0]; + tzz = GH->data[vi+5][0]; + 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); + } + + /* what a neat way to hack out the zombies for x < 0, y = 0 */ + if (dj == 0) {i -= s; ijk -= s; dj = jj*lnx;} } - if (pr) printf("cartoon_2d_tensorbound called for %s\n", name); - return(SetBoundAxisym(GH, vi, tensortype)); + /* syncs needed after interpolation (actually only for x direction) */ + for (n = 0; n < tensor_type_n[tensor_type]; n++) + CCTK_SyncGroup(GH, CCTK_GroupNameFromVarI[vi+n]); + + return(0); } -void FMODIFIER FORTRAN_NAME(Cartoon2DBCVarI)(int *retval, cGH *GH, int *vi, int *tensortype) { +void FMODIFIER FORTRAN_NAME(Cartoon2DBCVarI)(CCTK_INT *retval, cGH *GH, CCTK_INT *vi, CCTK_INT *tensortype) { *retval = Cartoon2DBCVarI(GH, *vi, *tensortype); } -int Cartoon2DBCVar(cGH *GH, const char *impvarname, int tensortype) +CCTK_INT Cartoon2DBCVar(cGH *GH, const char *impvarname, CCTK_INT tensortype) { - int vi; + CCTK_INT vi; + static CCTK_INT firstcall = 1, pr = 0; + + if (firstcall) { + firstcall = 0; + pr = CCTK_Equals(verbose, "yes"); + } + if (pr) printf("cartoon_2d_tensorbound called for %s\n", name); vi = CCTK_VarIndex(impvarname); return(Cartoon2DBCVarI(GH, vi, tensortype)); } -void FMODIFIER FORTRAN_NAME(Cartoon2DBCVar)(int *retval, cGH *GH, ONE_FORTSTRING_ARG, int *tensortype) { +void FMODIFIER FORTRAN_NAME(Cartoon2DBCVar)(CCTK_INT *retval, cGH *GH, ONE_FORTSTRING_ARG, CCTK_INT *tensortype) { ONE_FORTSTRING_CREATE(impvarname) *retval = Cartoon2DBCVar(GH, impvarname, *tensortype); free(impvarname); } + +/* interpolation on x-axis */ +CCTK_REAL Cartoon2DInterp(cGH *GH, CCTK_REAL *f, CCTK_INT i, CCTK_INT di, + CCTK_INT ijk, CCTK_REAL x) +{ + DECLARE_CCTK_PARAMETERS + + CCTK_REAL c, d, x0; + CCTK_REAL lx0, dx0; + CCTK_REAL y[6], r; + CCTK_INT n, offset; + + lnx = GH->cctk_lsh[0]; + lx0 = GH->data[CCTK_CoordIndex("x")][0][0]; + dx0 = GH->cctk_delta_space[0]/GH->cctk_levfac[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]; + } + + /* call interpolator */ + r = interpolate_local(order, x0, dx0, y, x); + + return r; +} diff --git a/src/SetBoundAxisym.c b/src/SetBoundAxisym.c deleted file mode 100644 index f5a4f7d..0000000 --- a/src/SetBoundAxisym.c +++ /dev/null @@ -1,226 +0,0 @@ -/* set_bound_axisym.c, BB */ -/*@@ - @file SetBoundAxisym.c - @date Sat Nov 6 16:35:46 MET 1999 - @author Sai Iyer - @desc - An implementation of Steve Brandt's idea for doing - axisymmetry with 3d stencils. - Cactus 4 version of Bernd Bruegmann's set_bound_axisym_GT. - @enddesc - @@*/ - -static char *rcsid = "$Id$" - -#include "cctk.h" -#include "cctk_arguments.h" -#include "cctk_parameters.h" - -#include "cartoon2D.h" - -#include "math.h" - -/* Number of components */ -int tensor_type_n[N_TENSOR_TYPES] = {1, 3, 6}; - -/* set boundaries of a grid tensor assuming axisymmetry - - handles lower boundary in x - - does not touch other boundaries - - coordinates and rotation coefficients are independent of z and should - be precomputed - - this is also true for the constants in interpolator, but this may - be too messy - - minimizes conceptual warpage, not computationally optimized -*/ -/* uses rotation matrix and its inverse as linear transformation on - arbitrary tensor indices -- I consider this a good compromise - between doing index loops versus using explicit formulas in - cos(phi) etc with messy signs -*/ -int SetBoundAxiSym(cGH *GH, int vi, int tensortype) { - -{ - int i, j, k, di, dj, dk, ijk; - int jj, n, s; - Double x, y, rho; - Double rxx, rxy, ryx, ryy; - Double sxx, sxy, syx, syy; - Double f[6], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz; - Double *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz; - - - /* stencil width: in x-direction use ghost zone width, - in y-direction use ny/2, see below */ - s = cctk_nghostzones[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 = GH->lnx; - dk = GH->lnx*GH->lny; - - /* y = 0 */ - j = GH->ny/2; - - /* z-direction: include lower and upper boundary */ - for (k = 0; k < GH->lnz; k++) - - /* y-direction: as many zombies as the evolution stencil needs */ - for (jj = 1, dj = jj*GH->lnx; jj <= GH->ny/2; jj++, dj = jj*GH->lnx) - - /* x-direction: zombie for x < 0, including upper boundary for extrapol */ - for (i = -s; i < GH->lnx; i++) { - - /* index into linear memory */ - ijk = DI(GH,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 = GH->lx0 + GH->dx0 * i; - y = (dj) ? GH->dy0 * jj : 0; - 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 */ - for (n = 0; n < GT->n; n++) - f[n] = axisym_interpolate(GH, GT->gf[n], i, di, ijk, rho); - - /* rotate grid tensor by matrix multiplication */ - if (GT->tensor_type == TENSOR_TYPE_SCALAR) { - t = GT->gf[0]->data; - t[ijk+dj] = f[0]; - t[ijk-dj] = f[0]; - } - else if (GT->tensor_type == TENSOR_TYPE_U) { - tx = GT->gf[0]->data; - ty = GT->gf[1]->data; - tz = GT->gf[2]->data; - 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 (GT->tensor_type == TENSOR_TYPE_DDSYM) { - txx = GT->gf[0]->data; - txy = GT->gf[1]->data; - txz = GT->gf[2]->data; - tyy = GT->gf[3]->data; - tyz = GT->gf[4]->data; - tzz = GT->gf[5]->data; - 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 { - /* oops */ - } - - /* what a neat way to hack out the zombies for x < 0, y = 0 */ - if (dj == 0) {i -= s; ijk -= s; dj = jj*GH->lnx;} - } - - /* syncs needed after interpolation (actually only for x direction) */ - for (n = 0; n < GT->n; n++) - SyncGF(GH, GT->gf[n]); -} - - - - - -/* interpolation on x-axis */ -Double axisym_interpolate(pGH *GH, pGF *GF, int i, int di, int ijk, Double x) -{ - Double *f = GF->data; - Double c, d, x0; - - /* 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 > GH->lx0 + GH->dx0 * (i+1)) {i++; ijk++;} - - /* the original demo */ - if (0) { - c = (f[ijk+di] - f[ijk])/GH->dx0; - d = f[ijk]; - x0 = GH->lx0 + GH->dx0 * i; - return c*(x-x0) + d; - } - - /* first attempt to interface to JT's interpolator */ - if (1) { - double interpolate_local(int order, double x0, double dx, - double y[], double x); - double y[6], r; - int n, offset; - int order = cartoon_order; - - /* 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 >= GH->lnx) offset = i+order-GH->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 = GH->lx0 + GH->dx0 * (i-offset); - for (n = 0; n <= order; n++) { - y[n] = f[ijk-offset+n]; - if (0) printf("y[%d]=%6.4f ", n, y[n]); - } - - /* call interpolator */ - r = interpolate_local(order, x0, GH->dx0, y, x); - - if (0) printf("i %d x %7.4f x0 %7.4f r %7.4f\n", i, x, x0, r); - return r; - } - - return 0; -} diff --git a/src/error_exit.c b/src/error_exit.c deleted file mode 100644 index 748de99..0000000 --- a/src/error_exit.c +++ /dev/null @@ -1,63 +0,0 @@ -/* error_exit.c -- print an error message and cleanly shutdown cactus */ -/* $Id$ */ - -#include <stdio.h> -#include <stdarg.h> -#include <stdlib.h> - -#include "cactus.h" -#include "thorn_cartoon_2d/src/error_exit.h" - -/* - * Jonathan Thornburg's favorite syntactic sugar to make the two branches - * of an if-else statement symmetrical: - * - * Also note (in)famous joke saying - * "syntactic sugar causes cancer of the semicolon" - * :) :) :) - */ -#define then /* empty */ - -/*****************************************************************************/ - -/* - * This function prints an error message (formatted via vfprintf(3S)), - * then cleanly shuts down cactus. It does *not* return to its caller. - * - * Despite not actually returning, this function is declared as returning - * an int , so it may be easily used in conditional expressions like - * foo = bar_ok ? baz(bar) : error_exit(...); - * - * Usage: - * error_exit(exit_code, message, args...) - * where exit_code is one of the codes defined in "error_exit.h" - * - * ***** THIS VERSION IS FOR ANSI/ISO C ***** - * - * Arguments: - * exit_code = (in) One of the codes defined in "error_exit.h", or any - * other integer to pass back to the shell as our exit - * status. - * format = (in) vprintf(3S) format string for error message to be printed. - * args... = (in) Any additional arguments are (presumably) used in formatting - * the error message string. - */ -#ifdef __cplusplus -extern "C" -#endif -/*VARARGS*/ -/* PUBLIC */ -int error_exit(const int exit_code, const char *const format, ...) -{ -va_list ap; - -va_start(ap, format); -vfprintf(stderr, format, ap); -va_end(ap); - - if (exit_code == 0) { - then STOP; /*NOTREACHED*/ - } else if (exit_code > 0) { - then ERROR_STOP(exit_code); /*NOTREACHED*/ - } else ERROR_ABORT(exit_code); /*NOTREACHED*/ -} diff --git a/src/interpolate.c b/src/interpolate.c index 346a1fa..fe8cbb0 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -1,247 +1,45 @@ /* interpolate.c -- interpolation functions */ /* -** <<<design notes>>> - - * interpolate_global - Lagrange polynomial interpolation (global) + * * interpolate_local - ... (local) ** interpolate_local_order1 - ... (order 1) (linear) (2-point) ** interpolate_local_order2 - ... (order 2) (3-point) ** interpolate_local_order3 - ... (order 3) (4-point) ** interpolate_local_order4 - ... (order 4) (5-point) ** interpolate_local_order5 - ... (order 5) (6-point) - - * interpolate_deriv_global - Lagrange polynomial differentiation (global) - * interpolate_deriv_local - ... (local) -** interpolate_deriv_local_order1 - ... (order 1) (linear) (2-point) -** interpolate_deriv_local_order2 - ... (order 2) (3-point) -** interpolate_deriv_local_order3 - ... (order 3) (4-point) -** interpolate_deriv_local_order4 - ... (order 4) (5-point) -** interpolate_deriv_local_order5 - ... (order 5) (6-point) - - * interpolate_nu_global - Lagrange polynomial interp (nonuniform/global) - * interpolate_nu_local - ... (nonuniform/local) -** interpolate_nu_local_order3 - ... (order 3) (4-point) -** interpolate_nu_local_order4 - ... (order 4) (5-point) -** interpolate_nu_local_order5 - ... (order 5) (6-point) - - * subarray_posn - compute subarray positions for interpolation etc */ +/* To make porting easy, only the interpolate_local* functions + * have been retained in this file. JT's original code (kept in + * the archive directory) also has functions for derivatives and + * interpolation on non-uniform grids. ---Sai. +*/ +#include "cctk.h" -/* minimal thing to do so that this files compiles in cactus, BB */ -#define CACTUS_STANDALONE - -#ifdef CACTUS_STANDALONE - -#include "cactus.h" -#include "thorn_cartoon_2d/src/error_exit.h" #define PRIVATE #define PUBLIC -#define then -#define boolean int -#define FALSE 0 -#define TRUE 1 -#define HOW_MANY_IN_RANGE(a,b) ((b) - (a) + 1) -#define floor_double_to_int(x) ((int)floor((double)x)) - -double interpolate_local(int order, double x0, double dx, - double y[], double x); -double interpolate_deriv_local(int order, double x0, double dx, - double y[], double x); -double interpolate_nu_local(int order, double xx[], double yy[], double x); - -#else - -#include <stdio.h> -#include <jt/stdc.h> -#include <jt/util.h> -#include <jt/interpolate.h> -#endif +CCTK_REAL interpolate_local(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], CCTK_REAL x); /* prototypes for private functions defined in this file */ -PRIVATE double interpolate_local_order1(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_local_order2(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_local_order3(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_local_order4(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_local_order5(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_deriv_local_order1(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_deriv_local_order2(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_deriv_local_order3(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_deriv_local_order4(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_deriv_local_order5(double x0, double dx, - double y[], - double x); -PRIVATE double interpolate_nu_local_order3(double xx[], - double yy[], - double x); -PRIVATE double interpolate_nu_local_order4(double xx[], - double yy[], - double x); -PRIVATE double interpolate_nu_local_order5(double xx[], - double yy[], - double x); +PRIVATE CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x); +PRIVATE CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x); +PRIVATE CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x); +PRIVATE CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x); +PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x); -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -/* - * ***** Design Notes ***** - */ - -/* - * The basic problem of interpolation is this: given a set of data points - * (x[i],y[i]) sampled from a smooth function y = y(x) , we wish to - * estimate y or one if its derivatives, at arbitrary x values. - * - * We classify interpolation functions along 3 dimensions: - * - local (the underlying routines, where all data points are used), - * vs global (search for the right place in a large array, then do a - * local interpolation), - * - uniform (the x[i] values are uniformly spaced, and may be implicitly - * specified as a linear function x[i] = x0 + i*dx ), vs non-uniform - * (the x[i] values are (potentially) non-uniformly spaced, and hence - * must be supplied to the interpolation function as an explicit array), - * - function (we want to estimate y ) vs derivative (we want to estimate - * dy/dx ; higher derivatives are also possible, but typically don't - * occur in an interpolation context). - */ - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -/* - * This function does global Lagrange polynomial interpolation on a - * uniform grid. That is, given N_points uniformly spaced data points - * (x0 + i*dx, y[i]) , i = 0...N_points-1, and given a specified x - * coordinate, it locates a subarray of order+1 data points centered - * around x , interpolates a degree-order (Lagrange) polynomial through - * them, and evaluates this at x . - * - * 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: - * N_points = (in) The number of data points. - * order = (in) The order of the local interpolation, i.e. the degree - * of the interpolated polynomial. - * end_action = (in) Specifies what action to take if the (centered) - * local interpolant would need data from outside - * the domain of the data points. - * 'e' ==> Do an error_exit() . - * 'o' ==> Off-center the local interpolant as - * necessary to keep it inside the domain - * of the data points. - * x0 = (in) The x coordinate corresponding to y[0] . - * dx = (in) The x spacing between the data points. - * y[0...N_points) = (in) The y data array. - * x = (in) The x coordinate for the interpolation. - */ -PUBLIC double interpolate_global(int N_points, int order, - char end_action, - double x0, double dx, - double y[], - double x) -{ -int N_interp, i, i_lo, i_hi, shift; - -/* - * locate subarray for local interpolation - */ -N_interp = order + 1; -i = floor_double_to_int( (x - x0)/dx ); -subarray_posn(N_interp, i, &i_lo, &i_hi); -if (! ((i_lo >= 0) && (i_hi < N_points)) ) - then { - /* centered local interpolation would need data */ - /* from outside y[] array */ - switch (end_action) - { - case 'e': - error_exit(ERROR_EXIT, -"***** interpolate_global: local interpolation needs data\n" -" from outside data array!\n" -" N_points=%d order=%d N_interp=%d\n" -" x0=%g dx=%g\n" -" x=%g\n" -" ==> i=%d i_[lo,hi]=[%d,%d]\n" -, - N_points, order, N_interp, - x0, dx, - x, - i, i_lo, i_hi); /*NOTREACHED*/ - - case 'o': - /* off-center the local interpolation as necessary */ - /* to keep it within the y[] array */ - shift = 0; - if (i_lo < 0) then shift = 0 - i_lo; - if (i_hi >= N_points) then shift = (N_points-1) - i_hi; - - i_lo += shift; - i_hi += shift; - - if (! ((i_lo >= 0) && (i_hi < N_points)) ) - then error_exit(ERROR_EXIT, -"***** interpolate_global: shifted local interpolation needs data\n" -" from outside data array on other side!\n" -" (this should only happen if N_interp > N_points!)\n" -" N_points=%d order=%d N_interp=%d\n" -" x0=%g dx=%g\n" -" x=%g\n" -" ==> i=%d\n" -" ==> shift=%d i_[lo,hi]=[%d,%d]\n" -, - N_points, order, N_interp, - x0, dx, - x, - i, - shift, i_lo, i_hi); /*NOTREACHED*/ - break; - - default: - error_exit(PANIC_EXIT, -"***** interpolate_global: unknown end_action=0x%02x='%c'\n", - (int) end_action, end_action); /*NOTREACHED*/ - } - } - -/* - * do local interpolation - */ -return( - interpolate_local(order, - x0 + i_lo*dx, dx, - y + i_lo, - x) - ); -} /*****************************************************************************/ @@ -269,10 +67,10 @@ return( * y[0...order] = (in) The y data array. * x = (in) The x coordinate for the interpolation. */ -PUBLIC double interpolate_local(int order, - double x0, double dx, - double y[], - double x) +PUBLIC CCTK_REAL interpolate_local(CCTK_INT order, + CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x) { switch (order) { @@ -295,11 +93,6 @@ case 4: case 5: return( interpolate_local_order5(x0, dx, y, x) ); break; - -default: - error_exit(ERROR_EXIT, -"***** interpolate_local: order=%d not implemented yet!\n", - order); /*NOTREACHED*/ } } @@ -309,15 +102,15 @@ default: * This function interpolates a 1st-order (linear) Lagrange polynomial * in a 2-point [0...1] data array centered about the [0.5] grid position. */ -PRIVATE double interpolate_local_order1(double x0, double dx, - double y[], - double x) +PRIVATE CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x) { -double c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1]; -double 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]; -double xc = x0 + 0.5*dx; -double xr = (x - xc) / dx; +CCTK_REAL xc = x0 + 0.5*dx; +CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*c1 ); } @@ -328,16 +121,16 @@ return( c0 + xr*c1 ); * This function interpolates a 2nd-order (quadratic) Lagrange polynomial * in a 3-point [0...2] data array centered about the [1] grid position. */ -PRIVATE double interpolate_local_order2(double x0, double dx, - double y[], - double x) +PRIVATE CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x) { -double c0 = y[1]; -double c1 = (-1.0/2.0)*y[0] + (+1.0/2.0)*y[2]; -double 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]; -double xc = x0 + 1.0*dx; -double xr = (x - xc) / dx; +CCTK_REAL xc = x0 + 1.0*dx; +CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*c2) ); } @@ -348,21 +141,21 @@ return( c0 + xr*(c1 + xr*c2) ); * This function interpolates a 3rd-order (cubic) Lagrange polynomial * in a 4-point [0...3] data array centered about the [1.5] grid position. */ -PRIVATE double interpolate_local_order3(double x0, double dx, - double y[], - double x) +PRIVATE CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x) { -double c0 = (-1.0/16.0)*y[0] + (+9.0/16.0)*y[1] +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]; -double c1 = (+1.0/24.0)*y[0] + (-9.0/ 8.0)*y[1] +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]; -double c2 = (+1.0/ 4.0)*y[0] + (-1.0/ 4.0)*y[1] +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]; -double c3 = (-1.0/ 6.0)*y[0] + (+1.0/ 2.0)*y[1] +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]; -double xc = x0 + 1.5*dx; -double xr = (x - xc) / dx; +CCTK_REAL xc = x0 + 1.5*dx; +CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*(c2 + xr*c3)) ); } @@ -373,22 +166,22 @@ return( c0 + xr*(c1 + xr*(c2 + xr*c3)) ); * This function interpolates a 4th-order (quartic) Lagrange polynomial * in a 5-point [0...4] data array centered about the [2] grid position. */ -PRIVATE double interpolate_local_order4(double x0, double dx, - double y[], - double x) +PRIVATE CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x) { -double c0 = y[2]; -double c1 = (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1] +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]; -double c2 = + (-1.0/24.0)*y[0] + (+2.0/ 3.0)*y[1] + (-5.0/ 4.0)*y[2] +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]; -double c3 = + (-1.0/12.0)*y[0] + (+1.0/ 6.0)*y[1] +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]; -double c4 = + (+1.0/24.0)*y[0] + (-1.0/ 6.0)*y[1] + (+1.0/ 4.0)*y[2] +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]; -double xc = x0 + 2.0*dx; -double xr = (x - xc) / dx; +CCTK_REAL xc = x0 + 2.0*dx; +CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4))) ); } @@ -399,751 +192,32 @@ return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4))) ); * This function interpolates a 5th-order (quintic) Lagrange polynomial * in a 6-point [0...5] data array centered about the [2.5] grid position. */ -PRIVATE double interpolate_local_order5(double x0, double dx, - double y[], - double x) +PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, + CCTK_REAL y[], + CCTK_REAL x) { -double c0 = (+ 3.0/256.0)*y[0] + (-25.0/256.0)*y[1] +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]; -double c1 = (- 3.0/640.0)*y[0] + (+25.0/384.0)*y[1] +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]; -double c2 = (- 5.0/ 96.0)*y[0] + (+13.0/ 32.0)*y[1] +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]; -double c3 = (+ 1.0/ 48.0)*y[0] + (-13.0/ 48.0)*y[1] +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]; -double c4 = (+ 1.0/ 48.0)*y[0] + (- 1.0/ 16.0)*y[1] +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]; -double c5 = (- 1.0/120.0)*y[0] + (+ 1.0/ 24.0)*y[1] +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]; -double xc = x0 + 2.5*dx; -double xr = (x - xc) / dx; +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)))) ); } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -/* - * This function does global Lagrange polynomial differentiation on a - * uniform grid. That is, given N_points uniformly spaced data points - * (x0 + i*dx, y[i]) , i = 0...N_points-1, and given a specified x - * coordinate, it locates a subarray of order+1 data points centered - * around x , interpolates a degree-order (Lagrange) polynomial through - * them, differentiates this, and evaluates the derivative at x . - * - * 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: - * N_points = (in) The number of data points. - * order = (in) The order of the local interpolation, i.e. the degree - * of the interpolated polynomial. - * end_action = (in) Specifies what action to take if the (centered) - * local interpolant would need data from outside - * the domain of the data points. - * 'e' ==> Do an error_exit() . - * 'o' ==> Off-center the local interpolant as - * necessary to keep it inside the domain - * of the data points. - * x0 = (in) The x coordinate corresponding to y[0]. - * dx = (in) The x spacing between the data points. - * y[0...N_points) = (in) The y data array. - * x = (in) The x coordinate for the interpolation. - */ -PUBLIC double interpolate_deriv_global(int N_points, int order, - char end_action, - double x0, double dx, - double y[], - double x) -{ -int N_interp, i, i_lo, i_hi, shift; - -/* - * locate subarray for local differentiation - */ -N_interp = order + 1; -i = floor_double_to_int( (x - x0)/dx ); -subarray_posn(N_interp, i, &i_lo, &i_hi); -if (! ((i_lo >= 0) && (i_hi < N_points)) ) - then { - /* centered local differentiation would need data */ - /* from outside y[] array */ - switch (end_action) - { - case 'e': - error_exit(ERROR_EXIT, -"***** interpolate_deriv_global: local differentiation needs data\n" -" from outside data array!\n" -" N_points=%d order=%d N_interp=%d\n" -" x0=%g dx=%g\n" -" x=%g\n" -" ==> i=%d i_[lo,hi]=[%d,%d]\n" -, - N_points, order, N_interp, - x0, dx, - x, - i, i_lo, i_hi); /*NOTREACHED*/ - - case 'o': - /* off-center the local differentiation as necessary */ - /* to keep it within the y[] array */ - shift = 0; - if (i_lo < 0) then shift = 0 - i_lo; - if (i_hi >= N_points) then shift = (N_points-1) - i_hi; - - i_lo += shift; - i_hi += shift; - - if (! ((i_lo >= 0) && (i_hi < N_points)) ) - then error_exit(ERROR_EXIT, -"***** interpolate_deriv_global: shifted local differentiation needs data\n" -" from outside data array on other side!\n" -" (this should only happen\n" -" if N_interp > N_points!)\n" -" N_points=%d order=%d N_interp=%d\n" -" x0=%g dx=%g\n" -" x=%g\n" -" ==> i=%d\n" -" ==> shift=%d i_[lo,hi]=[%d,%d]\n" -, - N_points, order, N_interp, - x0, dx, - x, - i, - shift, i_lo, i_hi); /*NOTREACHED*/ - break; - - default: - error_exit(PANIC_EXIT, -"***** interpolate_deriv_global: unknown end_action=0x%02x='%c'\n", - (int) end_action, end_action); /*NOTREACHED*/ - } - } - -/* - * do local differentiation - */ -return( - interpolate_deriv_local(order, - x0 + i_lo*dx, dx, - y + i_lo, - x) - ); -} - -/*****************************************************************************/ - -/* - * This function does local Lagrange polynomial differentiation 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, differentiates it, and evaluates - * the derivative at a specified x coordinate. In general the error - * in this computed function value is O(h^(order)) . - * - * 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. - * - * 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. - */ -PUBLIC double interpolate_deriv_local(int order, - double x0, double dx, - double y[], - double x) -{ -switch (order) - { -case 1: - return( interpolate_deriv_local_order1(x0, dx, y, x) ); - break; - -case 2: - return( interpolate_deriv_local_order2(x0, dx, y, x) ); - break; - -case 3: - return( interpolate_deriv_local_order3(x0, dx, y, x) ); - break; - -case 4: - return( interpolate_deriv_local_order4(x0, dx, y, x) ); - break; - -case 5: - return( interpolate_deriv_local_order5(x0, dx, y, x) ); - break; - -default: - error_exit(ERROR_EXIT, -"***** interpolate_deriv_local: order=%d not implemented yet!\n", - order); /*NOTREACHED*/ - } -} - -/*****************************************************************************/ - -/* - * This function interpolates a 1st-order (linear) Lagrange polynomial - * in a 2-point [0...1] data array centered about the [0.5] grid position, - * then evaluates the (0th-order) (constant) derivative of that polynomial. - */ -PRIVATE double interpolate_deriv_local_order1(double x0, double dx, - double y[], - double x) -{ -double c0 = y[1] - y[0]; - -return( c0 / dx ); -} - -/*****************************************************************************/ - -/* - * This function interpolates a 2nd-order (quadratic) Lagrange polynomial - * in a 3-point [0...2] data array centered about the [1] grid position, - * then evaluates the (1st-order) (linear) derivative of that polynomial. - */ -PRIVATE double interpolate_deriv_local_order2(double x0, double dx, - double y[], - double x) -{ -double c0 = (-1.0/2.0)*y[0] + (+1.0/2.0)*y[2]; -double c1 = y[0] - 2.0*y[1] + y[2]; - -double xc = x0 + 1.0*dx; -double xr = (x - xc) / dx; - -return( (c0 + xr*c1) / dx ); -} - -/*****************************************************************************/ - -/* - * This function interpolates a 3rd-order (cubic) Lagrange polynomial - * in a 4-point [0...3] data array centered about the [1.5] grid position, - * then evaluates the (2nd-order) (quadratic) derivative of that polynomial. - */ -PRIVATE double interpolate_deriv_local_order3(double x0, double dx, - double y[], - double x) -{ -double c0 = (+1.0/24.0)*y[0] + (-9.0/ 8.0)*y[1] - + (+9.0/ 8.0)*y[2] + (-1.0/24.0)*y[3]; -double c1 = (+1.0/ 2.0)*y[0] + (-1.0/ 2.0)*y[1] - + (-1.0/ 2.0)*y[2] + (+1.0/ 2.0)*y[3]; -double c2 = (-1.0/ 2.0)*y[0] + (+3.0/ 2.0)*y[1] - + (-3.0/ 2.0)*y[2] + ( 1.0/ 2.0)*y[3]; - -double xc = x0 + 1.5*dx; -double xr = (x - xc) / dx; - -return( (c0 + xr*(c1 + xr*c2)) / dx ); -} - -/*****************************************************************************/ - -/* - * This function interpolates a 4th-order (quartic) Lagrange polynomial - * in a 5-point [0...4] data array centered about the [2] grid position, - * then evaluates the (3rd-order) (cubic) derivative of that polynomial. - */ -PRIVATE double interpolate_deriv_local_order4(double x0, double dx, - double y[], - double x) -{ -double c0 = + (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1] - + (+2.0/ 3.0)*y[3] + (-1.0/12.0)*y[4]; -double c1 = + (-1.0/12.0)*y[0] + (+4.0/ 3.0)*y[1] + (-5.0/ 2.0)*y[2] - + (+4.0/ 3.0)*y[3] + (-1.0/12.0)*y[4]; -double c2 = + (-1.0/ 4.0)*y[0] + (+1.0/ 2.0)*y[1] - + (-1.0/ 2.0)*y[3] + (+1.0/ 4.0)*y[4]; -double c3 = + (+1.0/ 6.0)*y[0] + (-2.0/ 3.0)*y[1] + y[2] - + (-2.0/ 3.0)*y[3] + (+1.0/ 6.0)*y[4]; - -double xc = x0 + 2.0*dx; -double xr = (x - xc) / dx; - -return( (c0 + xr*(c1 + xr*(c2 + xr*c3))) / dx ); -} - -/*****************************************************************************/ - -/* - * This function interpolates a 5th-order (quartic) Lagrange polynomial - * in a 6-point [0...5] data array centered about the [2.5] grid position, - * then evaluates the (4th-order) (quartic) derivative of that polynomial. - */ -PRIVATE double interpolate_deriv_local_order5(double x0, double dx, - double y[], - double x) -{ -double c0 = (- 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]; -double c1 = (- 5.0/ 48.0)*y[0] + (+13.0/ 16.0)*y[1] - + (-17.0/ 24.0)*y[2] + (-17.0/ 24.0)*y[3] - + (+13.0/ 16.0)*y[4] + (- 5.0/ 48.0)*y[5]; -double c2 = (+ 1.0/ 16.0)*y[0] + (-13.0/ 16.0)*y[1] - + (+17.0/ 8.0)*y[2] + (-17.0/ 8.0)*y[3] - + (+13.0/ 16.0)*y[4] + (- 1.0/ 16.0)*y[5]; -double c3 = (+ 1.0/ 12.0)*y[0] + (- 1.0/ 4.0)*y[1] - + (+ 1.0/ 6.0)*y[2] + (+ 1.0/ 6.0)*y[3] - + (- 1.0/ 4.0)*y[4] + (+ 1.0/ 12.0)*y[5]; -double c4 = (- 1.0/ 24.0)*y[0] + (+ 5.0/ 24.0)*y[1] - + (- 5.0/ 12.0)*y[2] + (+ 5.0/ 12.0)*y[3] - + (- 5.0/ 24.0)*y[4] + (+ 1.0/ 24.0)*y[5]; - -double xc = x0 + 2.5*dx; -double xr = (x - xc) / dx; - -return( (c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4)))) / dx ); -} - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -/* - * This function does global Lagrange polynomial interpolation on a - * nonuniform grid. That is, given N_points uniformly spaced data - * points (x0 + i*dx, y[i]) , i = 0...N_points-1, and given a specified - * x coordinate, it locates a subarray of order+1 data points centered - * around x , interpolates a degree-order (Lagrange) polynomial through - * them, and evaluates this at x . - * - * The local interpolation is done in the Lagrange form, which is fairly - * well-conditioned as ways of writing the interpolating polynomial go. - * Even so, the interpolation is still highly ill-conditioned for small - * dx and/or large order and/or x outside the domain of the data - * points. - * - * Arguments: - * N_points = (in) The number of data points. - * order = (in) The order of the local interpolation, i.e. the degree - * of the interpolated polynomial. - * end_action = (in) Specifies what action to take if the (centered) - * local interpolant would need data from outside - * the domain of the data points. - * 'e' ==> Do an error_exit() . - * 'o' ==> Off-center the local interpolant as - * necessary to keep it inside the domain - * of the data points. - * xx[0...N_points) = (in) The x data array. This is assumed to be in - * (increasing) order. - * yy[0...N_points) = (in) The y data array. - * x = (in) The x coordinate for the interpolation. - */ -PUBLIC double interpolate_nu_global(int N_points, int order, - char end_action, - double xx[], - double yy[], - double x) -{ -/* - * Implementation note: - * - * We assume this function will typically be used to interpolate the - * same xx[] and yy[] values to a number of different x values, - * and that these x values will typically change monotonically. - * Therefore, we use a sequential search to locate x in the xx[] - * array, starting from a cached position and going up or down as - * appropriate. - */ -static int cache_valid = FALSE; -static int cache_N_points; /* cache valid only if N_points == this */ -static double *cache_xx; /* && xx == this */ -static int cache_i; - -int N_interp, i, i_lo, i_hi, shift; -boolean ok; - -if (cache_valid && (N_points == cache_N_points) && (xx == cache_xx)) - then i = cache_i; /* use cache */ - else i = 0; /* default if cache isn't available */ - -/* - * Find i such that (roughly speaking) xx[i] <= x < xx[i+1] : - * - * More precisely, there are 3 cases to consider: - * x < xx[0] ==> i = 0 - * xx[0] <= x <= xx[N_points-1] ==> i >= 0 && i+1 < N_points - * && xx[i] <= x < xx[i+1] - * x > xx[N_points-1] ==> i = N_points-1 - */ -if (x >= xx[i]) - then { - /* search up from starting position */ - for (++i ; ((i < N_points) && (x >= xx[i])) ; ++i) - { - } - --i; - } - else { - /* search down from starting position */ - for (--i ; ((i >= 0) && (x < xx[i])) ; --i) - { - } - ++i; - } - -/* - * firewall: sanity-check that we found i correctly - */ -if ((x >= xx[0]) && (x <= xx[N_points-1])) - then ok = (i >= 0) && (i+1 < N_points) && (x >= xx[i]) && (x < xx[i+1]); -else if (x < xx[0]) - then ok = (i == 0); -else if (x > xx[N_points-1]) - then ok = (i == N_points-1); -else ok = FALSE; -if (! ok) - then error_exit(PANIC_EXIT, -"***** interpolate_nu_global: search result failed consistency check!\n" -" (this should never happen!)\n" -" N_points=%d\n" -" xx[0]=%g xx[N_points-1]=%g\n" -" x=%g\n" -" ==> i=%d\n" -" xx[i]=%g xx[i+1]=%g\n" -, - N_points, - xx[0], xx[N_points-1], - x, - i, - xx[i], xx[i+1]); /*NOTREACHED*/ - -/* - * update cache for future use - */ -cache_valid = TRUE; -cache_N_points = N_points; -cache_xx = xx; -cache_i = i; - -/* - * locate subarray - */ -N_interp = order + 1; -subarray_posn(N_interp, i, &i_lo, &i_hi); -if (! ((i_lo >= 0) && (i_hi < N_points)) ) - then { - /* centered local interpolation would need data */ - /* from outside y[] array */ - switch (end_action) - { - case 'e': - error_exit(ERROR_EXIT, -"***** interpolate_nu_global: local interpolation needs data\n" -" from outside data array!\n" -" N_points=%d order=%d N_interp=%d\n" -" xx[0]=%g xx[N_points-1]=%g\n" -" x=%g\n" -" ==> i=%d i_[lo,hi]=[%d,%d]\n" -, - N_points, order, N_interp, - xx[0], xx[N_points-1], - x, - i, i_lo, i_hi); /*NOTREACHED*/ - - case 'o': - /* off-center the local interpolation as necessary */ - /* to keep it within the y[] array */ - shift = 0; - if (i_lo < 0) then shift = 0 - i_lo; - if (i_hi >= N_points) then shift = (N_points-1) - i_hi; - - i_lo += shift; - i_hi += shift; - - if (! ((i_lo >= 0) && (i_hi < N_points)) ) - then error_exit(ERROR_EXIT, -"***** interpolate_nu_global: shifted local interpolation needs data\n" -" from outside data array on other side!\n" -" (this should only happen\n" -" if N_interp > N_points!)\n" -" N_points=%d order=%d N_interp=%d\n" -" xx[0]=%g xx[N_points-1]=%g\n" -" x=%g\n" -" ==> i=%d\n" -" ==> shift=%d i_[lo,hi]=[%d,%d]\n" -, - N_points, order, N_interp, - xx[0], xx[N_points-1], - x, - i, - shift, i_lo, i_hi); /*NOTREACHED*/ - break; - - default: - error_exit(PANIC_EXIT, -"***** interpolate_nu_global: unknown end_action=0x%02x='%c'\n", - (int) end_action, end_action); /*NOTREACHED*/ - } - } - -/* - * do local interpolation - */ -return( - interpolate_nu_local(order, - xx + i_lo, - yy + i_lo, - x) - ); -} - -/*****************************************************************************/ - -/* - * This function does local Lagrange polynomial interpolation on a - * nonuniform grid. That is, given order+1 data points (xx[i],yy[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) . - * - * The interpolation is done in the Lagrange form, which is fairly - * well-conditioned as ways of writing the interpolating polynomial go. - * Even so, the interpolation is still highly ill-conditioned for small - * dx and/or large order and/or x outside the domain of the data - * points. - * - * Arguments: - * N_points = (in) The number of data points. - * order = (in) The order of the local interpolation, i.e. the degree - * of the interpolated polynomial. - * xx[0...N_points) = (in) The x data array. This is assumed to be in - * (increasing) order. - * yy[0...N_points) = (in) The y data array. - * x = (in) The x coordinate for the interpolation. - */ -PUBLIC double interpolate_nu_local(int order, - double xx[], - double yy[], - double x) -{ -switch (order) - { -case 3: - return( interpolate_nu_local_order3(xx, yy, x) ); - break; - -case 4: - return( interpolate_nu_local_order4(xx, yy, x) ); - break; - -case 5: - return( interpolate_nu_local_order5(xx, yy, x) ); - break; - -default: - error_exit(ERROR_EXIT, -"***** interpolate_nu_local: order=%d not implemented yet!\n", - order); /*NOTREACHED*/ - } -} - -/*****************************************************************************/ - -/* - * This function does 4-point local Lagrange polynomial interpolation - * on a nonuniform grid, interpolating a 3rd-order polynomial. - */ -PRIVATE double interpolate_nu_local_order3(double xx[], - double yy[], - double x) -{ -double num0 = ( x - xx[1]) * ( x - xx[2]) * ( x - xx[3]); -double den0 = (xx[0] - xx[1]) * (xx[0] - xx[2]) * (xx[0] - xx[3]); -double c0 = num0 / den0; - -double num1 = ( x - xx[0]) * ( x - xx[2]) * ( x - xx[3]); -double den1 = (xx[1] - xx[0]) * (xx[1] - xx[2]) * (xx[1] - xx[3]); -double c1 = num1 / den1; - -double num2 = ( x - xx[0]) * ( x - xx[1]) * ( x - xx[3]); -double den2 = (xx[2] - xx[0]) * (xx[2] - xx[1]) * (xx[2] - xx[3]); -double c2 = num2 / den2; - -double num3 = ( x - xx[0]) * ( x - xx[1]) * ( x - xx[2]); -double den3 = (xx[3] - xx[0]) * (xx[3] - xx[1]) * (xx[3] - xx[2]); -double c3 = num3 / den3; - -return( c0*yy[0] + c1*yy[1] + c2*yy[2] + c3*yy[3] ); -} - -/*****************************************************************************/ - -/* - * This function does 5-point local Lagrange polynomial interpolation - * on a nonuniform grid, interpolating a 4th-order polynomial. - */ -PRIVATE double interpolate_nu_local_order4(double xx[], - double yy[], - double x) -{ -double num0 = ( x -xx[1]) * ( x -xx[2]) * ( x -xx[3]) * ( x -xx[4]); -double den0 = (xx[0]-xx[1]) * (xx[0]-xx[2]) * (xx[0]-xx[3]) * (xx[0]-xx[4]); -double c0 = num0 / den0; - -double num1 = ( x -xx[0]) * ( x -xx[2]) * ( x -xx[3]) * ( x -xx[4]); -double den1 = (xx[1]-xx[0]) * (xx[1]-xx[2]) * (xx[1]-xx[3]) * (xx[1]-xx[4]); -double c1 = num1 / den1; - -double num2 = ( x -xx[0]) * ( x -xx[1]) * ( x -xx[3]) * ( x -xx[4]); -double den2 = (xx[2]-xx[0]) * (xx[2]-xx[1]) * (xx[2]-xx[3]) * (xx[2]-xx[4]); -double c2 = num2 / den2; - -double num3 = ( x -xx[0]) * ( x -xx[1]) * ( x -xx[2]) * ( x -xx[4]); -double den3 = (xx[3]-xx[0]) * (xx[3]-xx[1]) * (xx[3]-xx[2]) * (xx[3]-xx[4]); -double c3 = num3 / den3; - -double num4 = ( x -xx[0]) * ( x -xx[1]) * ( x -xx[2]) * ( x -xx[3]); -double den4 = (xx[4]-xx[0]) * (xx[4]-xx[1]) * (xx[4]-xx[2]) * (xx[4]-xx[3]); -double c4 = num4 / den4; - -return( - c0*yy[0] + c1*yy[1] + c2*yy[2] + c3*yy[3] + c4*yy[4] - ); -} - -/*****************************************************************************/ - -/* - * This function does 6-point local Lagrange polynomial interpolation - * on a nonuniform grid, interpolating a 5th-order polynomial. - */ -PRIVATE double interpolate_nu_local_order5(double xx[], - double yy[], - double x) -{ -double num0 - = ( x -xx[1])*( x -xx[2])*( x -xx[3])*( x -xx[4])*( x -xx[5]); -double den0 - = (xx[0]-xx[1])*(xx[0]-xx[2])*(xx[0]-xx[3])*(xx[0]-xx[4])*(xx[0]-xx[5]); -double c0 = num0 / den0; - -double num1 - = ( x -xx[0])*( x -xx[2])*( x -xx[3])*( x -xx[4])*( x -xx[5]); -double den1 - = (xx[1]-xx[0])*(xx[1]-xx[2])*(xx[1]-xx[3])*(xx[1]-xx[4])*(xx[1]-xx[5]); -double c1 = num1 / den1; - -double num2 - = ( x -xx[0])*( x -xx[1])*( x -xx[3])*( x -xx[4])*( x -xx[5]); -double den2 - = (xx[2]-xx[0])*(xx[2]-xx[1])*(xx[2]-xx[3])*(xx[2]-xx[4])*(xx[2]-xx[5]); -double c2 = num2 / den2; - -double num3 - = ( x -xx[0])*( x -xx[1])*( x -xx[2])*( x -xx[4])*( x -xx[5]); -double den3 - = (xx[3]-xx[0])*(xx[3]-xx[1])*(xx[3]-xx[2])*(xx[3]-xx[4])*(xx[3]-xx[5]); -double c3 = num3 / den3; - -double num4 - = ( x -xx[0])*( x -xx[1])*( x -xx[2])*( x -xx[3])*( x -xx[5]); -double den4 - = (xx[4]-xx[0])*(xx[4]-xx[1])*(xx[4]-xx[2])*(xx[4]-xx[3])*(xx[4]-xx[5]); -double c4 = num4 / den4; - -double num5 - = ( x -xx[0])*( x -xx[1])*( x -xx[2])*( x -xx[3])*( x -xx[4]); -double den5 - = (xx[5]-xx[0])*(xx[5]-xx[1])*(xx[5]-xx[2])*(xx[5]-xx[3])*(xx[5]-xx[4]); -double c5 = num5 / den5; - -return( - c0*yy[0] + c1*yy[1] + c2*yy[2] + c3*yy[3] + c4*yy[4] + c5 * yy[5] - ); -} - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -/* - * This function computes subarray positions for an interpolation or - * similar operation. That is, it encapsulates the tricky odd/even - * computations needed to find the exact boundaries of the subarray in - * which an interpolation or similar operation should be done. - * - * Suppose we have a data array [i_min ... i_max] from which we wish to - * select an N point subarray for interpolation, and the subarray is to - * be centered between [i] and [i+1]. Then if N is even, the centering - * is unique, but if N is odd, we have to choose either the i-greater or - * i-less side on which to add the last point. We choose the former here, - * as would be appropriate for functions varying most rapidly at smaller - * [i] and less rapidly at larger [i]. Thus, for example, we have - * - * center - * here - * N i_lo i_hi [i-2] [i-1] [i+0] | [i+1] [i+2] [i+3] [i+4] - * - ---- ---- ----- ----- ----- | ----- ----- ----- ----- - * 2 i-0 i+1 x | x - * 3 i-0 i+2 x | x x - * 4 i-1 i+2 x x | x x - * 5 i-1 i+3 x x | x x x - * 6 i-2 i+3 x x x | x x x - * 7 i-2 i+4 x x x | x x x x - * - * The purpose of this function is to compute the bounds of the subarray, - * i_lo and i_hi . This function intentionally doesn't check for these - * overflowing the data array, because our caller can presumably provide - * a more useful error message if this occurs. - * - * Arguments: - * N = (in) The number of points in the subarray. - * i = (in) The subarray is to be centered between [i] and [i+1]. - * pi_{lo,hi} --> (out) The inclusive {lower,upper} bound of the subarray. - * Either or both of these arguments may be NULL, in - * which case the corresponding result isn't stored. - * - * Results: - * The function returns i_lo, the inclusive lower bound of the subarray. - */ -PUBLIC int subarray_posn(int N, int i, int *pi_lo, int *pi_hi) -{ -int i_lo = i + 1 - (N >> 1); -int i_hi = i + ((N+1) >> 1); - -if (HOW_MANY_IN_RANGE(i_lo,i_hi) != N) - then error_exit(PANIC_EXIT, -"***** interpolate_posn: subarray has wrong number of points!\n" -" (this should never happen!)\n" -" N=%d i=%d i_lo=%d i_hi=%d\n" -" HOW_MANY_IN_RANGE(i_lo,i_hi)=%d\n" -, - N, i, i_lo, i_hi, - HOW_MANY_IN_RANGE(i_lo,i_hi)); /*NOTREACHED*/ - -if (pi_lo != NULL) - then *pi_lo = i_lo; -if (pi_hi != NULL) - then *pi_hi = i_hi; - -return(i_lo); -} diff --git a/src/make.code.defn b/src/make.code.defn index 0483ce6..f720acb 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,7 +2,7 @@ # $Header$ # Source files in this directory -SRCS = +SRCS = CheckParameters.c Cartoon2DBC.c interpolate.c # Subdirectories containing source files SUBDIRS = |