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, 266 insertions, 0 deletions

diff --git a/src/NoExcision_mod.F90 b/src/NoExcision_mod.F90
new file mode 100644
index 0000000..646d335
--- /dev/null
+++ b/src/NoExcision_mod.F90
@@ -0,0 +1,266 @@
+! $Header$
+
+#include "cctk.h"
+
+module NoExcision_mod
+
+  implicit none
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  CCTK_REAL, parameter :: one = 1.0_wp
+  CCTK_REAL, parameter :: two = 2.0_wp
+  CCTK_INT :: infnorm_handle, sum_handle
+  CCTK_INT :: sym_selector
+  CCTK_INT :: loop_counter
+  CCTK_REAL, dimension(16) :: infnormresid
+  CCTK_REAL, dimension(16) :: delta_new, delta_old, alpha, beta, absres
+  character(len=18), dimension(16) :: red_names = (/ 'NoExcision::redgxx', &
+                                                     'NoExcision::redgxy', &
+                                                     'NoExcision::redgxz', &
+                                                     'NoExcision::redgyy', &
+                                                     'NoExcision::redgyz', &
+                                                     'NoExcision::redgzz', &
+                                                     'NoExcision::redkxx', &
+                                                     'NoExcision::redkxy', &
+                                                     'NoExcision::redkxz', &
+                                                     'NoExcision::redkyy', &
+                                                     'NoExcision::redkyz', &
+                                                     'NoExcision::redkzz', &
+                                                     'NoExcision::red', &
+                                                     'NoExcision::redx', &
+                                                     'NoExcision::redy', &
+                                                     'NoExcision::redz' /)
+  character(len=5), dimension(16) :: var_names = (/ '  gxx', '  gxy', &
+                                                    '  gxz', '  gyy', &
+                                                    '  gyz', '  gzz', &
+                                                    '  kxx', '  kxy', &
+                                                    '  kxz', '  kyy', &
+                                                    '  kyz', '  kzz', &
+                                                    '  alp', 'betax', &
+                                                    'betay', 'betaz' /)
+  integer :: nx, ny, nz
+  integer :: ii
+  logical, dimension(16) :: cont = (/ (.true., ii=1, 16) /)
+
+  contains
+
+    ! This is the basic residual calculation routine...
+
+    subroutine residual (val, mask, res, si, order )
+      CCTK_REAL, dimension(:,:,:), intent(in) :: val
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+      CCTK_REAL, dimension(:,:,:), intent(out) :: res
+      CCTK_REAL, intent(in) :: si
+      CCTK_INT, intent(in) :: order
+      CCTK_REAL, parameter :: six = 6.0_wp
+      CCTK_REAL, parameter :: sixth = 1.0_wp / six
+      CCTK_REAL, parameter :: four = 4.0_wp
+      CCTK_REAL, parameter :: eighteen = 18.0_wp
+      CCTK_REAL, parameter :: eighteenth = 1.0_wp / eighteen
+      CCTK_REAL, parameter :: fifteen = 15.0_wp
+      CCTK_REAL, parameter :: sixty = 60.0_wp
+      CCTK_REAL, parameter :: sixtieth = 1.0_wp / sixty
+      integer :: i, j, k
+  
+      select case (order)
+      case (2) 
+
+        do k = 2, nz-1
+          do j = 2, ny-1
+            do i = 2, nx-1
+    
+              if (mask(i,j,k) > 0) then
+  
+                res(i,j,k) = val(i+1,j,k) + val(i-1,j,k) + val(i,j+1,k) + &
+                             val(i,j-1,k) + val(i,j,k+1) + val(i,j,k-1) - &
+                             six * val(i,j,k)
+                res(i,j,k) = sign(one,si) * res(i,j,k)
+  
+              end if
+    
+            end do
+          end do
+        end do
+
+      case (4)
+
+        do k = 3, nz-2
+          do j = 3, ny-2
+            do i = 3, nx-2
+
+              if (mask(i,j,k) > 0) then
+  
+                res(i,j,k) = - ( val(i+2,j,k) + val(i-2,j,k) + &
+                                 val(i,j+2,k) + val(i,j-2,k) + &
+                                 val(i,j,k+2) + val(i,j,k-2) ) + &
+                        four * ( val(i+1,j,k) + val(i-1,j,k) + &
+                                 val(i,j+1,k) + val(i,j-1,k) + &
+                                 val(i,j,k+1) + val(i,j,k-1) ) - &
+                      eighteen * val(i,j,k)
+                res(i,j,k) = sign(one,si) * res(i,j,k)
+
+              end if
+
+            end do
+          end do
+        end do
+
+      case (6)
+
+        do k = 4, nz-3
+          do j = 4, ny-3
+            do i = 4, nx-3
+
+              if (mask(i,j,k) > 0) then
+
+                res(i,j,k) =   ( val(i+3,j,k) + val(i-3,j,k) + &
+                                 val(i,j+3,k) + val(i,j-3,k) + &
+                                 val(i,j,k+3) + val(i,j,k-3) ) - &
+                         six * ( val(i+2,j,k) + val(i-2,j,k) + &
+                                 val(i,j+2,k) + val(i,j-2,k) + &
+                                 val(i,j,k+2) + val(i,j,k-2) ) + &
+                     fifteen * ( val(i+1,j,k) + val(i-1,j,k) + &
+                                 val(i,j+1,k) + val(i,j-1,k) + &
+                                 val(i,j,k+1) + val(i,j,k-1) ) - &
+                         sixty * val(i,j,k)
+                res(i,j,k) = sign(one,si) * res(i,j,k)
+
+              end if
+
+            end do
+          end do
+        end do
+
+      case default
+
+        call CCTK_WARN ( 0, 'Internal error: Order out of range' )
+  
+      end select
+
+    end subroutine residual
+
+
+    ! This is the routine calculating residuals for all variables, that
+    ! have not yet converged.
+
+    subroutine residual_all ( v1, v2, v3, v4, v5, v6, v7, v8, &
+                               v9, v10, v11, v12, v13, v14, v15, v16, &
+                               r1, r2, r3, r4, r5, r6, r7, r8, &
+                               r9, r10, r11, r12, r13, r14, r15, r16, &
+                               mask, si, order )
+
+      CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+                                                 v7, v8, v9, v10, v11, v12, &
+                                                 v13, v14, v15, v16
+      CCTK_REAL, dimension(:,:,:), intent(out) :: r1, r2, r3, r4, r5, r6, &
+                                                  r7, r8, r9, r10, r11, r12, &
+                                                  r13, r14, r15, r16
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+      CCTK_REAL, intent(in) :: si
+      CCTK_INT, intent(in) :: order
+
+      if ( cont(1) ) call residual ( v1, mask, r1, si, order )
+      if ( cont(2) ) call residual ( v2, mask, r2, si, order )
+      if ( cont(3) ) call residual ( v3, mask, r3, si, order )
+      if ( cont(4) ) call residual ( v4, mask, r4, si, order )
+      if ( cont(5) ) call residual ( v5, mask, r5, si, order )
+      if ( cont(6) ) call residual ( v6, mask, r6, si, order )
+      if ( cont(7) ) call residual ( v7, mask, r7, si, order )
+      if ( cont(8) ) call residual ( v8, mask, r8, si, order )
+      if ( cont(9) ) call residual ( v9, mask, r9, si, order )
+      if ( cont(10) ) call residual ( v10, mask, r10, si, order )
+      if ( cont(11) ) call residual ( v11, mask, r11, si, order )
+      if ( cont(12) ) call residual ( v12, mask, r12, si, order )
+      if ( cont(13) ) call residual ( v13, mask, r13, si, order )
+      if ( cont(14) ) call residual ( v14, mask, r14, si, order )
+      if ( cont(15) ) call residual ( v15, mask, r15, si, order )
+      if ( cont(16) ) call residual ( v16, mask, r16, si, order )
+
+    end subroutine residual_all
+
+
+    ! This routine multiplies two variables. Only does it according
+    ! to the mask and only for variables that have not converged yet.
+
+    subroutine multiply ( u1, u2, u3, u4, u5, u6, u7, u8, &
+                          u9, u10, u11, u12, u13, u14, u15, u16, &
+                          v1, v2, v3, v4, v5, v6, v7, v8, &
+                          v9, v10, v11, v12, v13, v14, v15, v16, &
+                          r1, r2, r3, r4, r5, r6, r7, r8, &
+                          r9, r10, r11, r12, r13, r14, r15, r16, &
+                          mask )
+
+      CCTK_REAL, dimension(:,:,:), intent(in) :: u1, u2, u3, u4, u5, u6, &
+                                                 u7, u8, u9, u10, u11, u12, &
+                                                 u13, u14, u15, u16
+      CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+                                                 v7, v8, v9, v10, v11, v12, &
+                                                 v13, v14, v15, v16
+      CCTK_REAL, dimension(:,:,:), intent(out) :: r1, r2, r3, r4, r5, r6, &
+                                                  r7, r8, r9, r10, r11, r12, &
+                                                  r13, r14, r15, r16
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+
+      if ( cont(1) ) where ( mask > 0 ) r1 = u1*v1
+      if ( cont(2) ) where ( mask > 0 ) r2 = u2*v2
+      if ( cont(3) ) where ( mask > 0 ) r3 = u3*v3
+      if ( cont(4) ) where ( mask > 0 ) r4 = u4*v4
+      if ( cont(5) ) where ( mask > 0 ) r5 = u5*v5
+      if ( cont(6) ) where ( mask > 0 ) r6 = u6*v6
+      if ( cont(7) ) where ( mask > 0 ) r7 = u7*v7
+      if ( cont(8) ) where ( mask > 0 ) r8 = u8*v8
+      if ( cont(9) ) where ( mask > 0 ) r9 = u9*v9
+      if ( cont(10) ) where ( mask > 0 ) r10 = u10*v10
+      if ( cont(11) ) where ( mask > 0 ) r11 = u11*v11
+      if ( cont(12) ) where ( mask > 0 ) r12 = u12*v12
+      if ( cont(13) ) where ( mask > 0 ) r13 = u13*v13
+      if ( cont(14) ) where ( mask > 0 ) r14 = u14*v14
+      if ( cont(15) ) where ( mask > 0 ) r15 = u15*v15
+      if ( cont(16) ) where ( mask > 0 ) r16 = u16*v16
+
+    end subroutine multiply
+                  
+
+    ! This routine takes the product of a scalar and a variable and
+    ! adds it to the product of another scalar and variable. Does it
+    ! according to the mask and only for variables that have not yet
+    ! converged.
+
+    subroutine multiply_sum ( u1, u2, u3, u4, u5, u6, u7, u8, &
+                              u9, u10, u11, u12, u13, u14, u15, u16, &
+                              v1, v2, v3, v4, v5, v6, v7, v8, &
+                              v9, v10, v11, v12, v13, v14, v15, v16, &
+                              c1, c2, mask )
+
+      CCTK_REAL, dimension(:,:,:), intent(inout) :: u1, u2, u3, u4, u5, u6, &
+                                                    u7, u8, u9, u10, u11, u12, &
+                                                    u13, u14, u15, u16
+      CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+                                                 v7, v8, v9, v10, v11, v12, &
+                                                 v13, v14, v15, v16
+      CCTK_REAL, dimension(16), intent(in) :: c1, c2
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+
+      
+      if ( cont(1) ) where ( mask > 0 ) u1 = c1(1)*u1 + c2(1)*v1
+      if ( cont(2) ) where ( mask > 0 ) u2 = c1(2)*u2 + c2(2)*v2
+      if ( cont(3) ) where ( mask > 0 ) u3 = c1(3)*u3 + c2(3)*v3
+      if ( cont(4) ) where ( mask > 0 ) u4 = c1(4)*u4 + c2(4)*v4
+      if ( cont(5) ) where ( mask > 0 ) u5 = c1(5)*u5 + c2(5)*v5
+      if ( cont(6) ) where ( mask > 0 ) u6 = c1(6)*u6 + c2(6)*v6
+      if ( cont(7) ) where ( mask > 0 ) u7 = c1(7)*u7 + c2(7)*v7
+      if ( cont(8) ) where ( mask > 0 ) u8 = c1(8)*u8 + c2(8)*v8
+      if ( cont(9) ) where ( mask > 0 ) u9 = c1(9)*u9 + c2(9)*v9
+      if ( cont(10) ) where ( mask > 0 ) u10 = c1(10)*u10 + c2(10)*v10
+      if ( cont(11) ) where ( mask > 0 ) u11 = c1(11)*u11 + c2(11)*v11
+      if ( cont(12) ) where ( mask > 0 ) u12 = c1(12)*u12 + c2(12)*v12
+      if ( cont(13) ) where ( mask > 0 ) u13 = c1(13)*u13 + c2(13)*v13
+      if ( cont(14) ) where ( mask > 0 ) u14 = c1(14)*u14 + c2(14)*v14
+      if ( cont(15) ) where ( mask > 0 ) u15 = c1(15)*u15 + c2(15)*v15
+      if ( cont(16) ) where ( mask > 0 ) u16 = c1(16)*u16 + c2(16)*v16
+
+    end subroutine multiply_sum
+
+
+end module NoExcision_mod