source: CIVL/examples/fortran/nek5000/core/eigsolv.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 13.3 KB
Line 
1C**********************************************************************
2C
3C ROUTINES FOR ESITMATING AND CALCULATING EIGENVALUES
4C USED IN NEKTON
5C
6C**********************************************************************
7C
8 SUBROUTINE ESTEIG
9C--------------------------------------------------------------
10C
11C Estimate eigenvalues
12C
13C-------------------------------------------------------------
14 INCLUDE 'SIZE'
15 INCLUDE 'GEOM'
16 INCLUDE 'INPUT'
17 INCLUDE 'EIGEN'
18 INCLUDE 'TSTEP'
19C
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
32C
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
50C
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
61C
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
72C
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
86C
87 RETURN
88 END
89C
90 SUBROUTINE EIGENV
91C-------------------------------------------------------------------------
92C
93C Compute the following eigenvalues:
94C EIGAA = minimum eigenvalue of the matrix A (=Laplacian)
95C EIGAE = minimum eigenvalue of the matrix E (=DB-1DT)
96C EIGAS = minimum eigenvalue of the matrix S (=DA-1DT)
97C EIGAST = minimum eigenvalue of the matrix St (=D(A+B/dt)-1DT
98C EIGGA = maximum eigenvalue of the matrix A
99C EIGGS = maximum eigenvalue of the matrix S
100C EIGGE = maximum eigenvalue of the matrix E
101C EIGGST = maximum eigenvalue of the matrix St
102C
103C Method : Power method/Inverse iteration & Rayleigh quotient wo shift
104C
105C-------------------------------------------------------------------------
106 INCLUDE 'SIZE'
107 INCLUDE 'EIGEN'
108 INCLUDE 'INPUT'
109 INCLUDE 'SOLN'
110 INCLUDE 'TSTEP'
111C
112 COMMON /SCRVH/ H1 (LX1,LY1,LZ1,LELT)
113 $ , H2 (LX1,LY1,LZ1,LELT)
114 COMMON /SCRHI/ H2INV (LX1,LY1,LZ1,LELV)
115C
116 NTOT1 = lx1*ly1*lz1*NELV
117C
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
131C
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
140C
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
149C
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
157C
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
166C
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
175C
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
183C
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
199C
200 RETURN
201 END
202C
203 SUBROUTINE ALPHAM1 (ALPHA,MASK,MULT,H1,H2,ISD)
204C---------------------------------------------------------------------------
205C
206C Compute minimum eigenvalue, ALPHA, of the discrete Helmholtz operator
207C
208C---------------------------------------------------------------------------
209 INCLUDE 'SIZE'
210 INCLUDE 'MASS'
211 INCLUDE 'TSTEP'
212C
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
220C
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'
226C
227 NXYZ1 = lx1*ly1*lz1
228 NTOT1 = NXYZ1*NEL
229 EVNEW = 0.
230 CALL STARTX1 (X1,Y1,MASK,MULT,NEL)
231C
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
253C
254 ALPHA = RQ
255 RETURN
256 END
257C
258 SUBROUTINE GAMMAM1 (GAMMA,MASK,MULT,H1,H2,ISD)
259C---------------------------------------------------------------------------
260C
261C Compute maximum eigenvalue of the discrete Helmholtz operator
262C
263C---------------------------------------------------------------------------
264 INCLUDE 'SIZE'
265 INCLUDE 'MASS'
266 INCLUDE 'TSTEP'
267C
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)
274C
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.
280c pff (2/15/96)
281 if (isd.eq.1) CALL STARTX1 (X1,Y1,MASK,MULT,NEL)
282C
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)
291C
292C HMT removed
293C
294C if (nio.eq.0) then
295C write(6,*) iter,' eig_max A:',evnew,crit,tolev
296C 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
304C
305 GAMMA = RQ
306 RETURN
307 END
308C
309 SUBROUTINE ALPHAM2 (ALPHA,H1,H2,H2INV,INLOC)
310C----------------------------------------------------------------------
311C
312C Compute minimum eigenvalue, ALPHA, of one of the matrices
313C defined on the pressure mesh:
314C INLOC = 0 : DA-1DT
315C INLOC = 1 : DB-1DT
316C INLOC = -1 : D(A+B/DT)-1DT
317C
318C----------------------------------------------------------------------
319 INCLUDE 'SIZE'
320 INCLUDE 'MASS'
321 INCLUDE 'TSTEP'
322C
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)
328C
329 NTOT2 = lx2*ly2*lz2*NELV
330 EVNEW = 0.
331 CALL STARTX2 (X2,Y2)
332C
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
338c 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
349C
350 ALPHA = RQ
351 RETURN
352 END
353C
354 SUBROUTINE GAMMAM2 (GAMMA,H1,H2,H2INV,INLOC)
355C-------------------------------------------------------------------
356C
357C Compute maximum eigenvalue, GAMMA, of one of the matrices
358C defined on the pressure mesh:
359C INLOC = 0 : DA-1DT
360C INLOC = 1 : DB-1DT
361C INLOC = -1 : D(A+B/DT)-1DT
362C
363C-------------------------------------------------------------------
364 INCLUDE 'SIZE'
365 INCLUDE 'MASS'
366 INCLUDE 'TSTEP'
367C
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)
373C
374 NTOT2 = lx2*ly2*lz2*NELV
375 EVNEW = 0.
376 CALL STARTX2 (X2,Y2)
377C
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
391C
392 GAMMA = RQ
393 RETURN
394 END
395C
396c-----------------------------------------------------------------------
397 SUBROUTINE STARTX1 (X1,Y1,MASK,MULT,NEL)
398
399c Compute startvector for finding an eigenvalue on mesh 1.
400c 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
428c-----------------------------------------------------------------------
429 SUBROUTINE STARTX2 (X2,Y2)
430C------------------------------------------------------------------
431C
432C Compute startvector for finding an eigenvalue on mesh 2.
433C
434C------------------------------------------------------------------
435 INCLUDE 'SIZE'
436 INCLUDE 'MASS'
437C
438 REAL X2 (LX2,LY2,LZ2,LELV)
439 REAL Y2 (LX2,LY2,LZ2,LELV)
440C
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
446C
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
457C
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)
463C
464 RETURN
465 END
Note: See TracBrowser for help on using the repository browser.