Take the static conformal factor into account.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/CalcK/trunk@5 f47d718b-0e4f-0410-8445-f2f96c8ccefb

6 files changed, 212 insertions, 22 deletions

diff --git a/interface.ccl b/interface.ccl
index 141167a..d4dc0ca 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -3,7 +3,7 @@
 
 IMPLEMENTS: CalcK
 
-INHERITS: ADMBase
+INHERITS: ADMBase StaticConformal
 
 CCTK_REAL metric_prev TYPE=gf
 {
diff --git a/par/calck.par b/par/calck.par
index bb3934f..1fc1926 100644
--- a/par/calck.par
+++ b/par/calck.par
@@ -6,8 +6,6 @@ Cactus::cctk_timer_output  = full
 
 Cactus::cctk_itlast = 0
 
-ActiveThorns = "always"
-
 ActiveThorns = "PUGH PUGHInterp PUGHReduce PUGHSlab"
 driver::global_nx = 80
 driver::global_ny = 80
@@ -19,16 +17,18 @@ ActiveThorns = "IOUtil"
 IO::out_dir = $parfile
 
 IO::filereader_ID_files = "
-	gxx_3d gxy_3d gxz_3d gyy_3d gyz_3d gzz_3d
-	alp_3d
-	betax_3d betay_3d betaz_3d
-        gxx_prev_3d gxy_prev_3d gxz_prev_3d gyy_prev_3d gyz_prev_3d gzz_prev_3d
-        gxx_next_3d gxy_next_3d gxz_next_3d gyy_next_3d gyz_next_3d gzz_next_3d
+	psi
+	gxx gxy gxz gyy gyz gzz
+	alp
+	betax betay betaz
+        gxx_prev gxy_prev gxz_prev gyy_prev gyz_prev gzz_prev
+        gxx_next gxy_next gxz_next gyy_next gyz_next gzz_next
 "
 IO::filereader_ID_vars = "
+	StaticConformal::confac
+	ADMBase::metric[cctk_iteration=101]
 	ADMBase::lapse[cctk_iteration=101]
 	ADMBase::shift[cctk_iteration=101]
-	ADMBase::metric[cctk_iteration=101]
 	CalcK::metric_prev
 	CalcK::metric_next
 "
@@ -43,7 +43,11 @@ CartGrid3D::avoid_origin = no
 
 ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConformal"
 
-ADMBase::initial_shift = "zero"
+ADMBase::metric_type = "static conformal"
+StaticConformal::conformal_storage = "factor"
+
+ADMBase::initial_data = none
+ADMBase::initial_shift = zero
 
 ActiveThorns = "CalcK"
 CalcK::calc_extcurv = yes
@@ -60,6 +64,7 @@ IOBasic::outScalar_every      = 1
 IOBasic::outScalar_reductions = "minimum maximum mean norm1 norm2 norm_inf"
 IOBasic::outScalar_style      = gnuplot
 IOBasic::outScalar_vars       = "
+	StaticConformal::confac
         ADMBase::metric
         ADMBase::curv
         ADMBase::lapse
@@ -70,6 +75,7 @@ ActiveThorns = "IOASCII"
 IOASCII::out1D_every = 1
 IOASCII::out1D_style = "gnuplot f(x)"
 IOASCII::out1D_vars  = "
+	StaticConformal::confac
         ADMBase::metric
         ADMBase::curv
         ADMBase::lapse
@@ -79,4 +85,9 @@ IOASCII::out1D_vars  = "
 IOHDF5::out_every = 1
 IOHDF5::out_vars  = "
         ADMBase::curv
+	StaticConformal::confac
+        ADMBase::metric
+        ADMBase::curv
+        ADMBase::lapse
+        ADMBase::shift
 "
diff --git a/par/copy_to_next.par b/par/copy_to_next.par
index ec68881..7cb1dea 100644
--- a/par/copy_to_next.par
+++ b/par/copy_to_next.par
@@ -6,8 +6,6 @@ Cactus::cctk_timer_output  = full
 
 Cactus::cctk_itlast = 0
 
-ActiveThorns = "always"
-
 ActiveThorns = "PUGH PUGHInterp PUGHReduce PUGHSlab"
 driver::global_nx = 80
 driver::global_ny = 80
@@ -19,7 +17,7 @@ ActiveThorns = "IOUtil"
 IO::out_dir = $parfile
 
 IO::filereader_ID_files = "
-	gxx_3d gxy_3d gxz_3d gyy_3d gyz_3d gzz_3d
+	gxx gxy gxz gyy gyz gzz
 "
 IO::filereader_ID_vars = "
 	ADMBase::metric[cctk_iteration=102]
@@ -35,7 +33,11 @@ CartGrid3D::avoid_origin = no
 
 ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConformal"
 
-ADMBase::initial_shift = "zero"
+ADMBase::metric_type = "static conformal"
+StaticConformal::conformal_storage = "factor"
+
+ADMBase::initial_data = none
+ADMBase::initial_shift = zero
 
 ActiveThorns = "CalcK"
 CalcK::copy_to_next = yes
diff --git a/par/copy_to_prev.par b/par/copy_to_prev.par
index 85113f0..f677fd9 100644
--- a/par/copy_to_prev.par
+++ b/par/copy_to_prev.par
@@ -6,8 +6,6 @@ Cactus::cctk_timer_output  = full
 
 Cactus::cctk_itlast = 0
 
-ActiveThorns = "always"
-
 ActiveThorns = "PUGH PUGHInterp PUGHReduce PUGHSlab"
 driver::global_nx = 80
 driver::global_ny = 80
@@ -19,7 +17,7 @@ ActiveThorns = "IOUtil"
 IO::out_dir = $parfile
 
 IO::filereader_ID_files = "
-	gxx_3d gxy_3d gxz_3d gyy_3d gyz_3d gzz_3d
+	gxx gxy gxz gyy gyz gzz
 "
 IO::filereader_ID_vars = "
 	ADMBase::metric[cctk_iteration=100]
@@ -35,7 +33,11 @@ CartGrid3D::avoid_origin = no
 
 ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConformal"
 
-ADMBase::initial_shift = "zero"
+ADMBase::metric_type = "static conformal"
+StaticConformal::conformal_storage = "factor"
+
+ADMBase::initial_data = none
+ADMBase::initial_shift = zero
 
 ActiveThorns = "CalcK"
 CalcK::copy_to_prev = yes
diff --git a/par/ih.par b/par/ih.par
new file mode 100644
index 0000000..baaf839
--- /dev/null
+++ b/par/ih.par
@@ -0,0 +1,168 @@
+# $Header$
+
+Cactus::cctk_run_title     = "Isolated Horizon"
+Cactus::cctk_full_warnings = yes
+Cactus::cctk_timer_output  = full
+
+Cactus::cctk_itlast = 0
+
+ActiveThorns = "LocalInterp AEILocalInterp"
+
+ActiveThorns = "Slab"
+
+ActiveThorns = "SphericalSurface"
+SphericalSurface::nsurfaces = 1
+SphericalSurface::symmetric_z[0] = yes
+SphericalSurface::ntheta[0]      = 11
+
+ActiveThorns = "TGRtensor"
+
+ActiveThorns = "PUGH PUGHInterp PUGHReduce PUGHSlab"
+driver::global_nx = 80
+driver::global_ny = 80
+driver::global_nz = 40
+
+ActiveThorns = "NaNChecker"
+
+ActiveThorns = "IOUtil"
+IO::out_dir = $parfile
+
+IO::filereader_ID_files = "
+	psi
+	gxx gxy gxz gyy gyz gzz
+	kxx kxy kxz kyy kyz kzz
+	alp
+	betax betay betaz
+"
+IO::filereader_ID_vars = "
+	StaticConformal::confac
+	ADMBase::lapse[cctk_iteration=101]
+	ADMBase::curv
+	ADMBase::shift[cctk_iteration=101]
+	ADMBase::metric[cctk_iteration=101]
+"
+
+ActiveThorns = "IOHDF5 IOHDF5Util"
+
+ActiveThorns = "Boundary CartGrid3D CoordBase"
+CartGrid3D::type         = byspacing
+CartGrid3D::domain       = bitant
+CartGrid3D::dxyz         = 0.08
+CartGrid3D::avoid_origin = no
+
+ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConformal"
+
+ADMBase::metric_type = "static conformal"
+StaticConformal::conformal_storage = "factor"
+
+ADMBase::initial_data = none
+ADMBase::initial_shift = zero
+
+ActiveThorns = "AHFinderDirect"
+AHFinderDirect::N_horizons = 2
+
+AHFinderDirect::origin_y[1]                                       = -1.168433
+AHFinderDirect::initial_guess_method[1]                           = "read from named file"
+AHFinderDirect::initial_guess__read_from_named_file__file_name[1] = "/data3/diener/pi0-384_dc_cou_20_3d2/h.t101.ah1.gp"
+AHFinderDirect::which_surface_to_store_info[1]                    =  0
+
+AHFinderDirect::origin_y[2]                                       = +1.168433
+AHFinderDirect::initial_guess_method[2]                           = "read from named file"
+AHFinderDirect::initial_guess__read_from_named_file__file_name[2] = "/data3/diener/pi0-384_dc_cou_20_3d2/h.t101.ah2.gp"
+
+#AHFinderDirect::initial_guess_method[3]                           = "read from named file"
+#AHFinderDirect::initial_guess__read_from_named_file__file_name[3]= "/data3/diener/pi0-384_dc_cou_20_3d2/h.t101.ah3.gp"
+
+# centroid 0.003181 1.168433 0.000000
+# radius 0.4446773218
+#
+# centroid -0.003181 -1.168433 0.000000
+# radius  0.4446773520
+#
+# centroid 0.000000 -0.000000 0.000000
+# radius 1.145238733 1.780038605 1.020373858
+
+ActiveThorns = "IsolatedHorizon"
+IsolatedHorizon::verbose       = yes
+IsolatedHorizon::veryverbose   = yes
+IsolatedHorizon::surface_index = 0
+
+ActiveThorns = "IOBasic"
+IOBasic::outInfo_every      = 1
+IOBasic::outInfo_reductions = "norm2"
+IOBasic::outInfo_vars       = "
+        IsolatedHorizon::ih_spin
+        IsolatedHorizon::ih_radius
+        IsolatedHorizon::ih_mass
+"
+
+IOBasic::outScalar_every      = 1
+IOBasic::outScalar_reductions = "minimum maximum mean norm1 norm2 norm_inf"
+IOBasic::outScalar_style      = gnuplot
+IOBasic::outScalar_vars       = "
+        ADMBase::gxx
+        ADMBase::kxx
+        ADMBase::lapse
+        SphericalSurface::sf_valid
+        SphericalSurface::sf_info
+        SphericalSurface::sf_radius
+        SphericalSurface::sf_origin
+        SphericalSurface::sf_coordinate_descriptors
+        IsolatedHorizon::ih_state
+        IsolatedHorizon::ih_time
+        IsolatedHorizon::ih_grid_int
+        IsolatedHorizon::ih_grid_real
+        IsolatedHorizon::ih_shape
+        IsolatedHorizon::ih_coordinates
+        IsolatedHorizon::ih_tetrad_l
+        IsolatedHorizon::ih_tetrad_n
+        IsolatedHorizon::ih_tetrad_m
+        IsolatedHorizon::ih_newman_penrose
+        IsolatedHorizon::ih_twometric
+        IsolatedHorizon::ih_killing_vector
+        IsolatedHorizon::ih_scalars
+        IsolatedHorizon::ih_variables
+"
+
+ActiveThorns = "IOASCII"
+IOASCII::out1D_every = 1
+IOASCII::out1D_style = "gnuplot f(x)"
+IOASCII::out1D_vars  = "
+        ADMBase::gxx
+        ADMBase::kxx
+        ADMBase::lapse
+        SphericalSurface::sf_radius
+        IsolatedHorizon::ih_shape
+        IsolatedHorizon::ih_coordinates
+        IsolatedHorizon::ih_tetrad_l
+        IsolatedHorizon::ih_tetrad_n
+        IsolatedHorizon::ih_tetrad_m
+        IsolatedHorizon::ih_newman_penrose
+        IsolatedHorizon::ih_twometric
+        IsolatedHorizon::ih_killing_vector
+        IsolatedHorizon::ih_variables
+"
+
+IOASCII::out2D_every = 1
+IOASCII::out2D_style = "gnuplot f(x,y)"
+IOASCII::out2D_vars  = "
+        SphericalSurface::sf_radius
+        IsolatedHorizon::ih_shape
+        IsolatedHorizon::ih_coordinates
+        IsolatedHorizon::ih_tetrad_l
+        IsolatedHorizon::ih_tetrad_n
+        IsolatedHorizon::ih_tetrad_m
+        IsolatedHorizon::ih_newman_penrose
+        IsolatedHorizon::ih_twometric
+        IsolatedHorizon::ih_killing_vector
+        IsolatedHorizon::ih_variables
+"
+
+IOHDF5::out_every = 1
+IOHDF5::out_vars  = "
+	StaticConformal::confac
+        ADMBase::metric
+        ADMBase::curv
+        ADMBase::lapse
+        ADMBase::shift
+"
diff --git a/src/CalcK.F90 b/src/CalcK.F90
index f77d500..2dacc50 100644
--- a/src/CalcK.F90
+++ b/src/CalcK.F90
@@ -22,6 +22,7 @@ subroutine CalcK (CCTK_ARGUMENTS)
   
   CCTK_REAL :: dt, dx(3)
   
+  CCTK_REAL :: ps, dps(3)
   CCTK_REAL :: gama(3,3), gama_dot(3,3), dgama(3,3,3)
   CCTK_REAL :: alfa
   CCTK_REAL :: beta(3), dbeta(3,3)
@@ -53,6 +54,8 @@ subroutine CalcK (CCTK_ARGUMENTS)
      do j = imin(2), imax(2)
         do i = imin(1), imax(1)
            
+           ps = psi(i,j,k)
+           
            gama(1,1) = gxx(i,j,k)
            gama(1,2) = gxy(i,j,k)
            gama(1,3) = gxz(i,j,k)
@@ -82,6 +85,10 @@ subroutine CalcK (CCTK_ARGUMENTS)
            gama_dot(3,1) = gama_dot(1,3)
            gama_dot(3,2) = gama_dot(2,3)
            
+           dps(1) = (psi(i+1,j,k) - psi(i-1,j,k)) / (2*dx(1))
+           dps(2) = (psi(i,j+1,k) - psi(i,j-1,k)) / (2*dx(2))
+           dps(3) = (psi(i,j,k+1) - psi(i,j,k-1)) / (2*dx(3))
+           
            dgama(1,1,1) = (gxx(i+1,j,k) - gxx(i-1,j,k)) / (2*dx(1))
            dgama(1,2,1) = (gxy(i+1,j,k) - gxy(i-1,j,k)) / (2*dx(1))
            dgama(1,3,1) = (gxz(i+1,j,k) - gxz(i-1,j,k)) / (2*dx(1))
@@ -127,11 +134,11 @@ subroutine CalcK (CCTK_ARGUMENTS)
            
            do a=1,3
               do b=1,3
-                 kk(a,b) = - gama_dot(a,b)
+                 kk(a,b) = - ps**4 * gama_dot(a,b)
                  do c=1,3
-                    kk(a,b) = kk(a,b) + gama(k,j) * dbeta(k,i) &
-                         &            + gama(i,k) * dbeta(k,j) &
-                         &            + beta(k) * dgama(i,j,k)
+                    kk(a,b) = kk(a,b) + ps**4 * gama(k,j) * dbeta(k,i) &
+                         &            + ps**4 * gama(i,k) * dbeta(k,j) &
+                         &            + beta(k) * (ps**4 * dgama(i,j,k) + 4*ps**3 * dps(k) * gama(i,j))
                  end do
                  kk(a,b) = kk(a,b) / (2*alfa)
               end do