V3FIT
spline_fields.f
1  SUBROUTINE spline_fields(arg1)
2  USE stel_kinds
3  USE stel_constants
4  USE bspline
5  USE spline3d_fit_coefs
6  USE mpi_params
7 #if defined(MPI_OPT)
8  USE mpi_inc
9 #endif
10  IMPLICIT NONE
11 
12  CHARACTER(LEN=60), INTENT(IN):: arg1
13  INTEGER:: i, j, k, istat, iphi0, itht0, is0, i1, j1,
14  & k1, isize
15  REAL(rprec), PARAMETER :: p5 = 0.5_dp
16  REAL(rprec):: maxabs, locabs, diff, c1
17  REAL(rprec), ALLOCATABLE, DIMENSION(:):: s, u, v
18  REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:):: r, z, phi,
19  & br, bz, bphi
20  REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:):: st, ut, vt
21  REAL(rprec), ALLOCATABLE, DIMENSION(:,:,:):: r_test,
22  & z_test, phi_test, br_test, bz_test, bphi_test
23 
24  OPEN(unit=4, file=trim(arg1), status="old", iostat=istat) ! File produced by SIESTA containing data
25  IF (istat .ne. 0) stop 'Cannot open specified file!'
26  READ(4,*) iphi0, is, itht0 ! IPHI0 = Number of toroidal nodes in file; IS = number o
27  iphi = iphi0 + 1; itht = itht0 + 1; is0 = is -1 ! Include U = 2*pi, V = 2*pi when splining. Thus, ITHT =
28 !
29  CALL alloc_splines ! Allocate quantities for splines which will be passed vi
30  ALLOCATE(s(is), u(itht), v(iphi), stat=istat) ! Allocate local quantities
31  IF (istat .NE. 0) stop 'Allocation error'
32  ALLOCATE(r(is,itht,iphi), z(is,itht,iphi),
33  & phi(is,itht,iphi), br(is,itht,iphi),
34  & bz(is,itht,iphi), bphi(is,itht,iphi), stat = istat)
35  IF (istat .NE. 0) stop 'Allocation error'
36  s = zero; v = zero; u = zero
37  r = zero; z = zero; phi = zero
38  br = zero; bz = zero; bphi = zero
39 !
40  DO k = 1, iphi ! Construct meshes. Assumed uniform.
41  v(k) = twopi*real(k-1,dp)/real(iphi0,dp) ! Extend last points in u, v to include 2pi:
42  ENDDO
43  DO i = 1, is
44  s(i) = real(i-1,dp)/real(is0,dp)
45  ENDDO
46  DO j = 1, itht
47  u(j) = twopi*real(j-1,dp)/real(itht0,dp)
48  ENDDO
49  if(myid .eq. 0) then
50  WRITE(*,*)
51  WRITE(*,*)' 3D fields being stored in file: ',
52  & trim(adjustl(arg1))
53  WRITE(*,'(a, i6, a8, i6, a8, i6)') ' Dimensions: Ns = ', is,
54  & ' Nu = ', itht0, ' Nv = ', iphi0
55  WRITE(*,*)' Reading ',trim(arg1),'..........'
56  endif !if(myid .eq. 0)
57 !OLD ORDER DO k = 1, iphi0 !EARLIER THAN 4-16-2013 - SPH
58 !OLD ORDER DO i = 1, is
59 !OLD ORDER DO j = 1, itht0
60  DO j = 1, itht0
61  DO k = 1, iphi0
62  DO i = 1, is
63  READ(4,*) r(i,j,k), z(i,j,k), phi(i,j,k),
64  & br(i,j,k), bz(i,j,k), bphi(i,j,k)
65  ENDDO
66  ENDDO
67  ENDDO
68  CLOSE(4)
69 
70  if(myid .eq. 0) then
71  WRITE(*,'(" .........completed")')
72  WRITE(*,*)
73  endif !if(myid .eq. 0) then
74 
75  CALL enforce_periodicity(is, itht, iphi, r, 0) ! Enforce periodicity of fields in U and V
76  CALL enforce_periodicity(is, itht, iphi, z, 0)
77  CALL enforce_periodicity(is, itht, iphi, phi, 1) ! Use 1 cause PHI is not periodic in V
78  CALL enforce_periodicity(is, itht, iphi, br, 0)
79  CALL enforce_periodicity(is, itht, iphi, bz, 0)
80  CALL enforce_periodicity(is, itht, iphi, bphi, 0)
81 
82  if(myid .eq. 0) then
83  WRITE(*,*) ' Creating spline knots'
84  endif !if(myid .eq. 0) then
85 
86  CALL dbsnak(is, s, ks, s_knots)
87  CALL dbsnak(iphi, v, kv, v_knots)
88  CALL dbsnak(itht, u, ku, u_knots)
89 
90  if(myid .eq. 0) then
91  WRITE(*,*) ' Computing spline coefficients'
92  endif !if(myid .eq. 0) then
93 
94  call dbs3in(is, s, itht, u, iphi, v, r, is, itht, ks, ku, kv,
95  > s_knots, u_knots, v_knots, r_coef)
96  call dbs3in(is, s, itht, u, iphi, v, z, is, itht, ks, ku, kv,
97  > s_knots, u_knots, v_knots, z_coef)
98  call dbs3in(is, s, itht, u, iphi, v, phi, is, itht, ks, ku, kv,
99  > s_knots, u_knots, v_knots, phi_coef)
100  call dbs3in(is, s, itht, u, iphi, v, br, is, itht, ks, ku, kv,
101  > s_knots, u_knots, v_knots, br_coef)
102  call dbs3in(is, s, itht, u, iphi, v, bz, is, itht, ks, ku, kv,
103  > s_knots, u_knots, v_knots, bz_coef)
104  call dbs3in(is, s, itht, u, iphi, v, bphi, is, itht, ks, ku, kv,
105  > s_knots, u_knots, v_knots, bphi_coef)
106 !
107 ! Test quality of spline at mid point of ALL 3D cells
108 ! by measuring the maximum relative difference between averaged
109 ! and splined value at center of 3D cell for ALL splined fields
110 !
111  IF (l_test_splines) THEN
112 
113  if(myid .eq. 0) then
114  WRITE(*,*) ' Testing quality of splines'
115  endif ! if(myid .eq. 0) then
116  isize = is0*itht0*iphi0
117  ALLOCATE(st(is0, itht0, iphi0), ut(is0, itht0, iphi0),
118  & vt(is0, itht0, iphi0),stat = istat)
119  IF (istat .NE. 0) stop 'Allocation error'
120  st = zero; ut = zero; vt = zero
121 
122  DO i = 1, is0
123  st(i,:,:) = (real(i,dp)-p5)/real(is,dp)
124  ENDDO
125  DO j = 1, itht0
126  ut(:,j,:) = twopi*(real(j,dp)-p5)/real(itht,dp)
127  ENDDO
128  DO k = 1, iphi0
129  vt(:,:,k) = twopi*(real(k,dp)-p5)/real(iphi,dp)
130  ENDDO
131 
132  ALLOCATE(br_test(is0, itht0, iphi0), stat=istat)
133  IF (istat .NE. 0) stop 'Allocation error'
134  br_test = zero
135  CALL dbs3vd(0, 0, 0, isize, st, ut, vt, ks, ku, kv,
136  & s_knots, u_knots, v_knots, is, itht, iphi,
137  & br_coef, br_test, is_uniformx, is_uniformy,
138  & is_uniformz)
139 
140  DO i = 1, is0
141  i1 = i + 1
142  DO j = 1, itht0
143  j1 = j + 1
144  DO k = 1, iphi0
145  k1 = k + 1
146  c1 = br(i,j,k) + br(i1,j,k)
147  & + br(i,j1,k) + br(i,j,k1)
148  & + br(i,j1,k1) + br(i1,j,k1)
149  & + br(i1,j1,k) + br(i1,j1,k1)
150  c1 = c1/8.d0 ! Average value at center of 3D cell
151  diff = abs(c1 - br_test(i,j,k))
152  IF (i==1 .AND. j == 1 .AND. k == 1) THEN
153  maxabs = diff
154  locabs = c1
155  ELSE
156  maxabs = max(maxabs, diff) ! Maximum (absolute) discrepancy
157  IF (maxabs == diff) locabs = c1 ! Store value at which max. discrepancy happens
158  ENDIF
159  ENDDO
160  ENDDO
161  ENDDO
162  if(myid .eq. 0) then
163  WRITE(*, 12) maxabs, locabs
164  endif ! if(myid .eq. 0) then
165 12 FORMAT(' - Max. absolute difference (Br):', f12.4,
166  & ' at value:', f12.4)
167 
168  DEALLOCATE(br_test)
169 
170  ALLOCATE(bz_test(is0, itht0, iphi0), stat = istat)
171  IF (istat .NE. 0) stop 'Allocation error'
172  bz_test = zero
173  CALL dbs3vd(0, 0, 0, isize, st, ut, vt, ks, ku, kv,
174  & s_knots, u_knots, v_knots, is, itht, iphi,
175  & bz_coef, bz_test, is_uniformx, is_uniformy,
176  & is_uniformz)
177 
178  DO i = 1, is0
179  i1 = i + 1
180  DO j = 1, itht0
181  j1 = j + 1
182  DO k = 1, iphi0
183  k1 = k + 1
184  c1 = bz(i,j,k) + bz(i1,j,k)
185  & + bz(i,j1,k) + bz(i,j,k1)
186  & + bz(i,j1,k1) + bz(i1,j,k1)
187  & + bz(i1,j1,k) + bz(i1,j1,k1)
188  c1 = c1/8.d0 ! Average value at center of 3D cell
189  diff = abs(c1 - bz_test(i,j,k))
190  IF (i==1 .AND. j == 1 .AND. k == 1) THEN
191  maxabs = diff
192  locabs = c1
193  ELSE
194  maxabs = max(maxabs, diff) ! Maximum (absolute) discrepancy
195  IF (maxabs == diff) locabs = c1 ! Store value at which max. discrepancy happens
196  ENDIF
197  ENDDO
198  ENDDO
199  ENDDO
200  if(myid .eq. 0) then
201  WRITE(*, 13) maxabs, locabs
202  endif ! if(myid .eq. 0) then
203 13 FORMAT(' - Max. absolute difference (Bz):', f12.4,
204  & ' at value:', f12.4)
205  DEALLOCATE(bz_test)
206 
207  ALLOCATE(bphi_test(is0, itht0, iphi0), stat = istat)
208  IF (istat .NE. 0) stop 'Allocation error'
209  bphi_test = zero
210  CALL dbs3vd(0, 0, 0, isize, st, ut, vt, ks, ku, kv,
211  & s_knots, u_knots, v_knots, is, itht, iphi,
212  & bphi_coef, bphi_test, is_uniformx, is_uniformy,
213  & is_uniformz)
214 
215  DO i = 1, is0
216  i1 = i + 1
217  DO j = 1, itht0
218  j1 = j + 1
219  DO k = 1, iphi0
220  k1 = k + 1
221  c1 = bphi(i,j,k) + bphi(i1,j,k)
222  & + bphi(i,j1,k) + bphi(i,j,k1)
223  & + bphi(i,j1,k1) + bphi(i1,j,k1)
224  & + bphi(i1,j1,k) + bphi(i1,j1,k1)
225  c1 = c1/8.d0
226  diff = abs(c1 - bphi_test(i,j,k))
227  IF (i==1 .AND. j == 1 .AND. k == 1) THEN
228  maxabs = diff
229  locabs = c1
230  ELSE
231  maxabs = max(maxabs, diff) ! Maximum (absolute) discrepancy
232  IF (maxabs == diff) locabs = c1 ! Store value at which max. discrepancy happens
233  ENDIF
234  ENDDO
235  ENDDO
236  ENDDO
237  if(myid .eq. 0) then
238  WRITE(*, 14) maxabs, locabs
239  endif ! if(myid .eq. 0) then
240 14 FORMAT(' - Max. absolute difference (Bphi):', f12.4,
241  & ' at value:', f12.4)
242  DEALLOCATE(bphi_test)
243 
244  ALLOCATE(r_test(is0, itht0, iphi0), stat = istat)
245  IF (istat .NE. 0) stop 'Allocation error'
246  r_test = zero
247  CALL dbs3vd(0, 0, 0, isize, st, ut, vt, ks, ku, kv,
248  & s_knots, u_knots, v_knots, is, itht, iphi,
249  & r_coef, r_test, is_uniformx, is_uniformy,
250  & is_uniformz)
251 
252  DO i = 1, is0
253  i1 = i + 1
254  DO j = 1, itht0
255  j1 = j + 1
256  DO k = 1, iphi0
257  k1 = k + 1
258  c1 = r(i,j,k) + r(i1,j,k)
259  & + r(i,j1,k) + r(i,j,k1)
260  & + r(i,j1,k1) + r(i1,j,k1)
261  & + r(i1,j1,k) + r(i1,j1,k1)
262  c1 = c1/8._dp
263  diff = abs(c1 - r_test(i,j,k))
264  IF (i==1 .AND. j == 1 .AND. k == 1) THEN
265  maxabs = diff
266  locabs = c1
267  ELSE
268  maxabs = max(maxabs, diff) ! Maximum (absolute) discrepancy
269  IF (maxabs == diff) locabs = c1 ! Store value at which max. discrepancy happens
270  ENDIF
271  ENDDO
272  ENDDO
273  ENDDO
274  if(myid .eq. 0) then
275  WRITE(*, 15) maxabs, locabs
276  endif ! if(myid .eq. 0) then
277 15 FORMAT(' - Max. absolute difference (R):', f12.4,
278  & ' at value:', f12.4)
279  DEALLOCATE(r_test)
280 
281  ALLOCATE(z_test(is0, itht0, iphi0), stat = istat)
282  IF (istat .NE. 0) stop 'Allocation error'
283  z_test = zero
284  CALL dbs3vd(0, 0, 0, isize, st, ut, vt, ks, ku, kv,
285  & s_knots, u_knots, v_knots, is, itht, iphi,
286  & z_coef, z_test, is_uniformx, is_uniformy,
287  & is_uniformz)
288 
289  DO i = 1, is0
290  i1 = i + 1
291  DO j = 1, itht0
292  j1 = j + 1
293  DO k = 1, iphi0
294  k1 = k + 1
295  c1 = z(i,j,k) + z(i1,j,k)
296  & + z(i,j1,k) + z(i,j,k1)
297  & + z(i,j1,k1) + z(i1,j,k1)
298  & + z(i1,j1,k) + z(i1,j1,k1)
299  c1 = c1/8._dp
300  diff = abs(c1 - z_test(i,j,k))
301  IF (i==1 .AND. j == 1 .AND. k == 1) THEN
302  maxabs = diff
303  locabs = c1
304  ELSE
305  maxabs = max(maxabs, diff) ! Maximum (absolute) discrepancy
306  IF (maxabs == diff) locabs = c1 ! Store value at which max. discrepancy happens
307  ENDIF
308  ENDDO
309  ENDDO
310  ENDDO
311  if(myid .eq. 0) then
312  WRITE(*, 16) maxabs, locabs
313  endif ! if(myid .eq. 0) then
314 16 FORMAT(' - Max. absolute difference (Z):', f12.4,
315  & ' at value:', f12.4)
316  DEALLOCATE(z_test)
317 
318  ALLOCATE(phi_test(is0, itht0, iphi0), stat = istat)
319  IF (istat .NE. 0) stop 'Allocation error'
320  phi_test = zero
321  CALL dbs3vd(0, 0, 0, isize, st, ut, vt, ks, ku, kv,
322  & s_knots, u_knots, v_knots, is, itht, iphi,
323  & phi_coef, phi_test, is_uniformx, is_uniformy,
324  & is_uniformz)
325 
326  DO i = 1, is0
327  i1 = i + 1
328  DO j = 1, itht0
329  j1 = j + 1
330  DO k = 1, iphi0
331  k1 = k + 1
332  c1 = phi(i,j,k) + phi(i1,j,k)
333  & + phi(i,j1,k) + phi(i,j,k1)
334  & + phi(i,j1,k1) + phi(i1,j,k1)
335  & + phi(i1,j1,k) + phi(i1,j1,k1)
336  c1 = c1/8._dp
337  diff = abs(c1 - phi_test(i,j,k))
338  IF (i==1 .AND. j == 1 .AND. k == 1) THEN
339  maxabs = diff
340  locabs = c1
341  ELSE
342  maxabs = max(maxabs, diff) ! Maximum (absolute) discrepancy
343  IF (maxabs == diff) locabs = c1 ! Store value at which max. discrepancy happens
344  ENDIF
345  ENDDO
346  ENDDO
347  ENDDO
348  if(myid .eq. 0) then
349  WRITE(*, 17) maxabs, locabs
350  endif ! if(myid .eq. 0) then
351 17 FORMAT(' - Max. absolute difference (PHI):', f12.4,
352  & ' at value:', f12.4)
353  DEALLOCATE(phi_test)
354 
355  DEALLOCATE(st, ut, vt)
356 
357  ENDIF
358  DEALLOCATE(s, u, v, r, z, phi, br, bz, bphi)
359 
360  END SUBROUTINE spline_fields
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17