Recreating code after origin /data crash.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@6 eec4d7dc-71c2-46d6-addf-10296150bf52

1 files changed, 229 insertions, 14 deletions

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;
+}