Changes in methodology for filling in the interior of Cook-Pfeiffer data as

used in the "turducken" paper. This includes some spline blending from
Erik and the most recent method of filling in the interior using an elliptic
equation with the "good" points as boundary values. This method uses a
conjugate gradient solver (built in), that should be robust (i.e. it should
always work) but may not be the fastest. However, since this is only done on
initial data this shouldn't be an issue. The elliptic method currently support
a second order, fourth order and sixth order equation, that will give C0, C1
and C2 solutions, respectively, across the boundary of the fill in region.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/NoExcision/trunk@7 4ec1db94-0e4f-0410-ada3-8bed251432c9

1 files changed, 150 insertions, 0 deletions

diff --git a/src/smooth.F90 b/src/smooth.F90
new file mode 100644
index 0000000..d6df156
--- /dev/null
+++ b/src/smooth.F90
@@ -0,0 +1,150 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_Smooth (CCTK_ARGUMENTS)
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero=0
+  CCTK_REAL :: mask (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz
+  integer   :: n
+  integer   :: iter
+  integer   :: ierr
+  
+  do iter = 1, smoothing_iterations
+     
+     do n = 1, num_regions
+        
+        cx = centre_x(n)
+        cy = centre_y(n)
+        cz = centre_z(n)
+        if (CCTK_EQUALS(region_shape(n), "sphere")) then
+           radx = radius(n)
+           rady = radius(n)
+           radz = radius(n)
+        else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+           radx = radius_x(n)
+           rady = radius_y(n)
+           radz = radius_z(n)
+        else
+           call CCTK_WARN (0, "internal error")
+        end if
+        
+        mask = 1 - sqrt (((x - cx) / radx)**2 + &
+             &           ((y - cy) / rady)**2 + &
+             &           ((z - cz) / radz)**2)
+        where (mask <= 0)
+           mask = 0             ! outside
+        else where (mask >= smoothing_zone_width(n))
+           mask = 1             ! far inside
+        elsewhere
+           mask = mask / smoothing_zone_width(n) ! a bit inside
+        end where
+        
+        if (overwrite_geometry(n) /= 0) then
+           
+           if (conformal_state >= 1) then
+              call smooth (psi)
+           end if
+           if (conformal_state >= 2) then
+              call smooth (psix)
+              call smooth (psiy)
+              call smooth (psiz)
+           end if
+           if (conformal_state >= 3) then
+              call smooth (psixx)
+              call smooth (psixy)
+              call smooth (psixz)
+              call smooth (psiyy)
+              call smooth (psiyz)
+              call smooth (psizz)
+           end if
+           
+           call smooth (gxx)
+           call smooth (gxy)
+           call smooth (gxz)
+           call smooth (gyy)
+           call smooth (gyz)
+           call smooth (gzz)
+           call smooth (kxx)
+           call smooth (kxy)
+           call smooth (kxz)
+           call smooth (kyy)
+           call smooth (kyz)
+           call smooth (kzz)
+           
+        end if
+        
+        if (overwrite_lapse(n) /= 0) then
+           
+           call smooth (alp)
+           
+        end if
+        
+        if (overwrite_shift(n) /= 0) then
+           
+           if (shift_state /= 0) then
+              call smooth (betax)
+              call smooth (betay)
+              call smooth (betaz)
+           end if
+           
+        end if
+        
+     end do
+     
+     if (overwrite_geometry(n) /= 0) then
+        
+        if (conformal_state >= 1) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac")
+        end if
+        if (conformal_state >= 2) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_1derivs")
+        end if
+        if (conformal_state >= 3) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_2derivs")
+        end if
+        
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::metric")
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::curv")
+        
+     end if
+     
+     if (overwrite_lapse(n) /= 0) then
+        
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::lapse")
+        
+     end if
+     
+     if (overwrite_shift(n) /= 0) then
+        
+        if (shift_state /= 0) then
+           call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::shift")
+        end if
+        
+     end if
+     
+  end do
+  
+contains
+  
+  subroutine smooth (val)
+    CCTK_REAL, intent(inout) :: val(:,:,:)
+    
+    where (mask > 0)
+       val =  (1 - mask * smoothing_factor) * val &
+            & + (mask * smoothing_factor / 6) &
+            &   * (+ eoshift(val, +1, dim=1) + eoshift(val, -1, dim=1) &
+            &      + eoshift(val, +1, dim=2) + eoshift(val, -1, dim=2) &
+            &      + eoshift(val, +1, dim=3) + eoshift(val, -1, dim=3))
+    end where
+  end subroutine smooth
+  
+end subroutine NoExcision_Smooth