diff options
Diffstat (limited to 'src/extrap.cc')
| -rw-r--r-- | src/extrap.cc | 26 |
1 files changed, 21 insertions, 5 deletions
diff --git a/src/extrap.cc b/src/extrap.cc index b752aa0..b44d0ae 100644 --- a/src/extrap.cc +++ b/src/extrap.cc @@ -18,6 +18,8 @@ using namespace std; // Adapted from BSSN_MoL's files Init.F +// Erik Schnetter: I assume this code was probably originally written +// by Miguel Alcubierre. static void extrap_kernel (cGH const* restrict const cctkGH, int const* restrict const bmin, @@ -55,6 +57,8 @@ void extrap_kernel (cGH const* restrict const cctkGH, } } + // Note: These loops are not parallel, since each iteration may + // access grid points set by previous iterations. 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]) { @@ -102,9 +106,9 @@ void extrap_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) { @@ -117,10 +121,15 @@ void extrap_loop (cGH const* restrict const cctkGH, } if (nnz == ifec) { - // one of tahe faces is a boundary + // one of the faces is a boundary 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) { @@ -130,6 +139,9 @@ void extrap_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 is_ipbnd[2*d+0]); break; case 0: bmin[d] = imin[d]; @@ -140,11 +152,15 @@ void extrap_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) { extrap_kernel (cctkGH, bmin, bmax, dir, var); } |