V3FIT
precondn.f
1  ! SKS-RANGE: All OUT arrays computed correctly between
2  ! [tlglob, trglob]
3  SUBROUTINE precondn_par(lu1, bsq, gsqrt, r12, xs, xu12, xue, xuo,
4  & xodd, axm, axd, bxm, bxd,
5  & cx, eqfactor, trigmult)
6  USE vmec_main
7  USE vmec_params, ONLY: signgs
8  USE realspace, ONLY: pshalf, pwint
9  USE parallel_include_module
10  IMPLICIT NONE
11 C-----------------------------------------------
12 C D u m m y A r g u m e n t s
13 C-----------------------------------------------
14  REAL(rprec), DIMENSION(nznt,ns), INTENT(in) ::
15  & lu1, bsq, gsqrt, r12, xs, xu12, xue, xuo, xodd
16  REAL(rprec), DIMENSION(ns+1,2), INTENT(out) ::
17  & axm, axd, bxm, bxd
18  REAL(rprec), DIMENSION(ns+1), INTENT(out) :: cx
19  REAL(rprec), DIMENSION(ns), INTENT(out) :: eqfactor
20  REAL(rprec), DIMENSION(nznt), INTENT(in) :: trigmult
21 C-----------------------------------------------
22 C L o c a l V a r i a b l e s
23 C-----------------------------------------------
24  INTEGER :: js, l, lk
25  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: ax, bx
26  !REAL(rprec) :: temp(ns+1)
27  REAL(rprec), ALLOCATABLE, DIMENSION(:) :: temp
28  REAL(rprec), DIMENSION(:), ALLOCATABLE :: ptau, ptau2
29  REAL(rprec) :: t1, t2, t3, pfactor
30 
31  INTEGER :: nsmin, nsmax, i, j, k, numjs
32  INTEGER, DIMENSION(:), ALLOCATABLE :: ldisps, lcounts
33  REAL(rprec), DIMENSION(:), ALLOCATABLE :: sbuf
34 
35  nsmin=tlglob; nsmax=t1rglob
36  ! Correct incoming lu1, bsq, gsqrt, , xue, xuo, xodd, trigmult
37  ! between [t1lglob, t1rglob]
38 
39  ALLOCATE (ax(ns+1,4), bx(ns+1,4),
40  & ptau(nznt), ptau2(nznt), temp(ns+1))
41  ax = 0
42  bx = 0
43  cx = 0
44  temp = 0
45  pfactor = -4*r0scale**2 !restored in v8.51
46 
47  nsmin = max(2,tlglob)
48  nsmax = t1rglob
49  DO js = nsmin, nsmax
50 !
51 ! COMPUTE DOMINANT (1/DELTA-S)**2 PRECONDITIONING
52 ! MATRIX ELEMENTS
53 !
54  lk = 0
55  DO k = 1, ntheta3
56  DO j = 1, nzeta
57  lk = lk + 1
58  t1 = pfactor*r12(lk,js)*bsq(lk,js)
59  ptau2(lk) = r12(lk,js)*t1/gsqrt(lk,js)
60  t1 = t1*pwint(lk,js)
61  temp(js) = temp(js) + t1*trigmult(lk)*xu12(lk,js)
62  ptau(lk) = r12(lk,js)*t1/gsqrt(lk,js)
63  t1 = xu12(lk,js)*ohs
64  t2 = cp25*(xue(lk,js)/pshalf(lk,js)+xuo(lk,js))
65  & / pshalf(lk,js)
66  t3 = cp25*(xue(lk,js-1)/pshalf(lk,js) +
67  & xuo(lk,js-1))/pshalf(lk,js)
68  ax(js,1) = ax(js,1) + ptau(lk)*t1*t1
69  ax(js,2) = ax(js,2) + ptau(lk)*(-t1+t3)*(t1+t2)
70  ax(js,3) = ax(js,3) + ptau(lk)*(t1+t2)*(t1+t2)
71  ax(js,4) = ax(js,4) + ptau(lk)*(-t1+t3)*(-t1+t3)
72  END DO
73  END DO
74 !
75 ! COMPUTE PRECONDITIONING MATRIX ELEMENTS FOR M**2, N**2 TERMS
76 !
77  lk = 0
78  DO k = 1, ntheta3
79  DO j = 1, nzeta
80  lk = lk+1
81  t1 = cp5*(xs(lk,js) + cp5*xodd(lk,js)/pshalf(lk,js))
82  t2 = cp5*(xs(lk,js) + cp5*xodd(lk,js-1)/pshalf(lk,js))
83  bx(js,1) = bx(js,1) + ptau(lk)*t1*t2
84  bx(js,2) = bx(js,2) + ptau(lk)*t1*t1
85  bx(js,3) = bx(js,3) + ptau(lk)*t2*t2
86  cx(js) = cx(js)
87  & + cp25*pfactor*lu1(lk,js)**2 *
88  & gsqrt(lk,js)*pwint(lk,js)
89  END DO
90  END DO
91 
92  END DO
93 
94  nsmin = max(2,tlglob)
95  nsmax = t1rglob
96  temp(1) = 0
97  temp(nsmin:nsmax) = temp(nsmin:nsmax)/vp(nsmin:nsmax);
98  temp(ns+1) = 0
99 
100  nsmin = t1lglob
101  nsmax = t1rglob
102  DO js = nsmin, nsmax
103  axm(js,1) = -ax(js,1)
104  axd(js,1) = ax(js,1) + ax(js+1,1)
105  axm(js,2) = ax(js,2) * sm(js) * sp(js-1)
106  axd(js,2) = ax(js,3)*sm(js)**2 + ax(js+1,4)*sp(js)**2
107  bxm(js,1) = bx(js,1)
108  bxm(js,2) = bx(js,1) * sm(js) * sp(js-1)
109  bxd(js,1) = bx(js,2) + bx(js+1,3)
110  bxd(js,2) = bx(js,2)*sm(js)**2 + bx(js+1,3)*sp(js)**2
111  cx(js) = cx(js) + cx(js+1)
112  temp(js) = signgs*(temp(js) + temp(js+1))
113  END DO
114 
115  nsmin = max(2,tlglob)
116  nsmax = min(ns-1,trglob)
117  eqfactor(nsmin:nsmax) = axd(nsmin:nsmax,2)*hs*hs/
118  & temp(nsmin:nsmax)
119  eqfactor(1) = 0
120  eqfactor(ns) = 0
121  axm(ns+1,:) = 0
122  axd(ns+1,:) = 0
123  bxm(ns+1,:) = 0
124  bxd(ns+1,:) = 0
125  DEALLOCATE (ax, bx, ptau, ptau2, temp)
126 
127  END SUBROUTINE precondn_par
128 
129  SUBROUTINE precondn(lu1, bsq, gsqrt, r12, xs, xu12, xue, xuo,
130  & xodd, axm, axd, bxm, bxd,
131  & cx, eqfactor, trigmult)
132  USE vmec_main
133  USE vmec_params, ONLY: signgs
134  USE realspace
135 
136  IMPLICIT NONE
137 C-----------------------------------------------
138 C D u m m y A r g u m e n t s
139 C-----------------------------------------------
140  REAL(rprec), DIMENSION(nrzt), INTENT(in) ::
141  1 lu1, bsq, gsqrt, r12, xs, xu12, xue, xuo, xodd
142  REAL(rprec), DIMENSION(ns+1,2), INTENT(out) ::
143  1 axm, axd, bxm, bxd
144  REAL(rprec), DIMENSION(ns+1), INTENT(out) :: cx
145  REAL(rprec), DIMENSION(ns), INTENT(out) :: eqfactor
146  REAL(rprec), DIMENSION(nznt), INTENT(in) :: trigmult
147 C-----------------------------------------------
148 C L o c a l V a r i a b l e s
149 C-----------------------------------------------
150  INTEGER :: js, l, lk
151  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: ax, bx
152  REAL(rprec) :: temp(ns+1)
153  REAL(rprec), DIMENSION(:), ALLOCATABLE :: ptau, ptau2
154  REAL(rprec) :: t1, t2, t3, pfactor
155 
156 C-----------------------------------------------
157 !
158 ! COMPUTE PRECONDITIONING MATRIX ELEMENTS FOR R,Z
159 ! FORCE. NOTE THAT THE NORMALIZATION IS:
160 !
161 ! AX(off-diag) ~ <(cosmui cosmu cosnv cosnv) 2(R**2*Xu**2*bsq/gsqrt)>
162 ! Factor of 2 arising from 1/gsqrt**2 in bsq
163 !
164 ! Now, cosmui cosmu ~ mscale(0)**2, cosnv**2 ~ nscale(0)**2
165 ! Therefore, AX ~ (mscale(0)*nscale(0))**2 2<R**2*Xu**2*bsq/gsqrt>
166 ! ~ 2*r0scale**2 <...>
167 !
168  ALLOCATE (ax(ns+1,4), bx(ns+1,4), ptau(nznt), ptau2(nznt))
169  ax = 0; bx = 0; cx = 0
170  temp = 0
171 ! pfactor = -2*r0scale**2 !v8.50
172  pfactor = -4*r0scale**2 !restored in v8.51
173 
174  DO 20 js = 2,ns
175 !
176 ! COMPUTE DOMINANT (1/DELTA-S)**2 PRECONDITIONING
177 ! MATRIX ELEMENTS
178 !
179  lk = 0
180  DO l = js,nrzt,ns
181  lk = lk + 1
182  t1 = pfactor*r12(l)*bsq(l)
183  ptau2(lk) = r12(l)*t1/gsqrt(l)
184  t1 = t1*wint(l)
185  temp(js) = temp(js) + t1*trigmult(lk)*xu12(l)
186  ptau(lk) = r12(l)*t1/gsqrt(l)
187  t1 = xu12(l)*ohs
188  t2 = cp25*(xue(l)/shalf(js) + xuo(l))/shalf(js)
189  t3 = cp25*(xue(l-1)/shalf(js) + xuo(l-1))/shalf(js)
190  ax(js,1) = ax(js,1) + ptau(lk)*t1*t1
191  ax(js,2) = ax(js,2) + ptau(lk)*(-t1+t3)*(t1+t2)
192  ax(js,3) = ax(js,3) + ptau(lk)*(t1+t2)*(t1+t2)
193  ax(js,4) = ax(js,4) + ptau(lk)*(-t1+t3)*(-t1+t3)
194  END DO
195 !
196 ! COMPUTE PRECONDITIONING MATRIX ELEMENTS FOR M**2, N**2 TERMS
197 !
198  lk = 0
199  DO l = js,nrzt,ns
200  lk = lk+1
201  t1 = cp5*(xs(l) + cp5*xodd(l)/shalf(js))
202  t2 = cp5*(xs(l) + cp5*xodd(l-1)/shalf(js))
203  bx(js,1) = bx(js,1) + ptau(lk)*t1*t2
204  bx(js,2) = bx(js,2) + ptau(lk)*t1*t1
205  bx(js,3) = bx(js,3) + ptau(lk)*t2*t2
206  cx(js) = cx(js) + cp25*pfactor*lu1(l)**2*gsqrt(l)*wint(l)
207  END DO
208  20 CONTINUE
209 
210  temp(1) = 0
211  temp(2:ns) = temp(2:ns)/vp(2:ns)
212  temp(ns+1) = 0
213  DO js = 1,ns
214  axm(js,1) =-ax(js,1)
215  axd(js,1) = ax(js,1) + ax(js+1,1)
216  axm(js,2) = ax(js,2) * sm(js) * sp(js-1)
217  axd(js,2) = ax(js,3)*sm(js)**2 + ax(js+1,4)*sp(js)**2
218  bxm(js,1) = bx(js,1)
219  bxm(js,2) = bx(js,1) * sm(js) * sp(js-1)
220  bxd(js,1) = bx(js,2) + bx(js+1,3)
221  bxd(js,2) = bx(js,2)*sm(js)**2 + bx(js+1,3)*sp(js)**2
222  cx(js) = cx(js) + cx(js+1)
223  temp(js) = signgs*(temp(js) + temp(js+1))
224  END DO
225 
226  eqfactor(2:ns-1) = axd(2:ns-1,2)*hs*hs/temp(2:ns-1)
227  eqfactor(1) = 0
228  eqfactor(ns) = 0
229  axm(ns+1,:) = 0
230  axd(ns+1,:) = 0
231  bxm(ns+1,:) = 0
232  bxd(ns+1,:) = 0
233  DEALLOCATE (ax, bx, ptau, ptau2)
234 
235  END SUBROUTINE precondn
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17