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

2 files changed, 60 insertions, 23 deletions

diff --git a/param.ccl b/param.ccl
index 3debb67..1fe7342 100644
--- a/param.ccl
+++ b/param.ccl
@@ -10,6 +10,8 @@ INT num_regions "Number of no-excision regions"
   0:10 :: ""
 } 0
 
+
+
 BOOLEAN overwrite_geometry[10] "Set the geometry to Minkowski"
 {
 } "yes"
@@ -19,6 +21,8 @@ REAL Minkowski_scale[10] "Scaling factor for Minkowski"
   (*:*) :: "Choose 1 for true Minkowski"
 } 1.0
 
+
+
 BOOLEAN overwrite_lapse[10] "Set the lapse to one"
 {
 } "yes"
@@ -28,10 +32,14 @@ REAL lapse_scale[10] "Scaling factor for lapse"
   (*:*) :: "Choose 1 for geodesic slicing, 0 to halt evolution"
 } 1.0
 
+
+
 BOOLEAN overwrite_shift[10] "Set the shift to zero"
 {
 } "yes"
 
+
+
 REAL centre_x[10] "X-coordinate of the centre of the region"
 {
   (*:*) :: ""
@@ -51,3 +59,8 @@ REAL radius[10] "Radius of the region"
 {
   0.0:*) :: ""
 } 1.0
+
+REAL smoothing_zone_width[10] "Width of smoothing zone inside the region"
+{
+  0.0:*) :: ""
+} 0.0
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