diff options
author | diener <diener@2a26948c-0e4f-0410-aee8-f1d3e353619c> | 2003-06-20 08:21:27 +0000 |
---|---|---|
committer | diener <diener@2a26948c-0e4f-0410-aee8-f1d3e353619c> | 2003-06-20 08:21:27 +0000 |
commit | 88a32348f469a016ee14a18ea4f52f3817943a6c (patch) | |
tree | 2e99eee10ffdd85584d0b9b9f47cafaacdfaa9bf /src | |
parent | 57e38b49b43b1501449900c68c1790ab0e57dc65 (diff) |
Added support for setting up generators on ellipsoidal initial data.
Added support for evolving the generators independently (used for testing).
Added support for detecting a peanut shaped excision region and take the
appropriate precautions.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@119 2a26948c-0e4f-0410-aee8-f1d3e353619c
Diffstat (limited to 'src')
-rw-r--r-- | src/EHFinder_Init.F90 | 94 | ||||
-rw-r--r-- | src/EHFinder_SetMask.F90 | 145 |
2 files changed, 203 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... diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90 index 21d7fc9..9a64ca0 100644 --- a/src/EHFinder_SetMask.F90 +++ b/src/EHFinder_SetMask.F90 @@ -570,3 +570,148 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) end if end subroutine EHFinder_SetMask2 + + +subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k + logical active, mask_modified + + active = .false. + +! If re-parametrization has just been done, we are in the process of +! resetting the mask so we set active=.true. + if ( mask_first ) active = .true. + if ( CCTK_EQUALS(re_param_method,'approx') .or. & + CCTK_EQUALS(re_param_method,'mixed') ) then + if ( reparametrize_every_approx .gt. 0 ) then + if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) then + active = .true. + end if + end if + end if + if ( CCTK_EQUALS(re_param_method,'pde') .or. & + CCTK_EQUALS(re_param_method,'mixed') ) then + if ( reparametrize_every_pde .gt. 0 ) then + if ( mod(cctk_iteration,reparametrize_every_pde) .eq. 0 ) active = .true. + end if + end if + +! If the reparametrization was not undone... + if ( active .and. .not.reparam_undone ) then + + mask_modified = .false. +# include "include/physical_part.h" + +! Make sure we are not on the physical outer boundary. + if ( ixl .eq. 1 ) ixl = 2 + if ( ixr .eq. nx ) ixr = nx - 1 + if ( jyl .eq. 1 ) jyl = 2 + if ( jyr .eq. ny ) jyr = ny - 1 + if ( kzl .eq. 1 ) kzl = 2 + if ( kzr .eq. nz ) kzr = nz - 1 + + tm_mask = eh_mask + + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k) .ge. 0 ) then + ! If we are on an excision boundary in the -x-direction... + if ( btest(eh_mask(i,j,k),0) ) then + ! If any of the two closest neighbours in the +x-direction is + ! excised we have to excise this point + if ( ( eh_mask(i+1,j,k) .eq. -1 ) .or. & + ( eh_mask(i+2,j,k) .eq. -1 ) ) then + tm_mask(i,j,k) = -1 + mask_modified = .true. + print*, i, j, k, 'x-' + print*, eh_mask(i:i+2,j,k) + end if + end if + + ! If we are on an excision boundary in the +x-direction... + if ( btest(eh_mask(i,j,k),1) ) then + ! If any of the two closest neighbours in the -x-direction is + ! excised we have to excise this point + if ( ( eh_mask(i-1,j,k) .eq. -1 ) .or. & + ( eh_mask(i-2,j,k) .eq. -1 ) ) then + tm_mask(i,j,k) = -1 + mask_modified = .true. + print*, i, j, k, 'x+' + print*, eh_mask(i:i-2,j,k) + end if + end if + + ! If we are on an excision boundary in the -y-direction... + if ( btest(eh_mask(i,j,k),2) ) then + ! If any of the two closest neighbours in the +y-direction is + ! excised we have to excise this point + if ( ( eh_mask(i,j+1,k) .eq. -1 ) .or. & + ( eh_mask(i,j+2,k) .eq. -1 ) ) then + tm_mask(i,j,k) = -1 + mask_modified = .true. + print*, i, j, k, 'y-' + print*, eh_mask(i,j:j+2,k) + end if + end if + + ! If we are on an excision boundary in the +y-direction... + if ( btest(eh_mask(i,j,k),3) ) then + ! If any of the two closest neighbours in the -y-direction is + ! excised we have to excise this point + if ( ( eh_mask(i,j-1,k) .eq. -1 ) .or. & + ( eh_mask(i,j-2,k) .eq. -1 ) ) then + tm_mask(i,j,k) = -1 + mask_modified = .true. + print*, i, j, k, 'y+' + print*, eh_mask(i,j:j-2,k) + end if + end if + + ! If we are on an excision boundary in the -z-direction... + if ( btest(eh_mask(i,j,k),4) ) then + ! If any of the two closest neighbours in the +z-direction is + ! excised we have to excise this point + if ( ( eh_mask(i,j,k+1) .eq. -1 ) .or. & + ( eh_mask(i,j,k+2) .eq. -1 ) ) then + tm_mask(i,j,k) = -1 + mask_modified = .true. + print*, i, j, k, 'z-' + print*, eh_mask(i,j,k:k+2) + end if + end if + + ! If we are on an excision boundary in the +z-direction... + if ( btest(eh_mask(i,j,k),5) ) then + ! If any of the two closest neighbours in the -z-direction is + ! excised we have to excise this point + if ( ( eh_mask(i,j,k-1) .eq. -1 ) .or. & + ( eh_mask(i,j,k-2) .eq. -1 ) ) then + tm_mask(i,j,k) = -1 + mask_modified = .true. + print*, i, j, k, 'z+' + print*, eh_mask(i,j,k:k-2) + end if + end if + end if + + end do + end do + end do + + if ( mask_modified ) then + call CCTK_WARN ( 1, "Mask modified because it was not convex" ) + end if + eh_mask = tm_mask + end if + +end subroutine EHFinder_SetMask3 |