diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/qlm_derivs.F90 | 103 | ||||
| -rw-r--r-- | src/qlm_twometric.F90 | 103 | ||||
| -rw-r--r-- | src/qlm_weyl_scalars.F90 | 47 |
3 files changed, 50 insertions, 203 deletions
diff --git a/src/qlm_derivs.F90 b/src/qlm_derivs.F90 index 24a3552..a2a73ef 100644 --- a/src/qlm_derivs.F90 +++ b/src/qlm_derivs.F90 @@ -10,9 +10,7 @@ module qlm_derivs public operator(.dot.) public deriv public deriv2 - public deriv3 public timederiv - public timederiv2 interface deriv module procedure rderiv @@ -22,18 +20,8 @@ module qlm_derivs 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.) @@ -163,47 +151,6 @@ contains 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 @@ -235,54 +182,4 @@ contains 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_twometric.F90 b/src/qlm_twometric.F90 index e450ebb..7a31c42 100644 --- a/src/qlm_twometric.F90 +++ b/src/qlm_twometric.F90 @@ -23,17 +23,14 @@ subroutine qlm_calc_twometric (CCTK_ARGUMENTS, hn) 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 :: gg(3,3), dgg(3,3,3) + CCTK_REAL :: ee(3,2), dee(3,2,2) + CCTK_REAL :: qq(2,2), dqq(2,2,2) CCTK_REAL :: delta_space(2) integer :: i, j - integer :: a, b, c, d, e, f, g, h + integer :: a, b, c, d, e, f if (veryverbose/=0) then call CCTK_INFO ("Calculating two-metric") @@ -77,49 +74,6 @@ subroutine qlm_calc_twometric (CCTK_ARGUMENTS, hn) 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) @@ -128,10 +82,6 @@ subroutine qlm_calc_twometric (CCTK_ARGUMENTS, hn) 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 @@ -160,52 +110,12 @@ subroutine qlm_calc_twometric (CCTK_ARGUMENTS, hn) 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) @@ -213,9 +123,6 @@ subroutine qlm_calc_twometric (CCTK_ARGUMENTS, hn) 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 @@ -229,7 +136,5 @@ subroutine qlm_calc_twometric (CCTK_ARGUMENTS, hn) 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_weyl_scalars.F90 b/src/qlm_weyl_scalars.F90 index ade2e31..73d2926 100644 --- a/src/qlm_weyl_scalars.F90 +++ b/src/qlm_weyl_scalars.F90 @@ -26,7 +26,7 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) 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 :: kk(3,3), dkk(3,3,3), kk_dot(3,3) CCTK_REAL :: tt(3,3) - CCTK_REAL :: dtg, gu(3,3), dgu(3,3,3), gamma(3,3,3), dgamma(3,3,3,3), ri(3,3) + CCTK_REAL :: dtg, gu(3,3), dgu(3,3,3), gamma(3,3,3), dgamma(3,3,3,3), ri(3,3), rsc 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) @@ -34,6 +34,10 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) 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) + CCTK_REAL :: ss(0:3) + CCTK_REAL :: nabla_ll(0:3,0:3), nabla_nn(0:3,0:3) + CCTK_REAL :: nabla_ss(0:3,0:3) + CCTK_REAL :: trkAH, kk_kk, tmpR integer :: i, j integer :: a, b, c, d @@ -182,6 +186,13 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) mm(2) = qlm_m2(i,j,hn) mm(3) = qlm_m3(i,j,hn) + ss = (ll - nn) / sqrt(two) + + nabla_ll = qlm_tetrad_derivs(i,j)%nabla_ll + nabla_nn = qlm_tetrad_derivs(i,j)%nabla_nn + + nabla_ss = (nabla_ll - nabla_nn) / sqrt(two) + ! Calculate 4-metric @@ -191,6 +202,7 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) call calc_connections (gu, dgg, gamma) call calc_connectionderivs (gu, dgg, dgu, ddgg, dgamma) call calc_ricci (gamma, dgamma, ri) + call calc_trace (gu, ri, rsc) call calc_3metricdot_simple (kk, gg_dot) call calc_3metricderivdot_simple (dkk, dgg_dot) @@ -278,6 +290,37 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) & + 2 * real (qlm_psi2(i,j,hn)) & & + 4 * qlm_lambda(i,j,hn) + + + trkAH = 0 + do a=1,3 + do b=1,3 + trkAH = trkAH + (gu(a,b) - ss(a) * ss(b)) * nabla_ss(a,b) + end do + end do + + kk_kk = 0 + do a=1,3 + do b=1,3 + do c=1,3 + do d=1,3 + kk_kk = kk_kk + & + (gu(a,c) - ss(a) * ss(c)) * nabla_ss(a,b) * & + (gu(b,d) - ss(b) * ss(d)) * nabla_ss(c,d) + end do + end do + end do + end do + + tmpR = 0 + do a=1,3 + do b=1,3 + tmpR = tmpR + ri(a,b) * ss(a) * ss(b) + end do + end do + + qlm_rsc(i,j,hn) = rsc - 2*tmpR + trkAH**2 - kk_kk + end do end do @@ -305,5 +348,7 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) 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) + + call set_boundary (CCTK_PASS_FTOF, hn, qlm_rsc(:,:,hn), +1) end subroutine qlm_calc_weyl_scalars |