| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | ! | ||
| 2 | module bminmax_mod | ||
| 3 | implicit none | ||
| 4 | |||
| 5 | ! Define real(dp) kind parameter | ||
| 6 | integer, parameter :: dp = kind(1.0d0) | ||
| 7 | ! | ||
| 8 | logical :: prop=.true. | ||
| 9 | ! | ||
| 10 | integer, parameter :: nsbmnx=100 | ||
| 11 | real(dp) :: hsbmnx | ||
| 12 | real(dp), dimension(0:nsbmnx) :: bmin_arr,bmax_arr | ||
| 13 | ! | ||
| 14 | end module bminmax_mod | ||
| 15 | ! | ||
| 16 | module find_bminmax_sub | ||
| 17 | implicit none | ||
| 18 | |||
| 19 | ! Define real(dp) kind parameter | ||
| 20 | integer, parameter :: dp = kind(1.0d0) | ||
| 21 | |||
| 22 | contains | ||
| 23 | |||
| 24 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 25 | ! | ||
| 26 | 808 | subroutine find_bminmax(s,bmin,bmax) | |
| 27 | ! | ||
| 28 | use new_vmec_stuff_mod, only : nper | ||
| 29 | ! | ||
| 30 | implicit none | ||
| 31 | ! | ||
| 32 | integer, parameter :: np=100, nt=100 | ||
| 33 | ! | ||
| 34 | integer :: it,ip,iter | ||
| 35 | real(dp) :: s,bmin,bmax,bmod,hp,ht,twopi | ||
| 36 | real(dp) :: tmin,tmax,pmin,pmax | ||
| 37 | real(dp), dimension(3) :: x,bder | ||
| 38 | ! | ||
| 39 | 808 | twopi=8.d0*atan(1.d0) | |
| 40 | 808 | ht=twopi/dble(nt) | |
| 41 | 808 | hp=twopi/dble(nt*nper) | |
| 42 | 808 | x(1)=s | |
| 43 |
2/2✓ Branch 0 taken 1616 times.
✓ Branch 1 taken 808 times.
|
2424 | x(2:3)=0.d0 |
| 44 | ! | ||
| 45 | 808 | call getbmod | |
| 46 | ! | ||
| 47 | 808 | tmin=0.d0 | |
| 48 | 808 | tmax=0.d0 | |
| 49 | 808 | pmin=0.d0 | |
| 50 | 808 | pmax=0.d0 | |
| 51 | 808 | bmin=bmod | |
| 52 | 808 | bmax=bmod | |
| 53 | ! | ||
| 54 |
2/2✓ Branch 0 taken 80800 times.
✓ Branch 1 taken 808 times.
|
81608 | do it=1,nt |
| 55 | 80800 | x(2)=ht*dble(it) | |
| 56 |
2/2✓ Branch 0 taken 8080000 times.
✓ Branch 1 taken 80800 times.
|
8161608 | do ip=1,np |
| 57 | 8080000 | x(3)=hp*dble(ip) | |
| 58 | ! | ||
| 59 | 8080000 | call getbmod | |
| 60 | ! | ||
| 61 |
2/2✓ Branch 0 taken 357206 times.
✓ Branch 1 taken 7722794 times.
|
8160800 | if(bmod.gt.bmax) then |
| 62 | 357206 | bmax=bmod | |
| 63 | 357206 | tmax=x(2) | |
| 64 | 357206 | pmax=x(3) | |
| 65 |
2/2✓ Branch 0 taken 45872 times.
✓ Branch 1 taken 7676922 times.
|
7722794 | elseif(bmod.lt.bmin) then |
| 66 | 45872 | bmin=bmod | |
| 67 | 45872 | tmin=x(2) | |
| 68 | 45872 | pmin=x(3) | |
| 69 | endif | ||
| 70 | ! | ||
| 71 | enddo | ||
| 72 | enddo | ||
| 73 | ! | ||
| 74 | 808 | call newtextr(tmin,pmin,bmin) | |
| 75 | ! | ||
| 76 | 808 | call newtextr(tmax,pmax,bmax) | |
| 77 | ! | ||
| 78 | !----------- | ||
| 79 | contains | ||
| 80 | !----------- | ||
| 81 | ! | ||
| 82 | 8161608 | subroutine getbmod | |
| 83 | use magfie_sub, only : magfie | ||
| 84 | ! | ||
| 85 | implicit none | ||
| 86 | ! | ||
| 87 | real(dp) :: sqrtg | ||
| 88 | real(dp), dimension(3) :: hcovar, hctrvr, hcurl | ||
| 89 | ! | ||
| 90 | 8161608 | call magfie(x,bmod,sqrtg,bder,hcovar,hctrvr,hcurl) | |
| 91 | ! | ||
| 92 | 8161608 | end subroutine getbmod | |
| 93 | ! | ||
| 94 | !----------- | ||
| 95 | ! | ||
| 96 | 1616 | subroutine newtextr(theta,phi,bextr) | |
| 97 | ! | ||
| 98 | implicit none | ||
| 99 | ! | ||
| 100 | integer, parameter :: niter=10 | ||
| 101 | real(dp), parameter :: hdt=1.d-3, hdp=1.d-3 | ||
| 102 | real(dp) :: bt,bp,btt,btp,bpt,bpp,theta,phi,bextr | ||
| 103 | real(dp) :: det,dtheta,dphi,w | ||
| 104 | ! | ||
| 105 |
2/2✓ Branch 0 taken 16160 times.
✓ Branch 1 taken 1616 times.
|
17776 | do iter=1,niter |
| 106 | 16160 | x(3)=phi | |
| 107 | 16160 | x(2)=theta | |
| 108 | ! | ||
| 109 | 16160 | call getbmod | |
| 110 | ! | ||
| 111 | 16160 | bt=bder(2)*bmod | |
| 112 | 16160 | bp=bder(3)*bmod | |
| 113 | ! | ||
| 114 | 16160 | x(3)=phi | |
| 115 | 16160 | x(2)=theta+hdt | |
| 116 | ! | ||
| 117 | 16160 | call getbmod | |
| 118 | ! | ||
| 119 | 16160 | btt=bder(2)*bmod | |
| 120 | 16160 | bpt=bder(3)*bmod | |
| 121 | 16160 | x(2)=theta-hdt | |
| 122 | ! | ||
| 123 | 16160 | call getbmod | |
| 124 | ! | ||
| 125 | 16160 | btt=(btt-bder(2)*bmod)/(2.d0*hdt) | |
| 126 | 16160 | bpt=(bpt-bder(3)*bmod)/(2.d0*hdt) | |
| 127 | 16160 | btp=bpt | |
| 128 | 16160 | x(2)=theta | |
| 129 | 16160 | x(3)=phi+hdp | |
| 130 | ! | ||
| 131 | 16160 | call getbmod | |
| 132 | ! | ||
| 133 | 16160 | bpp=bder(3)*bmod | |
| 134 | ! | ||
| 135 | 16160 | x(3)=phi-hdp | |
| 136 | ! | ||
| 137 | 16160 | call getbmod | |
| 138 | ! | ||
| 139 | 16160 | bpp=(bpp-bder(3)*bmod)/(2.d0*hdp) | |
| 140 | 16160 | det=btt*bpp-btp*bpt | |
| 141 | 16160 | dtheta=(bt*bpp-bp*btp)/det | |
| 142 | 16160 | dphi=(bp*btt-bt*bpt)/det | |
| 143 | 16160 | w=min(ht/max(abs(dtheta),ht),hp/max(abs(dphi),hp)) | |
| 144 | 16160 | theta=theta-w*dtheta | |
| 145 | 17776 | phi=phi-w*dphi | |
| 146 | enddo | ||
| 147 | ! | ||
| 148 | 1616 | bextr=bmod | |
| 149 | ! | ||
| 150 | 8161608 | end subroutine newtextr | |
| 151 | ! | ||
| 152 | !----------- | ||
| 153 | ! | ||
| 154 | end subroutine find_bminmax | ||
| 155 | ! | ||
| 156 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc | ||
| 157 | ! | ||
| 158 | 8 | subroutine init_bminmax_arrays | |
| 159 | use bminmax_mod, only: prop, nsbmnx, hsbmnx, bmin_arr, bmax_arr | ||
| 160 | |||
| 161 | integer :: k | ||
| 162 | real(dp) :: s0 | ||
| 163 | |||
| 164 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
|
8 | if (.not. prop) return |
| 165 | 8 | prop = .false. | |
| 166 | 8 | hsbmnx = 1.d0/dble(nsbmnx) | |
| 167 | ! The cache is defined in the coordinate system of the active magfie | ||
| 168 | ! backend. The driver initializes it serially after selecting the tracing | ||
| 169 | ! backend, before parallel orbit tracing reads it. | ||
| 170 |
2/2✓ Branch 0 taken 808 times.
✓ Branch 1 taken 8 times.
|
816 | do k = 0, nsbmnx |
| 171 | 808 | s0 = max(1.d-8, hsbmnx*dble(k)) | |
| 172 | 816 | call find_bminmax(s0, bmin_arr(k), bmax_arr(k)) | |
| 173 | end do | ||
| 174 | end subroutine init_bminmax_arrays | ||
| 175 | |||
| 176 | 45 | subroutine get_bminmax(s, bmin, bmax) | |
| 177 | use bminmax_mod, only: prop, nsbmnx, hsbmnx, bmin_arr, bmax_arr | ||
| 178 | |||
| 179 | real(dp), intent(in) :: s | ||
| 180 | real(dp), intent(out) :: bmin, bmax | ||
| 181 | integer :: k | ||
| 182 | real(dp) :: ws | ||
| 183 | |||
| 184 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
|
45 | if (prop) call init_bminmax_arrays |
| 185 | |||
| 186 | 45 | ws = s/hsbmnx | |
| 187 | 45 | k = min(nsbmnx - 1, max(0, int(ws))) | |
| 188 | 45 | ws = ws - dble(k) | |
| 189 | 45 | bmin = bmin_arr(k)*(1.d0 - ws) + bmin_arr(k + 1)*ws | |
| 190 | 45 | bmax = bmax_arr(k)*(1.d0 - ws) + bmax_arr(k + 1)*ws | |
| 191 | 45 | end subroutine get_bminmax | |
| 192 | |||
| 193 | end module find_bminmax_sub | ||
| 194 |