KTH framework for Nek5000 toolboxes; testing version  0.0.1
navier2.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine aspect_ratios(ar)
3 c
4 c 7/6/96
5 c
6 c This routine returns the aspect ratio of a
7 c conglomerate of a set of simplices defined by (tri,nt)
8 c
9  include 'SIZE'
10  include 'GEOM'
11  include 'INPUT'
12  include 'MASS'
13 c
14  real ar(1)
15  real xx(9)
16 c
17  nxyz = lx1*ly1*lz1
18 c
19  if (if3d) then
20  do ie=1,nelv
21  vol = vlsum(bm1(1,1,1,ie),nxyz)
22  x0 = vlsc2(xm1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
23  y0 = vlsc2(ym1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
24  z0 = vlsc2(zm1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
25 c
26  call rzero(xx,9)
27  do i=1,nxyz
28  x10 = xm1(i,1,1,ie) - x0
29  y10 = ym1(i,1,1,ie) - y0
30  z10 = zm1(i,1,1,ie) - z0
31 c
32  xx(1) = xx(1) + bm1(i,1,1,ie)*x10*x10
33  xx(2) = xx(2) + bm1(i,1,1,ie)*x10*y10
34  xx(3) = xx(3) + bm1(i,1,1,ie)*x10*z10
35  xx(4) = xx(4) + bm1(i,1,1,ie)*y10*x10
36  xx(5) = xx(5) + bm1(i,1,1,ie)*y10*y10
37  xx(6) = xx(6) + bm1(i,1,1,ie)*y10*z10
38  xx(7) = xx(7) + bm1(i,1,1,ie)*z10*x10
39  xx(8) = xx(8) + bm1(i,1,1,ie)*z10*y10
40  xx(9) = xx(9) + bm1(i,1,1,ie)*z10*z10
41  enddo
42  vi = 1./vol
43  call cmult(xx,vi,9)
44 c call eig3(xx,eign,eig1)
45  call eig2(xx,eign,eig1)
46  ar(ie) = sqrt(eign/eig1)
47  enddo
48  else
49  do ie=1,nelv
50  vol = vlsum(bm1(1,1,1,ie),nxyz)
51  x0 = vlsc2(xm1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
52  y0 = vlsc2(ym1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
53 c
54  call rzero(xx,4)
55  do i=1,nxyz
56  x10 = xm1(i,1,1,ie) - x0
57  y10 = ym1(i,1,1,ie) - y0
58  xx(1) = xx(1) + bm1(i,1,1,ie)*x10*x10
59  xx(2) = xx(2) + bm1(i,1,1,ie)*x10*y10
60  xx(3) = xx(3) + bm1(i,1,1,ie)*y10*x10
61  xx(4) = xx(4) + bm1(i,1,1,ie)*y10*y10
62  enddo
63  vi = 1./vol
64  call cmult(xx,vi,4)
65  call eig2(xx,eign,eig1)
66 c write(6,6) ie,vol,eign,eig1
67 c 6 format(i5,' veig:',1p3e16.6)
68  ar(ie) = sqrt(eign/eig1)
69  enddo
70  endif
71 c
72 c do ie=1,nelv
73 c write(6,*) ar(ie),ie,' aspect'
74 c enddo
75 c
76  return
77  end
78 c-----------------------------------------------------------------------
79  subroutine eig2(AA,eign,eig1)
80 c
81 c return max and min eigenvalues of a 2x2 matrix
82 c
83  real aa(2,2)
84 c
85  a = aa(1,1)
86  b = aa(1,2)
87  c = aa(2,1)
88  d = aa(2,2)
89 c
90  qa = 1.
91  qb = -(a+d)
92  qc = (a*d - c*b)
93  call quadratic_h(eign,eig1,qa,qb,qc)
94 c
95  return
96  end
97 c-----------------------------------------------------------------------
98  subroutine quadratic_h(x1,x2,a,b,c)
99 c
100 c Stable routine for computation of real roots of quadratic
101 c
102 c The "_h" denotes the hack below so we don't need to worry
103 c about complex arithmetic in the near-double root case. pff 1/22/97
104 c
105 c Upon return, | x1 | >= | x2 |
106 c
107  x1 = 0.
108  x2 = 0.
109 c
110  if (a.eq.0.) then
111  if (b.eq.0) then
112  write(6,10) x1,x2,a,b,c
113  return
114  endif
115  x1 = -c/b
116  write(6,11) x1,a,b,c
117  return
118  endif
119 c
120  d = b*b - 4.*a*c
121  if (d.lt.0) then
122 c write(6,12) a,b,c,d
123 c hack, for this application we'll assume d<0 by epsilon, and just set
124  d = 0
125  endif
126  if (d.gt.0) d = sqrt(d)
127 c
128  q = -0.5 * (b + b/abs(b) * d)
129  x1 = q/a
130  x2 = c/q
131 c
132  if (abs(x2).gt.abs(x1)) then
133  tmp = x2
134  x2 = x1
135  x1 = tmp
136  endif
137 c
138  10 format('ERROR: Both a & b zero in routine quadratic NO ROOTS.'
139  $ ,1p5e12.4)
140  11 format('ERROR: a = 0 in routine quadratic, only one root.'
141  $ ,1p5e12.4)
142  12 format('ERROR: negative discriminate in routine quadratic.'
143  $ ,1p5e12.4)
144 c
145  return
146  end
147 c-----------------------------------------------------------------------
148  subroutine out_sem(iel)
149 c
150  include 'SIZE'
151  include 'INPUT'
152  include 'GEOM'
153 c
154 c
155  open(unit=33,file='v1')
156  if (iel.eq.0) then
157  do ie=1,nelv
158  write(33,33) xm1( 1, 1,1,ie),ym1( 1, 1,1,ie)
159  write(33,33) xm1(lx1, 1,1,ie),ym1(lx1, 1,1,ie)
160  write(33,33) xm1(lx1,ly1,1,ie),ym1(lx1,ly1,1,ie)
161  write(33,33) xm1( 1,ly1,1,ie),ym1( 1,ly1,1,ie)
162  enddo
163  else
164  ie = iel
165  write(33,33) xm1( 1, 1,1,ie),ym1( 1, 1,1,ie)
166  write(33,33) xm1(lx1, 1,1,ie),ym1(lx1, 1,1,ie)
167  write(33,33) xm1(lx1,ly1,1,ie),ym1(lx1,ly1,1,ie)
168  write(33,33) xm1( 1,ly1,1,ie),ym1( 1,ly1,1,ie)
169  endif
170  33 format(f14.6)
171  close(unit=33)
172  return
173  end
174 c-----------------------------------------------------------------------
175  subroutine gradm11(ux,uy,uz,u,e)
176 c
177 c Compute gradient of T -- mesh 1 to mesh 1 (vel. to vel.)
178 c
179 c Single element case
180 c
181  include 'SIZE'
182  include 'DXYZ'
183  include 'GEOM'
184  include 'INPUT'
185  include 'TSTEP'
186 c
187  parameter(lxyz=lx1*ly1*lz1)
188  real ux(lxyz),uy(lxyz),uz(lxyz),u(lxyz,1)
189 c
190  common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz)
191 c
192  integer e
193 
194 c
195  n = lx1-1
196  if (if3d) then
197  call local_grad3(ur,us,ut,u,n,e,dxm1,dxtm1)
198  do i=1,lxyz
199  ux(i) = jacmi(i,e)*(ur(i)*rxm1(i,1,1,e)
200  $ + us(i)*sxm1(i,1,1,e)
201  $ + ut(i)*txm1(i,1,1,e) )
202  uy(i) = jacmi(i,e)*(ur(i)*rym1(i,1,1,e)
203  $ + us(i)*sym1(i,1,1,e)
204  $ + ut(i)*tym1(i,1,1,e) )
205  uz(i) = jacmi(i,e)*(ur(i)*rzm1(i,1,1,e)
206  $ + us(i)*szm1(i,1,1,e)
207  $ + ut(i)*tzm1(i,1,1,e) )
208  enddo
209  else
210  if (ifaxis) call setaxdy (ifrzer(e))
211  call local_grad2(ur,us,u,n,e,dxm1,dytm1)
212  do i=1,lxyz
213  ux(i) =jacmi(i,e)*(ur(i)*rxm1(i,1,1,e)
214  $ + us(i)*sxm1(i,1,1,e) )
215  uy(i) =jacmi(i,e)*(ur(i)*rym1(i,1,1,e)
216  $ + us(i)*sym1(i,1,1,e) )
217  enddo
218  endif
219 c
220  return
221  end
222 c-----------------------------------------------------------------------
vlsum
real function vlsum(vec, n)
Definition: math.f:417
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
rzero
subroutine rzero(a, n)
Definition: math.f:208
eig2
subroutine eig2(AA, eign, eig1)
Definition: navier2.f:80
aspect_ratios
subroutine aspect_ratios(ar)
Definition: navier2.f:3
quadratic_h
subroutine quadratic_h(x1, x2, a, b, c)
Definition: navier2.f:99
out_sem
subroutine out_sem(iel)
Definition: navier2.f:149
gradm11
subroutine gradm11(ux, uy, uz, u, e)
Definition: navier2.f:176
local_grad2
subroutine local_grad2(ur, us, u, N, e, D, Dt)
Definition: navier5.f:448
local_grad3
subroutine local_grad3(ur, us, ut, u, N, e, D, Dt)
Definition: navier5.f:429
setaxdy
subroutine setaxdy(ifaxdy)
Definition: subs1.f:2342