From cc166f4d2e715bdfdef04e217602a18acdbef1ae Mon Sep 17 00:00:00 2001 From: diener Date: Fri, 13 Feb 2004 18:13:42 +0000 Subject: Now the 2d grid array for the generators should work. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@166 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_Generator_Sources.F90 | 73 ++------------------------------------ src/EHFinder_Init.F90 | 53 ++++++++++++++++++++------- 2 files changed, 44 insertions(+), 82 deletions(-) diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90 index 3dee327..ab19d29 100644 --- a/src/EHFinder_Generator_Sources.F90 +++ b/src/EHFinder_Generator_Sources.F90 @@ -194,18 +194,7 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) dxg(i,l) = - betaxg(i) + ssign * factor * dfux dyg(i,l) = - betayg(i) + ssign * factor * dfuy dzg(i,l) = - betazg(i) + ssign * factor * dfuz - if ( i == 5 ) then - print*,'Debug : ', xg(i,l), yg(i,l), zg(i,l), & - dxg(i,l), dyg(i,l), dzg(i,l) - print*,'lapse : ', alp2 - print*,'shift : ', betaxg(i), betayg(i), betazg(i) - print*,'metric : ', gxxg(i), gxyg(i), gxzg(i), & - gyyg(i), gyzg(i), gzzg(i) - print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz - print*,'df : ', dfxg(i), dfyg(i), dfzg(i) - print*,'dfu : ', dfux, dfuy, dfuz - print*,'factor : ', factor - end if + end do else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then call Util_TableSetIntArray ( status, table_handle, 14, & @@ -306,8 +295,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) CCTK_INT, dimension(14) :: out_types CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz - CCTK_REAL, dimension(:,:), allocatable :: tmp, postmp - CCTK_INT, dimension(1) :: shape ! Set the sign correctly depending on the surface direction. if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one @@ -355,12 +342,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) end if - shape(1) = lsh(1)*lsh(2) - allocate ( tmp(lsh(1)*lsh(2),14), postmp(lsh(1)*lsh(2),3) ) - - do l = 1, 14 - out_arrays(l) = CCTK_PointerTo(tmp(:,l)) - end do out_arrays(1) = CCTK_PointerTo(alpg2) out_arrays(2) = CCTK_PointerTo(betaxg2) out_arrays(3) = CCTK_PointerTo(betayg2) @@ -376,18 +357,12 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) out_arrays(13) = CCTK_PointerTo(dfzg2) out_arrays(14) = CCTK_PointerTo(psig2) -!! Set the pointers to the output arrays. +! Set the pointers to the output arrays. ! Loop over the level sets do l = 1, eh_number_level_sets ! Set the pointers to the points to be interpolated to. -! postmp(:,1) = reshape(xg2(:,:,l),shape) -! postmp(:,2) = reshape(yg2(:,:,l),shape) -! postmp(:,3) = reshape(zg2(:,:,l),shape) -! interp_coords(1) = CCTK_PointerTo(postmp(:,1)) -! interp_coords(2) = CCTK_PointerTo(postmp(:,2)) -! interp_coords(3) = CCTK_PointerTo(postmp(:,3)) interp_coords(1) = CCTK_PointerTo(xg2(1,1,l)) interp_coords(2) = CCTK_PointerTo(yg2(1,1,l)) @@ -442,20 +417,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) call CCTK_INFO ( 'Interpolation failed.' ) end if -! alpg2 = reshape(tmp(:,1),lsh) -! betaxg2 = reshape(tmp(:,2),lsh) -! betayg2 = reshape(tmp(:,3),lsh) -! betazg2 = reshape(tmp(:,4),lsh) -! gxxg2 = reshape(tmp(:,5),lsh) -! gxyg2 = reshape(tmp(:,6),lsh) -! gxzg2 = reshape(tmp(:,7),lsh) -! gyyg2 = reshape(tmp(:,8),lsh) -! gyzg2 = reshape(tmp(:,9),lsh) -! gzzg2 = reshape(tmp(:,10),lsh) -! dfxg2 = reshape(tmp(:,11),lsh) -! dfyg2 = reshape(tmp(:,12),lsh) -! dfzg2 = reshape(tmp(:,13),lsh) - ! For each point on this processor calculate the right hand side of the ! characteristic evolution equation. do j = 1, lsh(2) @@ -494,19 +455,7 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) dxg2(i,j,l) = - betaxg2(i,j) + ssign * factor * dfux dyg2(i,j,l) = - betayg2(i,j) + ssign * factor * dfuy dzg2(i,j,l) = - betazg2(i,j) + ssign * factor * dfuz - if ( j == 1 .and. i == 5 ) then - - print*,'debug2 : ', xg2(i,j,l), yg2(i,j,l), zg2(i,j,l), & - dxg2(i,j,l), dyg2(i,j,l), dzg2(i,j,l) - print*,'lapse : ', alp2 - print*,'shift : ', betaxg2(i,j), betayg2(i,j), betazg2(i,j) - print*,'metric : ', gxxg2(i,j), gxyg2(i,j), gxzg2(i,j), & - gyyg2(i,j), gyzg2(i,j), gzzg2(i,j) - print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz - print*,'df : ', dfxg2(i,j), dfyg2(i,j), dfzg2(i,j) - print*,'dfu : ', dfux, dfuy, dfuz - print*,'factor : ', factor - end if + end do end do else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then @@ -534,21 +483,6 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) call CCTK_INFO ( 'Interpolation failed.' ) end if -! alpg2 = reshape(tmp(:,1),lsh) -! betaxg2 = reshape(tmp(:,2),lsh) -! betayg2 = reshape(tmp(:,3),lsh) -! betazg2 = reshape(tmp(:,4),lsh) -! gxxg2 = reshape(tmp(:,5),lsh) -! gxyg2 = reshape(tmp(:,6),lsh) -! gxzg2 = reshape(tmp(:,7),lsh) -! gyyg2 = reshape(tmp(:,8),lsh) -! gyzg2 = reshape(tmp(:,9),lsh) -! gzzg2 = reshape(tmp(:,10),lsh) -! dfxg2 = reshape(tmp(:,11),lsh) -! dfyg2 = reshape(tmp(:,12),lsh) -! dfzg2 = reshape(tmp(:,13),lsh) -! psig2 = reshape(tmp(:,14),lsh) - do j = 1, lsh(2) do i = 1, lsh(1) alp2 = alpg2(i,j)**2 @@ -590,6 +524,5 @@ subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) end if end do - deallocate ( tmp, postmp ) return end subroutine EHFinder_Generator_Sources_2D diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90 index c9dd60c..d9aa13e 100644 --- a/src/EHFinder_Init.F90 +++ b/src/EHFinder_Init.F90 @@ -23,6 +23,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) CCTK_REAL :: last_time CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el CCTK_REAL :: phi, dphi, phimin, phimax + CCTK_REAL :: costh, sinth, cosph, sinph CCTK_INT, dimension(1) :: lsh, lbnd CCTK_INT, dimension(2) :: lsh2, lbnd2 @@ -115,13 +116,13 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) if ( number_of_generators_theta .eq. 1 ) then theta = half * ( thetamax - thetamin ) + thetamin else - dtheta = ( thetamax - thetamin ) / ( number_of_generators_theta - 1 ) + dtheta = ( thetamax - thetamin ) / ( number_of_generators_theta ) end if if ( number_of_generators_phi .eq. 1 ) then phi = half * ( phimax - phimin ) + phimin else - dphi = ( phimax - phimin ) / ( number_of_generators_phi - 1 ) + dphi = ( phimax - phimin ) / ( number_of_generators_phi ) end if end if @@ -140,20 +141,21 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) if ( evolve_generators .gt. 0 ) then if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then do i = 1, lsh(1) - theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) + theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) + half * dtheta xg(i,l) = initial_rad(l) * sin(theta) + translate_x(l) yg(i,l) = translate_y(l) zg(i,l) = initial_rad(l) * cos(theta) + translate_z(l) end do else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then do j = 1, lsh2(2) - phi = phimin + dphi * ( j + lbnd2(2) - 1 ) + phi = phimin + dphi * ( j + lbnd2(2) - 1 ) + half * dphi + cosph = cos(phi); sinph = sin(phi) do i = 1, lsh2(1) theta = thetamin + dtheta * ( i + lbnd2(1) - 1 ) - xg2(i,j,l) = initial_rad(l) * sin(theta) * cos(phi) + translate_x(l) - yg2(i,j,l) = initial_rad(l) * sin(theta) * sin(phi) + translate_x(l) - zg2(i,j,l) = initial_rad(l) * cos(theta) + translate_z(l) - print*,xg2(i,j,l),yg2(i,j,l),zg2(i,j,l) + costh = cos(theta); sinth = sin(theta) + xg2(i,j,l) = initial_rad(l) * sinth * cosph + translate_x(l) + yg2(i,j,l) = initial_rad(l) * sinth * sinph + translate_x(l) + zg2(i,j,l) = initial_rad(l) * costh + translate_z(l) end do end do end if @@ -204,16 +206,43 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then do i = 1, lsh(1) theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) - r_el = sqrt ( one / ( sin(theta)**2 / initial_a(l)**2 + & - cos(theta)**2 / initial_c(l)**2 ) ) - xp(1) = r_el * sin(theta) + translate_x(l) + costh = cos(theta); sinth = sin(theta) + r_el = sqrt ( one / ( sinth**2 / initial_a(l)**2 + & + costh**2 / initial_c(l)**2 ) ) + xp(1) = r_el * sinth + translate_x(l) xp(2) = translate_y(l) - xp(3) = r_el * cos(theta) + translate_z(l) + xp(3) = r_el * costh + translate_z(l) xpt = matmul ( txyz, xp ) xg(i,l) = xpt(1) yg(i,l) = xpt(2) zg(i,l) = xpt(3) end do + else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then + do j = 1, lsh2(2) + phi = phimin + dphi * ( j + lbnd2(2) - 1 ) + half * dphi + cosph = cos(phi); sinph = sin(phi) + do i = 1, lsh2(1) + theta = thetamin + dtheta * ( i + lbnd2(1) - 1 ) + half * dtheta + costh = cos(theta); sinth = sin(theta) + r_el = sqrt ( one / & + ( sinth**2*cosph**2 / initial_a(l)**2 + & + sinth**2*sinph**2 / initial_b(l)**2 + & + costh**2 / initial_c(l)**2 ) ) + xp(1) = r_el * sinth * cosph + translate_x(l) + xp(2) = r_el * sinth * sinph + translate_y(l) + xp(3) = r_el * costh + translate_z(l) + xpt = matmul ( txyz, xp ) + xg2(i,j,l) = xpt(1) + yg2(i,j,l) = xpt(2) + zg2(i,j,l) = xpt(3) + end do + end do +! print*,xg2 +! print* +! print*,yg2 +! print* +! print*,zg2 +! stop end if end if end if -- cgit v1.2.3