diff options
Diffstat (limited to 'src/qlm_newman_penrose.F90')
| -rw-r--r-- | src/qlm_newman_penrose.F90 | 156 |
1 files changed, 156 insertions, 0 deletions
diff --git a/src/qlm_newman_penrose.F90 b/src/qlm_newman_penrose.F90 new file mode 100644 index 0000000..b6cc742 --- /dev/null +++ b/src/qlm_newman_penrose.F90 @@ -0,0 +1,156 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_calc_newman_penrose (CCTK_ARGUMENTS, hn) + use adm_metric + use cctk + use classify + use constants + use qlm_boundary + use qlm_derivs + use qlm_variables + use ricci4 + use tensor + use tensor4 + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: zero=0, one=1, two=2, half=1/two + CCTK_REAL :: ll(0:3), nn(0:3) + CCTK_COMPLEX :: mm(0:3) + CCTK_REAL :: nabla_ll(0:3,0:3), nabla_nn(0:3,0:3) + CCTK_COMPLEX :: nabla_mm(0:3,0:3) + + CCTK_REAL :: t0, t1, t2 + logical :: ce0, ce1, ce2 + CCTK_REAL :: delta_space(2) + + integer :: i, j + integer :: a, b + CCTK_REAL :: theta, phi + + if (veryverbose/=0) then + call CCTK_INFO ("Calculating Newman-Penrose quantities") + 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) /) + + if (qlm_nghoststheta(hn)<2 .or. qlm_nghostsphi(hn)<2) call CCTK_WARN (0, "internal error") + + ! 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 + 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) + + nabla_ll(:,:) = qlm_tetrad_derivs(i,j)%nabla_ll(:,:) + nabla_nn(:,:) = qlm_tetrad_derivs(i,j)%nabla_nn(:,:) + nabla_mm(:,:) = qlm_tetrad_derivs(i,j)%nabla_mm(:,:) + + + + ! kappa = - m^a l^b D_b l_a + ! tau = - m^a n^b D_b l_a + ! sigma = - m^a m^b D_b l_a + ! rho = - m^a ~m^b D_b l_a [ ~m = congj(m) ] + qlm_npkappa (i,j,hn) = 0 + qlm_nptau (i,j,hn) = 0 + qlm_npsigma (i,j,hn) = 0 + qlm_nprho (i,j,hn) = 0 + ! epsilon = - 1/2 ( n^a l^b D_b l_a - ~m^a l^b D_b m_a ) + ! gamma = - 1/2 ( n^a n^b D_b l_a - ~m^a n^b D_b m_a ) + ! beta = - 1/2 ( n^a m^b D_b l_a - ~m^a m^b D_b m_a ) + ! alpha = - 1/2 ( n^a ~m^b D_b l_a - ~m^a ~m^b D_b m_a ) + qlm_npepsilon(i,j,hn) = 0 + qlm_npgamma (i,j,hn) = 0 + qlm_npbeta (i,j,hn) = 0 + qlm_npalpha (i,j,hn) = 0 + ! pi = ~m^a l^b D_b n_a + ! nu = ~m^a n^b D_b n_a + ! mu = ~m^a m^b D_b n_a + ! lambda = ~m^a ~m^b D_b n_a + qlm_nppi (i,j,hn) = 0 + qlm_npnu (i,j,hn) = 0 + qlm_npmu (i,j,hn) = 0 + qlm_nplambda (i,j,hn) = 0 + +!!$ qlm_lie_l_npsigma(i,j,hn) = 0 +!!$ qlm_lie_n_npsigma(i,j,hn) = 0 + +!!$ ! Theta is the expansion +!!$ ! Theta_(l) is (- 2 Re rho) +!!$ ! Theta_(n) is (+ 2 Re mu) +!!$ ! Theta_(T) is [Theta_(l) + Theta_(n)] / sqrt(2) +!!$ ! Theta_(S) is [Theta_(l) - Theta_(n)] / sqrt(2) + +!!$ qlm_lie_l_theta_l(i,j,hn) = 0 +!!$ qlm_lie_l_theta_n(i,j,hn) = 0 +!!$ qlm_lie_n_theta_l(i,j,hn) = 0 +!!$ qlm_lie_n_theta_n(i,j,hn) = 0 + + do a=0,3 + do b=0,3 + qlm_npkappa (i,j,hn) = qlm_npkappa (i,j,hn) - mm(a) * ll(b) * nabla_ll(a,b) + qlm_nptau (i,j,hn) = qlm_nptau (i,j,hn) - mm(a) * nn(b) * nabla_ll(a,b) + qlm_npsigma (i,j,hn) = qlm_npsigma (i,j,hn) - mm(a) * mm(b) * nabla_ll(a,b) + qlm_nprho (i,j,hn) = qlm_nprho (i,j,hn) - mm(a) * conjg(mm(b)) * nabla_ll(a,b) + qlm_npepsilon(i,j,hn) = qlm_npepsilon(i,j,hn) - half * (nn(a) * ll(b) * nabla_ll(a,b) - conjg(mm(a)) * ll(b) * nabla_mm(a,b)) + qlm_npgamma (i,j,hn) = qlm_npgamma (i,j,hn) - half * (nn(a) * nn(b) * nabla_ll(a,b) - conjg(mm(a)) * nn(b) * nabla_mm(a,b)) + qlm_npbeta (i,j,hn) = qlm_npbeta (i,j,hn) - half * (nn(a) * mm(b) * nabla_ll(a,b) - conjg(mm(a)) * mm(b) * nabla_mm(a,b)) + qlm_npalpha (i,j,hn) = qlm_npalpha (i,j,hn) - half * (nn(a) * conjg(mm(b)) * nabla_ll(a,b) - conjg(mm(a)) * conjg(mm(b)) * nabla_mm(a,b)) + qlm_nppi (i,j,hn) = qlm_nppi (i,j,hn) + conjg(mm(a)) * ll(b) * nabla_nn(a,b) + qlm_npnu (i,j,hn) = qlm_npnu (i,j,hn) + conjg(mm(a)) * nn(b) * nabla_nn(a,b) + qlm_npmu (i,j,hn) = qlm_npmu (i,j,hn) + conjg(mm(a)) * mm(b) * nabla_nn(a,b) + qlm_nplambda (i,j,hn) = qlm_nplambda (i,j,hn) + conjg(mm(a)) * conjg(mm(b)) * nabla_nn(a,b) + + end do + end do + + end do + end do + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npkappa (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_nptau (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npsigma (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_nprho (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npepsilon(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npgamma (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npbeta (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npalpha (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_nppi (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npnu (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_npmu (:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_nplambda (:,:,hn), +1) + +end subroutine qlm_calc_newman_penrose |