| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module cut_detector | ||
| 2 | use, intrinsic :: iso_fortran_env, only: dp => real64 | ||
| 3 | use util, only: twopi | ||
| 4 | use simple, only: tstep | ||
| 5 | use orbit_symplectic, only: symplectic_integrator_t | ||
| 6 | use field_can_mod, only: field_can_t | ||
| 7 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 8 | use params, only: debug | ||
| 9 | #endif | ||
| 10 | |||
| 11 | implicit none | ||
| 12 | |||
| 13 | save | ||
| 14 | |||
| 15 | integer, parameter :: n_tip_vars = 6 | ||
| 16 | integer, parameter :: nplagr = 6 | ||
| 17 | integer, parameter :: nder = 0 | ||
| 18 | |||
| 19 | public | ||
| 20 | |||
| 21 | type :: cut_detector_t | ||
| 22 | real(dp) :: fper ! field period | ||
| 23 | |||
| 24 | ! for Poincare cuts | ||
| 25 | real(dp) :: alam_prev, par_inv | ||
| 26 | integer :: iper, itip, kper | ||
| 27 | |||
| 28 | real(dp) :: orb_sten(6,nplagr), coef(0:nder,nplagr) | ||
| 29 | integer :: ipoi(nplagr) | ||
| 30 | end type cut_detector_t | ||
| 31 | |||
| 32 | contains | ||
| 33 | |||
| 34 | ✗ | subroutine init(self, fper, z) | |
| 35 | type(cut_detector_t) :: self | ||
| 36 | real(dp), intent(in) :: fper | ||
| 37 | real(dp), intent(in) :: z(:) | ||
| 38 | integer :: i | ||
| 39 | |||
| 40 | ✗ | self%fper = fper | |
| 41 | |||
| 42 | !--------------------------------------------------------------------------- | ||
| 43 | ! Prepare calculation of orbit tip by interpolation and buffer for Poincare plot: | ||
| 44 | |||
| 45 | ✗ | do i=1,nplagr | |
| 46 | ✗ | self%ipoi(i)=i | |
| 47 | enddo | ||
| 48 | |||
| 49 | !-------------------------------- | ||
| 50 | ! Initialize tip detector | ||
| 51 | |||
| 52 | ✗ | self%itip=nplagr/2+1 | |
| 53 | ✗ | self%alam_prev=z(5) | |
| 54 | |||
| 55 | ! End initialize tip detector | ||
| 56 | !-------------------------------- | ||
| 57 | ! Initialize period crossing detector | ||
| 58 | |||
| 59 | ✗ | self%iper=nplagr/2+1 | |
| 60 | ✗ | self%kper=int(z(3)/self%fper) | |
| 61 | |||
| 62 | ! End initialize period crossing detector | ||
| 63 | !-------------------------------- | ||
| 64 | ✗ | self%par_inv = 0.0d0 | |
| 65 | ! | ||
| 66 | |||
| 67 | ! End prepare calculation of orbit tip by interpolation | ||
| 68 | !-------------------------------------------------------------------------- | ||
| 69 | ✗ | end subroutine init | |
| 70 | |||
| 71 | ✗ | subroutine trace_to_cut(self, si, f, z, var_cut, cut_type, ierr) | |
| 72 | use plag_coeff_sub, only : plag_coeff | ||
| 73 | |||
| 74 | type(cut_detector_t) :: self | ||
| 75 | type(symplectic_integrator_t) :: si | ||
| 76 | type(field_can_t) :: f | ||
| 77 | |||
| 78 | real(dp), intent(inout) :: z(:) | ||
| 79 | ! variables to evaluate at tip: z(1..5), par_inv | ||
| 80 | real(dp), dimension(:), intent(inout) :: var_cut | ||
| 81 | integer, intent(out) :: cut_type | ||
| 82 | integer, intent(out) :: ierr | ||
| 83 | |||
| 84 | integer, parameter :: nstep_max = 1000000000 | ||
| 85 | integer :: i | ||
| 86 | real(dp) :: phiper = 0.0d0 | ||
| 87 | |||
| 88 | ✗ | do i = 1, nplagr | |
| 89 | ✗ | call tstep(si, f, z, ierr) | |
| 90 | ✗ | if (ierr.ne.0) exit | |
| 91 | |||
| 92 | ✗ | self%par_inv = self%par_inv + z(5)**2 ! parallel adiabatic invariant | |
| 93 | |||
| 94 | ✗ | self%orb_sten(1:5, i) = z | |
| 95 | ✗ | self%orb_sten(6, i) = self%par_inv | |
| 96 | |||
| 97 | !------------------------------------------------------------------------- | ||
| 98 | ! Tip detection and interpolation | ||
| 99 | |||
| 100 | ✗ | if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) then | |
| 101 | ✗ | self%itip=0 !<=tip has been passed | |
| 102 | end if | ||
| 103 | ✗ | self%itip=self%itip+1 | |
| 104 | ✗ | self%alam_prev=z(5) | |
| 105 | |||
| 106 | !<=use only initialized stencil | ||
| 107 | if(i.gt.nplagr .and. self%itip.eq.nplagr/2) then | ||
| 108 | cut_type = 0 | ||
| 109 | call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), & | ||
| 110 | self%coef) | ||
| 111 | var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) | ||
| 112 | var_cut(2) = modulo(var_cut(2), twopi) | ||
| 113 | var_cut(3) = modulo(var_cut(3), twopi) | ||
| 114 | self%par_inv = self%par_inv - var_cut(6) | ||
| 115 | return | ||
| 116 | end if | ||
| 117 | |||
| 118 | ! End tip detection and interpolation | ||
| 119 | |||
| 120 | !------------------------------------------------------------------------- | ||
| 121 | ! Periodic boundary footprint detection and interpolation | ||
| 122 | |||
| 123 | ✗ | if(z(3).gt.dble(self%kper+1)*self%fper) then | |
| 124 | ✗ | self%iper=0 !<=periodic boundary has been passed | |
| 125 | ✗ | phiper=dble(self%kper+1)*self%fper | |
| 126 | ✗ | self%kper=self%kper+1 | |
| 127 | ✗ | elseif(z(3).lt.dble(self%kper)*self%fper) then | |
| 128 | ✗ | self%iper=0 !<=periodic boundary has been passed | |
| 129 | ✗ | phiper=dble(self%kper)*self%fper | |
| 130 | ✗ | self%kper=self%kper-1 | |
| 131 | endif | ||
| 132 | ✗ | self%iper=self%iper+1 | |
| 133 | |||
| 134 | !<=use only initialized stencil | ||
| 135 | ✗ | if(i.gt.nplagr .and. self%iper.eq.nplagr/2) then | |
| 136 | cut_type = 1 | ||
| 137 | !<=stencil around periodic boundary is complete, interpolate | ||
| 138 | call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - & | ||
| 139 | phiper, self%coef) | ||
| 140 | var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) | ||
| 141 | var_cut(2) = modulo(var_cut(2), twopi) | ||
| 142 | var_cut(3) = modulo(var_cut(3), twopi) | ||
| 143 | return | ||
| 144 | end if | ||
| 145 | |||
| 146 | ! End periodic boundary footprint detection and interpolation | ||
| 147 | !------------------------------------------------------------------------- | ||
| 148 | |||
| 149 | end do | ||
| 150 | |||
| 151 | ✗ | do i = nplagr + 1, nstep_max | |
| 152 | ✗ | call tstep(si, f, z, ierr) | |
| 153 | ✗ | if (ierr.ne.0) exit | |
| 154 | |||
| 155 | ✗ | self%par_inv = self%par_inv + z(5)**2 ! parallel adiabatic invariant | |
| 156 | |||
| 157 | ✗ | self%orb_sten(1:5, self%ipoi(1)) = z | |
| 158 | ✗ | self%orb_sten(6, self%ipoi(1)) = self%par_inv | |
| 159 | ✗ | self%ipoi = cshift(self%ipoi, 1) | |
| 160 | |||
| 161 | !------------------------------------------------------------------------- | ||
| 162 | ! Tip detection and interpolation | ||
| 163 | |||
| 164 | ✗ | if(self%alam_prev.lt.0.d0.and.z(5).gt.0.d0) then | |
| 165 | ✗ | self%itip=0 !<=tip has been passed | |
| 166 | end if | ||
| 167 | ✗ | self%itip=self%itip+1 | |
| 168 | ✗ | self%alam_prev=z(5) | |
| 169 | |||
| 170 | ✗ | if(self%itip.eq.nplagr/2) then | |
| 171 | ✗ | cut_type = 0 | |
| 172 | call plag_coeff(nplagr, nder, 0d0, self%orb_sten(5, self%ipoi), & | ||
| 173 | ✗ | self%coef) | |
| 174 | ✗ | var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) | |
| 175 | ✗ | var_cut(2) = modulo(var_cut(2), twopi) | |
| 176 | ✗ | var_cut(3) = modulo(var_cut(3), twopi) | |
| 177 | ✗ | self%par_inv = self%par_inv - var_cut(6) | |
| 178 | return | ||
| 179 | end if | ||
| 180 | |||
| 181 | ! End tip detection and interpolation | ||
| 182 | |||
| 183 | !------------------------------------------------------------------------- | ||
| 184 | ! Periodic boundary footprint detection and interpolation | ||
| 185 | |||
| 186 | ✗ | if(z(3).gt.dble(self%kper+1)*self%fper) then | |
| 187 | ✗ | self%iper=0 !<=periodic boundary has been passed | |
| 188 | ✗ | phiper=dble(self%kper+1)*self%fper | |
| 189 | ✗ | self%kper=self%kper+1 | |
| 190 | ✗ | elseif(z(3).lt.dble(self%kper)*self%fper) then | |
| 191 | ✗ | self%iper=0 !<=periodic boundary has been passed | |
| 192 | ✗ | phiper=dble(self%kper)*self%fper | |
| 193 | ✗ | self%kper=self%kper-1 | |
| 194 | endif | ||
| 195 | ✗ | self%iper=self%iper+1 | |
| 196 | |||
| 197 | ✗ | if(self%iper.eq.nplagr/2) then | |
| 198 | ✗ | cut_type = 1 | |
| 199 | !<=stencil around periodic boundary is complete, interpolate | ||
| 200 | call plag_coeff(nplagr, nder, 0d0, self%orb_sten(3, self%ipoi) - & | ||
| 201 | ✗ | phiper, self%coef) | |
| 202 | ✗ | var_cut = matmul(self%orb_sten(:, self%ipoi), self%coef(0,:)) | |
| 203 | ✗ | var_cut(2) = modulo(var_cut(2), twopi) | |
| 204 | ✗ | var_cut(3) = modulo(var_cut(3), twopi) | |
| 205 | return | ||
| 206 | end if | ||
| 207 | |||
| 208 | ! End periodic boundary footprint detection and interpolation | ||
| 209 | !------------------------------------------------------------------------- | ||
| 210 | |||
| 211 | end do | ||
| 212 | end subroutine trace_to_cut | ||
| 213 | |||
| 214 | ✗ | subroutine fract_dimension(ntr,rt,fraction) | |
| 215 | |||
| 216 | integer, parameter :: iunit=1003 | ||
| 217 | integer :: itr,ntr,ngrid,nrefine,irefine,kr,kt,nboxes | ||
| 218 | real(dp) :: fraction,rmax,rmin,tmax,tmin,hr,ht | ||
| 219 | real(dp), dimension(2,ntr) :: rt!0, rt | ||
| 220 | ✗ | logical, dimension(:,:), allocatable :: free | |
| 221 | |||
| 222 | ! TODO: check if this works better on tips only | ||
| 223 | ! rt(1,:) = rt0(1,:)*cos(rt0(2,:)) | ||
| 224 | ! rt(2,:) = rt0(1,:)*sin(rt0(2,:)) | ||
| 225 | |||
| 226 | ✗ | rmin=minval(rt(1,:)) | |
| 227 | ✗ | rmax=maxval(rt(1,:)) | |
| 228 | ✗ | tmin=minval(rt(2,:)) | |
| 229 | ✗ | tmax=maxval(rt(2,:)) | |
| 230 | |||
| 231 | ✗ | nrefine=int(log(dble(ntr))/log(4.d0)) | |
| 232 | |||
| 233 | ✗ | ngrid=1 | |
| 234 | ✗ | nrefine=nrefine+3 !<=add 3 for curiousity | |
| 235 | ✗ | do irefine=1,nrefine | |
| 236 | ✗ | ngrid=ngrid*2 | |
| 237 | ✗ | !$omp critical | |
| 238 | ✗ | allocate(free(0:ngrid,0:ngrid)) | |
| 239 | !$omp end critical | ||
| 240 | ✗ | free=.true. | |
| 241 | ✗ | hr=(rmax-rmin)/dble(ngrid) | |
| 242 | ✗ | ht=(tmax-tmin)/dble(ngrid) | |
| 243 | ✗ | nboxes=0 | |
| 244 | ✗ | do itr=1,ntr | |
| 245 | ✗ | kr=int((rt(1,itr)-rmin)/hr) | |
| 246 | ✗ | kr=min(ngrid-1,max(0,kr)) | |
| 247 | ✗ | kt=int((rt(2,itr)-tmin)/ht) | |
| 248 | ✗ | kt=min(ngrid-1,max(0,kt)) | |
| 249 | ✗ | if(free(kr,kt)) then | |
| 250 | ✗ | free(kr,kt)=.false. | |
| 251 | ✗ | nboxes=nboxes+1 | |
| 252 | endif | ||
| 253 | enddo | ||
| 254 | ✗ | !$omp critical | |
| 255 | ✗ | deallocate(free) | |
| 256 | !$omp end critical | ||
| 257 | #ifdef SIMPLE_ENABLE_DEBUG_OUTPUT | ||
| 258 | if(debug) then | ||
| 259 | !$omp critical | ||
| 260 | ! | ||
| 261 | ! Right now criterion for regular is at nboxes/ngrid**2 < 0.2 . | ||
| 262 | ! For fractal dimension d = log(nboxes)/log(ngrid) this means | ||
| 263 | ! d_thresh = 2 + log(0.2)/log(ngrid) | ||
| 264 | ! | ||
| 265 | write(iunit,*) irefine, nboxes, ngrid, dble(nboxes)/dble(ngrid**2), 0.2d0,& | ||
| 266 | log(1d0*nboxes)/log(1d0*ngrid), 2d0 + log(0.2d0)/log(1d0*ngrid) | ||
| 267 | !$omp end critical | ||
| 268 | end if | ||
| 269 | #endif | ||
| 270 | ✗ | if(irefine.eq.nrefine-3) fraction=dble(nboxes)/dble(ngrid**2) | |
| 271 | enddo | ||
| 272 | ✗ | close(iunit) | |
| 273 | |||
| 274 | ✗ | end subroutine fract_dimension | |
| 275 | |||
| 276 | ✗ | end module cut_detector | |
| 277 |