aboutsummaryrefslogtreecommitdiff
path: root/src/newrad.cc
diff options
context:
space:
mode:
Diffstat (limited to 'src/newrad.cc')
-rw-r--r--src/newrad.cc165
1 files changed, 102 insertions, 63 deletions
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 <math.h>
+#include <cassert>
+#include <cmath>
#include <cctk.h>
@@ -12,6 +13,8 @@ extern "C" {
#include <GenericFD.h>
}
+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]; k<bmax[2]; ++k) {
- for (int j=bmin[1]; j<bmax[1]; ++j) {
- for (int i=bmin[0]; i<bmax[0]; ++i) {
+ 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);
+
{
// 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 i<ni-1) {
+ if (svi==0) {
derivx = 0.5*(var[ind+di]-var[ind-di])*idx;
} else {
derivx = svi*0.5*(3*var[ind] - 4*var[ind-svi*di]
@@ -89,7 +115,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
// Find y derivative
CCTK_REAL derivy;
- if (j>0 and j<nj-1) {
+ if (svj==0) {
derivy = 0.5*(var[ind+dj]-var[ind-dj])*idy;
} else {
derivy = svj*0.5*(3*var[ind] - 4*var[ind-svj*dj]
@@ -98,7 +124,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
// Find z derivative
CCTK_REAL derivz;
- if (k>0 and k<nk-1) {
+ if (svk==0) {
derivz = 0.5*(var[ind+dk]-var[ind-dk])*idz;
} else {
derivz = svk*0.5*(3*var[ind] - 4*var[ind-svk*dk]
@@ -130,6 +156,9 @@ void newrad_kernel (cGH const* restrict const cctkGH,
int const ip = i-si;
int const jp = j-sj;
int const kp = k-sk;
+ assert (ip>=0 and ip<ni);
+ assert (jp>=0 and jp<nj);
+ assert (kp>=0 and kp<nk);
// Find local wave speeds
int const indp = CCTK_GFINDEX3D(cctkGH, ip,jp,kp);
@@ -141,13 +170,13 @@ void newrad_kernel (cGH const* restrict const cctkGH,
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;
+ 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 ip<ni-1) {
+ if (svi==0) {
derivx = 0.5*(var[indp+di]-var[indp-di])*idx;
} else {
derivx = svi*0.5*(3*var[indp] - 4*var[indp-svi*di]
@@ -156,7 +185,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
// Find y derivative
CCTK_REAL derivy;
- if (jp>0 and jp<nj-1) {
+ if (svj==0) {
derivy = 0.5*(var[indp+dj]-var[indp-dj])*idy;
} else {
derivy = svj*0.5*(3*var[indp] - 4*var[indp-svj*dj]
@@ -165,7 +194,7 @@ void newrad_kernel (cGH const* restrict const cctkGH,
// Find z derivative
CCTK_REAL derivz;
- if (kp>0 and kp<nk-1) {
+ if (svk==0) {
derivz = 0.5*(var[indp+dk]-var[indp-dk])*idz;
} else {
derivz = svk*0.5*(3*var[indp] - 4*var[indp-svk*dk]
@@ -201,8 +230,7 @@ void newrad_loop (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 width)
+ int const radpower)
{
int imin[3], imax[3], is_symbnd[6], is_physbnd[6], is_ipbnd[6];
GenericFD_GetBoundaryInfo
@@ -210,45 +238,57 @@ void newrad_loop (cGH const* restrict const cctkGH,
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;
+ // 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<cGH const*> (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;
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-12 23:54:32 +0000