diff options
Diffstat (limited to 'src/D3_extract.F')
| -rw-r--r-- | src/D3_extract.F | 221 |
1 files changed, 221 insertions, 0 deletions
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 + + + |