Use function aliasing to convert between physical and conformal

metric.

Fix severe indexing bug in calculating K_ij.

Use new ADMBase convention for not initialising the fields.


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

7 files changed, 78 insertions, 15 deletions

diff --git a/interface.ccl b/interface.ccl
index d4dc0ca..cd74515 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -30,3 +30,29 @@ CCTK_INT FUNCTION Boundary_SelectGroupForBC (	\
 	CCTK_STRING IN group_name,		\
 	CCTK_STRING IN bc_name)
 USES FUNCTION Boundary_SelectGroupForBC
+
+
+
+void FUNCTION ConfToPhysInPlace (CCTK_INT IN nx,            \
+                                 CCTK_INT IN ny,            \
+                                 CCTK_INT IN nz,            \
+                                 CCTK_REAL ARRAY IN psi,    \
+                                 CCTK_REAL ARRAY INOUT gxx, \
+                                 CCTK_REAL ARRAY INOUT gxy, \
+                                 CCTK_REAL ARRAY INOUT gxz, \
+                                 CCTK_REAL ARRAY INOUT gyy, \
+                                 CCTK_REAL ARRAY INOUT gyz, \
+                                 CCTK_REAL ARRAY INOUT gzz)
+USES FUNCTION ConfToPhysInPlace
+
+void FUNCTION PhysToConfInPlace (CCTK_INT IN nx,            \
+                                 CCTK_INT IN ny,            \
+                                 CCTK_INT IN nz,            \
+                                 CCTK_REAL ARRAY IN psi,    \
+                                 CCTK_REAL ARRAY INOUT gxx, \
+                                 CCTK_REAL ARRAY INOUT gxy, \
+                                 CCTK_REAL ARRAY INOUT gxz, \
+                                 CCTK_REAL ARRAY INOUT gyy, \
+                                 CCTK_REAL ARRAY INOUT gyz, \
+                                 CCTK_REAL ARRAY INOUT gzz)
+USES FUNCTION PhysToConfInPlace
diff --git a/par/calck.par b/par/calck.par
index 1fc1926..ec4b83d 100644
--- a/par/calck.par
+++ b/par/calck.par
@@ -46,7 +46,7 @@ ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConform
 ADMBase::metric_type = "static conformal"
 StaticConformal::conformal_storage = "factor"
 
-ADMBase::initial_data = none
+ADMBase::initial_data = unset
 ADMBase::initial_shift = zero
 
 ActiveThorns = "CalcK"
diff --git a/par/copy_to_next.par b/par/copy_to_next.par
index 7cb1dea..0aede90 100644
--- a/par/copy_to_next.par
+++ b/par/copy_to_next.par
@@ -36,7 +36,7 @@ ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConform
 ADMBase::metric_type = "static conformal"
 StaticConformal::conformal_storage = "factor"
 
-ADMBase::initial_data = none
+ADMBase::initial_data = unset
 ADMBase::initial_shift = zero
 
 ActiveThorns = "CalcK"
diff --git a/par/copy_to_prev.par b/par/copy_to_prev.par
index f677fd9..2d76df7 100644
--- a/par/copy_to_prev.par
+++ b/par/copy_to_prev.par
@@ -36,7 +36,7 @@ ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConform
 ADMBase::metric_type = "static conformal"
 StaticConformal::conformal_storage = "factor"
 
-ADMBase::initial_data = none
+ADMBase::initial_data = unset
 ADMBase::initial_shift = zero
 
 ActiveThorns = "CalcK"
diff --git a/par/ih.par b/par/ih.par
index baaf839..ca171be 100644
--- a/par/ih.par
+++ b/par/ih.par
@@ -55,7 +55,7 @@ ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConform
 ADMBase::metric_type = "static conformal"
 StaticConformal::conformal_storage = "factor"
 
-ADMBase::initial_data = none
+ADMBase::initial_data = unset
 ADMBase::initial_shift = zero
 
 ActiveThorns = "AHFinderDirect"
diff --git a/param.ccl b/param.ccl
index 71918e4..b46399d 100644
--- a/param.ccl
+++ b/param.ccl
@@ -29,3 +29,9 @@ STRING extcurv_boundary_options "Options for the boundry condition for the extri
 {
   ".*" :: "must have the options table syntax"
 } "decay_power=2.0"
+
+
+
+SHARES: ADMBase
+
+USES KEYWORD metric_type
diff --git a/src/CalcK.F90 b/src/CalcK.F90
index 2dacc50..7a64d5a 100644
--- a/src/CalcK.F90
+++ b/src/CalcK.F90
@@ -22,7 +22,6 @@ 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)
@@ -50,12 +49,29 @@ subroutine CalcK (CCTK_ARGUMENTS)
   where (cctk_bbox(1::2)/=0) imin(:) = 1+bndwidth
   where (cctk_bbox(2::2)/=0) imax(:) = cctk_lsh(:)-bndwidth
   
+  ! Convert to physical metric, if necessary
+  if (CCTK_EQUALS(metric_type, "static conformal")) then
+     call ConfToPhysInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx, gxy, gxz, gyy, gyz, gzz)
+     call ConfToPhysInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx_prev, gxy_prev, gxz_prev, gyy_prev, gyz_prev, gzz_prev)
+#if 0
+     call ConfToPhysInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx_prev2, gxy_prev2, gxz_prev2, gyy_prev2, gyz_prev2, gzz_prev2)
+#else
+     call ConfToPhysInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx_next, gxy_next, gxz_next, gyy_next, gyz_next, gzz_next)
+#endif
+  end if
+  
   do k = imin(3), imax(3)
      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)
@@ -85,10 +101,6 @@ 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))
@@ -134,11 +146,11 @@ subroutine CalcK (CCTK_ARGUMENTS)
            
            do a=1,3
               do b=1,3
-                 kk(a,b) = - ps**4 * gama_dot(a,b)
+                 kk(a,b) = - gama_dot(a,b)
                  do c=1,3
-                    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))
+                    kk(a,b) = kk(a,b) + gama(c,b) * dbeta(c,a) &
+                         &            + gama(a,c) * dbeta(c,b) &
+                         &            + beta(c) * dgama(a,b,c)
                  end do
                  kk(a,b) = kk(a,b) / (2*alfa)
               end do
@@ -161,4 +173,23 @@ subroutine CalcK (CCTK_ARGUMENTS)
        int(options,ik), "ADMBase::curv", str_boundary)
   if (ierr /= 0) call CCTK_WARN (0, "internal error")
   
+  ! Convert back to conformal metric, if necessary
+  if (CCTK_EQUALS(metric_type, "static conformal")) then
+     call PhysToConfInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx, gxy, gxz, gyy, gyz, gzz)
+     call PhysToConfInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx_prev, gxy_prev, gxz_prev, gyy_prev, gyz_prev, gzz_prev)
+#if 0
+     call PhysToConfInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx_prev2, gxy_prev2, gxz_prev2, gyy_prev2, gyz_prev2, gzz_prev2)
+#else
+     call PhysToConfInPlace &
+          (int(cctk_lsh(1),ik), int(cctk_lsh(2),ik), int(cctk_lsh(3),ik), &
+          psi, gxx_next, gxy_next, gxz_next, gyy_next, gyz_next, gzz_next)
+#endif
+  end if
+  
 end subroutine CalcK