diff options
Diffstat (limited to 'src/qlm_3determinant.F90')
-rw-r--r-- | src/qlm_3determinant.F90 | 112 |
1 files changed, 112 insertions, 0 deletions
diff --git a/src/qlm_3determinant.F90 b/src/qlm_3determinant.F90 new file mode 100644 index 0000000..d263009 --- /dev/null +++ b/src/qlm_3determinant.F90 @@ -0,0 +1,112 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +SUBROUTINE qlm_calc_3determinant (CCTK_ARGUMENTS, hn) + USE adm_metric + USE cctk + USE qlm_boundary + USE qlm_derivs + USE qlm_variables + USE tensor + + IMPLICIT NONE + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + INTEGER :: hn + + CCTK_REAL :: gg(3,3) + CCTK_REAL :: alfa + CCTK_REAL :: beta(3) + CCTK_REAL :: g4(0:3,0:3) + CCTK_COMPLEX :: mm(0:3) + CCTK_REAL :: rr(3), rr_p(3), rr_p_p(3) + CCTK_REAL :: Ttilde(0:3) + CCTK_REAL :: a + CCTK_COMPLEX :: d + + CCTK_REAL :: t0, t1, t2 + logical :: ce0, ce1, ce2 + CCTK_REAL :: delta_space(2) + + INTEGER :: i, j + CCTK_REAL :: theta, phi + + IF (veryverbose/=0) THEN + CALL CCTK_INFO ("Computing 3-Volume Element of Horizon") + END IF + + t0 = qlm_time(hn) + t1 = qlm_time_p(hn) + t2 = qlm_time_p_p(hn) + + ce0 = qlm_have_valid_data(hn) == 0 + ce1 = qlm_have_valid_data_p(hn) == 0 + ce2 = qlm_have_valid_data_p_p(hn) == 0 + + delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /) + + ! 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) + + alfa = qlm_alpha(i,j) + + beta(1) = qlm_betax(i,j) + beta(2) = qlm_betay(i,j) + beta(3) = qlm_betaz(i,j) + + 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) + + ! Build the four-metric + CALL calc_4metric (gg,alfa,beta, g4) + + ! Find a 3rd vector of triad + rr(1) = qlm_x(i,j,hn) + rr(2) = qlm_y(i,j,hn) + rr(3) = qlm_z(i,j,hn) + + rr_p(1) = qlm_x_p(i,j,hn) + rr_p(2) = qlm_y_p(i,j,hn) + rr_p(3) = qlm_z_p(i,j,hn) + + rr_p_p(1) = qlm_x_p_p(i,j,hn) + rr_p_p(2) = qlm_y_p_p(i,j,hn) + rr_p_p(3) = qlm_z_p_p(i,j,hn) + + Ttilde(0) = 1 + Ttilde(1:3) = timederiv (rr, rr_p, rr_p_p, t0,t1,t2, ce0,ce1,ce2) + + ! Compute some scalar products + a=SUM(MATMUL(g4, Ttilde)*Ttilde) + d=SUM(MATMUL(g4, mm)*Ttilde) + + ! This is the determinant + qlm_3det(i,j,hn)=a-2*CONJG(d)*d + + END DO + END DO + + CALL set_boundary (CCTK_PASS_FTOF, hn, qlm_3det(:,:,hn), +1) + +END SUBROUTINE qlm_calc_3determinant |