From b8e4ace240c743b05228d2a57ef317059ed6d525 Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 30 Sep 2003 08:57:28 +0000 Subject: Moved EHFinder_ReInitialize2 to EHFinder_ReInitialize. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@141 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_ReInitialize.F90 | 271 ++++++++++++++++++++++++++++++++---------- 1 file changed, 211 insertions(+), 60 deletions(-) diff --git a/src/EHFinder_ReInitialize.F90 b/src/EHFinder_ReInitialize.F90 index 70b8f47..3ff473a 100644 --- a/src/EHFinder_ReInitialize.F90 +++ b/src/EHFinder_ReInitialize.F90 @@ -19,12 +19,13 @@ subroutine EHFinder_ReInitializeControl(CCTK_ARGUMENTS) CCTK_REAL, dimension(eh_number_level_sets) :: fmin, fmax, fminloc, fmaxloc CCTK_INT :: i, j, k, l + character(len=128) :: info_message ! Initialize the control variables. The default means no re-initialization. re_init_control = 0 -! Check if it is time for a pde re-parametrization and set the control +! 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 @@ -42,26 +43,23 @@ subroutine EHFinder_ReInitializeControl(CCTK_ARGUMENTS) call CCTK_WARN(0,'Could not obtain a handle for minimum reduction') end if -! Set the courant factor for the 'Euler' re-parametrization scheme. +! Set the courant factor for the 'Euler' re-initialization scheme. if ( CCTK_EQUALS ( re_init_int_method, 'Euler' ) ) then -! hfac = one / eight hfac = one / four end if -! Set the courant factor for the 'rk2' re-parametrization scheme. +! Set the courant factor for the 'rk2' re-initialization scheme. if ( CCTK_EQUALS ( re_init_int_method, 'rk2' ) ) then hfac = half -! hfac = one / four end if ! 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-parametrization has +! 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 -! rep_mask = 0 do l = 1, eh_number_level_sets fminloc(l) = minval(f(:,:,:,l)) fmaxloc(l) = maxval(f(:,:,:,l)) @@ -80,21 +78,23 @@ subroutine EHFinder_ReInitializeControl(CCTK_ARGUMENTS) call CCTK_WARN(0,'Reduction of fmin failed') end if - reparam_this_level_set = .true. + re_init_this_level_set = .true. ! If fmin and fmax have the same sign, there is no surface present and -! re-parametrization will not be performed for the given surface. +! 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 - reparam_this_level_set(l) = .false. + re_init_this_level_set(l) = .false. end if end do -! If none of the surfaces should be reparametrized set the control variables +! If none of the surfaces should be re-initialized set the control variables ! to zero. - if ( all ( .not. reparam_this_level_set ) ) then + if ( all ( .not. re_init_this_level_set ) ) then re_init_control = 0 - call CCTK_INFO ( 'No zero-points of the level set functions. No re-parametrization performed.' ) + info_message = 'No zero-points of the level set functions. ' + info_message = info_message//'No re-initialization performed.' + call CCTK_INFO ( info_message ) end if end subroutine EHFinder_ReInitializeControl @@ -116,32 +116,83 @@ subroutine EHFinder_ReInitializeRK2_1(CCTK_ARGUMENTS) CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus +! Find the physical part of the computational domain. +#include "include/physical_part.h" + +! Set the step size from the courant factor and the minimum of the +! grid spacing. h = hfac*minval(cctk_delta_space) +! If centered finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then -# include "include/centered_second.h" - +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr +# include "include/centered_second2.h" + end do + end do + end do + end do end if +! If first order upwinded finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = 1, nz + do j = 1, ny + do i = 1, nx # include "include/upwind_first.h" - + end do + end do + end do + end do end if +! If second order upwinded finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then -# include "include/upwind_second.h" - +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr +# include "include/upwind_second2.h" + end do + end do + end do + end do end if +! For all active grid points take half an euler step. where ( eh_mask .ge. 0 ) sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) sftmp = - half * h * f / sqrt(f**2+one) * ( sftmp - one ) f = f + sftmp elsewhere - sftmp = one + sftmp = zero end where end subroutine EHFinder_ReInitializeRK2_1 @@ -166,34 +217,83 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus character(len=128) :: info_message +#include "include/physical_part.h" +! If centered finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then -# include "include/centered_second.h" - +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr +# include "include/centered_second2.h" + end do + end do + end do + end do end if +! If first order upwinded finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr # include "include/upwind_first.h" - + end do + end do + end do + end do end if +! If second order upwinded finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then -# include "include/upwind_second.h" - +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr +# include "include/upwind_second2.h" + end do + end do + end do + end do end if +! 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 ) - sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) f = ftmp + 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) @@ -201,6 +301,8 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) end do +! 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 ) @@ -208,6 +310,8 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) call CCTK_WARN(0,'Reduction of maxdf failed') 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 ) @@ -215,25 +319,31 @@ subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) call CCTK_WARN(0,'Reduction of mindf failed') 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 end if +! Update the number of calls counter. ncalls = ncalls + 1 -! print*,ncalls,maxdf,h*minval(cctk_delta_space)**2 - +! If we are finished write an info message. if ( re_init_control .eq. 0 ) then write(info_message,'(a35,i5,a12)') & - 'PDE re-parametrization complete in ',ncalls,' iterations.' + 'PDE re-initialization complete in ',ncalls,' iterations.' call CCTK_INFO(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 - call CCTK_WARN(1,'Re-parametrization failed to converge') + call CCTK_WARN(1,'Re-initialization failed to converge') re_init_control = 0 end if +! Copy the new values for the level set function into the temporary variable, +! so we are ready for the next step. ftmp = f end subroutine EHFinder_ReInitializeRK2_2 @@ -258,26 +368,77 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus character(len=128) :: info_message +! Find the physical part of the computational domain. +#include "include/physical_part.h" + +! Set the step size from the courant factor and the minimum of the +! grid spacing. h = hfac*minval(cctk_delta_space) +! If centered finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then -# include "include/centered_second.h" - +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr +# include "include/centered_second2.h" + end do + end do + end do + end do end if +! If first order upwinded finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr # include "include/upwind_first.h" - + end do + end do + end do + end do end if +! If second order upwinded finite differences are required... if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then -# include "include/upwind_second.h" - +! 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) + +! For all level sets and all physical points... + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr +# include "include/upwind_second2.h" + end do + end do + end do + end do end if +! For all active grid points take a full euler step using. where ( eh_mask .ge. 0 ) sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) @@ -285,25 +446,9 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) elsewhere sftmp = one end where -! print*,f(2:4,2:4,2:4,5) -! print*,sftmp(3,3,3,5) -! stop - -! do l = 1, eh_number_level_sets -! maxdfloc = zero -! mindfloc = 1.0d23 -! do k = kzl, kzr -! do j = jyl, jyr -! do i = ixl, ixr -! if ( eh_mask(i,j,k,l) .ge. 0 ) then -! maxdfloc = max ( maxdfloc, abs(sftmp(i,j,k,l)) ) -! mindfloc = min ( mindfloc, abs(sftmp(i,j,k,l)) ) -! end if -! end do -! end do -! end do -! end do +! 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) @@ -311,12 +456,17 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) end do +! 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') 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 ) @@ -324,26 +474,27 @@ subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) call CCTK_WARN(0,'Reduction of mindf failed') 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 endif +! Update the number of calls counter. ncalls = ncalls + 1 -! print*,ncalls,maxdf,h*minval(cctk_delta_space)**2 +! If we are finished write an info message. if ( re_init_control .eq. 0 ) then write(info_message,'(a35,i5,a12)') & - 'PDE re-parametrization complete in ',ncalls,' iterations.' + 'PDE re-initialization complete in ',ncalls,' iterations.' call CCTK_INFO(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 - call CCTK_WARN(1,'Re-parametrization failed to converge') + call CCTK_WARN(1,'Re-initialization failed to converge') re_init_control = 0 end if -! print*,ncalls -! print*,maxdf -! print*,mindf - end subroutine EHFinder_ReInitializeEuler -- cgit v1.2.3