V3FIT
symcheck.f
1  MODULE sym_check
2 
3  USE biotsavart
4  USE makegrid_global, ONLY: rprec
5  USE write_mgrid, ONLY: rmin, rmax, zmin, zmax
6 
7  IMPLICIT NONE
8 
9  PRIVATE
10 
11 !-----------------------------------------------
12 ! L o c a l V a r i a b l e s
13 !-----------------------------------------------
14  INTEGER :: sym_nfp, sym_nextcur
15  REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:) ::
16  1 sym_br, sym_bz, sym_bp, sym_b
17 !-----------------------------------------------
18 ! P u b l i c S u b r o u t i n e s A n d V a r i a b l e s
19 !-----------------------------------------------
20  PUBLIC :: check_symmetry
21  PUBLIC :: init_symmetry
22  PUBLIC :: cleanup_symmetry
23 
24  INTEGER, PUBLIC :: sym_ir, sym_jz, sym_kp
25  LOGICAL, PUBLIC :: sym_perform_tests
26 !-----------------------------------------------
27 
28  CONTAINS
29 
30 !----------------------------------------------------------------------
31 !*******************************************************************************
32 !----------------------------------------------------------------------
33 
34  SUBROUTINE init_symmetry
35 
36  sym_nfp = nfp_bs
37  sym_nextcur = SIZE(coil_group)
38 
39  ALLOCATE( sym_br(sym_ir, sym_jz, sym_kp*sym_nfp) )
40  ALLOCATE( sym_bz(sym_ir, sym_jz, sym_kp*sym_nfp) )
41  ALLOCATE( sym_bp(sym_ir, sym_jz, sym_kp*sym_nfp) )
42  ALLOCATE( sym_b(sym_ir, sym_jz, sym_kp*sym_nfp) )
43 
44  END SUBROUTINE init_symmetry
45 
46 !----------------------------------------------------------------------
47 !*******************************************************************************
48 !----------------------------------------------------------------------
49 
50  SUBROUTINE cleanup_symmetry
51 
52  DEALLOCATE(sym_bp)
53  DEALLOCATE(sym_bz)
54  DEALLOCATE(sym_br)
55  DEALLOCATE(sym_b)
56 
57  END SUBROUTINE cleanup_symmetry
58 
59 !----------------------------------------------------------------------
60 !*******************************************************************************
61 !----------------------------------------------------------------------
62 
63  SUBROUTINE check_symmetry(ig)
64 !-----------------------------------------------
65 ! D u m m y A r g u m e n t s
66 !-----------------------------------------------
67  INTEGER, INTENT(in) :: ig
68 !-----------------------------------------------
69 ! L o c a l V a r i a b l e s
70 !-----------------------------------------------
71  INTEGER :: i, j, ki, fp
72  REAL(rprec) :: maxerr, avgerr
73 !-----------------------------------------------
74 
75  IF(sym_perform_tests) THEN
76 
77  CALL sym_compute_bfield(ig)
78 
79  WRITE(*,225)
80  DO ki=1, sym_kp
81  CALL sym_checkrot(ki,maxerr,avgerr)
82  WRITE(*,300) ki, maxerr, avgerr
83  END DO
84 
85  WRITE(*,275)
86  DO fp=1, sym_nfp
87  CALL sym_checkstel(fp,maxerr,avgerr)
88  WRITE(*,300) fp, maxerr, avgerr
89  END DO
90 
91  END IF
92 
93  225 FORMAT(/" Error factors, rotational symmetry ",
94  1 "(corresponding to KI'th plane):",/,
95  2 " KI", 6x, "Max Error", 5x, "Avg Error")
96 
97  275 FORMAT(/" Error factors, stellarator symmetry ",
98  1 "(for FP'th field period):",/,
99  2 " FP", 6x, "Max Error", 5x, "Avg Error")
100 
101  300 FORMAT(i2,4x,es14.5,es14.5)
102 
103  END SUBROUTINE check_symmetry
104 
105 !----------------------------------------------------------------------
106 !*******************************************************************************
107 !----------------------------------------------------------------------
108  SUBROUTINE sym_compute_bfield(ig)
109 !-----------------------------------------------
110 ! D u m m y A r g u m e n t s
111 !-----------------------------------------------
112  INTEGER, INTENT(in) :: ig
113 !-----------------------------------------------
114 ! L o c a l V a r i a b l e s
115 !-----------------------------------------------
116  INTEGER :: numcoils, i, j, k
117  REAL(rprec) :: rrr,ppp,zzz,delr,delp,delz
118 !-----------------------------------------------
119 
120  delr = (rmax - rmin)/(sym_ir - 1)
121  delz = (zmax - zmin)/(sym_ir - 1)
122  delp = (8*datan(1.0d0))/(sym_nfp*sym_kp)
123 
124  sym_bz = 0d0
125  sym_bp = 0d0
126  sym_br = 0d0
127 
128  ppp = 0d0
129  DO k=1,sym_kp*sym_nfp
130  zzz = zmin
131 
132  DO j=1,sym_jz
133  rrr = rmin
134 
135  DO i=1,sym_ir
136  ! Compute br, bp, bz at each grid point
137  CALL bfield (rrr, ppp, zzz, sym_br(i,j,k),
138  1 sym_bp(i,j,k), sym_bz(i,j,k), ig)
139 
140  ! Compute magnitude, b, at each grid point
141  sym_b(i,j,k) =
142  1 dsqrt(sym_br(i,j,k)**2 + sym_bp(i,j,k)**2 + sym_bz(i,j,k)**2)
143 
144  rrr = rrr + delr
145  END DO
146  zzz = zzz + delz
147  END DO
148  ppp = ppp + delp
149  END DO
150 
151  END SUBROUTINE sym_compute_bfield
152 
153 !----------------------------------------------------------------------
154 !*******************************************************************************
155 !----------------------------------------------------------------------
156 
157  SUBROUTINE sym_checkrot(ki,maxerr,avgerr)
158 !-----------------------------------------------
159 ! D u m m y A r g u m e n t s
160 !-----------------------------------------------
161  INTEGER, INTENT(in) :: ki ! 1 <= ki <= sym_kp
162  REAL(rprec), INTENT(out) :: maxerr, avgerr
163 !-----------------------------------------------
164 ! L o c a l V a r i a b l e s
165 !-----------------------------------------------
166  INTEGER :: i,j,k,l
167  REAL(rprec) :: dbr, dbp, dbz, currerr, avgb
168 !-----------------------------------------------
169 
170  maxerr = 0d0
171  avgerr = 0d0
172 
173  ! Compute maximum and average rotational error factors
174  ! for the set of grids corresponding to k=ki
175  DO k=ki+sym_kp, sym_kp*sym_nfp, sym_kp
176  currerr = 0d0
177  DO j=1, sym_jz
178  DO i=1, sym_ir
179 
180  dbr = sym_br(i,j,k) - sym_br(i,j,ki)
181  dbp = sym_bp(i,j,k) - sym_bp(i,j,ki)
182  dbz = sym_bz(i,j,k) - sym_bz(i,j,ki)
183  avgb = (sym_b(i,j,k) + sym_b(i,j,ki))/2.0d0
184 
185  IF(avgb .NE. 0d0) THEN
186  currerr = currerr + dsqrt(dbr**2+dbp**2+dbz**2)/avgb
187  END IF
188  END DO
189  END DO
190 
191  avgerr = avgerr + currerr
192  IF(currerr .GT. maxerr) maxerr = currerr
193  END DO
194 
195  IF (sym_nfp .GT. 1) avgerr = avgerr/(sym_nfp-1)
196 
197  END SUBROUTINE sym_checkrot
198 
199 !----------------------------------------------------------------------
200 !*******************************************************************************
201 !----------------------------------------------------------------------
202 
203  SUBROUTINE sym_checkstel(fp,maxerr,avgerr)
204 !-----------------------------------------------
205 ! D u m m y A r g u m e n t s
206 !-----------------------------------------------
207  INTEGER, INTENT(in) :: fp ! 1 <= fp <= sym_nfp
208  REAL(rprec), INTENT(out) :: maxerr, avgerr
209 !-----------------------------------------------
210 !-----------------------------------------------
211 ! L o c a l V a r i a b l e s
212 !-----------------------------------------------
213  INTEGER :: i,j,k,h,l, startk, midk, endk
214  REAL(rprec) :: dbr, dbp, dbz, avgb, currerr
215 !-----------------------------------------------
216 
217  startk = sym_kp*(fp-1)+1
218  midk = startk+sym_kp/2-1
219  endk = sym_kp*fp+1
220 
221  maxerr = 0d0
222  avgerr = 0d0
223 
224  ! Compute maximum stellarator symmetry error factor in the
225  ! fp'th field period
226  DO k=startk, midk
227  l = mod(endk-(k-startk)-1,sym_kp*sym_nfp)+1
228  currerr = 0d0
229  DO i=1, sym_ir
230  DO j=1, sym_jz
231  h = sym_jz - (j-1)
232 
233  dbr = sym_br(i,j,k)+sym_br(i,h,l)
234  dbp = sym_bp(i,j,k)-sym_bp(i,h,l)
235  dbz = sym_bz(i,j,k)-sym_bz(i,h,l)
236  avgb = (sym_b(i,j,k)+sym_b(i,h,l))/2.0d0
237 
238  IF(avgb .NE. 0d0) THEN
239  currerr = currerr + dsqrt(dbr**2+dbp**2+dbz**2)/avgb
240  END IF
241  END DO
242  END DO
243 
244  avgerr = avgerr + currerr
245  IF(currerr .GT. maxerr) maxerr = currerr
246  END DO
247 
248  IF(midk-startk+1 .GT. 0) avgerr = avgerr/(midk-startk+1)
249 
250  END SUBROUTINE sym_checkstel
251 
252 !----------------------------------------------------------------------
253 !*******************************************************************************
254 !----------------------------------------------------------------------
255 
256  END MODULE sym_check
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17