From 8fe9c3a66c1dc44e137e176ff5018f4d8d65eb4d Mon Sep 17 00:00:00 2001 From: sai Date: Fri, 19 Nov 1999 10:25:13 +0000 Subject: Recreating code after origin /data crash. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@6 eec4d7dc-71c2-46d6-addf-10296150bf52 --- src/Cartoon2DBC.c | 243 ++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 229 insertions(+), 14 deletions(-) (limited to 'src/Cartoon2DBC.c') 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 +#include +#include +#include +#include +#include +#include + #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; +} -- cgit v1.2.3