From 76f702c98e1e81b0e81f5d3b7a39824b64799243 Mon Sep 17 00:00:00 2001 From: bentivegna Date: Sat, 10 Apr 2010 15:47:30 +0000 Subject: Changed label "ADM" -> "Landau-Lifshitz" for the existing definitions (from Weinberg); added ADM definitions (from Alcubierre's book). git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/QuasiLocalMeasures/trunk@9 ef6f4158-a8ee-47d1-ba14-cb73256398e0 --- interface.ccl | 3 ++ src/qlm_analyse.F90 | 144 ++++++++++++++++++++++++++++++++++++++++++---------- 2 files changed, 121 insertions(+), 26 deletions(-) diff --git a/interface.ccl b/interface.ccl index 0c01900..050d3d0 100644 --- a/interface.ccl +++ b/interface.ccl @@ -177,6 +177,9 @@ REAL qlm_scalars[num_surfaces] TYPE=scalar qlm_adm_energy qlm_adm_momentum_x qlm_adm_momentum_y qlm_adm_momentum_z qlm_adm_angular_momentum_x qlm_adm_angular_momentum_y qlm_adm_angular_momentum_z + qlm_ll_energy + qlm_ll_momentum_x qlm_ll_momentum_y qlm_ll_momentum_z + qlm_ll_angular_momentum_x qlm_ll_angular_momentum_y qlm_ll_angular_momentum_z } "Scalar quantities on the surface" REAL qlm_scalars_p[num_surfaces] TYPE=scalar diff --git a/src/qlm_analyse.F90 b/src/qlm_analyse.F90 index 44f0cba..887652e 100644 --- a/src/qlm_analyse.F90 +++ b/src/qlm_analyse.F90 @@ -34,7 +34,8 @@ subroutine qlm_analyse (CCTK_ARGUMENTS, hn) 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_REAL :: adm_energy, adm_mom(3), adm_amom(3) + CCTK_REAL :: ll_energy, ll_mom(3), ll_amom(3,3) CCTK_COMPLEX :: ev CCTK_REAL :: spin CCTK_REAL :: npspin @@ -95,6 +96,14 @@ subroutine qlm_analyse (CCTK_ARGUMENTS, hn) qlm_adm_angular_momentum_y(hn) = 0 qlm_adm_angular_momentum_z(hn) = 0 + qlm_ll_energy(hn) = 0 + qlm_ll_momentum_x(hn) = 0 + qlm_ll_momentum_y(hn) = 0 + qlm_ll_momentum_z(hn) = 0 + qlm_ll_angular_momentum_x(hn) = 0 + qlm_ll_angular_momentum_y(hn) = 0 + qlm_ll_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 @@ -318,68 +327,137 @@ subroutine qlm_analyse (CCTK_ARGUMENTS, hn) + coordspinz * sqrt(dtq) * weights(i) ! ADM quantities - ! Weinberg, chapter 7.6, pp. 165 ff, eqns. (7.6.22) - (7.6.24): + ! Alcubierre, Appendix A, p. 402 ff, eqns. (A.5) - (A.7): ! ADM energy - ! E_adm = (-1/16 pi) int_S(r) [h_jj,i - hij,j] n_i r^2 dOmega + ! E_adm = (1/16 pi) int_S(r) [delta^jk h_ik,j - (trace(h)),i] 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) + do c=1,3 + adm_energy = adm_energy + (delta3(a,b) * dh4(a,c,b) & + & - gu(a,b) * dh4(a,b,c) - h4(a,b) * dgu(a,b,c) ) * ss(c) + end do end do end do qlm_adm_energy(hn) = qlm_adm_energy(hn) & - & + adm_energy / (-16*pi) & - & * sqrt(dtq) * weights(i) + & + adm_energy / (16*pi) & + & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) ! ADM momentum - ! P_adm^j = (-1/16 pi) int_S(r) + ! P_adm^i = (1/8 pi) int_S(r) [K^i_j - K delta^i_j] n^j r^2 dOmega + do a=1,3 + adm_mom(a) = 0 + do b=1,3 + do c=1,3 + adm_mom(a) = adm_mom(a) & + & + gu(a,c) * kk(c,b) * ss(b) + end do + adm_mom(a) = adm_mom(a) & + + delta3(a,b) * trk * ss(b) + end do + end do + qlm_adm_momentum_x(hn) = qlm_adm_momentum_x(hn) & + & + adm_mom(1) / (8*pi) & + & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + qlm_adm_momentum_y(hn) = qlm_adm_momentum_y(hn) & + & + adm_mom(2) / (8*pi) & + & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + qlm_adm_momentum_z(hn) = qlm_adm_momentum_z(hn) & + & + adm_mom(3) / (8*pi) & + & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + + ! ADM angular momentum + ! J_adm^i = (1/8 pi) int_S(r) + ! epsilon^ijk x_j (K_kl - delta_kl K) * n^l dS + do a=1,3 + adm_amom(a) = 0 + do b=1,3 + do c=1,3 + do l=1,3 + do m=1,3 + adm_amom(a) = adm_amom(a) & + & + epsilon3(a,b,c) * xx(b) * gu(c,m) * kk(m,l) * ss(l) + end do + adm_amom(a) = adm_amom(a) & + & - epsilon3(a,b,c) * xx(b) * delta3(c,l) * trk * ss(l) + end do + end do + end do + end do + qlm_adm_angular_momentum_x(hn) = qlm_adm_angular_momentum_x(hn) & + & + adm_amom(1) / (8*pi) & + & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + qlm_adm_angular_momentum_y(hn) = qlm_adm_angular_momentum_y(hn) & + & + adm_amom(2) / (8*pi) & + & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + qlm_adm_angular_momentum_z(hn) = qlm_adm_angular_momentum_z(hn) & + & + adm_amom(3) / (8*pi) & + & * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) + + ! Landau-Lifshitz quantities + ! Weinberg, chapter 7.6, pp. 165 ff, eqns. (7.6.22) - (7.6.24): + + ! Landau-Lifshitz energy + ! E_ll = (-1/16 pi) int_S(r) [h_jj,i - hij,j] n_i r^2 dOmega + ll_energy = 0 + do a=1,3 + do b=1,3 + ll_energy = ll_energy + (dh4(b,b,a) - dh4(a,b,b)) * ss(a) + end do + end do + qlm_ll_energy(hn) = qlm_ll_energy(hn) & + & + ll_energy / (-16*pi) & + & * sqrt(dtq) * weights(i) + + ! Landau-Lifshitz momentum + ! P_ll^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 + ! P_ll^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 + ll_mom(a) = 0 do b=1,3 - adm_mom(a) = adm_mom(a) & + ll_mom(a) = ll_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) & + qlm_ll_momentum_x(hn) = qlm_ll_momentum_x(hn) & + & + ll_mom(1) / (-16*pi) & & * sqrt(dtq) * weights(i) - qlm_adm_momentum_y(hn) = qlm_adm_momentum_y(hn) & - & + adm_mom(2) / (-16*pi) & + qlm_ll_momentum_y(hn) = qlm_ll_momentum_y(hn) & + & + ll_mom(2) / (-16*pi) & & * sqrt(dtq) * weights(i) - qlm_adm_momentum_z(hn) = qlm_adm_momentum_z(hn) & - & + adm_mom(3) / (-16*pi) & + qlm_ll_momentum_z(hn) = qlm_ll_momentum_z(hn) & + & + ll_mom(3) / (-16*pi) & & * sqrt(dtq) * weights(i) - ! ADM angular momentum - ! J_adm^jk = (-1/16 pi) int_S(r) + ! Landau-Lifshitz angular momentum + ! J_ll^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 + ll_amom(a,b) = 0 do c=1,3 - adm_amom(a,b) = adm_amom(a,b) + ss(c) * ( & + ll_amom(a,b) = ll_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) & + qlm_ll_angular_momentum_x(hn) = qlm_ll_angular_momentum_x(hn) & + & + ll_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) & + qlm_ll_angular_momentum_y(hn) = qlm_ll_angular_momentum_y(hn) & + & + ll_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) & + qlm_ll_angular_momentum_z(hn) = qlm_ll_angular_momentum_z(hn) & + & + ll_amom(1,2) / (-16*pi) & & * sqrt(dtq) * weights(i) end do @@ -483,6 +561,20 @@ subroutine qlm_analyse (CCTK_ARGUMENTS, hn) call CCTK_INFO (msg) write (msg, '(" ADM angular momentum z: ",g16.6)') qlm_adm_angular_momentum_z(hn) call CCTK_INFO (msg) + write (msg, '(" Landau-Lifshitz energy: ",g16.6)') qlm_ll_energy(hn) + call CCTK_INFO (msg) + write (msg, '(" Landau-Lifshitz momentum x: ",g16.6)') qlm_ll_momentum_x(hn) + call CCTK_INFO (msg) + write (msg, '(" Landau-Lifshitz momentum y: ",g16.6)') qlm_ll_momentum_y(hn) + call CCTK_INFO (msg) + write (msg, '(" Landau-Lifshitz momentum z: ",g16.6)') qlm_ll_momentum_z(hn) + call CCTK_INFO (msg) + write (msg, '(" Landau-Lifshitz angular momentum x: ",g16.6)') qlm_ll_angular_momentum_x(hn) + call CCTK_INFO (msg) + write (msg, '(" Landau-Lifshitz angular momentum y: ",g16.6)') qlm_ll_angular_momentum_y(hn) + call CCTK_INFO (msg) + write (msg, '(" Landau-Lifshitz angular momentum z: ",g16.6)') qlm_ll_angular_momentum_z(hn) + call CCTK_INFO (msg) end if -- cgit v1.2.3