From ff6dcc8f76dce5197014401eba240c860af6eacf Mon Sep 17 00:00:00 2001 From: schnetter Date: Fri, 4 Sep 2009 23:47:56 +0000 Subject: Correct errors. Add extrapolator to set up the boundary conditions Gamma for initial data. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/NewRad/trunk@3 16f5cafb-89ad-4b9f-9d0b-9d7a0a37d03b --- interface.ccl | 9 ++- src/extrap.cc | 186 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 2 +- src/newrad.cc | 165 +++++++++++++++++++++++++++++------------------ 4 files changed, 296 insertions(+), 66 deletions(-) create mode 100644 src/extrap.cc diff --git a/interface.ccl b/interface.ccl index 57f6416..234262d 100644 --- a/interface.ccl +++ b/interface.ccl @@ -4,6 +4,12 @@ IMPLEMENTS: NewRad USES INCLUDE HEADER: GenericFD.h +CCTK_INT FUNCTION \ + ExtrapolateGammas \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_REAL ARRAY INOUT var) +PROVIDES FUNCTION ExtrapolateGammas WITH ExtrapolateGammas1 LANGUAGE C + CCTK_INT FUNCTION \ NewRad_Apply \ (CCTK_POINTER_TO_CONST IN cctkGH, \ @@ -11,6 +17,5 @@ CCTK_INT FUNCTION \ CCTK_REAL ARRAY INOUT rhs, \ CCTK_REAL IN var0, \ CCTK_REAL IN v0, \ - CCTK_REAL IN radpower, \ - CCTK_INT IN width) + CCTK_INT IN radpower) PROVIDES FUNCTION NewRad_Apply WITH NewRad_Apply1 LANGUAGE C diff --git a/src/extrap.cc b/src/extrap.cc new file mode 100644 index 0000000..b752aa0 --- /dev/null +++ b/src/extrap.cc @@ -0,0 +1,186 @@ +#include +#include + +#include + +#ifdef CCTK_CXX_RESTRICT +# undef restrict +# define restrict CCTK_CXX_RESTRICT +#endif + +#define KRANC_C +extern "C" { +#include +} + +using namespace std; + + + +// Adapted from BSSN_MoL's files Init.F +static +void extrap_kernel (cGH const* restrict const cctkGH, + int const* restrict const bmin, + int const* restrict const bmax, + int const* restrict const dir, + CCTK_REAL* restrict const var) +{ + int const ni = cctkGH->cctk_lsh[0]; + int const nj = cctkGH->cctk_lsh[1]; + int const nk = cctkGH->cctk_lsh[2]; + + int const si = -dir[0]; + int const sj = -dir[1]; + int const sk = -dir[2]; + + int const di = 1; + int const dj = ni; + int const dk = ni*nj; + int const np = ni*nj*nk; + + int const dind = si*di + sj*dj + sk*dk; + + int imin[3], imax[3], idir[3]; + for (int d=0; d<3; ++d) { + if (dir[d]<0) { + // lower boundary + imin[d] = bmax[d]-1; + imax[d] = bmin[d]-1; + idir[d] = -1; + } else { + // interior and upper boundary + imin[d] = bmin[d]; + imax[d] = bmax[d]; + idir[d] = +1; + } + } + + for (int k=imin[2]; k!=imax[2]; k+=idir[2]) { + for (int j=imin[1]; j!=imax[1]; j+=idir[1]) { + for (int i=imin[0]; i!=imax[0]; i+=idir[0]) { + int const ind = CCTK_GFINDEX3D(cctkGH, i,j,k); + + // Test looping directions + if (i==0) assert (idir[0]<0); + if (j==0) assert (idir[1]<0); + if (k==0) assert (idir[2]<0); + if (i==ni-1) assert (idir[0]>0); + if (j==nj-1) assert (idir[1]>0); + if (k==nk-1) assert (idir[2]>0); + + // Apply boundary conditions to get e.g. Gammas on physical + // boundaries. Notice that we only want to apply boundary + // conditions to the Gammas, all other variables have been + // defined point-wise all the way to the boundaries from the + // original ADM quantities (which SHOULD be correctly defined + // all the way to the boundaries by the initial data + // routines). Here I use cubic extrapolation (though rational + // extrapolation might be better). + + assert (ind >=0 and ind =0 and ind+4*dindcctk_lsh, cctkGH->cctk_lssh, cctkGH->cctk_bbox, + cctkGH->cctk_nghostzones, + imin, imax, is_symbnd, is_physbnd, is_ipbnd); + + // Loop over all faces: + // Loop over faces first, then corners, athen edges, so that the + // stencil only sees points that have already been treated + // ifec means: interior-face-edge-corner + for (int ifec=1; ifec<=3; ++ifec) { + for (int dir2=-1; dir2<=+1; ++dir2) { + for (int dir1=-1; dir1<=+1; ++dir1) { + for (int dir0=-1; dir0<=+1; ++dir0) { + int const dir[3] = { dir0, dir1, dir2 }; + + int nnz = 0; + for (int d=0; d<3; ++d) { + if (dir[d]) ++nnz; + } + if (nnz == ifec) { + + // one of tahe faces is a boundary + bool have_bnd = false; + // all boundary faces are physical boundaries + bool all_physbnd = true; + + int bmin[3], bmax[3]; + for (int d=0; d<3; ++d) { + switch (dir[d]) { + case -1: + bmin[d] = 0; + bmax[d] = imin[d]; + have_bnd = true; + all_physbnd = all_physbnd and is_physbnd[2*d+0]; + break; + case 0: + bmin[d] = imin[d]; + bmax[d] = imax[d]; + break; + case +1: + bmin[d] = imax[d]; + bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)]; + have_bnd = true; + all_physbnd = all_physbnd and is_physbnd[2*d+1]; + break; + } + } + + if (have_bnd and all_physbnd) { + extrap_kernel (cctkGH, bmin, bmax, dir, var); + } + + } + } // for dir0 dir1 dir2 + } + } + } +} + + + +extern "C" +CCTK_INT ExtrapolateGammas1 (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_REAL* restrict const var) +{ + cGH const* restrict const cctkGH = static_cast (cctkGH_); + if (not cctkGH) { + CCTK_WARN (CCTK_WARN_ABORT, + "cctkGH is NULL"); + } + +#if 0 + CCTK_REAL* restrict const var = CCTK_VarDataPtr (cctkGH, 0, varname); + if (not var) { + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot access variable \"%s\"", varname); + } +#endif + + if (not var) { + CCTK_WARN (CCTK_WARN_ABORT, + "Pointer to variable is NULL"); + } + + extrap_loop (cctkGH, var); + + return 0; +} diff --git a/src/make.code.defn b/src/make.code.defn index 244b721..40d1780 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn NewRad # Source files in this directory -SRCS = newrad.cc +SRCS = extrap.cc newrad.cc # Subdirectories containing source files SUBDIRS = diff --git a/src/newrad.cc b/src/newrad.cc index a1edcf0..3845075 100644 --- a/src/newrad.cc +++ b/src/newrad.cc @@ -1,4 +1,5 @@ -#include +#include +#include #include @@ -12,6 +13,8 @@ extern "C" { #include } +using namespace std; + // Adapted from BSSN_MoL's files NewRad.F and newrad.h @@ -28,15 +31,15 @@ void newrad_kernel (cGH const* restrict const cctkGH, CCTK_REAL const* restrict const r, CCTK_REAL const& var0, CCTK_REAL const& v0, - CCTK_REAL const& radpower) + int const radpower) { int const ni = cctkGH->cctk_lsh[0]; int const nj = cctkGH->cctk_lsh[1]; int const nk = cctkGH->cctk_lsh[2]; - int const si = -dir[0]; - int const sj = -dir[1]; - int const sk = -dir[2]; + int const si = dir[0]; + int const sj = dir[1]; + int const sk = dir[2]; int const di = 1; int const dj = ni; @@ -49,11 +52,34 @@ void newrad_kernel (cGH const* restrict const cctkGH, CCTK_REAL const idy = 1.0/dy; CCTK_REAL const idz = 1.0/dz; - for (int k=bmin[2]; k0); + if (j==nj-1) assert (idir[1]>0); + if (k==nk-1) assert (idir[2]>0); + { // The main part of the boundary condition assumes that we // have an outgoing radial wave with some speed v0: @@ -74,13 +100,13 @@ void newrad_kernel (cGH const* restrict const cctkGH, CCTK_REAL const vy = v0*y[ind]*rpi; CCTK_REAL const vz = v0*z[ind]*rpi; - int const svi = i==0 ? +1 : i==ni-1 ? -1 : vx>0 ? +1 : -1; - int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1; - int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1; + int const svi = dir[0]; + int const svj = dir[1]; + int const svk = dir[2]; // Find x derivative CCTK_REAL derivx; - if (i>0 and i0 and j0 and k=0 and ip=0 and jp=0 and kp0 ? +1 : -1; - int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1; - int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1; + int const svi = dir[0]; + int const svj = dir[1]; + int const svk = dir[2]; // Find x derivative CCTK_REAL derivx; - if (ip>0 and ip0 and jp0 and kpcctk_nghostzones, imin, imax, is_symbnd, is_physbnd, is_ipbnd); - // Loop over all faces - for (int dir2=-1; dir2<=+1; ++dir2) { - for (int dir1=-1; dir1<=+1; ++dir1) { - for (int dir0=-1; dir0<=+1; ++dir0) { - int const dir[3] = { dir0, dir1, dir2 }; - - // one of tahe faces is a boundary - bool have_bnd = false; - // all boundary faces are physical boundaries - bool all_physbnd = true; - - int bmin[3], bmax[3]; - for (int d=0; d<3; ++d) { - switch (dir[d]) { - case -1: - bmin[d] = 0; - bmax[d] = imin[d]; - have_bnd = true; - all_physbnd = all_physbnd and is_physbnd[2*d+0]; - break; - case 0: - bmin[d] = imin[d]; - bmax[d] = imax[d]; - break; - case +1: - bmin[d] = imax[d]; - bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)]; - have_bnd = true; - all_physbnd = all_physbnd and is_physbnd[2*d+1]; - break; + // Loop over all faces: + // Loop over faces first, then corners, athen edges, so that the + // stencil only sees points that have already been treated + // ifec means: interior-face-edge-corner + for (int ifec=1; ifec<=3; ++ifec) { + for (int dir2=-1; dir2<=+1; ++dir2) { + for (int dir1=-1; dir1<=+1; ++dir1) { + for (int dir0=-1; dir0<=+1; ++dir0) { + int const dir[3] = { dir0, dir1, dir2 }; + + int nnz = 0; + for (int d=0; d<3; ++d) { + if (dir[d]) ++nnz; } - } - - if (have_bnd and all_physbnd) { - newrad_kernel (cctkGH, bmin, bmax, dir, - var, rhs, x,y,z,r, var0, v0, radpower); - } - - } // for dir0 dir1 dir2 + if (nnz == ifec) { + + // one of tahe faces is a boundary + bool have_bnd = false; + // all boundary faces are physical boundaries + bool all_physbnd = true; + + int bmin[3], bmax[3]; + for (int d=0; d<3; ++d) { + switch (dir[d]) { + case -1: + bmin[d] = 0; + bmax[d] = imin[d]; + have_bnd = true; + all_physbnd = all_physbnd and is_physbnd[2*d+0]; + break; + case 0: + bmin[d] = imin[d]; + bmax[d] = imax[d]; + break; + case +1: + bmin[d] = imax[d]; + bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)]; + have_bnd = true; + all_physbnd = all_physbnd and is_physbnd[2*d+1]; + break; + } + } + + if (have_bnd and all_physbnd) { + newrad_kernel (cctkGH, bmin, bmax, dir, + var, rhs, x,y,z,r, var0, v0, radpower); + } + + } + } // for dir0 dir1 dir2 + } } } } @@ -261,8 +301,7 @@ CCTK_INT NewRad_Apply1 (CCTK_POINTER_TO_CONST const cctkGH_, CCTK_REAL * restrict const rhs, CCTK_REAL const var0, CCTK_REAL const v0, - CCTK_REAL const radpower, - CCTK_INT const width) + CCTK_INT const radpower) { cGH const* restrict const cctkGH = static_cast (cctkGH_); if (not cctkGH) { @@ -305,7 +344,7 @@ CCTK_INT NewRad_Apply1 (CCTK_POINTER_TO_CONST const cctkGH_, "Cannot access coordinate variables x, y, z, and r"); } - newrad_loop (cctkGH, var, rhs, x,y,z,r, var0, v0, radpower, width); + newrad_loop (cctkGH, var, rhs, x,y,z,r, var0, v0, radpower); return 0; } -- cgit v1.2.3