V3FIT
scalfor.f
1  SUBROUTINE scalfor_par(gcx, axm, bxm, axd, bxd, cx, iflag)
2  USE vmec_main
3  USE vmec_params
4  USE vmec_dim, ONLY: ns
5  USE realspace, ONLY: wint, ru0
6  USE parallel_include_module
7  USE parallel_vmec_module, ONLY: padsides1x
8  USE xstuff, ONLY: pxc, pgc
9  IMPLICIT NONE
10 C-----------------------------------------------
11 C Dummy Arguments
12 C-----------------------------------------------
13  INTEGER, INTENT(IN) :: iflag
14  REAL(dp), DIMENSION(0:ntor,0:mpol1,ns,ntmax),
15  1 INTENT(INOUT) :: gcx
16  REAL(dp), DIMENSION(ns+1,2), INTENT(INOUT) ::
17  1 axm, bxm, axd, bxd
18  REAL(dp), DIMENSION(ns), INTENT(IN) :: cx
19 C-----------------------------------------------
20 C L o c a l V a r i a b l e s
21 C-----------------------------------------------
22  REAL(dp), PARAMETER :: ftol_edge = 1.e-9_dp, c1p5 = 1.5_dp,
23  1 fac = 0.25_dp, edge_pedestal = 0.05_dp
24  INTEGER :: m , mp, n, js, jmax, jmin4(0:mnsize-1)
25  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: ax, bx, dx
26  REAL(dp) :: mult_fac
27  INTEGER :: nsmin, nsmax, i, j, k, l
28  INTEGER, ALLOCATABLE, DIMENSION(:) :: counts, disps
29  REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: send_buf2
30  INTEGER :: MPI_STAT(MPI_STATUS_SIZE)
31  REAL(dp) :: tridslvton, tridslvtoff
32  REAL(dp) :: scalforton, scalfortoff
33 C-----------------------------------------------
34  IF (.NOT.lactive) RETURN
35 
36  DO i = 1, 2
37  CALL padsides1x(axm(:,i))
38  CALL padsides1x(bxm(:,i))
39  END DO
40 
41  ALLOCATE (ax(0:ntor,0:mpol1,ns), bx(0:ntor,0:mpol1,ns),
42  1 dx(0:ntor,0:mpol1,ns))
43  ax(:,:,1) = 0; bx(:,:,1) = 0; dx(:,:,1) = 0
44  ax(:,:,ns) = 0; bx(:,:,ns) = 0; dx(:,:,ns) = 0
45 
46  jmax = ns
47  IF (ivac .lt. 1) THEN
48  jmax = ns1
49  END IF
50 
51 ! FOR SOME 3D PLASMAS, THIS SOMETIME HELPS (CHOOSE mult_fac =1 otherwise)
52 ! TO AVOID JACOBIAN RESETS BY GIVING A SMOOTH TRANSITION FROM FIXED TO FREE ITERATIONS
53 ! mult_fac = 1._dp/(1._dp + 10*(fsqr+fsqz))
54 ! gcx(ns,:,:,:) = mult_fac*gcx(ns,:,:,:)
55 
56  nsmax = min(trglob, jmax)
57  DO m = 0, mpol1
58  nsmin = max(jmin2(m), tlglob)
59  mp = mod(m,2) + 1
60  DO n = 0, ntor
61  ax(n,m,tlglob:trglob) = 0
62  bx(n,m,tlglob:trglob) = 0
63  dx(n,m,tlglob:trglob) = 0
64  DO js = nsmin, nsmax
65  ax(n,m,js) = -(axm(js+1,mp) + bxm(js+1,mp)*m**2)
66  bx(n,m,js) = -(axm(js,mp) + bxm(js,mp)*m**2)
67  dx(n,m,js) = -(axd(js,mp) + bxd(js,mp)*m**2
68  & + cx(js)*(n*nfp)**2)
69  END DO
70 
71  IF (m .eq. 1 .and. nsmin .eq. 2) THEN
72  dx(n,m,2) = dx(n,m,2) + bx(n,m,2)
73  END IF
74  END DO
75  END DO
76 
77  IF (jmax .GE. ns) THEN
78 !
79 ! SMALL EDGE PEDESTAL NEEDED TO IMPROVE CONVERGENCE
80 ! IN PARTICULAR, NEEDED TO ACCOUNT FOR POTENTIAL ZERO
81 ! EIGENVALUE DUE TO NEUMANN (GRADIENT) CONDITION AT EDGE
82 !
83  dx(:,0:1,ns) = (1 + edge_pedestal) *dx(:,0:1,ns)
84  dx(:,2:mpol1,ns) = (1 + 2*edge_pedestal)*dx(:,2:mpol1,ns)
85 !
86 ! STABILIZATION ALGORITHM FOR ZC_00(NS)
87 ! FOR UNSTABLE CASE, HAVE TO FLIP SIGN OF -FAC -> +FAC FOR CONVERGENCE
88 ! COEFFICIENT OF < Ru (R Pvac)> ~ -fac*(z-zeq) WHERE fac (EIGENVALUE, OR
89 ! FIELD INDEX) DEPENDS ON THE EQUILIBRIUM MAGNETIC FIELD AND CURRENT,
90 ! AND zeq IS THE EQUILIBRIUM EDGE VALUE OF Z00
91  mult_fac = min(fac, fac*hs*15)
92  IF (iflag .eq. 1) THEN
93 !
94 ! METHOD 1: SUBTRACT (INSTABILITY) Pedge ~ fac*z/hs FROM PRECONDITIONER AT EDGE
95 !
96  dx(0,0,ns) = dx(0,0,ns)*(1 - mult_fac)/(1 + edge_pedestal)
97  END IF
98 
99  ENDIF
100 !
101 ! ACCELERATE (IMPROVE) CONVERGENCE OF FREE BOUNDARY. THIS WAS ADDED
102 ! TO DEAL WITH CASES WHICH MIGHT OTHERWISE DIVERGE. BY DECREASING THE
103 ! FSQ TOLERANCE LEVEL WHERE THIS KICKS IN (FTOL_EDGE), THE USER CAN
104 ! TURN-OFF THIS FEATURE
105 !
106 ! DIAGONALIZE (DX DOMINANT) AND REDUCE FORCE (DX ENHANCED) AT EDGE
107 ! TO IMPROVE CONVERGENCE FOR N != 0 TERMS
108 !
109 
110 ! ledge = .false.
111 ! IF ((fsqr+fsqz) .lt. ftol_edge) ledge = .true.
112 ! IF ((iter2-iter1).lt.400 .or. ivac.lt.1) ledge = .false.
113 
114 ! IF (ledge) THEN
115 ! dx(ns,1:,1:) = 3*dx(ns,1:,1:)
116 ! END IF
117 
118 ! FOR DATA MATCHING MODE (0 <= IRESIDUE < 3),
119 ! MAGNETIC AXIS IS FIXED SO JMIN3(0) => 2 FOR M=0,N=0
120 
121  jmin4 = jmin3
122  IF (iresidue .GE. 0 .AND. iresidue .LT. 3) THEN
123  jmin4(0) = 2
124  END IF
125 
126 ! Padsides moved to SUBROUTINE residue AFTER this completes
127  CALL second0(tridslvton)
128  CALL bst_parallel_tridiag_solver(ax,dx,bx,gcx,jmin4,jmax,
129  1 mnsize - 1,ns,ntmax)
130  CALL second0(tridslvtoff)
131  tridslv_time = tridslv_time + (tridslvtoff-tridslvton)
132 
133  DEALLOCATE (ax, bx, dx)
134 
135  CALL second0(scalfortoff)
136 ! scalfor_time = scalfor_time + (scalfortoff-scalforton)
137 
138  END SUBROUTINE scalfor_par
139 
140  SUBROUTINE bst_parallel_tridiag_solver(a, d, b, c, jmin,
141  1 jmax, mnd1, ns, nrhs)
142  USE stel_kinds
143  USE parallel_include_module
144  USE blocktridiagonalsolver_bst, ONLY: setmatrixrowcoll_bst
145  USE blocktridiagonalsolver_bst, ONLY: setmatrixrowcold_bst
146  USE blocktridiagonalsolver_bst, ONLY: setmatrixrowcolu_bst
147  USE blocktridiagonalsolver_bst, ONLY: forwardsolve_bst
148  USE blocktridiagonalsolver_bst, ONLY: setmatrixrhs_bst
149  USE blocktridiagonalsolver_bst, ONLY: backwardsolve_bst
150  USE blocktridiagonalsolver_bst, ONLY: getsolutionvector_bst
151 
152  IMPLICIT NONE
153 C-----------------------------------------------
154 C D u m m y A r g u m e n t s
155 C-----------------------------------------------
156  INTEGER, INTENT(IN) :: jmax, mnd1, ns, nrhs
157  INTEGER, DIMENSION(0:mnd1), INTENT(IN) :: jmin
158  REAL(dp), DIMENSION(0:mnd1,ns) :: a, d, b
159  REAL(dp), DIMENSION(0:mnd1,ns,nrhs), INTENT(INOUT) :: c
160 C-----------------------------------------------
161 C L o c a l P a r a m e t e r s
162 C-----------------------------------------------
163  REAL(dp), PARAMETER :: zero = 0, one = 1
164 C-----------------------------------------------
165 C L o c a l V a r i a b l e s
166 C-----------------------------------------------
167  INTEGER :: mn, in0, in1, jrhs
168  INTEGER :: irow, icol, blklength, i, j
169  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: tmp
170  REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmpv
171  REAL(dp), DIMENSION(0:mnd1) :: psi0
172  REAL(dp) :: t1, t2
173 C-----------------------------------------------
174 ! SOLVES B(I)*X(I-1)+D(I)*X(I)+A(I)*X(I+1)=C(I), I=IN,JMAX
175 ! AND RETURNS ANSWER IN C(I)
176 
177  CALL second0(t1)
178  IF (jmax .GT. ns) THEN
179  stop 'jmax>ns in tridslv_par'
180  END IF
181  in0 = minval(jmin)
182  DO mn = 0, mnd1
183  in1 = jmin(mn)-1
184  IF (in1 .ge. in0) THEN
185  d(mn, in0:in1) = 1
186  b(mn, in0:in1) = 0
187  a(mn, in0:in1) = 0
188  c(mn, in0:in1, 1:nrhs) = 0
189  END IF
190  END DO
191 
192  blklength=mnd1+1
193  ALLOCATE(tmp(blklength,blklength))
194  tmp=zero
195  CALL second0(t2)
196  init_time = init_time + (t2-t1)
197 
198  CALL second0(t1)
199  DO irow = tlglob, trglob
200 
201  ! Set up L
202  IF (irow .EQ. ns .AND. jmax .LT. ns) THEN
203  b(:,irow) = 0
204  END IF
205  CALL setmatrixrowcoll_bst(irow,b(:,irow))
206 
207  ! Set up D
208  IF (irow .EQ. ns .AND. jmax .LT. ns) THEN
209  d(:,irow) = 1
210  END IF
211  CALL setmatrixrowcold_bst(irow,d(:,irow))
212 
213  ! Set up U
214  CALL setmatrixrowcolu_bst(irow,a(:,irow))
215  END DO
216  CALL second0(t2)
217  setup_time = setup_time + (t2-t1)
218 
219  CALL second0(t1)
220  CALL forwardsolve_bst
221  CALL second0(t2)
222  forwardsolve_time = forwardsolve_time + (t2-t1)
223 
224  ALLOCATE(tmpv(0:mnd1))
225  CALL second0(t1)
226  DO jrhs = 1, nrhs
227 
228  ! Set RHS
229  DO irow = tlglob, trglob
230  tmpv(0:mnd1)=c(:,irow,jrhs)
231  IF (irow.EQ.ns.AND.jmax.LT.ns) tmpv(0:mnd1)=0
232  CALL setmatrixrhs_bst(irow,tmpv)
233  END DO
234 
235  ! Backward solve
236  CALL backwardsolve_bst
237 
238  ! Get solution vector
239  DO irow = tlglob, trglob
240  CALL getsolutionvector_bst(irow, tmpv)
241  c(:,irow,jrhs)=tmpv(0:mnd1)
242  END DO
243 
244  END DO
245  DEALLOCATE(tmp, tmpv)
246  CALL second0(t2)
247  backwardsolve_time = backwardsolve_time + (t2-t1)
248 
249  END SUBROUTINE bst_parallel_tridiag_solver
250 
251  SUBROUTINE scalfor(gcx, axm, bxm, axd, bxd, cx, iflag)
252  USE vmec_main
253  USE vmec_params
254  USE vmec_dim, ONLY: ns
255  USE realspace, ONLY: wint, ru0
256 
257  IMPLICIT NONE
258 C-----------------------------------------------
259 C Dummy Arguments
260 C-----------------------------------------------
261  INTEGER, INTENT(in) :: iflag
262  REAL(dp), DIMENSION(ns,0:ntor,0:mpol1,ntmax),
263  1 INTENT(inout) :: gcx
264  REAL(dp), DIMENSION(ns+1,2), INTENT(in) ::
265  1 axm, bxm, axd, bxd
266  REAL(dp), DIMENSION(ns), INTENT(in) :: cx
267 C-----------------------------------------------
268 C L o c a l V a r i a b l e s
269 C-----------------------------------------------
270  REAL(dp), PARAMETER :: ftol_edge = 1.e-9_dp, c1p5 = 1.5_dp,
271  1 fac = 0.25_dp, edge_pedestal = 0.05_dp
272  INTEGER :: m , mp, n, js, jmax, jmin4(0:mnsize-1)
273  REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: ax, bx, dx
274  REAL(dp) :: mult_fac
275 
276 C-----------------------------------------------
277 
278  ALLOCATE (ax(ns,0:ntor,0:mpol1), bx(ns,0:ntor,0:mpol1),
279  & dx(ns,0:ntor,0:mpol1))
280  ax(1,:,:) = 0; bx(1,:,:) = 0; dx(1,:,:) = 0
281  ax(ns,:,:) = 0; bx(ns,:,:) = 0; dx(ns,:,:) = 0
282 
283  jmax = ns
284  IF (ivac .lt. 1) THEN
285  jmax = ns1
286  END IF
287 
288 ! FOR SOME 3D PLASMAS, THIS SOMETIME HELPS (CHOOSE mult_fac =1 otherwise)
289 ! TO AVOID JACOBIAN RESETS BY GIVING A SMOOTH TRANSITION FROM FIXED TO FREE ITERATIONS
290 ! mult_fac = 1._dp/(1._dp + 10*(fsqr+fsqz))
291 ! gcx(ns,:,:,:) = mult_fac*gcx(ns,:,:,:)
292 
293  DO m = 0, mpol1
294  mp = mod(m,2) + 1
295  DO n = 0, ntor
296  DO js = jmin2(m), jmax
297  ax(js,n,m) = -(axm(js+1,mp) + bxm(js+1,mp)*m**2)
298  bx(js,n,m) = -(axm(js,mp) + bxm(js,mp)*m**2)
299  dx(js,n,m) = -(axd(js,mp) + bxd(js,mp)*m**2
300  1 + cx(js)*(n*nfp)**2)
301  END DO
302 
303  IF (m .eq. 1) THEN
304  dx(2,n,m) = dx(2,n,m) + bx(2,n,m)
305 !OFF 050311 DO js = jmin2(m), jmax
306 !OFF 050311 ax(js,n,m) = c1p5*ax(js,n,m)
307 !OFF 050311 bx(js,n,m) = c1p5*bx(js,n,m)
308 !OFF 050311 dx(js,n,m) = c1p5*dx(js,n,m)
309 !OFF 050311 END DO
310  END IF
311  END DO
312  END DO
313 
314  IF (jmax .ge. ns) THEN
315 !
316 ! SMALL EDGE PEDESTAL NEEDED TO IMPROVE CONVERGENCE
317 ! IN PARTICULAR, NEEDED TO ACCOUNT FOR POTENTIAL ZERO
318 ! EIGENVALUE DUE TO NEUMANN (GRADIENT) CONDITION AT EDGE
319 !
320  dx(ns,:,0:1) = (1+edge_pedestal) *dx(ns,:,0:1)
321  dx(ns,:,2:mpol1) = (1+2*edge_pedestal)*dx(ns,:,2:mpol1)
322 !
323 ! STABILIZATION ALGORITHM FOR ZC_00(NS)
324 ! FOR UNSTABLE CASE, HAVE TO FLIP SIGN OF -FAC -> +FAC FOR CONVERGENCE
325 ! COEFFICIENT OF < Ru (R Pvac)> ~ -fac*(z-zeq) WHERE fac (EIGENVALUE, OR
326 ! FIELD INDEX) DEPENDS ON THE EQUILIBRIUM MAGNETIC FIELD AND CURRENT,
327 ! AND zeq IS THE EQUILIBRIUM EDGE VALUE OF Z00
328  mult_fac = min(fac, fac*hs*15)
329  IF (iflag .eq. 1) THEN
330 !
331 ! METHOD 1: SUBTRACT (INSTABILITY) Pedge ~ fac*z/hs FROM PRECONDITIONER AT EDGE
332 !
333  dx(ns,0,0) = dx(ns,0,0)*(1-mult_fac)/(1+edge_pedestal)
334  END IF
335  END IF
336 
337 
338 !
339 ! ACCELERATE (IMPROVE) CONVERGENCE OF FREE BOUNDARY. THIS WAS ADDED
340 ! TO DEAL WITH CASES WHICH MIGHT OTHERWISE DIVERGE. BY DECREASING THE
341 ! FSQ TOLERANCE LEVEL WHERE THIS KICKS IN (FTOL_EDGE), THE USER CAN
342 ! TURN-OFF THIS FEATURE
343 !
344 ! DIAGONALIZE (DX DOMINANT) AND REDUCE FORCE (DX ENHANCED) AT EDGE
345 ! TO IMPROVE CONVERGENCE FOR N != 0 TERMS
346 !
347 
348 ! ledge = .false.
349 ! IF ((fsqr+fsqz) .lt. ftol_edge) ledge = .true.
350 ! IF ((iter2-iter1).lt.400 .or. ivac.lt.1) ledge = .false.
351 
352 ! IF (ledge) THEN
353 ! dx(ns,1:,1:) = 3*dx(ns,1:,1:)
354 ! END IF
355 
356 ! FOR DATA MATCHING MODE (0 <= IRESIDUE < 3),
357 ! MAGNETIC AXIS IS FIXED SO JMIN3(0) => 2 FOR M=0,N=0
358 
359  jmin4 = jmin3
360  IF (iresidue.GE.0 .and. iresidue.LT.3) THEN
361  jmin4(0) = 2
362  END IF
363 
364 ! SOLVES BX(I)*X(I-1)+DX(I)*X(I)+AX(I)*X(I+1)=GCX(I), I=JMIN4,JMAX
365 ! AND RETURNS ANSWER IN GCX(I)
366  CALL serial_tridslv (ax, dx, bx, gcx, jmin4, jmax, mnsize - 1,
367  & ns, ntmax)
368  DEALLOCATE (ax, bx, dx)
369 
370  END SUBROUTINE scalfor
371 
372  SUBROUTINE serial_tridslv(a, d, b, c, jmin, jmax, mnd1, ns, nrhs)
373  USE stel_kinds
374  IMPLICIT NONE
375 C-----------------------------------------------
376 C D u m m y A r g u m e n t s
377 C-----------------------------------------------
378  INTEGER, INTENT(in) :: jmax, mnd1, ns, nrhs
379  INTEGER, DIMENSION(0:mnd1), INTENT(in) :: jmin
380  REAL(dp), DIMENSION(ns,0:mnd1) :: a, d, b
381  REAL(dp), DIMENSION(ns,0:mnd1, nrhs), INTENT(inout) :: c
382 C-----------------------------------------------
383 C L o c a l P a r a m e t e r s
384 C-----------------------------------------------
385  REAL(dp), PARAMETER :: zero = 0, one = 1
386 C-----------------------------------------------
387 C L o c a l V a r i a b l e s
388 C-----------------------------------------------
389  INTEGER :: mn, in, i0, in1, jrhs
390  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: alf
391  REAL(dp), DIMENSION(0:mnd1) :: psi0
392 C-----------------------------------------------
393 !
394 ! SOLVES B(I)*X(I-1)+D(I)*X(I)+A(I)*X(I+1)=C(I), I=IN,JMAX
395 ! AND RETURNS ANSWER IN C(I)
396 ! ADDED VECTORIZATION ON FOURIER MODE ARGUMENT (01-2000)
397 ! AND NEW ARGUMENT (NRHS) TO DO MULTIPLE RIGHT SIDES SIMULTANEOUSLY
398 !
399  IF (jmax .gt. ns) THEN
400  stop 'jmax>ns in tridslv'
401  END IF
402 
403  ALLOCATE (alf(ns,0:mnd1), stat = in)
404  IF (in .ne. 0) THEN
405  stop 'Allocation error in tridslv'
406  END IF
407 
408  in = minval(jmin)
409 !
410 ! FILL IN MN BELOW MAX(JMIN) WITH DUMMY VALUES
411 ! TO ALLOW VECTORIZATION ON MN INDEX
412 !
413  DO mn = 0, mnd1
414  in1 = jmin(mn)-1
415  IF (in1 .ge. in) THEN
416  d(in:in1, mn) = 1
417  c(in:in1, mn, 1:nrhs) = 0
418  b(in:in1, mn) = 0
419  a(in:in1, mn) = 0
420  END IF
421  END DO
422 
423  in1 = in + 1
424 
425  psi0(:)= d(in,:)
426  IF (any(psi0 .eq. zero)) THEN
427  stop 'psi0 == 0 error in tridslv'
428  END IF
429  psi0 = one/psi0
430  DO jrhs = 1, nrhs
431  c(in,:,jrhs) = c(in,:,jrhs)*psi0(:)
432  END DO
433 
434  DO i0 = in1,jmax
435  alf(i0-1,:) = a(i0-1,:)*psi0(:)
436  psi0 = d(i0,:) - b(i0,:)*alf(i0-1,:)
437  IF (any(abs(psi0) .le. 1.e-8_dp*abs(d(i0,:)))) THEN
438  stop 'psi0/d(i0) < 1.E-8: possible singularity in tridslv'
439  END IF
440  psi0 = one/psi0
441  DO jrhs = 1, nrhs
442  c(i0,:,jrhs) = (c(i0,:,jrhs) - b(i0,:)*c(i0-1,:,jrhs))*psi0
443  END DO
444  END DO
445 
446  DO i0 = jmax - 1, in, -1
447  DO jrhs = 1,nrhs
448  c(i0,:,jrhs) = c(i0,:,jrhs) - alf(i0,:)*c(i0+1,:,jrhs)
449  END DO
450  END DO
451 
452  DEALLOCATE (alf)
453 
454  END SUBROUTINE serial_tridslv
455