Erik's changes to make things timelevel aware.

Tom


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Boundary/trunk@85 6a38eb6e-646e-4a02-a296-d141613ad6c4

1 files changed, 81 insertions, 89 deletions

diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c
index 866954c..e48df15 100644
--- a/src/RadiationBoundary.c
+++ b/src/RadiationBoundary.c
@@ -46,18 +46,15 @@
             a wave speed of v0, so the boundaries should be out in
             the region where the characteristic speed is constant.
             Notice that this speed does not have to be 1.  For gauge
-            quantities {alpha,phi,trK} we can have a different
-asymptotic
-            speed, which is why the value of v0 is passed as a
-parameter.
+            quantities {alpha,phi,trK} we can have a different asymptotic
+            speed, which is why the value of v0 is passed as a parameter.
 
 	    Some remarks on the C version for Fortran programmers:
 
 	    This is 1:1 translation of the F code, which used the very
             nice looking F array assignments x(2,:,:), etc. In C this
             needs to be written as loops, but to my surprise it doesn't
-            look to bad actually. F statments like x(1,:,:) become
-x(xgp0)
+            look too bad actually. F statments like x(1,:,:) become x(xgp0)
             (the integer index "xgp0" meaning the zeroth gridpoint from
             the boundary), x(2,:,:) --> x(xgp1), where:
 
@@ -129,29 +126,24 @@ void RadiativeBCi(cGH *GH,
 
   dth = half*dt;
 
-  /* Print warning. */
-
-  if ((sw[0]!=1) || (sw[1]!=1) || (sw[2]!=1))
-     CCTK_WARN(0,"Radiative BC only with stencil_size=1");
-
   /* Lower x-bound: x(2,:,:) --> x[xgp2] */
 
   if (doBC[0]==1) {
     for (k=0;k<lsh[2];k++) {
       for (j=0;j<lsh[1];j++) {
-
-        xgp0 = CCTK_GFINDEX3D(GH,0,j,k);
-	xgp1 = CCTK_GFINDEX3D(GH,1,j,k);
-
-        var_n[xgp0] = ( dtv*var0* ( x[xgp0]/(SQR(r[xgp0])) 
-              + x[xgp1]/(SQR(r[xgp1])))
-	      - var_n[xgp1]*(rhox + x[xgp1]/r[xgp1]
-              * (one + v0 * dth/r[xgp1]))
-              + var_p[xgp0]*(rhox + x[xgp0]/r[xgp0]
-              * (one - v0*dth/r[xgp0]))
-              - var_p[xgp1]*(rhox - x[xgp1]/r[xgp1]
-              * (one -v0*dth/r[xgp1])))
-              / (-rhox +x[xgp0]/r[xgp0] * (one +v0*dth/r[xgp0]));
+	for (i=sw[0]-1;i>=0;i--) {
+	  
+	  xgp0 = CCTK_GFINDEX3D(GH,i  ,j,k);
+	  xgp1 = CCTK_GFINDEX3D(GH,i+1,j,k);
+	  
+	  var_n[xgp0] =
+	    (dtv*var0 * (x[xgp0]/(SQR(r[xgp0])) + x[xgp1]/(SQR(r[xgp1])))
+	     - var_n[xgp1] * ( rhox + x[xgp1]/r[xgp1] * (one+v0*dth/r[xgp1]))
+	     + var_p[xgp0] * ( rhox + x[xgp0]/r[xgp0] * (one-v0*dth/r[xgp0]))
+	     - var_p[xgp1] * ( rhox - x[xgp1]/r[xgp1] * (one-v0*dth/r[xgp1])))
+	    /                (-rhox + x[xgp0]/r[xgp0] * (one+v0*dth/r[xgp0]));
+	  
+	}
       }
     }
   }
@@ -161,19 +153,19 @@ void RadiativeBCi(cGH *GH,
   if (doBC[1]==1){
     for (k=0;k<lsh[2];k++) {
       for (j=0;j<lsh[1];j++) {
-
-	xgp0 = CCTK_GFINDEX3D(GH,lsh[0]-1,j,k);
-        xgp1 = CCTK_GFINDEX3D(GH,lsh[0]-2,j,k);
-
-        var_n[xgp0] = (dtv*var0*(x[xgp0]/(SQR(r[xgp0]))
-              + x[xgp1]/(SQR(r[xgp1])))
-              + var_n[xgp1]*(rhox - x[xgp1]/r[xgp1]
-              *(one + v0*dth/r[xgp1]))
-              + var_p[xgp0]*(-rhox + x[xgp0]/r[xgp0]
-              *(one - v0*dth/r[xgp0]))
-              + var_p[xgp1]*(rhox + x[xgp1]/r[xgp1]
-              *(one - v0*dth/r[xgp1])))
-              /(rhox + x[xgp0]/r[xgp0]*(one + v0*dth/r[xgp0]));
+	for (i=lsh[0]-sw[0];i<lsh[0];i++) {
+	  
+	  xgp0 = CCTK_GFINDEX3D(GH,i  ,j,k);
+	  xgp1 = CCTK_GFINDEX3D(GH,i-1,j,k);
+	  
+	  var_n[xgp0] =
+	    (dtv*var0 * (x[xgp0]/(SQR(r[xgp0])) + x[xgp1]/(SQR(r[xgp1])))
+	     + var_n[xgp1] * ( rhox - x[xgp1]/r[xgp1] * (one+v0*dth/r[xgp1]))
+	     + var_p[xgp0] * (-rhox + x[xgp0]/r[xgp0] * (one-v0*dth/r[xgp0]))
+	     + var_p[xgp1] * ( rhox + x[xgp1]/r[xgp1] * (one-v0*dth/r[xgp1])))
+	    /                ( rhox + x[xgp0]/r[xgp0] * (one+v0*dth/r[xgp0]));
+	  
+	}
       }
     }
   }
@@ -183,19 +175,19 @@ void RadiativeBCi(cGH *GH,
   if (doBC[2] == 1) {
     for (k=0;k<lsh[2];k++) {
       for (i=0;i<lsh[0];i++) {
-
-	ygp0 = CCTK_GFINDEX3D(GH,i,0,k);
-	ygp1 = CCTK_GFINDEX3D(GH,i,1,k);
-
-	var_n[ygp0] = (dtv*var0*(y[ygp0]/(SQR(r[ygp0]))
-              + y[ygp1]/(SQR(r[ygp1])))
-              - var_n[ygp1]*(rhoy + y[ygp1]/r[ygp1]
-              *(one + v0*dth/r[ygp1]))
-              + var_p[ygp0]*(rhoy + y[ygp0]/r[ygp0]
-              *(one - v0*dth/r[ygp0]))
-              - var_p[ygp1]*(rhoy - y[ygp1]/r[ygp1]
-              *(one - v0*dth/r[ygp1])))
-              /(-rhoy + y[ygp0]/r[ygp0]*(one + v0*dth/r[ygp0]));
+	for (j=sw[1]-1;j>=0;j--) {
+	  
+	  ygp0 = CCTK_GFINDEX3D(GH,i,j  ,k);
+	  ygp1 = CCTK_GFINDEX3D(GH,i,j+1,k);
+	  
+	  var_n[ygp0] =
+	    (dtv*var0 * (y[ygp0]/(SQR(r[ygp0])) + y[ygp1]/(SQR(r[ygp1])))
+	     - var_n[ygp1] * ( rhoy + y[ygp1]/r[ygp1] * (one+v0*dth/r[ygp1]))
+	     + var_p[ygp0] * ( rhoy + y[ygp0]/r[ygp0] * (one-v0*dth/r[ygp0]))
+	     - var_p[ygp1] * ( rhoy - y[ygp1]/r[ygp1] * (one-v0*dth/r[ygp1])))
+	    /                (-rhoy + y[ygp0]/r[ygp0] * (one+v0*dth/r[ygp0]));
+	  
+	}
       }
     }
   }
@@ -205,19 +197,19 @@ void RadiativeBCi(cGH *GH,
   if (doBC[3] == 1) {
     for (k=0;k<lsh[2];k++) {
       for (i=0;i<lsh[0];i++) {
-
-        ygp0 = CCTK_GFINDEX3D(GH,i,lsh[1]-1,k);
-	ygp1 = CCTK_GFINDEX3D(GH,i,lsh[1]-2,k);
-
-	var_n[ygp0] = (dtv*var0*(y[ygp0]/(SQR(r[ygp0]))
-              + y[ygp1]/(SQR(r[ygp1])))
-              + var_n[ygp1]*(rhoy - y[ygp1]/r[ygp1]
-              *(one + v0*dth/r[ygp1]))
-              + var_p[ygp0]*(-rhoy + y[ygp0]/r[ygp0]
-              *(one - v0*dth/r[ygp0]))
-              + var_p[ygp1]*(rhoy + y[ygp1]/r[ygp1]
-              *(one - v0*dth/r[ygp1])))
-              /(rhoy + y[ygp0]/r[ygp0]*(one + v0*dth/r[ygp0]));
+	for (j=lsh[1]-sw[1];j<lsh[1];j++) {
+	  
+	  ygp0 = CCTK_GFINDEX3D(GH,i,j  ,k);
+	  ygp1 = CCTK_GFINDEX3D(GH,i,j-1,k);
+	  
+	  var_n[ygp0] =
+	    (dtv*var0 * (y[ygp0]/(SQR(r[ygp0])) + y[ygp1]/(SQR(r[ygp1])))
+	     + var_n[ygp1] * ( rhoy - y[ygp1]/r[ygp1] * (one+v0*dth/r[ygp1]))
+	     + var_p[ygp0] * (-rhoy + y[ygp0]/r[ygp0] * (one-v0*dth/r[ygp0]))
+	     + var_p[ygp1] * ( rhoy + y[ygp1]/r[ygp1] * (one-v0*dth/r[ygp1])))
+	    /                ( rhoy + y[ygp0]/r[ygp0] * (one+v0*dth/r[ygp0]));
+	  
+	}
       }
     }
   }
@@ -227,19 +219,19 @@ void RadiativeBCi(cGH *GH,
   if (doBC[4]==1) {
     for (j=0;j<lsh[1];j++) {
       for (i=0;i<lsh[0];i++) {
-
-	zgp0 = CCTK_GFINDEX3D(GH,i,j,0);
-        zgp1 = CCTK_GFINDEX3D(GH,i,j,1);
-
-	var_n[zgp0] = (dtv*var0*(z[zgp0]/(SQR(r[zgp0]))
-              + z[zgp1]/(SQR(r[zgp1])))
-              - var_n[zgp1]*(rhoz + z[zgp1]/r[zgp1]
-              *(one + v0*dth/r[zgp1]))
-              + var_p[zgp0]*(rhoz + z[zgp0]/r[zgp0]
-              *(one - v0*dth/r[zgp0]))
-              - var_p[zgp1]*(rhoz - z[zgp1]/r[zgp1]
-              *(one - v0*dth/r[zgp1])))
-              /(-rhoz + z[zgp0]/r[zgp0]*(one + v0*dth/r[zgp0]));	
+	for (k=sw[2]-1;k>=0;k--) {
+	  
+	  zgp0 = CCTK_GFINDEX3D(GH,i,j,k  );
+	  zgp1 = CCTK_GFINDEX3D(GH,i,j,k+1);
+	  
+	  var_n[zgp0] =
+	    (dtv*var0 * (z[zgp0]/(SQR(r[zgp0])) + z[zgp1]/(SQR(r[zgp1])))
+	     - var_n[zgp1] * ( rhoz + z[zgp1]/r[zgp1] * (one+v0*dth/r[zgp1]))
+	     + var_p[zgp0] * ( rhoz + z[zgp0]/r[zgp0] * (one-v0*dth/r[zgp0]))
+	     - var_p[zgp1] * ( rhoz - z[zgp1]/r[zgp1] * (one-v0*dth/r[zgp1])))
+	    /                (-rhoz + z[zgp0]/r[zgp0] * (one+v0*dth/r[zgp0]));
+	  
+	}
       }
     }
   }
@@ -249,19 +241,19 @@ void RadiativeBCi(cGH *GH,
   if (doBC[5] == 1) {
     for (j=0;j<lsh[1];j++) {
       for (i=0;i<lsh[0];i++) {
-
-	zgp0 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-1);
-	zgp1 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-2);
-
-	var_n[zgp0] = (dtv*var0*(z[zgp0]/(SQR(r[zgp0]))
-              + z[zgp1]/(SQR(r[zgp1])))
-              + var_n[zgp1]*(rhoz - z[zgp1]/r[zgp1]
-              *(one + v0*dth/r[zgp1]))
-              + var_p[zgp0]*(-rhoz + z[zgp0]/r[zgp0]
-              *(one - v0*dth/r[zgp0]))
-              + var_p[zgp1]*(rhoz + z[zgp1]/r[zgp1]
-              *(one - v0*dth/r[zgp1])))
-              /(rhoz + z[zgp0]/r[zgp0]*(one + v0*dth/r[zgp0]));
+	for (k=lsh[2]-sw[2];k<lsh[2];k++) {
+	  
+	  zgp0 = CCTK_GFINDEX3D(GH,i,j,k  );
+	  zgp1 = CCTK_GFINDEX3D(GH,i,j,k-1);
+	  
+	  var_n[zgp0] =
+	    (dtv*var0 * (z[zgp0]/(SQR(r[zgp0])) + z[zgp1]/(SQR(r[zgp1])))
+	     + var_n[zgp1] * ( rhoz - z[zgp1]/r[zgp1] * (one+v0*dth/r[zgp1]))
+	     + var_p[zgp0] * (-rhoz + z[zgp0]/r[zgp0] * (one-v0*dth/r[zgp0]))
+	     + var_p[zgp1] * ( rhoz + z[zgp1]/r[zgp1] * (one-v0*dth/r[zgp1])))
+	    /                ( rhoz + z[zgp0]/r[zgp0] * (one+v0*dth/r[zgp0]));
+	  
+	}
       }
     }
   }