diff options
Diffstat (limited to 'src/boost.F')
| -rw-r--r-- | src/boost.F | 462 |
1 files changed, 462 insertions, 0 deletions
diff --git a/src/boost.F b/src/boost.F new file mode 100644 index 0000000..6b993ea --- /dev/null +++ b/src/boost.F @@ -0,0 +1,462 @@ +c This subroutine calculates the 4-metric and its inverse at an event, +c taking into account an optional Lorentz boost and an optional rotation. +c The model is first rotated and then boosted, such that the boost is +c applied to the rotated model. +c $Header$ +c +c The coordinates are +c Cx(a) = Cactus $x^a$ +c Mx(a) = Model $X^a$ +c The 4-metrics are +c Cgdd(a,b) = Cactus $g_{ab}$ Cguu(a,b) = Cactus $g^{ab}$ +c Mgdd(a,b) = Model $g_{ab}$ Mguu(a,b) = Model $g^{ab}$ +c +c For a definition of the Euler angles in the conventions used below, see +c http://mathworld.wolfram.com/EulerAngles.html +c Another useful resource may be +c http://en.wikipedia.org/wiki/Euler_angles +c although this uses (on 2006-11-29) different conventions. +c +c This file is copyright (c) 2003 by Jonathan Thornburg <jthorn@aei.mpg.de>. +c This file is covered by the GNU GPL license; see the files ../README +c and ../COPYING for details. +c + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#include "param_defs.inc" + + subroutine Exact__metric( + $ decoded_exact_model, + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi) + + implicit none + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + +c input arguments + CCTK_INT decoded_exact_model + CCTK_REAL x, y, z, t + +c output arguments + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi + +c static local variables describing Lorentz transformation + logical firstcall + data firstcall /.true./ + CCTK_REAL gamma + CCTK_REAL vv(3), nn(3) + CCTK_REAL parallel(3,3), perp(3,3) + CCTK_REAL Cx_par(3), Cx_perp(3) + CCTK_REAL partial_Mx_wrt_Cx(0:3,0:3) + CCTK_REAL partial_Cx_wrt_Mx(0:3,0:3) + CCTK_REAL R(0:3,0:3) + save firstcall + save gamma + save vv + save parallel, perp + save partial_Mx_wrt_Cx + save partial_Cx_wrt_Mx + save R +c$omp threadprivate (firstcall, gamma, vv, parallel, perp, +c$omp+ partial_Mx_wrt_Cx,partial_Cx_wrt_Mx, R) + +c coordinates and 4-metric + CCTK_REAL Cx(0:3) + CCTK_REAL Cgdd(0:3,0:3), Cguu(0:3,0:3) + CCTK_REAL Mx(0:3) + CCTK_REAL Mgdd(0:3,0:3), Mguu(0:3,0:3) + CCTK_REAL Nx(0:3) + CCTK_REAL Ngdd(0:3,0:3), Nguu(0:3,0:3) + +c misc temps + CCTK_REAL vnorm, vnormsq + CCTK_REAL delta_ij + CCTK_REAL Cx_par_i, Cx_perp_i + CCTK_REAL vdotCx + CCTK_REAL Cgdd_ab, Cguu_ab + CCTK_REAL cos_phi, sin_phi + CCTK_REAL cos_theta, sin_theta + CCTK_REAL cos_psi, sin_psi + CCTK_REAL R_phi(0:3,0:3), R_theta(0:3,0:3), R_psi(0:3,0:3) + character*1000 warn_buffer + +c flags, array indices, etc + logical Tmunu_flag + integer i, j, k, l + integer Ca, Cb, MA, MB + +c constants + integer n + parameter (n = 3) + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c optimized fast-path if no Lorentz boost and no rotation +c + if ( (boost_vx .eq. 0.0) + $ .and. (boost_vy .eq. 0.0) + $ .and. (boost_vz .eq. 0.0) + $ .and. (rotation_euler_phi .eq. 0.0) + $ .and. (rotation_euler_theta .eq. 0.0) + $ .and. (rotation_euler_psi .eq. 0.0)) then + call Exact__metric_for_model( + $ decoded_exact_model, + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ psi, + $ Tmunu_flag) + return + end if + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c the rest of this function is the Lorentz-boost case: +c - Lorentz-transform Cactus coordinates --> Model coordinates +c - compute Model 4-metric and inverse at Model coordinates +c - tensor-transform 4-metric and inverse from Model coordinates +c --> Cactus coordinates +c +c All the equations used are given in ../doc/documentation.tex +c + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c compute Lorentz transformation information on first call +c + if (firstcall) then + +c boost velocity + vv(1) = boost_vx + vv(2) = boost_vy + vv(3) = boost_vz + +c Lorentz gamma factor, unit vector in direction of boost velocity + vnormsq = 0.0 + do i = 1,n + vnormsq = vnormsq + vv(i)*vv(i) + end do + gamma = 1.0 / sqrt(1.0 - vnormsq) + vnorm = sqrt(vnormsq) + if ( (boost_vx .eq. 0.0) + $ .and. (boost_vy .eq. 0.0) + $ .and. (boost_vz .eq. 0.0)) then + nn(1) = 1 + nn(2) = 0 + nn(3) = 0 + else + do i = 1,n + nn(i) = vv(i) / vnorm + end do + end if + +c projection operators parallel(*,*) and perp(*,*) + do j = 1,n + do i = 1,n + parallel(i,j) = nn(i) * nn(j) + if (i .eq. j) then + delta_ij = 1.0 + else + delta_ij = 0.0 + end if + perp(i,j) = delta_ij - parallel(i,j) + end do + end do + +c partial derivatives of Model coordinates with respect to Cactus coordinates + partial_Mx_wrt_Cx(0,0) = gamma + do i = 1,n + partial_Mx_wrt_Cx(0,i) = -gamma*vv(i) + end do + do i = 1,n + partial_Mx_wrt_Cx(i,0) = -gamma*vv(i) + do j=1,n + partial_Mx_wrt_Cx(i,j) = gamma*parallel(i,j) + perp(i,j) + end do + end do + +c partial derivatives of Cactus coordinates with respect to Model coordinates + partial_Cx_wrt_Mx(0,0) = gamma + do i = 1,n + partial_Cx_wrt_Mx(0,i) = + gamma*vv(i) + end do + do i = 1,n + partial_Cx_wrt_Mx(i,0) = + gamma*vv(i) + do j=1,n + partial_Cx_wrt_Mx(i,j) = gamma*parallel(i,j) + perp(i,j) + end do + end do + +c Sines and cosines of rotation angles + cos_phi = cos (rotation_euler_phi) + sin_phi = sin (rotation_euler_phi) + cos_theta = cos (rotation_euler_theta) + sin_theta = sin (rotation_euler_theta) + cos_psi = cos (rotation_euler_psi) + sin_psi = sin (rotation_euler_psi) + +c Set up individual rotation matrices + R_phi(0,0) = 1 + R_phi(0,1) = 0 + R_phi(0,2) = 0 + R_phi(0,3) = 0 + R_phi(1,0) = 0 + R_phi(1,1) = + cos_phi + R_phi(1,2) = + sin_phi + R_phi(1,3) = 0 + R_phi(2,0) = 0 + R_phi(2,1) = - sin_phi + R_phi(2,2) = + cos_phi + R_phi(2,3) = 0 + R_phi(3,0) = 0 + R_phi(3,1) = 0 + R_phi(3,2) = 0 + R_phi(3,3) = 1 + + R_theta(0,0) = 1 + R_theta(0,1) = 0 + R_theta(0,2) = 0 + R_theta(0,3) = 0 + R_theta(1,0) = 0 + R_theta(1,1) = 1 + R_theta(1,2) = 0 + R_theta(1,3) = 0 + R_theta(2,0) = 0 + R_theta(2,1) = 0 + R_theta(2,2) = + cos_theta + R_theta(2,3) = + sin_theta + R_theta(3,0) = 0 + R_theta(3,1) = 0 + R_theta(3,2) = - sin_theta + R_theta(3,3) = + cos_theta + + R_psi(0,0) = 1 + R_psi(0,1) = 0 + R_psi(0,2) = 0 + R_psi(0,3) = 0 + R_psi(1,0) = 0 + R_psi(1,1) = + cos_psi + R_psi(1,2) = + sin_psi + R_psi(1,3) = 0 + R_psi(2,0) = 0 + R_psi(2,1) = - sin_psi + R_psi(2,2) = + cos_psi + R_psi(2,3) = 0 + R_psi(3,0) = 0 + R_psi(3,1) = 0 + R_psi(3,2) = 0 + R_psi(3,3) = 1 + +c Combine individual rotation matrices + do i = 0,n + do j = 0,n + R(i,j) = 0 + do k = 0,n + do l = 0,n + R(i,j) = R(i,j) + R_psi(i,k) * R_theta(k,l) * R_phi(l,j) + end do + end do + end do + end do + +c Notes that help me (Erik Schnetter) think: +c This considers a rotation with phi=0, theta=pi/2, psi=0. +c Nx(i) = Nx(i) + R(j,i) * Mx(j) +c Nx(1) = - Mx(3) +c Nx(3) = Mx(1) +c Mgxx(1,1) = Ngxx(3,3) +c Mgxx(1,3) = - Ngxx(1,3) +c Mgxx(3,3) = Ngxx(1,1) +c Mgxx(i,j) = R(i,k) R(j,l) Ngxx(k,l) [correct] +c Mgxx(i,j) = R(k,i) R(l,j) Ngxx(k,l) [correct] +c Mbetax(1) = - Nbetax(3) +c Mbetax(3) = Nbetax(1) +c Mbetax(i) + R(i,j) Nbetax(j) [wrong] +c Mbetax(i) + R(j,i) Nbetax(j) [correct] + + firstcall = .false. + end if + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c +c compute flat-space components of Cx(*) parallel and perpendicular to vv(*) +c + Cx(0) = t + Cx(1) = x + Cx(2) = y + Cx(3) = z + + do i=1,n + Cx_par_i = 0.0 + Cx_perp_i = 0.0 + do j=1,n + Cx_par_i = Cx_par_i + parallel(i,j)*Cx(j) + Cx_perp_i = Cx_perp_i + perp (i,j)*Cx(j) + end do + Cx_par (i) = Cx_par_i + Cx_perp(i) = Cx_perp_i + end do + +c +c Lorentz-transform and rotate the Cactus coordinate +c to get the Model coordinates +c + +c Boost + vdotCx = 0.0 + do i = 1,n + vdotCx = vdotCx + vv(i)*Cx(i) + end do + + Mx(0) = gamma * (Cx(0) - vdotCx) + do i=1,n + Mx(i) = gamma * (Cx_par(i) - vv(i)*Cx(0)) + Cx_perp(i) + end do + +c Rotation + do i=0,n + Nx(i) = 0 + do j = 0,n + Nx(i) = Nx(i) + R(i,j) * Mx(j) + end do + end do + + +c +c compute the Model 4-metric and inverse 4-metric at the Model coordinates +c + call Exact__metric_for_model( + $ decoded_exact_model, + $ Nx(1), Nx(2), Nx(3), Nx(0), + $ Ngdd(0,0), Ngdd(0,1), Ngdd(0,2), Ngdd(0,3), + $ Ngdd(1,1), Ngdd(2,2), Ngdd(3,3), + $ Ngdd(1,2), Ngdd(2,3), Ngdd(1,3), + $ Nguu(0,0), Nguu(0,1), Nguu(0,2), Nguu(0,3), + $ Nguu(1,1), Nguu(2,2), Nguu(3,3), + $ Nguu(1,2), Nguu(2,3), Nguu(1,3), + $ psi, Tmunu_flag) + + if (Tmunu_flag) then + write (warn_buffer, '(a,i8,a,a)') + $ 'exact_model = ', decoded_exact_model, + $ 'sets the stress-energy tensor', + $ ' ==> we cannot Lorentz-boost or rotate it!' + call CCTK_WARN(0, warn_buffer) + end if + +c +c symmetrize the Model 4-metric and inverse 4-metric arrays +c (the Exact__metric_for_model() call only set the upper triangles) +c + Ngdd(1,0) = Ngdd(0,1) + Ngdd(2,0) = Ngdd(0,2) + Ngdd(2,1) = Ngdd(1,2) + Ngdd(3,0) = Ngdd(0,3) + Ngdd(3,1) = Ngdd(1,3) + Ngdd(3,2) = Ngdd(2,3) + + Nguu(1,0) = Nguu(0,1) + Nguu(2,0) = Nguu(0,2) + Nguu(2,1) = Nguu(1,2) + Nguu(3,0) = Nguu(0,3) + Nguu(3,1) = Nguu(1,3) + Nguu(3,2) = Nguu(2,3) + +c +c tensor-transorm (the upper triangle of) the 4-metric and inverse 4-metric +c + +c Rotations + do i = 0,n + do j = 0,n + Mgdd(i,j) = 0 + Mguu(i,j) = 0 + do k = 0,n + do l = 0,n + Mgdd(i,j) = Mgdd(i,j) + R(k,i) * R(l,j) * Ngdd(k,l) +c The inverse of R is also its transpose. That means that the +c transpose of the inverse, which you would use for g^ij, is just R +c again. + Mguu(i,j) = Mguu(i,j) + R(k,i) * R(l,j) * Nguu(k,l) + end do + end do + end do + end do + +c Boost + do Ca = 0,n + do Cb = Ca,n + Cgdd_ab = 0.0 + do Ma = 0,n + do Mb = 0,n + Cgdd_ab = Cgdd_ab + $ + Mgdd(Ma,Mb) + $ * partial_Mx_wrt_Cx(Ma,Ca) + $ * partial_Mx_wrt_Cx(Mb,Cb) + end do + end do + Cgdd(Ca,Cb) = Cgdd_ab + end do + end do + + do Ca = 0,n + do Cb = Ca,n + Cguu_ab = 0.0 + do Ma = 0,n + do Mb = 0,n + Cguu_ab = Cguu_ab + $ + Mguu(Ma,Mb) + $ * partial_Cx_wrt_Mx(Ca,Ma) + $ * partial_Cx_wrt_Mx(Cb,Mb) + end do + end do + Cguu(Ca,Cb) = Cguu_ab + end do + end do + +c +c unpack the Cactus-coordinates 4-metric and inverse 4-metric +c into the corresponding output arguments +c + gdtt = Cgdd(0,0) + gdtx = Cgdd(0,1) + gdty = Cgdd(0,2) + gdtz = Cgdd(0,3) + gdxx = Cgdd(1,1) + gdxy = Cgdd(1,2) + gdxz = Cgdd(1,3) + gdyy = Cgdd(2,2) + gdyz = Cgdd(2,3) + gdzz = Cgdd(3,3) + + gutt = Cguu(0,0) + gutx = Cguu(0,1) + guty = Cguu(0,2) + gutz = Cguu(0,3) + guxx = Cguu(1,1) + guxy = Cguu(1,2) + guxz = Cguu(1,3) + guyy = Cguu(2,2) + guyz = Cguu(2,3) + guzz = Cguu(3,3) + + end |