aboutsummaryrefslogtreecommitdiff
path: root/src/D3_to_D2.F
diff options
context:
space:
mode:
Diffstat (limited to 'src/D3_to_D2.F')
-rw-r--r--src/D3_to_D2.F161
1 files changed, 141 insertions, 20 deletions
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
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 14:08:33 +0000