From 5167aed638fffc8014bf1fb4a46605574d01e31b Mon Sep 17 00:00:00 2001 From: diener Date: Mon, 1 Sep 2003 12:32:13 +0000 Subject: EHFinder_ReParametrize.F90 is now EHFinder_ReInitialize.F90 to be consistent with the exclusive use of reinitialization instead of reparametrization. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@132 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_ReInitialize.F90 | 349 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 349 insertions(+) create mode 100644 src/EHFinder_ReInitialize.F90 diff --git a/src/EHFinder_ReInitialize.F90 b/src/EHFinder_ReInitialize.F90 new file mode 100644 index 0000000..70b8f47 --- /dev/null +++ b/src/EHFinder_ReInitialize.F90 @@ -0,0 +1,349 @@ +! Re-Initialization of the level set function +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +! Control routine for re-initialization of the level set function. + +subroutine EHFinder_ReInitializeControl(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_REAL, dimension(eh_number_level_sets) :: fmin, fmax, fminloc, fmaxloc + CCTK_INT :: i, j, k, l + +! 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 +! variable to 1. + if ( re_initialize_every .gt. 0 ) then + if ( mod(cctk_iteration,re_initialize_every) .eq. 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 + +! Set the courant factor for the 'Euler' re-parametrization 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. + 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 +! 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)) + end do + + 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 + + reparam_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. + do l = 1, eh_number_level_sets + if ( fmin(l) * fmax(l) .gt. zero ) then + reparam_this_level_set(l) = .false. + end if + end do + +! If none of the surfaces should be reparametrized set the control variables +! to zero. + if ( all ( .not. reparam_this_level_set ) ) then + re_init_control = 0 + call CCTK_INFO ( 'No zero-points of the level set functions. No re-parametrization performed.' ) + end if + +end subroutine EHFinder_ReInitializeControl + + +subroutine EHFinder_ReInitializeRK2_1(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, l + CCTK_REAL :: idx, idy, idz + CCTK_REAL :: al, ar, bl, br, cl, cr + CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus + CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus + + h = hfac*minval(cctk_delta_space) + + if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then + +# include "include/centered_second.h" + + end if + + if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then + +# include "include/upwind_first.h" + + end if + + if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then + +# include "include/upwind_second.h" + + end if + + 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 + end where + +end subroutine EHFinder_ReInitializeRK2_1 + + +subroutine EHFinder_ReInitializeRK2_2(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, l + CCTK_REAL :: idx, idy, idz, maxf + CCTK_REAL, dimension(eh_number_level_sets) :: maxdfloc, maxdf, & + mindfloc, mindf + CCTK_REAL :: al, ar, bl, br, cl, cr + CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus + CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus + character(len=128) :: info_message + + + if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then + +# include "include/centered_second.h" + + end if + + if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then + +# include "include/upwind_first.h" + + end if + + if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then + +# include "include/upwind_second.h" + + end if + + 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 + end where + + 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 + + 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 + + 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') + end if + + if ( all ( maxdf .lt. h*minval(cctk_delta_space)**2 ) ) then + re_init_control = 0 + end if + + ncalls = ncalls + 1 + +! print*,ncalls,maxdf,h*minval(cctk_delta_space)**2 + + if ( re_init_control .eq. 0 ) then + write(info_message,'(a35,i5,a12)') & + 'PDE re-parametrization complete in ',ncalls,' iterations.' + call CCTK_INFO(info_message) + end if + + if ( ncalls .gt. re_init_max_iter ) then + call CCTK_WARN(1,'Re-parametrization failed to converge') + re_init_control = 0 + end if + + ftmp = f + +end subroutine EHFinder_ReInitializeRK2_2 + + +subroutine EHFinder_ReInitializeEuler(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, l + CCTK_REAL :: idx, idy, idz, maxf + CCTK_REAL, dimension(eh_number_level_sets) :: maxdfloc, maxdf, & + mindfloc, mindf + CCTK_REAL :: al, ar, bl, br, cl, cr + CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus + CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus + character(len=128) :: info_message + + h = hfac*minval(cctk_delta_space) + + if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then + +# include "include/centered_second.h" + + end if + + if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then + +# include "include/upwind_first.h" + + end if + + if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then + +# include "include/upwind_second.h" + + end if + + where ( eh_mask .ge. 0 ) + sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) + sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) + f = f + sftmp + 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 + + 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 + + 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 + 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') + end if + + if ( all ( maxdf .lt. h*minval(cctk_delta_space)**2 ) ) then + re_init_control = 0 + endif + + ncalls = ncalls + 1 + +! print*,ncalls,maxdf,h*minval(cctk_delta_space)**2 + if ( re_init_control .eq. 0 ) then + write(info_message,'(a35,i5,a12)') & + 'PDE re-parametrization complete in ',ncalls,' iterations.' + call CCTK_INFO(info_message) + end if + + if ( ncalls .gt. re_init_max_iter ) then + call CCTK_WARN(1,'Re-parametrization failed to converge') + re_init_control = 0 + end if + +! print*,ncalls +! print*,maxdf +! print*,mindf + +end subroutine EHFinder_ReInitializeEuler -- cgit v1.2.3