diff options
author | schnetter <schnetter@4ec1db94-0e4f-0410-ada3-8bed251432c9> | 2004-10-12 19:26:47 +0000 |
---|---|---|
committer | schnetter <schnetter@4ec1db94-0e4f-0410-ada3-8bed251432c9> | 2004-10-12 19:26:47 +0000 |
commit | 68f2a97a76104b96344e1a1a1c437a38a83db00b (patch) | |
tree | cf63acc473b3888737715525ab5515cc5ef52688 /src | |
parent | 5d921c7d6457f0392fec9ecee61b7e4372937cf3 (diff) |
Make it possible to smooth the transition from the exact data to flat
space
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/NoExcision/trunk@4 4ec1db94-0e4f-0410-ada3-8bed251432c9
Diffstat (limited to 'src')
-rw-r--r-- | src/overwrite.F90 | 70 |
1 files changed, 47 insertions, 23 deletions
diff --git a/src/overwrite.F90 b/src/overwrite.F90 index 87cd534..6c6089f 100644 --- a/src/overwrite.F90 +++ b/src/overwrite.F90 @@ -11,7 +11,9 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS DECLARE_CCTK_PARAMETERS - CCTK_REAL :: cx, cy, cz, rad + CCTK_REAL, parameter :: zero=0 + CCTK_REAL :: dist2 (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3)) + CCTK_REAL :: cx, cy, cz, rad, width integer :: n do n = 1, num_regions @@ -20,23 +22,26 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS) cy = centre_y(n) cz = centre_z(n) rad = radius(n) + width = smoothing_zone_width(n) + + dist2 = (x - cx)**2 + (y - cy)**2 + (z - cz)**2 if (overwrite_geometry(n) /= 0) then if (conformal_state >= 1) then - where ((x - cx)**2 + (y - cy)**2 + (z - cz)**2 <= rad**2) + where (dist2 <= rad**2) psi = 1 end where end if if (conformal_state >= 2) then - where ((x - cx)**2 + (y - cy)**2 + (z - cz)**2 <= rad**2) + where (dist2 <= rad**2) psix = 0 psiy = 0 psiz = 0 end where end if if (conformal_state >= 3) then - where ((x - cx)**2 + (y - cy)**2 + (z - cz)**2 <= rad**2) + where (dist2 <= rad**2) psixx = 0 psixy = 0 psixz = 0 @@ -46,27 +51,27 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS) end where end if - where ((x - cx)**2 + (y - cy)**2 + (z - cz)**2 <= rad**2) - gxx = Minkowski_scale(n) - gxy = 0 - gxz = 0 - gyy = Minkowski_scale(n) - gyz = 0 - gzz = Minkowski_scale(n) - kxx = 0 - kxy = 0 - kxz = 0 - kyy = 0 - kyz = 0 - kzz = 0 + where (dist2 <= rad**2) + gxx = smooth (gxx, Minkowski_scale(n), dist2) + gxy = smooth (gxy, zero , dist2) + gxz = smooth (gxz, zero , dist2) + gyy = smooth (gyy, Minkowski_scale(n), dist2) + gyz = smooth (gyz, zero , dist2) + gzz = smooth (gzz, Minkowski_scale(n), dist2) + kxx = smooth (kxx, zero , dist2) + kxy = smooth (kxy, zero , dist2) + kxz = smooth (kxz, zero , dist2) + kyy = smooth (kyy, zero , dist2) + kyz = smooth (kyz, zero , dist2) + kzz = smooth (kzz, zero , dist2) end where end if if (overwrite_lapse(n) /= 0) then - where ((x - cx)**2 + (y - cy)**2 + (z - cz)**2 <= rad**2) - alp = lapse_scale(n) + where (dist2 <= rad**2) + alp = smooth (alp, lapse_scale(n), dist2) end where end if @@ -74,10 +79,10 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS) if (overwrite_shift(n) /= 0) then if (shift_state /= 0) then - where ((x - cx)**2 + (y - cy)**2 + (z - cz)**2 <= rad**2) - betax = 0 - betay = 0 - betaz = 0 + where (dist2 <= rad**2) + betax = smooth (betax, zero, dist2) + betay = smooth (betay, zero, dist2) + betaz = smooth (betaz, zero, dist2) end where end if @@ -85,4 +90,23 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS) end do +contains + + elemental function smooth (oldval, goodval, dist2) result (newval) + CCTK_REAL :: newval + CCTK_REAL, intent(in) :: oldval + CCTK_REAL, intent(in) :: goodval + CCTK_REAL, intent(in) :: dist2 + + if (dist2 <= (rad-width)**2) then + newval = goodval + else if (width == 0 .or. dist2 >= rad**2) then + newval = oldval + else + newval = goodval & + & + (oldval - goodval) * (sqrt(dist2) - (rad-width)) / width + end if + + end function smooth + end subroutine NoExcision_Overwrite |