V3FIT
grid_extension.f90
Go to the documentation of this file.
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
9 !*******************************************************************************
10  MODULE grid_extension
11  USE vmec_info, ONLY: jcurrumnc, jcurrvmnc
12  USE stel_kinds
13  USE stel_constants
14  USE shared_data, ONLY: lasym, r1_i, z1_i, ru_i, zu_i, rv_i, zv_i, nsp
15  USE island_params, ns=>ns_i, mpol=>mpol_i, ntor=>ntor_i, &
16  ntheta=>nu_i, nzeta=>nv_i, mnmax=>mnmax_i, nsv=>ns_v, &
17  ohs=>ohs_i, nfp=>nfp_i
18  USE island_params, ONLY: phipf_i, chipf_i, presf_i, nsh
19  USE v3_utilities, ONLY: assert
20  USE siesta_namelist, ONLY: l_vessel
21  USE descriptor_mod, ONLY: iam
22 
23  IMPLICIT NONE
24 
25 !*******************************************************************************
26 ! grid extension module parameters.
27 !*******************************************************************************
29  INTEGER, PARAMETER :: lin = 1
31  INTEGER, PARAMETER :: quad = 2
33  INTEGER, PARAMETER :: cub = 3
35  REAL (dp), PARAMETER :: p5 = 0.5_dp
36 
37 !*******************************************************************************
38 ! grid extension module variables.
39 !*******************************************************************************
41  INTEGER :: interp_type
42  INTEGER :: nse
43 
44  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: ri, zi
45 
46  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1_vmec, z1_vmec
47  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: ru_vmec, zu_vmec
48  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: rv_vmec, zv_vmec
49 
50  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1s_vmec, r1ss_vmec
51  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: z1s_vmec, z1ss_vmec
52 
53  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1su_vmec, z1su_vmec
54  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1sv_vmec, z1sv_vmec
55 
56  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: rmnc_v, zmns_v
57  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: rmns_v, zmnc_v
58 
59  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: r_ext, z_ext
60  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: ru_ext, zu_ext
61  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: rv_ext, zv_ext
62 
63  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1_ves, z1_ves
64  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: ru_ves, zu_ves
65  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: rv_ves, zv_ves
66 
67  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: r_full, z_full
68  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: ru_full, zu_full
69  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: rv_full, zv_full
70 
71  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: sqrtg_full
72 
73  REAL(dp) :: s_edge, ming, maxg
74 
75  CONTAINS
76 
77 !*******************************************************************************
78 ! UTILITY SUBROUTINES
79 !*******************************************************************************
80  SUBROUTINE grid_extender(wout_file, itype, istat)
81  USE fourier, ONLY: f_cos, f_sin, f_none, f_sum, f_du, f_dv
82  USE v3_utilities, ONLY: assert
83  USE vmec_info, ONLY: rmnc_i, rmns_i, zmns_i, zmnc_i, &
84  jcurrumnc, jcurrvmnc, jcurrumns, jcurrvmns
85  USE descriptor_mod, ONLY: siesta_comm
86  USE siesta_namelist, ONLY: nsin_ext
87  USE shared_data, ONLY: lverbose
88  USE island_params, ONLY: fourier_context
89 
90  IMPLICIT NONE
91 
92 ! Declare Arguments
93  CHARACTER(*), INTENT(in) :: wout_file
94  CHARACTER(*), INTENT(in) :: itype
95  INTEGER, INTENT(out) :: istat
96 
97  CHARACTER(LEN=10) :: chfit
98  CHARACTER(LEN=100) :: device
99 
100  INTEGER :: j, m, n, u, v, i, k, js, js1, nsp1
101  REAL(dp) :: rho, vol0, vol1, vols, rho1, rho2, rho3, vp1, s_ext, delv
102  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: r12, rs12, zs12, ru12, zu12
103  REAL(dp), DIMENSION(:), ALLOCATABLE :: vp
104  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: tmp, tmp2
105  REAL(dp), DIMENSION(:), ALLOCATABLE :: tmp1
106 
107  class(fourier_class), POINTER :: four => null()
108 
109  IF (.NOT.l_vessel) RETURN ! FIXME: This check should be outsize the subroutine
110 
111  CALL init_extender
112 
113  i = index(wout_file, '.', back=.true.)
114  j = index(wout_file,'_')
115  IF (j .LT. i) THEN
116  device = wout_file(j+1:i-1)
117  ELSE
118  device = wout_file(1:i)
119  END IF
120 
121 ! chfit = 'quad'
122 ! chfit = TRIM(itype)
123  IF (index(itype,'lin').EQ.1 .OR. index(itype,'LIN').EQ.1) THEN
124  interp_type = lin; chfit = "linear"
125  ELSE IF (index(itype,'quad').EQ.1 .OR. index(itype,'QUAD').EQ.1) THEN
126  interp_type = quad; chfit = "quadratic"
127  ELSE IF (index(itype,'cub').EQ.1 .OR. index(itype,'CUB').EQ.1) THEN
128  interp_type = cub; chfit = "cubic"
129  ELSE
130  CALL assert(.false., 'Unknown interpolation <TYPE> specified')
131  END IF
132 
133  ALLOCATE(ri(ntheta,nzeta), zi(ntheta,nzeta))
134  ALLOCATE(r1_ves(ntheta,nzeta), z1_ves(ntheta,nzeta), &
135  ru_ves(ntheta,nzeta), zu_ves(ntheta,nzeta), &
136  rv_ves(ntheta,nzeta), zv_ves(ntheta,nzeta))
137 
138 ! FIXME: The vessel file could have a different periodicity than the equilbirum.
139  four => fourier_class(mpol_v, ntor_v, ntheta, nzeta, nfp, lasym)
140 
141 ! Set internal normalization
142  rmnc_v = rmnc_v/four%orthonorm
143  zmns_v = zmns_v/four%orthonorm
144 
145 ! FIXME: What about non stellarator symmetric vessels?
146  CALL four%toijsp(rmnc_v, r1_ves, f_none, f_cos)
147  CALL four%toijsp(zmns_v, z1_ves, f_none, f_sin)
148 
149  CALL four%toijsp(rmnc_v, ru_ves, f_du, f_cos)
150  CALL four%toijsp(zmns_v, zu_ves, f_du, f_sin)
151 
152  CALL four%toijsp(rmnc_v, rv_ves, f_dv, f_cos)
153  CALL four%toijsp(zmns_v, zv_ves, f_dv, f_sin)
154 
155  DEALLOCATE(four)
156 
157 ! Compute extended "s" coordinate for s>1, assuming Vp(s) = Vp(s=1)*s
158 ! where Vp(s=1) is obtained from extrapolating VMEC values. Results is
159 ! V(s) = Vp*(s^2-1)/2 + V0 => s_edge = SQRT(2(V1-V0)/Vp + 1)
160  vol0 = get_volume(zero) !rho = 0: s=1 vmec
161  vol1 = get_volume(one) !rho = 1: vessel
162 
163  vp1 = getvp1()
164  s_edge = sqrt(1 + 2*(vol1 - vol0)/vp1)
165 
166 ! NOTE: at js=nse, s(js) == (js-1)*hs = s_edge
167  nse = nint(s_edge*ohs) + 1
168 
169  nsv = nse - ns
170  s_edge = (nsv - 1)*hs_i + 1
171 
172 ! Rounding discrepancy in volume due to integer nse
173  delv = vol1 - (vol0 + p5*vp1*(s_edge*s_edge - 1))
174  delv = delv/(s_edge-1)**2
175 
176  ALLOCATE(r_ext(ntheta,nzeta,2:nsv), z_ext(ntheta,nzeta,2:nsv), &
177  ru_ext(ntheta,nzeta,2:nsv), zu_ext(ntheta,nzeta,2:nsv), &
178  rv_ext(ntheta,nzeta,2:nsv), zv_ext(ntheta,nzeta,2:nsv), &
179  stat=istat)
180 
181  DO j = 2, nsv
182 
183  s_ext = 1 + (j-1)*hs_i !at j=nsv, s_ext = s_edge
184 
185  !V(s) = Vp*(s^2-1)/2 + V0
186  vols = vol0 + p5*vp1*(s_ext**2-1) !matches vp exactly at s=1
187  vols = vols + delv*(s_ext-1)**2 !matches vol1 exactly at s=s_edge
188  vols = min(vols, vol1)
189  CALL remap_radial(rho, vols)
190 
191  rho1 = 1 - rho
192  IF (interp_type.EQ.lin) THEN
193  r_ext(:,:,j) = rho1*r1_vmec + rho*r1_ves
194  z_ext(:,:,j) = rho1*z1_vmec + rho*z1_ves
195  ru_ext(:,:,j) = rho1*ru_vmec + rho*ru_ves
196  zu_ext(:,:,j) = rho1*zu_vmec + rho*zu_ves
197  rv_ext(:,:,j) = rho1*rv_vmec + rho*rv_ves
198  zv_ext(:,:,j) = rho1*zv_vmec + rho*zv_ves
199  ELSE
200  rho2 = rho*rho
201  IF (interp_type .eq. quad) THEN
202 ! QUADRATIC (R' match at rho = 0)
203  r_ext(:,:,j) = p5*r1s_vmec*(1 - rho2 - rho1*rho1) &
204  + r1_vmec*(1-rho2) + r1_ves*rho2
205  z_ext(:,:,j) = p5*z1s_vmec*(1 - rho2 - rho1*rho1) &
206  + z1_vmec*(1-rho2) + z1_ves*rho2
207  ru_ext(:,:,j) = p5*r1su_vmec*(1 - rho2 - rho1*rho1) &
208  + ru_vmec*(1-rho2) + ru_ves*rho2
209  zu_ext(:,:,j) = p5*z1su_vmec*(1 - rho2 - rho1*rho1) &
210  + zu_vmec*(1-rho2) + zu_ves*rho2
211  rv_ext(:,:,j) = p5*r1sv_vmec*(1 - rho2 - rho1*rho1) &
212  + rv_vmec*(1-rho2) + rv_ves*rho2
213  zv_ext(:,:,j) = p5*z1sv_vmec*(1 - rho2 - rho1*rho1) &
214  + zv_vmec*(1-rho2) + zv_ves*rho2
215  ELSE IF (interp_type .eq. cub) THEN
216  rho3 = rho*rho2
217  CALL assert(.false., 'Unknown interpolation <TYPE> specified')
218  END IF
219  END IF
220  END DO
221 
222  nse = nse-1 !Not to double count s=1 surface
223  ALLOCATE(r_full(ntheta,nzeta,nse), z_full(ntheta,nzeta,nse), &
224  ru_full(ntheta,nzeta,nse), zu_full(ntheta,nzeta,nse), &
225  rv_full(ntheta,nzeta,nse), zv_full(ntheta,nzeta,nse))
226 
227  r_full(:,:,1:ns) = r1_i
228  r_full(:,:,ns + 1:nse) = r_ext
229  z_full(:,:,1:ns) = z1_i
230  z_full(:,:,ns + 1:nse) = z_ext
231  ru_full(:,:,1:ns) = ru_i
232  ru_full(:,:,ns + 1:nse) = ru_ext
233  zu_full(:,:,1:ns) = zu_i
234  zu_full(:,:,ns + 1:nse) = zu_ext
235  rv_full(:,:,1:ns) = rv_i
236  rv_full(:,:,ns + 1:nse) = rv_ext
237  zv_full(:,:,1:ns) = zv_i
238  zv_full(:,:,ns + 1:nse) = zv_ext
239 
240  nse = min(ns + nsin_ext, nse)
241 
242  IF (lverbose .and. iam .eq. 0) THEN
243  WRITE (*,1000) nse, ns, nsin_ext
244  END IF
245 1000 FORMAT('Using ',i4,' Total Surfaces.',/, &
246  'Number of VMEC surfaces ',i4,/, &
247  'Number of extended surfaces ',i4)
248 
249 ! Assign the new value of ns.
250 ! REASSIGN THE NEW VALUE OF NS
251  nsp = ns
252  ns = nse
253  nsh = ns - 1
254 
255  IF (nsp .NE. ns) THEN
256  DEALLOCATE(r1_i)
257  DEALLOCATE(z1_i)
258  DEALLOCATE(ru_i)
259  DEALLOCATE(zu_i)
260  DEALLOCATE(rv_i)
261  DEALLOCATE(zv_i)
262  ALLOCATE(r1_i(ntheta,nzeta,ns), z1_i(ntheta,nzeta,ns), &
263  ru_i(ntheta,nzeta,ns), zu_i(ntheta,nzeta,ns), &
264  rv_i(ntheta,nzeta,ns), zv_i(ntheta,nzeta,ns))
265  r1_i = r_full(:,:,1:ns)
266  z1_i = z_full(:,:,1:ns)
267  ru_i = ru_full(:,:,1:ns)
268  zu_i = zu_full(:,:,1:ns)
269  rv_i = rv_full(:,:,1:ns)
270  zv_i = zv_full(:,:,1:ns)
271  END IF
272 
273  DEALLOCATE(r_full, z_full)
274  DEALLOCATE(ru_full, zu_full)
275  DEALLOCATE(rv_full, zv_full)
276 
277  IF (nsp .NE. ns) THEN
278  nsp1 = nsp + 1
279  ALLOCATE(tmp(0:mpol,-ntor:ntor,ns), tmp2(0:mpol,-ntor:ntor, ns))
280  tmp(:,:,1:nsp) = rmnc_i
281  tmp2(:,:,1:nsp) = zmns_i
282  DEALLOCATE(rmnc_i, zmns_i)
283  ALLOCATE(rmnc_i(0:mpol,-ntor:ntor,ns), &
284  zmns_i(0:mpol,-ntor:ntor,ns), stat=istat)
285  rmnc_i(:,:,1:nsp) = tmp(:,:,1:nsp)
286  zmns_i(:,:,1:nsp) = tmp2(:,:,1:nsp)
287  CALL fourier_context%tomnsp(r1_i(:,:,nsp1:ns), rmnc_i(:,:,nsp1:ns), &
288  f_cos)
289  CALL fourier_context%tomnsp(z1_i(:,:,nsp1:ns), zmns_i(:,:,nsp1:ns), &
290  f_sin)
291  DO js = nsp1, ns
292  rmnc_i(:,:,js) = rmnc_i(:,:,js) &
293  / fourier_context%orthonorm(0:mpol,-ntor:ntor)
294  zmns_i(:,:,js) = zmns_i(:,:,js) &
295  / fourier_context%orthonorm(0:mpol,-ntor:ntor)
296  END DO
297 
298 ! PUT PROFILES ONTO EXTENDED MESH
299  tmp = 0
300  tmp(:,:,1:nsp) = jcurrumnc
301  DEALLOCATE(jcurrumnc)
302  ALLOCATE(jcurrumnc(0:mpol,-ntor:ntor,ns))
303  jcurrumnc = tmp
304 
305  tmp = 0
306  tmp(:,:,1:nsp) = jcurrvmnc
307  DEALLOCATE(jcurrvmnc)
308  ALLOCATE(jcurrvmnc(0:mpol,-ntor:ntor,ns))
309  jcurrvmnc = tmp
310 
311 ! DO THE SAME FOR THE phipf_i, chipf_i AND presf_i
312 ! phipf, chipf will be computed in PCHELMS in the vacuum region
313  ALLOCATE(tmp1(ns))
314  tmp1 = 0
315  tmp1(1:nsp) = phipf_i
316  DEALLOCATE(phipf_i)
317  ALLOCATE(phipf_i(ns))
318  phipf_i = tmp1
319 
320  tmp1 = 0
321  tmp1(1:nsp) = chipf_i
322  DEALLOCATE(chipf_i)
323  ALLOCATE(chipf_i(ns))
324  chipf_i = tmp1
325 
326  tmp1(1:nsp) = presf_i
327  tmp1(nsp + 1:) = presf_i(nsp)
328  DEALLOCATE(presf_i)
329  ALLOCATE(presf_i(ns))
330  presf_i = tmp1
331  DEALLOCATE(tmp1)
332 
333  IF (lasym) THEN
334  tmp(:,:,1:nsp) = rmns_i
335  tmp2(:,:,1:nsp) = zmnc_i
336  DEALLOCATE(rmns_i, zmnc_i)
337  ALLOCATE(rmns_i(0:mpol,-ntor:ntor,ns), &
338  zmnc_i(0:mpol,-ntor:ntor,ns), stat=istat)
339  rmns_i(:,:,1:nsp) = tmp(:,:,1:nsp)
340  zmnc_i(:,:,1:nsp) = tmp2(:,:,1:nsp)
341  CALL fourier_context%tomnsp(r1_i(:,:,nsp1:ns), rmns_i(:,:,nsp1:ns), &
342  f_sin)
343  CALL fourier_context%tomnsp(z1_i(:,:,nsp1:ns), zmnc_i(:,:,nsp1:ns), &
344  f_cos)
345  DO js = nsp1, ns
346  rmns_i(:,:,js) = rmns_i(:,:,js) &
347  / fourier_context%orthonorm(0:mpol,-ntor:ntor)
348  zmnc_i(:,:,js) = zmnc_i(:,:,js) &
349  / fourier_context%orthonorm(0:mpol,-ntor:ntor)
350  END DO
351 
352  tmp = 0
353  tmp(:,:,1:nsp)=jcurrumns
354  DEALLOCATE(jcurrumns)
355  ALLOCATE(jcurrumns(0:mpol,-ntor:ntor,ns))
356  jcurrumns=tmp
357 
358  tmp=0
359  tmp(:,:,1:nsp)=jcurrvmns
360  DEALLOCATE(jcurrvmns)
361  ALLOCATE(jcurrvmns(0:mpol,-ntor:ntor,ns))
362  jcurrvmns=tmp
363  END IF
364 
365  DEALLOCATE(tmp, tmp2)
366 
367  END IF
368 
369 #if defined(NOSKS)
370  DO u=1, ntheta
371  DO v=1, nzeta
372  DO js=1, ns
373  WRITE(1600,*) r1_i(u,v,js), z1_i(u,v,js)
374  CALL flush(1600)
375  END DO
376  END DO
377  END DO
378 
379  ! u = theta, v = zeta
380  ALLOCATE (r12(ntheta,nzeta), &
381  rs12(ntheta,nzeta), &
382  zs12(ntheta,nzeta), &
383  ru12(ntheta,nzeta), &
384  zu12(ntheta,nzeta), &
385  sqrtg_full(ntheta,nzeta,nse), &
386  stat = istat)
387  CALL assert(istat.eq.0,' Allocation error in grid_extender')
388  sqrtg_full(:,:,1) = 0
389  DO js = 2, nse
390  js1 = js-1
391  r12 = p5*(r1_i(:,:,js) + r1_i(:,:,js1))
392  rs12 = (r1_i(:,:,js) - r1_i(:,:,js1))*ohs
393  zs12 = (z1_i(:,:,js) - z1_i(:,:,js1))*ohs
394  ru12 = p5*(ru_i(:,:,js) + ru_i(:,:,js1))
395  zu12 = p5*(zu_i(:,:,js) + zu_i(:,:,js1))
396  sqrtg_full(:,:,js) = r12*(ru12*zs12 - rs12*zu12)
397  END DO
398 
399  DEALLOCATE (r12, rs12, zs12, ru12, zu12)
400 
401  ming = minval(sqrtg_full(:,:,ns:nse) )
402  maxg = maxval(sqrtg_full(:,:,ns:nse))
403  print *,' IAM: ', iam, 'MING: ', ming,' MAXG: ',maxg
404  CALL assert(ming*maxg.GT.zero,' Jacobian changed sign in grid_extender')
405 
406  ALLOCATE (vp(nse))
407  CALL surfavglocal(vp,sqrtg_full,1,nse)
408  vp(1) = 0
409  DO js=2, nse
410  ming = minval(sqrtg_full(:,:,js) )
411  maxg = maxval(sqrtg_full(:,:,js))
412  WRITE(5000,'(i4,1p,3E14.4)') js, -vp(js), ming, maxg
413  END DO
414  CALL flush(5000)
415  DEALLOCATE (vp,r_ext, z_ext, ru_ext, zu_ext)
416  DEALLOCATE (ri,zi)
417 #endif
418 
419  DEALLOCATE(r1s_vmec, r1su_vmec, r1sv_vmec, z1s_vmec, z1su_vmec, z1sv_vmec)
420 
421  END SUBROUTINE grid_extender
422 
423  SUBROUTINE init_extender
424  REAL(dp) :: fac
425 
426  ALLOCATE(r1_vmec(ntheta,nzeta),z1_vmec(ntheta,nzeta))
427  ALLOCATE(ru_vmec(ntheta,nzeta),zu_vmec(ntheta,nzeta))
428  ALLOCATE(rv_vmec(ntheta,nzeta),zv_vmec(ntheta,nzeta))
429  r1_vmec = r1_i(:,:,ns); z1_vmec = z1_i(:,:,ns)
430  ru_vmec = ru_i(:,:,ns); zu_vmec(:,:) = zu_i(:,:,ns)
431  rv_vmec = rv_i(:,:,ns); zv_vmec(:,:) = zv_i(:,:,ns)
432 
433  !First derivatives on s=1 surface
434  ALLOCATE(r1s_vmec(ntheta,nzeta),z1s_vmec(ntheta,nzeta))
435  ALLOCATE(r1su_vmec(ntheta,nzeta),z1su_vmec(ntheta,nzeta))
436  ALLOCATE(r1sv_vmec(ntheta,nzeta),z1sv_vmec(ntheta,nzeta))
437 
438  fac =.5 !COULD SCAN THIS TO MAKE MAX(sqrt(g) < 0) IF NECESSARY
439 !! fac =.17 !COULD SCAN THIS TO MAKE MAX(sqrt(g) < 0) IF NECESSARY
440  r1s_vmec=(r1_vmec - r1_i(:,:,ns-1))*ohs*fac
441  z1s_vmec=(z1_vmec - z1_i(:,:,ns-1))*ohs*fac
442 
443  r1su_vmec = (ru_vmec - ru_i(:,:,ns-1))*ohs*fac
444  z1su_vmec = (zu_vmec - zu_i(:,:,ns-1))*ohs*fac
445 
446  r1sv_vmec = (rv_vmec - rv_i(:,:,ns-1))*ohs*fac
447  z1sv_vmec = (zv_vmec - zv_i(:,:,ns-1))*ohs*fac
448 
449  END SUBROUTINE init_extender
450 
451 !-------------------------------------------------------------------------------
455 !-------------------------------------------------------------------------------
456  FUNCTION getvp1()
457 
458  IMPLICIT NONE
459 
460 ! Declare Arguments
461  REAL(dp) :: getvp1
462 
463 ! Local Variables
464  INTEGER :: js
465  INTEGER :: js1
466  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: r12
467  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: rs12
468  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zs12
469  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ru12
470  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zu12
471  REAL(dp), DIMENSION(:), ALLOCATABLE :: vp
472 
473 ! Start of executable code.
474  ALLOCATE(vp(ns - 1:ns), sqrtg_full(ntheta,nzeta,ns - 1:ns))
475  ALLOCATE(r12(ntheta,nzeta), ru12(ntheta,nzeta), rs12(ntheta,nzeta), &
476  zu12(ntheta,nzeta), zs12(ntheta,nzeta))
477 
478 ! Guess assuming V(s) = V0*s**2: vp1 = 2*vol0
479 ! NOTE: jacobian not computed yet
480 
481  DO js = ns - 1, ns
482  js1 = js - 1
483  r12 = p5*(r1_i(:,:,js) + r1_i(:,:,js1))
484  rs12 = (r1_i(:,:,js) - r1_i(:,:,js1))*ohs
485  zs12 = (z1_i(:,:,js) - z1_i(:,:,js1))*ohs
486  ru12= p5*(ru_i(:,:,js) + ru_i(:,:,js1))
487  zu12= p5*(zu_i(:,:,js) + zu_i(:,:,js1))
488  sqrtg_full(:,:,js) = r12*(ru12*zs12 - rs12*zu12)
489  END DO
490 
491  CALL surfavglocal(vp, sqrtg_full, ns - 1, ns)
492  vp = twopi*twopi*abs(vp)
493  getvp1 = p5*(3*vp(ns) - vp(ns - 1))
494 
495  DEALLOCATE(vp, sqrtg_full, r12, ru12, rs12, zu12, zs12)
496 
497  END FUNCTION
498 
499  SUBROUTINE remap_radial(rho, vrho)
500  IMPLICIT NONE
501  INTEGER :: it
502  REAL(dp), INTENT(INOUT) :: rho
503  REAL(dp), INTENT(IN) :: vrho
504  REAL(dp) :: rootl, rootr, vol, delv
505  REAL(dp), PARAMETER :: tol=1.e-6_dp
506  !
507  ! ROOT FINDER: FINDS THE NEW rho VALUE CORRESPONDING TO THE
508  ! INITIAL VOLUME VRHO
509  !
510  rho = p5
511  rootr = 1
512  rootl = 0
513 
514  DO it = 1, 100
515  vol = get_volume(rho)
516  delv = vol-vrho
517  IF (abs(delv) .LT. tol*vrho) EXIT
518  IF (delv .GT. 0) THEN
519  rootr = rho
520  rho = (rootl+rho)/2
521  ELSE
522  rootl = rho
523  rho = rootl+(rootr-rho)/2
524  END IF
525  END DO
526  END SUBROUTINE remap_radial
527 
528  REAL(dp) function get_volume (rho)
529  IMPLICIT NONE
530  !
531  ! USES FORMULA VOLUME(rho) = INTEGRAL[.5*(R^2 * ZU)]
532  ! WHERE INTEGRATION IS OVER U,V FOR A FIXED rho=X
533  ! FOR SIMPLICITY, WE'LL DO THE INTEGRAL PLANE BY PLANE (WHICH IS NOT QUITE RIGHT)
534  !
535  REAL(dp) :: rho
536  INTEGER :: i, j
537  REAL(dp) :: zu, r1, dnorm, odelu
538 
539  get_volume = 0
540  odelu = ntheta
541  odelu = odelu/twopi
542  dnorm = p5*twopi/(nzeta*odelu)
543  IF (.NOT. lasym) odelu = 2*odelu
544  CALL loadsurfacepoints(rho)
545  DO j = 1, nzeta
546  DO i = 1, ntheta-1
547  zu = (zi(i+1,j)-zi(i,j))*odelu
548  r1 = p5*(ri(i+1,j)+ri(i,j))
549  get_volume = get_volume + r1**2 * zu
550  END DO
551  IF (lasym) THEN
552  zu = (zi(1,j)-zi(ntheta,j))*odelu
553  r1 = p5*(ri(1,j)+ri(ntheta,j))
554  get_volume = get_volume + r1**2*zu
555  END IF
556  END DO
557  get_volume = dnorm*get_volume
558 
559  END FUNCTION get_volume
560 
561  SUBROUTINE loadsurfacepoints(rho)
562  USE island_params, ns=>ns_i, mpol=>mpol_i, ntor=>ntor_i, &
563  ntheta=>nu_i, nzeta=>nv_i, mnmax=>mnmax_i, nsv=>ns_v
564  IMPLICIT NONE
565  REAL(dp), INTENT(IN) :: rho
566  REAL(dp) :: rho1, rho2, rho3
567 
568  !LOADS THE THETA AND ZETA VALUES OF RI, ZI ARRAYS
569  !FOR THE INPUT VALUE OF rho ("radial" surface)
570  rho1 = 1 - rho
571  IF (interp_type .EQ. lin) THEN
572  ri = rho1*r1_vmec + rho*r1_ves
573  zi = rho1*z1_vmec + rho*z1_ves
574  ELSE
575  rho2 = rho*rho
576  IF (interp_type .EQ. quad) THEN
577  !QUADRATIC (R' match at rho = 0)
578  ri = p5*r1s_vmec*(1 - rho2 - rho1*rho1) &
579  + r1_vmec*(1-rho2) + r1_ves*rho2
580  zi = p5*z1s_vmec*(1 - rho2 - rho1*rho1) &
581  + z1_vmec*(1-rho2) + z1_ves*rho2
582  ELSE IF (interp_type .EQ. cub) THEN
583  !CUBIC (R', R'' match at rho=0)
584  rho3 = rho*rho2
585  CALL assert(.false., 'Unknown interpolation <TYPE> specified')
586  END IF
587  END IF
588  END SUBROUTINE loadsurfacepoints
589 
590 !-------------------------------------------------------------------------------
594 !-------------------------------------------------------------------------------
595  FUNCTION read_vessel_file()
596  USE island_params, ns=>ns_i, mpol=>mpol_i, ntor=>ntor_i, &
597  ntheta=>nu_i, nzeta=>nv_i, mnmax=>mnmax_i
598  USE siesta_namelist, ONLY: vessel_file
599 
600  IMPLICIT NONE
601 
602 ! Declare Arguments
603  INTEGER :: read_vessel_file
604 
605 ! Local variables
606  INTEGER :: m
607  INTEGER :: n
608  INTEGER :: mn
609  INTEGER :: mnmax_v
610  INTEGER :: iou
611  REAL (dp) :: rv1
612  REAL (dp) :: zv1
613 
614 ! Start of executable code.
615  mpol_v = 0
616  ntor_v = 0
617  mnmax_v = 0
618 
619  iou = 36
620  OPEN(unit=iou, file=vessel_file, status="old", form="formatted", &
621  iostat=read_vessel_file)
622  IF (read_vessel_file .NE. 0) THEN
623  WRITE (*,1000) vessel_file
624  RETURN
625  END IF
626 
627  DO WHILE(.true.)
628  READ(iou,*,END=100) n, m, rv1, zv1
629  mpol_v = max(m, mpol_v)
630  ntor_v = max(abs(n), ntor_v)
631  mnmax_v = mnmax_v + 1
632  END DO
633  100 CONTINUE
634 
635  rewind(unit=iou, iostat=read_vessel_file)
636  IF (read_vessel_file .NE. 0) THEN
637  WRITE (*,*) vessel_file
638  RETURN
639  END IF
640  ALLOCATE(rmnc_v(0:mpol_v,-ntor_v:ntor_v), &
641  zmns_v(0:mpol_v,-ntor_v:ntor_v))
642  rmnc_v = 0
643  zmns_v = 0
644  IF (lasym) THEN
645  ALLOCATE(rmns_v(0:mpol_v,-ntor_v:ntor_v), &
646  zmnc_v(0:mpol_v,-ntor_v:ntor_v))
647  rmns_v = 0
648  zmnc_v = 0
649  END IF
650 
651  DO mn = 1, mnmax_v
652  READ(iou,*,END=101) n, m, rv1, zv1
653  CALL assert(m .ge. 0 .and. m .le. mpol_v, &
654  'm out of bounds: read_vessel_file')
655  CALL assert(n .ge. -ntor_v .and. n .le. ntor_v, &
656  'n out of bounds: read_vessel_file')
657  rmnc_v(m,-n) = rv1
658  zmns_v(m,-n) = zv1
659  END DO
660  101 CONTINUE
661  CLOSE(unit=iou)
662 
663 1000 FORMAT('Warning: Failed to open vessel file ',a)
664 1001 FORMAT('Warning: Failed to read vessel file ',a)
665 
666  END FUNCTION
667 
668 !-------------------------------------------------------------------------------
670 !-------------------------------------------------------------------------------
671  SUBROUTINE surfavglocal(average, q3d, nsmin, nsmax)
672  USE island_params, ONLY: fourier_context
673 
674  IMPLICIT NONE
675 
676 ! Declare Arguments
677  REAL(dp), INTENT(out) :: average(nsmin:nsmax)
678  REAL(dp), INTENT(in) :: q3d(ntheta,nzeta,nsmin:nsmax)
679  INTEGER, INTENT(in) :: nsmin
680  INTEGER, INTENT(in) :: nsmax
681 
682 ! Local Variables
683  INTEGER :: i
684  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: wint
685 
686 ! Start of executable code.
687  ALLOCATE(wint(ntheta,nzeta,nsmin:nsmax))
688 
689  DO i = 1, ntheta
690  wint(i,:,:) = fourier_context%cosmui(i,0)
691  END DO
692  DO i = nsmin, nsmax
693  average(i) = sum(q3d(:,:,i)*wint(:,:,i))
694  END DO
695 
696  DEALLOCATE(wint)
697 
698  END SUBROUTINE
699 
700  END MODULE
shared_data::r1_i
real(dp), dimension(:,:,:), allocatable r1_i
R coordinates of the computational grid.
Definition: shared_data.f90:153
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
shared_data::zv_i
real(dp), dimension(:,:,:), allocatable zv_i
dZ/dv coordinates of the computational grid.
Definition: shared_data.f90:163
siesta_namelist::l_vessel
logical l_vessel
If extended grid is to be used using an available vessel file.
Definition: siesta_namelist.f90:140
siesta_namelist::vessel_file
character(len=siesta_namelist_name_length) vessel_file
Name of the restart file extension.
Definition: siesta_namelist.f90:195
grid_extension::surfavglocal
subroutine surfavglocal(average, q3d, nsmin, nsmax)
Get surface average quantity.
Definition: grid_extension.f90:672
island_params::phipf_i
real(dp), dimension(:), allocatable phipf_i
Radial toroidal flux derivative.
Definition: island_params.f90:68
grid_extension::lin
integer, parameter lin
Linear extension type.
Definition: grid_extension.f90:29
island_params::fourier_context
type(fourier_class), pointer fourier_context
Fourier transform object.
Definition: island_params.f90:76
fourier
Module contains subroutines for computing FFT on parallel radial intervals. Converts quantities betwe...
Definition: fourier.f90:13
v3_utilities::assert
Definition: v3_utilities.f:55
fourier::f_sin
integer, parameter f_sin
Sine parity.
Definition: fourier.f90:27
grid_extension::getvp1
real(dp) function getvp1()
UNKNOWN.
Definition: grid_extension.f90:457
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
shared_data::zu_i
real(dp), dimension(:,:,:), allocatable zu_i
dZ/du coordinates of the computational grid.
Definition: shared_data.f90:159
shared_data::lasym
logical lasym
Use non-stellarator symmetry.
Definition: shared_data.f90:230
island_params::chipf_i
real(dp), dimension(:), allocatable chipf_i
Radial poloidal flux derivative.
Definition: island_params.f90:70
island_params::presf_i
real(dp), dimension(:), allocatable presf_i
Radial pressure. FIXME: Check if this is really needed.
Definition: island_params.f90:72
grid_extension::quad
integer, parameter quad
Quadtratic extension type.
Definition: grid_extension.f90:31
shared_data::rv_i
real(dp), dimension(:,:,:), allocatable rv_i
dR/dv coordinates of the computational grid.
Definition: shared_data.f90:161
fourier::f_du
integer, parameter f_du
Poloidal derivative flag.
Definition: fourier.f90:36
shared_data::nsp
integer nsp
Total radial grid size in the VMEC region.
Definition: shared_data.f90:64
grid_extension::read_vessel_file
integer function read_vessel_file()
Read the vessel file.
Definition: grid_extension.f90:596
grid_extension::cub
integer, parameter cub
Cubic extension type.
Definition: grid_extension.f90:33
fourier::f_sum
integer, parameter f_sum
Sum fouier real space transformed quantity. This is used when a non-stellarator symmetric parity is b...
Definition: fourier.f90:43
grid_extension
This file contains utilities for extending the vmec grid.
Definition: grid_extension.f90:10
island_params::ohs_i
real(dp) ohs_i
Radial grid derivative factor. d X/ ds = (x_i+1 - x_i)/(s_i+1 - s_i) (1) Where ohs = 1/(s_i+1 - s_i) ...
Definition: island_params.f90:51
grid_extension::p5
real(dp), parameter p5
Divide by two constant.
Definition: grid_extension.f90:35
shared_data::z1_i
real(dp), dimension(:,:,:), allocatable z1_i
Z coordinates of the computational grid.
Definition: shared_data.f90:155
fourier::f_cos
integer, parameter f_cos
Cosine parity.
Definition: fourier.f90:25
siesta_namelist
This file contains all the variables and maximum sizes of the inputs for a SIESTA namelist input file...
Definition: siesta_namelist.f90:103
shared_data::ru_i
real(dp), dimension(:,:,:), allocatable ru_i
dR/du coordinates of the computational grid.
Definition: shared_data.f90:157
shared_data
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.f90:10
siesta_namelist::nsin_ext
integer nsin_ext
Radial size of the extended grid.
Definition: siesta_namelist.f90:164
fourier::f_dv
integer, parameter f_dv
Toroidal derivative flag.
Definition: fourier.f90:38
fourier::fourier_class
Base class containing fourier memory.
Definition: fourier.f90:77
fourier::f_none
integer, parameter f_none
Do not sum fouier real space transformed quantity. This is used when a stellarator symmetric parity i...
Definition: fourier.f90:32
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17