Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
navier5.F90
Go to the documentation of this file.
1 !-----------------------------------------------------------------------
3 subroutine q_filter(wght)
4  use kinds, only : dp
5  use size_m, only : lx1, lelt
6  use size_m, only : nx1, nz1, nid, ndim
7  use input, only : ifflow, ifbase, if3d, ifsplit, ifldmhd, ifpert, ifheat
8  use input, only : ifcvode, param, ifmhd, iflomach, npscal
9  use soln, only : vx, vy, vz, pr, t
10  use tstep, only : ifield, istep
11  use wz_m, only : zgm1
12  implicit none
13 
14  real(DP) :: wght
15 
16  real(DP), save :: intv(lx1,lx1)
17 
18  real(DP) :: intt(lx1,lx1)
19  real(DP), allocatable :: wk1 (:,:,:,:)
20  real(DP) :: wk2 (lx1,lx1,lx1)
21  real(DP) :: tmax(100),omax(103)
22 
23 ! outpost arrays
24  integer, save :: icalld = 0
25  real(DP), save :: old_wght = 0._dp
26 
27  integer :: imax, jmax, ncut, ifldt, mmax, nfldt, ifld, k
28  integer, external :: iglmax
29  real(DP) :: umax, vmax, wmax, pmax
30  real(DP), external :: glmax
31  logical :: if_fltv
32 
33  imax = nid
34  imax = iglmax(imax,1)
35  jmax = iglmax(imax,1)
36  if (icalld == 0 .or. abs(wght-old_wght) > 1.e-8 ) then
37  icalld = 1
38  old_wght = wght
39  ncut = int(param(101)+1)
40  call build_new_filter(intv,zgm1,nx1,ncut,wght,nid)
41  elseif (icalld < 0) then ! old (std.) filter
42  write(*,*) "Oops: icalld < 0 in qfilter"
43 #if 0
44  icalld = 1
45  call zwgll(zgmv,wgtv,lxv)
46  call igllm(intuv,intw,zgmv,zgm1,lxv,nx1,lxv,nx1)
47  call igllm(intdv,intw,zgm1,zgmv,nx1,lxv,nx1,lxv)
48 
49  call zwgl(zgmp,wgtp,lxp)
50  call iglm(intup,intw,zgmp,zgm2,lxp,nx2,lxp,nx2)
51  call iglm(intdp,intw,zgm2,zgmp,nx2,lxp,nx2,lxp)
52 
53  ! Multiply up and down interpolation into single operator
54 
55  call mxm(intup,nx2,intdp,lxp,intp,nx2)
56  call mxm(intuv,nx1,intdv,lxv,intv,nx1)
57 
58  ! Weight the filter to make it a smooth (as opposed to truncated)
59  ! decay in wave space
60 
61  w0 = 1.-wght
62  call ident(intup,nx2)
63  call add2sxy(intp,wght,intup,w0,nx2*nx2)
64 
65  call ident(intuv,nx1)
66  call add2sxy(intv ,wght,intuv,w0,nx1*nx1)
67 #endif
68  endif
69 
70  ifldt = ifield
71 ! ifield = 1
72 
73  if_fltv = .false.
74  if ( ifflow .AND. .NOT. ifmhd ) if_fltv = .true.
75  if ( ifield == 1 .AND. ifmhd ) if_fltv = .true.
76 
77 ! Adam Peplinski; to take into account freezing of base flow
78  if ( .NOT. ifbase ) if_fltv = .false. ! base-flow frozen
79 
80  allocate(wk1(lx1,lx1,lx1,lelt))
81  if ( if_fltv ) then
82  call filterq(vx,intv,nx1,nz1,wk1,wk2,intt,if3d,umax)
83  call filterq(vy,intv,nx1,nz1,wk1,wk2,intt,if3d,vmax)
84  if (if3d) &
85  call filterq(vz,intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
86  if (ifsplit .AND. .NOT. iflomach) &
87  call filterq(pr,intv,nx1,nz1,wk1,wk2,intt,if3d,pmax)
88  endif
89 
90  if (ifmhd .AND. ifield == ifldmhd) then
91  write(*,*) "Oops: ifmhd"
92 #if 0
93  call filterq(bx,intv,nx1,nz1,wk1,wk2,intt,if3d,umax)
94  call filterq(by,intv,nx1,nz1,wk1,wk2,intt,if3d,vmax)
95  if (if3d) &
96  call filterq(bz,intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
97 #endif
98  endif
99 
100  if (ifpert) then
101  write(*,*) "Oops: ifpert"
102 #if 0
103  do j=1,npert
104 
105  ifield = 1
106  call filterq(vxp(1,j),intv,nx1,nz1,wk1,wk2,intt,if3d,umax)
107  call filterq(vyp(1,j),intv,nx1,nz1,wk1,wk2,intt,if3d,vmax)
108  if (if3d) &
109  call filterq(vzp(1,j),intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
110 
111  ifield = 2
112  if (ifheat .AND. .NOT. ifcvode) &
113  call filterq(tp(1,j,1),intv,nx1,nz1,wk1,wk2,intt,if3d,wmax)
114 
115  enddo
116 #endif
117  endif
118 
119  mmax = 0
120  if (ifflow) then
121  ! pmax = glmax(pmax,1)
122  omax(1) = glmax(umax,1)
123  omax(2) = glmax(vmax,1)
124  omax(3) = glmax(wmax,1)
125  mmax = ndim
126  endif
127 
128 
129  nfldt = 1+npscal
130  if (ifheat .AND. .NOT. ifcvode) then
131  do ifld=1,nfldt
132  ifield = ifld + 1
133  call filterq(t(1,1,1,1,ifld),intv &
134  ,nx1,nz1,wk1,wk2,intt,if3d,tmax(ifld))
135  mmax = mmax+1
136  omax(mmax) = glmax(tmax(ifld),1)
137  enddo
138  endif
139  deallocate(wk1)
140 
141  if (nid == 0) then
142  if (npscal == 0) then
143  ! write(6,101) mmax
144  ! write(sfmt,101) mmax
145  ! 101 format('''(i8,1p',i1,'e12.4,a6)''')
146  ! write(6,sfmt) istep,(omax(k),k=1,mmax),' qfilt'
147  ! write(6,'(i8,1p4e12.4,a6)') istep,(omax(k),k=1,mmax),' qfilt'
148  else
149  if (if3d) then
150  write(6,1) istep,ifield,umax,vmax,wmax
151  else
152  write(6,1) istep,ifield,umax,vmax
153  endif
154  1 format(4x,i7,i3,' qfilt:',1p3e12.4)
155  if(ifheat .AND. .NOT. ifcvode) &
156  write(6,'(1p50e12.4)') (tmax(k),k=1,nfldt)
157  endif
158  endif
159 
160  ifield = ifldt ! RESTORE ifield
161 
162 
163  return
164 end subroutine q_filter
165 
166 !-----------------------------------------------------------------------
167 subroutine filterq(v,f,nx,nz,w1,w2,ft,if3d,dmax)
168  use kinds, only : dp
169  use size_m, only : nelt
170  use tstep, only : nelfld, ifield
171  implicit none
172 
173  integer :: nx, nz
174  real(DP) :: v(nx*nx*nz,nelt),w1(*),w2(*)
175  real(DP) :: f(nx,nx),ft(nx,nx)
176  real(DP) :: dmax
177  logical :: if3d
178 
179 
180  integer :: e, nxyz, nel, i, j, k
181  real(DP) :: smax
182  real(DP), external :: vlamax
183 
184  call transpose(ft,nx,f,nx)
185 
187  nxyz=nx*nx*nz
188  dmax = 0.
189 
190  nel = nelfld(ifield)
191 
192  if (if3d) then
193  do e=1,nel
194  ! Filter
195  call copy(w2,v(1,e),nxyz)
196  call mxm(f,nx,w2,nx,w1,nx*nx)
197  i=1
198  j=1
199  do k=1,nx
200  call mxm(w1(i),nx,ft,nx,w2(j),nx)
201  i = i+nx*nx
202  j = j+nx*nx
203  enddo
204  call mxm(w2,nx*nx,ft,nx,w1,nx)
205  w2(1:nxyz) = v(:,e) - w1(1:nxyz)
206  call copy(v(1,e),w1,nxyz)
207  smax = vlamax(w2,nxyz)
208  dmax = max(dmax,abs(smax))
209  enddo
210 
211  else
212  do e=1,nel
213  ! Filter
214  call copy(w1,v(1,e),nxyz)
215  call mxm(f ,nx,w1,nx,w2,nx)
216  call mxm(w2,nx,ft,nx,w1,nx)
217 
218  w2(1:nxyz) = v(:,e) - w1(1:nxyz)
219  call copy(v(1,e),w1,nxyz)
220  smax = vlamax(w2,nxyz)
221  dmax = max(dmax,abs(smax))
222  enddo
223  endif
224 
225  return
226 end subroutine filterq
227 
228 !-----------------------------------------------------------------------
230 subroutine local_grad3(ur,us,ut,u,N,e,D,Dt)
231  use kinds, only : dp
232  implicit none
233  integer :: N,e
234  real(DP) :: ur(0:n,0:n,0:n),us(0:n,0:n,0:n),ut(0:n,0:n,0:n)
235  real(DP) :: u (0:n,0:n,0:n,*)
236  real(DP) :: D (0:n,0:n),Dt(0:n,0:n)
237 
238  integer :: m1, m2, k
239  m1 = n+1
240  m2 = m1*m1
241 
242  call mxm(d ,m1,u(0,0,0,e),m1,ur,m2)
243  do k=0,n
244  call mxm(u(0,0,k,e),m1,dt,m1,us(0,0,k),m1)
245  enddo
246  call mxm(u(0,0,0,e),m2,dt,m1,ut,m1)
247 
248  return
249 end subroutine local_grad3
250 
251 !-----------------------------------------------------------------------
256 subroutine gaujordf(a,m,n,indr,indc,ipiv,ierr,rmult)
257  use kinds, only : dp
258  implicit none
259 
260  integer :: m, n, indr(m),indc(n),ipiv(n), ierr
261  real(DP) :: a(m,n),rmult(m)
262 
263  real(DP) :: eps
264  integer :: i, j, k, ir, jc
265  real(DP) :: amx, tmp, piv, work
266 
267 ! call outmat(a,m,n,'ab4',n)
268 ! do i=1,m
269 ! write(6,1) (a(i,j),j=1,n)
270 ! enddo
271 ! 1 format('mat: ',1p6e12.4)
272 
273  ierr = 0
274  eps = 1.e-9
275  ipiv(1:m) = 0
276 
277  ir = -1
278  do k=1,m
279  amx=0.
280  do i=1,m ! Pivot search
281  if (ipiv(i) /= 1) then
282  do j=1,m
283  if (ipiv(j) == 0) then
284  if (abs(a(i,j)) >= amx) then
285  amx = abs(a(i,j))
286  ir = i
287  jc = j
288  endif
289  elseif (ipiv(j) > 1) then
290  ierr = -ipiv(j)
291  return
292  endif
293  enddo
294  endif
295  enddo
296  ipiv(jc) = ipiv(jc) + 1
297 
298  ! Swap rows
299  if (ir /= jc) then
300  do j=1,n
301  tmp = a(ir,j)
302  a(ir,j) = a(jc,j)
303  a(jc,j) = tmp
304  enddo
305  endif
306  indr(k)=ir
307  indc(k)=jc
308  ! write(6 ,*) k,' Piv:',jc,a(jc,jc)
309  ! write(28,*) k,' Piv:',jc,a(jc,jc)
310  if (abs(a(jc,jc)) < eps) then
311  write(6 ,*) 'small Gauss Jordan Piv:',jc,a(jc,jc)
312  write(28,*) 'small Gauss Jordan Piv:',jc,a(jc,jc)
313  ierr = jc
314  call exitt
315  return
316  endif
317  piv = 1./a(jc,jc)
318  a(jc,jc)=1.
319  do j=1,n
320  a(jc,j) = a(jc,j)*piv
321  enddo
322 
323  do j=1,n
324  work = a(jc,j)
325  a(jc,j) = a(1 ,j)
326  a(1 ,j) = work
327  enddo
328  do i=2,m
329  rmult(i) = a(i,jc)
330  a(i,jc) = 0.
331  enddo
332 
333  do j=1,n
334  do i=2,m
335  a(i,j) = a(i,j) - rmult(i)*a(1,j)
336  enddo
337  enddo
338 
339  do j=1,n
340  work = a(jc,j)
341  a(jc,j) = a(1 ,j)
342  a(1 ,j) = work
343  enddo
344 
345  ! do i=1,m
346  ! if (i.ne.jc) then
347  ! rmult = a(i,jc)
348  ! a(i,jc) = 0.
349  ! do j=1,n
350  ! a(i,j) = a(i,j) - rmult*a(jc,j)
351  ! enddo
352  ! endif
353  ! enddo
354 
355  enddo
356 
357 ! Unscramble matrix
358  do j=m,1,-1
359  if (indr(j) /= indc(j)) then
360  do i=1,m
361  tmp=a(i,indr(j))
362  a(i,indr(j))=a(i,indc(j))
363  a(i,indc(j))=tmp
364  enddo
365  endif
366  enddo
367 
368  return
369 end subroutine gaujordf
370 
371 !-----------------------------------------------------------------------
373 subroutine legendre_poly(L,x,N)
374  use kinds, only : dp
375  implicit none
376  integer :: N
377  real(DP) :: L(0:n), x
378 
379  integer :: j
380  l(0) = 1.
381  l(1) = x
382 
383  do j=2,n
384  l(j) = ( (2*j-1) * x * l(j-1) - (j-1) * l(j-2) ) / j
385  enddo
386 
387  return
388 end subroutine legendre_poly
389 
390 !-----------------------------------------------------------------------
405 subroutine build_new_filter(intv,zpts,nx,kut,wght,nid)
406  use kinds, only : dp
407  implicit none
408 
409  integer :: nx, nid, kut
410  real(DP) :: intv(nx,nx),zpts(nx), wght
411 
412  integer, parameter :: lm=40
413  integer, parameter :: lm2=lm*lm
414  real(DP) :: phi(lm2),pht(lm2),diag(lm2),rmult(lm),Lj(lm)
415  integer :: indr(lm),indc(lm),ipiv(lm)
416  integer :: kj, n, j, k, k0, kk, ierr, np1
417  real(DP) :: z, amp
418 
419  if (nx > lm) then
420  write(6,*) 'ABORT in build_new_filter:',nx,lm
421  call exitt
422  endif
423 
424  kj = 0
425  n = nx-1
426  do j=1,nx
427  z = zpts(j)
428  call legendre_poly(lj,z,n)
429  kj = kj+1
430  pht(kj) = lj(1)
431  kj = kj+1
432  pht(kj) = lj(2)
433  do k=3,nx
434  kj = kj+1
435  pht(kj) = lj(k)-lj(k-2)
436  enddo
437  enddo
438  call transpose(phi,nx,pht,nx)
439  call copy(pht,phi,nx*nx)
440  call gaujordf(pht,nx,nx,indr,indc,ipiv,ierr,rmult)
441 
442 ! Set up transfer function
443 
444  call ident(diag,nx)
445 
446  k0 = nx-kut
447  do k=k0+1,nx
448  kk = k+nx*(k-1)
449  amp = wght*(k-k0)*(k-k0)/(kut*kut) ! quadratic growth
450  diag(kk) = 1.-amp
451  enddo
452 
453  call mxm(diag,nx,pht,nx,intv,nx) ! -1
454  call mxm(phi ,nx,intv,nx,pht,nx) ! V D V
455  call copy(intv,pht,nx*nx)
456 
457  do k=1,nx*nx
458  pht(k) = 1.-diag(k)
459  enddo
460  np1 = nx+1
461  if (nid == 0) then
462  write(6,6) 'filt amp',(pht(k),k=1,nx*nx,np1)
463  write(6,6) 'filt trn',(diag(k),k=1,nx*nx,np1)
464  6 format(a8,16f7.4,6(/,8x,16f7.4))
465  endif
466 
467  return
468 end subroutine build_new_filter
469 
470 !-----------------------------------------------------------------------
471 real(DP) function ran1(idum)
472  use kinds, only : dp
473  implicit none
474 
475  integer :: idum
476 
477  integer, parameter :: ia=16807,im=2147483647,iq=127773,ir=2836,ntab=32,ndiv=1+(im-1)/ntab
478  real(DP), parameter :: am=1./im,eps=1.2e-7,rnmx=1.-eps
479 
480 ! Numerical Rec. in Fortran, 2nd eD. P. 271
481 
482  integer :: j,k
483  integer :: iv(ntab),iy
484  save iv,iy
485  data iv,iy /ntab*0,0/
486 
487  if (idum <= 0 .OR. iy == 0) then
488  idum=max(-idum,1)
489  do j=ntab+8,1,-1
490  k = idum/iq
491  idum = ia*(idum-k*iq)-ir*k
492  if(idum < 0) idum = idum+im
493  if (j <= ntab) iv(j) = idum
494  enddo
495  iy = iv(1)
496  endif
497  k = idum/iq
498  idum = ia*(idum-k*iq)-ir*k
499  if(idum < 0) idum = idum+im
500  j = 1+iy/ndiv
501  iy = iv(j)
502  iv(j) = idum
503  ran1 = min(am*iy,rnmx)
504 ! ran1 = cos(ran1*1.e8)
505 
506  return
507 end function ran1
508 
509 !-----------------------------------------------------------------------
510 subroutine rand_fld_h1(x)
511  use kinds, only : dp
512  use size_m, only : nx1, ny1, nz1, nelt
513  implicit none
514 
515  real(DP) :: x(*)
516  real(DP), external :: ran1
517 
518  integer :: n, id, i
519  n=nx1*ny1*nz1*nelt
520  id = n
521  do i=1,n
522  x(i) = ran1(id)
523  enddo
524  call dsavg(x)
525 
526  return
527 end subroutine rand_fld_h1
528 
529 !-----------------------------------------------------------------------
tstep
cleaned
Definition: tstep_mod.F90:2
ident
subroutine ident(a, n)
Construct A = I_n (identity matrix)
Definition: math.F90:555
input
Input parameters from preprocessors.
Definition: input_mod.F90:11
soln
Definition: soln_mod.F90:1
wz_m
cleaned
Definition: wz_mod.F90:2
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Compute matrix-matrix product C = A*B.
Definition: mxm_wrapper.F90:7
speclib::iglm
subroutine iglm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
Definition: speclib.F90:1167
exitt
void exitt()
Definition: comm_mpi.F90:411
dsavg
subroutine dsavg(u)
Take direct stiffness avg of u.
Definition: ic.F90:1559
copy
subroutine copy(a, b, n)
Definition: math.F90:52
speclib::zwgll
subroutine, public zwgll(Z, W, NP)
Definition: speclib.F90:95
speclib::zwgl
subroutine, public zwgl(Z, W, NP)
Generate NP Gauss Legendre points (Z) and weights (W) associated with Jacobi polynomial P(N)(alpha=0...
Definition: speclib.F90:77
speclib::igllm
subroutine, public igllm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
Definition: speclib.F90:1197
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8