diff options
| -rw-r--r-- | param.ccl | 8 | ||||
| -rw-r--r-- | src/RadiationBoundary.c | 414 |
2 files changed, 341 insertions, 81 deletions
@@ -1,2 +1,10 @@ # Parameter definitions for thorn Boundary # $Header$ + + +private: + +INT radpower "Power of decay rate in extrapolation used in radiative boundaries" +{ + : :: "A negative value switches off this feature" +} -1
\ No newline at end of file diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c index f4610d6..ef760ad 100644 --- a/src/RadiationBoundary.c +++ b/src/RadiationBoundary.c @@ -17,37 +17,64 @@ /*@@ @routine RadiativeBC - @date Mon Mar 15 15:51:57 1999 + @date May 2000 @author Miguel Alcubierre @desc Implements radiative boundary condition. - The boundary condition that is implemented is: + The radiative boundary condition that is implemented is: - f = f0 + u(r - v0*t) / r + f = f0 + u(r - v*t) / r + h(r + v*t) / r - which leads to the differential equation: + That is, I assume outgoing radial waves with a 1/r + fall off, and the correct asymptotic value f0, plus + I include the possibility of incoming waves as well + (these incoming waves should be modeled somehow). - (x / r) d f + v0 d f + v0 x (f - f0) / r^2 = 0 - i t i i + The condition above leads to the differential equation: - where xi is the normal direction to the given boundaries. + + (x / r) d f + v d f + v x (f - f0) / r^2 = v x H / r^2 + i t i i i - That is, I assume outgoing radial waves with the correct 1/r - fall off, and the correct asymptotic value f0 (f0 = 1 for - diagonal metric components, and f0 = 0 for everything else), - but at a given boundary I only worry about derivatives in - the normal direction. This works very well. + where x_i is the normal direction to the given boundaries, + and H = 2 dh(s)/ds. - The only problem with this condition is that it does not + So at a given boundary I only worry about derivatives in + the normal direction. Notice that u(r-v*t) has dissapeared, + but we still do not know the value of H. + + To get H what I do is the following: I evaluate the + expression one point in from the boundary and solve for H + there. We now need a way of extrapolation H to the boundary. + For this I assume that H falls off as a power law: + + + H = k/r**n => d H = - n H/r + i + + The value of n is is defined by the parameter "radpower". + If this parameter is negative, H is forced to be zero (this + corresponds to pure outgoing waves and is the default). + + The behaviour I have observed is the following: Using H=0 + is very stable, but has a very bad initial transient. Taking + n to be 0 or positive improves the initial behaviour considerably, + but introduces a drift that can kill the evolution at very late + times. Empirically, the best value I have found is n=2, for + which the initial behaviour is very nice, and the late time drift + is quite small. + + Another problem with this condition is that it does not use the physical characteristic speed, but rather it assumes - a wave speed of v0, so the boundaries should be out in + a wave speed of v, 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. + speed, which is why the value of v is passed as a parameter. + Some remarks on the C version for Fortran programmers: @@ -90,11 +117,15 @@ int BndApplyRadiative3Di(cGH *GH, CCTK_REAL var0, CCTK_REAL v0) { + + DECLARE_CCTK_PARAMETERS + int i,j,k; - CCTK_REAL dtv,dth; + CCTK_REAL dtv,dtvh,dtvvar0; CCTK_REAL dx,dy,dz; CCTK_REAL rhox,rhoy,rhoz; + CCTK_REAL H,fac,aux; CCTK_REAL half,one; /* Linear grid point index, used to access the 0th/1st gridpoint away @@ -104,11 +135,12 @@ int BndApplyRadiative3Di(cGH *GH, x(2,:,:) --> x[xgp1], where xgp0/1 are calculated with CCTK_GFINDEX3D */ - int xgp0,xgp1,ygp0,ygp1,zgp0,zgp1; - - /* Inverses of the radius */ + int xgp0,xgp1,xgp2,ygp0,ygp1,ygp2,zgp0,zgp1,zgp2; - CCTK_REAL ri0, ri1; + CCTK_REAL u0,u1,u2; + CCTK_REAL ui1,ui2; + CCTK_REAL r0,r1,r2; + CCTK_REAL ri0,ri1; /* Grid sizes */ @@ -143,13 +175,22 @@ int BndApplyRadiative3Di(cGH *GH, /* Find Courant parameters. */ - dtv = v0*dt; + dtv = v0*dt; + dtvh = half*dtv; + + dtvvar0 = dtv*var0; rhox = dtv/dx; rhoy = dtv/dy; rhoz = dtv/dz; - dth = half*dt; + /* Extrapolation power */ + + if (radpower<0) { + fac = 0; + } else { + fac = radpower; + } /* Lower x-bound: x(2,:,:) --> x[xgp2] */ @@ -157,19 +198,55 @@ int BndApplyRadiative3Di(cGH *GH, for (k=0;k<lssh2;k++) { for (j=0;j<lssh1;j++) { for (i=sw0-1;i>=0;i--) { - + xgp0 = CCTK_GFINDEX3D(GH,i ,j,k); xgp1 = CCTK_GFINDEX3D(GH,i+1,j,k); - ri0 = 1.0/r[xgp0]; - ri1 = 1.0/r[xgp1]; + xgp2 = CCTK_GFINDEX3D(GH,i+2,j,k); + + u0 = x[xgp0]; + u1 = x[xgp1]; + u2 = x[xgp2]; + + r0 = r[xgp0]; + r1 = r[xgp1]; + r2 = r[xgp2]; + + ui1 = 1.0/u1; + ui2 = 1.0/u2; + + ri0 = 1.0/r0; + ri1 = 1.0/r1; + + if (fac==0) { + + H = 0; + + } else { + + H = 0; + + H = dtv*(-var0 + 0.25 + *(var_n[xgp1] + var_n[xgp2] + + var_p[xgp1] + var_p[xgp2])) + +(r1*(var_n[xgp1] - var_p[xgp1]) + + r2*(var_n[xgp2] - var_p[xgp2]))*half + + 0.25*rhox*(SQR(r1)*ui1 + SQR(r2)*ui2) + *(var_n[xgp2] - var_n[xgp1] + + var_p[xgp2] - var_p[xgp1]); + + aux = 0.25*radpower*dx*(u0*SQR(ri0) + u1*SQR(ri1)); + + H = H*(one + aux)/(one - aux); + + } var_n[xgp0] = - (dtv*var0 * (x[xgp0]*(SQR(ri0)) + x[xgp1]*(SQR(ri1))) - - var_n[xgp1] * ( rhox + x[xgp1]*ri1 * (one+v0*dth*ri1)) - + var_p[xgp0] * ( rhox + x[xgp0]*ri0 * (one-v0*dth*ri0)) - - var_p[xgp1] * ( rhox - x[xgp1]*ri1 * (one-v0*dth*ri1))) - / (-rhox + x[xgp0]*ri0 * (one+v0*dth*ri0)); - + ((dtvvar0 + H)*(u0*SQR(ri0) + u1*SQR(ri1)) + - var_n[xgp1]*( rhox + u1*ri1*(one + dtvh*ri1)) + + var_p[xgp0]*( rhox + u0*ri0*(one - dtvh*ri0)) + - var_p[xgp1]*( rhox - u1*ri1*(one - dtvh*ri1))) + / (-rhox + u0*ri0*(one + dtvh*ri0)); + } } } @@ -181,20 +258,55 @@ int BndApplyRadiative3Di(cGH *GH, for (k=0;k<lssh2;k++) { for (j=0;j<lssh1;j++) { for (i=lssh0-sw0;i<lssh0;i++) { - + xgp0 = CCTK_GFINDEX3D(GH,i ,j,k); xgp1 = CCTK_GFINDEX3D(GH,i-1,j,k); + xgp2 = CCTK_GFINDEX3D(GH,i-2,j,k); + + u0 = x[xgp0]; + u1 = x[xgp1]; + u2 = x[xgp2]; + + r0 = r[xgp0]; + r1 = r[xgp1]; + r2 = r[xgp2]; + + ui1 = 1.0/u1; + ui2 = 1.0/u2; + + ri0 = 1.0/r0; + ri1 = 1.0/r1; + + if (fac==0) { + + H = 0; + + } else { + + H = 0; + + H = dtv*(-var0 + 0.25 + *(var_n[xgp1] + var_n[xgp2] + + var_p[xgp1] + var_p[xgp2])) + +(r1*(var_n[xgp1] - var_p[xgp1]) + + r2*(var_n[xgp2] - var_p[xgp2]))*half + + 0.25*rhox*(SQR(r1)*ui1 + SQR(r2)*ui2) + *(var_n[xgp1] - var_n[xgp2] + + var_p[xgp1] - var_p[xgp2]); + + aux = 0.25*radpower*dx*(u0*SQR(ri0) + u1*SQR(ri1)); + + H = H*(one - aux)/(one + aux); + + } - ri0 = 1.0/r[xgp0]; - ri1 = 1.0/r[xgp1]; - var_n[xgp0] = - (dtv*var0 * (x[xgp0]*(SQR(ri0)) + x[xgp1]*(SQR(ri1))) - + var_n[xgp1] * ( rhox - x[xgp1]*ri1 * (one+v0*dth*ri1)) - + var_p[xgp0] * (-rhox + x[xgp0]*ri0 * (one-v0*dth*ri0)) - + var_p[xgp1] * ( rhox + x[xgp1]*ri1 * (one-v0*dth*ri1))) - / ( rhox + x[xgp0]*ri0 * (one+v0*dth*ri0)); - + ((dtvvar0 + H)*(u0*(SQR(ri0)) + u1*(SQR(ri1))) + + var_n[xgp1]*( rhox - u1*ri1*(one + dtvh*ri1)) + + var_p[xgp0]*(-rhox + u0*ri0*(one - dtvh*ri0)) + + var_p[xgp1]*( rhox + u1*ri1*(one - dtvh*ri1))) + / ( rhox + u0*ri0*(one + dtvh*ri0)); + } } } @@ -206,20 +318,55 @@ int BndApplyRadiative3Di(cGH *GH, for (k=0;k<lssh2;k++) { for (i=0;i<lssh0;i++) { for (j=sw1-1;j>=0;j--) { - + ygp0 = CCTK_GFINDEX3D(GH,i,j ,k); ygp1 = CCTK_GFINDEX3D(GH,i,j+1,k); + ygp2 = CCTK_GFINDEX3D(GH,i,j+2,k); + + u0 = y[ygp0]; + u1 = y[ygp1]; + u2 = y[ygp2]; + + r0 = r[ygp0]; + r1 = r[ygp1]; + r2 = r[ygp2]; + + ui1 = 1.0/u1; + ui2 = 1.0/u2; + + ri0 = 1.0/r0; + ri1 = 1.0/r1; + + if (fac==0) { + + H = 0; + + } else { + + H = 0; + + H = dtv*(-var0 + 0.25 + *(var_n[ygp1] + var_n[ygp2] + + var_p[ygp1] + var_p[ygp2])) + +(r1*(var_n[ygp1] - var_p[ygp1]) + + r2*(var_n[ygp2] - var_p[ygp2]))*half + + 0.25*rhoy*(SQR(r1)*ui1 + SQR(r2)*ui2) + *(var_n[ygp2] - var_n[ygp1] + + var_p[ygp2] - var_p[ygp1]); + + aux = 0.25*radpower*dy*(u0*SQR(ri0) + u1*SQR(ri1)); + + H = H*(one + aux)/(one - aux); + + } - ri0 = 1.0/r[ygp0]; - ri1 = 1.0/r[ygp1]; - var_n[ygp0] = - (dtv*var0 * (y[ygp0]*(SQR(ri0)) + y[ygp1]*(SQR(ri1))) - - var_n[ygp1] * ( rhoy + y[ygp1]*ri1 * (one+v0*dth*ri1)) - + var_p[ygp0] * ( rhoy + y[ygp0]*ri0 * (one-v0*dth*ri0)) - - var_p[ygp1] * ( rhoy - y[ygp1]*ri1 * (one-v0*dth*ri1))) - / (-rhoy + y[ygp0]*ri0 * (one+v0*dth*ri0)); - + ((dtvvar0 + H)*(u0*SQR(ri0) + u1*SQR(ri1)) + - var_n[ygp1]*( rhoy + u1*ri1*(one + dtvh*ri1)) + + var_p[ygp0]*( rhoy + u0*ri0*(one - dtvh*ri0)) + - var_p[ygp1]*( rhoy - u1*ri1*(one - dtvh*ri1))) + / (-rhoy + u0*ri0*(one + dtvh*ri0)); + } } } @@ -231,20 +378,55 @@ int BndApplyRadiative3Di(cGH *GH, for (k=0;k<lssh2;k++) { for (i=0;i<lssh0;i++) { for (j=lssh1-sw1;j<lssh1;j++) { - + ygp0 = CCTK_GFINDEX3D(GH,i,j ,k); ygp1 = CCTK_GFINDEX3D(GH,i,j-1,k); - - ri0 = 1.0/r[ygp0]; - ri1 = 1.0/r[ygp1]; + ygp2 = CCTK_GFINDEX3D(GH,i,j-2,k); + + u0 = y[ygp0]; + u1 = y[ygp1]; + u2 = y[ygp2]; + + r0 = r[ygp0]; + r1 = r[ygp1]; + r2 = r[ygp2]; + + ui1 = 1.0/u1; + ui2 = 1.0/u2; + + ri0 = 1.0/r0; + ri1 = 1.0/r1; + + if (fac==0) { + + H = 0; + + } else { + + H = 0; + + H = dtv*(-var0 + 0.25 + *(var_n[ygp1] + var_n[ygp2] + + var_p[ygp1] + var_p[ygp2])) + +(r1*(var_n[ygp1] - var_p[ygp1]) + + r2*(var_n[ygp2] - var_p[ygp2]))*half + + 0.25*rhoy*(SQR(r1)*ui1 + SQR(r2)*ui2) + *(var_n[ygp1] - var_n[ygp2] + + var_p[ygp1] - var_p[ygp2]); + + aux = 0.25*radpower*dy*(u0*SQR(ri0) + u1*SQR(ri1)); + + H = H*(one - aux)/(one + aux); + + } var_n[ygp0] = - (dtv*var0 * (y[ygp0]*(SQR(ri0)) + y[ygp1]*(SQR(ri1))) - + var_n[ygp1] * ( rhoy - y[ygp1]*ri1 * (one+v0*dth*ri1)) - + var_p[ygp0] * (-rhoy + y[ygp0]*ri0 * (one-v0*dth*ri0)) - + var_p[ygp1] * ( rhoy + y[ygp1]*ri1 * (one-v0*dth*ri1))) - / ( rhoy + y[ygp0]*ri0 * (one+v0*dth*ri0)); - + ((dtvvar0 + H)*(u0*(SQR(ri0)) + u1*(SQR(ri1))) + + var_n[ygp1]*( rhoy - u1*ri1*(one + dtvh*ri1)) + + var_p[ygp0]*(-rhoy + u0*ri0*(one - dtvh*ri0)) + + var_p[ygp1]*( rhoy + u1*ri1*(one - dtvh*ri1))) + / ( rhoy + u0*ri0*(one + dtvh*ri0)); + } } } @@ -256,20 +438,55 @@ int BndApplyRadiative3Di(cGH *GH, for (j=0;j<lssh1;j++) { for (i=0;i<lssh0;i++) { for (k=sw2-1;k>=0;k--) { - + zgp0 = CCTK_GFINDEX3D(GH,i,j,k ); zgp1 = CCTK_GFINDEX3D(GH,i,j,k+1); + zgp2 = CCTK_GFINDEX3D(GH,i,j,k+2); + + u0 = z[zgp0]; + u1 = z[zgp1]; + u2 = z[zgp2]; + + r0 = r[zgp0]; + r1 = r[zgp1]; + r2 = r[zgp2]; + + ui1 = 1.0/u1; + ui2 = 1.0/u2; + + ri0 = 1.0/r0; + ri1 = 1.0/r1; + + if (fac==0) { + + H = 0; + + } else { + + H = 0; + + H = dtv*(-var0 + 0.25 + *(var_n[zgp1] + var_n[zgp2] + + var_p[zgp1] + var_p[zgp2])) + +(r1*(var_n[zgp1] - var_p[zgp1]) + + r2*(var_n[zgp2] - var_p[zgp2]))*half + + 0.25*rhoz*(SQR(r1)*ui1 + SQR(r2)*ui2) + *(var_n[zgp2] - var_n[zgp1] + + var_p[zgp2] - var_p[zgp1]); + + aux = 0.25*radpower*dz*(u0*SQR(ri0) + u1*SQR(ri1)); + + H = H*(one + aux)/(one - aux); + + } - ri0 = 1.0/r[zgp0]; - ri1 = 1.0/r[zgp1]; - var_n[zgp0] = - (dtv*var0 * (z[zgp0]*(SQR(ri0)) + z[zgp1]*(SQR(ri1))) - - var_n[zgp1] * ( rhoz + z[zgp1]*ri1 * (one+v0*dth*ri1)) - + var_p[zgp0] * ( rhoz + z[zgp0]*ri0 * (one-v0*dth*ri0)) - - var_p[zgp1] * ( rhoz - z[zgp1]*ri1 * (one-v0*dth*ri1))) - / (-rhoz + z[zgp0]*ri0 * (one+v0*dth*ri0)); - + ((dtvvar0 + H)*(u0*SQR(ri0) + u1*SQR(ri1)) + - var_n[zgp1]*( rhoz + u1*ri1*(one + dtvh*ri1)) + + var_p[zgp0]*( rhoz + u0*ri0*(one - dtvh*ri0)) + - var_p[zgp1]*( rhoz - u1*ri1*(one - dtvh*ri1))) + / (-rhoz + u0*ri0*(one + dtvh*ri0)); + } } } @@ -281,20 +498,55 @@ int BndApplyRadiative3Di(cGH *GH, for (j=0;j<lssh1;j++) { for (i=0;i<lssh0;i++) { for (k=lssh2-sw2;k<lssh2;k++) { - + zgp0 = CCTK_GFINDEX3D(GH,i,j,k ); zgp1 = CCTK_GFINDEX3D(GH,i,j,k-1); + zgp2 = CCTK_GFINDEX3D(GH,i,j,k-2); + + u0 = z[zgp0]; + u1 = z[zgp1]; + u2 = z[zgp2]; + + r0 = r[zgp0]; + r1 = r[zgp1]; + r2 = r[zgp2]; + + ui1 = 1.0/u1; + ui2 = 1.0/u2; + + ri0 = 1.0/r0; + ri1 = 1.0/r1; + + if (fac==0) { + + H = 0; + + } else { + + H = 0; + + H = dtv*(-var0 + 0.25 + *(var_n[zgp1] + var_n[zgp2] + + var_p[zgp1] + var_p[zgp2])) + +(r1*(var_n[zgp1] - var_p[zgp1]) + + r2*(var_n[zgp2] - var_p[zgp2]))*half + + 0.25*rhoz*(SQR(r1)*ui1 + SQR(r2)*ui2) + *(var_n[zgp1] - var_n[zgp2] + + var_p[zgp1] - var_p[zgp2]); + + aux = 0.25*radpower*dz*(u0*SQR(ri0) + u1*SQR(ri1)); + + H = H*(one - aux)/(one + aux); + + } - ri0 = 1.0/r[zgp0]; - ri1 = 1.0/r[zgp1]; - var_n[zgp0] = - (dtv*var0 * (z[zgp0]*(SQR(ri0)) + z[zgp1]*(SQR(ri1))) - + var_n[zgp1] * ( rhoz - z[zgp1]*ri1 * (one+v0*dth*ri1)) - + var_p[zgp0] * (-rhoz + z[zgp0]*ri0 * (one-v0*dth*ri0)) - + var_p[zgp1] * ( rhoz + z[zgp1]*ri1 * (one-v0*dth*ri1))) - / ( rhoz + z[zgp0]*ri0 * (one+v0*dth*ri0)); - + ((dtvvar0 + H)*(u0*(SQR(ri0)) + u1*(SQR(ri1))) + + var_n[zgp1]*( rhoz - u1*ri1*(one + dtvh*ri1)) + + var_p[zgp0]*(-rhoz + u0*ri0*(one - dtvh*ri0)) + + var_p[zgp1]*( rhoz + u1*ri1*(one - dtvh*ri1))) + / ( rhoz + u0*ri0*(one + dtvh*ri0)); + } } } |