From 76d0fff25c2a024fe8a588543afcdf8c220c4ac7 Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 22 Jun 2004 13:26:25 +0000 Subject: Moved stuff around so that global mode and level mode operations are clearly separated, so that the routines can be scheduled with Carpet. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@178 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_ReInitialize.F90 | 360 +++++++++++++++++++++++++++--------------- 1 file changed, 231 insertions(+), 129 deletions(-) diff --git a/src/EHFinder_ReInitialize.F90 b/src/EHFinder_ReInitialize.F90 index 610d94d..5daf5fd 100644 --- a/src/EHFinder_ReInitialize.F90 +++ b/src/EHFinder_ReInitialize.F90 @@ -18,87 +18,141 @@ subroutine EHFinder_ReInitializeControl(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_REAL, dimension(eh_number_level_sets) :: fmin, fmax, fminloc, fmaxloc - CCTK_INT :: i, j, k, l - character(len=128) :: info_message + CCTK_REAL, dimension(eh_number_level_sets) :: fmin, fmax, fminl, fmaxl + CCTK_INT :: i, j, l !, re_init_control + character(len=128) :: info_message, variable_name + character(11) :: form(4) = (/ '(a12,i1,a1)', '(a12,i2,a1)', & + '(a12,i3,a1)', '(a12,i4,a1)' /) ! Initialize the control variables. The default means no re-initialization. re_init_control = 0 ! Check if it is time for a pde re-initialization and set the control ! variable to 1. - if ( re_initialize_every .gt. 0 ) then - if ( mod(cctk_iteration,re_initialize_every) .eq. 0 ) then + if ( re_initialize_every > 0 ) then + if ( mod(cctk_iteration,re_initialize_every) == 0 ) then re_init_control = 1 end if end if -! Get reduction handles for minimum and maximum reductions. - call CCTK_ReductionArrayHandle ( max_handle, 'maximum' ) - if ( max_handle .lt. 0 ) then - call CCTK_WARN(0,'Could not obtain a handle for maximum reduction') - end if - call CCTK_ReductionArrayHandle ( min_handle, 'minimum' ) - if ( min_handle .lt. 0 ) then - call CCTK_WARN(0,'Could not obtain a handle for minimum reduction') - end if + if ( re_init_control == 1 ) then -! Set the courant factor for the 'Euler' re-initialization scheme. - if ( CCTK_EQUALS ( re_init_int_method, 'Euler' ) ) then - hfac = one / four - end if + ! Get reduction handles for minimum and maximum reductions. + call CCTK_ReductionHandle ( max_handle, 'maximum' ) + if ( max_handle < 0 ) then + call CCTK_WARN(0,'Could not obtain a handle for maximum reduction') + end if + call CCTK_ReductionHandle ( min_handle, 'minimum' ) + if ( min_handle < 0 ) then + call CCTK_WARN(0,'Could not obtain a handle for minimum reduction') + end if -! Set the courant factor for the 'rk2' re-initialization scheme. - if ( CCTK_EQUALS ( re_init_int_method, 'rk2' ) ) then - hfac = half - end if + do l = 1, eh_number_level_sets + select case (l) + case (1:9) + j = 1 + case(10:99) + j = 2 + case(100:999) + j = 3 + case(1000:9999) + j = 4 + case(10000:) + call CCTK_WARN(0,'To many simultaneous level set functions') + end select + write (variable_name,form(j)) 'ehfinder::f[',l-1,']' + call CCTK_VarIndex ( i, variable_name ) + if ( i < 0 ) then + write (info_message,'(a9,a14,a16)') 'Variable ',variable_name, & + ' does not exist' + call CCTK_WARN(0,info_message) + end if + call CCTK_Reduce ( ierr, cctkGH, -1, max_handle, 1, & + CCTK_VARIABLE_REAL, fmax(l), 1, i) + if ( ierr .ne. 0 ) then + call CCTK_WARN(0,'Reduction of fmax failed') + end if + call CCTK_Reduce ( ierr, cctkGH, -1, min_handle, 1, & + CCTK_VARIABLE_REAL, fmin(l), 1, i) + if ( ierr .ne. 0 ) then + call CCTK_WARN(0,'Reduction of fmin failed') + end if + end do -! Set up counter for the number of iterations. Copy f into ftmp and -! fbak and eh_mask into eh_mask_bak (in case the re-initialization has -! to be undone) and find the min and max values of f. - ncalls = 0 - ftmp = f - fbak = f - eh_mask_bak = eh_mask - do l = 1, eh_number_level_sets - fminloc(l) = minval(f(:,:,:,l)) - fmaxloc(l) = maxval(f(:,:,:,l)) - end do + re_init_this_level_set = .true. + + ! If fmin and fmax have the same sign, there is no surface present and + ! re-initialization will not be performed for the given surface. + do l = 1, eh_number_level_sets + if ( fmin(l) * fmax(l) > zero ) then + re_init_this_level_set(l) = .false. + end if + end do + + ! If none of the surfaces should be re-initialized set the control variables + ! to zero. + if ( all ( .not. re_init_this_level_set ) ) then + re_init_control = 0 + info_message = 'No zero-points of the level set functions.' + info_message = trim(info_message)//' No re-initialization performed.' + call CCTK_INFO ( trim(info_message) ) + end if - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & - fmaxloc, fmax, eh_number_level_sets, & - CCTK_VARIABLE_REAL ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of fmax failed') - end if - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & - fminloc, fmin, eh_number_level_sets, & - CCTK_VARIABLE_REAL ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of fmin failed') end if - re_init_this_level_set = .true. + if ( re_init_control == 1 ) then -! If fmin and fmax have the same sign, there is no surface present and -! re-initialization will not be performed for the given surface. - do l = 1, eh_number_level_sets - if ( fmin(l) * fmax(l) .gt. zero ) then - re_init_this_level_set(l) = .false. +! If we are called more than once at the same iteration number (for example +! during the 3-level data initialization routine with Carpet) only do +! something the first time. + if ( cctk_iteration == last_called_at_iteration ) then + re_init_control = 0 + else + last_called_at_iteration = cctk_iteration end if - end do + end if -! If none of the surfaces should be re-initialized set the control variables -! to zero. - if ( all ( .not. re_init_this_level_set ) ) then - re_init_control = 0 - info_message = 'No zero-points of the level set functions. ' - info_message = trim(info_message)//'No re-initialization performed.' - call CCTK_INFO ( trim(info_message) ) + if ( re_init_control == 1 ) then + + ! Set the courant factor for the 'Euler' re-initialization scheme. + if ( CCTK_EQUALS ( re_init_int_method, 'Euler' ) ) then + hfac = one / four + end if + + ! Set the courant factor for the 'rk2' re-initialization scheme. + if ( CCTK_EQUALS ( re_init_int_method, 'rk2' ) ) then + hfac = half + end if + + niter_reinit = 0 + + call CCTK_INFO ('Re-Initialization started') end if end subroutine EHFinder_ReInitializeControl +subroutine EHFinder_ReInitializeInitialize(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + + if ( re_init_control == 1 ) then + + ! Set up counter for the number of iterations. Copy f into ftmp and + ! fbak and eh_mask into eh_mask_bak (in case the re-initialization has + ! to be undone) and find the min and max values of f. + ftmp = f + fbak = f + eh_mask_bak = eh_mask + + end if +end subroutine EHFinder_ReInitializeInitialize subroutine EHFinder_ReInitializeRK2_1(CCTK_ARGUMENTS) @@ -121,16 +175,16 @@ subroutine EHFinder_ReInitializeRK2_1(CCTK_ARGUMENTS) ! Set the step size from the courant factor and the minimum of the ! grid spacing. - h = hfac*minval(cctk_delta_space) + h = hfac*minval(CCTK_DELTA_SPACE) ! If centered finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then ! Caclulate the factor used in the finite differences used in the ! include file. - idx = half / cctk_delta_space(1) - idy = half / cctk_delta_space(2) - idz = half / cctk_delta_space(3) + idx = half / CCTK_DELTA_SPACE(1) + idy = half / CCTK_DELTA_SPACE(2) + idz = half / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -149,9 +203,9 @@ subroutine EHFinder_ReInitializeRK2_1(CCTK_ARGUMENTS) ! Caclulate the factor used in the finite differences used in the ! include file. - idx = one / cctk_delta_space(1) - idy = one / cctk_delta_space(2) - idz = one / cctk_delta_space(3) + idx = one / CCTK_DELTA_SPACE(1) + idy = one / CCTK_DELTA_SPACE(2) + idz = one / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -170,9 +224,9 @@ subroutine EHFinder_ReInitializeRK2_1(CCTK_ARGUMENTS) ! Caclulate the factor used in the finite differences used in the ! include file. - idx = half / cctk_delta_space(1) - idy = half / cctk_delta_space(2) - idz = half / cctk_delta_space(3) + idx = half / CCTK_DELTA_SPACE(1) + idy = half / CCTK_DELTA_SPACE(2) + idz = half / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -187,7 +241,7 @@ subroutine EHFinder_ReInitializeRK2_1(CCTK_ARGUMENTS) end if ! For all active grid points take half an euler step. - where ( eh_mask .ge. 0 ) + where ( eh_mask >= 0 ) sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) sftmp = - half * h * f / sqrt(f**2+one) * ( sftmp - one ) f = f + sftmp @@ -224,9 +278,9 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) ! Caclulate the factor used in the finite differences used in the ! include file. - idx = half / cctk_delta_space(1) - idy = half / cctk_delta_space(2) - idz = half / cctk_delta_space(3) + idx = half / CCTK_DELTA_SPACE(1) + idy = half / CCTK_DELTA_SPACE(2) + idz = half / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -245,9 +299,9 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) ! Caclulate the factor used in the finite differences used in the ! include file. - idx = one / cctk_delta_space(1) - idy = one / cctk_delta_space(2) - idz = one / cctk_delta_space(3) + idx = one / CCTK_DELTA_SPACE(1) + idy = one / CCTK_DELTA_SPACE(2) + idz = one / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -266,9 +320,9 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) ! Caclulate the factor used in the finite differences used in the ! include file. - idx = half / cctk_delta_space(1) - idy = half / cctk_delta_space(2) - idz = half / cctk_delta_space(3) + idx = half / CCTK_DELTA_SPACE(1) + idy = half / CCTK_DELTA_SPACE(2) + idz = half / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -284,7 +338,7 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) ! For all active grid points take a full euler step using the right hand ! side calculated at the half step. - where ( eh_mask .ge. 0 ) + where ( eh_mask >= 0 ) sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) f = ftmp + sftmp @@ -296,9 +350,9 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) ! done in this step for all active grid points on the local processor do l = 1, eh_number_level_sets maxdfloc(l) = maxval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), & - mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) + mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) >= 0) mindfloc(l) = minval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), & - mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) + mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) >= 0) end do ! Find the maximum of the changes to the level set functions done on all @@ -321,24 +375,27 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) ! If the maximum change on all processors are small enough signal that ! we are done with re-initialization by setting re_init_control = 0. - if ( all ( maxdf .lt. h*minval(cctk_delta_space)**2 ) ) then + if ( all ( maxdf < h*minval(CCTK_DELTA_SPACE)**2 ) ) then +! pugh_re_init_control = 0 re_init_control = 0 end if ! Update the number of calls counter. - ncalls = ncalls + 1 + niter_reinit = niter_reinit + 1 ! If we are finished write an info message. - if ( re_init_control .eq. 0 ) then +! if ( pugh_re_init_control == 0 ) then + if ( re_init_control == 0 ) then write(info_message,'(a35,i5,a12)') & - 'PDE re-initialization complete in ',ncalls,' iterations.' + 'PDE re-initialization complete in ',niter_reinit,' iterations.' call CCTK_INFO( trim(info_message) ) end if ! If we ran out of iterations without converging, print a level 1 warning ! message. - if ( ncalls .gt. re_init_max_iter ) then + if ( niter_reinit > re_init_max_iter ) then call CCTK_WARN(1,'Re-initialization failed to converge') +! pugh_re_init_control = 0 re_init_control = 0 end if @@ -373,16 +430,16 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) ! Set the step size from the courant factor and the minimum of the ! grid spacing. - h = hfac*minval(cctk_delta_space) + h = hfac*min(CCTK_DELTA_SPACE(1),CCTK_DELTA_SPACE(2),CCTK_DELTA_SPACE(3)) ! If centered finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then ! Caclulate the factor used in the finite differences used in the ! include file. - idx = half / cctk_delta_space(1) - idy = half / cctk_delta_space(2) - idz = half / cctk_delta_space(3) + idx = half / CCTK_DELTA_SPACE(1) + idy = half / CCTK_DELTA_SPACE(2) + idz = half / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -401,9 +458,9 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) ! Caclulate the factor used in the finite differences used in the ! include file. - idx = one / cctk_delta_space(1) - idy = one / cctk_delta_space(2) - idz = one / cctk_delta_space(3) + idx = one / CCTK_DELTA_SPACE(1) + idy = one / CCTK_DELTA_SPACE(2) + idz = one / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -422,9 +479,9 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) ! Caclulate the factor used in the finite differences used in the ! include file. - idx = half / cctk_delta_space(1) - idy = half / cctk_delta_space(2) - idz = half / cctk_delta_space(3) + idx = half / CCTK_DELTA_SPACE(1) + idy = half / CCTK_DELTA_SPACE(2) + idz = half / CCTK_DELTA_SPACE(3) ! For all level sets and all physical points... do l = 1, eh_number_level_sets @@ -439,62 +496,107 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) end if ! For all active grid points take a full euler step using. - where ( eh_mask .ge. 0 ) + where ( eh_mask >= 0 ) sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) f = f + sftmp elsewhere - sftmp = one + sftmp = zero end where -! Find the maximum and minimum of the changes to the level set functions -! done in this step for all active grid points on the local processor - do l = 1, eh_number_level_sets - maxdfloc(l) = maxval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), & - mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) - mindfloc(l) = minval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), & - mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) - end do +! Store the absulute value of the change to be used in the exit criteria, +! scaled by an appropriate factor. -! Find the maximum of the changes to the level set functions done on all -! processors - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & - maxdfloc, maxdf, eh_number_level_sets, & - CCTK_VARIABLE_REAL ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of maxdf failed') + sftmp = abs ( sftmp ) / ( h*min(CCTK_DELTA_SPACE(1),CCTK_DELTA_SPACE(2), & + CCTK_DELTA_SPACE(3))**2 ) + +end subroutine EHFinder_ReInitializeEuler + +subroutine EHFinder_ReInitializePostStep(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL, dimension(eh_number_level_sets) :: maxdf + CCTK_INT :: i,j,l + character(len=128) :: info_message, variable_name + character(15) :: form(4) = (/ '(a16,i1,a1)', '(a16,i2,a1)', & + '(a16,i3,a1)', '(a16,i4,a1)' /) + +! Get reduction handles for maximum reduction. + call CCTK_ReductionHandle ( max_handle, 'maximum' ) + if ( max_handle < 0 ) then + call CCTK_WARN(0,'Could not obtain a handle for maximum reduction') end if -! Find the minimum of the changes to the level set functions done on all -! processors - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & - mindfloc, mindf, eh_number_level_sets, & - CCTK_VARIABLE_REAL ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of mindf failed') + do l = 1, eh_number_level_sets + select case (l) + case (1:9) + j = 1 + case(10:99) + j = 2 + case(100:999) + j = 3 + case(1000:9999) + j = 4 + case(10000:) + call CCTK_WARN(0,'To many simultaneous level set functions') + end select + + write (variable_name,form(j)) 'ehfinder::sftmp[',l-1,']' + + call CCTK_VarIndex ( i, variable_name ) + + if ( i < 0 ) then + write (info_message,'(a9,a18,a16)') 'Variable ',variable_name, & + ' does not exist' + call CCTK_WARN(0,info_message) + end if + + call CCTK_Reduce ( ierr, cctkGH, -1, max_handle, 1, & + CCTK_VARIABLE_REAL, maxdf(l), 1, i) + + if ( ierr .ne. 0 ) then + call CCTK_WARN(0,'Reduction of fmax failed') + end if + + end do + + if ( re_init_verbose /= 0 ) then + write(info_message,'(a13,i5,a9,es12.5)') 'At iteration ', niter_reinit, & + ' maxdf = ', maxval(maxdf) + call CCTK_INFO ( trim(info_message) ) end if ! If the maximum change on all processors are small enough signal that -! we are done with re-initialization by setting re_init_control = 0. - if ( all ( maxdf .lt. h*minval(cctk_delta_space)**2 ) ) then - re_init_control = 0 +! we are done with re-initialization by setting *_re_init_control = 0. + if ( all ( maxdf < one ) ) then + re_init_control = 0 endif -! Update the number of calls counter. - ncalls = ncalls + 1 + if ( CCTK_IsThornACtive ( 'PUGH' ) /= 0 ) then +! Update the number of calls counter. + niter_reinit = niter_reinit + 1 + end if ! If we are finished write an info message. - if ( re_init_control .eq. 0 ) then +! if ( pugh_re_init_control == 0 .and. carpet_re_init_control == 0 ) then + if ( re_init_control == 0 ) then write(info_message,'(a35,i5,a12)') & - 'PDE re-initialization complete in ',ncalls,' iterations.' + 'PDE re-initialization complete in ',niter_reinit,' iterations.' call CCTK_INFO( trim(info_message) ) end if ! If we ran out of iterations without converging, print a level 1 warning ! message. - if ( ncalls .gt. re_init_max_iter ) then + if ( niter_reinit > re_init_max_iter ) then call CCTK_WARN(1,'Re-initialization failed to converge') re_init_control = 0 end if -end subroutine EHFinder_ReInitializeEuler +end subroutine EHFinder_ReInitializePostStep -- cgit v1.2.3