diff options
Diffstat (limited to 'src/qlm_weyl_scalars.F90')
| -rw-r--r-- | src/qlm_weyl_scalars.F90 | 252 |
1 files changed, 252 insertions, 0 deletions
diff --git a/src/qlm_weyl_scalars.F90 b/src/qlm_weyl_scalars.F90 new file mode 100644 index 0000000..d79751c --- /dev/null +++ b/src/qlm_weyl_scalars.F90 @@ -0,0 +1,252 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) + use adm_metric_simple + use cctk + use constants + use qlm_boundary + use qlm_derivs + use qlm_variables + use ricci + use ricci4 + use tensor + use tensor4 + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: two=2, four=4 + CCTK_REAL :: gg(3,3), dgg(3,3,3), ddgg(3,3,3,3), gg_dot(3,3), gg_dot2(3,3), dgg_dot(3,3,3) + CCTK_REAL :: g4(0:3,0:3), dg4(0:3,0:3,0:3), ddg4(0:3,0:3,0:3,0:3) + CCTK_REAL :: gu4(0:3,0:3), dgu4(0:3,0:3,0:3) + CCTK_REAL :: gamma4(0:3,0:3,0:3), dgamma4(0:3,0:3,0:3,0:3) + CCTK_REAL :: ri4(0:3,0:3), rsc4 + CCTK_REAL :: rm4(0:3,0:3,0:3,0:3), we4(0:3,0:3,0:3,0:3) + CCTK_REAL :: ll(0:3), nn(0:3) + CCTK_COMPLEX :: mm(0:3) + + integer :: i, j + integer :: a, b, c, d + CCTK_REAL :: theta, phi + + if (veryverbose/=0) then + call CCTK_INFO ("Calculating Weyl scalars") + end if + + ! Calculate the coordinates + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + theta = qlm_origin_theta(hn) + (i-1)*qlm_delta_theta(hn) + phi = qlm_origin_phi(hn) + (j-1)*qlm_delta_phi(hn) + + ! Get the stuff from the arrays + 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) + + ll(0) = qlm_l0(i,j,hn) + ll(1) = qlm_l1(i,j,hn) + ll(2) = qlm_l2(i,j,hn) + ll(3) = qlm_l3(i,j,hn) + + nn(0) = qlm_n0(i,j,hn) + nn(1) = qlm_n1(i,j,hn) + nn(2) = qlm_n2(i,j,hn) + nn(3) = qlm_n3(i,j,hn) + + mm(0) = qlm_m0(i,j,hn) + mm(1) = qlm_m1(i,j,hn) + mm(2) = qlm_m2(i,j,hn) + mm(3) = qlm_m3(i,j,hn) + + + + ! Calculate 4-metric + call calc_4metricderivs2_simple (gg, dgg, & + ddgg, gg_dot, gg_dot2, dgg_dot, g4,dg4,ddg4) + call calc_4inv (g4, gu4) + call calc_4invderiv (gu4, dg4, dgu4) + call calc_4connections (gu4,dg4, gamma4) + call calc_4connectionderivs (gu4, dg4, dgu4, ddg4, dgamma4) + call calc_4ricci (gamma4, dgamma4, ri4) + call calc_4riemann (g4, gamma4, dgamma4, rm4) + call calc_4trace (ri4, gu4, rsc4) + call calc_4weyl (g4, rm4, ri4, rsc4, we4) + + ! debugging +! qlm_rsc4(i,j,hn) = rsc4 + + qlm_psi0(i,j,hn) = 0 ! transverse radiation along n + qlm_psi1(i,j,hn) = 0 ! longitudinal radiation along n + qlm_psi2(i,j,hn) = 0 ! Coulomb field and spin + qlm_psi3(i,j,hn) = 0 ! longitudinal radiation along l + qlm_psi4(i,j,hn) = 0 ! transverse radiation along l + + do a=0,3 + do b=0,3 + do c=0,3 + do d=0,3 + qlm_psi0(i,j,hn) = qlm_psi0(i,j,hn) + we4(a,b,c,d) * ll(a) * mm(b) * ll(c) * mm(d) + qlm_psi1(i,j,hn) = qlm_psi1(i,j,hn) + we4(a,b,c,d) * ll(a) * mm(b) * ll(c) * nn(d) + qlm_psi2(i,j,hn) = qlm_psi2(i,j,hn) + we4(a,b,c,d) * ll(a) * mm(b) * conjg(mm(c)) * nn(d) + qlm_psi3(i,j,hn) = qlm_psi3(i,j,hn) + we4(a,b,c,d) * ll(a) * nn(b) * conjg(mm(c)) * nn(d) + qlm_psi4(i,j,hn) = qlm_psi4(i,j,hn) + we4(a,b,c,d) * conjg(mm(a)) * nn(b) * conjg(mm(c)) * nn(d) + end do + end do + end do + end do + + ! gr-qc/0104063, (3.3) + qlm_i(i,j,hn) = + 3 * qlm_psi2(i,j,hn)**2 & + & - 4 * qlm_psi1(i,j,hn) * qlm_psi3(i,j,hn) & + & + qlm_psi0(i,j,hn) * qlm_psi4(i,j,hn) + qlm_j(i,j,hn) = - qlm_psi2(i,j,hn)**3 & + & + qlm_psi0(i,j,hn) * qlm_psi2(i,j,hn) * qlm_psi4(i,j,hn) & + & + 2 * qlm_psi1(i,j,hn) * qlm_psi2(i,j,hn) * qlm_psi3(i,j,hn) & + & - qlm_psi1(i,j,hn)**2 * qlm_psi4(i,j,hn) & + & - qlm_psi0(i,j,hn) * qlm_psi3(i,j,hn)**2 + + ! gr-qc/0104063, (3.1) + qlm_s(i,j,hn) = 27 * qlm_j(i,j,hn)**2 / qlm_i(i,j,hn)**3 + qlm_sdiff(i,j,hn) = (27 * qlm_j(i,j,hn)**2 - qlm_i(i,j,hn)**3) / sqrt((abs2(qlm_psi0(i,j,hn)) + abs2(qlm_psi1(i,j,hn)) + abs2(qlm_psi2(i,j,hn)) + abs2(qlm_psi3(i,j,hn)) + abs2(qlm_psi4(i,j,hn))) / 5)**3 + + qlm_phi00(i,j,hn) = 0 + qlm_phi11(i,j,hn) = 0 + qlm_phi01(i,j,hn) = 0 + qlm_phi12(i,j,hn) = 0 + qlm_phi10(i,j,hn) = 0 + qlm_phi21(i,j,hn) = 0 + qlm_phi02(i,j,hn) = 0 + qlm_phi22(i,j,hn) = 0 + qlm_phi20(i,j,hn) = 0 + + do a=0,3 + do b=0,3 + + qlm_phi00(i,j,hn) = qlm_phi00(i,j,hn) - 1/two * ri4(a,b) * ll(a) * ll(b) + qlm_phi11(i,j,hn) = qlm_phi11(i,j,hn) - 1/two * ri4(a,b) * (ll(a) * nn(b) + mm(a) * conjg(mm(b))) / 2 + qlm_phi01(i,j,hn) = qlm_phi01(i,j,hn) - 1/two * ri4(a,b) * ll(a) * mm(b) + qlm_phi12(i,j,hn) = qlm_phi12(i,j,hn) - 1/two * ri4(a,b) * nn(a) * mm(b) + qlm_phi10(i,j,hn) = qlm_phi10(i,j,hn) - 1/two * ri4(a,b) * ll(a) * conjg(mm(b)) + qlm_phi21(i,j,hn) = qlm_phi21(i,j,hn) - 1/two * ri4(a,b) * nn(a) * conjg(mm(b)) + qlm_phi02(i,j,hn) = qlm_phi02(i,j,hn) - 1/two * ri4(a,b) * mm(a) * mm(b) + qlm_phi22(i,j,hn) = qlm_phi22(i,j,hn) - 1/two * ri4(a,b) * nn(a) * nn(b) + qlm_phi20(i,j,hn) = qlm_phi20(i,j,hn) - 1/two * ri4(a,b) * conjg(mm(a)) * conjg(mm(b)) + + end do + end do + + qlm_lambda(i,j,hn) = rsc4 / 24 + + qlm_lie_n_theta_l(i,j,hn) = & + & + 2 * real (qlm_npsigma(i,j,hn) * qlm_nplambda(i,j,hn)) & + & + 2 * real (qlm_psi2(i,j,hn)) & + & + 4 * qlm_lambda(i,j,hn) + + end do + end do + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_psi0(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_psi1(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_psi2(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_psi3(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_psi4(:,:,hn), +1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_i(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_j(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_s(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_sdiff(:,:,hn), +1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi00(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi11(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi01(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi12(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi10(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi21(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi02(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi22(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_phi20(:,:,hn), +1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_lambda(:,:,hn), +1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_lie_n_theta_l(:,:,hn), +1) + +end subroutine qlm_calc_weyl_scalars |