aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/D2_extract.F60
-rw-r--r--src/D2_extract_int.F24
-rw-r--r--src/D3_extract.F29
-rw-r--r--src/D3_extract_int.F22
-rw-r--r--src/D3_to_D2.F161
-rw-r--r--src/D3_to_D2_int.F32
-rw-r--r--src/Extract.F22
7 files changed, 274 insertions, 76 deletions
diff --git a/src/D2_extract.F b/src/D2_extract.F
index ca832f5..69f65f1 100644
--- a/src/D2_extract.F
+++ b/src/D2_extract.F
@@ -3,8 +3,8 @@
c ========================================================================
- SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11,
- & g12,g13,g22,g23,g33,dg22,dg23,dg33,
+ 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,
@@ -32,19 +32,22 @@ c Input variables
& 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,
+ INTEGER,INTENT(IN) :: do_ADMmass(2)
+ CCTK_REAL,INTENT(IN) :: eta
+ CCTK_REAL,INTENT(IN),DIMENSION (Nt) :: theta
+ CCTK_REAL,INTENT(IN),DIMENSION (Np) :: phi
+ CCTK_REAL,INTENT(IN),DIMENSION (Nt,Np) ::
+ & g00,g12,g13,g23,dg23,
& ADMmass_int1,ADMmass_int2,
& momentum_int1,momentum_int2,momentum_int3,
& spin_int1,spin_int2,spin_int3
+ CCTK_REAL,INTENT(INOUT),DIMENSION (Nt,Np) ::
+ & g11,g22,g33,dg22,dg33
+
c Output variables
CCTK_REAL,INTENT(OUT) ::
- & dtaudt,ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),
+ & ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
& momentum(3),spin(3)
c Local Variables
@@ -80,6 +83,9 @@ c 0. Initial things
c
c ------------------------------------------------------------------------
+c print*,'D2_extract Qodd'
+c print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
Pi = two*ASIN(one)
SELECT CASE (igrid)
@@ -112,6 +118,7 @@ c Calculate ADM mass
c ------------------
IF (do_ADMmass(1) == 1) THEN
+c print*,'1'
c Calculate ADM mass using standard formula
DO j = 2, Np+1
DO i = 2, Nt+1
@@ -123,6 +130,7 @@ c Calculate ADM mass using standard formula
ENDIF
IF (do_ADMmass(2) == 1) THEN
+c print*,'2'
c Calculate ADM mass using conformal formula
DO j = 2, Np+1
DO i = 2, Nt+1
@@ -137,6 +145,7 @@ c Calculate momentum
c ------------------
IF (do_momentum == 1) THEN
+c print*,'3'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = momentum_int1(i-1,j-1)*SIN(theta(i-1))
@@ -144,6 +153,7 @@ c ------------------
ENDDO
momentum(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'4'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = momentum_int2(i-1,j-1)*SIN(theta(i-1))
@@ -151,6 +161,7 @@ c ------------------
ENDDO
momentum(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'5'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = momentum_int3(i-1,j-1)*SIN(theta(i-1))
@@ -164,6 +175,7 @@ c Calculate spin
c --------------
IF (do_spin == 1) THEN
+c print*,'6'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = spin_int1(i-1,j-1)*SIN(theta(i-1))
@@ -171,6 +183,7 @@ c --------------
ENDDO
spin(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'7'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = spin_int2(i-1,j-1)*SIN(theta(i-1))
@@ -178,6 +191,7 @@ c --------------
ENDDO
spin(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'8'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = spin_int3(i-1,j-1)*SIN(theta(i-1))
@@ -194,6 +208,7 @@ c
c ------------------------------------------------------------------
c ... g00
+c print*,'9'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = g00(i-1,j-1)*SIN(theta(i-1))
@@ -202,6 +217,7 @@ c ... g00
dtaudt = sqrt(one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp))
c ... g11
+c print*,'10'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = g11(i-1,j-1)*SIN(theta(i-1))
@@ -210,6 +226,7 @@ c ... g11
sph_g11 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
c ... g22
+c print*,'11'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = g22(i-1,j-1)*SIN(theta(i-1))
@@ -218,6 +235,7 @@ c ... g22
sph_g22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
c ... dg22
+c print*,'12'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = dg22(i-1,j-1)*SIN(theta(i-1))
@@ -227,6 +245,7 @@ c ... dg22
c ... g33/sin(theta)**2
+c print*,'13'
DO j = 2, Np+1
DO i = 2, Nt+1
temp(i,j) = g33(i-1,j-1)/SIN(theta(i-1))
@@ -253,6 +272,7 @@ c 1. Compute the Schwarzschild radius
c
c ------------------------------------------------------------------
+c print*,'14'
DO j = 2, Np+1
DO i = 2, Nt+1
c temp(i,j) = g22(i-1,j-1)*SIN(theta(i-1))
@@ -273,6 +293,7 @@ c respect to the isotropic radius eta.
c
c ------------------------------------------------------------------
+c print*,'15'
DO j = 2, Np+1
DO i = 2, Nt+1
c temp(i,j) = (dg22(i-1,j-1))*SIN(theta(i-1))
@@ -305,6 +326,7 @@ c Do not know how important this is
c
c ------------------------------------------------------------------
+c print*,'16'
DO i=1,Nt
DO j=1,Np
st = SIN(theta(i))
@@ -360,12 +382,14 @@ c ------------------------------------------------------------------
c Calculate integrands
+c print*,'17'
loop_phi1: DO j = 1, Np
loop_theta1: DO i = 1, Nt
st = SIN(theta(i))
ist = one/st
+c print*,'Calling spher_harm_coms',i,j
CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2,Y3
& ,Y4)
@@ -390,20 +414,30 @@ c Calculate integrands
END DO loop_theta1
END DO loop_phi1
+c print*,'End loop'
c Integrations over the 2-sphere
h1(1) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'h1 done'
H2(1) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'H2 done'
G(1) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'G done'
K(1) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp)
& +DBLE(il*(il+1))*half*G(1)
+c print*,'K done'
c1(1) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'c1 done'
c2(1) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'c2 done'
dG(1) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'dG done'
dK(1) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp)
& +DBLE(il*(il+1))*half*dG(1)
+c print*,'dK done'
dc2(1) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp)
+c print*,'dc2 done'
c ------------------------------------------------------------------
@@ -412,10 +446,12 @@ c 5b. Imaginary parts of the Regge-Wheeler variables
c
c ------------------------------------------------------------------
+c print*,'17 1/2'
complex_parts: IF (igrid /= 1 .AND. igrid /=2) THEN
c Calculate integrands
+c print*,'18'
loop_phi2: DO j = 1, Np
loop_theta2: DO i = 1, Nt
@@ -486,11 +522,16 @@ c 6. Construct gauge invariant variables
c
c ------------------------------------------------------------------
+c print*,'Qodd',il,im
+c print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
Qodd(:,il,im) = SQRT(two*DBLE((il+2)*(il+1)*il*(il-1)))*S/
& rsch*(c1+half*(dc2-two/rsch*c2))
+c print*,'Lambda'
Lambda = DBLE((il-1)*(il+2))+3.0D0*(one-S)
+c print*,'Qeven'
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)
@@ -526,6 +567,7 @@ c ------------------------------------------------------------------
WRITE(302,*) eta,Qodd(:,lmax,mmax)
ENDIF
+c print*,'Exiting D2_extract'
END SUBROUTINE D2_extract
diff --git a/src/D2_extract_int.F b/src/D2_extract_int.F
index 43f6906..58a264d 100644
--- a/src/D2_extract_int.F
+++ b/src/D2_extract_int.F
@@ -7,8 +7,8 @@ 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,
+ 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,
@@ -19,19 +19,21 @@ c ------------------------------------------------------------------
& 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,
+ INTEGER,INTENT(IN) :: do_ADMmass(2)
+ CCTK_REAL,INTENT(IN) :: eta
+ CCTK_REAL,INTENT(IN),DIMENSION (Nt) :: theta
+ CCTK_REAL,INTENT(IN),DIMENSION (Np) :: phi
+ CCTK_REAL,INTENT(IN),DIMENSION (Nt,Np) ::
+ & g00,g12,g13,g23,dg23,
& ADMmass_int1, ADMmass_int2,
& momentum_int1,momentum_int2,momentum_int3,
& spin_int1,spin_int2,spin_int3
+
+ CCTK_REAL,INTENT(INOUT),DIMENSION (Nt,Np) ::
+ & g11,g22,g33,dg22,dg33
+
CCTK_REAL,INTENT(OUT) ::
- & ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt,
+ & ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
& momentum(3),spin(3)
END SUBROUTINE
diff --git a/src/D3_extract.F b/src/D3_extract.F
index 47a7fb3..7ef0f0f 100644
--- a/src/D3_extract.F
+++ b/src/D3_extract.F
@@ -4,7 +4,8 @@
c ==================================================================
SUBROUTINE D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,igrid,
- & origin,myproc,interpolation_operator,interpolation_order,Nt,Np,all_modes,
+ & origin,myproc,interpolation_operator,interpolation_order,Nt,Np,
+ & nx,ny,nz,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)
@@ -73,17 +74,18 @@ c Input variables
INTEGER,INTENT(IN) ::
& conformal_state,igrid,l,m,Psi_power,myproc
CCTK_INT,INTENT(IN) ::
- & Nt,Np,all_modes,do_momentum,do_spin,interpolation_order
+ & Nt,Np,nx,ny,nz,all_modes,do_momentum,do_spin,interpolation_order
INTEGER,INTENT(IN) ::
& do_ADMmass(2)
CCTK_REAL,INTENT(IN) ::
& origin(3),Dx,Dy,Dz,eta
- CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) ::
+ CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) ::
& 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(:,:,:) ::
+ CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+ CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+ CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+ CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
& Extract_temp3d
CCTK_STRING,INTENT(IN) ::
& interpolation_operator
@@ -122,6 +124,11 @@ c 1. Specify polar coordinates on chosen 2-surface
c
c ------------------------------------------------------------------
+c print*,'Entering D3_extract'
+c print*,nx,ny,nz
+c print*,'D3_extract Qodd'
+c print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
Pi = two*ASIN(one)
c Set polar gridspacing
@@ -163,14 +170,17 @@ c 2. Project quantities onto the 2-surface
c
c ------------------------------------------------------------------
+c print*,nx,ny,nz,Nt,Np
+c print*,'Calling D3_to_D2'
CALL D3_to_D2(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
& Psi_power,origin,myproc,interpolation_operator,interpolation_order,
& 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,
+ & x,y,z,eta,Nt,Np,nx,ny,nz,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 print*,'Exited D3_to_D2'
if (myproc .eq. 0) then
@@ -203,12 +213,17 @@ c 5. Extract gauge invariant quantities on 2-surface
c
c ------------------------------------------------------------------
+c print*,'Calling D2_extract'
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 print*,Qodd
+c print*,Qeven
+c print*,'Exiting D2_extract'
+
c End of myproc = 0 conditional
endif
diff --git a/src/D3_extract_int.F b/src/D3_extract_int.F
index ce14b37..e00b321 100644
--- a/src/D3_extract_int.F
+++ b/src/D3_extract_int.F
@@ -9,7 +9,7 @@ c ------------------------------------------------------------------
SUBROUTINE D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,igrid,
& origin,myproc,interpolation_operator,interpolation_order,
- & Nt,Np,all_modes,l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi,
+ & Nt,Np,nx,ny,nz,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)
@@ -19,26 +19,28 @@ c ------------------------------------------------------------------
CCTK_POINTER :: cctkGH
INTEGER,INTENT(IN) ::
- & igrid,l,m,Psi_power,myproc
+ & conformal_state,igrid,l,m,Psi_power,myproc
CCTK_INT,INTENT(IN) ::
- & conformal_state,Nt,Np,all_modes,do_momentum,do_spin,interpolation_order
+ & Nt,Np,nx,ny,nz,all_modes,do_momentum,do_spin,interpolation_order
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(:,:,:) ::
+ CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) ::
& Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
& hxx,hxy,hxz,hyy,hyz,hzz
- CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
+ CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+ CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+ CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+ CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
& Extract_Temp3d
- CCTK_REAL,INTENT(OUT) ::
- & ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt,
- & momentum(3),spin(3)
CCTK_STRING,INTENT(IN) ::
& interpolation_operator
+ CCTK_REAL,INTENT(OUT) ::
+ & ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
+ & momentum(3),spin(3)
+
END SUBROUTINE
diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F
index da2294e..1e71ff4 100644
--- a/src/D3_to_D2.F
+++ b/src/D3_to_D2.F
@@ -7,7 +7,7 @@ c ========================================================================
& Psi_power,origin,myproc,interpolation_operator,interpolation_order,
& 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,
+ & x,y,z,eta,Nt,Np,nx,ny,nz,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,
@@ -35,24 +35,26 @@ c Input variables
INTEGER,INTENT(IN) ::
& conformal_state,myproc,Psi_power
CCTK_INT, INTENT(IN) ::
- & Nt,Np,do_momentum,do_spin,interpolation_order
+ & Nt,Np,nx,ny,nz,do_momentum,do_spin,interpolation_order
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(:,:,:) ::
+ CCTK_REAL,INTENT(IN),DIMENSION(Nt) :: theta
+ CCTK_REAL,INTENT(IN),DIMENSION(Np) :: phi
+ CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+ CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+ CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+ CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) ::
& 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)
+ INTEGER,INTENT(IN) :: do_ADMmass(2)
CCTK_STRING,INTENT(IN) ::
& interpolation_operator
+ CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
+ & Extract_temp3d
c Output variables
- CCTK_REAL,INTENT(OUT),DIMENSION(:,:) ::
+ CCTK_REAL,INTENT(OUT),DIMENSION(Nt,Np) ::
& Psis,g00s,gxxs,gxys,gxzs,gyys,
& gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs,
& ADMmass_int1,ADMmass_int2,
@@ -65,7 +67,7 @@ c Local variables, passed on
LOGICAL ::
& err_flag
INTEGER ::
- & iorder,npoints,Nx,Ny,Nz
+ & iorder,npoints
INTEGER,DIMENSION(Nt,Np) ::
& ib,jb,kb
CCTK_REAL,DIMENSION(Nt,Np) ::
@@ -83,18 +85,14 @@ c Local variables, here only
character(128) options_string
character(128) operator
CCTK_INT nchars
+ character(len=20) :: error_code
c ------------------------------------------------------------------------
-c Initial Stuff
-c -------------
- Nx = SIZE(x)
- Ny = SIZE(y)
- Nz = SIZE(z)
-
-
c Compute Cartesian coordinates on the surface
c --------------------------------------------
+c print*,'Entering D3_to_D2'
+c print*,nx,ny,nz,Np,Nt
DO j = 1, Np
DO i = 1, Nt
@@ -105,6 +103,7 @@ c --------------------------------------------
ENDDO
ENDDO
+c print*,'Coordinates set'
c Only do interpolation on one processor
c --------------------------------------
@@ -118,6 +117,7 @@ c --------------------------------------
END SELECT
+c print*,'Number of points set'
c Get the interpolator, parameter table, and coordinate system handles
c --------------------------------------------------------------------
@@ -143,6 +143,7 @@ c --------------------------------------------------------------------
call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
endif
+c print*,'Handles obtained'
c fill in the input/output arrays for the interpolator
c ----------------------------------------------------
@@ -179,13 +180,19 @@ c ------------------------------------------------------------
num_arrays = 7
end if
+c print*,'1'
+c print*,in_array_indices
+c print*,out_arrays
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ num_arrays, in_array_indices,
$ num_arrays, out_array_type_codes, out_arrays)
+c print*,'1 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
@@ -194,57 +201,87 @@ c ----------------------------------------------------
call CCTK_VarIndex(in_array_indices(1), "extract::temp3d")
if (conformal_state > 0) then
+c print*,'1 1/2'
CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d)
+c print*,'1 1/2 done'
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
out_arrays(1) = CCTK_PointerTo(dPsis)
+c print*,'2'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'2 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
end if
+c print*,'3'
CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d)
+c print*,'3 done'
+c print*,'4'
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'4 done'
out_arrays(1) = CCTK_PointerTo(dgxxs)
+c print*,'5'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'5 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,Extract_temp3d)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
out_arrays(1) = CCTK_PointerTo(dgxys)
+c print*,'6'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'6 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,Extract_temp3d)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
out_arrays(1) = CCTK_PointerTo(dgxzs)
+c print*,'7'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'7 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,Extract_temp3d)
@@ -258,56 +295,80 @@ c ----------------------------------------------------
$ 1, out_array_type_codes, out_arrays(1))
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,Extract_temp3d)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
out_arrays(1) = CCTK_PointerTo(dgyzs)
+c print*,'8'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'8 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,Extract_temp3d)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
out_arrays(1) = CCTK_PointerTo(dgzzs)
+c print*,'9'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'9 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
c Calculate integrands for ADM masses
c -----------------------------------
c Standard equation
- IF (do_ADMmass(1)) THEN
+ IF (do_ADMmass(1) == 1) THEN
+c print*,'10'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
CALL ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,
& gxz,gyy,gyz,gzz,Extract_temp3d,Psi,Psi_power,conformal_state)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'10 done'
out_arrays(1) = CCTK_PointerTo(ADMmass_int1)
+c print*,'11'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'11 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
END IF
c Conformal equation
- IF (do_ADMmass(2)) THEN
+ IF (do_ADMmass(2) == 1) THEN
ADMmass_int2 = -eta**2/2D0/3.1416D0*dPsis
ENDIF
@@ -315,49 +376,79 @@ c Calculate integrands for momentum
c ---------------------------------
IF (do_momentum==1) THEN
+c print*,'12'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
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,conformal_state)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'12 done'
out_arrays(1) = CCTK_PointerTo(momentum_int1)
+c print*,'13'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'13 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
+c print*,'14'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
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,conformal_state)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'14 done'
out_arrays(1) = CCTK_PointerTo(momentum_int2)
+c print*,'15'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'15 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
+c print*,'16'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
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,conformal_state)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'16 done'
out_arrays(1) = CCTK_PointerTo(momentum_int3)
+c print*,'17'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'17 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
END IF
@@ -365,49 +456,79 @@ c Calculate integrands for spin
c -----------------------------
IF (do_spin==1) THEN
+c print*,'18'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
CALL spin_integrand3D(origin,x,y,z,
& gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
& Extract_temp3d,Psi,Psi_power,conformal_state)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'18 done'
out_arrays(1) = CCTK_PointerTo(spin_int1)
+c print*,'19'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'19 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
+c print*,'20'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
CALL spin_integrand3D(origin,x,y,z,
& gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
& Extract_temp3d,Psi,Psi_power,conformal_state)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'20 done'
out_arrays(1) = CCTK_PointerTo(spin_int2)
+c print*,'21'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'21 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
+c print*,'22'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
CALL spin_integrand3D(origin,x,y,z,
& gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
& Extract_temp3d,Psi,Psi_power,conformal_state)
CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c print*,'22 done'
out_arrays(1) = CCTK_PointerTo(spin_int3)
+c print*,'23'
+c print*,in_array_indices(1)
+c print*,out_arrays(1)
call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
$ param_table_handle, coord_system_handle,
$ npoints, CCTK_VARIABLE_REAL, interp_coords,
$ 1, in_array_indices(1),
$ 1, out_array_type_codes, out_arrays(1))
+c print*,'23 done'
if (ierror < 0) then
call CCTK_WARN (1, "interpolator call returned an error code");
+ write(error_code,'(a7,i13)') 'code = ',ierror
+ call CCTK_WARN (1, error_code );
endif
END IF
diff --git a/src/D3_to_D2_int.F b/src/D3_to_D2_int.F
index db013a2..4aeef0a 100644
--- a/src/D3_to_D2_int.F
+++ b/src/D3_to_D2_int.F
@@ -11,7 +11,7 @@ c ------------------------------------------------------------------
& Psi_power,origin,myproc,interpolation_operator,interpolation_order,
& 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,
+ & x,y,z,eta,Nt,Np,nx,ny,nz,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,
@@ -22,26 +22,30 @@ c ------------------------------------------------------------------
INTEGER,INTENT(IN) ::
& conformal_state,myproc,Psi_power
CCTK_INT, INTENT(IN) ::
- & Nt,Np,do_momentum,do_spin,interpolation_order
- INTEGER ::
- & do_ADMmass(2)
+ & Nt,Np,nx,ny,nz,do_momentum,do_spin,interpolation_order
CCTK_REAL,INTENT(IN) ::
& origin(3),Dx,Dy,Dz,eta
- CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) ::
+ CCTK_REAL,INTENT(IN),DIMENSION(Nt) :: theta
+ CCTK_REAL,INTENT(IN),DIMENSION(Np) :: phi
+ CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+ CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+ CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+ CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) ::
& 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
+ INTEGER,INTENT(IN) :: do_ADMmass(2)
CCTK_STRING,INTENT(IN) ::
& interpolation_operator
+ CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
+ & Extract_temp3d
+
+ CCTK_REAL,INTENT(OUT),DIMENSION(Nt,Np) ::
+ & 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
END SUBROUTINE
END INTERFACE
diff --git a/src/Extract.F b/src/Extract.F
index ace99fb..b4ba753 100644
--- a/src/Extract.F
+++ b/src/Extract.F
@@ -83,6 +83,7 @@ c Local variables
INTEGER myproc,ierr
CCTK_REAL dx,dy,dz
+ CCTK_INT :: nx, ny, nz
character*80 infoline
@@ -93,6 +94,10 @@ c ------------------------------------------------------------------
dy = cctk_delta_space(2)
dz = cctk_delta_space(3)
+ nx = cctk_lsh(1)
+ ny = cctk_lsh(2)
+ nz = cctk_lsh(3)
+
c ------------------------------------------------------------------
c
c 0. Check to see if Extract should have been called
@@ -149,13 +154,13 @@ c Set the value of igrid
c Create 1D coordinate arrays
- do ix = 1, cctk_lsh(1)
+ do ix = 1, nx
x_1d(ix) = x(ix,1,1)
enddo
- do iy = 1, cctk_lsh(2)
+ do iy = 1, ny
y_1d(iy) = y(1,iy,1)
enddo
- do iz = 1, cctk_lsh(3)
+ do iz = 1, nz
z_1d(iz) = z(1,1,iz)
enddo
@@ -218,6 +223,9 @@ c Get m_mode in this case too, or the T3E has a junk value (PW)
mmode = m_mode
ALLOCATE(Qodd(2,2:lmode,0:mmode),Qeven(2,2:lmode,0:mmode))
ENDIF
+c print*,'Extract Qodd'
+c print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
c Check do not have bad values for the modes
@@ -377,13 +385,17 @@ c Do extraction at each radius
extract_at_each_radius: DO WHILE (radius < r2)
+c print*,'Calling D3_extract'
+c print*,nx,ny,nz,Nt,Np
CALL D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
& igrid,orig,myproc,interpolation_operator,interpolation_order,
- & Nt,Np,all_modes,lmode,
+ & Nt,Np,nx,ny,nz,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)
+c print*,Qodd,Qeven
+c print*,'Exited D3_extract'
IF (verbose == 1) THEN
WRITE(*,*) "Extracted at r =",radius
@@ -551,7 +563,7 @@ c Cannot use the conformal equation for ADM mass now
CALL D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
& igrid,orig,myproc,interpolation_operator,interpolation_order,
- & Nt,Np,all_modes,
+ & Nt,Np,nx,ny,nz,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,
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 14:45:55 +0000