KTH framework for Nek5000 toolboxes; testing version  0.0.1
eigsolv.f
Go to the documentation of this file.
1 C**********************************************************************
2 C
3 C ROUTINES FOR ESITMATING AND CALCULATING EIGENVALUES
4 C USED IN NEKTON
5 C
6 C**********************************************************************
7 C
8  SUBROUTINE esteig
9 C--------------------------------------------------------------
10 C
11 C Estimate eigenvalues
12 C
13 C-------------------------------------------------------------
14  include 'SIZE'
15  include 'GEOM'
16  include 'INPUT'
17  include 'EIGEN'
18  include 'TSTEP'
19 C
20  ntot1=lx1*ly1*lz1*nelfld(ifield)
21  xmin = glmin(xm1,ntot1)
22  xmax = glmax(xm1,ntot1)
23  ymin = glmin(ym1,ntot1)
24  ymax = glmax(ym1,ntot1)
25  IF (if3d) THEN
26  zmin = glmin(zm1,ntot1)
27  zmax = glmax(zm1,ntot1)
28  ELSE
29  zmin = 0.0
30  zmax = 0.0
31  ENDIF
32 C
33  xx = xmax - xmin
34  yy = ymax - ymin
35  zz = zmax - zmin
36  rxy = xx/yy
37  ryx = yy/xx
38  rmin = rxy
39  IF (ryx .LT. rmin) rmin = ryx
40  IF (ldim .EQ. 3) THEN
41  rxz = xx/zz
42  rzx = zz/xx
43  ryz = yy/zz
44  rzy = zz/yy
45  IF (rxz .LT. rmin) rmin = rxz
46  IF (rzx .LT. rmin) rmin = rzx
47  IF (ryz .LT. rmin) rmin = ryz
48  IF (rzy .LT. rmin) rmin = rzy
49  ENDIF
50 C
51  xx2 = 1./xx**2
52  yy2 = 1./yy**2
53  xyzmin = xx2
54  xyzmax = xx2+yy2
55  IF (yy2 .LT. xyzmin) xyzmin = yy2
56  IF (ldim .EQ. 3) THEN
57  zz2 = 1./zz**2
58  xyzmax = xyzmax+zz2
59  IF (zz2 .LT. xyzmin) xyzmin = zz2
60  ENDIF
61 C
62  one = 1.
63  pi = 4.*atan(one)
64  ratio = xyzmin/xyzmax
65  eigae = pi*pi*xyzmin
66  eigge = eigga
67  IF (ldim .EQ. 2) eigaa = pi*pi*(xx2+yy2)/2.
68  IF (ldim .EQ. 3) eigaa = pi*pi*(xx2+yy2+zz2)/3.
69  IF (ifaxis) eigaa = .25*pi*pi*yy2
70  eigas = 0.25*ratio
71  eiggs = 2.0
72 C
73  IF (nio.EQ.0 .AND. istep.LE.0) THEN
74  WRITE (6,*) ' '
75  WRITE (6,*) 'Estimated eigenvalues'
76  WRITE (6,*) 'EIGAA = ',eigaa
77  WRITE (6,*) 'EIGGA = ',eigga
78  IF (ifflow) THEN
79  WRITE (6,*) 'EIGAE = ',eigae
80  WRITE (6,*) 'EIGAS = ',eigas
81  WRITE (6,*) 'EIGGE = ',eigge
82  WRITE (6,*) 'EIGGS = ',eiggs
83  ENDIF
84  WRITE (6,*) ' '
85  ENDIF
86 C
87  RETURN
88  END
89 C
90  SUBROUTINE eigenv
91 C-------------------------------------------------------------------------
92 C
93 C Compute the following eigenvalues:
94 C EIGAA = minimum eigenvalue of the matrix A (=Laplacian)
95 C EIGAE = minimum eigenvalue of the matrix E (=DB-1DT)
96 C EIGAS = minimum eigenvalue of the matrix S (=DA-1DT)
97 C EIGAST = minimum eigenvalue of the matrix St (=D(A+B/dt)-1DT
98 C EIGGA = maximum eigenvalue of the matrix A
99 C EIGGS = maximum eigenvalue of the matrix S
100 C EIGGE = maximum eigenvalue of the matrix E
101 C EIGGST = maximum eigenvalue of the matrix St
102 C
103 C Method : Power method/Inverse iteration & Rayleigh quotient wo shift
104 C
105 C-------------------------------------------------------------------------
106  include 'SIZE'
107  include 'EIGEN'
108  include 'INPUT'
109  include 'SOLN'
110  include 'TSTEP'
111 C
112  COMMON /scrvh/ h1(lx1,ly1,lz1,lelt)
113  $ , h2(lx1,ly1,lz1,lelt)
114  COMMON /scrhi/ h2inv(lx1,ly1,lz1,lelv)
115 C
116  ntot1 = lx1*ly1*lz1*nelv
117 C
118  IF (ifaa) THEN
119  ntot1 = lx1*ly1*lz1*nelv
120  CALL rone (h1,ntot1)
121  CALL rzero (h2,ntot1)
122  CALL alpham1 (eigaa1,v1mask,vmult,h1,h2,1)
123  CALL alpham1 (eigaa2,v2mask,vmult,h1,h2,2)
124  eigaa = min(eigaa1,eigaa2)
125  IF (ldim.EQ.3) THEN
126  CALL alpham1 (eigaa3,v3mask,vmult,h1,h2,3)
127  eigaa = min(eigaa,eigaa3)
128  ENDIF
129  IF (nio.EQ.0 .AND. istep.LE.0) WRITE (6,*) 'EIGAA = ',eigaa
130  ENDIF
131 C
132  IF (ifas) THEN
133  inloc = 0
134  CALL rone (h1,ntot1)
135  CALL rzero (h2,ntot1)
136  CALL rzero (h2inv,ntot1)
137  CALL alpham2 (eigas,h1,h2,h2inv,inloc)
138  IF (nio.EQ.0 .AND. istep.LE.0) WRITE (6,*) 'EIGAS = ',eigas
139  ENDIF
140 C
141  IF (ifae) THEN
142  inloc = 1
143  CALL rzero (h1,ntot1)
144  CALL rone (h2,ntot1)
145  CALL rone (h2inv,ntot1)
146  CALL alpham2 (eigae,h1,h2,h2inv,inloc)
147  IF (nio.EQ.0 .AND. istep.LE.0) WRITE (6,*) 'EIGAE = ',eigae
148  ENDIF
149 C
150  IF (ifast) THEN
151  inloc = -1
152  CALL sethlm (h1,h2,inloc)
153  CALL invers2 (h2inv,h2,ntot1)
154  CALL alpham2 (eigast,h1,h2,h2inv,inloc)
155  IF (nio.EQ.0 .AND. istep.LE.0) WRITE (6,*) 'EIGAST = ',eigast
156  ENDIF
157 C
158  IF (ifgs) THEN
159  inloc = 0
160  CALL rone (h1,ntot1)
161  CALL rzero (h2,ntot1)
162  CALL rzero (h2inv,ntot1)
163  CALL gammam2 (eiggs,h1,h2,h2inv,inloc)
164  IF (nio.EQ.0 .AND. istep.LE.0) WRITE (6,*) 'EIGGS = ',eiggs
165  ENDIF
166 C
167  IF (ifge) THEN
168  inloc = 1
169  CALL rzero (h1,ntot1)
170  CALL rone (h2,ntot1)
171  CALL rone (h2inv,ntot1)
172  CALL gammam2 (eigge,h1,h2,h2inv,inloc)
173  IF (nio.EQ.0 .AND. istep.LE.0) WRITE (6,*) 'EIGGE = ',eigge
174  ENDIF
175 C
176  IF (ifgst) THEN
177  inloc = -1
178  CALL sethlm (h1,h2,inloc)
179  CALL invers2 (h2inv,h2,ntot1)
180  CALL gammam2 (eiggst,h1,h2,h2inv,inloc)
181  IF (nio.EQ.0 .AND. istep.LE.0) WRITE (6,*) 'EIGGST = ',eiggst
182  ENDIF
183 C
184  IF (ifga) THEN
185  ntot1 = lx1*ly1*lz1*nelv
186  CALL rone (h1,ntot1)
187  CALL rzero (h2,ntot1)
188  IF (.NOT.ifstrs) THEN
189  CALL gammam1 (eigga1,v1mask,vmult,h1,h2,1)
190  CALL gammam1 (eigga2,v2mask,vmult,h1,h2,2)
191  eigga3 = 0.
192  IF (ldim.EQ.3)
193  $ CALL gammam1 (eigga3,v3mask,vmult,h1,h2,3)
194  eigga = max(eigga1,eigga2,eigga3)
195  ELSE
196  CALL gammasf (h1,h2)
197  ENDIF
198  ENDIF
199 C
200  RETURN
201  END
202 C
203  SUBROUTINE alpham1 (ALPHA,MASK,MULT,H1,H2,ISD)
204 C---------------------------------------------------------------------------
205 C
206 C Compute minimum eigenvalue, ALPHA, of the discrete Helmholtz operator
207 C
208 C---------------------------------------------------------------------------
209  include 'SIZE'
210  include 'MASS'
211  include 'TSTEP'
212 C
213  REAL MASK (LX1,LY1,LZ1,1)
214  REAL MULT (LX1,LY1,LZ1,1)
215  REAL H1 (LX1,LY1,LZ1,1)
216  REAL H2 (LX1,LY1,LZ1,1)
217  COMMON /screv/ x1(lx1,ly1,lz1,lelt)
218  $ , y1(lx1,ly1,lz1,lelt)
219  CHARACTER NAME*4
220 C
221  IF (imesh.EQ.1) nel = nelv
222  IF (imesh.EQ.2) nel = nelt
223  IF (isd .EQ.1) name = 'EVVX'
224  IF (isd .EQ.2) name = 'EVVX'
225  IF (isd .EQ.3) name = 'EVVX'
226 C
227  nxyz1 = lx1*ly1*lz1
228  ntot1 = nxyz1*nel
229  evnew = 0.
230  CALL startx1 (x1,y1,mask,mult,nel)
231 C
232  DO 1000 iter=1,nmxe
233  CALL axhelm (y1,x1,h1,h2,imesh,isd)
234  CALL col2 (y1,mask,ntot1)
235  CALL dssum (y1,lx1,ly1,lz1)
236  rq = glsc3(x1,y1,mult,ntot1)
237  evold = evnew
238  evnew = rq
239  write (6,*) 'alphaa = ',rq
240  crit = abs((evnew-evold)/evnew)
241  IF (crit.LT.tolev) GOTO 2000
242  CALL col2 (x1,bm1,ntot1)
243  CALL hmholtz ('NOMG',y1,x1,h1,h2,mask,mult,
244  $ imesh,tolhe,nmxe,isd)
245  CALL col3 (x1,bm1,y1,ntot1)
246  CALL dssum (x1,lx1,ly1,lz1)
247  yy = glsc3(x1,y1,mult,ntot1)
248  ynorm = 1./sqrt(yy)
249  CALL cmult (y1,ynorm,ntot1)
250  CALL copy (x1,y1,ntot1)
251  1000 CONTINUE
252  2000 CONTINUE
253 C
254  alpha = rq
255  RETURN
256  END
257 C
258  SUBROUTINE gammam1 (GAMMA,MASK,MULT,H1,H2,ISD)
259 C---------------------------------------------------------------------------
260 C
261 C Compute maximum eigenvalue of the discrete Helmholtz operator
262 C
263 C---------------------------------------------------------------------------
264  include 'SIZE'
265  include 'MASS'
266  include 'TSTEP'
267 C
268  REAL MASK (LX1,LY1,LZ1,1)
269  REAL MULT (LX1,LY1,LZ1,1)
270  REAL H1 (LX1,LY1,LZ1,1)
271  REAL H2 (LX1,LY1,LZ1,1)
272  COMMON /screv/ x1(lx1,ly1,lz1,lelt)
273  $ , y1(lx1,ly1,lz1,lelt)
274 C
275  IF (imesh.EQ.1) nel = nelv
276  IF (imesh.EQ.2) nel = nelt
277  nxyz1 = lx1*ly1*lz1
278  ntot1 = nxyz1*nel
279  evnew = 0.
280 c pff (2/15/96)
281  if (isd.eq.1) CALL startx1 (x1,y1,mask,mult,nel)
282 C
283  DO 1000 iter=1,nmxe
284  CALL axhelm (y1,x1,h1,h2,imesh,isd)
285  CALL col2 (y1,mask,ntot1)
286  CALL dssum (y1,lx1,ly1,lz1)
287  rq = glsc3(x1,y1,mult,ntot1)
288  evold = evnew
289  evnew = rq
290  crit = abs((evnew-evold)/evnew)
291 C
292 C HMT removed
293 C
294 C if (nio.eq.0) then
295 C write(6,*) iter,' eig_max A:',evnew,crit,tolev
296 C endif
297  IF (crit.LT.tolev) GOTO 2000
298  CALL col3 (x1,binvm1,y1,ntot1)
299  xx = glsc3(x1,y1,mult,ntot1)
300  xnorm = 1./sqrt(xx)
301  CALL cmult (x1,xnorm,ntot1)
302  1000 CONTINUE
303  2000 CONTINUE
304 C
305  gamma = rq
306  RETURN
307  END
308 C
309  SUBROUTINE alpham2 (ALPHA,H1,H2,H2INV,INLOC)
310 C----------------------------------------------------------------------
311 C
312 C Compute minimum eigenvalue, ALPHA, of one of the matrices
313 C defined on the pressure mesh:
314 C INLOC = 0 : DA-1DT
315 C INLOC = 1 : DB-1DT
316 C INLOC = -1 : D(A+B/DT)-1DT
317 C
318 C----------------------------------------------------------------------
319  include 'SIZE'
320  include 'MASS'
321  include 'TSTEP'
322 C
323  REAL H1 (LX1,LY1,LZ1,1)
324  REAL H2 (LX1,LY1,LZ1,1)
325  REAL H2INV(LX1,LY1,LZ1,1)
326  COMMON /screv/ x2(lx2,ly2,lz2,lelv)
327  $ , y2(lx2,ly2,lz2,lelv)
328 C
329  ntot2 = lx2*ly2*lz2*nelv
330  evnew = 0.
331  CALL startx2 (x2,y2)
332 C
333  DO 1000 iter=1,nmxe
334  CALL cdabdtp (y2,x2,h1,h2,h2inv,inloc)
335  rq = glsc2(x2,y2,ntot2)
336  evold = evnew
337  evnew = rq
338 c write (6,*) 'new eigenvalue ************* eigas = ',evnew
339  crit = abs((evnew-evold)/evnew)
340  IF (crit.LT.tolev) GOTO 2000
341  CALL col2 (x2,bm2,ntot2)
342  CALL uzawa (x2,h1,h2,h2inv,inloc,icg)
343  CALL col3 (y2,bm2,x2,ntot2)
344  xx = glsc2(x2,y2,ntot2)
345  xnorm = 1./sqrt(xx)
346  CALL cmult (x2,xnorm,ntot2)
347  1000 CONTINUE
348  2000 CONTINUE
349 C
350  alpha = rq
351  RETURN
352  END
353 C
354  SUBROUTINE gammam2 (GAMMA,H1,H2,H2INV,INLOC)
355 C-------------------------------------------------------------------
356 C
357 C Compute maximum eigenvalue, GAMMA, of one of the matrices
358 C defined on the pressure mesh:
359 C INLOC = 0 : DA-1DT
360 C INLOC = 1 : DB-1DT
361 C INLOC = -1 : D(A+B/DT)-1DT
362 C
363 C-------------------------------------------------------------------
364  include 'SIZE'
365  include 'MASS'
366  include 'TSTEP'
367 C
368  REAL H1 (LX1,LY1,LZ1,1)
369  REAL H2 (LX1,LY1,LZ1,1)
370  REAL H2INV (LX1,LY1,LZ1,1)
371  COMMON /screv/ x2(lx2,ly2,lz2,lelv)
372  $ , y2(lx2,ly2,lz2,lelv)
373 C
374  ntot2 = lx2*ly2*lz2*nelv
375  evnew = 0.
376  CALL startx2 (x2,y2)
377 C
378  DO 1000 iter=1,nmxe
379  CALL cdabdtp (y2,x2,h1,h2,h2inv,inloc)
380  rq = glsc2(x2,y2,ntot2)
381  evold = evnew
382  evnew = rq
383  crit = abs((evnew-evold)/evnew)
384  IF (crit.LT.tolev) GOTO 2000
385  CALL invcol3 (x2,y2,bm2,ntot2)
386  xx = glsc2(y2,x2,ntot2)
387  xnorm = 1./sqrt(xx)
388  CALL cmult (x2,xnorm,ntot2)
389  1000 CONTINUE
390  2000 CONTINUE
391 C
392  gamma = rq
393  RETURN
394  END
395 C
396 c-----------------------------------------------------------------------
397  SUBROUTINE startx1 (X1,Y1,MASK,MULT,NEL)
398 
399 c Compute startvector for finding an eigenvalue on mesh 1.
400 c Normalization: XT*B*X = 1
401 
402  include 'SIZE'
403  include 'MASS'
404 
405  REAL X1 (LX1,LY1,LZ1,1)
406  REAL Y1 (LX1,LY1,LZ1,1)
407  REAL MASK (LX1,LY1,LZ1,1)
408  REAL MULT (LX1,LY1,LZ1,1)
409 
410  ntot1 = lx1*ly1*lz1*nel
411  CALL copy (x1,bm1,ntot1)
412 
413 
414  call rand_fld_h1(y1) ! pff 3/21/12
415  small = 0.001*glamax(x1,ntot1)
416  call add2s2(x1,y1,small,ntot1)
417 
418 
419  CALL col2 (x1,mask,ntot1)
420  CALL col3 (y1,bm1,x1,ntot1)
421  CALL dssum (y1,lx1,ly1,lz1)
422  xx = glsc3(x1,y1,mult,ntot1)
423  xnorm = 1./sqrt(xx)
424  CALL cmult (x1,xnorm,ntot1)
425 
426  RETURN
427  END
428 c-----------------------------------------------------------------------
429  SUBROUTINE startx2 (X2,Y2)
430 C------------------------------------------------------------------
431 C
432 C Compute startvector for finding an eigenvalue on mesh 2.
433 C
434 C------------------------------------------------------------------
435  include 'SIZE'
436  include 'MASS'
437 C
438  REAL X2 (LX2,LY2,LZ2,LELV)
439  REAL Y2 (LX2,LY2,LZ2,LELV)
440 C
441  nxyz2 = lx2*ly2*lz2
442  ntot2 = nxyz2*nelv
443  iconst = 0
444  IF ((ldim .EQ. 2).AND.(nxyz2 .EQ. 4)) iconst = 1
445  IF ((ldim .EQ. 3).AND.(nxyz2 .EQ. 8)) iconst = 1
446 C
447  IF (iconst .EQ. 1) THEN
448  DO 1000 iel=1,nelv
449  DO 1000 k=1,lz2
450  DO 1000 j=1,ly2
451  DO 1000 i=1,lx2
452  x2(i,j,k,iel) = i*j*k
453  1000 CONTINUE
454  ELSE
455  CALL copy (x2,bm2,ntot2)
456  ENDIF
457 C
458  call ortho (x2)
459  CALL col3 (y2,bm2,x2,ntot2)
460  xx = glsc2(x2,y2,ntot2)
461  xnorm = 1./sqrt(xx)
462  CALL cmult (x2,xnorm,ntot2)
463 C
464  RETURN
465  END
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
esteig
subroutine esteig
Definition: eigsolv.f:9
alpham1
subroutine alpham1(ALPHA, MASK, MULT, H1, H2, ISD)
Definition: eigsolv.f:204
gammam1
subroutine gammam1(GAMMA, MASK, MULT, H1, H2, ISD)
Definition: eigsolv.f:259
startx1
subroutine startx1(X1, Y1, MASK, MULT, NEL)
Definition: eigsolv.f:398
startx2
subroutine startx2(X2, Y2)
Definition: eigsolv.f:430
eigenv
subroutine eigenv
Definition: eigsolv.f:91
gammam2
subroutine gammam2(GAMMA, H1, H2, H2INV, INLOC)
Definition: eigsolv.f:355
alpham2
subroutine alpham2(ALPHA, H1, H2, H2INV, INLOC)
Definition: eigsolv.f:310
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.f:3
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
glmin
function glmin(a, n)
Definition: math.f:973
glsc2
function glsc2(x, y, n)
Definition: math.f:794
glsc3
function glsc3(a, b, mult, n)
Definition: math.f:776
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
copy
subroutine copy(a, b, n)
Definition: math.f:260
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
glamax
real function glamax(a, n)
Definition: math.f:874
rone
subroutine rone(a, n)
Definition: math.f:230
invcol3
subroutine invcol3(a, b, c, n)
Definition: math.f:98
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
uzawa
subroutine uzawa(rcg, h1, h2, h2inv, intype, iter)
Definition: navier1.f:2827
ortho
subroutine ortho(respr)
Definition: navier1.f:224
cdabdtp
subroutine cdabdtp(ap, wp, h1, h2, h2inv, intype)
Definition: navier1.f:259
rand_fld_h1
subroutine rand_fld_h1(x)
Definition: navier5.f:2687
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
gammasf
subroutine gammasf(H1, H2)
Definition: subs2.f:256