Adding subroutine to calculate horizon locking shift.

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

2 files changed, 155 insertions, 1 deletions

diff --git a/src/AHFinder_shift.F b/src/AHFinder_shift.F
new file mode 100644
index 0000000..4029e57
--- /dev/null
+++ b/src/AHFinder_shift.F
@@ -0,0 +1,154 @@
+/*@@
+  @file      AHFinder_shift.F
+  @date      May 2000
+  @author    Miguel Alcubierre
+  @desc 
+             Set up horizon locking shift.
+  @enddesc 
+@@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+      subroutine AHFinder_shift(CCTK_FARGUMENTS)
+
+      use AHFinder_dat
+
+      implicit none
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+      logical shift_exp,shift_pow
+
+      CCTK_INT i,j,k
+
+      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
+
+
+!     ***************************************
+!     ***   FIND OUT WHAT KIND OF SHIFT   ***
+!     ***************************************
+
+      shift_exp = CCTK_EQUALS(shift,'locking1')
+      shift_pow = CCTK_EQUALS(shift,'locking2')
+
+      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
+
+      write(*,*)
+
+
+!     **************************************************
+!     ***   FIND HORIZON FUNCTION AND ITS GRADIENT   ***
+!     **************************************************
+
+      call AHFinder_fun(CCTK_FARGUMENTS)
+      call AHFinder_exp(CCTK_FARGUMENTS)
+
+
+!     *******************************
+!     ***   SET UP SHIFT VECTOR   ***
+!     *******************************
+
+      do k=1,nz
+         do j=1,ny
+            do i=1,nx
+
+!              Find spatial metric.
+
+               gdxx = gxx(i,j,k)
+               gdyy = gyy(i,j,k)
+               gdzz = gzz(i,j,k)
+               gdxy = gxy(i,j,k)
+               gdxz = gxz(i,j,k)
+               gdyz = gyz(i,j,k)
+
+!              Find determinant of spatial metric.
+
+               det = gdxx*gdyy*gdzz + 2.0D0*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
+
+!              Find inverse spatial metric.
+
+               guxx = idet*(gdyy*gdzz - gdyz**2)
+               guyy = idet*(gdxx*gdzz - gdxz**2)
+               guzz = idet*(gdxx*gdyy - gdxy**2)
+
+               guxy = idet*(gdxz*gdyz - gdzz*gdxy)
+               guxz = idet*(gdxy*gdyz - gdyy*gdxz)
+               guyz = idet*(gdxy*gdxz - gdxx*gdyz)
+
+!              Find covariant unit normal vector.
+
+               aux = ahfgradn(i,j,k)
+
+               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
+               else
+                  ndx = 0.0D0
+                  ndy = 0.0D0
+                  ndz = 0.0D0
+               end if
+
+!              Find contravariant unit normal vector.
+
+               nux = guxx*ndx + guxy*ndy + guxz*ndz
+               nuy = guxy*ndx + guyy*ndy + guyz*ndz
+               nuz = guxz*ndx + guyz*ndy + guzz*ndz
+
+!              Set up horizon locking shift.
+
+               if (shift_exp) then
+
+!                 Exponential shift.
+
+                  aux = alp(i,j,k)*dexp(-ahfgrid(i,j,k))
+
+               else if (shift_pow) then 
+
+!                 Power law shift.
+
+                  aux = alp(i,j,k)*rhmax/(rhmax + ahfgrid(i,j,k))
+
+               end if
+
+               betax(i,j,k) = aux*nux
+               betay(i,j,k) = aux*nuy
+               betaz(i,j,k) = aux*nuz
+
+            end do
+         end do
+      end do
+
+      call CartSymBCGroup(ierr,cctkGH,"einstein::shift")
+      call CCTK_SyncGroup(cctkGH,"einstein::shift")
+
+
+!     ***************
+!     ***   END   ***
+!     ***************
+
+      end subroutine AHFinder_shift
+
diff --git a/src/make.code.defn b/src/make.code.defn
index dadb6cb..15eb506 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -4,6 +4,6 @@
 SRCS += AHFinder_dat.F AHFinder.F \
         AHFinder_fun.F AHFinder_exp.F AHFinder_int.F \
         AHFinder_min.F AHFinder_pow.F AHFinder_leg.F \
-        AHFinder_flow.F AHFinder_mask.F \
+        AHFinder_flow.F AHFinder_mask.F AHFinder_shift.F \
 	AHFinder_find3.F AHFinder_calcsigma.F \
         AHFinder_initguess.F AHFinder_gau.F SetSym.F