1 SUBROUTINE scalfor_par(gcx, axm, bxm, axd, bxd, cx, iflag)
5 USE realspace,
ONLY: wint, ru0
6 USE parallel_include_module
7 USE parallel_vmec_module,
ONLY: padsides1x
8 USE xstuff,
ONLY: pxc, pgc
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) ::
18 REAL(dp),
DIMENSION(ns),
INTENT(IN) :: cx
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
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
34 IF (.NOT.lactive)
RETURN
37 CALL padsides1x(axm(:,i))
38 CALL padsides1x(bxm(:,i))
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
56 nsmax = min(trglob, jmax)
58 nsmin = max(jmin2(m), tlglob)
61 ax(n,m,tlglob:trglob) = 0
62 bx(n,m,tlglob:trglob) = 0
63 dx(n,m,tlglob:trglob) = 0
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)
71 IF (m .eq. 1 .and. nsmin .eq. 2)
THEN
72 dx(n,m,2) = dx(n,m,2) + bx(n,m,2)
77 IF (jmax .GE. ns)
THEN
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)
91 mult_fac = min(fac, fac*hs*15)
92 IF (iflag .eq. 1)
THEN
96 dx(0,0,ns) = dx(0,0,ns)*(1 - mult_fac)/(1 + edge_pedestal)
122 IF (iresidue .GE. 0 .AND. iresidue .LT. 3)
THEN
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)
133 DEALLOCATE (ax, bx, dx)
135 CALL second0(scalfortoff)
138 END SUBROUTINE scalfor_par
140 SUBROUTINE bst_parallel_tridiag_solver(a, d, b, c, jmin,
141 1 jmax, mnd1, ns, nrhs)
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
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
163 REAL(dp),
PARAMETER :: zero = 0, one = 1
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
178 IF (jmax .GT. ns)
THEN
179 stop
'jmax>ns in tridslv_par'
184 IF (in1 .ge. in0)
THEN
188 c(mn, in0:in1, 1:nrhs) = 0
193 ALLOCATE(tmp(blklength,blklength))
196 init_time = init_time + (t2-t1)
199 DO irow = tlglob, trglob
202 IF (irow .EQ. ns .AND. jmax .LT. ns)
THEN
205 CALL setmatrixrowcoll_bst(irow,b(:,irow))
208 IF (irow .EQ. ns .AND. jmax .LT. ns)
THEN
211 CALL setmatrixrowcold_bst(irow,d(:,irow))
214 CALL setmatrixrowcolu_bst(irow,a(:,irow))
217 setup_time = setup_time + (t2-t1)
220 CALL forwardsolve_bst
222 forwardsolve_time = forwardsolve_time + (t2-t1)
224 ALLOCATE(tmpv(0:mnd1))
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)
236 CALL backwardsolve_bst
239 DO irow = tlglob, trglob
240 CALL getsolutionvector_bst(irow, tmpv)
241 c(:,irow,jrhs)=tmpv(0:mnd1)
245 DEALLOCATE(tmp, tmpv)
247 backwardsolve_time = backwardsolve_time + (t2-t1)
249 END SUBROUTINE bst_parallel_tridiag_solver
251 SUBROUTINE scalfor(gcx, axm, bxm, axd, bxd, cx, iflag)
254 USE vmec_dim,
ONLY: ns
255 USE realspace,
ONLY: wint, ru0
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) ::
266 REAL(dp),
DIMENSION(ns),
INTENT(in) :: cx
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
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
284 IF (ivac .lt. 1)
THEN
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)
304 dx(2,n,m) = dx(2,n,m) + bx(2,n,m)
314 IF (jmax .ge. ns)
THEN
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)
328 mult_fac = min(fac, fac*hs*15)
329 IF (iflag .eq. 1)
THEN
333 dx(ns,0,0) = dx(ns,0,0)*(1-mult_fac)/(1+edge_pedestal)
360 IF (iresidue.GE.0 .and. iresidue.LT.3)
THEN
366 CALL serial_tridslv (ax, dx, bx, gcx, jmin4, jmax, mnsize - 1,
368 DEALLOCATE (ax, bx, dx)
370 END SUBROUTINE scalfor
372 SUBROUTINE serial_tridslv(a, d, b, c, jmin, jmax, mnd1, ns, nrhs)
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
385 REAL(dp),
PARAMETER :: zero = 0, one = 1
389 INTEGER :: mn, in, i0, in1, jrhs
390 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:) :: alf
391 REAL(dp),
DIMENSION(0:mnd1) :: psi0
399 IF (jmax .gt. ns)
THEN
400 stop
'jmax>ns in tridslv'
403 ALLOCATE (alf(ns,0:mnd1), stat = in)
405 stop
'Allocation error in tridslv'
415 IF (in1 .ge. in)
THEN
417 c(in:in1, mn, 1:nrhs) = 0
426 IF (any(psi0 .eq. zero))
THEN
427 stop
'psi0 == 0 error in tridslv'
431 c(in,:,jrhs) = c(in,:,jrhs)*psi0(:)
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'
442 c(i0,:,jrhs) = (c(i0,:,jrhs) - b(i0,:)*c(i0-1,:,jrhs))*psi0
446 DO i0 = jmax - 1, in, -1
448 c(i0,:,jrhs) = c(i0,:,jrhs) - alf(i0,:)*c(i0+1,:,jrhs)
454 END SUBROUTINE serial_tridslv