diff options
author | schnetter <schnetter@6a3ddf76-46e1-4315-99d9-bc56cac1ef84> | 2004-07-06 15:53:55 +0000 |
---|---|---|
committer | schnetter <schnetter@6a3ddf76-46e1-4315-99d9-bc56cac1ef84> | 2004-07-06 15:53:55 +0000 |
commit | 65cf9444b03ca7feff561eada96c9c8f9cbf0b56 (patch) | |
tree | 4c8ac70c2aaef8b194180e08add256a00098753b /src | |
parent | 8e9335ee54ecc33e8b07be4ce6fd815b316f768c (diff) |
Use a fudge factor to remove singularities when grid points fall onto
an axis.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAnalyticBH/trunk@152 6a3ddf76-46e1-4315-99d9-bc56cac1ef84
Diffstat (limited to 'src')
-rw-r--r-- | src/Schwarzschild.c | 19 |
1 files changed, 12 insertions, 7 deletions
diff --git a/src/Schwarzschild.c b/src/Schwarzschild.c index 850bd9a..331b10a 100644 --- a/src/Schwarzschild.c +++ b/src/Schwarzschild.c @@ -9,6 +9,7 @@ @version $Id$ @@*/ +#include <math.h> #include <string.h> #include "cctk.h" @@ -68,14 +69,16 @@ void Schwarzschild(CCTK_ARGUMENTS) for (i = 0; i < npoints; i++) { /* Compute conformal factor */ - psi[i] = ( one + mass/(two*r[i]+1e-20)); + CCTK_REAL const rr = pow(pow(r[i], 4) + pow(epsilon, 4), 0.25); + + psi[i] = (one + mass / (two*rr)); if(make_conformal_derivs) { /* derivatives of psi / psi */ - r_squared = r[i]*r[i]; - r_cubed = r[i]*r_squared; - tmp = mass/(two*r_cubed*psi[i]+1e-20); + r_squared = rr*rr; + r_cubed = rr*r_squared; + tmp = mass / (two*r_cubed*psi[i]); psix[i] = -x[i]*tmp; @@ -84,7 +87,7 @@ void Schwarzschild(CCTK_ARGUMENTS) if(*conformal_state > 2) { - tmp = mass/(two*r_squared*r_cubed*psi[i]+1e-20); + tmp = mass/(two*r_squared*r_cubed*psi[i]); psixy[i] = three*x[i]*y[i]*tmp; psixz[i] = three*x[i]*z[i]*tmp; psiyz[i] = three*y[i]*z[i]*tmp; @@ -106,7 +109,9 @@ void Schwarzschild(CCTK_ARGUMENTS) { for (i = 0; i < npoints; i++) { - tmp = one + mass/(two*r[i]+1e-20); + CCTK_REAL const rr = pow(pow(r[i], 4) + pow(epsilon, 4), 0.25); + + tmp = one + mass / (two*rr); gxx[i] = tmp*tmp*tmp*tmp; gyy[i] = gxx[i]; gzz[i] = gxx[i]; @@ -124,7 +129,7 @@ void Schwarzschild(CCTK_ARGUMENTS) for (i = 0; i < npoints; i++) { - alp[i] = (2.*r[i] - mass)/(2.*r[i]+mass); + alp[i] = (two*r[i] - mass) / (two*r[i] + mass); } } |