Many changes.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@89 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 264 insertions, 58 deletions

diff --git a/src/AHFinder_shift.F b/src/AHFinder_shift.F
index 4029e57..c72c016 100644
--- a/src/AHFinder_shift.F
+++ b/src/AHFinder_shift.F
@@ -3,13 +3,40 @@
   @date      May 2000
   @author    Miguel Alcubierre
   @desc 
-             Set up horizon locking shift.
+             Set up horizon (area) locking shift.  The area
+             locking shift is obtained from the following
+             expression:
+
+                 a        __a
+             beta   =  b  \/ f
+
+
+             where f is the "horizon function", and where b is given by:
+
+
+                         /             ab               2 \
+             b  =  alpha | trK   -    K      d f d f / u  |
+                         \                    a   b       /
+
+                         / __2      __a__b              2 \
+                     /   | \/ f  -  \/ \/ f  d f d f / u  |
+                         \                    a   b       /
+
+
+             with "alpha" the lapse function and:
+
+
+              2     mn
+             u  =  g   d f d f
+                        m   n
+
   @enddesc 
 @@*/
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
+#include "cctk_DefineThorn.h"
 
       subroutine AHFinder_shift(CCTK_FARGUMENTS)
 
@@ -21,51 +48,68 @@
       DECLARE_CCTK_PARAMETERS
       DECLARE_CCTK_FUNCTIONS
 
-      logical shift_exp,shift_pow
-
       CCTK_INT i,j,k
+      CCTK_INT order_save
 
       CCTK_REAL det,idet
       CCTK_REAL gdxx,gdyy,gdzz,gdxy,gdxz,gdyz
       CCTK_REAL guxx,guyy,guzz,guxy,guxz,guyz
-      CCTK_REAL ndx,ndy,ndz
-      CCTK_REAL nux,nuy,nuz
-      CCTK_REAL aux
+      CCTK_REAL kdxx,kdyy,kdzz,kdxy,kdxz,kdyz
+      CCTK_REAL ddxxx,ddxyy,ddxzz,ddxxy,ddxxz,ddxyz
+      CCTK_REAL ddyxx,ddyyy,ddyzz,ddyxy,ddyxz,ddyyz
+      CCTK_REAL ddzxx,ddzyy,ddzzz,ddzxy,ddzxz,ddzyz
+      CCTK_REAL dfx,dfy,dfz,dfux,dfuy,dfuz
+      CCTK_REAL d2fxx,d2fyy,d2fzz,d2fxy,d2fxz,d2fyz
+      CCTK_REAL c2fxx,c2fyy,c2fzz,c2fxy,c2fxz,c2fyz
+      CCTK_REAL idx,idy,idz,idx2,idy2,idz2,idxy,idxz,idyz
+      CCTK_REAL zero,half,one,two
+      CCTK_REAL aux,T0,T1,T2,T3,T4,T5
 
+      CCTK_REAL, DIMENSION(nx,ny,nz) :: mask,betax_old,betay_old,betaz_old
 
-!     ***************************************
-!     ***   FIND OUT WHAT KIND OF SHIFT   ***
-!     ***************************************
 
-      shift_exp = CCTK_EQUALS(shift,'locking1')
-      shift_pow = CCTK_EQUALS(shift,'locking2')
+!     *******************
+!     ***   NUMBERS   ***
+!     *******************
 
-      if (shift_exp) then
-         write(*,*) 'AHFinder: Setting up exponential horizon locking shift.'
-      else if (shift_pow) then
-         write(*,*) 'AHFinder: Setting up power law horizon locking shift.'
-      else
-         call CCTK_WARN(0,"AHFinder: Unkown type of horizon locking shift")
-      end if
+      zero = 0.0D0
+      half = 0.5D0
+      one  = 1.0D0
+      two  = 2.0D0
 
-      write(*,*)
 
+!     *********************************
+!     ***   SAY WHAT WE ARE DOING   ***
+!     *********************************
 
-!     **************************************************
-!     ***   FIND HORIZON FUNCTION AND ITS GRADIENT   ***
-!     **************************************************
+      write(*,*)
+      write(*,*) 'AHFinder: Setting up horizon locking shift.'
+      write(*,*)
 
-      call AHFinder_fun(CCTK_FARGUMENTS)
-      call AHFinder_exp(CCTK_FARGUMENTS)
+      betax_old = betax
+      betay_old = betay
+      betaz_old = betaz
 
 
 !     *******************************
 !     ***   SET UP SHIFT VECTOR   ***
 !     *******************************
 
-      do k=1,nz
-         do j=1,ny
-            do i=1,nx
+      idx = half/dx
+      idy = half/dy
+      idz = half/dz
+
+      idxy = idx*idy
+      idxz = idx*idz
+      idyz = idy*idz
+
+      idx2 = one/dx**2
+      idy2 = one/dy**2
+      idz2 = one/dz**2
+
+      do k=2,nz-1
+         do j=2,ny-1
+            do i=2,nx-1
 
 !              Find spatial metric.
 
@@ -76,16 +120,21 @@
                gdxz = gxz(i,j,k)
                gdyz = gyz(i,j,k)
 
+!              Find extrinsic curvature.
+
+               kdxx = kxx(i,j,k)
+               kdyy = kyy(i,j,k)
+               kdzz = kzz(i,j,k)
+               kdxy = kxy(i,j,k)
+               kdxz = kxz(i,j,k)
+               kdyz = kyz(i,j,k)
+
 !              Find determinant of spatial metric.
 
-               det = gdxx*gdyy*gdzz + 2.0D0*gdxy*gdxz*gdyz
+               det = gdxx*gdyy*gdzz + two*gdxy*gdxz*gdyz
      .             - gdxx*gdyz**2 - gdyy*gdxz**2 - gdzz*gdxy**2
 
-               if (det.ne.0.0D0) then
-                  idet = 1.0D0/det
-               else
-                  idet = 0.0D0
-               end if
+               idet = one/det
 
 !              Find inverse spatial metric.
 
@@ -97,51 +146,208 @@
                guxz = idet*(gdxy*gdyz - gdyy*gdxz)
                guyz = idet*(gdxy*gdxz - gdxx*gdyz)
 
-!              Find covariant unit normal vector.
+!              Find first spatial derivatives of f.
+
+               dfx = idx*(ahfgrid(i+1,j,k) - ahfgrid(i-1,j,k))
+               dfy = idy*(ahfgrid(i,j+1,k) - ahfgrid(i,j-1,k))
+               dfz = idz*(ahfgrid(i,j,k+1) - ahfgrid(i,j,k-1))
+
+!              Raise indices in first derivatives.
+
+               dfux = guxx*dfx + guxy*dfy + guxz*dfz
+               dfuy = guxy*dfx + guyy*dfy + guyz*dfz
+               dfuz = guxz*dfx + guyz*dfy + guzz*dfz
+
+!              Find derivatives of metric using finite differences.
+
+               ddxxx = idx*(gxx(i+1,j,k) - gxx(i-1,j,k))
+               ddxyy = idx*(gyy(i+1,j,k) - gyy(i-1,j,k))
+               ddxzz = idx*(gzz(i+1,j,k) - gzz(i-1,j,k))
+               ddxxy = idx*(gxy(i+1,j,k) - gxy(i-1,j,k))
+               ddxxz = idx*(gxz(i+1,j,k) - gxz(i-1,j,k))
+               ddxyz = idx*(gyz(i+1,j,k) - gyz(i-1,j,k))
+
+               ddyxx = idy*(gxx(i,j+1,k) - gxx(i,j-1,k))
+               ddyyy = idy*(gyy(i,j+1,k) - gyy(i,j-1,k))
+               ddyzz = idy*(gzz(i,j+1,k) - gzz(i,j-1,k))
+               ddyxy = idy*(gxy(i,j+1,k) - gxy(i,j-1,k))
+               ddyxz = idy*(gxz(i,j+1,k) - gxz(i,j-1,k))
+               ddyyz = idy*(gyz(i,j+1,k) - gyz(i,j-1,k))
+
+               ddzxx = idz*(gxx(i,j,k+1) - gxx(i,j,k-1))
+               ddzyy = idz*(gyy(i,j,k+1) - gyy(i,j,k-1))
+               ddzzz = idz*(gzz(i,j,k+1) - gzz(i,j,k-1))
+               ddzxy = idz*(gxy(i,j,k+1) - gxy(i,j,k-1))
+               ddzxz = idz*(gxz(i,j,k+1) - gxz(i,j,k-1))
+               ddzyz = idz*(gyz(i,j,k+1) - gyz(i,j,k-1))
+
+!              Find second spatial derivatives of f.
+
+               d2fxx = (ahfgrid(i+1,j,k) - two*ahfgrid(i,j,k)
+     .               +  ahfgrid(i-1,j,k))*idx2
+               d2fyy = (ahfgrid(i,j+1,k) - two*ahfgrid(i,j,k)
+     .               +  ahfgrid(i,j-1,k))*idy2
+               d2fzz = (ahfgrid(i,j,k+1) - two*ahfgrid(i,j,k)
+     .               +  ahfgrid(i,j,k-1))*idz2
+
+               d2fxy = (ahfgrid(i+1,j+1,k) + ahfgrid(i-1,j-1,k)
+     .               -  ahfgrid(i+1,j-1,k) - ahfgrid(i-1,j+1,k))*idxy
+               d2fxz = (ahfgrid(i+1,j,k+1) + ahfgrid(i-1,j,k-1)
+     .               -  ahfgrid(i+1,j,k-1) - ahfgrid(i-1,j,k+1))*idxz
+               d2fyz = (ahfgrid(i,j+1,k+1) + ahfgrid(i,j-1,k-1)
+     .               -  ahfgrid(i,j+1,k-1) - ahfgrid(i,j-1,k+1))*idyz
+
+!              Find second covariant derivatives of f.
+
+               c2fxx = d2fxx - half*(dfux*ddxxx
+     .               + dfuy*(two*ddxxy - ddyxx)
+     .               + dfuz*(two*ddxxz - ddzxx))
 
-               aux = ahfgradn(i,j,k)
+               c2fyy = d2fyy - half*(dfuy*ddyyy
+     .               + dfux*(two*ddyxy - ddxyy)
+     .               + dfuz*(two*ddyyz - ddzyy))
 
-               if (aux.ne.0.0D0) then
-                  aux = 1.0D0/aux
-                  ndx = ahfgradx(i,j,k)*aux
-                  ndy = ahfgrady(i,j,k)*aux
-                  ndz = ahfgradz(i,j,k)*aux
+               c2fzz = d2fzz - half*(dfuz*ddzzz
+     .               + dfux*(two*ddzxz - ddxzz)
+     .               + dfuy*(two*ddzyz - ddyzz))
+
+               c2fxy = d2fxy - half*(dfux*ddyxx + dfuy*ddxyy
+     .               + dfuz*(ddxyz + ddyxz - ddzxy))
+
+               c2fxz = d2fxz - half*(dfux*ddzxx + dfuz*ddxzz
+     .               + dfuy*(ddxyz + ddzxy - ddyxz))
+
+               c2fyz = d2fyz - half*(dfuy*ddzyy + dfuz*ddyzz
+     .               + dfux*(ddyxz + ddzxy - ddxyz))
+
+!                                 2    /  m      \ - 1
+!              Find:   T0 =  1 / u  =  | d f d f |
+!                                      \      m  /
+
+               T0 = dfx*dfux + dfy*dfuy + dfz*dfuz
+
+               if (T0.gt.zero) then
+                  T0 = one/T0
                else
-                  ndx = 0.0D0
-                  ndy = 0.0D0
-                  ndz = 0.0D0
+                  betax(i,j,k) = zero
+                  betay(i,j,k) = zero
+                  betaz(i,j,k) = zero
+                  goto 10
                end if
 
-!              Find contravariant unit normal vector.
+!                            __2
+!              Find:  T1  =  \/ f
+ 
+               T1 = guxx*c2fxx + guyy*c2fyy + guzz*c2fzz
+     .            + two*(guxy*c2fxy + guxz*c2fxz + guyz*c2fyz)
 
-               nux = guxx*ndx + guxy*ndy + guxz*ndz
-               nuy = guxy*ndx + guyy*ndy + guyz*ndz
-               nuz = guxz*ndx + guyz*ndy + guzz*ndz
+!                            __ __     a   b     2
+!              Find:  T2  =  \/ \/ f  d f d f / u
+!                              a  b 
 
-!              Set up horizon locking shift.
+               T2 = c2fxx*dfux**2 + c2fyy*dfuy**2 + c2fzz*dfuz**2
+     .            + two*(c2fxy*dfux*dfuy + c2fxz*dfux*dfuz
+     .            + c2fyz*dfuy*dfuz)
 
-               if (shift_exp) then
+               T2 = T2*T0
 
-!                 Exponential shift.
+!              Find:  T3  =  trK
 
-                  aux = alp(i,j,k)*dexp(-ahfgrid(i,j,k))
+               T3 = guxx*kdxx + guyy*kdyy + guzz*kdzz
+     .            + two*(guxy*kdxy + guxz*kdxz + guyz*kdyz)
 
-               else if (shift_pow) then 
+!                                 a   b     2
+!              Find:  T4  =  K   d f d f / u
+!                             ab
 
-!                 Power law shift.
+               T4 = kdxx*dfux**2 + kdyy*dfuy**2 + kdzz*dfuz**2
+     .            + two*(kdxy*dfux*dfuy + kdxz*dfux*dfuz
+     .            + kdyz*dfuy*dfuz)
 
-                  aux = alp(i,j,k)*rhmax/(rhmax + ahfgrid(i,j,k))
+               T4 = T4*T0
 
-               end if
+!              Find shift.
+
+               aux = alp(i,j,k)*(T3 - T4)/(T1 - T2)
+
+               betax(i,j,k) = aux*dfux
+               betay(i,j,k) = aux*dfuy
+               betaz(i,j,k) = aux*dfuz
 
-               betax(i,j,k) = aux*nux
-               betay(i,j,k) = aux*nuy
-               betaz(i,j,k) = aux*nuz
+ 10            continue
 
             end do
          end do
       end do
 
+!     Boundaries.
+
+      if (CCTK_EQUALS(ahf_shiftbound,'extrap')) then
+
+!        X direction.
+
+         betax(1,:,:) = two*betax(2,:,:) - betax(3,:,:)
+         betay(1,:,:) = two*betay(2,:,:) - betay(3,:,:)
+         betaz(1,:,:) = two*betaz(2,:,:) - betaz(3,:,:)
+
+         betax(nx,:,:) = two*betax(nx-1,:,:) - betax(nx-2,:,:)
+         betay(nx,:,:) = two*betay(nx-1,:,:) - betay(nx-2,:,:)
+         betaz(nx,:,:) = two*betaz(nx-1,:,:) - betaz(nx-2,:,:)
+
+!        Y direction.
+
+         betax(:,1,:) = two*betax(:,2,:) - betax(:,3,:)
+         betay(:,1,:) = two*betay(:,2,:) - betay(:,3,:)
+         betaz(:,1,:) = two*betaz(:,2,:) - betaz(:,3,:)
+
+         betax(:,ny,:) = two*betax(:,ny-1,:) - betax(:,ny-2,:)
+         betay(:,ny,:) = two*betay(:,ny-1,:) - betay(:,ny-2,:)
+         betaz(:,ny,:) = two*betaz(:,ny-1,:) - betaz(:,ny-2,:)
+
+!        Z direction.
+
+         betax(:,:,1) = two*betax(:,:,2) - betax(:,:,3)
+         betay(:,:,1) = two*betay(:,:,2) - betay(:,:,3)
+         betaz(:,:,1) = two*betaz(:,:,2) - betaz(:,:,3)
+
+         betax(:,:,nz) = two*betax(:,:,nz-1) - betax(:,:,nz-2)
+         betay(:,:,nz) = two*betay(:,:,nz-1) - betay(:,:,nz-2)
+         betaz(:,:,nz) = two*betaz(:,:,nz-1) - betaz(:,:,nz-2)
+
+      end if
+
+
+!     ********************
+!     ***   EXCISION   ***
+!     ********************
+
+      if (ahf_shiftexcise.eq.1) then
+
+#ifdef DEVELOPMENT_SIMPLEEXCISION
+
+         order_save = order
+         order = 1
+
+         call SimpleExcision(_CCTK_FARGUMENTS,x,y,z,r,mask,
+     .      betax,betax_old,zero,one,'extrapolation')
+         call SimpleExcision(_CCTK_FARGUMENTS,x,y,z,r,mask,
+     .      betay,betay_old,zero,one,'extrapolation')
+         call SimpleExcision(_CCTK_FARGUMENTS,x,y,z,r,mask,
+     .      betaz,betaz_old,zero,one,'extrapolation')
+
+         order = order_save
+
+#else
+         call CCTK_WARN(0,"You have not compiled with SimpleExcision")
+#endif
+
+      end if
+
+
+!     ********************************************
+!     ***   APPLY SYMMETRIES AND SYNCHRONIZE   ***
+!     ********************************************
+
       call CartSymBCGroup(ierr,cctkGH,"einstein::shift")
       call CCTK_SyncGroup(cctkGH,"einstein::shift")