From 91853ef68b64448bc7d902de45084407b7e18b77 Mon Sep 17 00:00:00 2001 From: diener Date: Wed, 25 Jun 2003 10:49:33 +0000 Subject: Changed assumed shape arrays into explicitly declared arrays to make the code run on seaborg (it was segmentation faulting before). This included passing the nx, ny, and nz of gridfunctions into the subroutines. Also made sure that the interfaces for the subroutine calls agree with the actual subroutines. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Extract/trunk@79 5301f0c2-dbc4-4cee-b2f5-8d7afba4d129 --- src/D2_extract.F | 60 ++++++++++++++++--- src/D2_extract_int.F | 24 ++++---- src/D3_extract.F | 29 +++++++--- src/D3_extract_int.F | 22 +++---- src/D3_to_D2.F | 161 ++++++++++++++++++++++++++++++++++++++++++++------- src/D3_to_D2_int.F | 32 +++++----- src/Extract.F | 22 +++++-- 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, -- cgit v1.2.3