aboutsummaryrefslogtreecommitdiff
path: root/src/EHFinder_Check.F90
blob: 4a6e07aca6259cda8a728608e22e54a5a5e55f36 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
!option opt(o(0)))

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"

subroutine EHFinder_ReParametrize_Check(CCTK_ARGUMENTS)

  use EHFinder_mod

  implicit none

  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_FUNCTIONS

  CCTK_INT :: i, j, k, l
  CCTK_REAL :: fmin_extremum, fmin_extremum_loc, fmax_loc, fmin_bound_loc
  character(len=128) :: info_message
  logical :: check_reparam
  CCTK_INT, dimension(3) :: ex_loc

  check_reparam = .false.
  reparam_undone = .false.
  if ( ( CCTK_EQUALS ( re_param_method, 'approx' ) ) .or. &
       ( CCTK_EQUALS ( re_param_method, 'mixed' ) ) ) then
    if ( reparametrize_every_approx .gt. 0 ) then
      if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) then
        check_reparam = .true.
      end if
    end if
  end if
  if ( ( CCTK_EQUALS ( re_param_method, 'pde' ) ) .or. &
       ( CCTK_EQUALS ( re_param_method, 'mixed' ) ) ) then
    if ( reparametrize_every_pde .gt. 0 ) then
      if ( mod(cctk_iteration,reparametrize_every_pde) .eq. 0 ) then
        check_reparam = .true.
      end if 
    end if 
  end if 
    
#include "include/physical_part2.h"

  if ( check_reparam ) then
    do l = 1, eh_number_level_sets
      fmin_bound_loc = huge
      if ( cctk_lbnd(1) .eq. 0 .and. .not.symx ) then
        fmin_bound_loc = min ( fmin_bound_loc, minval ( f(1,:,:,l) ) )
      end if
      if ( cctk_ubnd(1) .eq. cctk_gsh(1)-1 ) then
        fmin_bound_loc = min ( fmin_bound_loc, minval ( f(nx,:,:,l) ) )
      end if
      if ( cctk_lbnd(2) .eq. 0 .and. .not.symy ) then
        fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,1,:,l) ) )
      end if
      if ( cctk_ubnd(2) .eq. cctk_gsh(2)-1 ) then
        fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,ny,:,l) ) )
      end if
      if ( cctk_lbnd(3) .eq. 0 .and. .not.symz ) then
        fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,:,1,l) ) )
      end if
      if ( cctk_ubnd(3) .eq. cctk_gsh(3)-1 ) then
        fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,:,nz,l) ) )
      end if
  !    print*,fmin_bound_loc
    
      fmax_loc = zero
  
      fmin_extremum_loc = ex_value
      do k = kzl, kzr
        do j = jyl, jyr
          do i = ixl, ixr
            fmax_loc = max ( f(i,j,k,l), fmax_loc )
            if ( f(i,j,k,l) .lt. zero ) then
              if ( ( (f(i,j,k,l)-f(i-1,j,k,l))*(f(i+1,j,k,l)-f(i,j,k,l)).le.zero ).and. &
                   ( (f(i,j,k,l)-f(i,j-1,k,l))*(f(i,j+1,k,l)-f(i,j,k,l)).le.zero ).and. &
                   ( (f(i,j,k,l)-f(i,j,k-1,l))*(f(i,j,k+1,l)-f(i,j,k,l)).le.zero ) ) then
                fmin_extremum_loc = max ( f(i,j,k,l), fmin_extremum_loc )
              end if
            end if
          end do
        end do
      end do
  
      call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, &
                         fmin_extremum_loc, fmin_extremum, CCTK_VARIABLE_REAL )
  
      if ( ierr .ne. 0 ) then
        call CCTK_WARN(0,'Reduction of fmax_extremum failed')
      end if
  
      call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, &
                         fmin_bound_loc, fmin_bound, CCTK_VARIABLE_REAL )
  
      if ( ierr .ne. 0 ) then
        call CCTK_WARN(0,'Reduction of fmax failed')
      end if
  
      write(info_message,'(a39,i3,a6,f12.5)') &
          'Minimum f at the extrema for level set ', l, ' is : ',fmin_extremum
      call CCTK_INFO(info_message)
  
      write(info_message,'(a45,i3,a6,f12.5)') &
          'Limit for centered differences for level set ', l, ' is : ', &
                      fmin_bound - three*maxval ( cctk_delta_space )
      call CCTK_INFO(info_message)
  
      if ( reparam_undo .gt. 0 ) then
        if ( fmin_extremum .gt. -two * minval(cctk_delta_space) ) then
          write(info_message,'(a45,i3,a6,f12.5)') &
               'Re-parametrization of level set ', l, &
               'undone due to imminent pinch-off'
          call CCTK_INFO(info_message)
          f(:,:,:,l) = fbak(:,:,:,l)
          eh_mask(:,:,:,l) = eh_mask_bak(:,:,:,l)
          reparam_undone(l) = .true.
        end if
      end if
    end do
  end if

end subroutine EHFinder_ReParametrize_Check