Moved EHFinder_ReInitialize2 to EHFinder_ReInitialize.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@141 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files 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