! Re-Initialization of the level set function
! $Header$

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.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
  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-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
      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-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

! 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

  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 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.
    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

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

! 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

!   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

!   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 = zero
  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

#include "include/physical_part.h"

! 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)

!   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

!   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 = 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

! 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 )
  if ( ierr .ne. 0 ) then
    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

! If we are finished write an info message.
  if ( re_init_control .eq. 0 ) then
    write(info_message,'(a35,i5,a12)') &
        'PDE re-initialization complete in ',ncalls,' 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
    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


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

! 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

!   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

!   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 )
    f = f + sftmp
  elsewhere
    sftmp = one
  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

! 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 )
  if ( ierr .ne. 0 ) then
    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

! If we are finished write an info message.
  if ( re_init_control .eq. 0 ) then
    write(info_message,'(a35,i5,a12)') &
        'PDE re-initialization complete in ',ncalls,' 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
    call CCTK_WARN(1,'Re-initialization failed to converge')
    re_init_control = 0
  end if

end subroutine EHFinder_ReInitializeEuler