Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
eigsolv.F90
Go to the documentation of this file.
1 !**********************************************************************
4 !**********************************************************************
5 
6 !--------------------------------------------------------------
8 !-------------------------------------------------------------
9 SUBROUTINE esteig
10  use kinds, only : dp
11  use size_m, only : nx1, ny1, nz1, ndim, nid
12  use eigen, only : eigae, eigge, eigga, eigaa, eigas, eiggs
13  use geom, only : xm1, ym1, zm1
14  use input, only : if3d, ifaxis, ifflow
15  use tstep, only : nelfld, ifield, pi, istep
16  implicit none
17 
18  integer :: ntot1
19  real(DP) :: xmin, xmax, ymin, ymax, zmin, zmax
20  real(DP) :: xx, yy, zz, rxy, ryx, rxz, rzx, ryz, rzy, rmin
21  real(DP) :: xx2, yy2, zz2, xyzmin, xyzmax
22  real(DP) :: one, ratio
23  real(DP), external :: glmin, glmax
24 
25  ntot1=nx1*ny1*nz1*nelfld(ifield)
26  xmin = glmin(xm1,ntot1)
27  xmax = glmax(xm1,ntot1)
28  ymin = glmin(ym1,ntot1)
29  ymax = glmax(ym1,ntot1)
30  IF (if3d) THEN
31  zmin = glmin(zm1,ntot1)
32  zmax = glmax(zm1,ntot1)
33  ELSE
34  zmin = 0.0
35  zmax = 0.0
36  ENDIF
37 
38  xx = xmax - xmin
39  yy = ymax - ymin
40  zz = zmax - zmin
41  rxy = xx/yy
42  ryx = yy/xx
43  rmin = rxy
44  IF (ryx < rmin) rmin = ryx
45  IF (ndim == 3) THEN
46  rxz = xx/zz
47  rzx = zz/xx
48  ryz = yy/zz
49  rzy = zz/yy
50  IF (rxz < rmin) rmin = rxz
51  IF (rzx < rmin) rmin = rzx
52  IF (ryz < rmin) rmin = ryz
53  IF (rzy < rmin) rmin = rzy
54  ENDIF
55 
56  xx2 = 1./xx**2
57  yy2 = 1./yy**2
58  xyzmin = xx2
59  xyzmax = xx2+yy2
60  IF (yy2 < xyzmin) xyzmin = yy2
61  IF (ndim == 3) THEN
62  zz2 = 1./zz**2
63  xyzmax = xyzmax+zz2
64  IF (zz2 < xyzmin) xyzmin = zz2
65  ENDIF
66 
67  one = 1.
68  pi = 4.*atan(one)
69  ratio = xyzmin/xyzmax
70  eigae = pi*pi*xyzmin
71  eigge = eigga
72  IF (ndim == 2) eigaa = pi*pi*(xx2+yy2)/2.
73  IF (ndim == 3) eigaa = pi*pi*(xx2+yy2+zz2)/3.
74  IF (ifaxis) eigaa = .25*pi*pi*yy2
75  eigas = 0.25*ratio
76  eiggs = 2.0
77 
78  IF (nid == 0 .AND. istep <= 0) THEN
79  WRITE (6,*) ' '
80  WRITE (6,*) 'Estimated eigenvalues'
81  WRITE (6,*) 'EIGAA = ',eigaa
82  WRITE (6,*) 'EIGGA = ',eigga
83  IF (ifflow) THEN
84  WRITE (6,*) 'EIGAE = ',eigae
85  WRITE (6,*) 'EIGAS = ',eigas
86  WRITE (6,*) 'EIGGE = ',eigge
87  WRITE (6,*) 'EIGGS = ',eiggs
88  ENDIF
89  WRITE (6,*) ' '
90  ENDIF
91 
92  RETURN
93 END SUBROUTINE esteig
94 
95 !-------------------------------------------------------------------------
106 !-------------------------------------------------------------------------
107 SUBROUTINE eigenv
108  use kinds, only : dp
109  use size_m, only : lx1, ly1, lz1, lelt, nx1, ny1, nz1, nelv, ndim
110  use eigen, only : ifaa, ifas, ifae, ifast, ifgs, ifge, ifgst, ifga, eigga
111  use input, only : ifstrs
112  use soln, only : vmult, v1mask, v2mask, v3mask
113  implicit none
114 
115  real(DP), allocatable :: H1(:,:,:,:), H2(:,:,:,:)
116 ! C!OMMON /SCRHI/ H2INV (LX1,LY1,LZ1,LELV)
117  integer :: ntot1
118  real(DP) :: eigga1, eigga2, eigga3
119 
120  allocate(h1(lx1,ly1,lz1,lelt), h2(lx1,ly1,lz1,lelt))
121 
122  ntot1 = nx1*ny1*nz1*nelv
123 
124  IF (ifaa) THEN
125  write(*,*) "Oops: IFAA"
126 #if 0
127  ntot1 = nx1*ny1*nz1*nelv
128  CALL rone(h1,ntot1)
129  CALL rzero(h2,ntot1)
130  CALL alpham1(eigaa1,v1mask,vmult,h1,h2,1)
131  CALL alpham1(eigaa2,v2mask,vmult,h1,h2,2)
132  eigaa = min(eigaa1,eigaa2)
133  IF (ndim == 3) THEN
134  CALL alpham1(eigaa3,v3mask,vmult,h1,h2,3)
135  eigaa = min(eigaa,eigaa3)
136  ENDIF
137  IF (nid == 0 .AND. istep <= 0) WRITE (6,*) 'EIGAA = ',eigaa
138 #endif
139  ENDIF
140 
141  IF (ifas) THEN
142  write(*,*) "Oops: IFAA"
143 #if 0
144  inloc = 0
145  CALL rone(h1,ntot1)
146  CALL rzero(h2,ntot1)
147  CALL rzero(h2inv,ntot1)
148  CALL alpham2(eigas,h1,h2,h2inv,inloc)
149  IF (nid == 0 .AND. istep <= 0) WRITE (6,*) 'EIGAS = ',eigas
150 #endif
151  ENDIF
152 
153  IF (ifae) THEN
154  write(*,*) "Oops: IFAE"
155 #if 0
156  inloc = 1
157  CALL rzero(h1,ntot1)
158  CALL rone(h2,ntot1)
159  CALL rone(h2inv,ntot1)
160  CALL alpham2(eigae,h1,h2,h2inv,inloc)
161  IF (nid == 0 .AND. istep <= 0) WRITE (6,*) 'EIGAE = ',eigae
162 #endif
163  ENDIF
164 
165  IF (ifast) THEN
166  write(*,*) "Oops: IFAST"
167 #if 0
168  inloc = -1
169  CALL sethlm(h1,h2,inloc)
170  CALL invers2(h2inv,h2,ntot1)
171  CALL alpham2(eigast,h1,h2,h2inv,inloc)
172  IF (nid == 0 .AND. istep <= 0) WRITE (6,*) 'EIGAST = ',eigast
173 #endif
174  ENDIF
175 
176  IF (ifgs) THEN
177  write(*,*) "Oops: IFGS"
178 #if 0
179  inloc = 0
180  CALL rone(h1,ntot1)
181  CALL rzero(h2,ntot1)
182  CALL rzero(h2inv,ntot1)
183  CALL gammam2(eiggs,h1,h2,h2inv,inloc)
184  IF (nid == 0 .AND. istep <= 0) WRITE (6,*) 'EIGGS = ',eiggs
185 #endif
186  ENDIF
187 
188  IF (ifge) THEN
189  write(*,*) "Oops: IFGE"
190 #if 0
191  inloc = 1
192  CALL rzero(h1,ntot1)
193  CALL rone(h2,ntot1)
194  CALL rone(h2inv,ntot1)
195  CALL gammam2(eigge,h1,h2,h2inv,inloc)
196  IF (nid == 0 .AND. istep <= 0) WRITE (6,*) 'EIGGE = ',eigge
197 #endif
198  ENDIF
199 
200  IF (ifgst) THEN
201  write(*,*) "Oops: IFGST"
202 #if 0
203  inloc = -1
204  CALL sethlm(h1,h2,inloc)
205  CALL invers2(h2inv,h2,ntot1)
206  CALL gammam2(eiggst,h1,h2,h2inv,inloc)
207  IF (nid == 0 .AND. istep <= 0) WRITE (6,*) 'EIGGST = ',eiggst
208 #endif
209  ENDIF
210 
211  IF (ifga) THEN
212  ntot1 = nx1*ny1*nz1*nelv
213  h1 = 1._dp; h2 = 0._dp
214  IF ( .NOT. ifstrs) THEN
215  CALL gammam1(eigga1,v1mask,vmult,h1,h2,1)
216  CALL gammam1(eigga2,v2mask,vmult,h1,h2,2)
217  eigga3 = 0.
218  IF (ndim == 3) &
219  CALL gammam1(eigga3,v3mask,vmult,h1,h2,3)
220  eigga = max(eigga1,eigga2,eigga3)
221  ELSE
222 !max CALL GAMMASF (H1,H2)
223  ENDIF
224  ENDIF
225 
226  RETURN
227 END SUBROUTINE eigenv
228 
229 !---------------------------------------------------------------------------
231 !---------------------------------------------------------------------------
232 SUBROUTINE gammam1 (GAMMA,MASK,MULT,H1,H2,ISD)
233  use kinds, only : dp
234  use size_m, only : lx1, ly1, lz1, lelt
235  use size_m, only : nx1, ny1, nz1, nelt, nelv
236  use geom, only : binvm1
237  use tstep, only : imesh, nmxe, tolev
238  implicit none
239 
240  real(DP) :: gamma
241  REAL(DP) :: MASK (lx1,ly1,lz1,*)
242  REAL(DP) :: MULT (lx1,ly1,lz1,*)
243  REAL(DP) :: H1 (lx1,ly1,lz1,*)
244  REAL(DP) :: H2 (lx1,ly1,lz1,*)
245  integer :: isd
246 
247  real(DP), allocatable :: X1(:,:,:,:), Y1(:,:,:,:)
248 
249  integer :: nel, nxyz1, ntot1
250  integer :: iter
251  real(DP) :: evnew, rq, evold, crit, xx, xnorm
252  real(DP), external :: glsc3
253 
254  allocate(x1(lx1,ly1,lz1,lelt), y1(lx1,ly1,lz1,lelt))
255 
256 
257  IF (imesh == 1) nel = nelv
258  IF (imesh == 2) nel = nelt
259  nxyz1 = nx1*ny1*nz1
260  ntot1 = nxyz1*nel
261  evnew = 0.
262 ! pff (2/15/96)
263  if (isd == 1) CALL startx1(x1,y1,mask,mult,nel)
264 
265  rq = 0._dp
266  DO 1000 iter=1,nmxe
267  CALL axhelm(y1,x1,h1,h2,imesh,isd)
268  y1 = y1 * mask(:,:,:,1:nel)
269  CALL dssum(y1)
270  rq = glsc3(x1,y1,mult,ntot1)
271  evold = evnew
272  evnew = rq
273  crit = abs((evnew-evold)/evnew)
274 
275  ! HMT removed
276 
277  ! if (nid.eq.0) then
278  ! write(6,*) iter,' eig_max A:',evnew,crit,tolev
279  ! endif
280  IF (crit < tolev) goto 2000
281  x1 = binvm1 * y1
282  xx = glsc3(x1,y1,mult,ntot1)
283  xnorm = 1./sqrt(xx)
284  x1 = x1 * xnorm
285  1000 END DO
286  2000 CONTINUE
287 
288  gamma = rq
289  RETURN
290 END SUBROUTINE gammam1
291 
292 !-----------------------------------------------------------------------
295 SUBROUTINE startx1 (X1,Y1,MASK,MULT,NEL)
296  use kinds, only : dp
297  use size_m, only : lx1, ly1, lz1, nx1, ny1, nz1
298  use geom, only : bm1
299  implicit none
300 
301  integer :: nel
302  REAL(DP), intent(out) :: X1 (lx1,ly1,lz1,nel)
303  REAL(DP) :: Y1 (lx1,ly1,lz1,nel)
304  REAL(DP), intent(in) :: MASK (lx1,ly1,lz1,nel)
305  REAL(DP), intent(in) :: MULT (lx1,ly1,lz1,nel)
306 
307  integer :: ntot1
308  real(DP) :: small, xx, xnorm
309  real(DP), external :: glamax, glsc3
310 
311  ntot1 = nx1*ny1*nz1*nel
312  CALL copy(x1,bm1,ntot1)
313 
314 
315  call rand_fld_h1(y1) ! pff 3/21/12
316  small = 0.001*glamax(x1,ntot1)
317  x1 = x1 + small * y1
318 
319 
320  x1 = x1 * mask
321  y1 = bm1 * x1
322  CALL dssum(y1)
323  xx = glsc3(x1,y1,mult,ntot1)
324  xnorm = 1./sqrt(xx)
325  x1 = x1 * xnorm
326 
327  RETURN
328 END SUBROUTINE startx1
329 
330 !-----------------------------------------------------------------------
dssum
subroutine dssum(u)
Direct stiffness sum.
Definition: dssum.F90:54
tstep
cleaned
Definition: tstep_mod.F90:2
input
Input parameters from preprocessors.
Definition: input_mod.F90:11
soln
Definition: soln_mod.F90:1
esteig
subroutine esteig
Estimate eigenvalues.
Definition: eigsolv.F90:9
sethlm
subroutine sethlm(h1, h2, intloc)
Set the variable property arrays H1 and H2 in the Helmholtz equation. (associated with variable IFIEL...
Definition: subs1.F90:694
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Compute the (Helmholtz) matrix-vector product, AU = helm1*[A]u + helm2*[B]u, for NEL elements...
Definition: hmholtz.F90:79
startx1
subroutine startx1(X1, Y1, MASK, MULT, NEL)
Compute startvector for finding an eigenvalue on mesh 1. Normalization: XT*B*X = 1.
Definition: eigsolv.F90:295
copy
subroutine copy(a, b, n)
Definition: math.F90:52
eigenv
subroutine eigenv
Compute the following eigenvalues:. EIGAA = minimum eigenvalue of the matrix A (=Laplacian) EIGAE = m...
Definition: eigsolv.F90:107
gammam1
subroutine gammam1(GAMMA, MASK, MULT, H1, H2, ISD)
Compute maximum eigenvalue of the discrete Helmholtz operator.
Definition: eigsolv.F90:232
geom
Geometry arrays.
Definition: geom_mod.F90:2
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8