diff options
Diffstat (limited to 'src/newrad.cc')
-rw-r--r-- | src/newrad.cc | 70 |
1 files changed, 66 insertions, 4 deletions
diff --git a/src/newrad.cc b/src/newrad.cc index f87e78b..c39c360 100644 --- a/src/newrad.cc +++ b/src/newrad.cc @@ -18,6 +18,8 @@ using namespace std; // Adapted from BSSN_MoL's files NewRad.F and newrad.h +// Erik Schnetter: I assume this code was probably originally written +// by Miguel Alcubierre. static void newrad_kernel (cGH const* restrict const cctkGH, int const* restrict const bmin, @@ -93,6 +95,22 @@ void newrad_kernel (cGH const* restrict const cctkGH, if (j==nj-1) assert (idir[1]>0); if (k==nk-1) assert (idir[2]>0); + if (si==0) { + assert (i-1>=0 and i+1<ni); + } else { + assert (i-2*si>=0 and i-2*si<ni); + } + if (sj==0) { + assert (j-1>=0 and j+1<nj); + } else { + assert (j-2*sj>=0 and j-2*sj<nj); + } + if (sk==0) { + assert (k-1>=0 and k+1<nk); + } else { + assert (k-2*sk>=0 and k-2*sk<nk); + } + { // The main part of the boundary condition assumes that we // have an outgoing radial wave with some speed v0: @@ -105,6 +123,22 @@ void newrad_kernel (cGH const* restrict const cctkGH, // // where vi = v0 xi/r + if (si==0) { + assert (i-1>=0 and i+1<ni); + } else { + assert (i-2*si>=0 and i-2*si<ni); + } + if (sj==0) { + assert (j-1>=0 and j+1<nj); + } else { + assert (j-2*sj>=0 and j-2*sj<nj); + } + if (sk==0) { + assert (k-1>=0 and k+1<nk); + } else { + assert (k-2*sk>=0 and k-2*sk<nk); + } + // Find local wave speeds CCTK_REAL const rp = r[ind]; CCTK_REAL const rpi = 1.0/rp; @@ -169,6 +203,22 @@ void newrad_kernel (cGH const* restrict const cctkGH, assert (jp>=0 and jp<nj); assert (kp>=0 and kp<nk); + if (si==0) { + assert (ip-1>=0 and ip+1<ni); + } else { + assert (ip-2*si>=0 and ip-2*si<ni); + } + if (sj==0) { + assert (jp-1>=0 and jp+1<nj); + } else { + assert (jp-2*sj>=0 and jp-2*sj<nj); + } + if (sk==0) { + assert (kp-1>=0 and kp+1<nk); + } else { + assert (kp-2*sk>=0 and kp-2*sk<nk); + } + // Find local wave speeds int const indp = CCTK_GFINDEX3D(cctkGH, ip,jp,kp); @@ -244,9 +294,9 @@ void newrad_loop (cGH const* restrict const cctkGH, 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 + // Loop over faces first, then corners, and then 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) { @@ -263,6 +313,11 @@ void newrad_loop (cGH const* restrict const cctkGH, bool have_bnd = false; // all boundary faces are physical boundaries bool all_physbnd = true; + // at least one boundary face is a physical boundary + bool any_physbnd = false; + // all boundary faces are either physical boundaries or + // ghost zones + bool all_physbnd_or_ghostbnd = true; int bmin[3], bmax[3]; for (int d=0; d<3; ++d) { @@ -272,6 +327,9 @@ void newrad_loop (cGH const* restrict const cctkGH, bmax[d] = imin[d]; have_bnd = true; all_physbnd = all_physbnd and is_physbnd[2*d+0]; + any_physbnd = any_physbnd or is_physbnd[2*d+0]; + all_physbnd_or_ghostbnd = all_physbnd_or_ghostbnd and + (is_physbnd[2*d+0] or not is_ipbnd[2*d+0]); break; case 0: bmin[d] = imin[d]; @@ -282,11 +340,15 @@ void newrad_loop (cGH const* restrict const cctkGH, bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)]; have_bnd = true; all_physbnd = all_physbnd and is_physbnd[2*d+1]; + any_physbnd = any_physbnd or is_physbnd[2*d+1]; + all_physbnd_or_ghostbnd = all_physbnd_or_ghostbnd and + (is_physbnd[2*d+1] or not is_ipbnd[2*d+1]); break; } } + assert (have_bnd); // must be true since nnz>0 - if (have_bnd and all_physbnd) { + if (have_bnd and any_physbnd and all_physbnd_or_ghostbnd) { newrad_kernel (cctkGH, bmin, bmax, dir, var, rhs, x,y,z,r, var0, v0, radpower); } |