From 3c4efa5d0aee2452f15b77db71f0fc5949f42937 Mon Sep 17 00:00:00 2001 From: schnetter Date: Thu, 5 Feb 2004 13:02:19 +0000 Subject: 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 --- interface.ccl | 26 ++++++++++++++++++++++++++ par/calck.par | 2 +- par/copy_to_next.par | 2 +- par/copy_to_prev.par | 2 +- par/ih.par | 2 +- param.ccl | 6 ++++++ src/CalcK.F90 | 53 +++++++++++++++++++++++++++++++++++++++++----------- 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 -- cgit v1.2.3