1 SUBROUTINE residue_par (gcr, gcz, gcl)
2 USE vmec_main, p5 => cp5
3 USE vmec_params,
ONLY: rss, zcs, rsc, zcc, &
4 meven, modd, ntmax, signgs
5 USE realspace,
ONLY: phip
8 USE parallel_include_module
9 USE parallel_vmec_module,
ONLY: tlglob_arr, trglob_arr, &
17 REAL(dp),
DIMENSION(0:ntor,0:mpol1,ns,ntmax),
INTENT(INOUT) ::
22 INTEGER,
PARAMETER :: n0=0, m0=0, m1=1
23 INTEGER,
PARAMETER :: n3d=0, nasym=1
27 INTEGER :: nsfix, jedge, delIter
28 REAL(dp) :: r1, fac, tmp, tmp2(ns), ftotal
30 INTEGER :: i, j, k, l, m, blksize, left, right
31 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: counts, disps
32 INTEGER :: MPI_STAT(MPI_STATUS_SIZE)
33 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:,:,:) :: send_buf
34 REAL(dp),
ALLOCATABLE,
DIMENSION(:) :: recv_buf
35 REAL(dp) :: tredon, tredoff
41 IF (lfreeb .AND. lrfp)
THEN
43 IF (ictrl_prec2d .EQ. 0)
THEN
46 gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
64 CALL constrain_m1_par(gcr(:,m1,:,rss), gcz(:,m1,:,zcs))
67 CALL constrain_m1_par(gcr(:,m1,:,rsc), gcz(:,m1,:,zcc))
70 IF (lfreeb .AND. lrfp)
THEN
72 IF (ictrl_prec2d .EQ. 0)
THEN
75 gcr(0,m0,ns,:) = fac*gcr(0,m0,ns,:)
76 gcz(0,m0,ns,:) = fac*gcz(0,m0,ns,:)
81 IF (ictrl_prec2d .GE. 2 .OR. ictrl_prec2d .EQ. -1)
RETURN
86 r1 = one/(2*r0scale)**2
90 IF (deliter .lt. 50 .and.
91 (fsqr + fsqz) .LT. 1.e-6_dp)
THEN
95 CALL getfsq_par (gcr, gcz, fsqr, fsqz, r1*fnorm, jedge)
98 tmp = sum(gcl(:,:,tlglob:trglob,:)*gcl(:,:,tlglob:trglob,:))
99 CALL mpi_allreduce(tmp,ftotal,1,mpi_real8,mpi_sum,ns_comm,mpi_err)
100 CALL second0(tredoff)
101 allreduce_time = allreduce_time + (tredoff - tredon)
103 IF(rank .EQ. nranks-1)
THEN
104 fedge = r1*fnorm*sum(gcr(:,:,ns,:)**2 + gcz(:,:,ns,:)**2)
109 IF (ictrl_prec2d .EQ. 1)
THEN
111 IF (l_colscale .AND. lactive)
THEN
112 CALL saxlastntype(pgc, pcol_scale, pgc)
115 lresiduecall = .true.
116 CALL block_precond_par(pgc)
117 lresiduecall = .false.
119 IF (.NOT.lfreeb .AND. any(gcr(:,:,ns,:) .NE. zero))
THEN
120 stop
'gcr(ns) != 0 for fixed boundary in residue'
122 IF (.NOT.lfreeb .AND. any(gcz(:,:,ns,:) .NE. zero))
THEN
123 stop
'gcz(ns) != 0 for fixed boundary in residue'
125 IF (any(gcl(1:,m0,:,zsc) .NE. zero))
THEN
126 stop
'gcl(m=0,n>0,sc) != 0 in residue'
128 IF (lthreed .AND. any(gcl(n0,:,:,zcs) .NE. zero))
THEN
129 stop
'gcl(n=0,m,cs) != 0 in residue'
140 CALL scale_m1_par(gcr(:,m1,:,rss), gcz(:,m1,:,zcs))
143 CALL scale_m1_par(gcr(:,m1,:,rsc), gcz(:,m1,:,zcc))
147 CALL scalfor_par (gcr, arm, brm, ard, brd, crd, jedge)
149 CALL scalfor_par (gcz, azm, bzm, azd, bzd, crd, jedge)
151 CALL getfsq_par (gcr, gcz, fsqr1, fsqz1, fnorm1, m1)
153 DO l = tlglob, trglob
154 gcl(:,:,l,:) = pfaclam(:,:,l,:)*gcl(:,:,l,:)
155 tmp2(l) = sum(gcl(:,:,l,:)**2)
157 CALL gather1xarray(tmp2)
158 ftotal = sum(tmp2(2:ns))
165 CALL second0 (tresoff)
166 residue_time = residue_time + (tresoff-treson)
168 END SUBROUTINE residue_par
170 SUBROUTINE constrain_m1_par(gcr, gcz)
172 USE parallel_include_module
173 USE precon2d,
ONLY: ictrl_prec2d
178 REAL(dp),
DIMENSION(0:ntor,ns),
INTENT(INOUT) :: gcr, gcz
182 REAL(dp),
PARAMETER :: FThreshold = 1.e-6_dp
183 REAL(dp),
ALLOCATABLE,
DIMENSION(:,:) :: temp
190 ALLOCATE(temp(0:ntor,ns))
192 temp(:,tlglob:trglob) = gcr(:,tlglob:trglob)
193 gcr(:,tlglob:trglob) = osqrt2*(gcr(:,tlglob:trglob) +
194 gcz(:,tlglob:trglob))
195 gcz(:,tlglob:trglob) = osqrt2*(temp(:,tlglob:trglob) -
196 gcz(:,tlglob:trglob))
200 IF (fsqz .LT. fthreshold .OR.
202 & ictrl_prec2d .NE. 0)
THEN
203 gcz(:,tlglob:trglob) = 0
207 END SUBROUTINE constrain_m1_par
209 SUBROUTINE scale_m1_par(gcr, gcz)
211 USE parallel_include_module
216 REAL(dp),
DIMENSION(0:ntor,ns),
INTENT(inout) :: gcr, gcz
220 INTEGER,
PARAMETER :: nodd=2
224 IF (.not.lconm1)
RETURN
226 fac(tlglob:trglob) = (ard(tlglob:trglob,nodd) +
227 brd(tlglob:trglob,nodd))
228 / (ard(tlglob:trglob,nodd) +
229 brd(tlglob:trglob,nodd) +
230 azd(tlglob:trglob,nodd) +
231 bzd(tlglob:trglob,nodd))
233 gcr(n,tlglob:trglob) = fac(tlglob:trglob)*gcr(n,tlglob:trglob)
236 fac(tlglob:trglob) = (azd(tlglob:trglob,nodd) +
237 bzd(tlglob:trglob,nodd))
238 / (ard(tlglob:trglob,nodd) +
239 brd(tlglob:trglob,nodd) +
240 azd(tlglob:trglob,nodd) +
241 bzd(tlglob:trglob,nodd))
243 gcz(n,tlglob:trglob) = fac(tlglob:trglob)*gcz(n,tlglob:trglob)
246 END SUBROUTINE scale_m1_par
248 SUBROUTINE residue (gcr, gcz, gcl)
249 USE vmec_main, p5 => cp5
250 USE vmec_params,
ONLY: rss, zcs, rsc, zcc, &
251 meven, modd, ntmax, signgs
252 USE realspace,
ONLY: phip
256 USE angle_constraints,
ONLY: scalfor_rho
262 REAL(dp),
DIMENSION(ns,0:ntor,0:mpol1,ntmax),
INTENT(inout) :: &
267 INTEGER,
PARAMETER :: n0=0, m0=0, m1=1
268 INTEGER,
PARAMETER :: n3d=0, nasym=1
272 INTEGER :: nsfix, jedge, delIter
274 INTEGER :: i, j, k, l
284 IF (lfreeb .AND. lrfp)
THEN
286 IF (ictrl_prec2d .EQ. 0)
THEN
289 gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
307 CALL constrain_m1(gcr(:,:,m1,rss), gcz(:,:,m1,zcs))
310 CALL constrain_m1(gcr(:,:,m1,rsc), gcz(:,:,m1,zcc))
314 IF (lfreeb .AND. lrfp)
THEN
316 IF (ictrl_prec2d .EQ. 0)
THEN
319 gcr(ns,0,m0,:) = fac*gcr(ns,0,m0,:)
320 gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
325 IF (ictrl_prec2d .GE. 2 .OR.
326 ictrl_prec2d .EQ. -1)
THEN
333 r1 = one/(2*r0scale)**2
337 deliter = iter2 - iter1
339 IF (deliter .lt. 50 .and.
340 (fsqr + fsqz) .lt. 1.e-6_dp)
THEN
344 CALL getfsq (gcr, gcz, fsqr, fsqz, r1*fnorm, jedge)
346 fsql = fnorml*sum(gcl*gcl)
347 fedge = r1*fnorm*sum(gcr(ns,:,:,:)**2 + gcz(ns,:,:,:)**2)
353 IF (ictrl_prec2d .EQ. 1)
THEN
355 CALL block_precond(gc)
357 IF (.not.lfreeb .and. any(gcr(ns,:,:,:) .ne. zero))
THEN
358 stop
'gcr(ns) != 0 for fixed boundary in residue'
360 IF (.not.lfreeb .and. any(gcz(ns,:,:,:) .ne. zero))
THEN
361 stop
'gcz(ns) != 0 for fixed boundary in residue'
363 IF (any(gcl(:,1:,0,zsc) .ne. zero))
THEN
364 stop
'gcl(m=0,n>0,sc) != 0 in residue'
366 IF (lthreed .and. any(gcl(:,n0,:,zcs) .ne. zero))
THEN
367 stop
'gcl(n=0,m,cs) != 0 in residue'
376 CALL scalfor_rho(gcr, gcz)
380 CALL scale_m1(gcr(:,:,m1,rss), gcz(:,:,m1,zcs))
383 CALL scale_m1(gcr(:,:,m1,rsc), gcz(:,:,m1,zcc))
387 CALL scalfor (gcr, arm, brm, ard, brd, crd, jedge)
389 CALL scalfor (gcz, azm, bzm, azd, bzd, crd, jedge)
393 CALL getfsq (gcr, gcz, fsqr1, fsqz1, fnorm1, m1)
399 IF ((iter2 - iter1) .LT. 25 .AND.
400 (fsqr + fsqz) .GT. 1.e-2_dp)
THEN
401 fac = fac / sqrt(1.e2_dp*(fsqr+fsqz))
411 fsql1 = hs*sum(gcl*gcl)
416 END SUBROUTINE residue
418 SUBROUTINE constrain_m1(gcr, gcz)
420 USE precon2d,
ONLY: ictrl_prec2d
425 REAL(dp),
DIMENSION(ns,0:ntor),
INTENT(INOUT) :: gcr, gcz
429 REAL(dp),
PARAMETER :: FThreshold = 1.e-6_dp
430 REAL(dp) :: temp(ns,0:ntor)
439 gcr = osqrt2*(gcr + gcz)
440 gcz = osqrt2*(temp - gcz)
444 IF (fsqz .LT. fthreshold .OR.
446 & ictrl_prec2d .NE. 0)
THEN
450 END SUBROUTINE constrain_m1
452 SUBROUTINE scale_m1(gcr, gcz)
458 REAL(dp),
DIMENSION(ns,0:ntor),
INTENT(inout) :: gcr, gcz
462 INTEGER,
PARAMETER :: nodd=2
466 IF (.not.lconm1)
RETURN
468 fac = (ard(:,nodd) + brd(:,nodd))
469 / (ard(:,nodd) + brd(:,nodd) + azd(:,nodd) + bzd(:,nodd))
471 gcr(:,n) = fac*gcr(:,n)
474 fac = (azd(:,nodd) + bzd(:,nodd))
475 / (ard(:,nodd) + brd(:,nodd) + azd(:,nodd) + bzd(:,nodd))
477 gcz(:,n) = fac*gcz(:,n)
480 END SUBROUTINE scale_m1