From 661c56c847f4334feb53fd8d94551b2ab1efdb9a Mon Sep 17 00:00:00 2001 From: allen Date: Mon, 1 Nov 1999 14:56:51 +0000 Subject: This commit was generated by cvs2svn to compensate for changes in r2, which included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Extract/trunk@3 5301f0c2-dbc4-4cee-b2f5-8d7afba4d129 --- src/ADMmass_integrand3D.F | 170 ++++++++ src/ADMmass_integrand3D_int.F | 27 ++ src/D2_extract.F | 918 +++++++++++++++++++++++++++++++++++++++ src/D2_extract_int.F | 42 ++ src/D3_extract.F | 221 ++++++++++ src/D3_extract_int.F | 44 ++ src/D3_to_D2.F | 398 +++++++++++++++++ src/D3_to_D2_int.F | 45 ++ src/Extract.F | 814 ++++++++++++++++++++++++++++++++++ src/cartesian_to_spherical.F | 116 +++++ src/cartesian_to_spherical_int.F | 24 + src/energy.f | 205 +++++++++ src/make.code.defn | 27 ++ src/make.code.deps | 2 + src/met_rad_der.F | 91 ++++ src/met_rad_der_int.F | 24 + src/momentum_integrand3D.F | 188 ++++++++ src/momentum_integrand3D_int.F | 30 ++ src/spin_integrand3D.F | 193 ++++++++ src/spin_integrand3D_int.F | 30 ++ src/unphysical_to_physical.F | 56 +++ src/unphysical_to_physical_int.F | 21 + 22 files changed, 3686 insertions(+) create mode 100644 src/ADMmass_integrand3D.F create mode 100644 src/ADMmass_integrand3D_int.F create mode 100644 src/D2_extract.F create mode 100644 src/D2_extract_int.F create mode 100644 src/D3_extract.F create mode 100644 src/D3_extract_int.F create mode 100644 src/D3_to_D2.F create mode 100644 src/D3_to_D2_int.F create mode 100644 src/Extract.F create mode 100644 src/cartesian_to_spherical.F create mode 100644 src/cartesian_to_spherical_int.F create mode 100644 src/energy.f create mode 100644 src/make.code.defn create mode 100644 src/make.code.deps create mode 100644 src/met_rad_der.F create mode 100644 src/met_rad_der_int.F create mode 100644 src/momentum_integrand3D.F create mode 100644 src/momentum_integrand3D_int.F create mode 100644 src/spin_integrand3D.F create mode 100644 src/spin_integrand3D_int.F create mode 100644 src/unphysical_to_physical.F create mode 100644 src/unphysical_to_physical_int.F (limited to 'src') diff --git a/src/ADMmass_integrand3D.F b/src/ADMmass_integrand3D.F new file mode 100644 index 0000000..06b1198 --- /dev/null +++ b/src/ADMmass_integrand3D.F @@ -0,0 +1,170 @@ +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz, + & gyy,gyz,gzz,ADMmass_int,Psi,Psi_power) + +c ------------------------------------------------------------------------ +c +c Estimates the ADM mass at a given radius using Equation (7) from +c O Murchadha and York, Phys Rev D, 10, 1974 page 2345 +c +c ------------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & ADMmass_int + +c Local variables, here only + INTEGER :: + & i,j,k,ip + CCTK_REAL,PARAMETER :: + & half = 0.5D0 + CCTK_REAL :: + & rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy, + & uyz,uzz,term1,term2,term3,dxx_y,dxx_z,dxy_x,dxy_y,dxy_z, + & dyy_x,dxz_x,dxz_y,dxz_z,dyz_x,dzz_x,dyy_z,dyz_y,dyz_z,dzz_y, + & Pi,idet,p,pip,pim,pjp,pjm,pkp,pkm + +c ------------------------------------------------------------------------ + + Pi = ACOS(-1D0) + +c Because other codes evolve Psi**4 + SELECT CASE (Psi_power) + + CASE (1) + ip = 4 + + CASE (4) + ip = 1 + + CASE DEFAULT + WRITE(*,*) "This value of Psi_power is not supported" + + END SELECT + + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + ux = (x(i)-origin(1))/rad + uy = (y(j)-origin(2))/rad + uz = (z(k)-origin(3))/rad + +c Abbreviations for metric coefficients +c ------------------------------------- + p = psi(i,j,k)**ip + + dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k) + dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k) + dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k) + +c Determinant of 3-metric +c ----------------------- + det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz + & - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2)) + + idet = 1.0/det + +c Inverse 3-metric +c ---------------- + uxx = idet*(dyy*dzz - dyz**2) + uyy = idet*(dxx*dzz - dxz**2) + uzz = idet*(dxx*dyy - dxy**2) + uxy = idet*(dxz*dyz - dzz*dxy) + uxz = idet*(dxy*dyz - dyy*dxz) + uyz = idet*(dxy*dxz - dxx*dyz) + +c Derivatives of the 3-metric +c --------------------------- + pip = psi(i+1,j,k)**ip + pim = psi(i-1,j,k)**ip + pjp = psi(i,j+1,k)**ip + pjm = psi(i,j-1,k)**ip + pkp = psi(i,j,k+1)**ip + pkm = psi(i,j,k-1)**ip + + dxx_y = half/Dy*(pjp*gxx(i,j+1,k)-pjm*gxx(i,j-1,k)) + dxx_z = half/Dz*(pkp*gxx(i,j,k+1)-pkm*gxx(i,j,k-1)) + dxy_x = half/Dx*(pip*gxy(i+1,j,k)-pim*gxy(i-1,j,k)) + dxy_y = half/Dy*(pjp*gxy(i,j+1,k)-pjm*gxy(i,j-1,k)) + dxy_z = half/Dz*(pkp*gxy(i,j,k+1)-pkm*gxy(i,j,k-1)) + dyy_x = half/Dx*(pip*gyy(i+1,j,k)-pim*gyy(i-1,j,k)) + dyy_z = half/Dz*(pkp*gyy(i,j,k+1)-pkm*gyy(i,j,k-1)) + dxz_x = half/Dx*(pip*gxz(i+1,j,k)-pim*gxz(i-1,j,k)) + dxz_y = half/Dy*(pjp*gxz(i,j+1,k)-pjm*gxz(i,j-1,k)) + dxz_z = half/Dz*(pkp*gxz(i,j,k+1)-pkm*gxz(i,j,k-1)) + dyz_x = half/Dx*(pip*gyz(i+1,j,k)-pim*gyz(i-1,j,k)) + dyz_y = half/Dy*(pjp*gyz(i,j+1,k)-pjm*gyz(i,j-1,k)) + dyz_z = half/Dz*(pkp*gyz(i,j,k+1)-pkm*gyz(i,j,k-1)) + dzz_x = half/Dx*(pip*gzz(i+1,j,k)-pim*gzz(i-1,j,k)) + dzz_y = half/Dy*(pjp*gzz(i,j+1,k)-pjm*gzz(i,j-1,k)) + + term1 = uxy*(dxx_y-dxy_x)+uxz*(dxx_z-dxz_x) + & +uyy*(dxy_y-dyy_x)+uyz*(dxy_z-dyz_x) + & +uyz*(dxz_y-dyz_x)+uzz*(dxz_z-dzz_x) + + term2 = uyz*(dyy_z-dyz_y)+uxy*(dyy_x-dxy_y) + & +uzz*(dyz_z-dzz_y)+uxz*(dyz_x-dxz_y) + & +uxz*(dxy_z-dxz_y)+uxx*(dxy_x-dxx_y) + + term3 = uxz*(dzz_x-dxz_z)+uyz*(dzz_y-dyz_z) + & +uxx*(dxz_x-dxx_z)+uxy*(dxz_y-dxy_z) + & +uxy*(dyz_x-dxy_z)+uyy*(dyz_y-dyy_z) + + ADMmass_int(i,j,k) = 1.0D0/16.0D0/Pi*(ux*term1+ + & uy*term2+uz*term3)*SQRT(det)*rad**2 + + ELSE + + ADMmass_int(i,j,k) = 0.0D0 + + ENDIF + + ENDDO + ENDDO + ENDDO + + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + ADMmass_int(1,j,k) = ADMmass_int(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + ADMmass_int(i,1,k) = ADMmass_int(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + ADMmass_int(i,j,1) = ADMmass_int(i,j,2) + ENDDO + ENDDO + + + END SUBROUTINE ADMmass_integrand3D diff --git a/src/ADMmass_integrand3D_int.F b/src/ADMmass_integrand3D_int.F new file mode 100644 index 0000000..c103f55 --- /dev/null +++ b/src/ADMmass_integrand3D_int.F @@ -0,0 +1,27 @@ + +#include "cctk.h" + + MODULE ADMmass_integrand3D_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy, + & gxz,gyy,gyz,gzz,ADMmass_int,Psi,Psi_power) + IMPLICIT NONE + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & ADMmass_int + END SUBROUTINE ADMmass_integrand3D + + END INTERFACE + + END MODULE ADMmass_integrand3D_int diff --git a/src/D2_extract.F b/src/D2_extract.F new file mode 100644 index 0000000..3b24ae0 --- /dev/null +++ b/src/D2_extract.F @@ -0,0 +1,918 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11, + & g12,g13,g22,g23,g33,dg22,dg23,dg33, + & do_ADMmass,ADMmass_int1,ADMmass_int2, + & do_momentum,momentum_int1,momentum_int2,momentum_int3, + & do_spin,spin_int1,spin_int2,spin_int3, + & mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt) + +c ------------------------------------------------------------------------ +c +c Extract the odd and even parity gauge invariant variables +c at a given isotropic radius (eta) from the center of a +c slightly perturbed Schwarzschild spacetime. +c +c The polar coordinates and metric variables must be given on +c an equally spaced grid, which is offset from the axis. +c +c Notice that for octant symmetry we only need calculate the real +c parts of the gauge invariant variables for the even-l,m modes. +c This subroutine is optimised (!!) for octant symmetry. +c +c ------------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & igrid,l,m + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & eta,theta(:),phi(:) + CCTK_REAL,DIMENSION (:,:) :: + & g00,g11,g12,g13,g22,g23,g33,dg22,dg23,dg33, + & ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + +c Output variables + CCTK_REAL,INTENT(OUT) :: + & dtaudt,ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:), + & momentum(3),spin(3) + +c Local Variables + LOGICAL, PARAMETER :: + & debug = .FALSE. + INTEGER :: + & i,j,il,im,lmin,lmax,mmin,mmax,lstep,mstep + CCTK_REAL,PARAMETER :: + & zero = 0.0D0, + & half = 0.5D0, + & one = 1.0D0, + & two = 2.0D0, + & four = 4.0D0, + & six = 6.0D0, + & eight = 8.0D0 + CCTK_REAL :: + & Dt,Dp,st,ist,drschdri,dridrsch, + & g22comb,g23comb,g33comb,S,rsch2,Lambda, + & Pi, + & fac_h1,fac_H2,fac_K,fac_G,fac_c1,fac_c2, + & fac_dG,fac_dK,fac_dc2,fac1,fac2,fac3,max_g11, + & sph_g11,sph_g22,sph_g33,sph_dg22,sph_dg33,error, + & sphere_int + CCTK_REAL,DIMENSION(2) :: + & Y,Y1,Y2,Y3,Y4,h1,H2,K,G,c1,c2,dG,dK,dc2 + CCTK_REAL,DIMENSION(Nt+1,Np+1) :: + & temp,h1i,H2i,Gi,Ki,c1i,c2i,dGi,dKi,dc2i + + +c ------------------------------------------------------------------------ +c +c 0. Initial things +c +c ------------------------------------------------------------------------ + + Pi = two*ASIN(one) + + SELECT CASE (igrid) + CASE (0) ! full grid + Dt = Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + lstep = 1 + mstep = 1 + CASE (1) ! octant grid + Dt = half*Pi/DBLE(Nt) + Dp = half*Pi/DBLE(Np) + lstep = 2 + mstep = 2 + CASE (2) ! cartoon + Dt = Pi/DBLE(Nt) + Dp = 2D0*Pi + lstep = 2 + mstep = 2 + CASE (3) ! bitant + Dt = half*Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + lstep = 1 + mstep = 1 + CASE DEFAULT + WRITE (*,*) "Grid incorrectly set in D2_extract" + STOP + END SELECT + +c Calculate ADM mass +c ------------------ + IF (do_ADMmass(1) == 1) THEN + +c Calculate ADM mass using standard formula + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = ADMmass_int1(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + ADMmass(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + ENDIF + IF (do_ADMmass(2) == 1) THEN + +c Calculate ADM mass using conformal formula + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = ADMmass_int2(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + ADMmass(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + END IF + +c Calculate momentum +c ------------------ + IF (do_momentum == 1) THEN + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = momentum_int1(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + momentum(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = momentum_int2(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + momentum(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = momentum_int3(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + momentum(3) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + ENDIF + +c Calculate spin +c -------------- + IF (do_spin == 1) THEN + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = spin_int1(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + spin(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = spin_int2(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + spin(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = spin_int3(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + spin(3) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + ENDIF + +c ------------------------------------------------------------------ +c +c 1. Find Spherical Parts of Metric Components +c +c ------------------------------------------------------------------ + +c ... g00 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g00(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + dtaudt = sqrt(one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)) + +c ... g11 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g11(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + sph_g11 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + +c ... g22 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g22(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + sph_g22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + +c ... dg22 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = dg22(i-1,j-1)*DSIN(theta(i-1)) + ENDDO + ENDDO + sph_dg22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + +c ... g33/sin(theta)**2 + DO j = 2, Np+1 + DO i = 2, Nt+1 + temp(i,j) = g33(i-1,j-1)/DSIN(theta(i-1)) + ENDDO + ENDDO + sph_g33 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + +c Calculate error in orthonormality of spherical harmonic +c DO j = 2, Np+1 +c DO i = 2, Nt+1 +c CALL spher_harm_combs(theta(i),phi(j),l,m,Y,Y1,Y2,Y3,Y4) +c temp(i,j) = (Y(1)*Y(1)+Y(2)*Y(2))*sin(theta(i-1)) +c ENDDO +c ENDDO +c error = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) +c print *,"Error is",one-error + + + +c ------------------------------------------------------------------ +c +c 1. Compute the Schwarzschild radius +c +c ------------------------------------------------------------------ + + DO j = 2, Np+1 + DO i = 2, Nt+1 +c temp(i,j) = g22(i-1,j-1)*SIN(theta(i-1)) + temp(i,j) = (g22(i-1,j-1)+g33(i-1,j-1)/SIN(theta(i-1))**2) + & *SIN(theta(i-1))/(two) + ENDDO + ENDDO + + rsch2= one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + rsch = DSQRT(rsch2) + + +c ------------------------------------------------------------------ +c +c 2. Compute the derivative of the schwarzschild radius with +c respect to the isotropic radius eta. +c +c ------------------------------------------------------------------ + + DO j = 2, Np+1 + DO i = 2, Nt+1 +c temp(i,j) = (dg22(i-1,j-1))*SIN(theta(i-1)) + temp(i,j) = (dg22(i-1,j-1)+dg33(i-1,j-1)/SIN(theta(i-1))**2) + & *SIN(theta(i-1))/(two) + ENDDO + ENDDO + + drschdri = one/(eight*Pi*rsch)* + & sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp) + + dridrsch = one/drschdri + + +c ------------------------------------------------------------------ +c +c 3. Calculate the Schwarzschild mass and S factor +c +c ------------------------------------------------------------------ + + S = drschdri**2/sph_g11 + mass = rsch*(one-S)/two + + + +c ------------------------------------------------------------------ +c +c 4. Subtract spherical part from metric +c Do not know how important this is +c +c ------------------------------------------------------------------ + + DO i=1,Nt + DO j=1,Np + st = SIN(theta(i)) + g11(i,j) = g11(i,j) - sph_g11 + g22(i,j) = g22(i,j) - sph_g22 + g33(i,j) = g33(i,j) - sph_g22*st**2 + dg22(i,j) = dg22(i,j) - sph_dg22 + dg33(i,j) = dg33(i,j) - sph_dg22*st**2 + ENDDO + ENDDO + + + IF (all_modes == 0) THEN + lmin = l ; lmax = l + ELSE + lmin = 2 ; lmax = l + ENDIF + + loop_l: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = m ; mmax = m + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m: DO im = mmin,mmax,mstep + +c ------------------------------------------------------------------ +c +c 5. Calculate all Regge-Wheeler variables +c +c ------------------------------------------------------------------ + +c Factors independent of angular coordinates +c ------------------------------------------ + fac_h1 = dridrsch/DBLE(il*(il+1)) + fac_H2 = S*dridrsch**2 + fac_G = one/(rsch**2*DBLE(il*(il+1)*(il-1)*(il+2))) + fac_K = half/rsch**2 + fac_c1 = dridrsch/DBLE(il*(il+1)) + fac_c2 = two/DBLE(il*(il+1)*(il-1)*(il+2)) + fac_dG = fac_G + fac_dK = fac_K + fac_dc2 = fac_c2 + + +c ------------------------------------------------------------------ +c +c 5a. Real parts of the Regge-Wheeler variables +c +c ------------------------------------------------------------------ + +c Calculate integrands + + loop_phi1: DO j = 1, Np + loop_theta1: DO i = 1, Nt + + st = SIN(theta(i)) + ist = one/st + + CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2,Y3 + & ,Y4) + + h1i(i+1,j+1) = st*g12(i,j)*Y1(1)+ist*g13(i,j)*Y2(1) + H2i(i+1,j+1) = st*g11(i,j)*Y(1) + Gi(i+1,j+1) = (st*g22(i,j)-ist*g33(i,j))*Y3(1) + & +four*ist*g23(i,j)*Y4(1) + Ki(i+1,j+1) = (st*g22(i,j)+ist*g33(i,j))*Y(1) + c1i(i+1,j+1) = g13(i,j)*Y1(1)-g12(i,j)*Y2(1) + c2i(i+1,j+1) = (g22(i,j)-ist**2*g33(i,j))*Y4(1) + & -g23(i,j)*Y3(1) + + g22comb = dridrsch*dg22(i,j)-two/rsch*g22(i,j) + g23comb = dridrsch*dg23(i,j)-two/rsch*g23(i,j) + g33comb = dridrsch*dg33(i,j)-two/rsch*g33(i,j) + + dGi(i+1,j+1) = (st*g22comb-ist*g33comb)*Y3(1) + & +four*ist*g23comb*Y4(1) + dKi(i+1,j+1) = (st*g22comb+ist*g33comb)*Y(1) + dc2i(i+1,j+1) = (dg22(i,j)-ist**2*dg33(i,j))*Y4(1) + & -dg23(i,j)*Y3(1) + + END DO loop_theta1 + END DO loop_phi1 + +c Integrations over the 2-sphere + + h1(1) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp) + H2(1) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp) + G(1) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp) + K(1) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*G(1) + c1(1) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp) + c2(1) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp) + dG(1) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp) + dK(1) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*dG(1) + dc2(1) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp) + + +c ------------------------------------------------------------------ +c +c 5b. Imaginary parts of the Regge-Wheeler variables +c +c ------------------------------------------------------------------ + + complex_parts: IF (igrid /= 1 .AND. igrid /=2) THEN + +c Calculate integrands + + loop_phi2: DO j = 1, Np + + loop_theta2: DO i = 1, Nt + + st = SIN(theta(i)) + ist = one/st + + CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2 + & ,Y3,Y4) + + h1i(i+1,j+1) =-st*g12(i,j)*Y1(2)-ist*g13(i,j)*Y2(2) + H2i(i+1,j+1) = -st*g11(i,j)*Y(2) + Gi(i+1,j+1) = -(st*g22(i,j)-ist*g33(i,j))*Y3(2) + & -four*ist*g23(i,j)*Y4(2) + Ki(i+1,j+1) = -(st*g22(i,j)+ist*g33(i,j))*Y(2) + c1i(i+1,j+1) = -g13(i,j)*Y1(2)+g12(i,j)*Y2(2) + c2i(i+1,j+1) = -(g22(i,j)-ist**2*g33(i,j))*Y4(2) + & +g23(i,j)*Y3(2) + + g22comb = dridrsch*dg22(i,j)-two/rsch*g22(i,j) + g23comb = dridrsch*dg23(i,j)-two/rsch*g23(i,j) + g33comb = dridrsch*dg33(i,j)-two/rsch*g33(i,j) + + dGi(i+1,j+1) = -(st*g22comb-ist*g33comb)*Y3(2) + & -four*ist*g23comb*Y4(2) + dKi(i+1,j+1) = -(st*g22comb+ist*g33comb)*Y(2) + dc2i(i+1,j+1) = -(dg22(i,j)-ist**2*dg33(i,j))*Y4(2) + & +dg23(i,j)*Y3(2) + + END DO loop_theta2 + END DO loop_phi2 + + +c Integrate over 2-sphere + + h1(2) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp) + H2(2) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp) + G(2) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp) + K(2) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*G(2) + c1(2) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp) + c2(2) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp) + dG(2) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp) + dK(2) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp) + & +DBLE(il*(il+1))*half*dG(2) + dc2(2) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp) + + ELSE + + h1(2) = zero + H2(2) = zero + G(2) = zero + K(2) = zero + c1(2) = zero + c2(2) = zero + dG(2) = zero + dK(2) = zero + dc2(2) = zero + + + END IF complex_parts + + + +c ------------------------------------------------------------------ +c +c 6. Construct gauge invariant variables +c +c ------------------------------------------------------------------ + + Qodd(:,il,im) = SQRT(two*DBLE((il+2)*(il+1)*il*(il-1)))*S/ + & rsch*(c1+half*(dc2-two/rsch*c2)) + + Lambda = DBLE((il-1)*(il+2))+3.0D0*(one-S) + + Qeven(:,il,im) = one/Lambda*SQRT(two*DBLE((il-1)*(il+2))/ + & DBLE(il*(il+1)))*(DBLE(il*(il+1))*S*(rsch**2*dG- + & two*h1)+two*rsch*S*(H2-rsch*dK) + & +Lambda*rsch*K) + + + END DO loop_m + + END DO loop_l + + + +c ------------------------------------------------------------------ +c +c 7. Write some diagnostics if required (for last l,m to be used) +c +c ------------------------------------------------------------------ + + IF (debug) THEN + WRITE(101,*) eta,h1 + WRITE(102,*) eta,H2 + WRITE(103,*) eta,G + WRITE(104,*) eta,K + WRITE(105,*) eta,c1 + WRITE(106,*) eta,c2 + WRITE(107,*) eta,dG + WRITE(108,*) eta,dK + WRITE(109,*) eta,dc2 + WRITE(201,*) eta,mass + WRITE(202,*) eta,rsch + WRITE(203,*) eta,drschdri + WRITE(301,*) eta,Qeven(:,lmax,mmax) + WRITE(302,*) eta,Qodd(:,lmax,mmax) + ENDIF + + + END SUBROUTINE D2_extract + +c =============================================================== +c +c +c +c +c +c =============================================================== + + SUBROUTINE simpsons(n,Dx,f,intf) + +c Simpsons rule to evaluate 1D integral equation + +c Variables In: +c n number of points [must be odd] +c Dx step size +c f(n) integrand at each abscissa + +c Variables Out: +c intf evaluated integral + +c Variables Local: +c i index + +c --------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & n + CCTK_REAL,INTENT(IN) :: + & Dx,f(n) + +c Output variables + CCTK_REAL,INTENT(OUT) :: + & intf + +c Local variables + INTEGER :: + & i + CCTK_REAL,PARAMETER :: + & two = 2.0D0, + & three = 3.0D0, + & four = 4.0D0 + +c --------------------------------------------------------------- + + IF (MOD(n,2) == 0) THEN + WRITE(*,*) "Call to -simpsons- with even number of points" + STOP + END IF + + intf = f(1) + four*f(n-1) + f(n) + + DO i = 2, n-3,2 + + intf = intf + four*f(i) + two*f(i+1) + + ENDDO + + intf = intf*Dx/three + + END SUBROUTINE simpsons + +c =============================================================== +c +c +c +c +c +c =============================================================== + + FUNCTION sphere_int(temp,igrid,Nt,Np,Dt,Dp) + +c --------------------------------------------------------------- +c +c Integration the function held in temp over the sphere. +c If the full grid is being used then Simpsons rule is used, +c if an octant is used then the symmetry means Simpsons rule +c becomes a simple sum. +c +c --------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Nt,Np,igrid + CCTK_REAL,INTENT(IN) :: + & temp(Nt,Np),Dt,Dp + +c Output variables + CCTK_REAL :: + & sphere_int + +c Local variables + INTEGER :: + & i,j,Np_start + CCTK_REAL :: + & ft(Nt),fp(Np) + CCTK_REAL,PARAMETER :: + & two = 2.0D0, + & four = 4.0D0 + +c --------------------------------------------------------------- + + DO j = 2, Np + +c Integrate with respect to theta + + DO i = 2, Nt + ft(i) = temp(i,j) + ENDDO + ft(1) = ft(2) ! from symmetry across axis + + SELECT CASE (igrid) + CASE (0) + CALL simpsons(Nt,Dt,ft,fp(j)) + CASE (1) + fp(j)=two*Dt*SUM(ft(2:Nt)) + CASE (2) + CALL simpsons(Nt,Dt,ft,fp(j)) + CASE (3) + fp(j)=two*Dt*SUM(ft(2:Nt)) + END SELECT + + ENDDO + +c Integrate with respect to phi + + fp(1) = fp(2) ! from symmetry across axis + + + SELECT CASE (igrid) + CASE (0) + CALL simpsons(Np,Dp,fp,sphere_int) + CASE (1) + sphere_int=four*Dp*SUM(fp(2:Np)) + CASE (2) + sphere_int=Dp*fp(1) + CASE (3) + CALL simpsons(Np,Dp,fp,sphere_int) + END SELECT + + + END FUNCTION sphere_int + +c ================================================================== +c +c +c +c +c +c ================================================================== + + SUBROUTINE spherical_harmonic(l,m,theta,phi,Ylm) + +c ------------------------------------------------------------------ +c +c Calculate the (l,m) spherical harmonic at given angular +c coordinates. This number is in general complex, and +c +c Ylm = Ylm(1) + i Ylm(2) +c +c where +c +c a ( 2 l + 1 (l-|m|)! ) i m phi +c Ylm = (-1) SQRT( ------- ------- ) P_l|m| (cos(theta)) e +c ( 4 Pi (l+|m|)! ) +c +c and where +c +c a = m/2 (sign(m)+1) +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER :: + & l,m + CCTK_REAL :: + & theta,phi + +c Output variables + CCTK_REAL :: + & Ylm(2) + +c Local variables + INTEGER :: + & i + CCTK_REAL,PARAMETER :: + & one = 1.0D0, + & two = 2.0D0, + & four = 4.0D0 + CCTK_REAL :: + & a,Pi,fac,plgndr + +c ------------------------------------------------------------------ + + Pi = ACOS(-one) + + fac = one + DO i = l-ABS(m)+1,l+ABS(m) + fac = fac*DBLE(i) + ENDDO + fac = one/fac + +c a = (-one)**((m*ISIGN(m,1)/ABS(m)+m)/2)*SQRT(DBLE(2*l+1) +c & /four/Pi*fac)*plgndr(l,ABS(m),cos(theta)) + + a = (-one)**MAX(m,0)*SQRT(DBLE(2*l+1) + & /four/Pi*fac)*plgndr(l,ABS(m),cos(theta)) + Ylm(1) = a*COS(DBLE(m)*phi) + Ylm(2) = a*SIN(DBLE(m)*phi) + + END SUBROUTINE spherical_harmonic + +c ================================================================== +c +c +c +c +c +c ================================================================== + + FUNCTION plgndr(l,m,x) + +c ------------------------------------------------------------------ +c +c From Numerical Recipes +c +c Calculates the associated Legendre polynomial Plm(x). +c Here m and l are integers satisfying 0 <= m <= l, +c while x lies in the range -1 <= x <= 1 +c +c ------------------------------------------------------------------ + +c Input variables + INTEGER,INTENT(IN) :: + & l,m + CCTK_REAL,INTENT(IN) :: + & x + +c Output variables + CCTK_REAL :: + & plgndr + +c Local Variables + INTEGER :: + & i,ll + CCTK_REAL,PARAMETER :: + & one = 1.0D0, + & two = 2.0D0 + CCTK_REAL :: + & pmm,somx2,fact,pmmp1,pll + +c ------------------------------------------------------------------ + + pmm = one + + IF (m.GT.0) THEN + somx2=SQRT((one-x)*(one+x)) + fact=one + DO i=1,m + pmm = -pmm*fact*somx2 + fact = fact+two + END DO + ENDIF + + IF (l.EQ.m) THEN + plgndr = pmm + ELSE + pmmp1 = x*(two*m+one)*pmm + IF (l.EQ.m+1) THEN + plgndr=pmmp1 + ELSE + DO ll=m+2,l + pll = (x*DBLE(2*ll-1)*pmmp1- + & DBLE(ll+m-1)*pmm)/DBLE(ll-m) + pmm = pmmp1 + pmmp1 = pll + END DO + plgndr = pll + ENDIF + ENDIF + + END FUNCTION plgndr + +c ================================================================== +c +c +c +c +c +c ================================================================== + + SUBROUTINE spher_harm_combs(theta,phi,l,m,Y,Y1,Y2,Y3,Y4) +c ------------------------------------------------------------------ +c +c Calculates the various combinations of spherical harmonics needed +c for the extraction (all are complex): +c +c Y = Ylm +c Y1 = Ylm,theta +c Y2 = Ylm,phi +c Y3 = Ylm,theta,theta-cot theta Ylm,theta-Ylm,phi,phi/sin^2 theta +c Y4 = Ylm,theta,phi-cot theta Ylm,phi +c +c The local variables Yplus is the spherical harmonic at (l+1,m) +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER :: + & l,m + CCTK_REAL :: + & theta,phi + +c Output variables + CCTK_REAL,DIMENSION(2) :: + & Y,Y1,Y2,Y3,Y4 + +c Local variables + INTEGER :: + & i + CCTK_REAL :: + & Yplus(2),rl,rm,cot_theta,Pi + CCTK_REAL,PARAMETER :: + & half = 0.5D0, + & one = 1.0D0, + & two = 2.0D0 + +c ------------------------------------------------------------------ + + Pi = DACOS(-one) + + rl = DBLE(l) + rm = DBLE(m) + + cot_theta = DCOS(theta)/DSIN(theta) + + CALL spherical_harmonic(l+1,m,theta,phi,Yplus) + + +c Find Y ... + + CALL spherical_harmonic(l,m,theta,phi,Y) + + +c Find Y1 ... + + DO i = 1,2 + Y1(i) = -(rl+one)*cot_theta*Y(i)+Yplus(i)/DSIN(theta)* + & DSQRT(((rl+one)**2-rm**2)*(rl+half)/(rl+one+half)) + ENDDO + + +c Find Y2 ... + + Y2(1) = -rm*Y(2) + Y2(2) = rm*Y(1) + + +c Find Y3 ... + + DO i = 1,2 + Y3(i) = -two*cot_theta*Y1(i)+(two*rm*rm/(DSIN(theta)**2) + & -rl*(rl+one))*Y(i) + ENDDO + + +c Find Y4 ... + + Y4(1) = rm*(cot_theta*Y(2)-Y1(2)) + Y4(2) = rm*(Y1(1)-cot_theta*Y(1)) + + + END SUBROUTINE spher_harm_combs + +c ================================================================== + + + + diff --git a/src/D2_extract_int.F b/src/D2_extract_int.F new file mode 100644 index 0000000..43f6906 --- /dev/null +++ b/src/D2_extract_int.F @@ -0,0 +1,42 @@ + +#include "cctk.h" + + MODULE D2_extract_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11, + & g12,g13,g22,g23,g33,dg22,dg23,dg33, + & do_ADMmass,ADMmass_int1,ADMmass_int2, + & do_momentum,momentum_int1,momentum_int2,momentum_int3, + & do_spin,spin_int1,spin_int2,spin_int3, + & mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt) + + IMPLICIT NONE + INTEGER,INTENT(IN) :: + & igrid,l,m + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & eta + CCTK_REAL,INTENT(IN),DIMENSION (:) :: + & theta,phi + CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: + & g00,g11,g12,g13,g22,g23,g33,dg22,dg23,dg33, + & ADMmass_int1, ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + CCTK_REAL,INTENT(OUT) :: + & ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt, + & momentum(3),spin(3) + END SUBROUTINE + + END INTERFACE + + END MODULE D2_extract_int + + diff --git a/src/D3_extract.F b/src/D3_extract.F new file mode 100644 index 0000000..f5f6850 --- /dev/null +++ b/src/D3_extract.F @@ -0,0 +1,221 @@ + +#include "cctk.h" + +c ================================================================== + + SUBROUTINE D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,igrid, + & origin,myproc,Nt,Np,all_modes, + & l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven,Extract_temp3d,dtaudt) + +c ------------------------------------------------------------------ +c +c Extract odd and even parity gauge invariant variables on a +c 2-surface defined by constant isotropic radius ETA from ORIGIN, +c also returns the corresponding Schwarzschild radius RSCH. +c +c Works for a full grid (IGRID=0) or an octant (IGRID=1), given +c non-conformal (PSI=1,Gij) or conformal (PSI^PSI_POWER,Gij). +c +c IMPORTANT : Nt and Np must both be even numbers (so that Simpsons +c rule works later) +c +c Variables in: +c ------------ +c igrid Grid type (0=full,1=octant) +c [Could extract in an octant from a full grid, but why +c would you want to] +c origin The position of the origin of extraction in the x,y,z +c coordinate system +c myproc Cactus variable, must be zero if not using Cactus +c Nt,Np Number of theta,phi grid DIVISIONS to use [even numbers] +c all_modes if true the extract all l,m modes up to l +c l,m Spherical harmonic to extract if ALL_MODES is false. +c otherwise m is ignored and l is the maximum l mode +c x,y,z Cartesian coordinates in the 3-space, hopefully +c isotropic about ORIGIN +c Dx,Dy,Dz Grid spacings +c Psi_power The power of Psi which is being passed in, usually one +c but may be four +c Psi The conformal factor, set it to one if non-conformal +c data is passed in +c gij The 3-metric components in cartesian coordinates +c eta The coordinate radius from ORIGIN for extraction, must +c obviously give a 2-surface which is contained in the +c grid +c +c Variables out: +c ------------- +c ADMmass Estimate of ADM mass +c momentum Estimate of momentum +c spin Estimate of spin +c mass The extracted mass +c rsch The extracted Schwarzschild radius +c Qodd The extracted odd parity gauge invariant variable +c Qeven The extracted even parity gauge invariant variable +c +c ------------------------------------------------------------------ + + USE D3_to_D2_int + USE cartesian_to_spherical_int + USE unphysical_to_physical_int + USE D2_extract_int + +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + + CCTK_POINTER :: cctkGH + + INTEGER,INTENT(IN) :: + & igrid,l,m,Psi_power,myproc + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & x,y,z + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_temp3d + +c Output variables + CCTK_REAL,INTENT(OUT) :: + & ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt, + & momentum(3),spin(3) + +c Local variables passed on + CCTK_REAL :: + & theta(Nt),phi(Np) + CCTK_REAL,DIMENSION(Nt,Np) :: + & Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs, + & dgyys,dgyzs,dgzzs,grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp, + & ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + +c Local variables here only + INTEGER :: + & i,j + CCTK_REAL,PARAMETER :: + & half = 0.5D0, + & one = 1.0D0, + & two = 2.0D0 + CCTK_REAL :: + & Pi,Dt,Dp + +c ------------------------------------------------------------------ + + +c ------------------------------------------------------------------ +c +c 1. Specify polar coordinates on chosen 2-surface +c +c ------------------------------------------------------------------ + + Pi = two*ASIN(one) + +c Set polar gridspacing + + SELECT CASE (igrid) + CASE (0) ! full grid + Dt = Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + CASE (1) ! octant grid + Dt = half*Pi/DBLE(Nt) + Dp = half*Pi/DBLE(Np) + CASE (2) ! cartoon + Dt = Pi/DBLE(Nt) + Dp = 2D0*Pi + CASE (3) ! bitant + Dt = half*Pi/DBLE(Nt) + Dp = two*Pi/DBLE(Np) + CASE DEFAULT + WRITE (*,*) "Grid incorrectly set in D3_extract" + STOP + END SELECT + + +c Set coordinates + + DO i = 1, Nt + theta(i) = (DBLE(i) - half)* Dt + ENDDO + DO j = 1, Np + phi(j) = (DBLE(j) - half)* Dp + ENDDO + + IF (igrid == 2) phi(1) = 0D0 + + +c ------------------------------------------------------------------ +c +c 2. Project quantities onto the 2-surface +c +c ------------------------------------------------------------------ + + CALL D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin, + & Psi_power,origin,myproc,Dx,Dy,Dz,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys,gxzs, + & gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs, + & ADMmass_int1,ADMmass_int2,momentum_int1,momentum_int2, + & momentum_int3,spin_int1,spin_int2,spin_int3,Extract_temp3d) + +c ------------------------------------------------------------------ +c Now this is done the grid on processor 0 has correct values and the +c grid on all other processors have junk. This means that processor +c 0 needs to calculate the wave forms and the other processors have +c to sit and do nothing. So we shall put this in a conditional +c to only happen on processor 0. +c PW Jan 31 98 + + if (myproc .eq. 0) then +c ------------------------------------------------------------------ +c +c 3. Transform tensor components to polar coordinates on 2-surface +c +c ------------------------------------------------------------------ + + CALL cartesian_to_spherical(theta,phi,eta,gxxs,gxys,gxzs,gyys, + & gyzs,gzzs,grr,grt,grp,gtt,gtp,gpp,dgxxs,dgxys,dgxzs,dgyys, + & dgyzs,dgzzs,dgtt,dgtp,dgpp) + + +c ------------------------------------------------------------------ +c +c 4. Calculate physical quantities on 2-surface +c +c ------------------------------------------------------------------ + + CALL unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp, + & dgpp,Psis,dPsis,Psi_power) + + +c ------------------------------------------------------------------ +c +c 5. Extract gauge invariant quantities on 2-surface +c +c ------------------------------------------------------------------ + + CALL D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00s,grr, + & grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp, + & do_ADMmass,ADMmass_int1,ADMmass_int2, + & do_momentum,momentum_int1,momentum_int2,momentum_int3, + & do_spin,spin_int1,spin_int2,spin_int3, + & mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt) + +c End of myproc = 0 conditional + endif + + END SUBROUTINE D3_extract + + + diff --git a/src/D3_extract_int.F b/src/D3_extract_int.F new file mode 100644 index 0000000..3bb617a --- /dev/null +++ b/src/D3_extract_int.F @@ -0,0 +1,44 @@ + +#include "cctk.h" + + MODULE D3_extract_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,igrid, + & origin,myproc, + & Nt,Np,all_modes,l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven, + & Extract_temp3d,dtaudt) + + IMPLICIT NONE + + CCTK_POINTER :: cctkGH + + INTEGER,INTENT(IN) :: + & igrid,l,m,Psi_power,myproc + CCTK_INT,INTENT(IN) :: + & Nt,Np,all_modes,do_momentum,do_spin + INTEGER,INTENT(IN) :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & x,y,z + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_Temp3d + CCTK_REAL,INTENT(OUT) :: + & ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt, + & momentum(3),spin(3) + + END SUBROUTINE + + END INTERFACE + + END MODULE D3_extract_int diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F new file mode 100644 index 0000000..884f8b8 --- /dev/null +++ b/src/D3_to_D2.F @@ -0,0 +1,398 @@ + +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin, + & Psi_power,origin,myproc,Dx,Dy,Dz,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & x,y,z,eta,Nt,Np,theta,phi, + & Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs, + & dgyys,dgyzs,dgzzs,ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3,Extract_temp3d) + +c ------------------------------------------------------------------------ +c +c Project the 3-metric and its 1st radial derivatives onto a +c 2-sphere. +c +c ------------------------------------------------------------------------ + + USE met_rad_der_int + USE ADMmass_integrand3D_int + USE momentum_integrand3D_int + USE spin_integrand3D_int + + IMPLICIT NONE + +c Input variables + + CCTK_POINTER :: cctkGH + + INTEGER,INTENT(IN) :: + & myproc,Psi_power + CCTK_INT, INTENT(IN) :: + & Nt,Np,do_momentum,do_spin + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & theta,phi,x,y,z + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_temp3d + LOGICAL :: + & do_ADMmass(2) + +c Output variables + + CCTK_REAL,INTENT(OUT),DIMENSION(:,:) :: + & Psis,g00s,gxxs,gxys,gxzs,gyys, + & gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs, + & ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + + +c Local variables, passed on + + LOGICAL :: + & err_flag + INTEGER :: + & iorder,npoints,Nx,Ny,Nz + INTEGER,DIMENSION(Nt,Np) :: + & ib,jb,kb + CCTK_REAL,DIMENSION(Nt,Np) :: + & xs,ys,zs,ux,uy,uz,xb,yb,zb + +c Local variables, here only + + INTEGER :: + & i,j,handle,ierror + CCTK_REAL :: + & xmin,xmax,ymin,ymax,zmin,zmax + +c ------------------------------------------------------------------------ + +c Initial Stuff +c ------------- + Nx = SIZE(x) + Ny = SIZE(y) + Nz = SIZE(z) + + +c Compute Cartesian coordinates on the surface +c -------------------------------------------- + DO j = 1, Np + DO i = 1, Nt + + xs(i,j) = origin(1)+eta*SIN(theta(i))*COS(phi(j)) + ys(i,j) = origin(2)+eta*SIN(theta(i))*SIN(phi(j)) + zs(i,j) = origin(3)+eta*COS(theta(i)) + + ENDDO + ENDDO + + +c Only do interpolation on one processor +c -------------------------------------- + SELECT CASE (myproc) + + CASE (0) + npoints = Np*Nt + + CASE DEFAULT + npoints = 0 + + END SELECT + + +c Choose the interpolator +c ----------------------- + call CCTK_GetInterpHandle (handle, "simple") + +c Find the origin of the spatial coordinates +c ------------------------------------------ + call CCTK_CoordRange(ierror,cctkGH,xmin,xmax,"x") + call CCTK_CoordRange(ierror,cctkGH,ymin,ymax,"y") + call CCTK_CoordRange(ierror,cctkGH,zmin,zmax,"z") + + +c Project un-physical metric and conformal factor onto sphere +c ------------------------------------------------------------ + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 8, 8, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL + $ ) + + +c Calculate radial derivatives and project onto sphere +c ---------------------------------------------------- + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dPsis, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgxxs, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgxys, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgxzs, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgyys, + $ CCTK_VARIABLE_REAL + $ ) + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgyzs, + $ CCTK_VARIABLE_REAL + $ ) + + + CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,Extract_temp3d) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ dgzzs, + $ CCTK_VARIABLE_REAL + $ ) + + +c Calculate integrands for ADM masses +c ----------------------------------- +c Standard equation + IF (do_ADMmass(1)) THEN + CALL ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy, + & gxz,gyy,gyz,gzz,Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ ADMmass_int1, + $ CCTK_VARIABLE_REAL + $ ) + + END IF + +c Conformal equation + IF (do_ADMmass(2)) THEN + ADMmass_int2 = -eta**2/2D0/3.1416D0*dPsis + ENDIF + +c Calculate integrands for momentum +c --------------------------------- + IF (do_momentum==1) THEN + + CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ momentum_int1, + $ CCTK_VARIABLE_REAL + $ ) + + CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ momentum_int2, + $ CCTK_VARIABLE_REAL + $ ) + + + CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ momentum_int3, + $ CCTK_VARIABLE_REAL + $ ) + + END IF + +c Calculate integrands for spin +c ----------------------------- + IF (do_momentum==1) THEN + + CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ spin_int1, + $ CCTK_VARIABLE_REAL + $ ) + + CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ spin_int2, + $ CCTK_VARIABLE_REAL + $ ) + + + CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z, + & gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & Extract_temp3d,Psi,Psi_power) + CALL CCTK_SyncGroup(cctkGH,"extract::temps") + call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1, + $ nx, ny, nz, + $ xs, ys, zs, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ xmin, ymin, zmin, + $ dx, dy, dz, + $ Extract_temp3d, + $ CCTK_VARIABLE_REAL, + $ spin_int3, + $ CCTK_VARIABLE_REAL + $ ) + + END IF + + END SUBROUTINE D3_to_D2 + + + + + + + + + + + + + + + + + diff --git a/src/D3_to_D2_int.F b/src/D3_to_D2_int.F new file mode 100644 index 0000000..297e22e --- /dev/null +++ b/src/D3_to_D2_int.F @@ -0,0 +1,45 @@ + +#include "cctk.h" + + MODULE D3_to_D2_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin, + & Psi_power,origin,myproc,Dx,Dy,Dz,Psi, + & g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz, + & x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys, + & gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs, + & dgzzs,ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3,Extract_temp3d) + + IMPLICIT NONE + CCTK_POINTER :: cctkGH + INTEGER,INTENT(IN) :: + & myproc,Psi_power + CCTK_INT, INTENT(IN) :: + & Nt,Np,do_momentum,do_spin + INTEGER :: + & do_ADMmass(2) + CCTK_REAL,INTENT(IN) :: + & origin(3),Dx,Dy,Dz,eta + CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: + & Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,INTENT(IN),DIMENSION(:) :: + & x,y,z,theta,phi + CCTK_REAL,INTENT(OUT),DIMENSION(:,:) :: + & Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs, + & dgyys,dgyzs,dgzzs,ADMmass_int1,ADMmass_int2, + & momentum_int1,momentum_int2,momentum_int3, + & spin_int1,spin_int2,spin_int3 + CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) :: + & Extract_temp3d + END SUBROUTINE + + END INTERFACE + + END MODULE D3_to_D2_int diff --git a/src/Extract.F b/src/Extract.F new file mode 100644 index 0000000..1aa8534 --- /dev/null +++ b/src/Extract.F @@ -0,0 +1,814 @@ +c/*@@ +c @file Extract.F +c @date 28th October 1997 +c @author Gab Allen +c @desc Waveform extraction +c +c @enddesc +c@@*/ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +c/*@@ +c @routine Extract +c @date 28th October 1997 +c @author Gab Allen +c @desc Entry point for waveform extraction routines +c +c @enddesc +c @calls D3_extract +c @calledby No idea +c @history +c +c @endhistory +c@@*/ + + + SUBROUTINE Extract(CCTK_FARGUMENTS) + + USE D3_extract_int + + IMPLICIT NONE + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + +c Non-Cactus input variables for D3_extract + + INTEGER :: + & igrid,lmode,mmode,Psi_power=1 + INTEGER :: + & do_ADMmass(2) + CCTK_REAL :: + & orig(3),radius + CCTK_REAL :: + & x_1d(cctk_lsh(1)),y_1d(cctk_lsh(2)),z_1d(cctk_lsh(3)) + + +c Output variables from D3_extract + + CCTK_REAL :: + & dtaudt,ADMmass(2),mass,rsch,momentum(3),spin(3) + CCTK_REAL,ALLOCATABLE :: + & Qodd(:,:,:), Qeven(:,:,:) + + +c Local variables + + INTEGER :: + & ix,iy,iz,fn1,fn2,lmin,lmax,mmin,mmax,lstep,mstep, + & il,im,it,ndet,idet,nx_parfile,nchar,iconv,ioutput + INTEGER,SAVE :: openfiles = 1 + INTEGER,PARAMETER :: + & max_detectors = 9,number_timecoords = 2 + CCTK_REAL,PARAMETER :: + & two = 2.0D0 + CCTK_REAL :: + & time,r1,r2,dr,r_max,xmin,xmax,ymin,ymax,zmin,zmax,gf_min,gf_max + CCTK_REAL :: + & detector(max_detectors) + CHARACTER*100 :: + & out1,out2,massfile,rschfile,ADMmassfile1,ADMmassfile2, + & momentumfile1,momentumfile2,momentumfile3,spinfile1,spinfile2,spinfile3 + CHARACTER*3 :: + & timestring + + LOGICAL :: use_proper,use_coordinate,do_output + LOGICAL,SAVE :: firsttime + CCTK_REAL,DIMENSION(10),SAVE :: proper_time + CCTK_REAL,SAVE :: last_time + INTEGER, save :: open_file_level(10) = 0 + + INTEGER myproc,ierr + CCTK_REAL dx,dy,dz + character*200 filestr + character*80 infoline + +c ------------------------------------------------------------------ + +c Get the output directory sorted out + call CCTK_FortranString(nchar,outdir,filestr) + if (firsttime) then + call CCTK_mkdir(ierr,filestr); + firsttime = .FALSE. + end if + +c ------------------------------------------------------------------ + + myproc = CCTK_MyProc(cctkGH) + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + +c ------------------------------------------------------------------ +c +c 0. Check to see if Extract should have been called +c +c ------------------------------------------------------------------ + + IF (MODULO(int(cctk_iteration),int(itout)) .NE. 0) RETURN + + if (verbose == 1) then + write(infoline,'(A24,G12.7)') + & 'Calling Extract at time ',cctk_time + call CCTK_INFO(infoline) + end if + +c ------------------------------------------------------------------ +c +c 1. Initial Stuff +c +c ------------------------------------------------------------------ + + +c Get the number of polar divisions + + IF (MOD(int(Nt),2) == 1 .OR. MOD(int(Np),2) == 1 .OR. + & Nt < 0 .OR. Np < 0) THEN + call CCTK_WARN(0,"Error in Nt or Np in Extract") + END IF + + +c Get the origin of spherical symmetry + + orig(1) = origin_x + orig(2) = origin_y + orig(3) = origin_z + + +c Set the value of igrid + + IF (CCTK_EQUALS(domain,"octant")) THEN + igrid = 1 + ELSEIF (CCTK_EQUALS(domain,"bitant")) THEN + igrid = 3 + ELSE + igrid = 0 + END IF + +#ifdef THORN_CARTOON_2D + IF (contains("cartoon_active","yes") .NE. 0) THEN + igrid = 2 + Np = 1 + END IF +#endif + +c Create 1D coordinate arrays + + do ix = 1, cctk_lsh(1) + x_1d(ix) = x(ix,1,1) + enddo + do iy = 1, cctk_lsh(2) + y_1d(iy) = y(1,iy,1) + enddo + do iz = 1, cctk_lsh(3) + z_1d(iz) = z(1,1,iz) + enddo + +c See if the ADM mass should be calculated and output + IF (doADMmass == 1) THEN + do_ADMmass(1) = 1 + IF (use_conformal == 1) THEN + do_ADMmass(2) = 1 + ELSE + do_ADMmass(2) = 0 + ENDIF + ELSE + do_ADMmass(:) = 0 + END IF + +c Set array containing gtt ... for the moment I am lazy and +c just make it the lapse and do not include the shift! + g00 = alp + +c What kind of timecoordinate do I want to use + if (CCTK_EQUALS(timecoord,"both")) then + use_proper = .true. + use_coordinate = .true. + if (cctk_iteration == 0) then + proper_time = 0 + last_time = cctk_time + end if + elseif (CCTK_EQUALS(timecoord,"proper")) then + use_proper = .true. + use_coordinate = .false. + if (cctk_iteration == 0) then + proper_time = 0 + last_time = cctk_time + end if + elseif (CCTK_EQUALS(timecoord,"coord")) then + use_proper = .false. + use_coordinate = .true. + else + use_proper = .false. + use_coordinate = .true. + endif + +c ------------------------------------------------------------------ +c +c 2. Sort out which modes to extract +c +c ------------------------------------------------------------------ + + lmode = l_mode + + IF (all_modes == 1) THEN + lmin = 2 + lmax = lmode +c Get m_mode in this case too, or the T3E has a junk value (PW) + mmode = 0 + ALLOCATE(Qodd(2,2:lmode,0:lmode),Qeven(2,2:lmode,0:lmode)) + ELSE + lmin = lmode + lmax = lmode + mmode = m_mode + ALLOCATE(Qodd(2,2:lmode,0:mmode),Qeven(2,2:lmode,0:mmode)) + ENDIF + + +c Check do not have bad values for the modes + + IF (lmode < 2 .OR. mmode > lmode .OR. mmode < 0) THEN + call CCTK_WARN(0,"Error in lmode,mmode in Extract") + END IF + IF (lmode > 9) THEN + call CCTK_WARN(4,"Filenames will probably be crazy in Extract") + END IF + +c If using an octant, do not do odd modes + + IF (igrid == 1 .OR. igrid == 2) THEN + lstep = 2 ; mstep = 2 + ELSE + lstep = 1 ; mstep = 1 + END IF + +c If we have cartoon, then jump past all the m-modes in loop + IF (igrid == 2) mstep = lmode+1 + + +c ------------------------------------------------------------------ +c +c 3. Find maximum radius of extraction on grid +c +c ------------------------------------------------------------------ + + call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,"x") + call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,"y") + call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z") + + IF (igrid == 2) THEN + r_max = MIN(xmax-two*Dx-orig(1), + & zmax-two*Dz-orig(3)) + ELSE IF (igrid == 1) THEN + r_max = MIN(zmax-two*Dx-orig(1), + & ymax-two*Dy-orig(2), + & zmax-two*Dz-orig(3)) + ELSE IF (igrid == 3) THEN + r_max = MIN(xmax-two*Dx-orig(1), + & ymax-two*Dy-orig(2), + & zmax-two*Dz-orig(3), + & DABS(xmin+two*Dx-orig(1)), + & DABS(ymin+two*Dy-orig(2))) + ELSE + r_max = MIN(xmax-two*Dx-orig(1), + & ymax-two*Dy-orig(2), + & zmax-two*Dz-orig(3), + & DABS(xmin+two*Dx-orig(1)), + & DABS(ymin+two*Dy-orig(2)), + & DABS(zmin+two*Dz-orig(3))) + END IF + + + + +c ------------------------------------------------------------------ +c +c 4. Extract Cauchy initial data for a linear wave equation +c +c ------------------------------------------------------------------ + + test_Cauchy: IF (Cauchy == 1) THEN + + it = Cauchy_timestep + + +c check_timestep: IF (cctk_iteration == it) THEN + + check_timestep: IF (open_file_level(1) .eq. 0) THEN + +c Open output files + write (*,*) "Open output file on level ", 1 + open_file_level(1)=1; + + test_myproc1: IF (CCTK_MyProc(cctkGH) == 0 ) THEN + + out1 = filestr(1:nchar)//"/rsch_ini.rl" + OPEN(UNIT= 77,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/mass_ini.rl" + OPEN(UNIT= 78,file = out1,STATUS="unknown") + IF (do_ADMmass(1) == 1) THEN + out1 = filestr(1:nchar)//"/ADMmass_ini.rl" + OPEN(UNIT= 79,file = out1,STATUS="unknown") + ENDIF + IF (do_ADMmass(2) == 1) THEN + out1 = filestr(1:nchar)//"/ADMmassc_ini.rl" + OPEN(UNIT= 80,file = out1,STATUS="unknown") + ENDIF + + IF (do_momentum == 1) THEN + out1 = filestr(1:nchar)//"/momentum_x_ini.rl" + OPEN(UNIT=781,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/momentum_y_ini.rl" + OPEN(UNIT=782,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/momentum_z_ini.rl" + OPEN(UNIT=783,file = out1,STATUS="unknown") + ENDIF + IF (do_spin == 1) THEN + out1 = filestr(1:nchar)//"/spin_x_ini.rl" + OPEN(UNIT=784,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/spin_y_ini.rl" + OPEN(UNIT=785,file = out1,STATUS="unknown") + out1 = filestr(1:nchar)//"/spin_z_ini.rl" + OPEN(UNIT=786,file = out1,STATUS="unknown") + ENDIF + + loop_l1: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m1: DO im = mmin,mmax,lstep + + + fn1 = il*10+im + out1 = filestr(1:nchar)//"/Qeven_ini_"// + & CHAR(il+48)//CHAR(im+48)//".rl" + OPEN(UNIT=fn1,FILE=out1,STATUS="unknown") + +c Only print odd-parity if full grid + IF (igrid == 0) THEN + fn2 = 100+il*10+im + out2 = filestr(1:nchar)//"/Qodd_ini_"// + & CHAR(il+48)//CHAR(im+48)//".rl" + OPEN(UNIT=fn2,FILE=out2,STATUS="unknown") + END IF + + END DO loop_m1 + + END DO loop_l1 + + END IF test_myproc1 + + +c Find range of extraction radii + + r1 = Cauchy_r1 + r2 = r_max + dr = Cauchy_dr + + + +c Do extraction at each radius + + IF (verbose == 1) THEN + WRITE(*,*) + WRITE(*,*) "Extracting Cauchy initial data" + WRITE(*,*) "------------------------------" + WRITE(*,*) "From r =",r1 + WRITE(*,*) "To r =",r2 + WRITE(*,*) "In steps of ",dr + WRITE(*,*) + END IF + + radius = r1 + + extract_at_each_radius: DO WHILE (radius < r2) + + CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin, + & igrid,orig,myproc,Nt,Np,all_modes,lmode, + & mmode,x_1d,y_1d,z_1d,Dx,Dy,Dz,Psi_power,Psi,g00, + & gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz, + & radius,ADMmass,momentum,spin,mass,rsch, + & Qodd,Qeven,temp3d,dtaudt) + + IF (verbose == 1) THEN + WRITE(*,*) "Extracted at r =",radius + WRITE(*,*) "Schwarzschild radius =",rsch + WRITE(*,*) "Schwarzschild mass =",mass + if (do_ADMmass(1) == 1) then + WRITE(*,*) "ADM mass =",ADMmass(1) + end if + if (do_ADMmass(2) == 1) then + WRITE(*,*) "ADM mass (conformal) =",ADMmass(2) + end if + if (do_momentum == 1) then + WRITE(*,*) "momentum = (",momentum(1),",", + & momentum(2),",",momentum(3),")" + end if + if (do_spin == 1) then + WRITE(*,*) "spin = (",spin(1),",", + & spin(2),",",spin(3),")" + end if + ENDIF + + +c Write to file at each radius + + test_myproc2: IF (CCTK_MyProc(cctkGH) == 0) THEN + + WRITE( 77,*) radius,rsch + WRITE( 78,*) radius,mass + IF (do_ADMmass(1) == 1) WRITE( 79,*) radius,ADMmass(1) + IF (do_ADMmass(2) == 1) WRITE( 80,*) radius,ADMmass(2) + IF (do_momentum == 1) then + WRITE(781,*) radius,momentum(1) + WRITE(782,*) radius,momentum(2) + WRITE(783,*) radius,momentum(3) + end if + + IF (do_spin == 1) THEN + WRITE(784,*) radius,spin(1) + WRITE(785,*) radius,spin(2) + WRITE(786,*) radius,spin(3) + END IF + + loop_l2: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m2: DO im = mmin,mmax,mstep + + fn1 = il*10+im + WRITE(fn1,*) rsch,Qeven(:,il,im) + +c Only print odd-parity if full grid + IF (igrid == 0) THEN + fn2 = 100+il*10+im + WRITE(fn2,*) rsch,Qeven(:,il,im) + END IF + + END DO loop_m2 + + END DO loop_l2 + + END IF test_myproc2 + + radius = radius+dr + + END DO extract_at_each_radius + + +c Now close all the files + + test_myproc3: IF (myproc == 0) THEN + + CLOSE( 77) ! Schwarzschild radius + CLOSE( 78) ! Schwarzschild mass + IF (do_ADMmass(1) == 1) CLOSE( 79) ! ADM mass + IF (do_ADMmass(2) == 1) CLOSE( 80) ! ADM mass + + IF (do_momentum == 1) THEN + CLOSE(781) ! momentum + CLOSE(782) ! momentum + CLOSE(783) ! momentum + end if + + IF (do_spin == 1) THEN + CLOSE(784) ! spin + CLOSE(785) ! spin + CLOSE(786) ! spin + END IF + + loop_l3: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m3: DO im = mmin,mmax,mstep + + fn1 = il*10+im + CLOSE(fn1) ! Qeven_ini_lm.dat + + IF (igrid == 0) THEN + fn2 = 100+il*10+im + CLOSE(fn2) ! Qodd_ini_lm.dat + END IF + + END DO loop_m3 + + END DO loop_l3 + + END IF test_myproc3 + + END IF check_timestep + + END IF test_Cauchy + + +c ------------------------------------------------------------------ +c +c 5. Extract waveforms at detectors +c +c ------------------------------------------------------------------ + +c Cannot use the conformal equation for ADM mass now + do_ADMmass(2) = 0 + + ndet = num_detectors + + test_detectors: IF (ndet > 0) THEN + + IF (ndet > 9) THEN + call CCTK_WARN(0,"Too many detectors in Extract") + END IF + + detector(1) = detector1 + detector(2) = detector2 + detector(3) = detector3 + detector(4) = detector4 + detector(5) = detector5 + detector(6) = detector6 + detector(7) = detector7 + detector(8) = detector8 + detector(9) = detector9 + + DO idet = ndet+1, max_detectors + detector(idet) = 0.0D0 + END DO + + DO idet = 1, ndet + IF (detector(idet) > r_max) THEN + IF (openfiles == 1 .OR. cctk_iteration == it) THEN + call CCTK_WARN(8,"Removing detectors outside coordinate range") +c IF (iconv == 0) STOP + END IF + ndet =idet-1 + GOTO 404 + ENDIF + END DO + 404 CONTINUE + + IF (verbose == 1) THEN + IF (openfiles == 1) THEN + WRITE(*,*) + WRITE(*,*) "Extracting at detectors" + WRITE(*,*) "-----------------------" + WRITE(*,*) "Using ",ndet," detectors at" + DO idet=1,ndet + WRITE(*,*) "r = ",detector(idet) + ENDDO + END IF + END IF + + detector_loop: DO idet = 1,ndet + + radius = detector(idet) + + IF (verbose == 1) THEN + WRITE(*,*) "Calling extract at radius ... ",radius + END IF + + CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin, + & igrid,orig,myproc,Nt, + & Np,all_modes,lmode,mmode,x_1d,y_1d,z_1d, + & Dx,Dy,Dz,Psi_power,Psi,g00, + & gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz, + & radius,ADMmass,momentum,spin,mass,rsch, + & Qodd,Qeven,temp3d,dtaudt) + + test_verbose: + & IF (verbose == 1) THEN + + WRITE(*,*) + WRITE(*,*) "Detector ",idet," ..." + WRITE(*,*) "Schwarzschild radius =",rsch + WRITE(*,*) "Schwarzschild mass =",mass + + loop_l4: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m4: DO im = mmin,mmax,mstep + + WRITE(*,*) filestr(1:nchar)//"/Qeven_"// + & CHAR(il+48)//CHAR(im+48),Qeven(:,il,im) + WRITE(*,*) filestr(1:nchar)//"/Qodd_"// + & CHAR(il+48)//CHAR(im+48),Qodd(:,il,im) + + END DO loop_m4 + + END DO loop_l4 + + END IF test_verbose + +c Output to files + + test_myproc4: IF (myproc == 0) THEN + + DO ioutput=1,number_timecoords + + do_output = .false. + timestring = "ERR" + time = 0D0 + IF (ioutput == 1 .AND. use_coordinate) THEN + timestring = ".tl" + do_output = .true. + time = cctk_time + ELSEIF (ioutput == 2 .AND. use_proper) THEN + timestring = ".ul" + do_output = .true. + time = proper_time(idet) + (cctk_time-last_time)*dtaudt + proper_time(idet) = time + if (idet == ndet) last_time = cctk_time + ENDIF + + IF (do_output) THEN + +c Output extracted radius and mass + rschfile = filestr(1:nchar)//"/rsch_"//"R"// + & CHAR(idet+48)//timestring + massfile = filestr(1:nchar)//"/mass_"//"R"// + & CHAR(idet+48)//timestring + ADMmassfile1 = filestr(1:nchar)//"/ADMmass_"//"R"// + & CHAR(idet+48)//timestring + ADMmassfile2 = filestr(1:nchar)//"/ADMmassc_"//"R"// + & CHAR(idet+48)//timestring + momentumfile1 = filestr(1:nchar)//"/momentum_x_"//"R"// + & CHAR(idet+48)//timestring + momentumfile2 = filestr(1:nchar)//"/momentum_y_"//"R"// + & CHAR(idet+48)//timestring + momentumfile3 = filestr(1:nchar)//"/momentum_z_"//"R"// + & CHAR(idet+48)//timestring + spinfile1 = filestr(1:nchar)//"/spin_x_"//"R"// + & CHAR(idet+48)//timestring + spinfile2 = filestr(1:nchar)//"/spin_y_"//"R"// + & CHAR(idet+48)//timestring + spinfile3 = filestr(1:nchar)//"/spin_z_"//"R"// + & CHAR(idet+48)//timestring + + IF (openfiles == 1) THEN + OPEN(UNIT= 77,FILE=rschfile,STATUS="unknown") + OPEN(UNIT= 78,FILE=massfile,STATUS="unknown") + IF (do_ADMmass(1) == 1) + & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="unknown") + IF (do_ADMmass(2) == 1) + & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="unknown") + IF (do_momentum == 1) THEN + OPEN(UNIT=781,FILE=momentumfile1,STATUS="unknown") + OPEN(UNIT=782,FILE=momentumfile2,STATUS="unknown") + OPEN(UNIT=783,FILE=momentumfile3,STATUS="unknown") + END IF + + IF (do_spin == 1) THEN + OPEN(UNIT=784,FILE=spinfile1,STATUS="unknown") + OPEN(UNIT=785,FILE=spinfile2,STATUS="unknown") + OPEN(UNIT=786,FILE=spinfile3,STATUS="unknown") + END IF + ELSE + OPEN(UNIT= 77,FILE=rschfile,STATUS="old", + & POSITION="append") + OPEN(UNIT= 78,FILE=massfile,STATUS="old", + & POSITION="append") + IF (do_ADMmass(1) == 1) + & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="old", + & POSITION="append") + IF (do_ADMmass(2) == 1) + & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="old", + & POSITION="append") + IF (do_momentum == 1) THEN + OPEN(UNIT=781,FILE=momentumfile1,STATUS="old", + & POSITION="append") + OPEN(UNIT=782,FILE=momentumfile2,STATUS="old", + & POSITION="append") + OPEN(UNIT=783,FILE=momentumfile3,STATUS="old", + & POSITION="append") + END IF + + IF (do_spin == 1)THEN + OPEN(UNIT=784,FILE=spinfile1,STATUS="old", + & POSITION="append") + OPEN(UNIT=785,FILE=spinfile2,STATUS="old", + & POSITION="append") + OPEN(UNIT=786,FILE=spinfile3,STATUS="old", + & POSITION="append") + END IF + + ENDIF + WRITE( 77,21) time,rsch + WRITE( 78,21) time,mass + IF (do_ADMmass(1) == 1) WRITE( 79,21) time,ADMmass(1) + IF (do_ADMmass(2) == 1) WRITE( 80,21) time,ADMmass(2) + IF (do_momentum == 1) THEN + WRITE(781,*) cctk_time,momentum(1) + WRITE(782,*) cctk_time,momentum(2) + WRITE(783,*) cctk_time,momentum(3) + END IF + IF (do_spin == 1) THEN + WRITE(784,*) cctk_time,spin(1) + WRITE(785,*) cctk_time,spin(2) + WRITE(786,*) cctk_time,spin(3) + END IF + + CLOSE( 77) + CLOSE( 78) + IF (do_ADMmass(1) == 1) CLOSE( 79) + IF (do_ADMmass(2) == 1) CLOSE( 80) + + IF (do_momentum == 1) THEN + CLOSE(781) + CLOSE(782) + CLOSE(783) + END IF + IF (do_spin == 1) THEN + CLOSE(782) + CLOSE(783) + CLOSE(784) + END IF + + +c Output gauge invariant variables + + loop_l5: DO il = lmin,lmax,lstep + + IF (all_modes == 0) THEN + mmin = mmode ; mmax = mmode + ELSE + mmin = 0 ; mmax = il + END IF + + loop_m5: DO im = mmin,mmax,mstep + + fn1 = 11 + fn2 = 12 + + out1 = filestr(1:nchar)//"/Qeven_"//"R" + & //CHAR(idet+48)//"_"//CHAR(il+48)// + & CHAR(im+48)//timestring + out2 = filestr(1:nchar)//"/Qodd_"//"R" + & //CHAR(idet+48)//"_"//CHAR(il+48)// + & CHAR(im+48)//timestring + +c Write even parity waveforms + IF (openfiles == 1) THEN + OPEN(UNIT=fn1,FILE=out1,STATUS="unknown") + ELSE + OPEN(UNIT=fn1,FILE=out1,STATUS="old", + & POSITION="append") + ENDIF + + WRITE(fn1,20) time,Qeven(:,il,im) + + CLOSE(fn1) + +c Only write odd parity waveforms if full grid + IF (igrid == 0) THEN + IF (openfiles == 1) THEN + OPEN(UNIT=fn2,FILE=out2,STATUS="unknown") + ELSE + OPEN(UNIT=fn2,FILE=out2,STATUS="old", + & POSITION="append") + ENDIF + + WRITE(fn2,20) time,Qodd(:,il,im) + + CLOSE(fn2) + END IF + + END DO loop_m5 + + END DO loop_l5 + + ENDIF + END DO + END IF test_myproc4 + + END DO detector_loop + + END IF test_detectors + + DEALLOCATE(Qodd,Qeven) + + openfiles = 0 + + 20 format(e24.15,e24.15,e24.15) + 21 format(e24.15,e24.15) + + END SUBROUTINE Extract diff --git a/src/cartesian_to_spherical.F b/src/cartesian_to_spherical.F new file mode 100644 index 0000000..d55f84f --- /dev/null +++ b/src/cartesian_to_spherical.F @@ -0,0 +1,116 @@ + +#include "cctk.h" + +c ================================================================== + + SUBROUTINE cartesian_to_spherical(theta,phi,r,gxxs,gxys,gxzs, + & gyys,gyzs, gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys, + & dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps) + +c ------------------------------------------------------------------ +c +c Convert components of a symmetric 2nd rank 3D tensor on a 2-sphere +c from Cartesian to spherical polar coordinates +c +c (x,y,z) ----> (r,t,p) +c +c ( gxx gxy gxz ) ( grr grt grp ) +c ( . gyy gyz ) ----> ( . gpp gtp ) +c ( . . gzz ) ( . . gpp ) +c +c Also convert radial derivatives of the tensor components +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + CCTK_REAL,INTENT(IN) :: + & r + CCTK_REAL,INTENT(IN),DIMENSION (:) :: + & theta,phi + CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: + & gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys,dgyzs, + & dgzzs + +c Output variables + CCTK_REAL,INTENT(OUT),DIMENSION (:,:) :: + & grrs,grts,grps,gtts,gtps, + & gpps,dgtts,dgtps,dgpps + +c Local variables + INTEGER :: + & it,ip,Nt,Np + CCTK_REAL :: + & ct,st,cp,sp,ct2,st2,cp2,sp2,r2,gxx,gxy,gxz,gyy,gyz,gzz,dgxx, + & dgxy,dgxz,dgyy,dgyz,dgzz + CCTK_REAL,PARAMETER :: + & two = 2.0D0 + +c ------------------------------------------------------------------ + + Nt = SIZE(theta) + Np = SIZE(phi) + + r2 = r*r + + DO it = 1, Nt + + ct = COS(theta(it)) ; ct2 = ct*ct + st = SIN(theta(it)) ; st2 = st*st + + DO ip = 1, Np + + cp = COS(phi(ip)) ; cp2 = cp*cp + sp = SIN(phi(ip)) ; sp2 = sp*sp + + gxx = gxxs(it,ip) ; dgxx = dgxxs(it,ip) + gxy = gxys(it,ip) ; dgxy = dgxys(it,ip) + gxz = gxzs(it,ip) ; dgxz = dgxzs(it,ip) + gyy = gyys(it,ip) ; dgyy = dgyys(it,ip) + gyz = gyzs(it,ip) ; dgyz = dgyzs(it,ip) + gzz = gzzs(it,ip) ; dgzz = dgzzs(it,ip) + + grrs(it,ip) = st2*cp2*gxx+st2*sp2*gyy+ct2*gzz + & +two*st2*cp*sp*gxy+two*st*cp*ct*gxz+two*st*ct*sp*gyz + + grts(it,ip) = r*(st*cp2*ct*gxx+two*st*ct*sp*cp*gxy + & +cp*(ct2-st2)*gxz+st*sp2*ct*gyy+sp*(ct2-st2)*gyz-ct*st*gzz) + + grps(it,ip) = r*st*(-st*sp*cp*gxx-st*(sp2-cp2)*gxy-sp*ct*gxz + & +st*sp*cp*gyy+ct*cp*gyz) + + gtts(it,ip) = r2*(ct2*cp2*gxx+two*ct2*sp*cp*gxy + & -two*st*ct*cp*gxz+ct2*sp2*gyy-two*st*sp*ct*gyz+st2*gzz) + + gtps(it,ip) = r2*st*(-cp*sp*ct*gxx-ct*(sp2-cp2)*gxy + & +st*sp*gxz+cp*sp*ct*gyy-st*cp*gyz) + + gpps(it,ip) = r2*st2*(sp2*gxx-two*cp*sp*gxy+cp2*gyy) + + dgtts(it,ip) = two/r*gtts(it,ip)+r2*(ct2*cp2*dgxx + & +two*ct2*sp*cp*dgxy-two*st*ct*cp*dgxz+ct2*sp2*dgyy + & -two*st*sp*ct*dgyz+st2*dgzz) + + dgtps(it,ip) = two/r*gtps(it,ip)+r2*st*(-cp*sp*ct*dgxx + & -ct*(sp2-cp2)*dgxy+st*sp*dgxz+cp*sp*ct*dgyy-st*cp*dgyz) + + dgpps(it,ip) = two/r*gpps(it,ip) +r2*st2*(sp2*dgxx + & -two*cp*sp*dgxy+cp2*dgyy) + + END DO + + END DO + + + END SUBROUTINE cartesian_to_spherical + + + + + + + + + + diff --git a/src/cartesian_to_spherical_int.F b/src/cartesian_to_spherical_int.F new file mode 100644 index 0000000..a4db3b2 --- /dev/null +++ b/src/cartesian_to_spherical_int.F @@ -0,0 +1,24 @@ + +#include "cctk.h" + + MODULE cartesian_to_spherical_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE cartesian_to_spherical(theta,phi,rad,gxxs,gxys,gxzs, + & gyys,gyzs,gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys, + & dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps) + IMPLICIT NONE + CCTK_REAL :: + & rad,theta(:),phi(:) + CCTK_REAL, DIMENSION (:,:) :: + & gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys, + & dgyzs,dgzzs,grrs,grts,grps,gtts,gtps,gpps,dgtts,dgtps, + & dgpps + END SUBROUTINE + + END INTERFACE + + END MODULE cartesian_to_spherical_int diff --git a/src/energy.f b/src/energy.f new file mode 100644 index 0000000..19ef5e2 --- /dev/null +++ b/src/energy.f @@ -0,0 +1,205 @@ + PROGRAM energy + +c ----------------------------------------------------------------- +c +c Calculate energy flux from even parity variable (only real part) +c +c Now copes with having unequally spaced timesteps, it works out +c the smallest timestep, and uses this to interpolate to a new +c array with equally spaced timesteps. +c +c ----------------------------------------------------------------- + + IMPLICIT NONE + + INTEGER,PARAMETER :: nmax=5000 + INTEGER :: i,itotal,it1,it2,Nt + CHARACTER*50 filename + REAL*8 t(nmax),Qplus(nmax),dEdt(nmax),E,ddQ(nmax),Qnew(nmax) + REAL*8 Dtnew,fac,time,t1,t2,dQdt_l,dQdt_r,dtmin,dt + + WRITE(6,*) + WRITE(6,*) "Starting PROGRAM energy ..." + WRITE(6,*) "---------------------------" + + WRITE(6,131) + 131 FORMAT("Filename containing Q ? ",$) + READ(5,*) filename + WRITE(6,*) "Filename is ... ",filename + + OPEN(UNIT=12,FILE=filename,STATUS="old") + + DO i=1,nmax + READ(12,*,ERR=101) t(i),Qplus(i) + itotal = i + END DO + WRITE(6,*) "Didn't read all of data file" + STOP + + 101 CLOSE(12) + +c Give information about the data + WRITE(6,*) "We have data from t=",t(1)," to ",t(itotal) + +c Find minimum time step + dtmin = 100000.0 + do i=2,itotal + dt = t(i)-t(i-1) + if (dt t(itotal)) GOTO 301 + +c Calculate final timeslice + do i=1,itotal + if (t(i), Bowen and York, Phys. Rev. D, 21 2047(1980) +c +c ----------------------------------------------------------------------- + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi, + & hxx,hxy,hxz,hyy,hyz,hzz + + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & spin_int + +c Local variables, here only + INTEGER :: + & i,j,k,ip,count + + CCTK_REAL,PARAMETER :: + & half = 0.5D0 + CCTK_REAL :: + & rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy, + & uyz,uzz, + & Pi,idet,p, + & term1,term2 + + data count / 1 / + save count +c ------------------------------------------------------------------------ + + Pi = ACOS(-1D0) + + SELECT CASE (Psi_power) + + CASE (1) + ip = 4 + + CASE (4) + ip = 1 + + CASE DEFAULT + WRITE(*,*) "This value of Psi_power is not supported" + + END SELECT + + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + ux = (x(i)-origin(1))/rad + uy = (y(j)-origin(2))/rad + uz = (z(k)-origin(3))/rad + +c Abbreviations for metric coefficients +c ------------------------------------- + p = psi(i,j,k)**ip + + dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k) + dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k) + dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k) + +c Determinant of 3-metric +c ----------------------- + det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz + & - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2)) + + idet = 1.0/det + +c Inverse 3-metric +c ---------------- + uxx = idet*(dyy*dzz - dyz**2) + uyy = idet*(dxx*dzz - dxz**2) + uzz = idet*(dxx*dyy - dxy**2) + uxy = idet*(dxz*dyz - dzz*dxy) + uxz = idet*(dxy*dyz - dyy*dxz) + uyz = idet*(dxy*dxz - dxx*dyz) + + +c Integrands +c ---------- + + if (count.eq.1) then + term1 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k)) + & +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxy*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uyz*(ux*hxz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (uy*term1-uz*term2)*rad**3 + + else if (count.eq.2) then + term1 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uxy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uxz*(uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k)) + & +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (uz*term1-ux*term2)*rad**3 + + else if (count.eq.3) then + term1 = uxy*(ux*hxx(i,j,k)+uz*hxz(i,j,k)) + & +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uyz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + term2 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)) + & +uxy*(uy*hyy(i,j,k)+uz*hyz(i,j,k)) + & +uxz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)) + + spin_int(i,j,k) = 1.0D0/8.0D0/Pi* + & (ux*term1-uy*term2)*rad**3 + + end if + + + ELSE + + spin_int(i,j,k) = 0.0D0 + + ENDIF + + ENDDO + ENDDO + ENDDO + + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + spin_int(1,j,k) = spin_int(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + spin_int(i,1,k) = spin_int(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + spin_int(i,j,1) = spin_int(i,j,2) + ENDDO + ENDDO + + +c setting counter +c --------------- + + if (count.eq.3) then + count = 1 + else + count = count + 1 + end if + + +c end +c --- + + END SUBROUTINE spin_integrand3D diff --git a/src/spin_integrand3D_int.F b/src/spin_integrand3D_int.F new file mode 100644 index 0000000..cb2a1ca --- /dev/null +++ b/src/spin_integrand3D_int.F @@ -0,0 +1,30 @@ + +#include "cctk.h" + + MODULE spin_integrand3D_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx, + & gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,spin_int, + & Psi,Psi_power) + + IMPLICIT NONE + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi, + & hxx,hxy,hxz,hyy,hyz,hzz + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & spin_int + END SUBROUTINE spin_integrand3D + + END INTERFACE + + END MODULE spin_integrand3D_int diff --git a/src/unphysical_to_physical.F b/src/unphysical_to_physical.F new file mode 100644 index 0000000..74b4ad8 --- /dev/null +++ b/src/unphysical_to_physical.F @@ -0,0 +1,56 @@ + +#include "cctk.h" + +c ================================================================== + + SUBROUTINE unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt, + & dgtp,dgpp,Psis,dPsis,Psi_power) + +c ------------------------------------------------------------------ +c +c Convert unphysical metric components and radial (eta) derivatives +c on the 2-sphere into physical quantities, using the conformal factor +c on the sphere +c +c ------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER :: + & Psi_power + CCTK_REAL,INTENT(INOUT),DIMENSION (:,:) :: + & grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp + CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: + & Psis,dPsis + +c Output variables +c WE ARE CHANGING THE INPUT VARIABLES !!!! + +c ------------------------------------------------------------------ + + dgtt = Psis**4*dgtt + 4.0d0*Psis**3*dPsis*gtt + dgtp = Psis**4*dgtp + 4.0d0*Psis**3*dPsis*gtp + dgpp = Psis**4*dgpp + 4.0d0*Psis**3*dPsis*gpp + + grr = Psis**4*grr + grt = Psis**4*grt + grp = Psis**4*grp + gtt = Psis**4*gtt + gtp = Psis**4*gtp + gpp = Psis**4*gpp + + + END SUBROUTINE unphysical_to_physical + +c ================================================================== + + + + + + + + + + diff --git a/src/unphysical_to_physical_int.F b/src/unphysical_to_physical_int.F new file mode 100644 index 0000000..b53fa0b --- /dev/null +++ b/src/unphysical_to_physical_int.F @@ -0,0 +1,21 @@ + +#include "cctk.h" + + MODULE unphysical_to_physical_int + +c ------------------------------------------------------------------ + + INTERFACE + + SUBROUTINE unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt, + & dgtp,dgpp,Psis,dPsis,Psi_power) + IMPLICIT NONE + INTEGER :: + & Psi_power + CCTK_REAL, DIMENSION (:,:) :: + & grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp,Psis,dPsis + END SUBROUTINE + + END INTERFACE + + END MODULE unphysical_to_physical_int -- cgit v1.2.3