V3FIT
vacfield_mod.f
1  MODULE vacfield_mod
2  USE stel_kinds
3  IMPLICIT NONE
4  INTEGER, PARAMETER :: ncur = 100, msmax = 64, nresmax = 10
5  INTEGER, PARAMETER :: iout = 10, jout = 11, kout = 12, lout = 13
6  INTEGER :: nfp, ngroups, ntmax, npmax, lopt, nres, numres
7  INTEGER, DIMENSION(0:nresmax) :: mres, nimax, nrp
8  REAL(rprec) :: drsurf, residsum, raxiscon
9  REAL(rprec), DIMENSION(msmax) :: ri, zi, di, drho, iota
10  REAL(rprec), DIMENSION(0:nresmax) :: vres, rres, dres
11  REAL(rprec), DIMENSION(0:nresmax) :: trace, deter, resid
12  !trace, determinant, residue of map
13  REAL(rprec), DIMENSION(0:nresmax) :: rres_min, rres_max, reswt
14  REAL(rprec), DIMENSION(1:nresmax) :: rmax_isl, rmin_isl
15  REAL(rprec), DIMENSION(1:nresmax) :: zmax_isl, zmin_isl
16  REAL(rprec), DIMENSION(1:nresmax) :: shear_r, delu_r
17  REAL(rprec), DIMENSION(1:nresmax) :: shear_p, delu_p
18  REAL(rprec), DIMENSION(1:nresmax) :: area_isl, flux_isl
19  REAL(rprec), DIMENSION(1:nresmax) :: alpha
20  REAL(rprec), DIMENSION(1:nresmax) :: sigma_resid
21  LOGICAL, DIMENSION(1:nresmax) :: lsymm
22  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: ropnts, zopnts
23  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: r_isl, z_isl
24  LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lli, lri
25  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: rt, zt
26  REAL(rprec), DIMENSION(msmax) :: rmag, zmag
27  INTEGER :: ierr, mc, msurf, nvar, nstell_coils
28  INTEGER :: ierr_mpi, myid, numprocs
29  REAL(rprec), DIMENSION(ncur) :: vac_extcur, xc
30  REAL(rprec) :: raxis, zaxis, r0, z0, dsmin, dsurf, tsurf, v0
31  REAL(rprec), DIMENSION(nresmax) :: atol, rtol
32  LOGICAL :: laxis, lprt, lisize, lsurf, ldisplace, lscreen
33  LOGICAL, DIMENSION(:), ALLOCATABLE :: lopnt
34 
35  TYPE angles_shifts
36  REAL(rprec) :: as_array(3)
37  END TYPE angles_shifts
38 
39 
40  TYPE (angles_shifts), DIMENSION(20) :: angles, angles_lo,
41  1 angles_hi, sigma_angles_lo, sigma_angles_hi
42  TYPE (angles_shifts), DIMENSION(20) :: shifts, shifts_lo,
43  1 shifts_hi, sigma_shifts_lo, sigma_shifts_hi
44  TYPE (angles_shifts), DIMENSION(20,10,10) ::
45  1 control_pt_init, control_pt_final
46  CHARACTER(LEN=100) :: extension, coils_file_extension
47 
48  namelist /nlst_vacopt/ vac_extcur, msurf, ntmax, ri, zi, lopt,
49  1 r0, z0, v0, dsmin, dsurf, tsurf, nfp, laxis, raxis, zaxis,
50  2 vres, rres, rres_min, rres_max, reswt, nrp, dres, mres, numres,
51  3 raxiscon, rmin_isl, rmax_isl, zmin_isl, zmax_isl, shear_r,
52  4 shear_p, alpha, lsymm, lisize, lsurf, lprt, ldisplace,
53  5 nstell_coils, nvar, rtol, atol, coils_file_extension,
54  6 rmag, zmag, sigma_resid,
55  7 angles, angles_lo, angles_hi, sigma_angles_lo, sigma_angles_hi,
56  8 shifts, shifts_lo, shifts_hi, sigma_shifts_lo, sigma_shifts_hi,
57  9 control_pt_init
58 
59  INTERFACE write_invac
60  MODULE PROCEDURE write_invac, write_invac2
61  END INTERFACE
62 
63  CONTAINS
64 
65  SUBROUTINE initialize_vacfield (arg1)
66  USE safe_open_mod
67  USE biotsavart
68 
69  INTEGER :: n, np, istat, invac
70  CHARACTER*(*), INTENT(in) :: arg1
71  CHARACTER :: invac_file*100
72 
73 ! Get extension of invac file, used for writing output files
74  extension = trim(arg1)
75  n = index(extension, ".", back=.true.)
76  extension = extension(n+1:)
77  coils_file_extension = trim(extension)
78 
79 ! Read namelist input file
80  invac = 20
81  vac_extcur = 0
82  invac_file = arg1
83 
84  CALL safe_open (invac, istat, invac_file, 'old', 'formatted')
85  IF (istat .ne. 0) THEN
86  print *, 'Error opening namelist input file', invac_file
87  stop 02
88  END IF
89 
90  CALL read_vacopt_namelist(invac, istat)
91 
92  IF (istat .ne. 0) THEN
93  print *,' Error reading ', trim(invac_file),
94  1 ' iostat = ', istat
95  print *,' Be sure EXTCUR -> VAC_EXTUR!'
96  stop 03
97  END IF
98 
99  CLOSE(unit=invac)
100 
101 ! Number of toroidal symmetry planes
102  npmax = 2*nfp
103 
104 ! Read in coils-dot file
105  CALL initialize_biotsavart (vac_extcur, coils_file_extension,
106  1 scaled=.false.)
107 
108  END SUBROUTINE initialize_vacfield
109 
110  SUBROUTINE read_vacopt_namelist (iunit, istat)
111  INTEGER :: iunit, istat
112  INTEGER :: n
113 !
114 ! READS IN nls_vacopt NAMELIST
115 !
116  ldisplace = .false. !.true, will apply rotations, shifts to produce new coils-dot file
117  laxis = .true.
118  lprt = .false.
119  lisize = .false.
120  lsurf = .false.
121  msurf = 16
122  ntmax = 200
123  nfp = 2
124  lopt = 0
125  nvar = 10
126  raxis = 0.7173_dp
127  raxiscon = 0.7173_dp
128  zaxis = 0
129  rmag(1) = 0.7173_dp
130  rmag(2) = 1.2945_dp
131  rmag(3) = 0.7173_dp
132  rmag(4) = 1.2945_dp
133  zmag(1) = 0.0_dp
134  zmag(2) = 0.0_dp
135  zmag(3) = 0.0_dp
136  zmag(4) = 0.0_dp
137  r0 = 0.7173_dp
138  z0 = 0.0_dp
139  v0 = 0.5_dp
140  numres = 3 ! number of targeted resonances
141  sigma_resid = 1 ! sigma (inverse weight) associated with targeted resonances
142  DO n = 1, SIZE(angles,1)
143  angles(n) = angles_shifts( (/0._dp, 0._dp, 0._dp/) )
144  shifts(n) = angles_shifts( (/0._dp, 0._dp, 0._dp/) )
145  sigma_angles_lo(n) = angles_shifts( (/1._dp, 1._dp, 1._dp/) )
146  sigma_shifts_lo(n) = angles_shifts( (/1._dp, 1._dp, 1._dp/) )
147  sigma_angles_hi(n) = angles_shifts( (/1._dp, 1._dp, 1._dp/) )
148  sigma_shifts_hi(n) = angles_shifts( (/1._dp, 1._dp, 1._dp/) )
149 
150  !control_pt arrays: fixed pts on coils; init is the initial value (x,y,z)
151  ! final is the rotated+shifted pt (x,y,z) written to chisq file
152  control_pt_init(n,:,:) =
153  1 angles_shifts( (/0._dp, 0._dp, 0._dp/) )
154  END DO
155  mres(0) = 1 ! order of fixed point (magnetic axis)
156  mres(1) = 10 ! mode number for targeted resonance
157  mres(2) = 9 ! mode number for targeted resonance
158  mres(3) = 8 ! mode number for targeted resonance
159  vres(0) = 0.5_dp ! toroidal plane for axis search
160  vres(1) = 0.0_dp ! toroidal plane for O-point search
161  vres(2) = 0.5_dp ! toroidal plane for O-point search
162  vres(3) = 0.0_dp ! toroidal plane for O-point search
163  rres(0) = 0.72_dp ! approx. radius of O-point
164  rres(1) = 1.33_dp ! approx. radius of O-point
165  rres(2) = 0.90_dp ! approx. radius of O-point
166  rres(3) = 1.37_dp ! approx. radius of O-point
167  nrp(0) = 1 ! no. subintervals in fixed pt. search
168  nrp(1) = 1 ! no. subintervals in fixed pt. search
169  nrp(2) = 1 ! no. subintervals in fixed pt. search
170  nrp(3) = 1 ! no. subintervals in fixed pt. search
171  rres_min(0) = 0.70_dp ! min. R for O-point search
172  rres_max(0) = 0.73_dp ! max. R for O-point search
173  rres_min(1) = 1.33_dp ! min. R for O-point search
174  rres_max(1) = 1.36_dp ! max. R for O-point search
175  rres_min(2) = 0.88_dp ! min. R for O-point search
176  rres_max(2) = 0.94_dp ! max. R for O-point search
177  rres_min(3) = 1.37_dp ! min. R for O-point search
178  rres_max(3) = 1.40_dp ! max. R for O-point search
179  reswt(0) = 1.0_dp ! weight for residue minimization
180  reswt(1) = 1.0_dp ! weight for residue minimization
181  reswt(2) = 1.0_dp ! weight for residue minimization
182  reswt(3) = 1.0_dp ! weight for residue minimization
183  dres(0) = 0.010_dp ! step size for tangent map approx.
184  dres(1) = 0.005_dp ! step size for tangent map approx.
185  dres(2) = 0.010_dp ! step size for tangent map approx.
186  dres(3) = 0.005_dp ! step size for tangent map approx.
187  vac_extcur(1) = 0 ! vacuum external current array
188  rmin_isl(1) = 1.33_dp ! min. R for island size calc.
189  rmax_isl(1) = 1.36_dp ! max. R for island size calc.
190  zmin_isl(1) = -0.40_dp ! min. Z for island size calc.
191  zmax_isl(1) = 0.40_dp ! max. Z for island size calc.
192  rmin_isl(2) = 0.88_dp ! min. R for island size calc.
193  rmax_isl(2) = 0.94_dp ! max. R for island size calc.
194  zmin_isl(2) = -0.10_dp ! min. Z for island size calc.
195  zmax_isl(2) = 0.10_dp ! max. Z for island size calc.
196  rmin_isl(3) = 1.37_dp ! min. R for island size calc.
197  rmax_isl(3) = 1.40_dp ! max. R for island size calc.
198  zmin_isl(3) = -0.40_dp ! min. Z for island size calc.
199  zmax_isl(3) = 0.40_dp ! max. Z for island size calc.
200  shear_r(1) = 0.56_dp ! shear (1/m) of unperturbed flow at res.
201  shear_r(2) = 0.56_dp ! shear (1/m) of unperturbed flow at res.
202  shear_r(3) = 0.53_dp ! shear (1/m) of unperturbed flow at res.
203  shear_p(1) = -0.3875_dp ! shear (1/Wb) of unperturbed flow at res.
204  shear_p(2) = -0.3875_dp ! shear (1/Wb) of unperturbed flow at res.
205  shear_p(3) = -0.2813_dp ! shear (1/Wb) of unperturbed flow at res.
206  lsymm(1) = .false. ! = true if starting in elongated plane
207  lsymm(2) = .false. ! = true if starting in elongated plane
208  lsymm(3) = .false. ! = true if starting in elongated plane
209  alpha(1) = 1.0_dp ! factor to adjust du(m) for starting pt.
210  alpha(2) = 1.0_dp ! factor to adjust du(m) for starting pt.
211  alpha(3) = 1.0_dp ! factor to adjust du(m) for starting pt.
212  dsmin = 0.00_dp
213  dsurf = 0.01_dp
214  tsurf = 0.00_dp
215  rtol = 1.e-8_dp
216  atol = 1.e-8_dp
217 
218  DO n = 1, msmax
219  di(n) = dsmin + n*dsurf
220  ri(n) = raxis + di(n)*cos(tsurf)
221  zi(n) = zaxis + di(n)*sin(tsurf)
222  END DO
223 
224  nstell_coils = -1 !No. of unique "stellarator-symmetric" coils to rotate
225 
226  READ (iunit, nml=nlst_vacopt, iostat=istat)
227 
228  END SUBROUTINE read_vacopt_namelist
229 
230  SUBROUTINE write_invac (iunit, istat)
231  USE stel_constants
232  USE write_array_generic
233  IMPLICIT NONE
234 !-----------------------------------------------
235 ! D u m m y A r g u m e n t s
236 !-----------------------------------------------
237  INTEGER :: iunit, istat
238 !-----------------------------------------------
239 ! L o c a l V a r i a b l e s
240 !-----------------------------------------------
241  INTEGER :: i, j, k
242  CHARACTER(LEN=100) :: form1, form2
243 !-----------------------------------------------
244 
245  WRITE(iunit,'(a12)')'&NLST_VACOPT'
246  WRITE(iunit,100) 'LDISPLACE = ', ldisplace
247  WRITE(iunit,100) 'LAXIS = ', laxis
248  WRITE(iunit,100) 'LPRT = ', lprt
249  WRITE(iunit,100) 'LISIZE = ', lisize
250  WRITE(iunit,200) 'LOPT = ', lopt
251  WRITE(iunit,200) 'MSURF = ', msurf
252  WRITE(iunit,200) 'NTMAX = ', ntmax
253  WRITE(iunit,200) 'NFP = ', nfp
254  WRITE(iunit,200) 'NSTELL_COILS = ', nstell_coils
255  WRITE(iunit,200) 'NVAR = ', nvar
256  WRITE(iunit,300) 'RAXIS = ', raxis
257  WRITE(iunit,300) 'R0 = ', r0
258  WRITE(iunit,300) 'Z0 = ', z0
259  WRITE(iunit,300) 'V0 = ', v0
260  WRITE(iunit,300) 'DSURF = ', dsurf
261  WRITE(iunit,300) 'DSMIN = ', dsmin
262  WRITE(iunit,300) 'TSURF = ', tsurf
263  WRITE(iunit,200) 'NUMRES = ', numres
264  WRITE(iunit,'(1x,3a)') "COILS_FILE_EXTENSION = '",
265  1 trim(coils_file_extension),"'"
266  IF (msmax .gt. 0) THEN
267  CALL write_array(iunit,'RMAG', rmag, msmax)
268  CALL write_array(iunit,'ZMAG', zmag, msmax)
269  CALL write_array(iunit,'RI', ri, msurf)
270  CALL write_array(iunit,'ZI', zi, msurf)
271  END IF
272  IF (nresmax .gt. 0) THEN
273  CALL write_array(iunit,'MRES', mres, nresmax+1, low_index=0)
274  CALL write_array(iunit,'VRES', vres, nresmax+1, low_index=0)
275  CALL write_array(iunit,'RRES', rres, nresmax+1, low_index=0)
276  CALL write_array(iunit,'NRP', nrp, nresmax+1, low_index=0)
277  CALL write_array(iunit,'RRES_MIN', rres_min, nresmax+1,
278  1 low_index=0)
279  CALL write_array(iunit,'RRES_MAX', rres_max, nresmax+1,
280  1 low_index=0)
281  CALL write_array(iunit,'RESWT', reswt, nresmax+1, low_index=0)
282  CALL write_array(iunit,'DRES', dres, nresmax+1, low_index=0)
283  END IF
284  IF (ncur .gt. 0) THEN
285  CALL write_array(iunit,'VAC_EXTCUR', vac_extcur, ncur)
286  END IF
287  IF (nresmax .gt. 0) THEN
288  CALL write_array(iunit,'RMIN_ISL', rmin_isl, nresmax)
289  CALL write_array(iunit,'RMAX_ISL', rmax_isl, nresmax)
290  CALL write_array(iunit,'ZMIN_ISL', zmin_isl, nresmax)
291  CALL write_array(iunit,'ZMAX_ISL', zmax_isl, nresmax)
292  CALL write_array(iunit,'SHEAR_R', shear_r, nresmax)
293  CALL write_array(iunit,'SHEAR_P', shear_p, nresmax)
294  CALL write_array(iunit,'LSYMM', lsymm, nresmax)
295  CALL write_array(iunit,'ALPHA', alpha, nresmax)
296  CALL write_array(iunit,'SIGMA_RESID', sigma_resid, nresmax)
297  END IF
298 
299  DO i = 1, nstell_coils
300  IF (i .lt. 10) THEN
301  form1 = '(2x,a,i1,a,1p3e15.6)'
302  form2 = '(2x,a,i1,2(a1,i2.2),a,1p3e15.6)'
303  ELSE
304  form1 = '(2x,a,i2,a,1p3e15.6)'
305  form2 = '(2x,a,i2,2(a1,i2.2),a,1p3e15.6)'
306  END IF
307  WRITE(iunit, form1, iostat=istat)
308  1 'ANGLES(',i,') = ', angles(i)
309  WRITE(iunit, form1, iostat=istat)
310  1 'SHIFTS(',i,') = ', shifts(i)
311  WRITE(iunit, form1, iostat=istat)
312  1 'ANGLES_LO(',i,') = ', angles_lo(i)
313  WRITE(iunit, form1, iostat=istat)
314  1 'SHIFTS_LO(',i,') = ', shifts_lo(i)
315  WRITE(iunit, form1, iostat=istat)
316  1 'ANGLES_HI(',i,') = ', angles_hi(i)
317  WRITE(iunit, form1, iostat=istat)
318  1 'SHIFTS_HI(',i,') = ', shifts_hi(i)
319  WRITE(iunit, form1, iostat=istat)
320  1 'SIGMA_ANGLES_LO(',i,') = ', sigma_angles_lo(i)
321  WRITE(iunit, form1, iostat=istat)
322  1 'SIGMA_SHIFTS_LO(',i,') = ', sigma_shifts_lo(i)
323  WRITE(iunit, form1, iostat=istat)
324  1 'SIGMA_ANGLES_HI(',i,') = ', sigma_angles_hi(i)
325  WRITE(iunit, form1, iostat=istat)
326  1 'SIGMA_SHIFTS_HI(',i,') = ', sigma_shifts_hi(i)
327 
328  DO j = 1, SIZE(control_pt_init,2)
329  DO k = 1, SIZE(control_pt_init,3)
330  WRITE(iunit, form2, iostat=istat)
331  1 'CONTROL_PT_INIT(',i,',',j,',',k,') = ',
332  2 control_pt_init(i,j,k)
333  END DO
334  END DO
335 
336  END DO
337 
338  WRITE(iunit,'(a)')'/'
339 
340  100 FORMAT(4(1x,a,l2,','))
341  200 FORMAT(4(1x,a,i6,','))
342  300 FORMAT(3(1x,a,1pe12.4,','))
343  350 FORMAT(2(1x,a,1pe12.4,','))
344  450 FORMAT(1x,a,1pe12.4,',')
345 
346  END SUBROUTINE write_invac
347 
348  SUBROUTINE write_invac2 (filename, istat)
349  USE stel_constants
350  USE write_array_generic
351  USE safe_open_mod
352  IMPLICIT NONE
353 !-----------------------------------------------
354 ! D u m m y A r g u m e n t s
355 !-----------------------------------------------
356  INTEGER :: iunit, istat
357  CHARACTER(LEN=*) :: filename
358 !-----------------------------------------------
359  CALL safe_open(iunit, istat, filename, 'replace', 'formatted')
360  IF (istat .ne. 0) RETURN
361 
362  CALL write_invac (iunit, istat)
363 
364  CLOSE (iunit)
365 
366  END SUBROUTINE write_invac2
367 
368  END MODULE vacfield_mod
vacfield_mod::angles_shifts
Definition: vacfield_mod.f:35
vacfield_mod::write_invac
Definition: vacfield_mod.f:59
write_array_generic::write_array
Definition: write_array_generic.f:11
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17