diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/make.code.defn | 7 | ||||
-rw-r--r-- | src/newrad.cc | 311 |
2 files changed, 318 insertions, 0 deletions
diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..244b721 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,7 @@ +# Main make.code.defn file for thorn NewRad + +# Source files in this directory +SRCS = newrad.cc + +# Subdirectories containing source files +SUBDIRS = diff --git a/src/newrad.cc b/src/newrad.cc new file mode 100644 index 0000000..a1edcf0 --- /dev/null +++ b/src/newrad.cc @@ -0,0 +1,311 @@ +#include <math.h> + +#include <cctk.h> + +#ifdef CCTK_CXX_RESTRICT +# undef restrict +# define restrict CCTK_CXX_RESTRICT +#endif + +#define KRANC_C +extern "C" { +#include <GenericFD.h> +} + + + +// Adapted from BSSN_MoL's files NewRad.F and newrad.h +static +void newrad_kernel (cGH const* restrict const cctkGH, + int const* restrict const bmin, + int const* restrict const bmax, + int const* restrict const dir, + CCTK_REAL const* restrict const var, + CCTK_REAL * restrict const rhs, + CCTK_REAL const* restrict const x, + CCTK_REAL const* restrict const y, + CCTK_REAL const* restrict const z, + CCTK_REAL const* restrict const r, + CCTK_REAL const& var0, + CCTK_REAL const& v0, + CCTK_REAL 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 di = 1; + int const dj = ni; + int const dk = ni*nj; + + CCTK_REAL const dx = cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0]; + CCTK_REAL const dy = cctkGH->cctk_delta_space[1] / cctkGH->cctk_levfac[1]; + CCTK_REAL const dz = cctkGH->cctk_delta_space[2] / cctkGH->cctk_levfac[2]; + CCTK_REAL const idx = 1.0/dx; + CCTK_REAL const idy = 1.0/dy; + CCTK_REAL const idz = 1.0/dz; + + for (int k=bmin[2]; k<bmax[2]; ++k) { + for (int j=bmin[1]; j<bmax[1]; ++j) { + for (int i=bmin[0]; i<bmax[0]; ++i) { + int const ind = CCTK_GFINDEX3D(cctkGH, i,j,k); + + { + // The main part of the boundary condition assumes that we + // have an outgoing radial wave with some speed v0: + // + // var = var0 + u(r-v0*t)/r + // + // This implies the following differential equation: + // + // d_t var = - v^i d_i var - v0 (var - var0) / r + // + // where vi = v0 xi/r + + // Find local wave speeds + CCTK_REAL const rp = r[ind]; + CCTK_REAL const rpi = 1.0/rp; + + CCTK_REAL const vx = v0*x[ind]*rpi; + 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; + + // Find x derivative + CCTK_REAL derivx; + if (i>0 and i<ni-1) { + derivx = 0.5*(var[ind+di]-var[ind-di])*idx; + } else { + derivx = svi*0.5*(3*var[ind] - 4*var[ind-svi*di] + + var[ind-2*svi*di])*idx; + } + + // Find y derivative + CCTK_REAL derivy; + if (j>0 and j<nj-1) { + derivy = 0.5*(var[ind+dj]-var[ind-dj])*idy; + } else { + derivy = svj*0.5*(3*var[ind] - 4*var[ind-svj*dj] + + var[ind-2*svj*dj])*idy; + } + + // Find z derivative + CCTK_REAL derivz; + if (k>0 and k<nk-1) { + derivz = 0.5*(var[ind+dk]-var[ind-dk])*idz; + } else { + derivz = svk*0.5*(3*var[ind] - 4*var[ind-svk*dk] + + var[ind-2*svk*dk])*idz; + } + + // Calculate source term + rhs[ind] = + - vx*derivx - vy*derivy - vz*derivz - v0*(var[ind] - var0)*rpi; + + } + + if (radpower >= 0) { + // ***************************************** + // *** EXTRAPOLATION OF MISSING PART *** + // ***************************************** + // + // Here we try to extrapolate for the part of the boundary + // that does not behave as a pure wave (i.e. Coulomb type + // terms caused by infall of the coordinate lines). + // + // This we do by comparing the source term one grid point + // away from the boundary (which we already have), to what + // we would have obtained if we had used the boundary + // condition there. The difference gives us an idea of the + // missing part and we extrapolate that to the boundary + // assuming a power-law decay. + + int const ip = i-si; + int const jp = j-sj; + int const kp = k-sk; + + // Find local wave speeds + int const indp = CCTK_GFINDEX3D(cctkGH, ip,jp,kp); + + CCTK_REAL const rp = r[indp]; + CCTK_REAL const rpi = 1.0/rp; + + CCTK_REAL const vx = v0*x[indp]*rpi; + CCTK_REAL const vy = v0*y[indp]*rpi; + CCTK_REAL const vz = v0*z[indp]*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; + + // Find x derivative + CCTK_REAL derivx; + if (ip>0 and ip<ni-1) { + derivx = 0.5*(var[indp+di]-var[indp-di])*idx; + } else { + derivx = svi*0.5*(3*var[indp] - 4*var[indp-svi*di] + + var[indp-2*svi*di])*idx; + } + + // Find y derivative + CCTK_REAL derivy; + if (jp>0 and jp<nj-1) { + derivy = 0.5*(var[indp+dj]-var[indp-dj])*idy; + } else { + derivy = svj*0.5*(3*var[indp] - 4*var[indp-svj*dj] + + var[indp-2*svj*dj])*idy; + } + + // Find z derivative + CCTK_REAL derivz; + if (kp>0 and kp<nk-1) { + derivz = 0.5*(var[indp+dk]-var[indp-dk])*idz; + } else { + derivz = svk*0.5*(3*var[indp] - 4*var[indp-svk*dk] + + var[indp-2*svk*dk])*idz; + } + + // Find difference in sources + CCTK_REAL const aux = + rhs[indp] + + vx*derivx + vy*derivy + vz*derivz + v0*(var[indp] - var0)*rpi; + + // Extrapolate difference and add it to source in boundary + rhs[ind] += aux*pow(rp/r[ind],radpower); + + } // if radpower>=0 + + } // for i j k + } + } +} + + + +// Adapted from Kranc's KrancNumericalTools/GenericFD's file +// GenericFD.c +static +void newrad_loop (cGH const* restrict const cctkGH, + CCTK_REAL const* restrict const var, + CCTK_REAL * restrict const rhs, + CCTK_REAL const* restrict const x, + CCTK_REAL const* restrict const y, + CCTK_REAL const* restrict const z, + CCTK_REAL const* restrict const r, + CCTK_REAL const& var0, + CCTK_REAL const& v0, + CCTK_REAL const& radpower, + int const width) +{ + int imin[3], imax[3], is_symbnd[6], is_physbnd[6], is_ipbnd[6]; + GenericFD_GetBoundaryInfo + (cctkGH, cctkGH->cctk_lsh, cctkGH->cctk_lssh, cctkGH->cctk_bbox, + cctkGH->cctk_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; + } + } + + 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 + } + } +} + + + +extern "C" +CCTK_INT NewRad_Apply1 (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_REAL const* restrict const var, + CCTK_REAL * restrict const rhs, + CCTK_REAL const var0, + CCTK_REAL const v0, + CCTK_REAL const radpower, + CCTK_INT const width) +{ + cGH const* restrict const cctkGH = static_cast<cGH const*> (cctkGH_); + if (not cctkGH) { + CCTK_WARN (CCTK_WARN_ABORT, + "cctkGH is NULL"); + } + +#if 0 + CCTK_REAL const* 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); + } + CCTK_REAL * restrict const rhs = CCTK_VarDataPtr (cctkGH, 0, rhsname); + if (not rhs) { + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot access RHS variable \"%s\"", rhsname); + } +#endif + + if (not var) { + CCTK_WARN (CCTK_WARN_ABORT, + "Pointer to variable is NULL"); + } + if (not rhs) { + CCTK_WARN (CCTK_WARN_ABORT, + "Pointer to RHS is NULL"); + } + + CCTK_REAL const* restrict const x = + static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::x")); + CCTK_REAL const* restrict const y = + static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::y")); + CCTK_REAL const* restrict const z = + static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::z")); + CCTK_REAL const* restrict const r = + static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::r")); + if (not x or not y or not z or not z) { + CCTK_WARN (CCTK_WARN_ABORT, + "Cannot access coordinate variables x, y, z, and r"); + } + + newrad_loop (cctkGH, var, rhs, x,y,z,r, var0, v0, radpower, width); + + return 0; +} |