diff options
Diffstat (limited to 'src/EHFinder_Init.F90')
| -rw-r--r-- | src/EHFinder_Init.F90 | 94 |
1 files changed, 58 insertions, 36 deletions
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90 index 3036a1d..2bf3d7d 100644 --- a/src/EHFinder_Init.F90 +++ b/src/EHFinder_Init.F90 @@ -20,7 +20,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) CCTK_REAL, dimension(3,3) :: txyz CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc CCTK_REAL :: last_time - CCTK_REAL :: theta, dtheta, thetamin, thetamax + CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el CCTK_INT, dimension(1) :: lsh, lbnd ! Initialize the current_iteration variable. @@ -35,6 +35,46 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) cctk_time = last_time + if ( evolve_generators .gt. 0 ) then + + call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::generators" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "cannot get lower bounds for generator arrays" ) + end if + call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "cannot get local size for generator arrays" ) + end if + + if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then + + if ( CCTK_EQUALS( domain, 'full' ) ) then + thetamin = zero; thetamax = pi + else if ( CCTK_EQUALS( domain, 'bitant') ) then + if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then + thetamin = zero; thetamax = half * pi + else + thetamin = zero; thetamax = pi + end if + else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then + if ( CCTK_EQUALS( quadrant_direction, 'x' ) .or. & + CCTK_EQUALS( quadrant_direction, 'y' ) ) then + thetamin = zero; thetamax = half * pi + else + thetamin = zero; thetamax = pi + end if + else if ( CCTK_EQUALS( domain, 'octant' ) ) then + thetamin = zero; thetamax = half * pi + end if + + if ( number_of_generators .eq. 1 ) then + theta = half * ( thetamax - thetamin ) + thetamin + else + dtheta = ( thetamax - thetamin ) / ( number_of_generators - 1 ) + end if + end if + end if + ! If a sphere is requested... if ( CCTK_EQUALS( initial_f, 'sphere' ) ) then @@ -44,42 +84,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) ( y - translate_y )**2 + & ( z - translate_z )**2 ) - initial_rad if ( evolve_generators .gt. 0 ) then - - call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::generators" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get lower bounds for generator arrays" ) - end if - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get local size for generator arrays" ) - end if - if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then - - if ( CCTK_EQUALS( domain, 'full' ) ) then - thetamin = zero; thetamax = pi - else if ( CCTK_EQUALS( domain, 'bitant') ) then - if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then - thetamin = zero; thetamax = half * pi - else - thetamin = zero; thetamax = pi - end if - else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then - if ( CCTK_EQUALS( quadrant_direction, 'x' ) .or. & - CCTK_EQUALS( quadrant_direction, 'y' ) ) then - thetamin = zero; thetamax = half * pi - else - thetamin = zero; thetamax = pi - end if - else if ( CCTK_EQUALS( domain, 'octant' ) ) then - thetamin = zero; thetamax = half * pi - end if - - if ( number_of_generators .eq. 1 ) then - theta = half * ( thetamax - thetamin ) + thetamin - else - dtheta = ( thetamax - thetamin ) / ( number_of_generators - 1 ) - end if do i = 1, lsh(1) theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) xg(i) = initial_rad * sin(theta) + translate_x @@ -130,6 +135,23 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) end do end do end do + + 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 ) + r_el = sqrt ( one / ( sin(theta)**2 / initial_a**2 + & + cos(theta)**2 / initial_c**2 ) ) + xp(1) = r_el * sin(theta) + translate_x + xp(2) = translate_y + xp(3) = r_el * cos(theta) + translate_z + xpt = matmul ( txyz, xp ) + xg(i) = xpt(1) + yg(i) = xpt(2) + zg(i) = xpt(3) + end do + end if + end if end if ! if an ovaloid of Cassini is requested... |