V3FIT
check3d.f90
1  SUBROUTINE check3d_symmetry(a,b,c,ap,bp,cp,mblk,mblkc,nblocks, &
2  asym_index)
3  USE descriptor_mod, mm_loc=>mm, csrc_loc=>csrc
4 ! USE prof_mod
5 ! USE ptrd_mod
6  IMPLICIT NONE
7 
8 !Written 03-12-10 by Eduardo D'Azevedo
9 
10 !fake integer, parameter :: dp = kind(1.0d0)
11 !fake integer, parameter :: dlen_ = 9, m_=2,n_=3
12 !fake integer, dimension(DLEN_) :: descA, descA_1xp
13 !fake integer, parameter :: Locq = 10, Locp = 10
14 
15 !-----------------------------------------------
16 ! D u m m y A r g u m e n t s
17 !-----------------------------------------------
18  INTEGER, INTENT(in) :: nblocks, mblk, mblkc
19  REAL(dp), TARGET, DIMENSION(mblk,mblkc,nblocks), INTENT(inout) :: &
20  a, b, c
21  REAL(dp), TARGET, DIMENSION(Locq,Locp,nblocks), INTENT(inout) :: &
22  ap
23  REAL(dp), TARGET, DIMENSION(Locq,Locp,nblocks-1), INTENT(inout) ::&
24  bp, cp
25 
26  REAL(dp), INTENT(out) :: asym_index
27 #if defined(MPI_OPT)
28 !-----------------------------------------------
29 ! L o c a l P a r a m e t e r s
30 !-----------------------------------------------
31  INTEGER :: k
32 
33  REAL(dp) :: w1, w2
34 
35  REAL(dp), parameter :: tol = 1.0d-6
36 
37  REAL(dp), POINTER :: amat(:,:), bmat(:,:), cmat(:,:)
38  REAL(dp), POINTER :: amatp(:,:), bmatp(:,:), cmatp(:,:)
39  REAL(dp), DIMENSION(:,:), POINTER :: Asrc, Bsrc, Csrc
40  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Atmp
41  INTEGER, DIMENSION(DLEN_) :: descAtmp
42 
43 ! ---------------------------
44 ! check a(:,:,:) - true - or ap(:,:,:) - false
45 ! ---------------------------
46  LOGICAL, PARAMETER :: check_matp = .false.
47 
48  INTEGER :: mm,nn
49  REAL(dp) :: dnorm, anorm,bnorm,cnorm, alpha,beta
50  REAL(dp), EXTERNAL :: pdlange
51  REAL(dp), DIMENSION(1) :: work
52 !-----------------------------------------------
53  IF (.NOT.lscalapack) RETURN
54 
55  w1 = 0
56  w2 = 0
57 ! -----------------------------
58 ! setup storage for temp matrix
59 ! -----------------------------
60  IF (check_matp) THEN
61  ALLOCATE( atmp(locq,locp) )
62  descatmp = desca
63  ELSE
64  ALLOCATE( atmp(mblk,mblkc) )
65  descatmp = desca_1xp
66  ENDIF
67 
68 ! ---------------------
69 ! check diagonal blocks
70 ! ---------------------
71  DO k=1,nblocks
72  amat => a(:,:,k)
73  amatp => ap(:,:,k)
74  IF (check_matp) THEN
75  asrc => amatp
76  ELSE
77  asrc => amat
78  ENDIF
79 
80 
81 ! ---------------------------
82 ! copy diagonal block to Atmp
83 ! ---------------------------
84  atmp = 0.0d0
85 
86  alpha = 1.0d0
87  beta = 0.0d0
88  mm = descatmp(m_)
89  nn = descatmp(n_)
90  CALL pdgeadd('Notranspose', mm,nn, &
91  alpha, asrc,1,1,descatmp, &
92  beta, atmp,1,1,descatmp )
93 
94  mm = descatmp(m_)
95  nn = descatmp(n_)
96  anorm = pdlange('F',mm,nn,atmp,1,1,descatmp,work)
97 ! ----------------------
98 ! compute Atmp <- Atmp - transpose(Asrc)
99 ! ----------------------
100 
101  alpha = -1.0d0
102  beta = 1.0d0
103  mm = descatmp(m_)
104  nn = descatmp(n_)
105  CALL pdgeadd('Transpose', mm,nn, &
106  alpha, asrc,1,1,descatmp, &
107  beta, atmp,1,1,descatmp )
108 
109 
110  mm = descatmp(m_)
111  nn = descatmp(n_)
112  dnorm = pdlange( 'F', mm,nn, atmp,1,1,descatmp, work)
113  w1 = w1 + dnorm*dnorm
114  w2 = w2 + anorm*anorm
115 
116 ! -------------------------
117 ! check off-diagonal blocks
118 !
119 ! Bsrc = transpose(Csrc)
120 ! -----------------------------
121  IF (k < nblocks) THEN
122  bmat => b(:,:,k+1)
123  bmatp => bp(:,:,k)
124 
125  cmat => c(:,:,k)
126  cmatp => cp(:,:,k)
127 
128  IF (check_matp) THEN
129  bsrc => bmatp
130  csrc => cmatp
131  ELSE
132  bsrc => bmat
133  csrc => cmat
134  ENDIF
135 
136  mm = descatmp(m_)
137  nn = descatmp(n_)
138  bnorm = pdlange('F',mm,nn,bsrc,1,1,descatmp,work)
139  cnorm = pdlange('F',mm,nn,csrc,1,1,descatmp,work)
140 
141 ! ---------------------------
142 ! copy off-diagonal B block to Atmp
143 ! ---------------------------
144  atmp = 0.0d0
145 
146  alpha = 1.0d0
147  beta = 0.0d0
148  mm = descatmp(m_)
149  nn = descatmp(n_)
150  CALL pdgeadd('Notranspose', mm,nn, &
151  alpha,bsrc,1,1,descatmp, &
152  beta,atmp,1,1,descatmp )
153 
154 ! ----------------------
155 ! compute Atmp <- Atmp - transpose(Csrc)
156 ! ----------------------
157 
158  alpha = -1.0d0
159  beta = 1.0d0
160  mm = descatmp(m_)
161  nn = descatmp(n_)
162  CALL pdgeadd('Transpose', mm,nn, &
163  alpha,csrc,1,1,descatmp, &
164  beta,atmp,1,1,descatmp )
165 
166 
167  mm = descatmp(m_)
168  nn = descatmp(n_)
169  dnorm = pdlange( 'F', mm,nn, atmp,1,1,descatmp, work)
170 
171  w1 = w1 + dnorm*dnorm
172  w2 = w2 + bnorm*bnorm + cnorm*cnorm
173 
174  ENDIF
175 
176  ENDDO
177 
178  DEALLOCATE( atmp )
179 
180  IF (w2 .ne. 0) THEN
181  asym_index = sqrt(w1/w2)
182  ELSE
183  asym_index = -1
184  END IF
185 #else
186  asym_index = -1
187 #endif
188 
189  END SUBROUTINE check3d_symmetry
190 
191 
192  SUBROUTINE check3d_symmetry_serial(dblk,lblk,ublk,mpol,ntor, &
193  mblk,nblocks,asym_index)
194  USE stel_kinds
195  USE shared_data, ONLY: ndims
196  IMPLICIT NONE
197 !-----------------------------------------------
198 ! D u m m y A r g u m e n t s
199 !-----------------------------------------------
200  INTEGER, INTENT(in) :: mblk, nblocks, mpol, ntor
201  REAL(dp), INTENT(out) :: asym_index
202  REAL(dp), INTENT(in), DIMENSION(mblk,mblk,nblocks) :: &
203  dblk, lblk, ublk
204 !-----------------------------------------------
205 ! L o c a l P a r a m e t e r s
206 !-----------------------------------------------
207  INTEGER :: ntype, m, n, icol
208  REAL(dp) :: w1, w2
209 !-----------------------------------------------
210 !
211 ! COMPUTES asym_index = SQRT(||A - A^T||**2/||A||**2)
212 !
213  icol = 0
214 
215  w1 = 0; w2 = 0
216  DO ntype = 1, ndims
217  DO n = -ntor, ntor
218  DO m = 0, mpol
219  icol = icol+1
220  IF (m.EQ.0 .AND. n.LT.0) cycle
221  w1 = w1 + &
222  sum((ublk(icol,:,1:nblocks-1) &
223  - lblk(:,icol,2:nblocks))**2)
224  w2 = w2 + sum(ublk(icol,:,1:nblocks-1)**2 &
225  + lblk(:,icol,2:nblocks)**2)
226  w1 = w1 + sum((dblk(icol,:,:) - &
227  dblk(:,icol,:))**2)
228  w2 = w2 + sum(dblk(icol,:,:)**2)
229  END DO
230  END DO
231  END DO
232 
233  IF (w2 .NE. 0) THEN
234  asym_index = sqrt(w1/w2)
235  ELSE
236  asym_index = -1
237  END IF
238 
239  END SUBROUTINE check3d_symmetry_serial
240 
241 
242  SUBROUTINE checkforces(xc, gc)
243  USE stel_kinds, ONLY: dp
244  USE stel_constants, ONLY: zero, twopi
245  USE descriptor_mod, ONLY: iam
249  ndims
250  USE shared_functions, ONLY: funct_island
251  USE quantities, ONLY: fsubsmncf, fsubumnsf, fsubvmnsf, &
252  mpol, ntor, ns, wp0, hs_i
253  USE hessian, ONLY: l_compute_hessian
254  USE nscalingtools, ONLY: startglobrow, endglobrow
255  USE v3_utilities, ONLY: assert
256  IMPLICIT NONE
257 !-----------------------------------------------
258 ! D u m m y A r g u m e n t s
259 !-----------------------------------------------
260  REAL(dp) :: xc(ndims*(1+mpol)*(2*ntor+1),ns), &
261  gc(ndims*(1+mpol)*(2*ntor+1),ns)
262 !-----------------------------------------------
263 ! L o c a l P a r a m e t e r s
264 !-----------------------------------------------
265  REAL(dp), ALLOCATABLE :: gc1(:,:)
266  INTEGER :: ntype, n, m, icol, js, mblk, nt1, icount
267  REAL(dp) :: wplus, wmins, force, eps, fnorm
268  LOGICAL :: lsave, lsave2, lwrite, l_edge, l_s, l_u, l_v
269  CHARACTER*(10) :: force_array(3) = (/' FORCE-S: ',' FORCE-U: ',' FORCE-V: '/)
270 !-----------------------------------------------
271 !
272 ! COMPUTES F = -GRAD(WMHD) NUMERICALLY
273 !
274  mblk = SIZE(xc,1)
275  icol = 0
276  lsave = l_linearize
277  lsave2= l_compute_hessian
278  l_edge = l_push_edge
279  l_s = l_push_s
280  l_u = l_push_u
281  l_v = l_push_v
282  l_linearize = .false.
283  l_compute_hessian = .false.
284  l_init_state = .false.
285  l_applyprecon = .false.
286  l_getfsq = .false.
287  l_getwmhd = .true.
288  l_push_edge = .true.
289  l_push_s = .true.; l_push_u = .true.; l_push_v = .true.
290 
291  xc = 0
292  CALL funct_island
293  ALLOCATE (gc1(mblk,ns))
294 
295  fnorm = (twopi*twopi)*hs_i
296  gc1 = gc*fnorm
297 
298  IF (iam .EQ. 0) print *,'WRITING CHECK-FORCES FILES FORT.(3000+iam)'
299 
300  WRITE (3000+iam, *) 'fnorm: ', fnorm,' iam: ', iam
301 
302 ! eps = 10*SQRT(epsilon(eps))
303  eps = 1.e-3_dp
304 
305  CALL assert(wtotal.GT.zero,'WTOTAL=0!')
306 
307  !Need l_push_s,u,v = T for js=1; l_push_edge = T for js=ns
308 
309  radial: DO icount = 1, 4
310 
311  IF (icount .EQ. 1) js = 1
312  IF (icount .EQ. 2) js = 2
313  IF (icount .EQ. 3) js = ns/2
314  IF (icount .EQ. 4) js = ns
315 
316  lwrite = (startglobrow.LE.js .AND. endglobrow.GE.js)
317  IF (lwrite) &
318  WRITE (3000+iam, *)'POINT: ', js,' START: ', &
319  startglobrow,' END: ', endglobrow
320 
321  icol = 0
322 
323  DO ntype = 1, ndims
324  DO n = -ntor, ntor
325  DO m = 0, mpol
326  icol = icol+1
327 
328  xc(icol,js) = eps
329  CALL funct_island
330  wplus = wtotal
331 
332  xc(icol,js) = -eps
333  CALL funct_island
334  wmins = wtotal
335  force = -(wplus-wmins)/(2*eps)
336 
337  CALL assert(l_getwmhd,'L_GETWMHD = FALSE!')
338 
339  xc(icol,js) = 0
340 
341  IF (lwrite) THEN
342  nt1 = ntype
343  IF (ntype .GT. 3) nt1 = ntype-3
344  WRITE (3000+iam,100) m, n, ntype, force, force_array(nt1), &
345  gc1(icol,js)
346  END IF
347 
348  END DO
349  END DO
350  END DO
351 
352  IF (lwrite) WRITE (3000+iam,*)
353  CALL flush(3000+iam)
354 
355  END DO radial
356 
357 
358  100 FORMAT(1x,'M: ',i4,' N: ',i4,' NTYPE: ',i4,' GRAD-W: ',1pe14.4, a, 1pe14.4)
359 
360  CALL assert(icol.EQ.mblk,'icol != mblk in CheckForces')
361  l_linearize = lsave
362  l_compute_hessian = lsave2
363  l_push_edge = l_edge
364  l_push_s = l_s; l_push_u = l_u; l_push_v = l_v
365  xc = 0
366  DEALLOCATE (gc1)
367 
368  END SUBROUTINE checkforces
369 
370 
371  SUBROUTINE tridiag_test (xc, gc, eps)
372 !
373 ! TESTS Tridiagonal STRUCTURE OF HESSIAN at js=ns/2
374 !
375  USE stel_kinds, ONLY: dp
376  USE stel_constants, ONLY: zero
377  USE quantities, ONLY: fsubsmncf, fsubumnsf, fsubvmnsf, mpol, ntor, ns
378  USE hessian, ONLY: l_compute_hessian
379  USE shared_data, ONLY: ndims, l_linearize, l_init_state, &
381  USE shared_functions, ONLY: funct_island
382  USE descriptor_mod, ONLY: iam
383  USE v3_utilities, ONLY: assert
384  IMPLICIT NONE
385 !-----------------------------------------------
386 ! D u m m y A r g u m e n t s
387 !-----------------------------------------------
388  REAL(dp), DIMENSION(ndims*(1+mpol)*(2*ntor+1),ns) :: &
389  xc, gc
390  REAL(dp), INTENT(IN) :: eps
391 !-----------------------------------------------
392 ! L o c a l P a r a m e t e r s
393 !-----------------------------------------------
394  INTEGER :: js, icol, ntype, n, m
395  REAL(dp), ALLOCATABLE :: gc0(:,:)
396  LOGICAL :: l_save, l_AppSave
397 !-----------------------------------------------
398 
399  js = ns/2
400  l_save=l_linearize
401  l_linearize=.false.
402  l_init_state = .true.
403  l_appsave = l_applyprecon
404  l_applyprecon = .false.
405  l_getwmhd = .true.
406 
407  ALLOCATE (gc0(SIZE(gc,1), SIZE(gc,2)))
408  xc = 0
409  CALL funct_island
410  gc0 = gc
411  l_getwmhd = .false.
412 
413  icol = 0
414  DO ntype = 1, ndims
415  DO n = -ntor, ntor
416  DO m = 0, mpol
417  icol = icol+1
418  IF (m .eq. 0 .and. n.lt.0) cycle
419  xc(icol,js) = eps
420  CALL funct_island
421  xc(icol,js) = 0
422  gc = gc - gc0
423  CALL assert(all(gc(:,1:js-2) .EQ. zero), &
424  'TRI_DIAGONAL FAILED FOR LOWER BANDS')
425  CALL assert(all(gc(:,js+2:ns) .EQ. zero), &
426  'TRI_DIAGONAL FAILED FOR UPPER BANDS')
427  END DO
428  END DO
429  END DO
430 
431  IF (iam .EQ. 0) print *,'TRI_DIAGONAL TEST PASSED'
432  l_linearize=l_save
433  l_applyprecon = l_appsave
434  DEALLOCATE (gc0)
435 
436  END SUBROUTINE tridiag_test
shared_data::l_push_edge
logical l_push_edge
Solve u,v components at s=1.
Definition: shared_data.f90:202
shared_data::l_getfsq
logical l_getfsq
Compute |F|^2.
Definition: shared_data.f90:216
shared_functions
Module contained subroutines and functions updating MHD forces and Wmhd.
Definition: shared_functions.f90:4
shared_data::l_push_u
logical l_push_u
Solve for u component at origin.
Definition: shared_data.f90:206
quantities
This file contains subroutines for allocating and initializing curvilinear magnetic covariant and pre...
Definition: quantities.f90:11
shared_data::l_push_s
logical l_push_s
Solve for s component at origin.
Definition: shared_data.f90:204
shared_data::l_applyprecon
logical l_applyprecon
Apply preconditioner.
Definition: shared_data.f90:218
shared_data::l_init_state
logical l_init_state
Store initial field/pressure state.
Definition: shared_data.f90:222
v3_utilities::assert
Definition: v3_utilities.f:55
shared_data::wtotal
real(dp) wtotal
MHD energy sum of magnetic and thermal.
Definition: shared_data.f90:121
shared_data::l_linearize
logical l_linearize
Use linearized forces.
Definition: shared_data.f90:210
shared_data::l_getwmhd
logical l_getwmhd
Compute MHD energy.
Definition: shared_data.f90:214
shared_data::ndims
integer ndims
Number of independent variables.
Definition: shared_data.f90:58
shared_data
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.f90:10
shared_data::l_push_v
logical l_push_v
Solve for v component at origin.
Definition: shared_data.f90:208