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

1 files changed, 349 insertions, 0 deletions

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