diff options
Diffstat (limited to 'src')
27 files changed, 5267 insertions, 0 deletions
diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..3205ee1 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,32 @@ +# Main make.code.defn file for thorn QuasiLocalMeasures -*-Makefile-*- + +# Source files in this directory +SRCS = qlm_3determinant.F90 \ + qlm_analyse.F90 \ + qlm_boundary.F90 \ + qlm_broadcast.c \ + qlm_calculate.F90 \ + qlm_coordinates.F90 \ + qlm_derivs.F90 \ + qlm_gram_schmidt.F90 \ + qlm_import_surface.F90 \ + qlm_init.F90 \ + qlm_interpolate.F90 \ + qlm_killing_axial.F90 \ + qlm_killing_gradient.F90 \ + qlm_killing_normalisation.F90 \ + qlm_killing_normalise.F90 \ + qlm_killing_test.F90 \ + qlm_killing_transport.F90 \ + qlm_killing_transportation.F90 \ + qlm_multipoles.F90 \ + qlm_newman_penrose.F90 \ + qlm_paramcheck.F90 \ + qlm_set_coordinates.F90 \ + qlm_tetrad.F90 \ + qlm_twometric.F90 \ + qlm_variables.F90 \ + qlm_weyl_scalars.F90 + +# Subdirectories containing source files +SUBDIRS = 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 diff --git a/src/qlm_analyse.F90 b/src/qlm_analyse.F90 new file mode 100644 index 0000000..44f0cba --- /dev/null +++ b/src/qlm_analyse.F90 @@ -0,0 +1,525 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_analyse (CCTK_ARGUMENTS, hn) + use adm_metric_simple + use cctk + use constants + use qlm_derivs + use qlm_variables + use tensor + use tensor2 + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: zero=0, two=2 + integer, parameter :: rk = kind(zero) + + CCTK_REAL :: theta, phi + CCTK_REAL :: ll(0:3), nn(0:3) + CCTK_COMPLEX :: mm(0:3) + CCTK_REAL :: ss(0:3), tt(0:3) + CCTK_COMPLEX :: npalpha, npbeta + CCTK_REAL :: gg(3,3), dgg(3,3,3), kk(3,3), gg_dot(3,3) + CCTK_REAL :: g4(0:3,0:3), dg4(0:3,0:3,0:3) + CCTK_REAL :: h4(0:3,0:3), dh4(0:3,0:3,0:3) + CCTK_REAL :: dtg, gu(3,3), trk + CCTK_REAL :: xx(3), ee(3,2) + CCTK_REAL :: xi(2), xi1(3) + CCTK_REAL :: qq(2,2), dtq + CCTK_REAL :: adm_energy, adm_mom(3), adm_amom(3,3) + CCTK_COMPLEX :: ev + CCTK_REAL :: spin + CCTK_REAL :: npspin + CCTK_REAL :: wsspin + CCTK_REAL :: tx, ty, tz + CCTK_REAL :: xi1_x(3), xi1_y(3), xi1_z(3) + CCTK_REAL :: coordspinx, coordspiny, coordspinz + + CCTK_REAL :: delta_space(2) + + CCTK_REAL, allocatable :: weights(:) + CCTK_REAL :: sum1 + + integer :: ntheta_inner, nphi_inner + + integer :: i_eq, j_p0, j_p2 + + integer :: i, j, l + integer :: a, b, c + + character :: msg*1000 + + if (veryverbose/=0) then + call CCTK_INFO ("Calculating spin") + end if + + delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /) + + ! Equatorial circumference + i_eq = (qlm_ntheta(hn) + 1) / 2 + + ! Polar circumference at phi=0 + j_p0 = 1+qlm_nghostsphi(hn) + + ! Polar circumference at phi=pi/2 + j_p2 = 1+qlm_nghostsphi(hn) + (qlm_nphi(hn) - 2*qlm_nghostsphi(hn) - 1) / 4 + + + + ! Initial values + qlm_equatorial_circumference(hn) = 0 + qlm_polar_circumference_0(hn) = 0 + qlm_polar_circumference_pi_2(hn) = 0 + qlm_area(hn) = 0 + qlm_spin(hn) = 0 + qlm_cvspin(hn) = 0 + qlm_npspin(hn) = 0 + qlm_wsspin(hn) = 0 + qlm_coordspinx(hn) = 0 + qlm_coordspiny(hn) = 0 + qlm_coordspinz(hn) = 0 + + qlm_adm_energy(hn) = 0 + qlm_adm_momentum_x(hn) = 0 + qlm_adm_momentum_y(hn) = 0 + qlm_adm_momentum_z(hn) = 0 + qlm_adm_angular_momentum_x(hn) = 0 + qlm_adm_angular_momentum_y(hn) = 0 + qlm_adm_angular_momentum_z(hn) = 0 + + ! Compute weights for spherical integration (see Driscoll and Healy). + ! These are the correct weights in a Gauss-Legendre-Senc. + ! Also compare qlm_area with AHFinderDirect's area, they are in much + ! better agreement than with the previous method using + ! delta_theta*delta_phi. + + ntheta_inner = qlm_ntheta(hn) - 2*qlm_nghoststheta(hn) + nphi_inner = qlm_nphi(hn) - 2*qlm_nghostsphi(hn) + + allocate (weights(1+qlm_nghoststheta(hn) : qlm_ntheta(hn)-qlm_nghoststheta(hn))) + + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + theta = qlm_origin_theta(hn) + (i-1)*qlm_delta_theta(hn) + sum1 = 0 + do l = 0, (ntheta_inner-1)/2 + sum1 = sum1 + sin((2*l+1)*theta)/(2*l+1) + end do + weights(i) = 8*pi * sum1 / (nphi_inner * ntheta_inner) + end do + + 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) + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + call calc_2det (qq, dtq) + + 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) + + tt = (ll + nn) / sqrt(two) + ss = (ll - nn) / sqrt(two) + + npalpha = qlm_npalpha(i,j,hn) + npbeta = qlm_npbeta(i,j,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(2,1,1) = dgg(1,2,1) + dgg(3,1,1) = dgg(1,3,1) + dgg(3,2,1) = dgg(2,3,1) + 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(2,1,2) = dgg(1,2,2) + dgg(3,1,2) = dgg(1,3,2) + dgg(3,2,2) = dgg(2,3,2) + 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,3) = dgg(1,2,3) + dgg(3,1,3) = dgg(1,3,3) + dgg(3,2,3) = dgg(2,3,3) + + kk(1,1) = qlm_kxx(i,j) + kk(1,2) = qlm_kxy(i,j) + kk(1,3) = qlm_kxz(i,j) + kk(2,2) = qlm_kyy(i,j) + kk(2,3) = qlm_kyz(i,j) + kk(3,3) = qlm_kzz(i,j) + kk(2,1) = kk(1,2) + kk(3,1) = kk(1,3) + kk(3,2) = kk(2,3) + + + + call calc_det (gg, dtg) + call calc_inv (gg, dtg, gu) + call calc_trace (gu, kk, trk) + + ! Calculate time derivatives + call calc_3metricdot_simple (kk, gg_dot) + + ! Calculate 4-metric + call calc_4metricderivs_simple (gg, dgg, gg_dot, g4,dg4) + + ! "Perturbative" 4-metric + h4 = g4 - eta4 + dh4 = dg4 + + xx(1) = qlm_x(i,j,hn) + xx(2) = qlm_y(i,j,hn) + xx(3) = qlm_z(i,j,hn) + + 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) + + xi(1) = qlm_xi_t(i,j,hn) + xi(2) = qlm_xi_p(i,j,hn) + + if (i == i_eq) then + qlm_equatorial_circumference(hn) = qlm_equatorial_circumference(hn) & + + sqrt(qq(2,2)) * qlm_delta_phi(hn) + end if + + if (j == j_p0) then + qlm_polar_circumference_0(hn) = qlm_polar_circumference_0(hn) & + + sqrt(qq(1,1)) * qlm_delta_theta(hn) + end if + + if (j == j_p2) then + qlm_polar_circumference_pi_2(hn) = qlm_polar_circumference_pi_2(hn) & + + sqrt(qq(1,1)) * qlm_delta_theta(hn) + end if + + qlm_area(hn) = qlm_area(hn) & + + sqrt(dtq) * weights(i) + + ! s^i: outward spacelike normal + ! K_ij: extrinsic curvature + ! phi^i: rotational Killing vector + + do a=1,3 + xi1(a) = 0 + do b=1,2 + xi1(a) = xi1(a) + ee(a,b) * xi(b) + end do + end do + + ! phi^i omega_i = - phi^i s^j K_ij + + spin = 0 + do a=1,3 + do b=1,3 + spin = spin + xi1(a) * ss(b) * kk(a,b) + end do + end do + + ! phi^i omega_i = (alpha + ~beta) phi^i m_i + complex conjugate + npspin = 0 + do a=1,3 + do b=1,3 + npspin = npspin + 2 * real((npalpha + conjg(npbeta)) * xi1(a) * gg(a,b) * mm(b)) + end do + end do + + ! phi^i omega_i = 1/2 f Im Psi_2 + ! (or is it phi^i omega_i = - f Im Psi_2 ?) + wsspin = - qlm_inv_z(i,j,hn) * aimag(qlm_psi2(i,j,hn)) + + ! x = sin theta cos phi + ! y = sin theta sin phi + ! z = cos theta + ! xi_x = (0,-z,y) + ! xi_y = (z,0,-x) + ! xi_z = (-y,x,0) + tx = qlm_x(i,j,hn) - qlm_origin_x(hn) + ty = qlm_y(i,j,hn) - qlm_origin_y(hn) + tz = qlm_z(i,j,hn) - qlm_origin_z(hn) + xi1_x(:) = (/ zero, -tz, ty /) + xi1_y(:) = (/ tz, zero, -tx /) + xi1_z(:) = (/ -ty, tx, zero /) + coordspinx = 0 + coordspiny = 0 + coordspinz = 0 + do a=1,3 + do b=1,3 + coordspinx = coordspinx + xi1_x(a) * ss(b) * kk(a,b) + coordspiny = coordspiny + xi1_y(a) * ss(b) * kk(a,b) + coordspinz = coordspinz + xi1_z(a) * ss(b) * kk(a,b) + end do + end do + + xi1_x(:) = (/ 1, 0, 0 /) + xi1_y(:) = (/ 0, 1, 0 /) + xi1_z(:) = (/ 0, 0, 1 /) + + qlm_spin(hn) = qlm_spin(hn) & + + spin * sqrt(dtq) * weights(i) + + qlm_npspin(hn) = qlm_npspin(hn) & + + npspin * sqrt(dtq) * weights(i) + + qlm_wsspin(hn) = qlm_wsspin(hn) & + + wsspin * sqrt(dtq) * weights(i) + + qlm_coordspinx(hn) = qlm_coordspinx(hn) & + + coordspinx * sqrt(dtq) * weights(i) + qlm_coordspiny(hn) = qlm_coordspiny(hn) & + + coordspiny * sqrt(dtq) * weights(i) + qlm_coordspinz(hn) = qlm_coordspinz(hn) & + + coordspinz * sqrt(dtq) * weights(i) + + ! ADM quantities + ! Weinberg, chapter 7.6, pp. 165 ff, eqns. (7.6.22) - (7.6.24): + + ! ADM energy + ! E_adm = (-1/16 pi) int_S(r) [h_jj,i - hij,j] n_i r^2 dOmega + adm_energy = 0 + do a=1,3 + do b=1,3 + adm_energy = adm_energy + (dh4(b,b,a) - dh4(a,b,b)) * ss(a) + end do + end do + qlm_adm_energy(hn) = qlm_adm_energy(hn) & + & + adm_energy / (-16*pi) & + & * sqrt(dtq) * weights(i) + + ! ADM momentum + ! P_adm^j = (-1/16 pi) int_S(r) + ! [- h_kk,0 delta_ij + h_k0,k delta_ij + ! - h_j0,i + h_ij,0 ] n_i r^2 dOmega + ! + ! P_adm^j = (1/8 pi) int_S(r) [K_ij - K delta_ij] n_i r^2 dOmega + do a=1,3 + adm_mom(a) = 0 + do b=1,3 + adm_mom(a) = adm_mom(a) & + + (- dh4(b,b,0) + dh4(b,0,b)) * ss(a) & + + (- dh4(a,0,b) + dh4(b,a,0)) * ss(b) + end do + end do + qlm_adm_momentum_x(hn) = qlm_adm_momentum_x(hn) & + & + adm_mom(1) / (-16*pi) & + & * sqrt(dtq) * weights(i) + qlm_adm_momentum_y(hn) = qlm_adm_momentum_y(hn) & + & + adm_mom(2) / (-16*pi) & + & * sqrt(dtq) * weights(i) + qlm_adm_momentum_z(hn) = qlm_adm_momentum_z(hn) & + & + adm_mom(3) / (-16*pi) & + & * sqrt(dtq) * weights(i) + + ! ADM angular momentum + ! J_adm^jk = (-1/16 pi) int_S(r) + ! [- x_j h_0k,i + x_k h_0j,i + ! + x_j h_ki,0 - x_k h_ji,0 + ! + h_0k delta_ij - h_0j delta_ik] n_i r^2 dOmega + do a=1,3 + do b=1,3 + adm_amom(a,b) = 0 + do c=1,3 + adm_amom(a,b) = adm_amom(a,b) + ss(c) * ( & + + (- xx(a) * dh4(0,b,c) + xx(b) * dh4(0,a,c)) & + + (+ xx(a) * dh4(b,c,0) - xx(b) * dh4(a,c,0)) & + + (h4(0,b) * delta4(c,a) - h4(0,a) * delta4(c,b))) + end do + end do + end do + qlm_adm_angular_momentum_x(hn) = qlm_adm_angular_momentum_x(hn) & + & + adm_amom(2,3) / (-16*pi) & + & * sqrt(dtq) * weights(i) + qlm_adm_angular_momentum_y(hn) = qlm_adm_angular_momentum_y(hn) & + & + adm_amom(3,1) / (-16*pi) & + & * sqrt(dtq) * weights(i) + qlm_adm_angular_momentum_z(hn) = qlm_adm_angular_momentum_z(hn) & + & + adm_amom(1,2) / (-16*pi) & + & * sqrt(dtq) * weights(i) + + end do + end do + + deallocate(weights) + + qlm_polar_circumference_0(hn) = qlm_polar_circumference_0(hn) * 2 + qlm_polar_circumference_pi_2(hn) = qlm_polar_circumference_pi_2(hn) * 2 + + ! A = 4 pi R^2 + ! R = 2 M + qlm_radius(hn) = sqrt(qlm_area(hn) / (4*pi)) + qlm_irreducible_mass(hn) = qlm_radius(hn) / 2 + + qlm_spin(hn) = qlm_spin(hn) / (8*pi) + qlm_mass(hn) = 1/(2*qlm_radius(hn)) * sqrt(qlm_radius(hn)**4 + 4*qlm_spin(hn)**2) + qlm_cvspin(hn) = qlm_cvspin(hn) / (8*pi) + qlm_npspin(hn) = qlm_npspin(hn) / (-8*pi) + qlm_wsspin(hn) = qlm_wsspin(hn) / (-4*pi) + qlm_coordspinx(hn) = qlm_coordspinx(hn) / (8*pi) + qlm_coordspiny(hn) = qlm_coordspiny(hn) / (8*pi) + qlm_coordspinz(hn) = qlm_coordspinz(hn) / (8*pi) + + ! The event horizon is at r = M + sqrt (M^2 - a^2) + ! with x^2 + y^2 + z^2 = rho^2 = r^2 + a^2 (1 - z^2 / r^2) + + call guess_mass_spin & + (qlm_area(hn), qlm_equatorial_circumference(hn), qlm_mass_guess(hn), qlm_spin_guess(hn)) + + + + if (verbose/=0) then + + write (msg, '("Geometric quantities for surface ",i4,":")') hn-1 + call CCTK_INFO (msg) + write (msg, '(" Area A: ",g16.6)') qlm_area(hn) + call CCTK_INFO (msg) + write (msg, '(" Irreducible mass M = R/2: ",g16.6)') qlm_irreducible_mass(hn) + call CCTK_INFO (msg) + write (msg, '(" Areal radius R = sqrt(A/4pi): ",g16.6)') qlm_radius(hn) + call CCTK_INFO (msg) + + write (msg, '("Coordinate-dependent quantities for surface ",i4,":")') hn-1 + call CCTK_INFO (msg) + write (msg, '(" Equatorial circumference: ",g16.6)') & + qlm_equatorial_circumference(hn) + call CCTK_INFO (msg) + write (msg, '(" Polar circumference at phi=0: ",g16.6)') & + qlm_polar_circumference_0(hn) + call CCTK_INFO (msg) + write (msg, '(" Polar circumference at phi=pi/2: ",g16.6)') & + qlm_polar_circumference_pi_2(hn) + call CCTK_INFO (msg) + if (qlm_spin_guess(hn) >= 0) then + write (msg, '(" Spin guess J from distortion: ",g16.6)') & + qlm_spin_guess(hn) + call CCTK_INFO (msg) + else + call CCTK_INFO (" No valid spin guess from distortion.") + call CCTK_INFO (" (Spin guess is imaginary.)") + write (msg, '(" Magnitude of invalid spin guess: ",g16.6)') & + abs(qlm_spin_guess(hn)) + call CCTK_INFO (msg) + end if + write (msg, '(" Mass guess M from distortion: ",g16.6)') & + qlm_mass_guess(hn) + call CCTK_INFO (msg) + + write (msg, '("Isolated Horizon quantities for surface ",i4,":")') hn-1 + call CCTK_INFO (msg) + ev = cmplx(qlm_killing_eigenvalue_re(hn), qlm_killing_eigenvalue_im(hn),rk) + write (msg, '(" Killing vector field eigenvalue norm: ",g14.6)') abs(ev) + call CCTK_INFO (msg) + write (msg, '(" Spin J: ",g14.6)') qlm_spin(hn) + call CCTK_INFO (msg) + write (msg, '(" Kerr spin parameter a = J/M: ",g14.6)') qlm_spin(hn) / qlm_mass(hn) + call CCTK_INFO (msg) + write (msg, '(" Dimensionless spin parameter a = J/M^2: ",g14.6)') qlm_spin(hn) / qlm_mass(hn)**2 + call CCTK_INFO (msg) + write (msg, '(" Spin J from NP: ",g14.6)') qlm_npspin(hn) + call CCTK_INFO (msg) + write (msg, '(" Spin J from phi-coordinate-vector: ",g14.6)') qlm_cvspin(hn) + call CCTK_INFO (msg) + write (msg, '(" Mass M: ",g14.6)') qlm_mass(hn) + call CCTK_INFO (msg) + + write (msg, '("Global quantities for surface ",i4,":")') hn-1 + call CCTK_INFO (msg) + write (msg, '(" ADM energy: ",g16.6)') qlm_adm_energy(hn) + call CCTK_INFO (msg) + write (msg, '(" ADM momentum x: ",g16.6)') qlm_adm_momentum_x(hn) + call CCTK_INFO (msg) + write (msg, '(" ADM momentum y: ",g16.6)') qlm_adm_momentum_y(hn) + call CCTK_INFO (msg) + write (msg, '(" ADM momentum z: ",g16.6)') qlm_adm_momentum_z(hn) + call CCTK_INFO (msg) + write (msg, '(" ADM angular momentum x: ",g16.6)') qlm_adm_angular_momentum_x(hn) + call CCTK_INFO (msg) + write (msg, '(" ADM angular momentum y: ",g16.6)') qlm_adm_angular_momentum_y(hn) + call CCTK_INFO (msg) + write (msg, '(" ADM angular momentum z: ",g16.6)') qlm_adm_angular_momentum_z(hn) + call CCTK_INFO (msg) + + end if + +contains + + subroutine guess_mass_spin (area, circumference, mass, spin) + CCTK_REAL, intent(in) :: area, circumference + CCTK_REAL, intent(out) :: mass, spin + CCTK_REAL :: radius, amom2, amom + + ! equatorial circumference L, area A + + ! L = 2 pi (r^2 + a^2) / r + ! A = 4 pi (r^2 + a^2) + ! r = M + sqrt (M^2 - a^2) + + ! r = A / (2 L) + ! a^2 = A / (4 pi) - r^2 ("spin" a = J/M = specific angular momentum) + ! M = (r^2 + a^2) / (2 r) + + ! J = a M (angular momentum) + + radius = area / (2*circumference) + amom2 = area / (4*pi) - radius**2 + amom = sign(sqrt(abs(amom2)), amom2) + mass = (radius**2 + amom2) / (2*radius) + spin = amom * mass + + ! equatorial circumference L, area A + ! mass M, areal radius R, spin a = J/M^2 + + ! A = 4 pi R^2 + ! L = 4 pi M + + ! a^2 = (R/M)^2 - 1/4 (R/M)^4 + ! J = a M^2 + + end subroutine guess_mass_spin + +end subroutine qlm_analyse diff --git a/src/qlm_boundary.F90 b/src/qlm_boundary.F90 new file mode 100644 index 0000000..33196bc --- /dev/null +++ b/src/qlm_boundary.F90 @@ -0,0 +1,154 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +module qlm_boundary + use cctk + implicit none + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + private + public set_boundary + + interface set_boundary + module procedure set_boundary_real + module procedure set_boundary_complex + end interface + +contains + + subroutine set_boundary_real (CCTK_ARGUMENTS, hn, f, parity) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + CCTK_REAL :: f(:,:) + integer, intent(in) :: parity + + CCTK_REAL, allocatable :: ff(:,:) + + integer :: ni, nj, gi, gj + + integer :: i, j, jj + + if (abs(parity) /= 1) then + call CCTK_WARN (0, "Parity must be +1 or -1") + end if + + + + ni = qlm_ntheta(hn) + nj = qlm_nphi(hn) + gi = qlm_nghoststheta(hn) + gj = qlm_nghostsphi(hn) + + + + f( :gi, :nj) = -42 + f(ni-gi+1:ni, :nj) = -42 + f(:ni, :gj) = -42 + f(:ni, nj-gi+1:nj) = -42 + + allocate (ff(gi, nj)) + + ! lower theta + ff(:, :nj) = f(1+gi:2*gi, :nj) + + ! polar boundary condition at north pole + do j=1,nj/2 + do i=1,gi + jj = j + (nj-2*gj) / 2 + f(i,j) = parity * ff(gi-i+1,jj) + end do + end do + do j=nj/2+1,nj + do i=1,gi + jj = j - (nj-2*gj) / 2 + f(i,j) = parity * ff(gi-i+1,jj) + end do + end do + + ! upper theta + ff(:, :nj) = f(ni-2*gi+1:ni-gi, :nj) + + ! polar boundary condition at south pole + do j=1,nj/2 + do i=1,gi + jj = j + (nj-2*gj) / 2 + f(ni-i+1,j) = parity * ff(i,jj) + end do + end do + do j=nj/2+1,nj + do i=1,gi + jj = j - (nj-2*gj) / 2 + f(ni-i+1,j) = parity * ff(i,jj) + end do + end do + + deallocate (ff) + + + + allocate (ff(ni,gj)) + + ! lower phi + ff(:ni, :) = f(:ni, nj-2*gj+1:nj-gj) + + ! periodic boundary at null meridian + do j=1,gj + do i=1,ni + f(i,j) = ff(i,j) + end do + end do + + ! upper phi + ff(:ni, :) = f(:ni, gj+1:2*gj) + + ! periodic boundary at null meridian + do j=1,gj + do i=1,ni + f(i,nj-j+1) = ff(i,gj-j+1) + end do + end do + + deallocate (ff) + + + + f(ni+1:, :nj) = 0 + f(:, nj+1:) = 0 + + end subroutine set_boundary_real + + + + subroutine set_boundary_complex (CCTK_ARGUMENTS, hn, f, parity) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + CCTK_COMPLEX :: f(:,:) + integer, intent(in) :: parity + + CCTK_REAL :: fre(size(f,1),size(f,2)) + CCTK_REAL :: fim(size(f,1),size(f,2)) + + integer :: ni, nj, gi, gj + + ni = qlm_ntheta(hn) + nj = qlm_nphi(hn) + gi = qlm_nghoststheta(hn) + gj = qlm_nghostsphi(hn) + + f( :gi, :nj) = 0 + f(ni-gi+1:ni, :nj) = 0 + f(:ni, :gj) = 0 + f(:ni, nj-gj+1:nj) = 0 + + fre = real(f) + fim = aimag(f) + call set_boundary_real (CCTK_PASS_FTOF, hn, fre, parity) + call set_boundary_real (CCTK_PASS_FTOF, hn, fim, parity) + f = cmplx(fre, fim, kind(f)) + end subroutine set_boundary_complex + +end module qlm_boundary diff --git a/src/qlm_broadcast.c b/src/qlm_broadcast.c new file mode 100644 index 0000000..a491fb5 --- /dev/null +++ b/src/qlm_broadcast.c @@ -0,0 +1,169 @@ +#include <assert.h> +#include <stdio.h> +#include <stdlib.h> + +#include <cctk.h> +#include <cctk_Arguments.h> +#include <cctk_Parameters.h> + +#ifdef CCTK_MPI +# include <mpi.h> +#endif + + + +/* Broadcast a vector element of a grid group */ +static void +bcast (cGH const * restrict const cctkGH, + char const * restrict const group, + int const vecind, + int const root) +{ +#ifdef CCTK_MPI + + int ierr; + + assert (group); + + MPI_Comm comm; + if (CCTK_IsFunctionAliased ("GetMPICommWorld")) + { + comm = * (MPI_Comm const *) GetMPICommWorld (cctkGH); + } + else + { + comm = MPI_COMM_WORLD; + } + + int num_procs; + MPI_Comm_size (comm, & num_procs); + + assert (root >= 0 && root < num_procs); + + int const gi = CCTK_GroupIndex (group); + assert (gi >= 0); + + cGroup data; + ierr = CCTK_GroupData (gi, & data); + assert (! ierr); + + assert (vecind >= 0 && vecind < data.vectorlength); + + if (data.numvars == 0) return; + + MPI_Datatype mpitype; + int items; + switch (data.vartype) + { + case CCTK_VARIABLE_INT: + assert (sizeof (CCTK_INT) == sizeof (int)); + mpitype = MPI_INT; + items = 1; + break; + case CCTK_VARIABLE_REAL: + assert (sizeof (CCTK_REAL) == sizeof (double)); + mpitype = MPI_DOUBLE; + items = 1; + break; + case CCTK_VARIABLE_COMPLEX: + assert (sizeof (CCTK_COMPLEX) == 2 * sizeof (double)); + mpitype = MPI_DOUBLE; + items = 2; + break; + default: + assert (0); + } + assert (items > 0); + + cGroupDynamicData dyndata; + ierr = CCTK_GroupDynamicData (cctkGH, gi, & dyndata); + assert (! ierr); + + int elems = 1; + { + int d; + for (d = 0; d < dyndata.dim; ++ d) + { + elems *= dyndata.lsh[d]; + } + } + assert (elems >= 0); + + int const v0 = CCTK_FirstVarIndexI (gi); + assert (v0 >= 0); + + int size; + MPI_Type_size (mpitype, & size); + assert (items * size == CCTK_VarTypeSize (data.vartype)); + + int const numvectors = data.numvars / data.vectorlength; + + { + int vec; + for (vec = 0; vec < numvectors; ++ vec) + { + int const vi = v0 + vecind + data.vectorlength * vec; + void * restrict const ptr = CCTK_VarDataPtrI (cctkGH, 0, vi); + assert (ptr); + + MPI_Bcast (ptr, elems * items, mpitype, root, comm); + } + } + +#endif /* ifdef CCTK_MPI */ +} + + +void CCTK_FCALL +CCTK_FNAME(qlm_broadcast) (CCTK_POINTER_TO_CONST * restrict cctkGH_); + +void CCTK_FCALL +CCTK_FNAME(qlm_broadcast) (CCTK_POINTER_TO_CONST * restrict const cctkGH_) +{ + cGH const * restrict const cctkGH = * (cGH const * const *) cctkGH_; + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (veryverbose) + { + CCTK_INFO ("Broadcasting quasi-local measures"); + } + + int const num_procs = CCTK_nProcs (cctkGH); + int hn; + for (hn = 0; hn < num_surfaces; ++ hn) + { + int const root = hn % num_procs; + + bcast (cctkGH, "QuasiLocalMeasures::qlm_state" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_state_p" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_grid_int" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_grid_real" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_grid_real_p", hn, root); + + bcast (cctkGH, "QuasiLocalMeasures::qlm_shapes" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_shapes_p" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_coordinates" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_l" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_n" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_m" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_l_p" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_n_p" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_m_p" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_newman_penrose" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_weyl_scalars" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_ricci_scalars" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_twometric" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_killing_vector" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_killed_twometric" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_invariant_coordinates", hn, root); + + bcast (cctkGH, "QuasiLocalMeasures::qlm_multipole_moments", hn, root); + + bcast (cctkGH, "QuasiLocalMeasures::qlm_3determinant", hn, root); + + bcast (cctkGH, "QuasiLocalMeasures::qlm_scalars" , hn, root); + bcast (cctkGH, "QuasiLocalMeasures::qlm_scalars_p", hn, root); + + } /* for hn */ +} diff --git a/src/qlm_calculate.F90 b/src/qlm_calculate.F90 new file mode 100644 index 0000000..eb4c4ba --- /dev/null +++ b/src/qlm_calculate.F90 @@ -0,0 +1,130 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_calculate (CCTK_ARGUMENTS) + use cctk + use qlm_variables + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + integer :: num_procs, my_proc + integer :: pass + integer :: h0, hn + + character :: msg*1000 + + logical :: did_allocate + + did_allocate = .false. + + num_procs = CCTK_nProcs (cctkGH) + my_proc = CCTK_MyProc (cctkGH) + + do pass = 1, (num_surfaces + num_procs - 1) / num_procs + + ! Calculate the range of horizons for this pass + h0 = (pass - 1) * num_procs + 1 + + ! This processor's horizon + hn = h0 + my_proc + + ! If there is nothing to do for this processor, set hn to zero + if (hn > num_surfaces) hn = 0 + + ! start calculations already? + if (hn > 0) then + if (cctk_time < begin_qlm_calculations_after(hn)) hn = 0 + end if + + if (verbose/=0 .or. veryverbose/=0) then + if (hn > 0) then + write (msg, '("Calculating Isolated and Dynamical Horizon quantities for horizon ",i4)') hn-1 + else + write (msg, '("Performing dummy calculation")') + end if + call CCTK_INFO (msg) + end if + + if (hn > 0) then + if (surface_index(hn) == -1) then + qlm_calc_error(hn) = 1 + qlm_have_valid_data(hn) = 0 + qlm_have_killing_vector(hn) = 0 + hn = 0 + end if + end if + + if (hn > 0) then + call qlm_import_surface (CCTK_PASS_FTOF, hn) + if (qlm_calc_error(hn) /= 0) hn = 0 + endif + + if (hn > 0) then + call qlm_set_coordinates (CCTK_PASS_FTOF, hn) + end if + + if (hn > 0) then + if (.not. did_allocate) then + ! allocate 2D arrays + call allocate_variables (int(maxntheta), int(maxnphi)) + did_allocate = .true. + end if + end if + + call qlm_interpolate (CCTK_PASS_FTOF, hn) + + if (hn > 0) then + if (qlm_calc_error(hn) /= 0) goto 9999 + + call qlm_calc_tetrad (CCTK_PASS_FTOF, hn) + call qlm_calc_newman_penrose (CCTK_PASS_FTOF, hn) + call qlm_calc_weyl_scalars (CCTK_PASS_FTOF, hn) + call qlm_calc_twometric (CCTK_PASS_FTOF, hn) + if (CCTK_EQUALS(killing_vector_method, "axial")) then + call qlm_killing_axial (CCTK_PASS_FTOF, hn) + else if (CCTK_EQUALS(killing_vector_method, "eigenvector")) then + call qlm_killing_transport (CCTK_PASS_FTOF, hn) + if (qlm_calc_error(hn) /= 0) goto 9999 + call qlm_killing_normalise (CCTK_PASS_FTOF, hn) + else if (CCTK_EQUALS(killing_vector_method, "gradient")) then + call qlm_killing_gradient (CCTK_PASS_FTOF, hn) + call qlm_killing_normalise (CCTK_PASS_FTOF, hn) + else + call CCTK_WARN (0, "internal error") + end if + if (qlm_calc_error(hn) /= 0) goto 9999 + if (qlm_have_killing_vector(hn) /= 0) then + call qlm_killing_test (CCTK_PASS_FTOF, hn) + call qlm_calc_coordinates (CCTK_PASS_FTOF, hn) + call qlm_multipoles (CCTK_PASS_FTOF, hn) + end if + call qlm_calc_3determinant (CCTK_PASS_FTOF, hn) + call qlm_analyse (CCTK_PASS_FTOF, hn) + +9999 continue + + if (qlm_timederiv_order(hn) < 2) then + call CCTK_WARN (2, "There were not enough past time levels available for accurate calculations") + end if + end if + + end do + + if (did_allocate) then + ! release 2D arrays + call deallocate_variables + end if + + call qlm_broadcast (cctkGH) + + if (veryverbose/=0) then + call CCTK_INFO ("Done.") + end if + +end subroutine qlm_calculate diff --git a/src/qlm_coordinates.F90 b/src/qlm_coordinates.F90 new file mode 100644 index 0000000..6807bcd --- /dev/null +++ b/src/qlm_coordinates.F90 @@ -0,0 +1,181 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_calc_coordinates (CCTK_ARGUMENTS, hn) + use cctk + use constants + use qlm_boundary + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: one=1, two=2, half=one/two + CCTK_REAL :: z0, z0dot, z1, z1dot + CCTK_REAL :: qq(2,2), dtq + CCTK_REAL :: integral_z, area, radius + integer :: i0,j0 + integer :: i,j + + if (veryverbose/=0) then + call CCTK_INFO ("Finding invariant coordinates") + end if + + ! latitude of "equator" + i0 = (qlm_ntheta(hn)+1)/2 + ! longitude of zero meridian + j0 = 1+qlm_nghostsphi(hn) + + + + ! calculate area + area = 0 + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1) + + area = area + sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + + end do + end do + radius = sqrt(area / (4*pi)) + + + + ! initial value + qlm_inv_z(i0,j0,hn) = 0 + + ! d_a z = 1/R^2 xi^b eps_ba + + ! transport along equator + do j = j0+1, qlm_nphi(hn)-qlm_nghostsphi(hn) + i = i0 + + z0 = qlm_inv_z(i,j-1,hn) + z0dot = rhs(i,j-1,(/0,1/)) + z1 = z0 + qlm_delta_phi(hn) * z0dot + z1dot = rhs(i,j,(/0,1/)) + qlm_inv_z(i,j,hn) = z0 + half * qlm_delta_phi(hn) * (z0dot + z1dot) + end do + + ! transport along meridians + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = i0-1, 1+qlm_nghoststheta(hn), -1 + z0 = qlm_inv_z(i+1,j,hn) + z0dot = rhs(i+1,j,(/-1,0/)) + z1 = z0 + qlm_delta_theta(hn) * z0dot + z1dot = rhs(i,j,(/-1,0/)) + qlm_inv_z(i,j,hn) = z0 + half * qlm_delta_theta(hn) * (z0dot + z1dot) + end do + + do i = i0+1, qlm_ntheta(hn)-qlm_nghoststheta(hn) + z0 = qlm_inv_z(i-1,j,hn) + z0dot = rhs(i-1,j,(/1,0/)) + z1 = z0 + qlm_delta_theta(hn) * z0dot + z1dot = rhs(i,j,(/1,0/)) + qlm_inv_z(i,j,hn) = z0 + half * qlm_delta_theta(hn) * (z0dot + z1dot) + end do + end do + + ! normalise + integral_z = 0 + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1) + + integral_z = integral_z + qlm_inv_z(i,j,hn) * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + + end do + end do + qlm_inv_z(:,:,hn) = qlm_inv_z(:,:,hn) - integral_z / area + + ! boundary conditions + call set_boundary (CCTK_PASS_FTOF, hn, qlm_inv_z(:,:,hn), +1) + + + +#if 0 + + ! initial value + qlm_inv_phi(i0,j0,hn) = 0 + + ! xi^a d_a phi = C + ! z^a d_a phi = 0 + ! z^a = (R^4 / q_bc xi^b xi^c) q^ab d_b z + ! v^a = A z^a + B xi^a + ! v^a d_a phi = (A z^a + B xi^a) d_a phi + ! = B C + ! (choose C=1, and normalise later) + +#error "replace z by phi" + ! transport along equator + do j = j0+1, qlm_nphi(hn)-qlm_nghostsphi(hn) + i = i0 + + z0 = qlm_inv_z(i,j-1,hn) + z0dot = rhs(i,j-1,(/0,1/)) + z1 = z0 + qlm_delta_phi(hn) * z0dot + z1dot = rhs(i,j,(/0,1/)) + qlm_inv_z(i,j,hn) = z0 + half * qlm_delta_phi(hn) * (z0dot + z1dot) + end do + + ! transport along meridians + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = i0-1, 1+qlm_nghoststheta(hn), -1 + z0 = qlm_inv_z(i+1,j,hn) + z0dot = rhs(i+1,j,(/-1,0/)) + z1 = z0 + qlm_delta_theta(hn) * z0dot + z1dot = rhs(i,j,(/-1,0/)) + qlm_inv_z(i,j,hn) = z0 + half * qlm_delta_theta(hn) * (z0dot + z1dot) + end do + + do i = i0+1, qlm_ntheta(hn)-qlm_nghoststheta(hn) + z0 = qlm_inv_z(i-1,j,hn) + z0dot = rhs(i-1,j,(/1,0/)) + z1 = z0 + qlm_delta_theta(hn) * z0dot + z1dot = rhs(i,j,(/1,0/)) + qlm_inv_z(i,j,hn) = z0 + half * qlm_delta_theta(hn) * (z0dot + z1dot) + end do + end do + +#error "normalise" + +#endif + +contains + + function rhs (i, j, vv) result (zdot) + integer, intent(in) :: i, j + integer, intent(in) :: vv(2) + CCTK_REAL :: zdot + CCTK_REAL :: qq(2,2), dtq + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1) + + zdot = vv(1) * (- (1/radius**2) * qlm_xi_p(i,j,hn) * sqrt(dtq)) & + + vv(2) * (+ (1/radius**2) * qlm_xi_t(i,j,hn) * sqrt(dtq)) + end function rhs + +end subroutine qlm_calc_coordinates diff --git a/src/qlm_derivs.F90 b/src/qlm_derivs.F90 new file mode 100644 index 0000000..24a3552 --- /dev/null +++ b/src/qlm_derivs.F90 @@ -0,0 +1,288 @@ +#include "cctk.h" +#include "cctk_Parameters.h" + +module qlm_derivs + use classify + implicit none + private + public abs2 + public operator(.outer.) + public operator(.dot.) + public deriv + public deriv2 + public deriv3 + public timederiv + public timederiv2 + + interface deriv + module procedure rderiv + end interface + + interface deriv2 + module procedure rderiv2 + end interface + + interface deriv3 + module procedure rderiv3 + end interface + + interface timederiv + module procedure rtimederiv + module procedure ctimederiv + end interface + + interface timederiv2 + module procedure rtimederiv2 + module procedure ctimederiv2 + end interface + + interface operator(.outer.) + module procedure router + module procedure couter + end interface + + interface operator(.dot.) + module procedure rdot + module procedure cdot + end interface + +contains + + ! abs(c)**2 for complex c without a square root + pure elemental function abs2 (a) + CCTK_COMPLEX, intent(in) :: a + CCTK_REAL :: abs2 + abs2 = a * conjg(a) + end function abs2 + + + + function router (left, right) result (result) + CCTK_REAL, intent(in) :: left(:), right(:) + CCTK_REAL :: result(size(left,1),size(right,1)) + integer :: i, j + forall (i=1:size(left,1), j=1:size(right,1)) + result(i,j) = left(i) * right(j) + end forall + end function router + + function couter (left, right) result (result) + CCTK_COMPLEX, intent(in) :: left(:), right(:) + CCTK_COMPLEX :: result(size(left,1),size(right,1)) + integer :: i, j + forall (i=1:size(left,1), j=1:size(right,1)) + result(i,j) = left(i) * right(j) + end forall + end function couter + + + + function rdot (left, right) result (result) + CCTK_REAL, intent(in) :: left(:), right(:) + CCTK_REAL :: result + integer :: i + if (size(left,1) /= size(right,1)) then + call CCTK_WARN (0, "Array sizes must have the same sizes") + end if +!!$ result = sum((/( left(i) * right(i), i=1,size(left,1) )/)) + result = 0 + do i=1,size(left,1) + result = result + left(i) * right(i) + end do + end function rdot + + function cdot (left, right) result (result) + CCTK_COMPLEX, intent(in) :: left(:), right(:) + CCTK_COMPLEX :: result + integer :: i + if (size(left,1) /= size(right,1)) then + call CCTK_WARN (0, "Array sizes must have the same sizes") + end if +!!$ result = sum((/( left(i) * right(i), i=1,size(left,1) )/)) + result = 0 + do i=1,size(left,1) + result = result + left(i) * right(i) + end do + end function cdot + + + + ! Calculate spatial derivatives + + pure function rderiv (f, i, j, dx) result (df) + DECLARE_CCTK_PARAMETERS + CCTK_REAL, intent(in) :: f(:,:) + integer, intent(in) :: i, j + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: df(2) + + select case (spatial_order) + case (2) + df(1) = (+ f(i+1,j) - f(i-1,j)) / (2 * dx(1)) + df(2) = (+ f(i,j+1) - f(i,j-1)) / (2 * dx(2)) + case (4) + df(1) = (- f(i+2,j) + 8*f(i+1,j) - 8*f(i-1,j) + f(i-2,j)) / (12 * dx(1)) + df(2) = (- f(i,j+2) + 8*f(i,j+1) - 8*f(i,j-1) + f(i,j-2)) / (12 * dx(2)) + case default + ! call CCTK_WARN (0, "internal error") + df = TAT_nan() + end select + end function rderiv + + pure function rderiv2 (f, i, j, dx) result (ddf) + DECLARE_CCTK_PARAMETERS + CCTK_REAL, intent(in) :: f(:,:) + integer, intent(in) :: i, j + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: ddf(2,2) + + select case (spatial_order) + case (2) + ddf(1,1) = (+ f(i+1,j) - 2*f(i,j) + f(i-1,j)) / dx(1)**2 + ddf(2,2) = (+ f(i,j+1) - 2*f(i,j) + f(i,j-1)) / dx(2)**2 + ddf(1,2) = (+ f(i-1,j-1) - f(i+1,j-1) - f(i-1,j+1) + f(i+1,j+1)) & + & / (4 * dx(1) * dx(2)) + ddf(2,1) = ddf(1,2) + case (4) + ddf(1,1) & + = (- f(i+2,j) + 16*f(i+1,j) - 30*f(i,j) + 16*f(i-1,j) - f(i-2,j)) & + & / (12 * dx(1)**2) + ddf(2,2) & + = (- f(i,j+2) + 16*f(i,j+1) - 30*f(i,j) + 16*f(i,j-1) - f(i,j-2)) & + & / (12 * dx(2)**2) + ddf(1,2) & + = (+ f(i+2,j+2) - 8*f(i+1,j+2) + 8*f(i-1,j+2) - f(i-2,j+2) & + & - 8*f(i+2,j+1) + 64*f(i+1,j+1) - 64*f(i-1,j+1) + 8*f(i-2,j+1) & + & + 8*f(i+2,j-1) - 64*f(i+1,j-1) + 64*f(i-1,j-1) - 8*f(i-2,j-1) & + & - f(i+2,j-2) + 8*f(i+1,j-2) - 8*f(i-1,j-2) + f(i-2,j-2)) & + & / (144 * dx(1) * dx(2)) + ddf(2,1) = ddf(1,2) + case default + ! call CCTK_WARN (0, "internal error") + ddf = TAT_nan() + end select + end function rderiv2 + + pure function rderiv3 (f, i, j, dx) result (dddf) + DECLARE_CCTK_PARAMETERS + CCTK_REAL, intent(in) :: f(:,:) + integer, intent(in) :: i, j + CCTK_REAL, intent(in) :: dx(2) + CCTK_REAL :: dddf(2,2,2) + + select case (spatial_order) + case (2, 4) + ! No separate 4th order stencil, since that would need 3 ghost + ! zones + dddf(1,1,1) = (- f(i-2,j ) & + & + 2*f(i-1,j ) & + & - 2*f(i+1,j ) & + & + f(i+2,j )) / (2*dx(1)**3) + dddf(1,1,2) = ( f(i+1,j+1) & + & - 2*f(i ,j+1) & + & + f(i-1,j+1) & + & - f(i+1,j-1) & + & + 2*f(i ,j-1) & + & - f(i-1,j-1)) / (2 * dx(1)**2 * dx(2)) + dddf(1,2,2) = ( f(i+1,j+1) & + & - 2*f(i+1,j ) & + & + f(i+1,j-1) & + & - f(i-1,j+1) & + & + 2*f(i-1,j ) & + & - f(i-1,j-1)) / (2 * dx(1) * dx(2)**2) + dddf(2,2,2) = (- f(i ,j-2) & + & + 2*f(i ,j-1) & + & - 2*f(i ,j+1) & + & + f(i ,j+2)) / (2*dx(2)**3) + dddf(1,2,1) = dddf(1,1,2) + dddf(2,1,1) = dddf(1,1,2) + dddf(2,1,2) = dddf(1,2,2) + dddf(2,2,1) = dddf(1,2,2) + case default + ! call CCTK_WARN (0, "internal error") + dddf = TAT_nan() + end select + end function rderiv3 + + + + ! Calculate a time derivate from several time levels + pure elemental function rtimederiv (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdot) + CCTK_REAL, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + logical, intent(in) :: ce0, ce1, ce2 + CCTK_REAL :: fdot + CCTK_REAL :: dt1, dt2 + CCTK_REAL :: fdot1, fdot2 + +!!$ dt1 = qlm_time(hn) - qlm_time_p(hn) +!!$ dt2 = qlm_time(hn) - qlm_time_p_p(hn) + dt1 = t0 - t1 + dt2 = t0 - t2 + + if (ce0 .or. ce1) then + fdot = 0 + else if (ce2) then + fdot = (f0 - f1) / dt1 + else + ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0) + ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0) + ! f'(0) = [f(dt1) - f(0)] / dt1 - dt1/2 f''(0) + ! f'(0) = [f(dt2) - f(0)] / dt2 - dt2/2 f''(0) + fdot1 = (f0 - f1) / dt1 + fdot2 = (f0 - f2) / dt2 + fdot = (fdot1 * dt2 - fdot2 * dt1) / (dt2 - dt1) + end if + end function rtimederiv + + pure elemental function ctimederiv (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdot) + CCTK_COMPLEX, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + logical, intent(in) :: ce0, ce1, ce2 + CCTK_COMPLEX :: fdot + + fdot = cmplx(timederiv(real(f0),real(f1),real(f2), t0,t1,t2, ce0,ce1,ce2), & + & timederiv(aimag(f0),aimag(f1),aimag(f2), t0,t1,t2, ce0,ce1,ce2), & + & kind(fdot)) + end function ctimederiv + + pure elemental function rtimederiv2 (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdotdot) + CCTK_REAL, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + logical, intent(in) :: ce0, ce1, ce2 + CCTK_REAL :: fdotdot + CCTK_REAL :: dt1, dt2 + CCTK_REAL :: fdotdot1, fdotdot2 + +!!$ dt1 = qlm_time(hn) - qlm_time_p(hn) +!!$ dt2 = qlm_time(hn) - qlm_time_p_p(hn) + dt1 = t0 - t1 + dt2 = t0 - t2 + + if (ce0 .or. ce1) then + fdotdot = 0 + else if (ce2) then + fdotdot = 0 + else + ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0) + ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0) + ! f''(0) = [f(dt1) - f(0)] / [dt1^2/2] - f'(0) / [dt1/2] + ! f''(0) = [f(dt2) - f(0)] / [dt2^2/2] - f'(0) / [dt2/2] + fdotdot1 = (f1 - f0) / (dt1**2/2) + fdotdot2 = (f2 - f0) / (dt2**2/2) + fdotdot = (fdotdot1 * dt1 - fdotdot2 * dt2) / (dt1 - dt2) + end if + end function rtimederiv2 + + pure elemental function ctimederiv2 (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdotdot) + CCTK_COMPLEX, intent(in) :: f0, f1, f2 + CCTK_REAL, intent(in) :: t0, t1, t2 + logical, intent(in) :: ce0, ce1, ce2 + CCTK_COMPLEX :: fdotdot + + fdotdot = cmplx(timederiv2(real(f0),real(f1),real(f2), t0,t1,t2, ce0,ce1,ce2), & + & timederiv2(aimag(f0),aimag(f1),aimag(f2), t0,t1,t2, ce0,ce1,ce2), & + & kind(fdotdot)) + end function ctimederiv2 + +end module qlm_derivs diff --git a/src/qlm_gram_schmidt.F90 b/src/qlm_gram_schmidt.F90 new file mode 100644 index 0000000..fb20293 --- /dev/null +++ b/src/qlm_gram_schmidt.F90 @@ -0,0 +1,91 @@ +#include "cctk.h" + +module qlm_gram_schmidt + implicit none + private + public gram_schmidt_project + public gram_schmidt_normalise + +contains + + subroutine gram_schmidt_project (g, dg, y, dy, y2, x, dx) + CCTK_REAL, intent(in) :: g(0:3,0:3), dg(0:3,0:3,0:3) + CCTK_REAL, intent(in) :: y(0:3), dy(0:3,0:3), y2 + CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3) + CCTK_REAL :: z(0:3), dz(0:3,0:3) + + integer :: a, b, d, e + + do a=0,3 + z(a) = x(a) + do d=0,3 + do e=0,3 + z(a) = z(a) - g(d,e) * x(d) * y(e) * y(a) / y2 + end do + end do + end do + + do a=0,3 + do b=0,3 + dz(a,b) = dx(a,b) + do d=0,3 + do e=0,3 + dz(a,b) = dz(a,b) & + - dg(d,e,b) * x(d) * y(e) * y(a) / y2 & + - g(d,e) * dx(d,b) * y(e) * y(a) / y2 & + - g(d,e) * x(d) * dy(e,b) * y(a) / y2 & + - g(d,e) * x(d) * y(e) * dy(a,b) / y2 + end do + end do + end do + end do + + x = z + dx = dz + end subroutine gram_schmidt_project + + subroutine gram_schmidt_normalise (g, dg, x, dx, x2) + CCTK_REAL, intent(in) :: g(0:3,0:3), dg(0:3,0:3,0:3) + CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3) + CCTK_REAL, intent(in) :: x2 + CCTK_REAL, parameter :: two=2, half=1/two + CCTK_REAL :: z2, dz2(0:3) + CCTK_REAL :: z(0:3), dz(0:3,0:3) + + integer :: a, b, d, e + + z2 = 0 + do d=0,3 + do e=0,3 + z2 = z2 + g(d,e) * x(d) * x(e) + end do + end do + + do a=0,3 + dz2(a) = 0 + do d=0,3 + do e=0,3 + dz2(a) = dz2(a) & + + dg(d,e,a) * x(d) * x(e) & + + g(d,e) * dx(d,a) * x(e) & + + g(d,e) * x(d) * dx(e,a) + end do + end do + end do + + do a=0,3 + z(a) = x(a) / sqrt(z2 / x2) + end do + + do a=0,3 + do b=0,3 + dz(a,b) = (+ dx(a,b) * z2 / x2 & + & - half * x(a) * dz2(b) / x2) / sqrt(z2 / x2)**3 + end do + end do + + x = z + dx = dz + end subroutine gram_schmidt_normalise + +end module qlm_gram_schmidt diff --git a/src/qlm_import_surface.F90 b/src/qlm_import_surface.F90 new file mode 100644 index 0000000..8eec453 --- /dev/null +++ b/src/qlm_import_surface.F90 @@ -0,0 +1,190 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_import_surface (CCTK_ARGUMENTS, hn) + use cctk + use qlm_boundary + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + integer :: i, j, sn + character :: msg*1000 + + if (veryverbose/=0) then + call CCTK_INFO ("Importing surface shape") + end if + + if (surface_index(hn) < 0 .or. surface_index(hn) >= nsurfaces) then + call CCTK_WARN (0, "Illegal spherical surface index specified") + end if + + if (verbose/=0 .or. veryverbose/=0) then + write (msg, '("Importing from spherical surface ",i4)') surface_index(hn) + call CCTK_INFO (msg) + end if + + sn = surface_index(hn) + 1 + + + + if (qlm_calc_error(hn) == 0 .and. cctk_iteration > qlm_iteration(hn)) then + + ! Cycle time levels + qlm_have_valid_data_p_p(hn) = qlm_have_valid_data_p(hn) + qlm_have_valid_data_p (hn) = qlm_have_valid_data (hn) + qlm_have_killing_vector_p_p(hn) = qlm_have_killing_vector_p(hn) + qlm_have_killing_vector_p (hn) = qlm_have_killing_vector (hn) + qlm_iteration(hn) = cctk_iteration + + qlm_time_p_p(hn) = qlm_time_p(hn) + qlm_time_p (hn) = qlm_time (hn) + qlm_radius_p_p(hn) = qlm_radius_p(hn) + qlm_radius_p (hn) = qlm_radius (hn) + + qlm_origin_x_p_p(hn) = qlm_origin_x_p(hn) + qlm_origin_x_p (hn) = qlm_origin_x (hn) + qlm_origin_y_p_p(hn) = qlm_origin_y_p(hn) + qlm_origin_y_p (hn) = qlm_origin_y (hn) + qlm_origin_z_p_p(hn) = qlm_origin_z_p(hn) + qlm_origin_z_p (hn) = qlm_origin_z (hn) + + qlm_shape_p_p(:,:,hn) = qlm_shape_p(:,:,hn) + qlm_shape_p (:,:,hn) = qlm_shape (:,:,hn) + + qlm_l0_p_p(:,:,hn) = qlm_l0_p(:,:,hn) + qlm_l0_p (:,:,hn) = qlm_l0 (:,:,hn) + qlm_l1_p_p(:,:,hn) = qlm_l1_p(:,:,hn) + qlm_l1_p (:,:,hn) = qlm_l1 (:,:,hn) + qlm_l2_p_p(:,:,hn) = qlm_l2_p(:,:,hn) + qlm_l2_p (:,:,hn) = qlm_l2 (:,:,hn) + qlm_l3_p_p(:,:,hn) = qlm_l3_p(:,:,hn) + qlm_l3_p (:,:,hn) = qlm_l3 (:,:,hn) + + qlm_n0_p_p(:,:,hn) = qlm_n0_p(:,:,hn) + qlm_n0_p (:,:,hn) = qlm_n0 (:,:,hn) + qlm_n1_p_p(:,:,hn) = qlm_n1_p(:,:,hn) + qlm_n1_p (:,:,hn) = qlm_n1 (:,:,hn) + qlm_n2_p_p(:,:,hn) = qlm_n2_p(:,:,hn) + qlm_n2_p (:,:,hn) = qlm_n2 (:,:,hn) + qlm_n3_p_p(:,:,hn) = qlm_n3_p(:,:,hn) + qlm_n3_p (:,:,hn) = qlm_n3 (:,:,hn) + + qlm_m0_p_p(:,:,hn) = qlm_m0_p(:,:,hn) + qlm_m0_p (:,:,hn) = qlm_m0 (:,:,hn) + qlm_m1_p_p(:,:,hn) = qlm_m1_p(:,:,hn) + qlm_m1_p (:,:,hn) = qlm_m1 (:,:,hn) + qlm_m2_p_p(:,:,hn) = qlm_m2_p(:,:,hn) + qlm_m2_p (:,:,hn) = qlm_m2 (:,:,hn) + qlm_m3_p_p(:,:,hn) = qlm_m3_p(:,:,hn) + qlm_m3_p (:,:,hn) = qlm_m3 (:,:,hn) + + qlm_xi_t_p_p(:,:,hn) = qlm_xi_t_p(:,:,hn) + qlm_xi_t_p (:,:,hn) = qlm_xi_t (:,:,hn) + qlm_xi_p_p_p(:,:,hn) = qlm_xi_p_p(:,:,hn) + qlm_xi_p_p (:,:,hn) = qlm_xi_p (:,:,hn) + + end if + + + + ! Check for valid horizon data + if (sf_valid(sn) <= 0) then + if (verbose/=0) then + call CCTK_INFO ("No valid horizon data found") + end if + qlm_calc_error(hn) = 1 + qlm_have_valid_data(hn) = 0 + qlm_have_killing_vector(hn) = 0 + return + end if + + + + ! Import surface description + qlm_calc_error(hn) = 0 + qlm_have_valid_data(hn) = 0 + qlm_have_killing_vector(hn) = 1 + + if (qlm_have_valid_data_p(hn) == 0) then + qlm_timederiv_order(hn) = 0 + else if (qlm_have_valid_data_p_p(hn) == 0) then + qlm_timederiv_order(hn) = 1 + else + qlm_timederiv_order(hn) = 2 + end if + + qlm_time(hn) = cctk_time + +#warning "TODO: Ensure that the surface parameters don't change" + qlm_nghoststheta(hn) = sf_nghoststheta(sn) + qlm_nghostsphi (hn) = sf_nghostsphi(sn) + qlm_ntheta (hn) = sf_ntheta(sn) + qlm_nphi (hn) = sf_nphi(sn) + + qlm_origin_x (hn) = sf_origin_x(sn) + qlm_origin_y (hn) = sf_origin_y(sn) + qlm_origin_z (hn) = sf_origin_z(sn) + qlm_origin_theta(hn) = sf_origin_theta(sn) + qlm_origin_phi (hn) = sf_origin_phi(sn) + qlm_delta_theta (hn) = sf_delta_theta(sn) + qlm_delta_phi (hn) = sf_delta_phi(sn) + + if (veryverbose /= 0) then + write (msg, '("calc error : ",i6)') qlm_calc_error(hn) + call CCTK_INFO (msg) + write (msg, '("time deriv order: ",i6)') qlm_timederiv_order(hn) + call CCTK_INFO (msg) + write (msg, '("time : ",g16.6)') qlm_time(hn) + call CCTK_INFO (msg) + write (msg, '("nghosts : ",2i6)') qlm_nghoststheta(hn), qlm_nghostsphi(hn) + call CCTK_INFO (msg) + write (msg, '("n : ",2i6)') qlm_ntheta(hn), qlm_nphi(hn) + call CCTK_INFO (msg) + write (msg, '("origin : ",3g16.6)') & + qlm_origin_x(hn), qlm_origin_y(hn), qlm_origin_z(hn) + call CCTK_INFO (msg) + write (msg, '("origin : ",2g16.6)') & + qlm_origin_theta(hn), qlm_origin_phi(hn) + call CCTK_INFO (msg) + write (msg, '("delta : ",2g16.6)') & + qlm_delta_theta(hn), qlm_delta_phi(hn) + call CCTK_INFO (msg) + end if + + if (qlm_ntheta(hn) > maxntheta .or. qlm_nphi(hn) > maxnphi) then + call CCTK_WARN (0, "Surface is too large") + end if + + if (qlm_nghoststheta(hn)<1 .or. qlm_nghostsphi(hn)<1) then + call CCTK_WARN (0, "Not enough ghost zones for the horizon surface -- need at least 1") + end if + + + ! Import the surface + ! Calculate the coordinates + do j = 1, qlm_nphi(hn) + do i = 1, qlm_ntheta(hn) + qlm_shape(i,j,hn) = sf_radius(i,j,sn) + end do + end do + + + + if (mod(int(qlm_ntheta(hn) - 2*qlm_nghoststheta(hn)),2) /= 1) then + ! We need a grid point on the equator + call CCTK_WARN (0, "The number of interior grid points in theta direction must be odd") + end if + + if (mod(int(qlm_nphi(hn) - 2*qlm_nghostsphi(hn)),4) /= 0) then + ! We need grid points on the four major meridians + call CCTK_WARN (0, "The number of interior grid points in phi direction must be a multiple of four") + end if + +end subroutine qlm_import_surface diff --git a/src/qlm_init.F90 b/src/qlm_init.F90 new file mode 100644 index 0000000..0106ad2 --- /dev/null +++ b/src/qlm_init.F90 @@ -0,0 +1,34 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_init (CCTK_ARGUMENTS) + use cctk + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + integer :: hn + + if (verbose/=0 .or. veryverbose/=0) then + call CCTK_INFO ("Initialising Quasi-Local Measure calculations") + end if + + do hn = 1, num_surfaces + + qlm_calc_error(hn) = 1 + qlm_have_valid_data(hn) = 0 + qlm_have_valid_data_p(hn) = 0 + qlm_have_valid_data_p_p(hn) = 0 + qlm_have_killing_vector(hn) = 0 + qlm_have_killing_vector_p(hn) = 0 + qlm_have_killing_vector_p_p(hn) = 0 + qlm_iteration(hn) = -1 + + end do + +end subroutine qlm_init diff --git a/src/qlm_interpolate.F90 b/src/qlm_interpolate.F90 new file mode 100644 index 0000000..82ec27e --- /dev/null +++ b/src/qlm_interpolate.F90 @@ -0,0 +1,602 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +! TODO: +! instead of interpolating to the symmetry points, copy them + + + +! A convenient shortcut +#define P(x) CCTK_PointerTo(x) + + + +subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) + use cctk + use qlm_variables + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_INT, parameter :: izero = 0 + integer, parameter :: ik = kind(izero) + integer, parameter :: sk = kind(interpolator) + CCTK_REAL, parameter :: one = 1 + + integer :: len_coordsystem + integer :: len_interpolator + integer :: len_interpolator_options + + character :: fort_coordsystem*100 + character :: fort_interpolator*100 + character :: fort_interpolator_options*1000 + + integer :: nvars + + integer :: coord_handle + integer :: interp_handle + integer :: options_table + + integer :: ninputs + integer :: noutputs + + CCTK_REAL, allocatable :: xcoord(:,:) + CCTK_REAL, allocatable :: ycoord(:,:) + CCTK_REAL, allocatable :: zcoord(:,:) + + integer :: ind_gxx, ind_gxy, ind_gxz, ind_gyy, ind_gyz, ind_gzz + integer :: ind_kxx, ind_kxy, ind_kxz, ind_kyy, ind_kyz, ind_kzz + integer :: ind_alpha + integer :: ind_betax, ind_betay, ind_betaz + + integer :: coord_type + CCTK_POINTER :: coords(3) + CCTK_INT :: inputs(16) + CCTK_INT :: output_types(88) + CCTK_POINTER :: outputs(88) + CCTK_INT :: operand_indices(88) + CCTK_INT :: operation_codes(88) + integer :: npoints + + character :: msg*1000 + + integer :: ni, nj + + integer :: ierr + + + + if (veryverbose/=0) then + call CCTK_INFO ("Interpolating 3d grid functions") + end if + + + + if (shift_state==0) then + call CCTK_WARN (0, "The shift must have storage") + end if + + + + ! Get coordinate system + call CCTK_FortranString & + (len_coordsystem, int(coordsystem,sk), fort_coordsystem) + call CCTK_CoordSystemHandle (coord_handle, fort_coordsystem) + if (coord_handle<0) then + write (msg, '("The coordinate system """, a, """ does not exist")') & + trim(fort_coordsystem) + call CCTK_WARN (0, msg) + end if + + ! Get interpolator + call CCTK_FortranString & + (len_interpolator, int(interpolator,sk), fort_interpolator) + call CCTK_InterpHandle (interp_handle, fort_interpolator) + if (interp_handle<0) then + write (msg, '("The interpolator """,a,""" does not exist")') & + trim(fort_interpolator) + call CCTK_WARN (0, msg) + end if + + ! Get interpolator options + call CCTK_FortranString & + (len_interpolator_options, int(interpolator_options,sk), & + fort_interpolator_options) + call Util_TableCreateFromString (options_table, fort_interpolator_options) + if (options_table<0) then + write (msg, '("The interpolator_options """,a,""" have a wrong syntax")') & + trim(fort_interpolator_options) + call CCTK_WARN (0, msg) + end if + + + + if (hn > 0) then + + ni = qlm_ntheta(hn) + nj = qlm_nphi(hn) + + allocate (xcoord(ni,nj)) + allocate (ycoord(ni,nj)) + allocate (zcoord(ni,nj)) + + xcoord(:,:) = qlm_x(:ni,:nj,hn) + ycoord(:,:) = qlm_y(:ni,:nj,hn) + zcoord(:,:) = qlm_z(:ni,:nj,hn) + + end if + + + + ! TODO: check the excision mask + + ! Get variable indices + call CCTK_VarIndex (ind_gxx , "ADMBase::gxx" ) + call CCTK_VarIndex (ind_gxy , "ADMBase::gxy" ) + call CCTK_VarIndex (ind_gxz , "ADMBase::gxz" ) + call CCTK_VarIndex (ind_gyy , "ADMBase::gyy" ) + call CCTK_VarIndex (ind_gyz , "ADMBase::gyz" ) + call CCTK_VarIndex (ind_gzz , "ADMBase::gzz" ) + call CCTK_VarIndex (ind_kxx , "ADMBase::kxx" ) + call CCTK_VarIndex (ind_kxy , "ADMBase::kxy" ) + call CCTK_VarIndex (ind_kxz , "ADMBase::kxz" ) + call CCTK_VarIndex (ind_kyy , "ADMBase::kyy" ) + call CCTK_VarIndex (ind_kyz , "ADMBase::kyz" ) + call CCTK_VarIndex (ind_kzz , "ADMBase::kzz" ) + call CCTK_VarIndex (ind_alpha, "ADMBase::alp" ) + call CCTK_VarIndex (ind_betax, "ADMBase::betax") + call CCTK_VarIndex (ind_betay, "ADMBase::betay") + call CCTK_VarIndex (ind_betaz, "ADMBase::betaz") + + + + ! Set up the interpolator arguments + coord_type = CCTK_VARIABLE_REAL + if (hn > 0) then + npoints = ni * nj + coords(:) = (/ P(xcoord), P(ycoord), P(zcoord) /) + else + npoints = 0 + coords(:) = CCTK_NullPointer() + end if + + inputs = (/ & + ind_gxx, ind_gxy, ind_gxz, ind_gyy, ind_gyz, ind_gzz, & + ind_kxx, ind_kxy, ind_kxz, ind_kyy, ind_kyz, ind_kzz, & + ind_alpha, & + ind_betax, ind_betay, ind_betaz /) + call CCTK_NumVars (nvars) + if (nvars < 0) call CCTK_WARN (0, "internal error") + if (any(inputs /= -1 .and. (inputs < 0 .or. inputs >= nvars))) then + call CCTK_WARN (0, "internal error") + end if + + operand_indices = (/ & + 00, 01, 02, 03, 04, 05, & ! g_ij + 00, 01, 02, 03, 04, 05, & ! g_ij,k + 00, 01, 02, 03, 04, 05, & + 00, 01, 02, 03, 04, 05, & + 00, 01, 02, 03, 04, 05, & ! g_ij,kl + 00, 01, 02, 03, 04, 05, & + 00, 01, 02, 03, 04, 05, & + 00, 01, 02, 03, 04, 05, & + 00, 01, 02, 03, 04, 05, & + 00, 01, 02, 03, 04, 05, & + 06, 07, 08, 09, 10, 11, & ! K_ij + 06, 07, 08, 09, 10, 11, & ! K_ij,k + 06, 07, 08, 09, 10, 11, & + 06, 07, 08, 09, 10, 11, & + 12, & ! alp + 13, 14, 15 /) ! beta^i + + operation_codes = (/ & + 0, 0, 0, 0, 0, 0, & ! g_ij + 1, 1, 1, 1, 1, 1, & ! g_ij,k + 2, 2, 2, 2, 2, 2, & + 3, 3, 3, 3, 3, 3, & + 11, 11, 11, 11, 11, 11, & ! g_ij,kl + 12, 12, 12, 12, 12, 12, & + 13, 13, 13, 13, 13, 13, & + 22, 22, 22, 22, 22, 22, & + 23, 23, 23, 23, 23, 23, & + 33, 33, 33, 33, 33, 33, & + 0, 0, 0, 0, 0, 0, & ! K_ij + 1, 1, 1, 1, 1, 1, & ! K_ij,k + 2, 2, 2, 2, 2, 2, & + 3, 3, 3, 3, 3, 3, & + 0, & ! alp + 0, 0, 0 /) ! beta^i + + output_types(:) = CCTK_VARIABLE_REAL + if (hn > 0) then + outputs = (/ & + P(qlm_gxx), P(qlm_gxy), P(qlm_gxz), P(qlm_gyy), P(qlm_gyz), P(qlm_gzz), & + P(qlm_dgxxx), P(qlm_dgxyx), P(qlm_dgxzx), P(qlm_dgyyx), P(qlm_dgyzx), P(qlm_dgzzx), & + P(qlm_dgxxy), P(qlm_dgxyy), P(qlm_dgxzy), P(qlm_dgyyy), P(qlm_dgyzy), P(qlm_dgzzy), & + P(qlm_dgxxz), P(qlm_dgxyz), P(qlm_dgxzz), P(qlm_dgyyz), P(qlm_dgyzz), P(qlm_dgzzz), & + P(qlm_ddgxxxx), P(qlm_ddgxyxx), P(qlm_ddgxzxx), P(qlm_ddgyyxx), P(qlm_ddgyzxx), P(qlm_ddgzzxx), & + P(qlm_ddgxxxy), P(qlm_ddgxyxy), P(qlm_ddgxzxy), P(qlm_ddgyyxy), P(qlm_ddgyzxy), P(qlm_ddgzzxy), & + P(qlm_ddgxxxz), P(qlm_ddgxyxz), P(qlm_ddgxzxz), P(qlm_ddgyyxz), P(qlm_ddgyzxz), P(qlm_ddgzzxz), & + P(qlm_ddgxxyy), P(qlm_ddgxyyy), P(qlm_ddgxzyy), P(qlm_ddgyyyy), P(qlm_ddgyzyy), P(qlm_ddgzzyy), & + P(qlm_ddgxxyz), P(qlm_ddgxyyz), P(qlm_ddgxzyz), P(qlm_ddgyyyz), P(qlm_ddgyzyz), P(qlm_ddgzzyz), & + P(qlm_ddgxxzz), P(qlm_ddgxyzz), P(qlm_ddgxzzz), P(qlm_ddgyyzz), P(qlm_ddgyzzz), P(qlm_ddgzzzz), & + P(qlm_kxx), P(qlm_kxy), P(qlm_kxz), P(qlm_kyy), P(qlm_kyz), P(qlm_kzz), & + P(qlm_dkxxx), P(qlm_dkxyx), P(qlm_dkxzx), P(qlm_dkyyx), P(qlm_dkyzx), P(qlm_dkzzx), & + P(qlm_dkxxy), P(qlm_dkxyy), P(qlm_dkxzy), P(qlm_dkyyy), P(qlm_dkyzy), P(qlm_dkzzy), & + P(qlm_dkxxz), P(qlm_dkxyz), P(qlm_dkxzz), P(qlm_dkyyz), P(qlm_dkyzz), P(qlm_dkzzz), & + P(qlm_alpha), & + P(qlm_betax), P(qlm_betay), P(qlm_betaz) /) + else + outputs(:) = CCTK_NullPointer() + end if + + + + ninputs = size(inputs) + noutputs = size(outputs) + + + +#if 0 + ! Poison the output variables +#define poison -42 + qlm_gxx = poison + qlm_gxy = poison + qlm_gxz = poison + qlm_gyy = poison + qlm_gyz = poison + qlm_gzz = poison + qlm_dgxxx = poison + qlm_dgxyx = poison + qlm_dgxzx = poison + qlm_dgyyx = poison + qlm_dgyzx = poison + qlm_dgzzx = poison + qlm_dgxxy = poison + qlm_dgxyy = poison + qlm_dgxzy = poison + qlm_dgyyy = poison + qlm_dgyzy = poison + qlm_dgzzy = poison + qlm_dgxxz = poison + qlm_dgxyz = poison + qlm_dgxzz = poison + qlm_dgyyz = poison + qlm_dgyzz = poison + qlm_dgzzz = poison + qlm_ddgxxxx = poison + qlm_ddgxyxx = poison + qlm_ddgxzxx = poison + qlm_ddgyyxx = poison + qlm_ddgyzxx = poison + qlm_ddgzzxx = poison + qlm_ddgxxxy = poison + qlm_ddgxyxy = poison + qlm_ddgxzxy = poison + qlm_ddgyyxy = poison + qlm_ddgyzxy = poison + qlm_ddgzzxy = poison + qlm_ddgxxxz = poison + qlm_ddgxyxz = poison + qlm_ddgxzxz = poison + qlm_ddgyyxz = poison + qlm_ddgyzxz = poison + qlm_ddgzzxz = poison + qlm_ddgxxyy = poison + qlm_ddgxyyy = poison + qlm_ddgxzyy = poison + qlm_ddgyyyy = poison + qlm_ddgyzyy = poison + qlm_ddgzzyy = poison + qlm_ddgxxyz = poison + qlm_ddgxyyz = poison + qlm_ddgxzyz = poison + qlm_ddgyyyz = poison + qlm_ddgyzyz = poison + qlm_ddgzzyz = poison + qlm_ddgxxzz = poison + qlm_ddgxyzz = poison + qlm_ddgxzzz = poison + qlm_ddgyyzz = poison + qlm_ddgyzzz = poison + qlm_ddgzzzz = poison + qlm_kxx = poison + qlm_kxy = poison + qlm_kxz = poison + qlm_kyy = poison + qlm_kyz = poison + qlm_kzz = poison + qlm_dkxxx = poison + qlm_dkxyx = poison + qlm_dkxzx = poison + qlm_dkyyx = poison + qlm_dkyzx = poison + qlm_dkzzx = poison + qlm_dkxxy = poison + qlm_dkxyy = poison + qlm_dkxzy = poison + qlm_dkyyy = poison + qlm_dkyzy = poison + qlm_dkzzy = poison + qlm_dkxxz = poison + qlm_dkxyz = poison + qlm_dkxzz = poison + qlm_dkyyz = poison + qlm_dkyzz = poison + qlm_dkzzz = poison + qlm_alpha = poison + qlm_betax = poison + qlm_betay = poison + qlm_betaz = poison +#endif + + + + ! Call the interpolator + call Util_TableSetIntArray & + (ierr, options_table, noutputs, & + operand_indices, "operand_indices") + if (ierr /= 0) call CCTK_WARN (0, "internal error") + call Util_TableSetIntArray & + (ierr, options_table, noutputs, & + operation_codes, "operation_codes") + if (ierr /= 0) call CCTK_WARN (0, "internal error") + + call CCTK_InterpGridArrays & + (ierr, cctkGH, 3, & + interp_handle, options_table, coord_handle, & + npoints, coord_type, coords, & + ninputs, inputs, & + noutputs, output_types, outputs) + + if (ierr /= 0) then + if (hn > 0) then + qlm_calc_error(hn) = 1 + end if + call CCTK_WARN (1, "Interpolator failed") + return + end if + + + + ! Unpack the variables + if (hn > 0) then + + call unpack (qlm_gxx , ni, nj) + call unpack (qlm_gxy , ni, nj) + call unpack (qlm_gxz , ni, nj) + call unpack (qlm_gyy , ni, nj) + call unpack (qlm_gyz , ni, nj) + call unpack (qlm_gzz , ni, nj) + call unpack (qlm_dgxxx , ni, nj) + call unpack (qlm_dgxyx , ni, nj) + call unpack (qlm_dgxzx , ni, nj) + call unpack (qlm_dgyyx , ni, nj) + call unpack (qlm_dgyzx , ni, nj) + call unpack (qlm_dgzzx , ni, nj) + call unpack (qlm_dgxxy , ni, nj) + call unpack (qlm_dgxyy , ni, nj) + call unpack (qlm_dgxzy , ni, nj) + call unpack (qlm_dgyyy , ni, nj) + call unpack (qlm_dgyzy , ni, nj) + call unpack (qlm_dgzzy , ni, nj) + call unpack (qlm_dgxxz , ni, nj) + call unpack (qlm_dgxyz , ni, nj) + call unpack (qlm_dgxzz , ni, nj) + call unpack (qlm_dgyyz , ni, nj) + call unpack (qlm_dgyzz , ni, nj) + call unpack (qlm_dgzzz , ni, nj) + call unpack (qlm_ddgxxxx , ni, nj) + call unpack (qlm_ddgxyxx , ni, nj) + call unpack (qlm_ddgxzxx , ni, nj) + call unpack (qlm_ddgyyxx , ni, nj) + call unpack (qlm_ddgyzxx , ni, nj) + call unpack (qlm_ddgzzxx , ni, nj) + call unpack (qlm_ddgxxxy , ni, nj) + call unpack (qlm_ddgxyxy , ni, nj) + call unpack (qlm_ddgxzxy , ni, nj) + call unpack (qlm_ddgyyxy , ni, nj) + call unpack (qlm_ddgyzxy , ni, nj) + call unpack (qlm_ddgzzxy , ni, nj) + call unpack (qlm_ddgxxxz , ni, nj) + call unpack (qlm_ddgxyxz , ni, nj) + call unpack (qlm_ddgxzxz , ni, nj) + call unpack (qlm_ddgyyxz , ni, nj) + call unpack (qlm_ddgyzxz , ni, nj) + call unpack (qlm_ddgzzxz , ni, nj) + call unpack (qlm_ddgxxyy , ni, nj) + call unpack (qlm_ddgxyyy , ni, nj) + call unpack (qlm_ddgxzyy , ni, nj) + call unpack (qlm_ddgyyyy , ni, nj) + call unpack (qlm_ddgyzyy , ni, nj) + call unpack (qlm_ddgzzyy , ni, nj) + call unpack (qlm_ddgxxyz , ni, nj) + call unpack (qlm_ddgxyyz , ni, nj) + call unpack (qlm_ddgxzyz , ni, nj) + call unpack (qlm_ddgyyyz , ni, nj) + call unpack (qlm_ddgyzyz , ni, nj) + call unpack (qlm_ddgzzyz , ni, nj) + call unpack (qlm_ddgxxzz , ni, nj) + call unpack (qlm_ddgxyzz , ni, nj) + call unpack (qlm_ddgxzzz , ni, nj) + call unpack (qlm_ddgyyzz , ni, nj) + call unpack (qlm_ddgyzzz , ni, nj) + call unpack (qlm_ddgzzzz , ni, nj) + call unpack (qlm_kxx , ni, nj) + call unpack (qlm_kxy , ni, nj) + call unpack (qlm_kxz , ni, nj) + call unpack (qlm_kyy , ni, nj) + call unpack (qlm_kyz , ni, nj) + call unpack (qlm_kzz , ni, nj) + call unpack (qlm_dkxxx , ni, nj) + call unpack (qlm_dkxyx , ni, nj) + call unpack (qlm_dkxzx , ni, nj) + call unpack (qlm_dkyyx , ni, nj) + call unpack (qlm_dkyzx , ni, nj) + call unpack (qlm_dkzzx , ni, nj) + call unpack (qlm_dkxxy , ni, nj) + call unpack (qlm_dkxyy , ni, nj) + call unpack (qlm_dkxzy , ni, nj) + call unpack (qlm_dkyyy , ni, nj) + call unpack (qlm_dkyzy , ni, nj) + call unpack (qlm_dkzzy , ni, nj) + call unpack (qlm_dkxxz , ni, nj) + call unpack (qlm_dkxyz , ni, nj) + call unpack (qlm_dkxzz , ni, nj) + call unpack (qlm_dkyyz , ni, nj) + call unpack (qlm_dkyzz , ni, nj) + call unpack (qlm_dkzzz , ni, nj) + call unpack (qlm_alpha , ni, nj) + call unpack (qlm_betax , ni, nj) + call unpack (qlm_betay , ni, nj) + call unpack (qlm_betaz , ni, nj) + + + +#if 0 + ! Check for poison + if (any(qlm_gxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_gxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_gxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_gyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_gyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_gzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxyx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxzx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgyyx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgyzx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgzzx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxzy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgyyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgyzy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgzzy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgxzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgyyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgyzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dgzzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxxxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxyxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxzxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyyxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyzxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgzzxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxxxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxyxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxzxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyyxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyzxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgzzxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxxxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxyxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxzxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyyxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyzxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgzzxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxxyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxyyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxzyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyyyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyzyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgzzyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxxyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxyyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxzyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyyyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyzyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgzzyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxxzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxyzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgxzzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyyzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgyzzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_ddgzzzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_kxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_kxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_kxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_kyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_kyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_kzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxxx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxyx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxzx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkyyx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkyzx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkzzx == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxxy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxzy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkyyy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkyzy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkzzy == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxxz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkxzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkyyz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkyzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_dkzzz == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_alpha == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_betax == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_betay == poison)) call CCTK_WARN (0, "poison found") + if (any(qlm_betaz == poison)) call CCTK_WARN (0, "poison found") +#endif + + end if + + + + ! Free interpolator options + call Util_TableDestroy (ierr, options_table) + + + + if (hn > 0) then + + qlm_have_valid_data(hn) = 1 + + deallocate (xcoord) + deallocate (ycoord) + deallocate (zcoord) + + end if + + + +contains + + subroutine pack (arr, ni, nj) + integer, intent(in) :: ni, nj + CCTK_REAL, intent(inout) :: arr(:,:) + CCTK_REAL :: tmp(ni,nj) + tmp(:,:) = arr(:ni, :nj) + call copy (arr, tmp, size(tmp)) + end subroutine pack + + subroutine unpack (arr, ni, nj) + integer, intent(in) :: ni, nj + CCTK_REAL, intent(inout) :: arr(:,:) + CCTK_REAL :: tmp(ni,nj) + call copy (tmp, arr, size(tmp)) + arr(:ni, :nj) = tmp(:,:) + arr(ni+1:, :nj) = 0 + arr(:, nj+1:) = 0 + end subroutine unpack + + subroutine copy (a, b, n) + integer, intent(in) :: n + CCTK_REAL, intent(out) :: a(n) + CCTK_REAL, intent(in) :: b(n) + a = b + end subroutine copy + +end subroutine qlm_interpolate diff --git a/src/qlm_killing_axial.F90 b/src/qlm_killing_axial.F90 new file mode 100644 index 0000000..feda999 --- /dev/null +++ b/src/qlm_killing_axial.F90 @@ -0,0 +1,35 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_killing_axial (CCTK_ARGUMENTS, hn) + use cctk + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL :: xi(2), chi + integer :: i, j + + do j = 1, qlm_nphi(hn) + do i = 1, qlm_ntheta(hn) + + xi(1) = 0 + xi(2) = 1 + + chi = 0 + + qlm_xi_t(i,j,hn) = xi(1) + qlm_xi_p(i,j,hn) = xi(2) + + qlm_chi(i,j,hn) = chi + + end do + end do + +end subroutine qlm_killing_axial diff --git a/src/qlm_killing_gradient.F90 b/src/qlm_killing_gradient.F90 new file mode 100644 index 0000000..fe970b5 --- /dev/null +++ b/src/qlm_killing_gradient.F90 @@ -0,0 +1,169 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_killing_gradient (CCTK_ARGUMENTS, hn) + use cctk + use constants + use qlm_boundary + use qlm_derivs + use tensor2 + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: two=2, half=1/two + CCTK_REAL :: qq(2,2), dqq(2,2,2), dtq, qu(2,2), dqu(2,2,2) + CCTK_REAL :: dpsi2(2), ddpsi2(2,2), ndpsi2, dndpsi2(2) + CCTK_REAL :: xi(2), dxi(2,2), chi + integer :: i, j + integer :: a, b + CCTK_REAL :: delta_space(2) + + delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /) + + ! Calculate the gradient of a scalar + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + +#if 0 + dqq(1,1,1) = qlm_dqttt(i,j,hn) + dqq(1,1,2) = qlm_dqttp(i,j,hn) + dqq(1,2,1) = qlm_dqtpt(i,j,hn) + dqq(1,2,2) = qlm_dqtpp(i,j,hn) + dqq(2,2,1) = qlm_dqppt(i,j,hn) + dqq(2,2,2) = qlm_dqppp(i,j,hn) + dqq(2,1,:) = dqq(1,2,:) +#endif + + call calc_2det (qq, dtq) +#if 0 + call calc_2inv (qq, dtq, qu) + call calc_2invderiv (qu, dqq, dqu) +#endif + + dpsi2(1) = (abs2(qlm_psi2(i+1,j,hn)) - abs2(qlm_psi2(i-1,j,hn))) / (2*qlm_delta_theta(hn)) + dpsi2(2) = (abs2(qlm_psi2(i,j+1,hn)) - abs2(qlm_psi2(i,j-1,hn))) / (2*qlm_delta_phi(hn)) + +#if 0 + ddpsi2(1,1) = (abs2(qlm_psi2(i+1,j,hn)) - 2*abs2(qlm_psi2(i,j,hn)) + abs2(qlm_psi2(i-1,j,hn))) / qlm_delta_theta(hn)**2 + ddpsi2(2,2) = (abs2(qlm_psi2(i,j+1,hn)) - 2*abs2(qlm_psi2(i,j,hn)) + abs2(qlm_psi2(i,j-1,hn))) / qlm_delta_phi(hn)**2 + ddpsi2(1,1) = (abs2(qlm_psi2(i-1,j-1,hn)) - abs2(qlm_psi2(i+1,j-1,hn)) - abs2(qlm_psi2(i-1,j+1,hn)) + abs2(qlm_psi2(i+1,j+1,hn))) / (4*qlm_delta_theta(hn)*qlm_delta_phi(hn)) + ddpsi2(2,1) = ddpsi2(1,2) + + ! ndpsi2 = ||grad |Psi_2|^2|| + ndpsi2 = 0 + do a=1,2 + do b=1,2 + ndpsi2 = ndpsi2 + qu(a,b) * dpsi2(a) * dpsi2(b) + end do + end do + ndpsi2 = sqrt(ndpsi2) + + ! dndpsi2 = grad ||grad |Psi_2|^2|| + do a=1,2 + dndpsi2(a) = 0 + do b=1,2 + do c=1,2 + dndpsi2(a) = dndpsi2(a) + 1 / (2*ndpsi2) * (qu(b,c) * ddpsi2(b,a) * dpsi2(c) + qu(b,c) * dpsi2(b) * ddpsi2(c,a) + dqu(b,c,a) * dpsi2(b) * dpsi2(c)) + end do + end do + end do +#endif + + ! xi^a = eps^ab D_b |Psi_2|^2 + do a=1,2 + xi(a) = 0 + do b=1,2 + xi(a) = xi(a) + sqrt(dtq) * epsilon2(a,b) * dpsi2(b) + end do + end do + + qlm_xi_t(i,j,hn) = xi(1) + qlm_xi_p(i,j,hn) = xi(2) + +#if 0 + ! xi^a = eps^ab D_b ||D_c |Psi_2|^2|| + do a=1,2 + xi(a) = 0 + do b=1,2 + xi(a) = xi(a) + sqrt(dtq) * epsilon2(a,b) * dndpsi2(b) + end do + end do + + qlm_xi_t(i,j,hn) = xi(1) + qlm_xi_p(i,j,hn) = xi(2) +#endif + + end do + end do + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_t(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_p(:,:,hn), -1) + + + + ! fix up xi (which must not be zero) + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + + xi(1) = qlm_xi_t(i,j,hn) + xi(2) = qlm_xi_p(i,j,hn) + + if (sum(xi**2) < 1.0d-4**2) then + + qlm_xi_t(i,j,hn) = sum(qlm_xi_t(i:i+1,j:j+1,hn)) / 4 + qlm_xi_p(i,j,hn) = sum(qlm_xi_p(i:i+1,j:j+1,hn)) / 4 + + end if + + end do + end do + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_t(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_p(:,:,hn), -1) + + + + ! set up the derivative of xi (which is not really needed) + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + call calc_2det (qq, dtq) + + dxi(1,1:2) = deriv (qlm_xi_t(:,:,hn), i, j, delta_space) + dxi(2,1:2) = deriv (qlm_xi_p(:,:,hn), i, j, delta_space) + + ! eps_ab sqrt(q) chi = D_b xi_a + ! sqrt(q) chi = -1/2 eps^ab D_a xi_b + chi = 0 + do a=1,2 + do b=1,2 + chi = chi - half * sqrt(dtq) * epsilon2(a,b) * dxi(b,a) + end do + end do + + qlm_chi(i,j,hn) = chi + + end do + end do + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_chi (:,:,hn), +1) + +end subroutine qlm_killing_gradient diff --git a/src/qlm_killing_normalisation.F90 b/src/qlm_killing_normalisation.F90 new file mode 100644 index 0000000..608fff8 --- /dev/null +++ b/src/qlm_killing_normalisation.F90 @@ -0,0 +1,226 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +module qlm_killing_normalisation + use cctk + use constants + implicit none + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + private + public killing_factor + +contains + + subroutine killing_factor (CCTK_ARGUMENTS, hn, theta, factor, nsteps) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + CCTK_REAL, intent(in) :: theta + CCTK_REAL, intent(out) :: factor + integer, intent(out) :: nsteps + + CCTK_REAL :: lambda0, theta0, phi0 + CCTK_REAL :: lambda1, theta1, phi1 + character :: msg*1000 + + lambda0 = 0 + theta0 = theta + phi0 = 0 + phi1 = 2*pi + call killing_geodesic & + (CCTK_PASS_FTOF, hn, lambda0, theta0, phi0, phi1, lambda1, theta1, & + nsteps) + + factor = (lambda1 - lambda0) / (2*pi) + + if (veryverbose/=0) then + write (msg, '("Integrated at theta=",g16.6)') theta + call CCTK_INFO (msg) + write (msg, '(" Integrated in ",i4," steps")') nsteps + call CCTK_INFO (msg) + write (msg, '(" Theta error is ",g16.6)') theta1 - theta0 + call CCTK_INFO (msg) + write (msg, '(" Normalisation factor is ",g16.6)') factor + call CCTK_INFO (msg) + end if + end subroutine killing_factor + + + + subroutine killing_geodesic & + (CCTK_ARGUMENTS, hn, lambda0, theta0, phi0, phi1, lambda1, theta1, nsteps) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + CCTK_REAL, intent(in) :: lambda0, theta0, phi0, phi1 + CCTK_REAL, intent(out) :: lambda1, theta1 + integer, intent(out) :: nsteps + + CCTK_REAL :: org_theta, org_phi, del_theta, del_phi + + CCTK_REAL :: lambda + CCTK_REAL :: theta, phi + CCTK_REAL :: dlambda + CCTK_REAL :: dtheta, dphi + CCTK_REAL :: theta2, phi2 + + integer :: ierr1, ierr2 + + org_theta = qlm_origin_theta(hn) + org_phi = qlm_origin_phi(hn) + del_theta = qlm_delta_theta(hn) + del_phi = qlm_delta_phi(hn) + + nsteps = 0 + lambda = lambda0 + theta = theta0 + phi = phi0 + + dtheta = killing_interp (qlm_xi_t(:,:,hn), & + org_theta, org_phi, del_theta, del_phi, theta, phi, ierr1) + dphi = killing_interp (qlm_xi_p(:,:,hn), & + org_theta, org_phi, del_theta, del_phi, theta, phi, ierr2) + + if (ierr1/=0 .or. ierr2/=0) then + call CCTK_WARN (2, "Integration path leaves the domain") + lambda1 = lambda0 + theta1 = theta0 + nsteps = -1 + return + end if + + if (abs(dphi) < 1.0d-8 .or. abs(dtheta) > abs(dphi)) then + call CCTK_WARN (2, "Integration path starts out too steep") + lambda1 = lambda0 + theta1 = theta0 + nsteps = -1 + return + end if + + dlambda = (qlm_delta_phi(hn) / dphi) / 2 + + do + + dtheta = killing_interp (qlm_xi_t(:,:,hn), & + org_theta, org_phi, del_theta, del_phi, theta, phi, ierr1) + dphi = killing_interp (qlm_xi_p(:,:,hn), & + org_theta, org_phi, del_theta, del_phi, theta, phi, ierr2) + + if (ierr1/=0 .or. ierr2/=0) then + call CCTK_WARN (2, "Integration path leaves the domain") + lambda1 = lambda0 + theta1 = theta0 + nsteps = -1 + return + end if + + theta2 = theta + dlambda * dtheta / 2 + phi2 = phi + dlambda * dphi / 2 + + dtheta = killing_interp (qlm_xi_t(:,:,hn), & + org_theta, org_phi, del_theta, del_phi, theta2, phi2, ierr1) + dphi = killing_interp (qlm_xi_p(:,:,hn), & + org_theta, org_phi, del_theta, del_phi, theta2, phi2, ierr2) + + if (ierr1/=0 .or. ierr2/=0) then + call CCTK_WARN (2, "Integration path leaves the domain") + lambda1 = lambda0 + theta1 = theta0 + nsteps = -1 + return + end if + + if (dphi<=0) then + call CCTK_WARN (2, "Integration path does not enclose the pole") + lambda1 = lambda0 + theta1 = theta0 + nsteps = -1 + return + end if + + theta2 = theta + dlambda * dtheta + phi2 = phi + dlambda * dphi + + if (phi2 >= phi1) exit + + if (nsteps > 100000) then + call CCTK_WARN (2, "Integration takes too many steps") + lambda1 = lambda0 + theta1 = theta0 + nsteps = -1 + return + exit + end if + + nsteps = nsteps + 1 + lambda = lambda + dlambda + theta = theta2 + phi = phi2 + + end do + + dlambda = (phi1 - phi) / dphi + + nsteps = nsteps + 1 + lambda = lambda + dlambda + theta = theta + dlambda * dtheta + phi = phi + dlambda * dphi + + lambda1 = lambda + theta1 = theta + + end subroutine killing_geodesic + + + + function killing_interp & + (array, origin_x1, origin_x2, delta_x1, delta_x2, x1, x2, ierr) + CCTK_REAL :: killing_interp + CCTK_REAL, intent(in) :: array(:,:) + CCTK_REAL, intent(in) :: origin_x1, origin_x2 + CCTK_REAL, intent(in) :: delta_x1, delta_x2 + CCTK_REAL, intent(in) :: x1, x2 + integer, intent(out) :: ierr + CCTK_REAL, parameter :: eps = 1.0d-10 + CCTK_REAL :: xx1, xx2 + CCTK_REAL :: dx1, dx2 + CCTK_REAL :: f1a, f1b, f1c + CCTK_REAL :: f2a, f2b, f2c + CCTK_REAL :: interp2a, interp2b, interp2c + integer :: i, j + xx1 = x1 + xx2 = x2 + i = nint((xx1 - origin_x1) / delta_x1) + 1 + j = nint((xx2 - origin_x2) / delta_x2) + 1 + i = max(2, min(size(array,1)-1, i)) + if (j>size(array,2)-1) then + ! periodicity in phi-direction + xx2 = xx2 - (size(array,2)-2) * delta_x2 + j = j - size(array,2)-2 + end if + j = max(2, min(size(array,2)-1, j)) + dx1 = xx1 - (origin_x1 + (i-1) * delta_x1) + dx2 = xx2 - (origin_x2 + (j-1) * delta_x2) + if (abs(dx1)>(1+eps)*delta_x1/2 .or. abs(dx2)>(1+eps)*delta_x2/2) then + call CCTK_WARN (2, "interpolating out of bounds") + killing_interp = 0 + ierr = 1 + return + end if + f1a = (-dx1) * (1-dx1) / ( 2 * 1 ) + f1b = (-1-dx1) * (1-dx1) / ( 1 * (-1)) + f1c = (-1-dx1) * (-dx1) / ((-1) * (-2)) + f2a = (-dx2) * (1-dx2) / ( 2 * 1 ) + f2b = (-1-dx2) * (1-dx2) / ( 1 * (-1)) + f2c = (-1-dx2) * (-dx2) / ((-1) * (-2)) + interp2a = f1a * array(i-1,j-1) + f1b * array(i,j-1) + f1c * array(i+1,j-1) + interp2b = f1a * array(i-1,j ) + f1b * array(i,j ) + f1c * array(i+1,j ) + interp2c = f1a * array(i-1,j+1) + f1b * array(i,j+1) + f1c * array(i+1,j+1) + killing_interp = f2a * interp2a + f2b * interp2b + f2c * interp2c + ierr = 0 + end function killing_interp + +end module qlm_killing_normalisation diff --git a/src/qlm_killing_normalise.F90 b/src/qlm_killing_normalise.F90 new file mode 100644 index 0000000..0d044f1 --- /dev/null +++ b/src/qlm_killing_normalise.F90 @@ -0,0 +1,141 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_killing_normalise (CCTK_ARGUMENTS, hn) + use cctk + use qlm_killing_normalisation + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: one = 1 + + integer, parameter :: ngeodesics = 7 + + CCTK_REAL :: theta + CCTK_REAL :: factor, factors(ngeodesics) + logical :: found(ngeodesics) + + integer :: nsteps + integer :: ii, iii + integer :: i, j + + character :: msg*1000 + + if (veryverbose/=0) then + call CCTK_INFO ("Normalising Killing vector field") + end if + + + + ! Starting meridian + j = 1+qlm_nghostsphi(hn) + + + + if (veryverbose/=0) call CCTK_INFO ("Calculating normalisation factor:") + +#if 0 + nsteps = -1 + do iii=1,ngeodesics + ii = (ngeodesics/2+1) + (2*mod(iii,2)-1) * (iii/2) + if (ii<1 .or. ii>ngeodesics) call CCTK_WARN (0, "internal error") + i = 1 + qlm_nghoststheta(hn) & + + ii * (qlm_ntheta(hn) - 2*qlm_nghoststheta(hn) - 1) / (ngeodesics+1) + if (qlm_xi_p(i,j,hn) < 0) then + qlm_xi_t(:,:,hn) = -qlm_xi_t(:,:,hn) + qlm_xi_p(:,:,hn) = -qlm_xi_p(:,:,hn) + qlm_chi (:,:,hn) = -qlm_chi (:,:,hn) + end if + theta = qlm_origin_theta(hn) + (i-1) * qlm_delta_theta(hn) + call killing_factor (CCTK_PASS_FTOF, hn, theta, factor, nsteps) + if (nsteps>0) exit + end do + if (nsteps < 0) then + call CCTK_WARN (1, "Did not manage to integrate along a Killing vector field line loop") + ! qlm_calc_error(hn) = 1 + qlm_have_killing_vector(hn) = 0 + factor = 1 + goto 9999 + end if +#else + do ii=1,ngeodesics + i = 1 + qlm_nghoststheta(hn) & + + ii * (qlm_ntheta(hn) - 2*qlm_nghoststheta(hn) - 1) / (ngeodesics+1) + if (qlm_xi_p(i,j,hn) < 0) then + qlm_xi_t(:,:,hn) = -qlm_xi_t(:,:,hn) + qlm_xi_p(:,:,hn) = -qlm_xi_p(:,:,hn) + qlm_chi (:,:,hn) = -qlm_chi (:,:,hn) + end if + theta = qlm_origin_theta(hn) + (i-1) * qlm_delta_theta(hn) + call killing_factor (CCTK_PASS_FTOF, hn, theta, factors(ii), nsteps) + found(ii) = nsteps>0 + end do + if (any(found)) then + if (CCTK_EQUALS(killing_vector_normalisation, "average")) then + if (count(found) >= 3) then + ! ignore smallest and largest factors + ii = minloc(factors, 1, found) + if (ii<=0) call CCTK_WARN (0, "internal error") + found(ii) = .false. + ii = maxloc(factors, 1, found) + if (ii<=0) call CCTK_WARN (0, "internal error") + found(ii) = .false. + end if + else if (CCTK_EQUALS(killing_vector_normalisation, "median")) then + do while (count(found) >= 3) + ! ignore smallest and largest factors + ii = minloc(factors, 1, found) + if (ii<=0) call CCTK_WARN (0, "internal error") + found(ii) = .false. + ii = maxloc(factors, 1, found) + if (ii<=0) call CCTK_WARN (0, "internal error") + found(ii) = .false. + end do + else + call CCTK_WARN (0, "internal error") + end if + ! average factors + factor = product(factors, found) ** (one / count(found)) + else + call CCTK_WARN (1, "Did not manage to integrate along a Killing vector field line loop") + ! qlm_calc_error(hn) = 1 + qlm_have_killing_vector(hn) = 0 + factor = 1 + goto 9999 + end if +#endif + + + + if (veryverbose/=0) then + write (msg, '("Normalising xi with the factor ",g16.6)') factor + call CCTK_INFO (msg) + end if + qlm_xi_t(:,:,hn) = qlm_xi_t(:,:,hn) * factor + qlm_xi_p(:,:,hn) = qlm_xi_p(:,:,hn) * factor + qlm_chi (:,:,hn) = qlm_chi (:,:,hn) * factor + + + + if (veryverbose/=0) then + call CCTK_INFO ("Checking normalisation factors for various Killing vector field line loops:") + do ii=1,ngeodesics + i = 1 + qlm_nghoststheta(hn) & + + ii * (qlm_ntheta(hn) - 2*qlm_nghoststheta(hn) - 1) / (ngeodesics+1) + theta = qlm_origin_theta(hn) + (i-1) * qlm_delta_theta(hn) + call killing_factor (CCTK_PASS_FTOF, hn, theta, factor, nsteps) + end do + end if + + ! qlm_have_killing_vector(hn) = 1 + +9999 continue + +end subroutine qlm_killing_normalise diff --git a/src/qlm_killing_test.F90 b/src/qlm_killing_test.F90 new file mode 100644 index 0000000..d268824 --- /dev/null +++ b/src/qlm_killing_test.F90 @@ -0,0 +1,83 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_killing_test (CCTK_ARGUMENTS, hn) + use cctk + use qlm_boundary + use qlm_derivs + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL :: qq(2,2), dqq(2,2,2) + CCTK_REAL :: xi(2), dxi(2,2) + CCTK_REAL :: lqq(2,2) + + CCTK_REAL :: delta_space(2) + + integer :: i,j + integer :: a, b, c + CCTK_REAL :: theta, phi + + if (veryverbose/=0) then + call CCTK_INFO ("Testing Killing vector field") + end if + + 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) + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + + dqq(1,1,1) = qlm_dqttt(i,j,hn) + dqq(1,2,1) = qlm_dqtpt(i,j,hn) + dqq(2,2,1) = qlm_dqppt(i,j,hn) + dqq(1,1,2) = qlm_dqttp(i,j,hn) + dqq(1,2,2) = qlm_dqtpp(i,j,hn) + dqq(2,2,2) = qlm_dqppp(i,j,hn) + dqq(2,1,:) = dqq(1,2,:) + + xi(1) = qlm_xi_t(i,j,hn) + xi(2) = qlm_xi_p(i,j,hn) + + dxi(1,1:2) = deriv (qlm_xi_t(:,:,hn), i, j, delta_space) + dxi(2,1:2) = deriv (qlm_xi_p(:,:,hn), i, j, delta_space) + + ! L_xi q_ab = q_cb d_a xi^c + q_ac d_b xi^c + xi^c d_c q_ab + do a=1,2 + do b=1,2 + lqq(a,b) = 0 + do c=1,2 + lqq(a,b) = lqq(a,b) + qq(c,b) * dxi(c,a) & + & + qq(a,c) * dxi(c,b) & + & + xi(c) * dqq(a,b,c) + end do + end do + end do + + qlm_lqtt(i,j,hn) = lqq(1,1) + qlm_lqtp(i,j,hn) = lqq(1,2) + qlm_lqpp(i,j,hn) = lqq(2,2) + + end do + end do + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_lqtt(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_lqtp(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_lqpp(:,:,hn), +1) + +end subroutine qlm_killing_test diff --git a/src/qlm_killing_transport.F90 b/src/qlm_killing_transport.F90 new file mode 100644 index 0000000..3ec4f7e --- /dev/null +++ b/src/qlm_killing_transport.F90 @@ -0,0 +1,130 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_killing_transport (CCTK_ARGUMENTS, hn) + use cctk + use constants + use qlm_boundary + use qlm_killing_transportation + use lapack + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + interface + integer function TAT_isnan (x) + implicit none + CCTK_REAL x + end function TAT_isnan + end interface + + CCTK_REAL :: xi(2,3), chi(3) + CCTK_REAL :: vec(3,3) + CCTK_REAL :: wr(3), wi(3), vl(3,3), vr(3,3) + integer :: i0, j0 + integer :: n + integer :: info + character :: msg*1000 + + integer, parameter :: lwork = 100 + CCTK_REAL :: work(lwork) + + + + if (veryverbose/=0) then + call CCTK_INFO ("Transporting Killing vector field") + end if + + + + ! latitude of "equator" + i0 = (qlm_ntheta(hn)+1)/2 + ! longitude of zero meridian + j0 = 1+qlm_nghostsphi(hn) + + vec = delta3 + + if (veryverbose/=0) then + write (msg, '("Initial vectors (xi^t xi^p chi):")') + call CCTK_INFO (msg) + do n=1,3 + write (msg, '(i2,3g18.8)') n, vec(:,n) + call CCTK_INFO (msg) + end do + end if + + xi(1,:) = vec(1,:) + xi(2,:) = vec(2,:) + chi(:) = vec(3,:) + + do n=1,3 + call transport_along_equator (CCTK_PASS_FTOF, hn, i0, xi(:,n), chi(n)) + end do + + vec(1,:) = xi(1,:) + vec(2,:) = xi(2,:) + vec(3,:) = chi(:) + + if (veryverbose/=0) then + write (msg, '("Final vectors (xi^t xi^p chi):")') + call CCTK_INFO (msg) + do n=1,3 + write (msg, '(i2,3g18.8)') n, vec(:,n) + call CCTK_INFO (msg) + end do + end if + + if (TAT_isnan(sum(vec)) /= 0) then + ! qlm_calc_error(hn) = 1 + qlm_have_killing_vector(hn) = 0 + call CCTK_WARN (3, "There are nans in the final vectors") + goto 9999 + end if + + call geev ('n', 'v', 3, vec, 3, wr, wi, vl, 3, vr, 3, work, lwork, info) + if (info/=0) then + ! qlm_calc_error(hn) = 1 + qlm_have_killing_vector(hn) = 0 + write (msg, '("Error in call to GEEV, info=",i2)') info + call CCTK_WARN (3, msg) + goto 9999 + end if + + if (veryverbose/=0) then + write (msg, '("Eigenvectors and eigenvalues:")') + call CCTK_INFO (msg) + do n=1,3 + write (msg, '(i2,3g18.8,(" (",g18.8,",",g18.8,")"))') & + n, vr(:,n), wr(n), wi(n) + call CCTK_INFO (msg) + end do + end if + + xi(1,:) = vr(1,:) + xi(2,:) = vr(2,:) + chi(:) = vr(3,:) + + ! TODO: + ! This transport scheme is not ideal. + ! It leads to large fluctuations in the phi direction. + n=3 + qlm_killing_eigenvalue_re(hn) = wr(n) + qlm_killing_eigenvalue_im(hn) = wi(n) + if (abs(cmplx(wr(n),wi(n),kind(wr)) - (1,0)) > 1.0d-4) then + call CCTK_WARN (3, "Did not manage to find an eigenvector with the eigenvalue 1") + end if + call transport_along_equator (CCTK_PASS_FTOF, hn, i0, xi(:,n), chi(n)) + call transport_along_meridians (CCTK_PASS_FTOF, hn, i0) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_t(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_p(:,:,hn), -1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_chi (:,:,hn), +1) + +9999 continue +end subroutine qlm_killing_transport diff --git a/src/qlm_killing_transportation.F90 b/src/qlm_killing_transportation.F90 new file mode 100644 index 0000000..779e524 --- /dev/null +++ b/src/qlm_killing_transportation.F90 @@ -0,0 +1,250 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +module qlm_killing_transportation + use cctk + use constants + use ricci2 + use tensor2 + implicit none + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + private + public transport_along_equator + public transport_along_meridians + +contains + + subroutine transport_along_equator (CCTK_ARGUMENTS, hn, i0, xi, chi) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + integer, intent(in) :: i0 + CCTK_REAL, intent(inout) :: xi(2), chi + integer :: j0 + integer :: nsteps + + j0 = 1+qlm_nghostsphi(hn) + nsteps = qlm_nphi(hn) - 2*qlm_nghostsphi(hn) + + call transport (CCTK_PASS_FTOF, hn, i0, j0, 0, 1, nsteps, xi, chi) + + end subroutine transport_along_equator + + + + subroutine transport_along_meridians (CCTK_ARGUMENTS, hn, i0) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + integer, intent(in) :: i0 + CCTK_REAL :: xi(2), chi + integer :: j0 + integer :: dir + integer :: nsteps + + do j0 = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + + do dir=-1,+1,2 + + if (dir==-1) nsteps = i0 - (1+qlm_nghoststheta(hn)) + if (dir==+1) nsteps = (qlm_ntheta(hn)-qlm_nghoststheta(hn)) - i0 + + xi(1) = qlm_xi_t(i0,j0,hn) + xi(2) = qlm_xi_p(i0,j0,hn) + chi = qlm_chi(i0,j0,hn) + + call transport (CCTK_PASS_FTOF, hn, i0, j0, dir, 0, nsteps, xi, chi) + + end do + + end do + + end subroutine transport_along_meridians + + + + subroutine transport (CCTK_ARGUMENTS, hn, i0, j0, di, dj, nsteps, xi, chi) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + integer, intent(in) :: i0, j0 + integer, intent(in) :: di, dj + integer, intent(in) :: nsteps + CCTK_REAL,intent(inout) :: xi(2), chi + CCTK_REAL :: vv(2) + CCTK_REAL :: xi_dot(2), chi_dot + CCTK_REAL :: xi1(2), chi1 + CCTK_REAL :: xi1_dot(2), chi1_dot + integer :: n + integer :: i, j + + i = i0 + j = j0 + + vv(1) = di * qlm_delta_theta(hn) + vv(2) = dj * qlm_delta_phi(hn) + + do n=1,nsteps + + call transport_rhs (CCTK_PASS_FTOF, hn, i, j, xi, chi, vv, xi_dot, chi_dot) + xi1 = xi + xi_dot + chi1 = chi + chi_dot + i = i + di + j = j + dj + if (j < 1+qlm_nghostsphi(hn)) j = j + (qlm_nphi(hn)-2*qlm_nghostsphi(hn)) + if (j > qlm_nphi(hn)-qlm_nghostsphi(hn)) j = j - (qlm_nphi(hn)-2*qlm_nghostsphi(hn)) + + call transport_rhs & + (CCTK_PASS_FTOF, hn, i, j, xi1, chi1, vv, xi1_dot, chi1_dot) + xi = xi + 0.5d0 * (xi_dot + xi1_dot) + chi = chi + 0.5d0 * (chi_dot + chi1_dot) + + qlm_xi_t(i,j,hn) = xi(1) + qlm_xi_p(i,j,hn) = xi(2) + qlm_chi(i,j,hn) = chi + + end do + + end subroutine transport + + + + subroutine transport_rhs (CCTK_ARGUMENTS, hn, i, j, xi, chi, vv, xi_dot, chi_dot) + DECLARE_CCTK_ARGUMENTS + integer, intent(in) :: hn + integer, intent(in) :: i, j + CCTK_REAL, intent(in) :: xi(2), chi + CCTK_REAL, intent(in) :: vv(2) + CCTK_REAL, intent(out) :: xi_dot(2), chi_dot + CCTK_REAL :: qq(2,2), dqq(2,2,2), dtq, qu(2,2), gamma(2,2,2), rsc + + if (i<1+qlm_nghoststheta(hn) .or. i>qlm_ntheta(hn)-qlm_nghoststheta(hn) & + .or. j<1+qlm_nghostsphi(hn) .or. j>qlm_nphi(hn)-qlm_nghostsphi(hn)) then + call CCTK_WARN (0, "internal error") + end if + if (i-1<1 .or. i+1>qlm_ntheta(hn) .or. j-1<1 .or. j+1>qlm_nphi(hn)) then + call CCTK_WARN (0, "internal error") + end if + + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + + dqq(1,1,1) = qlm_dqttt(i,j,hn) + dqq(1,2,1) = qlm_dqtpt(i,j,hn) + dqq(2,2,1) = qlm_dqppt(i,j,hn) + dqq(1,1,2) = qlm_dqttp(i,j,hn) + dqq(1,2,2) = qlm_dqtpp(i,j,hn) + dqq(2,2,2) = qlm_dqppp(i,j,hn) + dqq(2,1,:) = dqq(1,2,:) + + rsc = qlm_rsc(i,j,hn) + + call calc_2det (qq, dtq) + call calc_2inv (qq, dtq, qu) + + call calc_2connections (qu, dqq, gamma) + + call killing_transport_rhs & + (xi, chi, qq, dtq, qu, gamma, rsc, vv, xi_dot, chi_dot) + + end subroutine transport_rhs + + + + subroutine killing_transport_rhs & + (xi, chi, qq, dtq, qu, gamma2, rsc2, vv, xi_dot, chi_dot) + CCTK_REAL, intent(in) :: xi(2), chi + CCTK_REAL, intent(in) :: qq(2,2), dtq, qu(2,2), gamma2(2,2,2), rsc2 + CCTK_REAL, intent(in) :: vv(2) + CCTK_REAL, intent(out) :: xi_dot(2), chi_dot + integer :: i, k, l + + ! Wald eqn (C.3.6): + ! D_k D_j xi^i = R^i_jkl xi^l + + ! define: + ! L^i_j = D_j x^i + + ! then: see Wald eqns (C.3.7) and (C.3.8): + ! v^k D_k xi^i = L^i_k v^k + ! v^k D_k L^i_j = R^i_jkl v^k xi^l + + ! in 2D we have: + ! R_ijkl = 1/2 q R epsilon2_ij epsilon2_kl + ! L_ij = epsilon2_ij sqrt(q) chi + + ! then: + ! v^k D_k xi^i = epsilon2^i_k sqrt(q) chi v^k + ! v^k D_k epsilon2^i_j sqrt(q) chi = R^i_jkl v^k xi^l + ! v^k D_k chi = 1/2 sqrt(q) R epsilon2_kl v^k xi^l + + ! define: + ! X_dot = v^i d/dx^i X (partial derivatives) + + do i=1,2 + xi_dot(i) = 0 + do k=1,2 + do l=1,2 + xi_dot(i) = xi_dot(i) & + + qu(i,l) * epsilon2(l,k) * sqrt(dtq) * chi * vv(k) & + - vv(k) * gamma2(i,l,k) * xi(l) + end do + end do + end do + + chi_dot = 0 + do k=1,2 + do l=1,2 + chi_dot = chi_dot & + + 0.5d0 * sqrt(dtq) * rsc2 * epsilon2(k,l) * vv(k) * xi(l) + end do + end do + + end subroutine killing_transport_rhs + + + +#if 0 + subroutine killing_equation (qq, qu, xi, gxi, zeta, trzeta, zetasq) + CCTK_REAL, intent(in) :: qq(2,2), qu(2,2) + CCTK_REAL, intent(in) :: xi(2), gxi(2,2) + CCTK_REAL, intent(out) :: zeta(2,2), trzeta, zetasq + integer :: i, j, k, l + + do i=1,2 + do j=1,2 + ! zeta_ij = D_i xi_j + D_j xi_i + zeta(i,j) = 0 + do k=1,2 + do l=1,2 + zeta(i,j) = zeta(i,j) + qq(j,k) * gxi(k,i) + qq(i,k) * gxi(k,j) + end do + end do + end do + end do + + trzeta = 0 + do i=1,2 + do j=1,2 + trzeta = trzeta + qu(i,j) * zeta(i,j) + end do + end do + + zetasq = 0 + do i=1,2 + do j=1,2 + do k=1,2 + do l=1,2 + zetasq = zetasq + qu(i,k) * qu(j,l) * zeta(i,j) * zeta(k,l) + end do + end do + end do + end do + end subroutine killing_equation +#endif + +end module qlm_killing_transportation diff --git a/src/qlm_multipoles.F90 b/src/qlm_multipoles.F90 new file mode 100644 index 0000000..79b75d5 --- /dev/null +++ b/src/qlm_multipoles.F90 @@ -0,0 +1,323 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_multipoles (CCTK_ARGUMENTS, hn) + use cctk + use constants + use qlm_boundary + use qlm_derivs + use qlm_variables + use ricci2 + use tensor2 + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL, parameter :: one=1, two=2 + CCTK_REAL, parameter :: half=one/2, fourth=one/4, eighth=one/8 + CCTK_REAL, parameter :: sixteenth=one/16 + CCTK_REAL, parameter :: o128=one/128 + + CCTK_REAL :: qq(2,2), dtq, rsc + CCTK_REAL :: ee(3,2) + CCTK_REAL :: kk(3,3) + CCTK_REAL :: ll(0:3), nn(0:3), ss(0:3) + CCTK_COMPLEX :: psi2 + CCTK_REAL :: zz, dzz(2) + CCTK_REAL :: area, mass, spin + + CCTK_REAL :: delta_space(2) + + integer :: i, j +!!$ integer :: a, b, c, d + + if (veryverbose/=0) then + call CCTK_INFO ("Calculating multipole moments") + end if + + delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /) + + qlm_mp_m0(hn) = 0 + qlm_mp_m1(hn) = 0 + qlm_mp_m2(hn) = 0 + qlm_mp_m3(hn) = 0 + qlm_mp_m4(hn) = 0 + qlm_mp_m5(hn) = 0 + qlm_mp_m6(hn) = 0 + qlm_mp_m7(hn) = 0 + qlm_mp_m8(hn) = 0 + + qlm_mp_j0(hn) = 0 + qlm_mp_j1(hn) = 0 + qlm_mp_j2(hn) = 0 + qlm_mp_j3(hn) = 0 + qlm_mp_j4(hn) = 0 + qlm_mp_j5(hn) = 0 + qlm_mp_j6(hn) = 0 + qlm_mp_j7(hn) = 0 + qlm_mp_j8(hn) = 0 + + do j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn) + do i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn) + + ! 2-metric on the horizon + qq(1,1) = qlm_qtt(i,j,hn) + qq(1,2) = qlm_qtp(i,j,hn) + qq(2,2) = qlm_qpp(i,j,hn) + qq(2,1) = qq(1,2) + + call calc_2det (qq, dtq) + + 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) + + kk(1,1) = qlm_kxx(i,j) + kk(1,2) = qlm_kxy(i,j) + kk(1,3) = qlm_kxz(i,j) + kk(2,2) = qlm_kyy(i,j) + kk(2,3) = qlm_kyz(i,j) + kk(3,3) = qlm_kzz(i,j) + kk(2,1) = kk(1,2) + kk(3,1) = kk(1,3) + kk(3,2) = kk(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) + + ss = (ll - nn) / sqrt(two) + + rsc = qlm_rsc(i,j,hn) + psi2 = qlm_psi2(i,j,hn) + + zz = qlm_inv_z(i,j,hn) + + dzz(:) = deriv (qlm_inv_z(:,:,hn), i, j, delta_space) + + area = sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + + mass = fourth * rsc + + qlm_mp_m0(hn) = qlm_mp_m0(hn) + mass * p0(zz) * dzz(1) * area + qlm_mp_m1(hn) = qlm_mp_m1(hn) + mass * p1(zz) * dzz(1) * area + qlm_mp_m2(hn) = qlm_mp_m2(hn) + mass * p2(zz) * dzz(1) * area + qlm_mp_m3(hn) = qlm_mp_m3(hn) + mass * p3(zz) * dzz(1) * area + qlm_mp_m4(hn) = qlm_mp_m4(hn) + mass * p4(zz) * dzz(1) * area + qlm_mp_m5(hn) = qlm_mp_m5(hn) + mass * p5(zz) * dzz(1) * area + qlm_mp_m6(hn) = qlm_mp_m6(hn) + mass * p6(zz) * dzz(1) * area + qlm_mp_m7(hn) = qlm_mp_m7(hn) + mass * p7(zz) * dzz(1) * area + qlm_mp_m8(hn) = qlm_mp_m8(hn) + mass * p8(zz) * dzz(1) * area + + spin = - aimag(psi2) + + qlm_mp_j0(hn) = qlm_mp_j0(hn) + spin * p0(zz) * dzz(1) * area + qlm_mp_j1(hn) = qlm_mp_j1(hn) + spin * p1(zz) * dzz(1) * area + qlm_mp_j2(hn) = qlm_mp_j2(hn) + spin * p2(zz) * dzz(1) * area + qlm_mp_j3(hn) = qlm_mp_j3(hn) + spin * p3(zz) * dzz(1) * area + qlm_mp_j4(hn) = qlm_mp_j4(hn) + spin * p4(zz) * dzz(1) * area + qlm_mp_j5(hn) = qlm_mp_j5(hn) + spin * p5(zz) * dzz(1) * area + qlm_mp_j6(hn) = qlm_mp_j6(hn) + spin * p6(zz) * dzz(1) * area + qlm_mp_j7(hn) = qlm_mp_j7(hn) + spin * p7(zz) * dzz(1) * area + qlm_mp_j8(hn) = qlm_mp_j8(hn) + spin * p8(zz) * dzz(1) * area + +!!$ spin = 0 +!!$ do a=1,2 +!!$ do b=1,2 +!!$ do c=1,3 +!!$ do d=1,3 +!!$ spin = spin - half * area * epsilon2(a,b) * dzz(b) * ee(d,a) * kk(d,c) * ss(c) +!!$ end do +!!$ end do +!!$ end do +!!$ end do +!!$ +!!$ qlm_mp_j0(hn) = qlm_mp_j0(hn) + spin * dp0(zz) * area +!!$ qlm_mp_j1(hn) = qlm_mp_j1(hn) + spin * dp1(zz) * area +!!$ qlm_mp_j2(hn) = qlm_mp_j2(hn) + spin * dp2(zz) * area +!!$ qlm_mp_j3(hn) = qlm_mp_j3(hn) + spin * dp3(zz) * area +!!$ qlm_mp_j4(hn) = qlm_mp_j4(hn) + spin * dp4(zz) * area +!!$ qlm_mp_j5(hn) = qlm_mp_j5(hn) + spin * dp5(zz) * area +!!$ qlm_mp_j6(hn) = qlm_mp_j6(hn) + spin * dp6(zz) * area +!!$ qlm_mp_j7(hn) = qlm_mp_j7(hn) + spin * dp7(zz) * area +!!$ qlm_mp_j8(hn) = qlm_mp_j8(hn) + spin * dp8(zz) * area + + end do + end do + + ! Normalise + +!!$ ! This is the normalisation for I_n and L_n +!!$ qlm_mp_m0(hn) = qlm_mp_m0(hn) / sqrt(4*pi/ 1) +!!$ qlm_mp_m1(hn) = qlm_mp_m1(hn) / sqrt(4*pi/ 3) +!!$ qlm_mp_m2(hn) = qlm_mp_m2(hn) / sqrt(4*pi/ 5) +!!$ qlm_mp_m3(hn) = qlm_mp_m3(hn) / sqrt(4*pi/ 7) +!!$ qlm_mp_m4(hn) = qlm_mp_m4(hn) / sqrt(4*pi/ 9) +!!$ qlm_mp_m5(hn) = qlm_mp_m5(hn) / sqrt(4*pi/11) +!!$ qlm_mp_m6(hn) = qlm_mp_m6(hn) / sqrt(4*pi/13) +!!$ qlm_mp_m7(hn) = qlm_mp_m7(hn) / sqrt(4*pi/15) +!!$ qlm_mp_m8(hn) = qlm_mp_m8(hn) / sqrt(4*pi/17) +!!$ +!!$ qlm_mp_j0(hn) = qlm_mp_j0(hn) / sqrt(4*pi/ 1) +!!$ qlm_mp_j1(hn) = qlm_mp_j1(hn) / sqrt(4*pi/ 3) +!!$ qlm_mp_j2(hn) = qlm_mp_j2(hn) / sqrt(4*pi/ 5) +!!$ qlm_mp_j3(hn) = qlm_mp_j3(hn) / sqrt(4*pi/ 7) +!!$ qlm_mp_j4(hn) = qlm_mp_j4(hn) / sqrt(4*pi/ 9) +!!$ qlm_mp_j5(hn) = qlm_mp_j5(hn) / sqrt(4*pi/11) +!!$ qlm_mp_j6(hn) = qlm_mp_j6(hn) / sqrt(4*pi/13) +!!$ qlm_mp_j7(hn) = qlm_mp_j7(hn) / sqrt(4*pi/15) +!!$ qlm_mp_j8(hn) = qlm_mp_j8(hn) / sqrt(4*pi/17) + + ! This is the normalisation for M_n and J_n + qlm_mp_m0(hn) = qlm_mp_m0(hn) * qlm_mass(hn) * qlm_radius(hn)**0 / (2*pi) + qlm_mp_m1(hn) = qlm_mp_m1(hn) * qlm_mass(hn) * qlm_radius(hn)**1 / (2*pi) + qlm_mp_m2(hn) = qlm_mp_m2(hn) * qlm_mass(hn) * qlm_radius(hn)**2 / (2*pi) + qlm_mp_m3(hn) = qlm_mp_m3(hn) * qlm_mass(hn) * qlm_radius(hn)**3 / (2*pi) + qlm_mp_m4(hn) = qlm_mp_m4(hn) * qlm_mass(hn) * qlm_radius(hn)**4 / (2*pi) + qlm_mp_m5(hn) = qlm_mp_m5(hn) * qlm_mass(hn) * qlm_radius(hn)**5 / (2*pi) + qlm_mp_m6(hn) = qlm_mp_m6(hn) * qlm_mass(hn) * qlm_radius(hn)**6 / (2*pi) + qlm_mp_m7(hn) = qlm_mp_m7(hn) * qlm_mass(hn) * qlm_radius(hn)**7 / (2*pi) + qlm_mp_m8(hn) = qlm_mp_m8(hn) * qlm_mass(hn) * qlm_radius(hn)**8 / (2*pi) + + qlm_mp_j0(hn) = qlm_mp_j0(hn) * qlm_radius(hn)**1 / (4*pi) + qlm_mp_j1(hn) = qlm_mp_j1(hn) * qlm_radius(hn)**2 / (4*pi) + qlm_mp_j2(hn) = qlm_mp_j2(hn) * qlm_radius(hn)**3 / (4*pi) + qlm_mp_j3(hn) = qlm_mp_j3(hn) * qlm_radius(hn)**4 / (4*pi) + qlm_mp_j4(hn) = qlm_mp_j4(hn) * qlm_radius(hn)**5 / (4*pi) + qlm_mp_j5(hn) = qlm_mp_j5(hn) * qlm_radius(hn)**6 / (4*pi) + qlm_mp_j6(hn) = qlm_mp_j6(hn) * qlm_radius(hn)**7 / (4*pi) + qlm_mp_j7(hn) = qlm_mp_j7(hn) * qlm_radius(hn)**8 / (4*pi) + qlm_mp_j8(hn) = qlm_mp_j8(hn) * qlm_radius(hn)**9 / (4*pi) + +contains + + ! Legendre polynomials + + function p0 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p0 + p0 = 1 + end function p0 + + function p1 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p1 + p1 = z + end function p1 + + function p2 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p2 + p2 = 3*half * z**2 - half + end function p2 + + function p3 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p3 + p3 = 5*half * z**3 - 3*half * z + end function p3 + + function p4 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p4 + p4 = 35*eighth * z**4 - 15*fourth * z**2 + 3*eighth + end function p4 + + function p5 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p5 + p5 = 63*eighth * z**5 - 35*fourth * z**3 + 15*eighth * z + end function p5 + + function p6 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p6 + p6 = 231*sixteenth * z**6 - 315*sixteenth * z**4 + 105*sixteenth * z**2 & + - 5*sixteenth + end function p6 + + function p7 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p7 + p7 = 429*sixteenth * z**7 - 693*sixteenth * z**5 + 315*sixteenth * z**3 & + - 35*sixteenth * z + end function p7 + + function p8 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: p8 + p8 = 6435*o128 * z**8 - 12012*o128 * z**6 + 6930*o128 * z**4 & + - 1260*o128 * z**2 + 35*o128 + end function p8 + + ! Derivatives of the Legendre polynomials + + function dp0 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp0 + dp0 = 0 + end function dp0 + + function dp1 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp1 + dp1 = 1 + end function dp1 + + function dp2 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp2 + dp2 = 3 * z + end function dp2 + + function dp3 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp3 + dp3 = 15*half * z**2 - 3*half + end function dp3 + + function dp4 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp4 + dp4 = 35*half * z**3 - 15*half*z + end function dp4 + + function dp5 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp5 + dp5 = 315*eighth * z**4 - 105*fourth * z**2 - 15*eighth + end function dp5 + + function dp6 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp6 + dp6 = 693*eighth * z**5 - 315*fourth * z**3 + 105*eighth * z + end function dp6 + + function dp7 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp7 + dp7 = 3003*sixteenth * z**6 - 3465*sixteenth * z**4 + 945*sixteenth * z**2 & + - 35*sixteenth + end function dp7 + + function dp8 (z) + CCTK_REAL, intent(in) :: z + CCTK_REAL :: dp8 + dp8 = 51480*o128 * z**7 - 72072*o128 * z**5 + 27720*o128 * z**3 & + - 2520*o128 * z + end function dp8 + +end subroutine qlm_multipoles 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 diff --git a/src/qlm_paramcheck.F90 b/src/qlm_paramcheck.F90 new file mode 100644 index 0000000..b8a8b8e --- /dev/null +++ b/src/qlm_paramcheck.F90 @@ -0,0 +1,76 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_paramcheck (CCTK_ARGUMENTS) + use cctk + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + integer :: sn + integer :: eq_theta, new_ntheta, int_nphi, new_nphi + character :: msg*1000 + + if (veryverbose/=0) then + call CCTK_INFO ("Checking parameters") + end if + + do hn = 1, num_surfaces + + if (surface_index(hn) == -1) then + ! no surface selected; everything is fine + goto 9999 + end if + + if (surface_index(hn) < 0 .or. surface_index(hn) >= nsurfaces) then + write (msg, '("Illegal surface index specified for surface ",i4," (index is ",i4,", must be less than ",i4,")")') hn-1, surface_index(hn), nsurfaces + call CCTK_PARAMWARN (msg) + goto 9999 + end if + + sn = surface_index(hn) + 1 + + ! Import surface description + qlm_nghoststheta(hn) = nghoststheta(sn) + qlm_nghostsphi (hn) = nghostsphi(sn) + qlm_ntheta (hn) = ntheta(sn) + qlm_nphi (hn) = nphi(sn) + + ! Symmetries + if (symmetric_x(sn) /= 0 .or. & + symmetric_y(sn) /= 0 .or. & + symmetric_z(sn) /= 0) then + call CCTK_WARN (CCTK_WARN_ABORT, "SphericalSurface symmetries are not supported") + end if + + if (auto_res(sn) /= 1 .and. (qlm_ntheta(hn) > maxntheta .or. qlm_nphi(hn) > maxnphi)) then + write (msg, '("Surface ",i4," is too large: shape is (",2i6,"), maximum is (",2i6,")")') hn-1, qlm_ntheta(hn), qlm_nphi(hn), maxntheta, maxnphi + call CCTK_PARAMWARN (msg) + end if + + if (qlm_nghoststheta(hn)<1 .or. qlm_nghostsphi(hn)<1) then + write (msg, '("Not enough ghost zones for surface ",i4,": nghosts=",2i4,", minimum is ",2i4)') hn-1, qlm_nghoststheta(hn), qlm_nghostsphi(hn), 1, 1 + call CCTK_PARAMWARN (msg) + end if + + if (auto_res(sn) /= 1 .and. (mod(int(qlm_ntheta(hn) - 2*qlm_nghoststheta(hn)),2) /= 1)) then + ! We need a grid point on the equator + write (msg, '("The number of interior grid points in the theta direction of surface ",i4," must be odd after the symmetries have been removed, but it is ",i6)') hn-1, qlm_ntheta(hn) - 2*qlm_nghoststheta(hn) + call CCTK_PARAMWARN (msg) + end if + + if (auto_res(sn) /= 1 .and. (mod(int(qlm_nphi(hn) - 2*qlm_nghostsphi(hn)),4) /= 0)) then + ! We need grid points on the four major meridians + write (msg, '("The number of interior grid points in the phi direction of surface ",i4," must be a multiple of four after the symmetries have been removed, but it is ",i6)') hn-1, qlm_nphi(hn) - 2*qlm_nghostsphi(hn) + call CCTK_PARAMWARN (msg) + end if + +9999 continue + end do + +end subroutine qlm_paramcheck diff --git a/src/qlm_set_coordinates.F90 b/src/qlm_set_coordinates.F90 new file mode 100644 index 0000000..309a361 --- /dev/null +++ b/src/qlm_set_coordinates.F90 @@ -0,0 +1,46 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_set_coordinates (CCTK_ARGUMENTS, hn) + use cctk + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + CCTK_REAL :: theta, phi + CCTK_REAL :: sin_theta, cos_theta, sin_phi, cos_phi + + integer :: i, j + + ! Calculate the coordinates + do j = 1, qlm_nphi(hn) + do i = 1, qlm_ntheta(hn) + theta = qlm_origin_theta(hn) + (i-1)*qlm_delta_theta(hn) + phi = qlm_origin_phi(hn) + (j-1)*qlm_delta_phi(hn) + + sin_theta = sin(theta) + cos_theta = cos(theta) + sin_phi = sin(phi) + cos_phi = cos(phi) + + qlm_x(i,j,hn) = qlm_origin_x(hn) + qlm_shape(i,j,hn) * sin_theta * cos_phi + qlm_y(i,j,hn) = qlm_origin_y(hn) + qlm_shape(i,j,hn) * sin_theta * sin_phi + qlm_z(i,j,hn) = qlm_origin_z(hn) + qlm_shape(i,j,hn) * cos_theta + + qlm_x_p(i,j,hn) = qlm_origin_x_p(hn) + qlm_shape_p(i,j,hn) * sin_theta * cos_phi + qlm_y_p(i,j,hn) = qlm_origin_y_p(hn) + qlm_shape_p(i,j,hn) * sin_theta * sin_phi + qlm_z_p(i,j,hn) = qlm_origin_z_p(hn) + qlm_shape_p(i,j,hn) * cos_theta + + qlm_x_p_p(i,j,hn) = qlm_origin_x_p_p(hn) + qlm_shape_p_p(i,j,hn) * sin_theta * cos_phi + qlm_y_p_p(i,j,hn) = qlm_origin_y_p_p(hn) + qlm_shape_p_p(i,j,hn) * sin_theta * sin_phi + qlm_z_p_p(i,j,hn) = qlm_origin_z_p_p(hn) + qlm_shape_p_p(i,j,hn) * cos_theta + end do + end do + +end subroutine qlm_set_coordinates diff --git a/src/qlm_tetrad.F90 b/src/qlm_tetrad.F90 new file mode 100644 index 0000000..85837dc --- /dev/null +++ b/src/qlm_tetrad.F90 @@ -0,0 +1,324 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "cctk_Parameters.h" + + + +subroutine qlm_calc_tetrad (CCTK_ARGUMENTS, hn) + use adm_metric_simple + use cctk + use classify + use matinv + use pointwise2 + use qlm_boundary + use qlm_derivs + use qlm_gram_schmidt + 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 :: one=1, two=2 + CCTK_REAL, parameter :: gg4(0:3,0:3,0:3) = 0 + CCTK_REAL :: gg(3,3), dgg(3,3,3), gg_dot(3,3) + CCTK_REAL :: kk(3,3) + CCTK_REAL :: g4(0:3,0:3), gu4(0:3,0:3), dg4(0:3,0:3,0:3) + CCTK_REAL :: gamma4(0:3,0:3,0:3) + CCTK_REAL :: ee(0:3,0:3), ee_p(0:3,1:3), ee_p_p(0:3,1:3) + CCTK_REAL :: dee_spher(1:3,1:3,1:3), ee_inv(1:3,1:3) + CCTK_REAL :: dee(0:3,0:3,0:3), gee(0:3,0:3,0:3) + CCTK_REAL :: m1(0:3), m2(0:3) ! temporary variables to calculate mm + CCTK_REAL :: ss(0:3) ! spacelike outward normal to horizon + CCTK_REAL :: tt(0:3) ! timelike unit normal to hypersurface + CCTK_REAL :: ll(0:3) ! future null vector on the horizon + CCTK_REAL :: nn(0:3) ! future inward null vector + CCTK_COMPLEX :: mm(0:3) ! vector on horizon within the hypersurface + CCTK_REAL :: gtt(0:3,0:3), gss(0:3,0:3) + CCTK_REAL :: gm1(0:3,0:3), gm2(0:3,0:3) + CCTK_REAL :: gll(0:3,0:3), gnn(0:3,0:3) + CCTK_COMPLEX :: gmm(0:3,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) + + CCTK_REAL :: count, accuracy + + integer :: lsh(2) + integer :: i, j + integer :: a, b, c, d + CCTK_REAL :: theta, phi + + logical :: lerr + + character*2, parameter :: crlf = achar(13) // achar(10) + character :: msg*1000 + + if (veryverbose/=0) then + call CCTK_INFO ("Setting tetrad") + end if + + lsh(:) = (/ qlm_ntheta(hn), qlm_nphi(hn) /) + delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /) + + count = 0 + accuracy = 0 + + ! 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 variables 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(2,1,1) = dgg(1,2,1) + dgg(3,1,1) = dgg(1,3,1) + dgg(3,2,1) = dgg(2,3,1) + 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(2,1,2) = dgg(1,2,2) + dgg(3,1,2) = dgg(1,3,2) + dgg(3,2,2) = dgg(2,3,2) + 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,3) = dgg(1,2,3) + dgg(3,1,3) = dgg(1,3,3) + dgg(3,2,3) = dgg(2,3,3) + + kk(1,1) = qlm_kxx(i,j) + kk(1,2) = qlm_kxy(i,j) + kk(1,3) = qlm_kxz(i,j) + kk(2,2) = qlm_kyy(i,j) + kk(2,3) = qlm_kyz(i,j) + kk(3,3) = qlm_kzz(i,j) + kk(2,1) = kk(1,2) + kk(3,1) = kk(1,3) + kk(3,2) = kk(2,3) + + + + ! Calculate 4-metric + call calc_3metricdot_simple (kk, gg_dot) + call calc_4metricderivs_simple (gg,dgg,gg_dot, g4,dg4) + call calc_4inv (g4, gu4) + call calc_4connections (gu4,dg4, gamma4) + + + + ! The following must be consistent with qlm_tetrad.F90 + + ee = TAT_nan() + dee_spher = TAT_nan() + dee = TAT_nan() + + ee(0,:) = 0 + dee(0,:,:) = 0 + + + + ee(1,0) = 0 + ee(1,1) = qlm_x(i,j,hn) - qlm_origin_x(hn) + ee(1,2) = qlm_y(i,j,hn) - qlm_origin_y(hn) + ee(1,3) = qlm_z(i,j,hn) - qlm_origin_z(hn) + + ee_p(1,1) = qlm_x_p(i,j,hn) - qlm_origin_x_p(hn) + ee_p(1,2) = qlm_y_p(i,j,hn) - qlm_origin_y_p(hn) + ee_p(1,3) = qlm_z_p(i,j,hn) - qlm_origin_z_p(hn) + + ee_p_p(1,1) = qlm_x_p_p(i,j,hn) - qlm_origin_x_p_p(hn) + ee_p_p(1,2) = qlm_y_p_p(i,j,hn) - qlm_origin_y_p_p(hn) + ee_p_p(1,3) = qlm_z_p_p(i,j,hn) - qlm_origin_z_p_p(hn) + + dee(1,0,:) = 0 + dee(1,1:3,0) = timederiv (ee(1,1:3), ee_p(1,1:3), ee_p_p(1,1:3), t0,t1,t2, ce0,ce1,ce2) + dee_spher(1,:,1) = 0 ! this is a choice + dee_spher(1,1,2:3) = deriv (qlm_x(:,:,hn), i,j, delta_space) + dee_spher(1,2,2:3) = deriv (qlm_y(:,:,hn), i,j, delta_space) + dee_spher(1,3,2:3) = deriv (qlm_z(:,:,hn), i,j, delta_space) + + + + ee(2:3,0) = 0 + ee(2:3,1) = deriv (qlm_x(:,:,hn), i,j, delta_space) + ee(2:3,2) = deriv (qlm_y(:,:,hn), i,j, delta_space) + ee(2:3,3) = deriv (qlm_z(:,:,hn), i,j, delta_space) + + ee_p(2:3,1) = deriv (qlm_x_p(:,:,hn), i,j, delta_space) + ee_p(2:3,2) = deriv (qlm_y_p(:,:,hn), i,j, delta_space) + ee_p(2:3,3) = deriv (qlm_z_p(:,:,hn), i,j, delta_space) + + ee_p_p(2:3,1) = deriv (qlm_x_p_p(:,:,hn), i,j, delta_space) + ee_p_p(2:3,2) = deriv (qlm_y_p_p(:,:,hn), i,j, delta_space) + ee_p_p(2:3,3) = deriv (qlm_z_p_p(:,:,hn), i,j, delta_space) + + dee(2:3,0,:) = 0 + dee(2:3,1:3,0) = timederiv (ee(2:3,1:3), ee_p(2:3,1:3), ee_p_p(2:3,1:3), t0,t1,t2, ce0,ce1,ce2) + dee_spher(2:3,:,1) = 0 ! this is a choice + dee_spher(2:3,1,2:3) = deriv2 (qlm_x(:,:,hn), i,j, delta_space) + dee_spher(2:3,2,2:3) = deriv2 (qlm_y(:,:,hn), i,j, delta_space) + dee_spher(2:3,3,2:3) = deriv2 (qlm_z(:,:,hn), i,j, delta_space) + + ! ee_a^i + ! dee_spher_a^i,b + ! dee_a^i,j = ee_j^b dee_spher_a^i,b + call calc_inv3 (ee(1:3,1:3), ee_inv, lerr) + if (lerr) then + call CCTK_WARN (3, "Could not invert matrix") + end if + do a=1,3 + do b=1,3 + do c=1,3 + dee(a,b,c) = 0 + do d=1,3 + dee(a,b,c) = dee(a,b,c) + ee_inv(c,d) * dee_spher(a,b,d) + end do + end do + end do + end do + + do a=0,3 + do b=0,3 + gee(:,a,b) = dee(:,a,b) + do c=0,3 + gee(:,a,b) = gee(:,a,b) + ee(:,c) * gamma4(a,c,b) + end do + end do + end do + + + + ! tt + tt(:) = ee(0,:) + gtt(:,:) = gee(0,:,:) + + ! m1 = ep + m1(:) = ee(3,:) + gm1(:,:) = gee(3,:,:) + call gram_schmidt_normalise (g4,gg4, m1,gm1, one) + + ! m2 = et + m2(:) = ee(2,:) + gm2(:,:) = gee(2,:,:) + call gram_schmidt_project (g4,gg4, m1,gm1, one, m2,gm2) + call gram_schmidt_normalise (g4,gg4, m2,gm2, one) + + ! ss = er + ss(:) = ee(1,:) + gss(:,:) = gee(1,:,:) + call gram_schmidt_project (g4,gg4, m1,gm1, one, ss,gss) + call gram_schmidt_project (g4,gg4, m2,gm2, one, ss,gss) + call gram_schmidt_normalise (g4,gg4, ss,gss, one) + + + + ! ll = (tt + ss) / sqrt(two) + ll = (tt + ss) / sqrt(two) + gll = (gtt + gss) / sqrt(two) + + ! nn = (tt - ss) / sqrt(two) + nn = (tt - ss) / sqrt(two) + gnn = (gtt - gss) / sqrt(two) + + ! mm = cmplx(m1, m2, kind(mm)) / sqrt(two) + mm = cmplx(m1, m2, kind(mm)) / sqrt(two) + gmm = cmplx(gm1, gm2, kind(gmm)) / sqrt(two) + + + + ! Store the stuff into the arrays + do a=0,3 + do b=0,3 + nabla_ll(a,b) = 0 + nabla_nn(a,b) = 0 + nabla_mm(a,b) = 0 + do c=0,3 + nabla_ll(a,b) = nabla_ll(a,b) + g4(a,c) * gll(c,b) + nabla_nn(a,b) = nabla_nn(a,b) + g4(a,c) * gnn(c,b) + nabla_mm(a,b) = nabla_mm(a,b) + g4(a,c) * gmm(c,b) + end do + qlm_tetrad_derivs(i,j)%nabla_ll(a,b) = nabla_ll(a,b) + qlm_tetrad_derivs(i,j)%nabla_nn(a,b) = nabla_nn(a,b) + qlm_tetrad_derivs(i,j)%nabla_mm(a,b) = nabla_mm(a,b) + end do + end do + + qlm_l0(i,j,hn) = ll(0) + qlm_l1(i,j,hn) = ll(1) + qlm_l2(i,j,hn) = ll(2) + qlm_l3(i,j,hn) = ll(3) + + qlm_n0(i,j,hn) = nn(0) + qlm_n1(i,j,hn) = nn(1) + qlm_n2(i,j,hn) = nn(2) + qlm_n3(i,j,hn) = nn(3) + + qlm_m0(i,j,hn) = mm(0) + qlm_m1(i,j,hn) = mm(1) + qlm_m2(i,j,hn) = mm(2) + qlm_m3(i,j,hn) = mm(3) + + end do + end do + + if (count > 0) then + accuracy = sqrt(accuracy / count) + end if + + if (veryverbose/=0) then + write (msg, '("Tetrad accuracy L2 norm: ",g12.4)') accuracy + call CCTK_INFO (msg) + end if + +!!$ if (accuracy > 1.0d-12) then + if (accuracy > 1.0d-8) then + call CCTK_WARN (1, "Tetrad is not accurate") + end if + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_l0(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_l1(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_l2(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_l3(:,:,hn), +1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_n0(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_n1(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_n2(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_n3(:,:,hn), +1) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_m0(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_m1(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_m2(:,:,hn), +1) + call set_boundary (CCTK_PASS_FTOF, hn, qlm_m3(:,:,hn), +1) + +end subroutine qlm_calc_tetrad 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 diff --git a/src/qlm_variables.F90 b/src/qlm_variables.F90 new file mode 100644 index 0000000..d21842c --- /dev/null +++ b/src/qlm_variables.F90 @@ -0,0 +1,313 @@ +#include "cctk.h" + +module qlm_variables + use classify + implicit none + public + save + CCTK_REAL, allocatable, dimension(:,:) :: & + qlm_gxx, qlm_gxy, qlm_gxz, qlm_gyy, qlm_gyz, qlm_gzz, & + qlm_dgxxx, qlm_dgxyx, qlm_dgxzx, qlm_dgyyx, qlm_dgyzx, qlm_dgzzx, & + qlm_dgxxy, qlm_dgxyy, qlm_dgxzy, qlm_dgyyy, qlm_dgyzy, qlm_dgzzy, & + qlm_dgxxz, qlm_dgxyz, qlm_dgxzz, qlm_dgyyz, qlm_dgyzz, qlm_dgzzz, & + qlm_ddgxxxx, qlm_ddgxyxx, qlm_ddgxzxx, qlm_ddgyyxx, qlm_ddgyzxx, qlm_ddgzzxx, & + qlm_ddgxxxy, qlm_ddgxyxy, qlm_ddgxzxy, qlm_ddgyyxy, qlm_ddgyzxy, qlm_ddgzzxy, & + qlm_ddgxxxz, qlm_ddgxyxz, qlm_ddgxzxz, qlm_ddgyyxz, qlm_ddgyzxz, qlm_ddgzzxz, & + qlm_ddgxxyy, qlm_ddgxyyy, qlm_ddgxzyy, qlm_ddgyyyy, qlm_ddgyzyy, qlm_ddgzzyy, & + qlm_ddgxxyz, qlm_ddgxyyz, qlm_ddgxzyz, qlm_ddgyyyz, qlm_ddgyzyz, qlm_ddgzzyz, & + qlm_ddgxxzz, qlm_ddgxyzz, qlm_ddgxzzz, qlm_ddgyyzz, qlm_ddgyzzz, qlm_ddgzzzz, & + qlm_kxx, qlm_kxy, qlm_kxz, qlm_kyy, qlm_kyz, qlm_kzz, & + qlm_dkxxx, qlm_dkxyx, qlm_dkxzx, qlm_dkyyx, qlm_dkyzx, qlm_dkzzx, & + qlm_dkxxy, qlm_dkxyy, qlm_dkxzy, qlm_dkyyy, qlm_dkyzy, qlm_dkzzy, & + qlm_dkxxz, qlm_dkxyz, qlm_dkxzz, qlm_dkyyz, qlm_dkyzz, qlm_dkzzz, & + qlm_alpha, & + qlm_betax, qlm_betay, qlm_betaz + + type tetrad_derivs + ! nabla_ll(a,b) = D_b l_a + CCTK_REAL :: nabla_ll(0:3,0:3) + CCTK_REAL :: nabla_nn(0:3,0:3) + CCTK_COMPLEX :: nabla_mm(0:3,0:3) + end type tetrad_derivs + type(tetrad_derivs), allocatable, dimension(:,:) :: qlm_tetrad_derivs + +contains + + subroutine allocate_variables(ntheta, nphi) + integer, intent(in) :: ntheta, nphi + + allocate(qlm_gxx(ntheta,nphi)) + allocate(qlm_gxy(ntheta,nphi)) + allocate(qlm_gxz(ntheta,nphi)) + allocate(qlm_gyy(ntheta,nphi)) + allocate(qlm_gyz(ntheta,nphi)) + allocate(qlm_gzz(ntheta,nphi)) + allocate(qlm_dgxxx(ntheta,nphi)) + allocate(qlm_dgxyx(ntheta,nphi)) + allocate(qlm_dgxzx(ntheta,nphi)) + allocate(qlm_dgyyx(ntheta,nphi)) + allocate(qlm_dgyzx(ntheta,nphi)) + allocate(qlm_dgzzx(ntheta,nphi)) + allocate(qlm_dgxxy(ntheta,nphi)) + allocate(qlm_dgxyy(ntheta,nphi)) + allocate(qlm_dgxzy(ntheta,nphi)) + allocate(qlm_dgyyy(ntheta,nphi)) + allocate(qlm_dgyzy(ntheta,nphi)) + allocate(qlm_dgzzy(ntheta,nphi)) + allocate(qlm_dgxxz(ntheta,nphi)) + allocate(qlm_dgxyz(ntheta,nphi)) + allocate(qlm_dgxzz(ntheta,nphi)) + allocate(qlm_dgyyz(ntheta,nphi)) + allocate(qlm_dgyzz(ntheta,nphi)) + allocate(qlm_dgzzz(ntheta,nphi)) + allocate(qlm_ddgxxxx(ntheta,nphi)) + allocate(qlm_ddgxyxx(ntheta,nphi)) + allocate(qlm_ddgxzxx(ntheta,nphi)) + allocate(qlm_ddgyyxx(ntheta,nphi)) + allocate(qlm_ddgyzxx(ntheta,nphi)) + allocate(qlm_ddgzzxx(ntheta,nphi)) + allocate(qlm_ddgxxxy(ntheta,nphi)) + allocate(qlm_ddgxyxy(ntheta,nphi)) + allocate(qlm_ddgxzxy(ntheta,nphi)) + allocate(qlm_ddgyyxy(ntheta,nphi)) + allocate(qlm_ddgyzxy(ntheta,nphi)) + allocate(qlm_ddgzzxy(ntheta,nphi)) + allocate(qlm_ddgxxxz(ntheta,nphi)) + allocate(qlm_ddgxyxz(ntheta,nphi)) + allocate(qlm_ddgxzxz(ntheta,nphi)) + allocate(qlm_ddgyyxz(ntheta,nphi)) + allocate(qlm_ddgyzxz(ntheta,nphi)) + allocate(qlm_ddgzzxz(ntheta,nphi)) + allocate(qlm_ddgxxyy(ntheta,nphi)) + allocate(qlm_ddgxyyy(ntheta,nphi)) + allocate(qlm_ddgxzyy(ntheta,nphi)) + allocate(qlm_ddgyyyy(ntheta,nphi)) + allocate(qlm_ddgyzyy(ntheta,nphi)) + allocate(qlm_ddgzzyy(ntheta,nphi)) + allocate(qlm_ddgxxyz(ntheta,nphi)) + allocate(qlm_ddgxyyz(ntheta,nphi)) + allocate(qlm_ddgxzyz(ntheta,nphi)) + allocate(qlm_ddgyyyz(ntheta,nphi)) + allocate(qlm_ddgyzyz(ntheta,nphi)) + allocate(qlm_ddgzzyz(ntheta,nphi)) + allocate(qlm_ddgxxzz(ntheta,nphi)) + allocate(qlm_ddgxyzz(ntheta,nphi)) + allocate(qlm_ddgxzzz(ntheta,nphi)) + allocate(qlm_ddgyyzz(ntheta,nphi)) + allocate(qlm_ddgyzzz(ntheta,nphi)) + allocate(qlm_ddgzzzz(ntheta,nphi)) + allocate(qlm_kxx(ntheta,nphi)) + allocate(qlm_kxy(ntheta,nphi)) + allocate(qlm_kxz(ntheta,nphi)) + allocate(qlm_kyy(ntheta,nphi)) + allocate(qlm_kyz(ntheta,nphi)) + allocate(qlm_kzz(ntheta,nphi)) + allocate(qlm_dkxxx(ntheta,nphi)) + allocate(qlm_dkxyx(ntheta,nphi)) + allocate(qlm_dkxzx(ntheta,nphi)) + allocate(qlm_dkyyx(ntheta,nphi)) + allocate(qlm_dkyzx(ntheta,nphi)) + allocate(qlm_dkzzx(ntheta,nphi)) + allocate(qlm_dkxxy(ntheta,nphi)) + allocate(qlm_dkxyy(ntheta,nphi)) + allocate(qlm_dkxzy(ntheta,nphi)) + allocate(qlm_dkyyy(ntheta,nphi)) + allocate(qlm_dkyzy(ntheta,nphi)) + allocate(qlm_dkzzy(ntheta,nphi)) + allocate(qlm_dkxxz(ntheta,nphi)) + allocate(qlm_dkxyz(ntheta,nphi)) + allocate(qlm_dkxzz(ntheta,nphi)) + allocate(qlm_dkyyz(ntheta,nphi)) + allocate(qlm_dkyzz(ntheta,nphi)) + allocate(qlm_dkzzz(ntheta,nphi)) + allocate(qlm_alpha(ntheta,nphi)) + allocate(qlm_betax(ntheta,nphi)) + allocate(qlm_betay(ntheta,nphi)) + allocate(qlm_betaz(ntheta,nphi)) + + allocate(qlm_tetrad_derivs(ntheta,nphi)) + + qlm_gxx = TAT_nan() + qlm_gxy = TAT_nan() + qlm_gxz = TAT_nan() + qlm_gyy = TAT_nan() + qlm_gyz = TAT_nan() + qlm_gzz = TAT_nan() + qlm_dgxxx = TAT_nan() + qlm_dgxyx = TAT_nan() + qlm_dgxzx = TAT_nan() + qlm_dgyyx = TAT_nan() + qlm_dgyzx = TAT_nan() + qlm_dgzzx = TAT_nan() + qlm_dgxxy = TAT_nan() + qlm_dgxyy = TAT_nan() + qlm_dgxzy = TAT_nan() + qlm_dgyyy = TAT_nan() + qlm_dgyzy = TAT_nan() + qlm_dgzzy = TAT_nan() + qlm_dgxxz = TAT_nan() + qlm_dgxyz = TAT_nan() + qlm_dgxzz = TAT_nan() + qlm_dgyyz = TAT_nan() + qlm_dgyzz = TAT_nan() + qlm_dgzzz = TAT_nan() + qlm_ddgxxxx = TAT_nan() + qlm_ddgxyxx = TAT_nan() + qlm_ddgxzxx = TAT_nan() + qlm_ddgyyxx = TAT_nan() + qlm_ddgyzxx = TAT_nan() + qlm_ddgzzxx = TAT_nan() + qlm_ddgxxxy = TAT_nan() + qlm_ddgxyxy = TAT_nan() + qlm_ddgxzxy = TAT_nan() + qlm_ddgyyxy = TAT_nan() + qlm_ddgyzxy = TAT_nan() + qlm_ddgzzxy = TAT_nan() + qlm_ddgxxxz = TAT_nan() + qlm_ddgxyxz = TAT_nan() + qlm_ddgxzxz = TAT_nan() + qlm_ddgyyxz = TAT_nan() + qlm_ddgyzxz = TAT_nan() + qlm_ddgzzxz = TAT_nan() + qlm_ddgxxyy = TAT_nan() + qlm_ddgxyyy = TAT_nan() + qlm_ddgxzyy = TAT_nan() + qlm_ddgyyyy = TAT_nan() + qlm_ddgyzyy = TAT_nan() + qlm_ddgzzyy = TAT_nan() + qlm_ddgxxyz = TAT_nan() + qlm_ddgxyyz = TAT_nan() + qlm_ddgxzyz = TAT_nan() + qlm_ddgyyyz = TAT_nan() + qlm_ddgyzyz = TAT_nan() + qlm_ddgzzyz = TAT_nan() + qlm_ddgxxzz = TAT_nan() + qlm_ddgxyzz = TAT_nan() + qlm_ddgxzzz = TAT_nan() + qlm_ddgyyzz = TAT_nan() + qlm_ddgyzzz = TAT_nan() + qlm_ddgzzzz = TAT_nan() + qlm_kxx = TAT_nan() + qlm_kxy = TAT_nan() + qlm_kxz = TAT_nan() + qlm_kyy = TAT_nan() + qlm_kyz = TAT_nan() + qlm_kzz = TAT_nan() + qlm_dkxxx = TAT_nan() + qlm_dkxyx = TAT_nan() + qlm_dkxzx = TAT_nan() + qlm_dkyyx = TAT_nan() + qlm_dkyzx = TAT_nan() + qlm_dkzzx = TAT_nan() + qlm_dkxxy = TAT_nan() + qlm_dkxyy = TAT_nan() + qlm_dkxzy = TAT_nan() + qlm_dkyyy = TAT_nan() + qlm_dkyzy = TAT_nan() + qlm_dkzzy = TAT_nan() + qlm_dkxxz = TAT_nan() + qlm_dkxyz = TAT_nan() + qlm_dkxzz = TAT_nan() + qlm_dkyyz = TAT_nan() + qlm_dkyzz = TAT_nan() + qlm_dkzzz = TAT_nan() + qlm_alpha = TAT_nan() + qlm_betax = TAT_nan() + qlm_betay = TAT_nan() + qlm_betaz = TAT_nan() + end subroutine allocate_variables + + subroutine deallocate_variables + deallocate(qlm_gxx) + deallocate(qlm_gxy) + deallocate(qlm_gxz) + deallocate(qlm_gyy) + deallocate(qlm_gyz) + deallocate(qlm_gzz) + deallocate(qlm_dgxxx) + deallocate(qlm_dgxyx) + deallocate(qlm_dgxzx) + deallocate(qlm_dgyyx) + deallocate(qlm_dgyzx) + deallocate(qlm_dgzzx) + deallocate(qlm_dgxxy) + deallocate(qlm_dgxyy) + deallocate(qlm_dgxzy) + deallocate(qlm_dgyyy) + deallocate(qlm_dgyzy) + deallocate(qlm_dgzzy) + deallocate(qlm_dgxxz) + deallocate(qlm_dgxyz) + deallocate(qlm_dgxzz) + deallocate(qlm_dgyyz) + deallocate(qlm_dgyzz) + deallocate(qlm_dgzzz) + deallocate(qlm_ddgxxxx) + deallocate(qlm_ddgxyxx) + deallocate(qlm_ddgxzxx) + deallocate(qlm_ddgyyxx) + deallocate(qlm_ddgyzxx) + deallocate(qlm_ddgzzxx) + deallocate(qlm_ddgxxxy) + deallocate(qlm_ddgxyxy) + deallocate(qlm_ddgxzxy) + deallocate(qlm_ddgyyxy) + deallocate(qlm_ddgyzxy) + deallocate(qlm_ddgzzxy) + deallocate(qlm_ddgxxxz) + deallocate(qlm_ddgxyxz) + deallocate(qlm_ddgxzxz) + deallocate(qlm_ddgyyxz) + deallocate(qlm_ddgyzxz) + deallocate(qlm_ddgzzxz) + deallocate(qlm_ddgxxyy) + deallocate(qlm_ddgxyyy) + deallocate(qlm_ddgxzyy) + deallocate(qlm_ddgyyyy) + deallocate(qlm_ddgyzyy) + deallocate(qlm_ddgzzyy) + deallocate(qlm_ddgxxyz) + deallocate(qlm_ddgxyyz) + deallocate(qlm_ddgxzyz) + deallocate(qlm_ddgyyyz) + deallocate(qlm_ddgyzyz) + deallocate(qlm_ddgzzyz) + deallocate(qlm_ddgxxzz) + deallocate(qlm_ddgxyzz) + deallocate(qlm_ddgxzzz) + deallocate(qlm_ddgyyzz) + deallocate(qlm_ddgyzzz) + deallocate(qlm_ddgzzzz) + deallocate(qlm_kxx) + deallocate(qlm_kxy) + deallocate(qlm_kxz) + deallocate(qlm_kyy) + deallocate(qlm_kyz) + deallocate(qlm_kzz) + deallocate(qlm_dkxxx) + deallocate(qlm_dkxyx) + deallocate(qlm_dkxzx) + deallocate(qlm_dkyyx) + deallocate(qlm_dkyzx) + deallocate(qlm_dkzzx) + deallocate(qlm_dkxxy) + deallocate(qlm_dkxyy) + deallocate(qlm_dkxzy) + deallocate(qlm_dkyyy) + deallocate(qlm_dkyzy) + deallocate(qlm_dkzzy) + deallocate(qlm_dkxxz) + deallocate(qlm_dkxyz) + deallocate(qlm_dkxzz) + deallocate(qlm_dkyyz) + deallocate(qlm_dkyzz) + deallocate(qlm_dkzzz) + deallocate(qlm_alpha) + deallocate(qlm_betax) + deallocate(qlm_betay) + deallocate(qlm_betaz) + + deallocate(qlm_tetrad_derivs) + end subroutine deallocate_variables + +end module qlm_variables 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 |