aboutsummaryrefslogtreecommitdiff
path: root/src/newrad.cc
diff options
context:
space:
mode:
Diffstat (limited to 'src/newrad.cc')
-rw-r--r--src/newrad.cc311
1 files changed, 311 insertions, 0 deletions
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;
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-11 19:25:57 +0000