diff options
Diffstat (limited to 'src/qlm_twometric.F90')
| -rw-r--r-- | src/qlm_twometric.F90 | 235 |
1 files changed, 235 insertions, 0 deletions
diff --git a/src/qlm_twometric.F90 b/src/qlm_twometric.F90 new file mode 100644 index 0000000..e450ebb --- /dev/null +++ b/src/qlm_twometric.F90 @@ -0,0 +1,235 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_calc_twometric (CCTK_ARGUMENTS, hn) + use adm_metric + use cctk + use qlm_boundary + use qlm_derivs + use qlm_variables + use ricci + use ricci2 + use tensor + use tensor2 + use tensor4 + + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: two = 2, half = 1d0/2d0 + CCTK_REAL :: gg(3,3), dgg(3,3,3), ddgg(3,3,3,3) + CCTK_REAL :: gamma(3,3,3), dgamma(3,3,3,3), ri(3,3), rsc + CCTK_REAL :: ee(3,2), dee(3,2,2), ddee(3,2,2,2) + CCTK_REAL :: qq(2,2), dqq(2,2,2), ddqq(2,2,2,2) + CCTK_REAL :: dtq, qu(2,2), dqu(2,2,2) + + CCTK_REAL :: delta_space(2) + + integer :: i, j + integer :: a, b, c, d, e, f, g, h + + if (veryverbose/=0) then + call CCTK_INFO ("Calculating two-metric") + end if + + delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /) + + ! Calculate the two-metric + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + + gg(1,1) = qlm_gxx(i,j) + gg(1,2) = qlm_gxy(i,j) + gg(1,3) = qlm_gxz(i,j) + gg(2,2) = qlm_gyy(i,j) + gg(2,3) = qlm_gyz(i,j) + gg(3,3) = qlm_gzz(i,j) + gg(2,1) = gg(1,2) + gg(3,1) = gg(1,3) + gg(3,2) = gg(2,3) + + dgg(1,1,1) = qlm_dgxxx(i,j) + dgg(1,2,1) = qlm_dgxyx(i,j) + dgg(1,3,1) = qlm_dgxzx(i,j) + dgg(2,2,1) = qlm_dgyyx(i,j) + dgg(2,3,1) = qlm_dgyzx(i,j) + dgg(3,3,1) = qlm_dgzzx(i,j) + dgg(1,1,2) = qlm_dgxxy(i,j) + dgg(1,2,2) = qlm_dgxyy(i,j) + dgg(1,3,2) = qlm_dgxzy(i,j) + dgg(2,2,2) = qlm_dgyyy(i,j) + dgg(2,3,2) = qlm_dgyzy(i,j) + dgg(3,3,2) = qlm_dgzzy(i,j) + dgg(1,1,3) = qlm_dgxxz(i,j) + dgg(1,2,3) = qlm_dgxyz(i,j) + dgg(1,3,3) = qlm_dgxzz(i,j) + dgg(2,2,3) = qlm_dgyyz(i,j) + dgg(2,3,3) = qlm_dgyzz(i,j) + dgg(3,3,3) = qlm_dgzzz(i,j) + dgg(2,1,:) = dgg(1,2,:) + dgg(3,1,:) = dgg(1,3,:) + dgg(3,2,:) = dgg(2,3,:) + + ddgg(1,1,1,1) = qlm_ddgxxxx(i,j) + ddgg(1,2,1,1) = qlm_ddgxyxx(i,j) + ddgg(1,3,1,1) = qlm_ddgxzxx(i,j) + ddgg(2,2,1,1) = qlm_ddgyyxx(i,j) + ddgg(2,3,1,1) = qlm_ddgyzxx(i,j) + ddgg(3,3,1,1) = qlm_ddgzzxx(i,j) + ddgg(1,1,1,2) = qlm_ddgxxxy(i,j) + ddgg(1,2,1,2) = qlm_ddgxyxy(i,j) + ddgg(1,3,1,2) = qlm_ddgxzxy(i,j) + ddgg(2,2,1,2) = qlm_ddgyyxy(i,j) + ddgg(2,3,1,2) = qlm_ddgyzxy(i,j) + ddgg(3,3,1,2) = qlm_ddgzzxy(i,j) + ddgg(1,1,1,3) = qlm_ddgxxxz(i,j) + ddgg(1,2,1,3) = qlm_ddgxyxz(i,j) + ddgg(1,3,1,3) = qlm_ddgxzxz(i,j) + ddgg(2,2,1,3) = qlm_ddgyyxz(i,j) + ddgg(2,3,1,3) = qlm_ddgyzxz(i,j) + ddgg(3,3,1,3) = qlm_ddgzzxz(i,j) + ddgg(1,1,2,2) = qlm_ddgxxyy(i,j) + ddgg(1,2,2,2) = qlm_ddgxyyy(i,j) + ddgg(1,3,2,2) = qlm_ddgxzyy(i,j) + ddgg(2,2,2,2) = qlm_ddgyyyy(i,j) + ddgg(2,3,2,2) = qlm_ddgyzyy(i,j) + ddgg(3,3,2,2) = qlm_ddgzzyy(i,j) + ddgg(1,1,2,3) = qlm_ddgxxyz(i,j) + ddgg(1,2,2,3) = qlm_ddgxyyz(i,j) + ddgg(1,3,2,3) = qlm_ddgxzyz(i,j) + ddgg(2,2,2,3) = qlm_ddgyyyz(i,j) + ddgg(2,3,2,3) = qlm_ddgyzyz(i,j) + ddgg(3,3,2,3) = qlm_ddgzzyz(i,j) + ddgg(1,1,3,3) = qlm_ddgxxzz(i,j) + ddgg(1,2,3,3) = qlm_ddgxyzz(i,j) + ddgg(1,3,3,3) = qlm_ddgxzzz(i,j) + ddgg(2,2,3,3) = qlm_ddgyyzz(i,j) + ddgg(2,3,3,3) = qlm_ddgyzzz(i,j) + ddgg(3,3,3,3) = qlm_ddgzzzz(i,j) + ddgg(2,1,:,:) = ddgg(1,2,:,:) + ddgg(3,1,:,:) = ddgg(1,3,:,:) + ddgg(3,2,:,:) = ddgg(2,3,:,:) + ddgg(:,:,2,1) = ddgg(:,:,1,2) + ddgg(:,:,3,1) = ddgg(:,:,1,3) + ddgg(:,:,3,2) = ddgg(:,:,2,3) + + ee(1,1:2) = deriv (qlm_x(:,:,hn), i, j, delta_space) + ee(2,1:2) = deriv (qlm_y(:,:,hn), i, j, delta_space) + ee(3,1:2) = deriv (qlm_z(:,:,hn), i, j, delta_space) + + dee(1,1:2,1:2) = deriv2 (qlm_x(:,:,hn), i, j, delta_space) + dee(2,1:2,1:2) = deriv2 (qlm_y(:,:,hn), i, j, delta_space) + dee(3,1:2,1:2) = deriv2 (qlm_z(:,:,hn), i, j, delta_space) + + ddee(1,1:2,1:2,1:2) = deriv3 (qlm_x(:,:,hn), i, j, delta_space) + ddee(2,1:2,1:2,1:2) = deriv3 (qlm_y(:,:,hn), i, j, delta_space) + ddee(3,1:2,1:2,1:2) = deriv3 (qlm_z(:,:,hn), i, j, delta_space) + + do a=1,2 + do b=1,2 + qq(a,b) = 0 + do c=1,3 + do d=1,3 + qq(a,b) = qq(a,b) + gg(c,d) * ee(c,a) * ee(d,b) + end do + end do + end do + end do + + do a=1,2 + do b=1,2 + do c=1,2 + dqq(a,b,c) = 0 + do d=1,3 + do e=1,3 + do f=1,3 + dqq(a,b,c) = dqq(a,b,c) + dgg(d,e,f) * ee(d,a) * ee(e,b) * ee(f,c) + end do + dqq(a,b,c) = dqq(a,b,c) + gg(d,e) * dee(d,a,c) * ee(e,b) + dqq(a,b,c) = dqq(a,b,c) + gg(d,e) * ee(d,a) * dee(e,b,c) + end do + end do + end do + end do + end do + + do a=1,2 + do b=1,2 + do c=1,2 + do d=1,2 + ddqq(a,b,c,d) = 0 + do e=1,3 + do f=1,3 + do g=1,3 + do h=1,3 + ddqq(a,b,c,d) = ddqq(a,b,c,d) + ddgg(e,f,g,h) * ee(e,a) * ee(f,b) * ee(g,c) * ee(h,d) + end do + ddqq(a,b,c,d) = ddqq(a,b,c,d) + dgg(e,f,g) * dee(e,a,d) * ee(f,b) * ee(g,c) + ddqq(a,b,c,d) = ddqq(a,b,c,d) + dgg(e,f,g) * ee(e,a) * dee(f,b,d) * ee(g,c) + ddqq(a,b,c,d) = ddqq(a,b,c,d) + dgg(e,f,g) * ee(e,a) * ee(f,b) * dee(g,c,d) + ddqq(a,b,c,d) = ddqq(a,b,c,d) + dgg(e,f,g) * dee(e,a,c) * ee(f,b) * ee(g,d) + ddqq(a,b,c,d) = ddqq(a,b,c,d) + dgg(e,f,g) * ee(e,a) * dee(f,b,c) * ee(g,d) + end do + ddqq(a,b,c,d) = ddqq(a,b,c,d) + gg(e,f) * ddee(e,a,c,d) * ee(f,b) + ddqq(a,b,c,d) = ddqq(a,b,c,d) + gg(e,f) * dee(e,a,c) * dee(f,b,d) + ddqq(a,b,c,d) = ddqq(a,b,c,d) + gg(e,f) * dee(e,a,d) * dee(f,b,c) + ddqq(a,b,c,d) = ddqq(a,b,c,d) + gg(e,f) * ee(e,a) * ddee(f,b,c,d) + end do + end do + end do + end do + end do + end do + + call calc_2det (qq, dtq) + call calc_2inv (qq, dtq, qu) + call calc_2invderiv (qu, dqq, dqu) + + call calc_2connections (qu, dqq, gamma) + call calc_2connectionderivs (qu, dqq, dqu, ddqq, dgamma) + call calc_2ricci (gamma, dgamma, ri) + + call calc_2trace (qu, ri, rsc) + + ! Could also calculate this as: + ! q^ab = m^a mbar^b + mbar^a m^b + qlm_qtt(i,j,hn) = qq(1,1) + qlm_qtp(i,j,hn) = qq(1,2) + qlm_qpp(i,j,hn) = qq(2,2) + + ! Could also calculate this from NP coefficients alpha and + ! beta (and epsilon and gamma?) + qlm_dqttt(i,j,hn) = dqq(1,1,1) + qlm_dqtpt(i,j,hn) = dqq(1,2,1) + qlm_dqppt(i,j,hn) = dqq(2,2,1) + qlm_dqttp(i,j,hn) = dqq(1,1,2) + qlm_dqtpp(i,j,hn) = dqq(1,2,2) + qlm_dqppp(i,j,hn) = dqq(2,2,2) + + ! Could also calculate this from 3Rsc and extrinsic curvature + qlm_rsc(i,j,hn) = rsc + + end do + end do + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_qtt(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_qtp(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_qpp(:,:,hn), +1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_dqttt(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_dqtpt(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_dqppt(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_dqttp(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_dqtpp(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_dqppp(:,:,hn), -1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_rsc(:,:,hn), +1) + +end subroutine qlm_calc_twometric |