V3FIT
residue.f90
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
6  USE xstuff
7  USE precon2d
8  USE parallel_include_module
9  USE parallel_vmec_module, ONLY: tlglob_arr, trglob_arr, &
10  lactive, saxlastntype
11  USE blocktridiagonalsolver, ONLY: l_colscale
12 
13  IMPLICIT NONE
14 !-----------------------------------------------
15 ! D u m m y A r g u m e n t s
16 !-----------------------------------------------
17  REAL(dp), DIMENSION(0:ntor,0:mpol1,ns,ntmax), INTENT(INOUT) :: &
18  gcr, gcz, gcl
19 !-----------------------------------------------
20 ! L o c a l P a r a m e t e r s
21 !-----------------------------------------------
22  INTEGER, PARAMETER :: n0=0, m0=0, m1=1
23  INTEGER, PARAMETER :: n3d=0, nasym=1
24 !-----------------------------------------------
25 ! L o c a l V a r i a b l e s
26 !-----------------------------------------------
27  INTEGER :: nsfix, jedge, delIter
28  REAL(dp) :: r1, fac, tmp, tmp2(ns), ftotal
29 
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
36 !-----------------------------------------------
37  CALL second0 (treson)
38 
39 #ifdef _HBANGLE
40 !FREE-BDY RFP MAY NEED THIS TO IMPROVE CONVERGENCE (SPH 022514)
41  IF (lfreeb .AND. lrfp) THEN
42  fac = 0
43  IF (ictrl_prec2d .EQ. 0) THEN
44  fac = 1.e-1_dp
45  END IF
46  gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
47  END IF
48 #else
49 !
50 ! SYMMETRIC PERTURBATIONS (BASED ON POLAR RELATIONS):
51 ! Rss(n) = Zcs(n), n != 0
52 ! ASYMMETRIC PERTURBATIONS:
53 ! Rsc(n) = Zcc(n), ALL n
54 !
55 ! INTERNALLY:
56 ! XC(rss) = .5*(Rss + Zcs), XC(zcs) = .5*(Rss - Zcs) -> 0
57 ! XC(rsc) = .5*(Rsc + Zcc), XC(zcc) = .5*(Rsc - Zcc) -> 0
58 ! THIS IMPLIES THE CONSTRAINT
59 ! 3D ONLY : GC(zcs) = 0;
60 ! ASYM: GC(zcc) = 0
61 !
62 
63  IF (lthreed) THEN
64  CALL constrain_m1_par(gcr(:,m1,:,rss), gcz(:,m1,:,zcs))
65  END IF
66  IF (lasym) THEN
67  CALL constrain_m1_par(gcr(:,m1,:,rsc), gcz(:,m1,:,zcc))
68  END IF
69 
70  IF (lfreeb .AND. lrfp) THEN
71  fac = 0
72  IF (ictrl_prec2d .EQ. 0) THEN
73  fac = 1.e-1_dp
74  END IF
75  gcr(0,m0,ns,:) = fac*gcr(0,m0,ns,:)
76  gcz(0,m0,ns,:) = fac*gcz(0,m0,ns,:)
77  END IF
78 #endif
79 
80 ! PRECONDITIONER MUST BE CALCULATED USING RAW (UNPRECONDITIONED) FORCES
81  IF (ictrl_prec2d .GE. 2 .OR. ictrl_prec2d .EQ. -1) RETURN
82 
83 !
84 ! COMPUTE INVARIANT RESIDUALS
85 !
86  r1 = one/(2*r0scale)**2
87  jedge = 0
88  deliter = iter2-iter1
89 
90  IF (deliter .lt. 50 .and. &
91  (fsqr + fsqz) .LT. 1.e-6_dp) THEN
92  jedge = 1
93  ENDIF
94 
95  CALL getfsq_par (gcr, gcz, fsqr, fsqz, r1*fnorm, jedge)
96 
97  CALL second0(tredon)
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)
102  fsql = fnorml*ftotal
103  IF(rank .EQ. nranks-1) THEN
104  fedge = r1*fnorm*sum(gcr(:,:,ns,:)**2 + gcz(:,:,ns,:)**2)
105  END IF
106 !
107 ! PERFORM PRECONDITIONING AND COMPUTE RESIDUES
108 !
109  IF (ictrl_prec2d .EQ. 1) THEN
110 
111  IF (l_colscale .AND. lactive) THEN
112  CALL saxlastntype(pgc, pcol_scale, pgc)
113  END IF
114 
115  lresiduecall = .true.
116  CALL block_precond_par(pgc)
117  lresiduecall = .false.
118 
119  IF (.NOT.lfreeb .AND. any(gcr(:,:,ns,:) .NE. zero)) THEN
120  stop 'gcr(ns) != 0 for fixed boundary in residue'
121  END IF
122  IF (.NOT.lfreeb .AND. any(gcz(:,:,ns,:) .NE. zero)) THEN
123  stop 'gcz(ns) != 0 for fixed boundary in residue'
124  END IF
125  IF (any(gcl(1:,m0,:,zsc) .NE. zero)) THEN
126  stop 'gcl(m=0,n>0,sc) != 0 in residue'
127  END IF
128  IF (lthreed .AND. any(gcl(n0,:,:,zcs) .NE. zero)) THEN
129  stop 'gcl(n=0,m,cs) != 0 in residue'
130  END IF
131 
132  fsqr1 = sum(gcr*gcr)
133  fsqz1 = sum(gcz*gcz)
134  fsql1 = sum(gcl*gcl)
135 
136  ELSE
137 ! m = 1 constraint scaling
138 
139  IF (lthreed) THEN
140  CALL scale_m1_par(gcr(:,m1,:,rss), gcz(:,m1,:,zcs))
141  END IF
142  IF (lasym) THEN
143  CALL scale_m1_par(gcr(:,m1,:,rsc), gcz(:,m1,:,zcc))
144  END IF
145 
146  jedge = 0
147  CALL scalfor_par (gcr, arm, brm, ard, brd, crd, jedge)
148  jedge = 1
149  CALL scalfor_par (gcz, azm, bzm, azd, bzd, crd, jedge)
150 
151  CALL getfsq_par (gcr, gcz, fsqr1, fsqz1, fnorm1, m1)
152 
153  DO l = tlglob, trglob
154  gcl(:,:,l,:) = pfaclam(:,:,l,:)*gcl(:,:,l,:)
155  tmp2(l) = sum(gcl(:,:,l,:)**2)
156  END DO
157  CALL gather1xarray(tmp2)
158  ftotal = sum(tmp2(2:ns))
159  fsql1 = hs*ftotal
160 
161  CALL padsides(pgc)
162 
163  ENDIF
164 
165  CALL second0 (tresoff)
166  residue_time = residue_time + (tresoff-treson)
167 
168  END SUBROUTINE residue_par
169 
170  SUBROUTINE constrain_m1_par(gcr, gcz)
171  USE vmec_main
172  USE parallel_include_module
173  USE precon2d, ONLY: ictrl_prec2d
174  IMPLICIT NONE
175 !-----------------------------------------------
176 ! D u m m y A r g u m e n t s
177 !-----------------------------------------------
178  REAL(dp), DIMENSION(0:ntor,ns), INTENT(INOUT) :: gcr, gcz
179 !-----------------------------------------------
180 ! L o c a l P a r a m e t e r s
181 !-----------------------------------------------
182  REAL(dp), PARAMETER :: FThreshold = 1.e-6_dp
183  REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: temp
184 !-----------------------------------------------
185 !
186 ! COMPUTE INTERNAL gr, gz
187 ! NOTE: internal gz => 0 for both values of lconm1 (although gz is different)
188 ! FOR lconm1=T, gcr(internal) = gcr+gcz, gcz(internal) = gcr-gcz->0
189 !
190  ALLOCATE(temp(0:ntor,ns))
191  IF (lconm1) THEN
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))
197  END IF
198 
199 !v8.50: ADD iter2<2 so reset=<WOUT_FILE> works
200  IF (fsqz .LT. fthreshold .OR. &
201  & iter2 .LT. 2 .OR. &
202  & ictrl_prec2d .NE. 0) THEN
203  gcz(:,tlglob:trglob) = 0
204  END IF
205 
206  DEALLOCATE(temp)
207  END SUBROUTINE constrain_m1_par
208 
209  SUBROUTINE scale_m1_par(gcr, gcz)
210  USE vmec_main
211  USE parallel_include_module
212  IMPLICIT NONE
213 !-----------------------------------------------
214 ! D u m m y A r g u m e n t s
215 !-----------------------------------------------
216  REAL(dp), DIMENSION(0:ntor,ns), INTENT(inout) :: gcr, gcz
217 !-----------------------------------------------
218 ! L o c a l P a r a m e t e r s
219 !-----------------------------------------------
220  INTEGER, PARAMETER :: nodd=2
221  INTEGER :: n
222  REAL(dp) :: fac(ns)
223 !-----------------------------------------------
224  IF (.not.lconm1) RETURN
225 
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))
232  DO n = 0, ntor
233  gcr(n,tlglob:trglob) = fac(tlglob:trglob)*gcr(n,tlglob:trglob)
234  END DO
235 
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))
242  DO n = 0, ntor
243  gcz(n,tlglob:trglob) = fac(tlglob:trglob)*gcz(n,tlglob:trglob)
244  END DO
245 
246  END SUBROUTINE scale_m1_par
247 
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
253  USE xstuff
254  USE precon2d
255 #ifdef _HBANGLE
256  USE angle_constraints, ONLY: scalfor_rho
257 #endif
258  IMPLICIT NONE
259 !-----------------------------------------------
260 ! D u m m y A r g u m e n t s
261 !-----------------------------------------------
262  REAL(dp), DIMENSION(ns,0:ntor,0:mpol1,ntmax), INTENT(inout) :: &
263  gcr, gcz, gcl
264 !-----------------------------------------------
265 ! L o c a l P a r a m e t e r s
266 !-----------------------------------------------
267  INTEGER, PARAMETER :: n0=0, m0=0, m1=1
268  INTEGER, PARAMETER :: n3d=0, nasym=1
269 !-----------------------------------------------
270 ! L o c a l V a r i a b l e s
271 !-----------------------------------------------
272  INTEGER :: nsfix, jedge, delIter
273  REAL(dp) :: r1, fac
274  INTEGER :: i, j, k, l
275 !-----------------------------------------------
276 !
277 ! IMPOSE M=1 MODE CONSTRAINT TO MAKE THETA ANGLE
278 ! INVARIANT TO PHI-SHIFTS (AND THETA SHIFTS FOR ASYMMETRIC CASE)
279 ! (ZCS = RSS, ZSS = RCS ARE THE CORRECT POLAR RELATIONS)
280 !
281 
282 #ifdef _HBANGLE
283 !FREE-BDY RFP MAY NEED THIS TO IMPROVE CONVERGENCE (SPH 022514)
284  IF (lfreeb .AND. lrfp) THEN
285  fac = 0
286  IF (ictrl_prec2d .EQ. 0) THEN
287  fac = 1.e-1_dp
288  END IF
289  gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
290  END IF
291 #else
292 !
293 ! SYMMETRIC PERTURBATIONS (BASED ON POLAR RELATIONS):
294 ! RSS(n) = ZCS(n), n != 0
295 ! ASYMMETRIC PERTURBATIONS:
296 ! RSC(n) = ZCC(n), ALL n
297 !
298 ! INTERNALLY:
299 ! XC(rss) = .5*(Rss + Zcs), XC(zcs) = .5*(Rss - Zcs) -> 0
300 ! XC(rsc) = .5*(Rsc + Zcc), XC(zcc) = .5*(Rsc - Zcc) -> 0
301 ! THIS IMPLIES THE CONSTRAINT
302 ! 3D ONLY : GC(zcs) = 0;
303 ! ASYM: GC(zcc) = 0
304 !
305 
306  IF (lthreed) THEN
307  CALL constrain_m1(gcr(:,:,m1,rss), gcz(:,:,m1,zcs))
308  END IF
309  IF (lasym) THEN
310  CALL constrain_m1(gcr(:,:,m1,rsc), gcz(:,:,m1,zcc))
311  END IF
312 
313 !FREE-BDY RFP MAY NEED THIS TO IMPROVE CONVERGENCE (SPH 022514)
314  IF (lfreeb .AND. lrfp) THEN
315  fac = 0
316  IF (ictrl_prec2d .EQ. 0) THEN
317  fac = 1.e-1_dp
318  END IF
319  gcr(ns,0,m0,:) = fac*gcr(ns,0,m0,:)
320  gcz(ns,0,m0,:) = fac*gcz(ns,0,m0,:)
321  END IF
322 #endif
323 
324 ! PRECONDITIONER MUST BE CALCULATED USING RAW (UNPRECONDITIONED) FORCES
325  IF (ictrl_prec2d .GE. 2 .OR. &
326  ictrl_prec2d .EQ. -1) THEN
327  RETURN
328  END IF
329 
330 !
331 ! COMPUTE INVARIANT RESIDUALS
332 !
333  r1 = one/(2*r0scale)**2
334  jedge = 0
335 !SPH-JAH013108: MUST INCLUDE EDGE FORCE (INITIALLY) FOR V3FITA TO WORK
336 !ADD A V3FIT RELATED FLAG? ADD fsq criterion first
337  deliter = iter2 - iter1
338 
339  IF (deliter .lt. 50 .and. &
340  (fsqr + fsqz) .lt. 1.e-6_dp) THEN
341  jedge = 1
342  END IF
343 
344  CALL getfsq (gcr, gcz, fsqr, fsqz, r1*fnorm, jedge)
345 
346  fsql = fnorml*sum(gcl*gcl)
347  fedge = r1*fnorm*sum(gcr(ns,:,:,:)**2 + gcz(ns,:,:,:)**2)
348 
349 !
350 ! PERFORM PRECONDITIONING AND COMPUTE RESIDUES
351 !
352 
353  IF (ictrl_prec2d .EQ. 1) THEN
354 
355  CALL block_precond(gc)
356 
357  IF (.not.lfreeb .and. any(gcr(ns,:,:,:) .ne. zero)) THEN
358  stop 'gcr(ns) != 0 for fixed boundary in residue'
359  END IF
360  IF (.not.lfreeb .and. any(gcz(ns,:,:,:) .ne. zero)) THEN
361  stop 'gcz(ns) != 0 for fixed boundary in residue'
362  END IF
363  IF (any(gcl(:,1:,0,zsc) .ne. zero)) THEN
364  stop 'gcl(m=0,n>0,sc) != 0 in residue'
365  END IF
366  IF (lthreed .and. any(gcl(:,n0,:,zcs) .ne. zero)) THEN
367  stop 'gcl(n=0,m,cs) != 0 in residue'
368  END IF
369 
370  fsqr1 = sum(gcr*gcr)
371  fsqz1 = sum(gcz*gcz)
372  fsql1 = sum(gcl*gcl)
373 
374  ELSE
375 #ifdef _HBANGLE
376  CALL scalfor_rho(gcr, gcz)
377 #else
378 ! m = 1 constraint scaling
379  IF (lthreed) THEN
380  CALL scale_m1(gcr(:,:,m1,rss), gcz(:,:,m1,zcs))
381  END IF
382  IF (lasym) THEN
383  CALL scale_m1(gcr(:,:,m1,rsc), gcz(:,:,m1,zcc))
384  END IF
385 
386  jedge = 0
387  CALL scalfor (gcr, arm, brm, ard, brd, crd, jedge)
388  jedge = 1
389  CALL scalfor (gcz, azm, bzm, azd, bzd, crd, jedge)
390 #endif
391 
392 !SPH: add fnorm1 ~ 1/R**2, since preconditioned forces gcr,gcz ~ Rmn or Zmn
393  CALL getfsq (gcr, gcz, fsqr1, fsqz1, fnorm1, m1)
394 #ifdef _HBANGLE
395 !
396 ! TO IMPROVE CONVERGENCE, REDUCE FORCES INITIALLY IF THEY ARE TOO LARGE
397 !
398  fac = .5_dp
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))
402  END IF
403  gcr = fac*gcr
404  gcz = fac*gcz
405 #endif
406 
407 !SPH: THIS IS NOT INVARIANT UNDER PHIP->A*PHIP, AM->A**2*AM IN PROFIL1D
408 ! (EXTCUR -> A*EXTCUR for FREE BOUNDARY)
409 
410  gcl = faclam*gcl
411  fsql1 = hs*sum(gcl*gcl)
412 !030514 fsql1 = hs*lamscale**2*SUM(gcl*gcl)
413 
414  END IF
415 
416  END SUBROUTINE residue
417 
418  SUBROUTINE constrain_m1(gcr, gcz)
419  USE vmec_main
420  USE precon2d, ONLY: ictrl_prec2d
421  IMPLICIT NONE
422 !-----------------------------------------------
423 ! D u m m y A r g u m e n t s
424 !-----------------------------------------------
425  REAL(dp), DIMENSION(ns,0:ntor), INTENT(INOUT) :: gcr, gcz
426 !-----------------------------------------------
427 ! L o c a l P a r a m e t e r s
428 !-----------------------------------------------
429  REAL(dp), PARAMETER :: FThreshold = 1.e-6_dp
430  REAL(dp) :: temp(ns,0:ntor)
431 !-----------------------------------------------
432 !
433 ! COMPUTE INTERNAL gr, gz
434 ! NOTE: internal gz => 0 for both values of lconm1 (although gz is different)
435 ! FOR lconm1=T, gcr(internal) = gcr+gcz, gcz(internal) = gcr-gcz->0
436 !
437  IF (lconm1) THEN
438  temp = gcr
439  gcr = osqrt2*(gcr + gcz)
440  gcz = osqrt2*(temp - gcz)
441  END IF
442 
443 !v8.50: ADD iter2<2 so reset=<WOUT_FILE> works
444  IF (fsqz .LT. fthreshold .OR. &
445  & iter2 .LT. 2 .OR. &
446  & ictrl_prec2d .NE. 0) THEN
447  gcz = 0
448  END IF
449 
450  END SUBROUTINE constrain_m1
451 
452  SUBROUTINE scale_m1(gcr, gcz)
453  USE vmec_main
454  IMPLICIT NONE
455 !-----------------------------------------------
456 ! D u m m y A r g u m e n t s
457 !-----------------------------------------------
458  REAL(dp), DIMENSION(ns,0:ntor), INTENT(inout) :: gcr, gcz
459 !-----------------------------------------------
460 ! L o c a l P a r a m e t e r s
461 !-----------------------------------------------
462  INTEGER, PARAMETER :: nodd=2
463  INTEGER :: n
464  REAL(dp) :: fac(ns)
465 !-----------------------------------------------
466  IF (.not.lconm1) RETURN
467 
468  fac = (ard(:,nodd) + brd(:,nodd)) &
469  / (ard(:,nodd) + brd(:,nodd) + azd(:,nodd) + bzd(:,nodd))
470  DO n = 0, ntor
471  gcr(:,n) = fac*gcr(:,n)
472  END DO
473 
474  fac = (azd(:,nodd) + bzd(:,nodd)) &
475  / (ard(:,nodd) + brd(:,nodd) + azd(:,nodd) + bzd(:,nodd))
476  DO n = 0, ntor
477  gcz(:,n) = fac*gcz(:,n)
478  END DO
479 
480  END SUBROUTINE scale_m1
481 
blocktridiagonalsolver
Solver for block tri-diagonal matrices. [Kalyan S. Perumalla, ORNL, 2009-2011].
Definition: blocktridiagonalsolver.f90:8