V3FIT
vertex_list.f
1  MODULE vertex_list
2  USE stel_kinds
3  IMPLICIT NONE
4  INTEGER, PARAMETER :: nvertex = 4
5 
6  PRIVATE :: find_vertex
7 
8  TYPE point
9  INTEGER :: xpt, ypt, phi_plane
10  END TYPE point
11 
12  TYPE magfield
13  REAL(rprec) :: br, bp, bz
14  END TYPE magfield
15 
16 
17  TYPE pbfield
18  TYPE (MAGFIELD), POINTER :: pbfld !POINTER to BFIELD
19  END TYPE pbfield
20 
21 
22  TYPE vertex
23  TYPE (POINT) :: p1 !integer-INDEXed position of vertex (lower-left point)
24  TYPE (VERTEX), POINTER :: pprev !POINTER to previously ALLOCATED vertex (next in list)
25  TYPE (PBFIELD) :: bvector(nvertex) !pointers to the bfield vector at nvertex nodes of square
26  LOGICAL :: lalloc(nvertex) !true IF bvector(k) is ALLOCATED
27  END TYPE vertex
28 
29  TYPE pvertex !POINTER to VERTEX
30  TYPE (VERTEX), POINTER :: pverx
31  END TYPE pvertex
32 
33  INTEGER :: nlist !number of lists
34  REAL(rprec) :: hR, hZ, hP !scale factors for converting xpt to R, etc.
35  INTEGER, ALLOCATABLE :: ibs_cnt(:) !counts number of biot-savart calls made
36  TYPE (PVERTEX), ALLOCATABLE, TARGET :: tail(:) !array of root vertices for separate lists
37  TYPE (VERTEX), POINTER :: root !last vertex most recently added to list
38  TYPE (VERTEX), POINTER :: new_root, prev_vertex !stores previously computed nodes
39  LOGICAL :: lsearch(nvertex) !true IF bfield at ith node already computed
40 
41  CONTAINS
42 
43  SUBROUTINE initialize_list (number_lists, hR_in, hZ_in, hP_in)
44  IMPLICIT NONE
45  INTEGER, INTENT(in) :: number_lists
46  REAL(rprec), INTENT(in) :: hR_in, hZ_in, hP_in
47  INTEGER :: ilist
48 
49  nlist = number_lists
50  IF (nlist .le. 0) stop ' Must specify nlist > 0'
51 
52  ALLOCATE (tail(nlist), ibs_cnt(nlist))
53 
54  DO ilist = 1, nlist
55  NULLIFY (tail(ilist)%pverx)
56  END DO
57 
58  ibs_cnt(1:nlist) = 0
59 
60  hr = hr_in; hz = hz_in; hp = hp_in
61 
62  END SUBROUTINE initialize_list
63 
64 
65  SUBROUTINE add_vertex (rin, zin, phi_pt, list_index, bfield)
66  IMPLICIT NONE
67  TYPE (VERTEX), POINTER :: new_vertex, swap
68  TYPE (MAGFIELD), POINTER :: MagFld
69  INTEGER, INTENT(in) :: phi_pt, list_index
70  REAL(rprec), INTENT(in) :: rin, zin
71  REAL(rprec) :: Rreal, Zreal, PhiReal
72  INTEGER :: rpt, zpt, xpt, ypt
73  INTEGER :: icount
74  LOGICAL :: lexists
75  EXTERNAL bfield !Could be magnetic field or something ELSE...
76 !
77 ! adds new_vertex to tail of list (if it does not exist already), assigning it to root
78 ! If it exists already, moves vertex to tail of list, assigning it to root
79 ! thus, this vertex can be identified with "root" for use by external programs
80 
81  IF (list_index .gt. nlist) stop 'list_index > nlist in add_vertex'
82  IF (list_index .le. 0) stop 'list_index <= 0 in add_vertex'
83 
84 ! Get lower left INDEX of nearest vertex. Z can be positive or negative so USE FLOOR
85  rpt = floor(rin/hr)
86  zpt = floor(zin/hz)
87 
88  root => tail(list_index)%pverx
89 
90  ALLOCATE (new_vertex)
91  IF (ASSOCIATED(root)) THEN
92  new_vertex%pprev => root !point to previous node unless this is the first in list
93  ELSE
94  NULLIFY(new_vertex%pprev)
95  END IF
96 
97  new_vertex%p1 = point(rpt, zpt, phi_pt) !store INTEGER INDEX of lower-left vertex on grid
98 
99 ! Assign bfield values to this vertex from exist vertex values
100  CALL find_vertex (new_vertex, lexists)
101 
102  IF (lexists) THEN
103 ! new_node already existed!
104 ! move it to root position. and assign root to new_root
105 ! prev_vertex is the preceding vertex in list
106  DEALLOCATE (new_vertex)
107  IF (ASSOCIATED(prev_vertex)) THEN
108  swap => root
109  root => new_root
110  root%pprev => swap
111  END IF
112 ! PRINT *,' NODE ALREADY EXISTED IN LIST #', list_index
113 
114  ELSE
115  root => new_vertex
116 
117  DO icount = 1, nvertex
118  IF (.not.lsearch(icount)) THEN
119 ! Node is one of four corners of square, depending on value of icount (1-4)
120  ALLOCATE (new_vertex%bvector(icount)%pbfld)
121  new_vertex%lalloc(icount) = .true.
122  magfld => new_vertex%bvector(icount)%pbfld
123  !Simulate bio-savart CALL; must add toroidal plane in arg list ...
124  xpt = rpt; ypt = zpt
125  IF (icount.eq.2 .or. icount.eq.4) xpt = xpt+1
126  IF (icount.ge.3) ypt = ypt+1
127  rreal = xpt*hr; zreal = ypt*hz
128  phireal = (phi_pt-1)*hp
129  CALL bfield (rreal, phireal, zreal,
130  1 magfld%br, magfld%bp, magfld%bz)
131  ibs_cnt(list_index) = ibs_cnt(list_index) + 1
132  ELSE
133  new_vertex%lalloc(icount) = .false.
134  END IF
135  END DO
136  END IF
137 
138  tail(list_index)%pverx => root
139 
140  END SUBROUTINE add_vertex
141 
142 
143  SUBROUTINE find_vertex (new_node, lexists)
144  IMPLICIT NONE
145  TYPE (VERTEX), POINTER :: new_node
146  LOGICAL, INTENT(out) :: lexists
147  TYPE (VERTEX), POINTER :: stepit, previous
148  INTEGER :: xpt(nvertex), ypt(nvertex), xvert,
149  1 yvert, xstep, ystep, icount, ivertex
150 
151 ! Scans vertex (assumed in a given phi plane) ASSOCIATED with new_node and assigns bfield values
152 ! IF they already exist in list. If it finds the bfield ASSOCIATED with the new_node
153 ! or ANY of the other 3 corners, the global pointers new_root, prev_vertex become ASSOCIATED.
154 
155 ! Find nodes associate with the vertex new_node WHERE Bfield is required
156  xpt(1) = new_node%p1%xpt; ypt(1) = new_node%p1%ypt
157  xpt(2) = xpt(1)+1; ypt(2) = ypt(1)
158  xpt(3) = xpt(1); ypt(3) = ypt(1)+1
159  xpt(4) = xpt(2); ypt(4) = ypt(3)
160 
161  stepit => root;
162  NULLIFY(previous); NULLIFY(new_root); NULLIFY(prev_vertex)
163  lexists = .false.
164 
165  lsearch(1:nvertex) = .false.
166 
167  icount = 0
168  DO WHILE (icount.lt.nvertex .and. ASSOCIATED(stepit))
169  xstep = stepit%p1%xpt; ystep = stepit%p1%ypt
170 !
171 ! Logic to catch case whether vertex already exists
172 ! This may miss the case if all nodes are filled first by adjacent vertices, but does not matter
173 !
174  IF (xstep.eq.xpt(1) .and. ystep.eq.ypt(1)) THEN
175  new_root => stepit
176  IF (ASSOCIATED(previous)) previous%pprev => stepit%pprev
177  prev_vertex => previous
178  lexists = .true.
179  EXIT
180  END IF
181  DO ivertex = 1, nvertex
182  xvert = xpt(ivertex); yvert = ypt(ivertex)
183 
184  IF ((xstep.gt.xvert) .or. (ystep.gt.yvert) .or.
185  1 (xstep.lt.xvert-1) .or. (ystep.lt.yvert-1) .or.
186  2 lsearch(ivertex) ) cycle
187 
188  SELECT CASE (ivertex)
189  CASE (1) !!Lower left vertex
190  IF (xstep.eq.xvert) THEN
191  new_node%bvector(1) = stepit%bvector(3)
192  ELSE
193  IF (ystep .eq. yvert) THEN
194  new_node%bvector(1) = stepit%bvector(2)
195  ELSE !!ystep < yvert
196  new_node%bvector(1) = stepit%bvector(4)
197  END IF
198  END IF
199 
200  CASE (2) !!Lower right vertex
201  IF (xstep.lt.xvert) THEN
202  new_node%bvector(2) = stepit%bvector(4)
203  ELSE !!xstep == xvert
204  IF (ystep .eq. yvert) THEN
205  new_node%bvector(2) = stepit%bvector(1)
206  ELSE !!ystep < yvert
207  new_node%bvector(2) = stepit%bvector(3)
208  END IF
209  END IF
210 
211  CASE (3) !!Upper left vertex
212  IF (xstep .eq. xvert) THEN
213  new_node%bvector(3) = stepit%bvector(1)
214  ELSE !!xstep < xvert
215  IF (ystep .eq. yvert) THEN
216  new_node%bvector(3) = stepit%bvector(2)
217  ELSE !!ystep < yvert
218  new_node%bvector(3) = stepit%bvector(4)
219  END IF
220  END IF
221 
222  CASE (4) !!Upper right vertex
223  IF (xstep .lt. xvert) THEN
224  new_node%bvector(4) = stepit%bvector(2)
225  ELSE !!xstep == xvert
226  IF (ystep .eq. yvert) THEN
227  new_node%bvector(4) = stepit%bvector(1)
228  ELSE !!ystep < yvert
229  new_node%bvector(4) = stepit%bvector(3)
230  END IF
231  END IF
232 
233  END SELECT
234 
235  icount = icount + 1
236  IF (icount .ge. nvertex) EXIT
237  lsearch(ivertex) = .true.
238 
239  END DO
240 
241  previous => stepit
242  stepit => stepit%pprev
243 
244  END DO
245 
246  END SUBROUTINE find_vertex
247 
248 
249  SUBROUTINE b_at_point (rin, zin, br, bp, bz)
250  IMPLICIT NONE
251 C-----------------------------------------------
252 C D u m m y A r g u m e n t s
253 C-----------------------------------------------
254  REAL(rprec), INTENT(in) :: rin, zin
255  REAL(rprec), INTENT(out) :: br, bp, bz
256  REAL(rprec), DIMENSION(nvertex) :: brvac, bpvac, bzvac
257 C-----------------------------------------------
258 C L o c a l V a r i a b l e s
259 C-----------------------------------------------
260  TYPE (MAGFIELD), POINTER :: bvac
261  REAL(rprec), PARAMETER :: zero = 0, one = 1
262  INTEGER :: ir, jz, i
263  REAL(rprec) :: ri, zj, dri, dzj, fint(nvertex)
264 C-----------------------------------------------
265 !
266 ! THIS ROUTINE FINDS THE CYLINDRICAL COMPONENTS OF THE EXTERNAL
267 ! MAGNETIC FIELD AT A FIXED PHI PLANE AT THE POINT (rin, zin) BY 2-D INTERPOLATION
268 !
269 ! THIS CALL MUST BE MADE JUST AFTER "ADD_VERTEX" CALL, SO THAT
270 ! "ROOT" CONTAINS THE CORRECT LOWER-LEFT VERTEX IN THE PROPER PHI-PLANE
271 !
272 
273 !
274 ! DETERMINE R, Z VERTEX CORNER INDICES (IR,IZ)
275 !
276  ir = root%p1%xpt
277  jz = root%p1%ypt
278 
279 !
280 ! COMPUTE R , Z AND DR , DZ WITH RESPECT TO LOWER-LEFT VERTEX
281 !
282 
283  ri = ir*hr
284  zj = jz*hz
285 
286  dri = (rin - ri)/hr
287  dzj = (zin - zj)/hz
288 
289 !
290 ! RETRIEVE B FIELD AT 4 VERTICES OF ROOT GRID ELEMENT
291 !
292 
293  DO i = 1, nvertex
294  bvac => root%bvector(i)%pbfld
295  brvac(i) = bvac%br
296  bpvac(i) = bvac%bp
297  bzvac(i) = bvac%bz
298  END DO
299 
300 !
301 ! COMPUTE INTERPOLATED B FIELD
302 !
303  fint(4) = dri*dzj
304  fint(3) = dzj - fint(4)
305  fint(2) = dri - fint(4)
306  fint(1) = 1 + fint(4) - (dri + dzj)
307 
308  IF (dri.lt.zero .or. dri.gt.one .or.
309  1 dzj.lt.zero .or. dzj.gt.one) THEN
310  WRITE (6, '(2(a,1pe12.4))')
311  1 'Possible interp. error in b_from_coil: dri = ',
312  2 dri,' dzj = ', dzj
313  END IF
314 
315  br = sum(fint*brvac)
316  bp = sum(fint*bpvac)
317  bz = sum(fint*bzvac)
318 
319  END SUBROUTINE b_at_point
320 
321  SUBROUTINE cleanup_list
322  IMPLICIT NONE
323  TYPE (VERTEX), POINTER :: stepit
324  INTEGER :: icount, ilist
325 
326  DO ilist = 1, nlist
327  root => tail(ilist)%pverx
328  DO WHILE (ASSOCIATED(root))
329  stepit => root%pprev
330  DO icount = 1, nvertex
331  IF (root%lalloc(icount))
332  1 DEALLOCATE (root%bvector(icount)%pbfld)
333  END DO
334  DEALLOCATE (root)
335  root => stepit
336  END DO
337  END DO
338 
339  DEALLOCATE (tail, ibs_cnt)
340 
341  END SUBROUTINE cleanup_list
342 
343  END MODULE vertex_list
vertex_list::magfield
Definition: vertex_list.f:12
vertex_list::pbfield
Definition: vertex_list.f:17
vertex_list::pvertex
Definition: vertex_list.f:29
vertex
A vertex.
Definition: vertex.hpp:21
vertex_list::point
Definition: vertex_list.f:8